## ANOCOVA
## example: Neter et al. page 871

y<-c(38,39,36,45,33, 43,38,38,27,34, 24,32,31,21,28)
x<-c(21,26,22,28,19, 34,26,29,18,25, 23,29,30,16,29)
%trt<-c(rep(1,5),rep(2,5),rep(3,5))
trt<-rep(1:3,c(5,5,5))

##############
# plot y vs. x
##############
plot(x[trt==1],y[trt==1],xlim=c(15,35),ylim=c(20,50),xlab='x',ylab='y')
points(x[trt==2],y[trt==2],pch=20)
points(x[trt==3],y[trt==3],pch=22)
legend(15,50,c('treatment 1', 'treatment 2', 'treatment 3'),
     pch=c(21,20,22))
## from the figure we see that the relationship between y and x is linear.
## the slopes may be different for the 3 trts, but are more or less parallel.

##########################
# testing for treat effect
########################## 
> contrasts(factor(trt))
  2 3
1 0 0
2 1 0
3 0 1
## treatment contrasts, not the desired form of design matrix

ind1<-rep(-1,15)
ind1<-ifelse(trt==1,1,ind1)
ind1<-ifelse(trt==2,0,ind1)
ind2<-rep(-1,15)
ind2<-ifelse(trt==1,0,ind2)
ind2<-ifelse(trt==2,1,ind2)
xnew<-x-mean(x)
data<-cbind(y,ind1,ind2,xnew)
       y ind1 ind2  xnew
 [1,] 38    1    0    -4
 [2,] 39    1    0     1
 [3,] 36    1    0    -3
 [4,] 45    1    0     3
 [5,] 33    1    0    -6
 [6,] 43    0    1     9
 [7,] 38    0    1     1
 [8,] 38    0    1     4
 [9,] 27    0    1    -7
[10,] 34    0    1     0
[11,] 24   -1   -1    -2
[12,] 32   -1   -1     4
[13,] 31   -1   -1     5
[14,] 21   -1   -1    -9
[15,] 28   -1   -1     4

fit1<-lm(y~xnew+ind1+ind2)   ## full model
summary(fit1,corr=T)
Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  33.8000     0.4835  69.908 6.37e-16 ***
ind1          6.0174     0.7083   8.496 3.67e-06 ***
ind2          0.9420     0.6987   1.348    0.205    
xnew          0.8986     0.1026   8.759 2.73e-06 ***

Correlation of Coefficients:
      (Intercept) ind1  ind2 
ind1   0.00                  
ind2   0.00       -0.53      
xnew   0.00        0.26 -0.21

anova(fit1)
Analysis of Variance Table
          Df Sum Sq Mean Sq  F value    Pr(>F)    
xnew       1 190.68  190.68  54.3786 1.405e-05 ***
ind1       1 410.78  410.78 117.1478 3.335e-07 ***
ind2       1   6.37    6.37   1.8178    0.2047    
Residuals 11  38.57    3.51                       

fit2<-lm(y~xnew)   ## reduced model
Analysis of Variance Table
          Df Sum Sq Mean Sq F value  Pr(>F)  
xnew       1 190.68  190.68  5.4393 0.03641 *
Residuals 13 455.72   35.06                  

#######################
# test for equal slopes 
#######################
## For this purpose, the full model above is the reduced model
## and model with interaction is the full model
fit3<-lm(y~ind1+ind2+xnew+ind1:xnew+ind2:xnew)
Analysis of Variance Table
          Df  Sum Sq Mean Sq F value    Pr(>F)    
ind1       1 302.500 302.500 86.3714 6.561e-06 ***
ind2       1  36.300  36.300 10.3646    0.0105 *  
xnew       1 269.029 269.029 76.8145 1.060e-05 ***
ind1:xnew  1   6.595   6.595  1.8830    0.2032    
ind2:xnew  1   0.456   0.456  0.1301    0.7267    
Residuals  9  31.521   3.502                      

###############
# ANOCOVA table
###############
fit4<-lm(y~factor(trt))
anova(fit4)
Analysis of Variance Table
            Df Sum of Sq  Mean Sq  F Value      Pr(F) 
factor(trt)  2     338.8 169.4000 6.608583 0.01161208
  Residuals 12     307.6  25.6333                    
## get SSTO and SSE for y here
> 338.8+307.6
[1] 646.4

fit5<-lm(x~factor(trt))
            Df Sum of Sq  Mean Sq  F Value     Pr(F) 
factor(trt)  2      26.8 13.40000 0.482593 0.6286587
  Residuals 12     333.2 27.76667                   
> 26.8+333.2
[1] 360
## get SSTO and SSE for x here

## SPTO 
> sum( (x-mean(x))*(y-mean(y)) )
[1] 262
## SPE
meanx<-rep( c(mean(x[trt==1]),mean(x[trt==2]),mean(x[trt==3])), c(5,5,5) )
meany<-rep( c(mean(y[trt==1]),mean(y[trt==2]),mean(y[trt==3])), c(5,5,5) )
> sum( (x-meanx)*(y-meany) )
[1] 299.4