Rachael Blake
2015-Jun-09 21:14 UTC
[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-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
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