#### Stratification to control for prognostic variables
##stratum 1
y11<-c(42,44,36,13,19,22)
y12<-c(28,23,34,42,13)
y13<-c(1,29,19)
y14<-c(24,9,22,-2,15)
##stratum 2
y21<-c(33,26,33,21)
y22<-c(34,33,31,36)
y23<-c(11,9,7,1,-6)
y24<-c(27,12,12,-5,16,15)
##stratum 3
y31<-c(31,-3,25,25,24)
y32<-c(3,26,28,32,4,16)
y33<-c(21,1,9,3)
y34<-c(22,7,25,5,12)
## table 6.10
tb610<-rbind(c(length(y11),length(y12),length(y13),length(y14)),
c(mean(y11),mean(y12),mean(y13),mean(y14)),
c(sqrt(var(y11)),sqrt(var(y12)),sqrt(var(y13)),sqrt(var(y14))),
c(length(y21),length(y22),length(y23),length(y24)),
c(mean(y21),mean(y22),mean(y23),mean(y24)),
c(sqrt(var(y21)),sqrt(var(y22)),sqrt(var(y23)),sqrt(var(y24))),
c(length(y31),length(y32),length(y33),length(y34)),
c(mean(y31),mean(y32),mean(y33),mean(y34)),
c(sqrt(var(y31)),sqrt(var(y32)),sqrt(var(y33)),sqrt(var(y34))))
> round(tb610,4)
        [,1]    [,2]    [,3]    [,4] 
[1,]  6.0000  5.0000  3.0000  5.0000
[2,] 29.3333 28.0000 16.3333 13.6000
[3,] 13.0179 10.9772 14.1892 10.5499
[4,]  4.0000  4.0000  5.0000  6.0000
[5,] 28.2500 33.5000  4.4000 12.8333
[6,]  5.8523  2.0817  6.9138 10.3425
[7,]  5.0000  6.0000  4.0000  5.0000
[8,] 20.4000 18.1667  8.5000 14.2000
[9,] 13.3716 12.5286  9.0000  8.9275

y<-c(y11,y12,y13,y14,y21,y22,y23,y24,y31,y32,y33,y34)
trt<-c(rep(1:4,c(length(y11),length(y12),length(y13),length(y14))),
rep(1:4,c(length(y21),length(y22),length(y23),length(y24))),
rep(1:4,c(length(y31),length(y32),length(y33),length(y34))))
strata<-c(rep(1,length(c(y11,y12,y13,y14))),
rep(2,length(c(y21,y22,y23,y24))),
rep(3,length(c(y31,y32,y33,y34))))

fit1<-lm(y~factor(strata)*factor(trt),subset=(trt<=2))
                           Df  Sum Sq Mean Sq F value  Pr(>F)  
factor(strata)              2  782.11  391.05  3.1867 0.05924 .
factor(trt)                 1    0.08    0.08  0.0007 0.97983  
factor(strata):factor(trt)  2   73.50   36.75  0.2995 0.74394  
Residuals                  24 2945.12  122.71                  

fit2<-lm(y~factor(strata)+factor(trt),subset=(trt<=2))
               Df  Sum Sq Mean Sq F value  Pr(>F)  
factor(strata)  2  782.11  391.05  3.3682 0.05003 .
factor(trt)     1    0.08    0.08  0.0007 0.97925  
Residuals      26 3018.61  116.10              
## for trt effect, trt should be entered after strata.

fit3<-lm(y~factor(strata)*factor(trt))
                           Df Sum Sq Mean Sq F value    Pr(>F)    
factor(strata)              2  488.6   244.3  2.2120    0.1210    
factor(trt)                 3 3063.4  1021.1  9.2451 6.748e-05 ***
factor(strata):factor(trt)  6  707.3   117.9  1.0672    0.3958    
Residuals                  46 5080.8   110.5                    
## compare to the results in table 6.12
## If interested in testing for trt effect, strata needs to be entered first.