[R] A-priori contrasts with type III sums of squares in R

[Rachael Blake]

I am analyzing data using a factorial three-way ANOVA with a-priori 
contrasts and type III sums of squares. (Please don't comment about type 
I SS vs. type III SS. That's not the point of my question.  I have read 
at length about the choice between types of SS and have made my 
decision.) I get the contrasts like I need using summary.aov(), however 
that uses type I SS. When I use the Anova() function from library(car) 
to get type III SS, I don't get the contrasts. I have also tried using 
drop1() with the lm() model, but I get the same results as Anova() 
(without the contrasts).

Please advise on a statistical method in R to analyze data using 
factorial ANOVA with a-priori contrasts and type III SS as shown in my 
example below.

Sample data:
     DF <- structure(list(Code = structure(c(1L, 1L, 1L, 2L, 2L, 2L, 3L, 
3L,
     3L, 4L, 4L, 4L, 5L, 5L, 5L, 6L, 6L, 6L, 7L, 7L, 7L, 8L, 8L, 8L, 9L, 
9L,
     9L, 10L, 10L, 10L, 11L, 11L, 11L, 12L, 12L, 12L), .Label = c("A",
     "B", "C", "D", "E", "F", "G", "H", "I", "J", "K", "L"), class =
     "factor"), GzrTreat = structure(c(3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
     3L, 3L, 3L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,  2L, 2L,
     2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), contrasts = structure(c(1,
     -2, 1, 1, 0, -1), .Dim = c(3L, 2L), .Dimnames = list(c("I",
     "N", "R"), NULL)), .Label = c("I", "N", "R"), class = "factor"),
     BugTreat = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
     1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
     3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L), .Label =
     c("Immigration", "Initial", "None"), class = "factor"), TempTreat =
     structure(c(2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 
2L,
     2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 1L,
     1L, 1L, 1L, 1L, 1L), .Label = c("Not Warm", "Warmed"), class =
     "factor"), ShadeTreat = structure(c(2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L,
     2L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L,
     1L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L), .Label = 
c("Light",
     "Shaded"), class = "factor"), EpiChla = c(0.268482353, 0.423119608,
     0.579507843, 0.738839216, 0.727856863, 0.523960784, 0.405801961,
     0.335964706, 0.584441176, 0.557543137, 0.436456863, 0.563909804,
     0.432398039, 0.344956863, 0.340309804, 0.992884314, 0.938390196,
     0.663270588, 0.239833333, 0.62875098, 0.466011765, 0.536182353,
     0.340309804, 0.721172549, 0.752082353, 0.269372549, 0.198180392,
     1.298882353, 0.298354902, 0.913139216, 0.846129412, 0.922317647,
     0.727033333, 1.187662745, 0.35622549, 0.073547059), log_EpiChla =
     c(0.10328443, 0.153241402, 0.198521787, 0.240259426, 0.237507762,
     0.182973791, 0.147924145, 0.125794985, 0.19987612, 0.192440084,
     0.157292589, 0.194211702, 0.156063718, 0.128708355, 0.127205194,
     0.299482089, 0.287441205, 0.220962908, 0.093363308, 0.21185469,
     0.166137456, 0.186442772, 0.127205194, 0.235824411, 0.243554515,
     0.103589102, 0.078522208, 0.361516746, 0.113393422, 0.281746574,
     0.266262141, 0.283825153, 0.23730072, 0.339980371, 0.132331903,
     0.030821087), MeanZGrowthAFDM_g = c(0.00665, 0.003966667, 0.004466667,
     0.01705, 0.0139, 0.0129, 0.0081, 0.003833333, 0.00575, 0.011266667,
     0.0103, 0.009, 0.0052, 0.00595, 0.0105, 0.0091, 0.00905, 0.0045, 
0.0031,
     0.006466667, 0.0053, 0.009766667, 0.0181, 0.00725, 0, 0.0012, 5e-04,
     0.0076, 0.00615, 0.0814, NA, 0.0038, 0.00165, 0.0046, 0, 0.0015)),
     .Names = c("Code", "GzrTreat", "BugTreat", "TempTreat", "ShadeTreat",
     "EpiChla", "log_EpiChla", "MeanZGrowthAFDM_g"), class = "data.frame",
     row.names = c(NA, -36L))


Code:

     ## a-priori contrasts
     library(stats)
     contrasts(DF$GzrTreat) <- cbind(c(1,-2,1), c(1,0,-1))
     round(crossprod(contrasts(DF$GzrTreat)))
     c_labels <- list(GzrTreat=list('presence'=1, 'immigration'=2))

     ## model
     library(car)
     EpiLM <- lm(log_EpiChla~TempTreat*GzrTreat*ShadeTreat, DF)
     summary.aov(EpiLM, split=c_labels) ### MUST USE summary.aov(), to get
     #contrast results, but sadly this uses Type I SS
     Anova(EpiLM, split=c_labels, type="III") # Uses Type III SS, but NO
     #CONTRASTS!!!!!
     drop1(EpiLM, ~., test="F") # again, this does not print contrasts

     # I need contrast results like from summary.aov(), AND Type III SS
     # like from Anova()



-- 
Rachael E. Blake, PhD
Post-doctoral Associate