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

John Fox jfox at mcmaster.ca
Wed Jun 10 00:18:46 CEST 2015


Dear Rachel,

How about this (using the data and model you sent originally)?

> linearHypothesis(EpiLM, "GzrTreatpresence = 0")
Linear hypothesis test

Hypothesis:
GzrTreatpresence = 0

Model 1: restricted model
Model 2: log_EpiChla ~ TempTreat * GzrTreat * ShadeTreat

  Res.Df     RSS Df  Sum of Sq      F Pr(>F)
1     25 0.12665                            
2     24 0.12623  1 0.00042195 0.0802 0.7794

> linearHypothesis(EpiLM, "GzrTreatimmigration = 0")
Linear hypothesis test

Hypothesis:
GzrTreatimmigration = 0

Model 1: restricted model
Model 2: log_EpiChla ~ TempTreat * GzrTreat * ShadeTreat

  Res.Df     RSS Df  Sum of Sq     F Pr(>F)
1     25 0.12623                           
2     24 0.12623  1 5.0931e-06 0.001 0.9754

Note that this tests main-effect contrasts in a model that includes
interactions to which the main effect is marginal. You should probably think
about whether you really want to do that.

BTW, the slides to which you refer are for *multivariate* linear models
(including repeated measures); you're using a univariate linear model.

Best,
 John

> -----Original Message-----
> From: Rachael Blake [mailto:blake at nceas.ucsb.edu]
> Sent: June-09-15 5:14 PM
> To: John Fox; r-help at r-project.org
> Subject: Re: [R] A-priori contrasts with type III sums of squares in R
> 
> 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/17Au
> g_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


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



More information about the R-help mailing list