I have some doubts about the validity of my procedure to estimeate linear contrasts ina a factorial design. For sake of semplicity, let's imagine a one way ANOVA with three levels. I am interested to test the significance of the difference between the first and third level (called here contrast C1) and between the first and the seconda level (called here contrast C2). I used the following procedure: ------------------- reading data from a text file -----------------------> ar <-read.table("C:/Programmi/R/myworks/contrasti/cont1.txt",header=TRUE)> arCC GROUP 1 3.0 0 2 3.0 0 3 4.0 0 4 5.0 0 5 6.0 0 6 7.0 0 7 3.0 0 8 2.0 0 9 1.0 1 10 6.0 1 11 5.0 1 12 7.0 1 13 2.0 1 14 3.0 1 15 1.5 1 16 1.7 1 17 17.0 2 18 12.0 2 19 15.0 2 20 16.0 2 21 12.0 2 22 23.0 2 23 19.0 2 24 21.0 2 ------------------- creating a new array of data-----------------------> ar<-data.frame(GROUP=factor(ar$GROUP),DIP=ar$CC)> arGROUP DIP 1 0 3.0 2 0 3.0 3 0 4.0 4 0 5.0 5 0 6.0 6 0 7.0 7 0 3.0 8 0 2.0 9 1 1.0 10 1 6.0 11 1 5.0 12 1 7.0 13 1 2.0 14 1 3.0 15 1 1.5 16 1 1.7 17 2 17.0 18 2 12.0 19 2 15.0 20 2 16.0 21 2 12.0 22 2 23.0 23 2 19.0 24 2 21.0 ------------------- creating two dummy variables (C1 and C2) for linear contrasts-----------------------> ar<-data.frame(GROUP=factor(ar$GROUP),C1=factor(ar$GROUP),C2=factor(ar$GROUP),DIP=ar$DIP)> arGROUP C1 C2 DIP 1 0 0 0 3.0 2 0 0 0 3.0 3 0 0 0 4.0 4 0 0 0 5.0 5 0 0 0 6.0 6 0 0 0 7.0 7 0 0 0 3.0 8 0 0 0 2.0 9 1 1 1 1.0 10 1 1 1 6.0 11 1 1 1 5.0 12 1 1 1 7.0 13 1 1 1 2.0 14 1 1 1 3.0 15 1 1 1 1.5 16 1 1 1 1.7 17 2 2 2 17.0 18 2 2 2 12.0 19 2 2 2 15.0 20 2 2 2 16.0 21 2 2 2 12.0 22 2 2 2 23.0 23 2 2 2 19.0 24 2 2 2 21.0 ------------------- selecting the contrast levels-----------------------> ar$C1 <- C(ar$C1, c(1,0,-1), how.many = 1)> ar$C2 <- C(ar$C2, c(1,-1,0), how.many = 1)------------------- contrast analysis of C2 -----------------------> r.aov8 <-aov(DIP ~ C2 + GROUP , data = ar)> anova(r.aov8)Analysis of Variance Table Response: DIP Df Sum Sq Mean Sq F value Pr(>F) C2 1 2.10 2.10 0.2622 0.614 GROUP 1 917.00 917.00 114.3460 5.915e-10 *** Residuals 21 168.41 8.02 --- Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 ------------------- contrast analysis of C1 -----------------------> r.aov9 <-aov(DIP ~ C1 + GROUP , data = ar)> anova(r.aov9)Analysis of Variance Table Response: DIP Df Sum Sq Mean Sq F value Pr(>F) C1 1 650.25 650.25 81.083 1.175e-08 *** GROUP 1 268.85 268.85 33.525 9.532e-06 *** Residuals 21 168.41 8.02 --- Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 ------------------- anova of the global design -----------------------> r.aov10 <-aov(DIP ~ GROUP , data = ar)> anova(r.aov10)Analysis of Variance Table Response: DIP Df Sum Sq Mean Sq F value Pr(>F) GROUP 2 919.10 459.55 57.304 3.121e-09 *** Residuals 21 168.41 8.02 --- Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 I would like to know if there is a more economic procedure with R to do linear contrasts. Every comments will be well accepted. Thank you very much and best regards Marco Tommasi [[alternative HTML version deleted]]
group=factor(rep(c(0:2), each = 8))
ar = data.frame(group, dip)
con = matrix(c(1, -1, 0, 1, 0, -1), nrow=3, ncol=2, byrow=F)
contrasts(ar$group)=con
aovRes = aov(dip~group, ar)
> summary.aov(aovRes, split=list(group = list("0 vs 1" = 1,
"0 vs 3" = 2)))
Df Sum Sq Mean Sq F value Pr(>F)
group 2 919.10 459.55 57.3041 3.121e-09 ***
group: 0 vs 1 1 2.10 2.10 0.2622 0.614
group: 0 vs 3 1 917.00 917.00 114.3460 5.915e-10 ***
Residuals 21 168.41 8.02
---
Signif. codes: 0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1
I only don't know why it does not work with within-subject designs.
I have posted this question before, does anybody know?
-steffen
> I have some doubts about the validity of my procedure to estimeate linear
contrasts ina a factorial design.
> For sake of semplicity, let's imagine a one way ANOVA with three
levels. I am interested to test the significance of the difference between the
first and third level (called here contrast C1) and between the first and the
seconda level (called here contrast C2). I used the following procedure:
>
>
> ------------------- reading data from a text file -----------------------
>
>> ar
<-read.table("C:/Programmi/R/myworks/contrasti/cont1.txt",header=TRUE)
>
>> ar
>
> CC GROUP
>
> 1 3.0 0
>
> 2 3.0 0
>
> 3 4.0 0
>
> 4 5.0 0
>
> 5 6.0 0
>
> 6 7.0 0
>
> 7 3.0 0
>
> 8 2.0 0
>
> 9 1.0 1
>
> 10 6.0 1
>
> 11 5.0 1
>
> 12 7.0 1
>
> 13 2.0 1
>
> 14 3.0 1
>
> 15 1.5 1
>
> 16 1.7 1
>
> 17 17.0 2
>
> 18 12.0 2
>
> 19 15.0 2
>
> 20 16.0 2
>
> 21 12.0 2
>
> 22 23.0 2
>
> 23 19.0 2
>
> 24 21.0 2
>
>
>
> ------------------- creating a new array of data-----------------------
>
>> ar<-data.frame(GROUP=factor(ar$GROUP),DIP=ar$CC)
>
>> ar
>
> GROUP DIP
>
> 1 0 3.0
>
> 2 0 3.0
>
> 3 0 4.0
>
> 4 0 5.0
>
> 5 0 6.0
>
> 6 0 7.0
>
> 7 0 3.0
>
> 8 0 2.0
>
> 9 1 1.0
>
> 10 1 6.0
>
> 11 1 5.0
>
> 12 1 7.0
>
> 13 1 2.0
>
> 14 1 3.0
>
> 15 1 1.5
>
> 16 1 1.7
>
> 17 2 17.0
>
> 18 2 12.0
>
> 19 2 15.0
>
> 20 2 16.0
>
> 21 2 12.0
>
> 22 2 23.0
>
> 23 2 19.0
>
> 24 2 21.0
>
>
>
> ------------------- creating two dummy variables (C1 and C2) for linear
contrasts-----------------------
>
>>
ar<-data.frame(GROUP=factor(ar$GROUP),C1=factor(ar$GROUP),C2=factor(ar$GROUP),DIP=ar$DIP)
>
>> ar
>
> GROUP C1 C2 DIP
>
> 1 0 0 0 3.0
>
> 2 0 0 0 3.0
>
> 3 0 0 0 4.0
>
> 4 0 0 0 5.0
>
> 5 0 0 0 6.0
>
> 6 0 0 0 7.0
>
> 7 0 0 0 3.0
>
> 8 0 0 0 2.0
>
> 9 1 1 1 1.0
>
> 10 1 1 1 6.0
>
> 11 1 1 1 5.0
>
> 12 1 1 1 7.0
>
> 13 1 1 1 2.0
>
> 14 1 1 1 3.0
>
> 15 1 1 1 1.5
>
> 16 1 1 1 1.7
>
> 17 2 2 2 17.0
>
> 18 2 2 2 12.0
>
> 19 2 2 2 15.0
>
> 20 2 2 2 16.0
>
> 21 2 2 2 12.0
>
> 22 2 2 2 23.0
>
> 23 2 2 2 19.0
>
> 24 2 2 2 21.0
>
>
>
> ------------------- selecting the contrast levels-----------------------
>
>> ar$C1 <- C(ar$C1, c(1,0,-1), how.many = 1)
>
>> ar$C2 <- C(ar$C2, c(1,-1,0), how.many = 1)
>
>
>
>
>
> ------------------- contrast analysis of C2 -----------------------
>
>> r.aov8 <-aov(DIP ~ C2 + GROUP , data = ar)
>
>> anova(r.aov8)
>
> Analysis of Variance Table
>
>
>
> Response: DIP
>
> Df Sum Sq Mean Sq F value Pr(>F)
>
> C2 1 2.10 2.10 0.2622 0.614
>
> GROUP 1 917.00 917.00 114.3460 5.915e-10 ***
>
> Residuals 21 168.41 8.02
>
> ---
>
> Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 `
' 1
>
>
>
> ------------------- contrast analysis of C1 -----------------------
>
>> r.aov9 <-aov(DIP ~ C1 + GROUP , data = ar)
>
>> anova(r.aov9)
>
> Analysis of Variance Table
>
>
>
> Response: DIP
>
> Df Sum Sq Mean Sq F value Pr(>F)
>
> C1 1 650.25 650.25 81.083 1.175e-08 ***
>
> GROUP 1 268.85 268.85 33.525 9.532e-06 ***
>
> Residuals 21 168.41 8.02
>
> ---
>
> Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 `
' 1
>
>
>
> ------------------- anova of the global design -----------------------
>
>> r.aov10 <-aov(DIP ~ GROUP , data = ar)
>
>> anova(r.aov10)
>
> Analysis of Variance Table
>
>
>
> Response: DIP
>
> Df Sum Sq Mean Sq F value Pr(>F)
>
> GROUP 2 919.10 459.55 57.304 3.121e-09 ***
>
> Residuals 21 168.41 8.02
>
> ---
>
> Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 `
' 1
>
>
>
>
>
>
>
>
>
> I would like to know if there is a more economic procedure with R to do
linear contrasts.
>
> Every comments will be well accepted.
>
>
>
> Thank you very much and best regards
>
>
>
> Marco Tommasi
>
> [[alternative HTML version deleted]]
>
Dear Marco
If you are interested in a comparison of a reference level
against each other level (in your case level 1 against level 2
and level 1 against level 3), you can use the summary.lm(),
since this is the default contrast (see ?contr.treatment)
ar <- data.frame(GROUP = factor(rep(1:3, each = 8)),
DIP = c(3.0, 3.0, 4.0, 5.0, 6.0, 7.0, 3.0, 2.0, 1.0, 6.0, 5.0,
7.0, 2.0, 3.0, 1.5, 1.7, 17.0, 12.0, 15.0, 16.0, 12.0, 23.0,
19.0, 21.0))
r.aov10 <- aov(DIP ~ GROUP, data = ar)
anova(r.aov10)
summary.lm(r.aov10)
As result you will get the comparison GROUP 2 against GROUP 1,
denoted by GROUP2 and the comparison GROUP 3 against GROUP 1,
denoted by GROUP3.
Be careful. In your example you include both GROUP and C1 or C2,
respectively in your model. This results in a over parameterized
model and you get a warning that not all coefficients have been
estimated, due to singularities.
It is possible to use other contrasts than contr.treatment,
contr.sum, contr.helmert, contr.poly, but then you have to
specify the correctly in the model.
Regards,
Christoph Buser
--------------------------------------------------------------
Christoph Buser <buser at stat.math.ethz.ch>
Seminar fuer Statistik, LEO C13
ETH (Federal Inst. Technology) 8092 Zurich SWITZERLAND
phone: x-41-44-632-4673 fax: 632-1228
http://stat.ethz.ch/~buser/
--------------------------------------------------------------
Posta Univ. Cagliari writes:
> I have some doubts about the validity of my procedure to estimeate linear
contrasts ina a factorial design.
> For sake of semplicity, let's imagine a one way ANOVA with three
levels. I am interested to test the significance of the difference between the
first and third level (called here contrast C1) and between the first and the
seconda level (called here contrast C2). I used the following procedure:
>
>
> ------------------- reading data from a text file -----------------------
>
> > ar
<-read.table("C:/Programmi/R/myworks/contrasti/cont1.txt",header=TRUE)
>
> > ar
>
> CC GROUP
>
> 1 3.0 0
>
> 2 3.0 0
>
> 3 4.0 0
>
> 4 5.0 0
>
> 5 6.0 0
>
> 6 7.0 0
>
> 7 3.0 0
>
> 8 2.0 0
>
> 9 1.0 1
>
> 10 6.0 1
>
> 11 5.0 1
>
> 12 7.0 1
>
> 13 2.0 1
>
> 14 3.0 1
>
> 15 1.5 1
>
> 16 1.7 1
>
> 17 17.0 2
>
> 18 12.0 2
>
> 19 15.0 2
>
> 20 16.0 2
>
> 21 12.0 2
>
> 22 23.0 2
>
> 23 19.0 2
>
> 24 21.0 2
>
>
>
> ------------------- creating a new array of data-----------------------
>
> > ar<-data.frame(GROUP=factor(ar$GROUP),DIP=ar$CC)
>
> > ar
>
> GROUP DIP
>
> 1 0 3.0
>
> 2 0 3.0
>
> 3 0 4.0
>
> 4 0 5.0
>
> 5 0 6.0
>
> 6 0 7.0
>
> 7 0 3.0
>
> 8 0 2.0
>
> 9 1 1.0
>
> 10 1 6.0
>
> 11 1 5.0
>
> 12 1 7.0
>
> 13 1 2.0
>
> 14 1 3.0
>
> 15 1 1.5
>
> 16 1 1.7
>
> 17 2 17.0
>
> 18 2 12.0
>
> 19 2 15.0
>
> 20 2 16.0
>
> 21 2 12.0
>
> 22 2 23.0
>
> 23 2 19.0
>
> 24 2 21.0
>
>
>
> ------------------- creating two dummy variables (C1 and C2) for linear
contrasts-----------------------
>
> >
ar<-data.frame(GROUP=factor(ar$GROUP),C1=factor(ar$GROUP),C2=factor(ar$GROUP),DIP=ar$DIP)
>
> > ar
>
> GROUP C1 C2 DIP
>
> 1 0 0 0 3.0
>
> 2 0 0 0 3.0
>
> 3 0 0 0 4.0
>
> 4 0 0 0 5.0
>
> 5 0 0 0 6.0
>
> 6 0 0 0 7.0
>
> 7 0 0 0 3.0
>
> 8 0 0 0 2.0
>
> 9 1 1 1 1.0
>
> 10 1 1 1 6.0
>
> 11 1 1 1 5.0
>
> 12 1 1 1 7.0
>
> 13 1 1 1 2.0
>
> 14 1 1 1 3.0
>
> 15 1 1 1 1.5
>
> 16 1 1 1 1.7
>
> 17 2 2 2 17.0
>
> 18 2 2 2 12.0
>
> 19 2 2 2 15.0
>
> 20 2 2 2 16.0
>
> 21 2 2 2 12.0
>
> 22 2 2 2 23.0
>
> 23 2 2 2 19.0
>
> 24 2 2 2 21.0
>
>
>
> ------------------- selecting the contrast levels-----------------------
>
> > ar$C1 <- C(ar$C1, c(1,0,-1), how.many = 1)
>
> > ar$C2 <- C(ar$C2, c(1,-1,0), how.many = 1)
>
>
>
>
>
> ------------------- contrast analysis of C2 -----------------------
>
> > r.aov8 <-aov(DIP ~ C2 + GROUP , data = ar)
>
> > anova(r.aov8)
>
> Analysis of Variance Table
>
>
>
> Response: DIP
>
> Df Sum Sq Mean Sq F value Pr(>F)
>
> C2 1 2.10 2.10 0.2622 0.614
>
> GROUP 1 917.00 917.00 114.3460 5.915e-10 ***
>
> Residuals 21 168.41 8.02
>
> ---
>
> Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 `
' 1
>
>
>
> ------------------- contrast analysis of C1 -----------------------
>
> > r.aov9 <-aov(DIP ~ C1 + GROUP , data = ar)
>
> > anova(r.aov9)
>
> Analysis of Variance Table
>
>
>
> Response: DIP
>
> Df Sum Sq Mean Sq F value Pr(>F)
>
> C1 1 650.25 650.25 81.083 1.175e-08 ***
>
> GROUP 1 268.85 268.85 33.525 9.532e-06 ***
>
> Residuals 21 168.41 8.02
>
> ---
>
> Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 `
' 1
>
>
>
> ------------------- anova of the global design -----------------------
>
> > r.aov10 <-aov(DIP ~ GROUP , data = ar)
>
> > anova(r.aov10)
>
> Analysis of Variance Table
>
>
>
> Response: DIP
>
> Df Sum Sq Mean Sq F value Pr(>F)
>
> GROUP 2 919.10 459.55 57.304 3.121e-09 ***
>
> Residuals 21 168.41 8.02
>
> ---
>
> Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 `
' 1
>
>
>
>
>
>
>
>
>
> I would like to know if there is a more economic procedure with R to do
linear contrasts.
>
> Every comments will be well accepted.
>
>
>
> Thank you very much and best regards
>
>
>
> Marco Tommasi
>
> [[alternative HTML version deleted]]
>
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide!
http://www.R-project.org/posting-guide.html
> Date: Mon, 23 Jan 2006 13:25:33 +0100 > From: Christoph Buser <buser at stat.math.ethz.ch> > Subject: Re: [R] linear contrasts with anova > To: "Posta Univ. Cagliari" <mtommasi at unica.it> > Cc: r-help at stat.math.ethz.ch > Message-ID: <17364.52029.251153.507164 at stat.math.ethz.ch> > Content-Type: text/plain; charset=us-ascii > > Dear Marco > > If you are interested in a comparison of a reference level > against each other level (in your case level 1 against level 2 > and level 1 against level 3), you can use the summary.lm(), > since this is the default contrast (see ?contr.treatment) > > ar <- data.frame(GROUP = factor(rep(1:3, each = 8)), > DIP = c(3.0, 3.0, 4.0, 5.0, 6.0, 7.0, 3.0, 2.0, 1.0, 6.0, 5.0, > 7.0, 2.0, 3.0, 1.5, 1.7, 17.0, 12.0, 15.0, 16.0, 12.0, 23.0, > 19.0, 21.0)) > > > r.aov10 <- aov(DIP ~ GROUP, data = ar) > anova(r.aov10) > summary.lm(r.aov10) > > As result you will get the comparison GROUP 2 against GROUP 1, > denoted by GROUP2 and the comparison GROUP 3 against GROUP 1, > denoted by GROUP3. > > Be careful. In your example you include both GROUP and C1 or C2, > respectively in your model. This results in a over parameterized > model and you get a warning that not all coefficients have been > estimated, due to singularities. > > It is possible to use other contrasts than contr.treatment, > contr.sum, contr.helmert, contr.poly, but then you have to > specify the correctly in the model. > > Regards, > > Christoph Buser > > -------------------------------------------------------------- > Christoph Buser <buser at stat.math.ethz.ch> > Seminar fuer Statistik, LEO C13 > ETH (Federal Inst. Technology) 8092 Zurich SWITZERLAND > phone: x-41-44-632-4673 fax: 632-1228 > http://stat.ethz.ch/~buser/ > --------------------------------------------------------------Dear Marco Try also the the below: # Loading packages library(gplots) library(gmodels) # Testing the nature of dF is.data.frame(dF) is.factor(dF$GROUP) is.numeric(dF$DIP) #Plotting GROUPS win.graph(w = 4, h = 5) plotmeans(DIP ~ GROUP, data = dF, mean.labels = TRUE, digits = 3, col = 'blue', connect = FALSE, ylab = 'DIP', xlab = 'GROUP', pch='') # Contrasts attach(dF) #1 2 3 GROUP cmat = rbind('1 vs. 3' = c( 1, 0, -1,), '1 vs. 2' = c( 1, -1, 0,)) av = aov(DIP ~ GROUP, data = dF, contrasts = list(GROUP = make.contrasts(cmat))) sav = summary(av1, split = list(GROUP=list('1 vs. 3'=1, '1 vs. 2'=2))) sav detach(dF) # Another option attach(dF) #A B C fit.contrast(av, GROUP, c( 1, 0, -1)) # from gmodels fit.contrast(av, GROUP, c( 1, -1, 0)) detach(dF) HTH, []s, -- Jose Claudio Faria Brasil/Bahia/Ilh??us/UESC/DCET Estat??stica Experimental/Prof. Adjunto mails: joseclaudio.faria at terra.com.br jc_faria at uesc.br jc_faria at uol.com.br