thr3ads.net - R help - [R] A-priori contrasts with type III sums of squares in R [Jun 2015]

If this information is useful, please help other people find it:
Share via:

Rachael Blake

2015-Jun-05 22:32 UTC

[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

John Fox

2015-Jun-06 19:35 UTC

head link

[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-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
> 
> ______________________________________________
> R-help at 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

Rachael Blake

2015-Jun-09 21:14 UTC

head link

[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

R help - Jun 2015 - A-priori contrasts with type III sums of squares in R

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

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

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