[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
