Anova-Oneway

# oneway.r
# code for 1-way ANOVA examples in R

diet <- read.csv('diet.csv', as.is=T)
diet <- diet[diet$diet !='lopro',]

# I use as.is=T to force me to explicitly declare 
#   factors. default conversion does not convert 
#   factors with numeric levels.

diet$diet.f <- factor(diet$diet)

# anova using lm:
diet.lm <- lm(longevity ~ diet.f, data=diet)

# lm() has lots of helper / reporting functions:
coef(diet.lm)         
   # coefficients
vcov(diet.lm)         
   # and their variance-covariance matrix
sqrt(diag(vcov(diet.lm)))  
   # se's of the coefficients
anova(diet.lm)        
   # ANOVA table using type I = sequential SS
summary(diet.lm)      
   # lots of information
plot(predict(diet.lm), resid(diet.lm))   
   # plot of residuals vs predicted values

# model comparison by hand:
diet.m0 <- lm(longevity ~ +1, data=diet)   
   # intercept only model
anova(diet.m0, diet.lm)    
   # change in SS from 1st to 2nd model

drop1(diet.lm)         
   # drop one term at a time = type II SS

# orthogonal contrasts
diet.means <- tapply(diet$longevity, diet$diet, mean)  
   # mean for each diet
diet.helm <- contr.helmert(5)   
   # matrix of Helmert coeff. for 5 groups
diet.c1 <- t(diet.helm) %*% diet.means   
   # estimate of each contrast
diet.ss1 <- 49*diet.c1^2/apply(diet.helm^2,2,sum)  
   # SS for each contrast
sum(diet.ss1)

# second set of contrasts were hand entered
#  using diet.2nd <- rbind(c(-2,-1,0,1,2), ...)

# third set: 
diet.trt <- -contr.treatment(5)
diet.trt[1,] <- 1

diet.c3 <- t(diet.trt) %*% diet.means   
   # estimate of each contrast
diet.ss3 <- 49*diet.c3^2/apply(diet.trt^2,2,sum)  
   # SS for each contrast
sum(diet.ss3)

# C beta tests
diet.c1 <- t(diet.helm) %*% diet.means   
   # estimate of contrast

X1 <- model.matrix(~ -1 + diet.f, data=diet)
   # X matrix using cell means parameterization
diet.cm1 <- t(diet.helm) %*% solve(t(X1) %*% X1) %*% 
   diet.helm

t(diet.c1) %*% solve(diet.cm1) %*% diet.c1

# how do the explicit contrasts of cell means 
#   compare to R using those contrasts?

contrasts(diet$diet.f) <- contr.helmert  
#  tell R to use helmert contrasts 
#  default is contr.treatment, 
#    which is drop first level of the factor
#  contr.SAS is SAS-like (drop last level)

diet.lmh <- lm(longevity ~ diet.f,data=diet)
coef(diet.lmh)
# coefficients are not the same as the 
#   hand-computed contrasts!

# the R coefficients are related to the
#   contrasts among cell means

apply(diet.helm^2,2,sum)
   # sum of squared coefficients
   
coef(diet.lmh)[2:5] * c(2,6,12,20)
# these are the same as diet.c1

# although R calls them contrasts, they really are not.
# they are columns of the X matrix for a regression.
model.matrix(diet.lmh)

# can do C beta from the lm() information

diet.mse <- 4687.7/240   # MSE = est of sigma^2
diet.clm <- coef(diet.lm)[2:5]
  # extract coefficients for diet factor
diet.vclm <- vcov(diet.lm)[2:5,2:5]
  # and their VC matrix
diet.mse * t(diet.clm) %*% solve(diet.vclm) %*% 
    diet.clm
  # SS for diet
t(diet.clm) %*% solve(diet.vclm) %*% diet.clm/5
  # F statistic for diet (5 is # rows in C)