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

[Rachael Blake]

Thank you for replying, John!

I am not using treatment contrasts in this analysis.  I am specifying
           options(contrasts=c("contr.sum", "contr.poly"))
earlier in my code in order to get interpretable results from the Type 
III SS.  However, I did not include that code in the example because it 
is not related to my initial question, and those contrasts are not of 
interest to me.  My interest is in my a-priori specified contrasts:
          contrasts(All09$GzrTreat) <- cbind('presence'=c(1,-2,1), 
'immigration'=c(1,0,-1))

I have made a valiant attempt to use linearHypothesis(), based on the 
example provided here
https://web.warwick.ac.uk/statsdept/user2011/TalkSlides/Contributed/17Aug_1705_FocusV_4-Multivariate_1-Fox.pdf
as well as other places.   I have tried two different ways of specifying 
my contrast matrix, but I keep getting error messages that I can not 
resolve.   My code based on that powerpoint presentation is as follows 
(still using the data included in my initial question):

         options(contrasts=c("contr.sum", "contr.poly"))
         EpiLM <- lm(log_EpiChla~TempTreat*GzrTreat*ShadeTreat, All09)
         Anova(EpiLM, type="III")
         class(EpiLM)
         contrasts(All09$GzrTreat) <- cbind('presence'=c(1,-2,1), 
'immigration'=c(1,0,-1))
         con <- contrasts(All09$GzrTreat) ; con
         EpiLM2 <- update(EpiLM)
         rownames(coef(EpiLM2))
         linearHypothesis(model=EpiLM2, 
hypothesis.matrix=c("presence","immigration"), verbose=T)  # first 
attempt to implement
         linearHypothesis(model=EpiLM2, hypothesis.matrix=con, 
verbose=T)                                              # second attempt 
to implement


Thanks again for your reply.

-Rachael


On 6/6/2015 12:35 PM, John Fox wrote:
> Dear Rachel,
>
> Anova() won't give you a breakdown of the SS for each term into 1 df
> components (there is no split argument, as you can see if you look at
> ?Anova). Because, with the exception of GzrTreat, your contrasts are not
> orthogonal in the row basis of the design (apparently you're using the
> default "contr.treatment" coding), you also won't get sensible type-III
> tests from Anova(). If you formulated the contrasts for the other factors
> properly (using, e.g., contr.sum), you could get single df tests from
> linearHypothesis() in the car package.
>
> I hope this helps,
>   John
>
> -----------------------------------------------
> John Fox, Professor
> McMaster University
> Hamilton, Ontario, Canada
> http://socserv.socsci.mcmaster.ca/jfox/
>
>
>
>
>> -----Original Message-----
>> From: R-help [mailto:r-help-bounces at r-project.org] On Behalf Of Rachael
>> Blake
>> Sent: June-05-15 6:32 PM
>> To: r-help at r-project.org
>> Subject: [R] A-priori contrasts with type III sums of squares in R
>>
>> 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
>>
>>

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