[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.23983, 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.72703, 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.00397, 0.00447, 0.01705, 0.0139, 0.0129, 0.0081, 0.00383, 0.00575, 0.01127, 0.0103, 0.009, 0.0052, 0.00595, 0.0105, 0.0091, 0.00905, 0.0045, 0.0031, 0.00647, 0.0053, 0.00977, 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.
Re: [R] A-priori contrasts with type III sums of squares in R
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.23983, 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.72703, 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.00397, > 0.00447, > 0.01705, 0.0139, 0.0129, 0.0081, 0.00383, 0.00575, 0.01127, > 0.0103, 0.009, 0.0052, 0.00595, 0.0105, 0.0091, 0.
Re: [R] A-priori contrasts with type III sums of squares in R
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:bl...@nceas.ucsb.edu] > Sent: June-09-15 5:14 PM > To: John Fox; r-help@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-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
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/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-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, 1