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 > > ______________________________________________ > 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. --- This email has been checked for viruses by Avast antivirus software. https://www.avast.com/antivirus ______________________________________________ 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.
