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

[John Fox]

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
> 
> ______________________________________________
> R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-
> guide.html
> and provide commented, minimal, self-contained, reproducible code.


---
This email has been checked for viruses by Avast antivirus software.
https://www.avast.com/antivirus