#### repeated measures

####################
# data and table 8.2
####################
y11<-c(30.9,31.9,31.3,32.1,30.9,31.3,31.3,32.1,30.3,32.2)
y12<-c(30.7,31.6,31.1,31.0,31.2,31.7,31.8,33.0,30.9,32.1)
y13<-c(30.9,31.6,31.0,31.7,30.5,31.4,31.8,31.7,30.8,32.2)
y14<-c(30.9,31.7,31.3,31.3,30.8,31.2,31.7,31.5,30.6,32.4)
y1<-cbind(y11,y12,y13,y14)
> apply(y1,1,mean)
 [1] 30.850 31.700 31.175 31.525 30.850 31.400 31.650 32.075 30.650 32.225
> apply(y1,2,mean)
  y11   y12   y13   y14 
31.43 31.51 31.36 31.34 
> round(apply(y1,2,sd),2)
 y11  y12  y13  y14 
0.63 0.69 0.54 0.52 
> mean(y1)
[1] 31.41

y21<-c(31.5,31.2,31.3,30.4,30.7,29.8,31.4,30.9,31.1,31.3)
y22<-c(30.6,31.2,31.3,30.8,30.9,30.8,32.0,32.4,31.3,31.5)
y23<-c(30.8,31.1,31.5,30.4,30.9,30.9,31.7,31.8,31.2,31.6)
y24<-c(31.0,31.3,31.4,30.2,30.9,30.8,31.6,31.9,31.2,31.7)
y2<-cbind(y21,y22,y23,y24)
> apply(y2,1,mean)
 [1] 30.975 31.200 31.375 30.450 30.850 30.575 31.675 31.750 31.200 31.525
> apply(y2,2,mean)
  y21   y22   y23   y24 
30.96 31.28 31.19 31.20 
> round(apply(y2,2,sd),2)
 y21  y22  y23  y24 
0.53 0.57 0.45 0.50 
> mean(y2)
[1] 31.1575

> mean(c(y11,y21))
[1] 31.195
> mean(c(y12,y22))
[1] 31.395
> mean(c(y13,y23))
[1] 31.275
> mean(c(y14,y24))
[1] 31.27
> mean(c(y1,y2))
[1] 31.28375

########################
# ANOVA table: Table 8.3
########################
y<-c(y1,y2)
grp<-rep(1:2,c(40,40))
subj<-c(rep(1:10,4),rep(11:20,4))
time<-rep(rep(1:4,rep(10,4)),2)

fit1<-lm(y~factor(grp)*factor(time)+factor(subj))
> anova(fit1)
                         Df  Sum Sq Mean Sq F value 
factor(grp)               1  1.2751  1.2751 13.1638 
factor(time)              3  0.4104  0.1368  1.4122   # F2
factor(subj)             18 17.1763  0.9542  9.8511 
factor(grp):factor(time)  3  0.3364  0.1121  1.1575   # F3 
Residuals                54  5.2307  0.0969                      
> 1.2751/0.9542
[1] 1.336303     # F1

## averaging over time
yy<-c(apply(y1,1,mean),apply(y2,1,mean))
trt<-rep(1:2,c(10,10))
fit2<-lm(yy~factor(trt))
> anova(fit2)
            Df Sum Sq Mean Sq F value Pr(>F)
factor(trt)  1 0.3188  0.3188  1.3363 0.2628  # same as F1, correct F test
Residuals   18 4.2941  0.2386    

#############################
# expression 8.1, Problem 8.2
#############################
s1<-round(var(y1),4)
> s1
       y11    y12    y13    y14
y11 0.4001 0.2697 0.2813 0.2664
y12 0.2697 0.4721 0.2349 0.2207
y13 0.2813 0.2349 0.2871 0.2518
y14 0.2664 0.2207 0.2518 0.2738
s2<-round(var(y2),4)
> s2
       y21    y22    y23    y24
y21 0.2804 0.0969 0.1173 0.1544
y22 0.0969 0.3218 0.2276 0.2333
y23 0.1173 0.2276 0.2054 0.2178
y24 0.1544 0.2333 0.2178 0.2489

##pooled
rr<-round((s1+s2)/2,4)   ## simpler formula because n1=n2
       y11    y12    y13    y14
y11 0.3402 0.1833 0.1993 0.2104
y12 0.1833 0.3970 0.2313 0.2270
y13 0.1993 0.2313 0.2463 0.2348
y14 0.2104 0.2270 0.2348 0.2613

> 0.1833/sqrt(0.3402*0.3970)
[1] 0.4987697
> 0.1993/sqrt(0.3402*0.2463)
[1] 0.6885063
> 0.2104/sqrt(0.3402*0.2613)
[1] 0.7056817
> 0.2313/sqrt(0.3970*0.2463)
[1] 0.7396873
> 0.2270/sqrt(0.3970*0.2613)
[1] 0.7047919
> 0.2348/sqrt(0.2463*0.2613)
[1] 0.9255421