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

Rachael Blake Tue, 09 Jun 2015 23:32:44 -0700

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 TypeIII SS. However, I did not include that code in the example because itis not related to my initial question, and those contrasts are not ofinterest 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 theexample 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 specifyingmy contrast matrix, but I keep getting error messages that I can notresolve. 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) # firstattempt to implementlinearHypothesis(model=EpiLM2, hypothesis.matrix=con,verbose=T) # second attemptto 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-boun...@r-project.org] On Behalf Of Rachael
Blake
Sent: June-05-15 6:32 PM
To: r-help@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

______________________________________________
R-help@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.

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

Reply via email to