[R] Contrast specified with C() - R vs S-Plus problem
Pascal A. Niklaus
Pascal.Niklaus at unibas.ch
Wed Oct 8 15:59:01 CEST 2003
Thanks for the reply. Below is the solution and the S-Plus and R code
that does the same (for documentation).
>I can't reproduce that in S-PLUS 6.1, and it is not as documented:
>
>
In S-Plus 2000, C() complements the contrast matrix with orthogonal
contrasts if only the first is given.
> CO2 <- factor(rep(c("A","C","E"),each=8))
> m <- rep(seq(0,3,length=8),3) + rep(c(0,1,2),each=8)
> summary.aov(aov(m ~ C(CO2,c(1,0,-1))))
Df Sum of Sq Mean Sq F Value Pr(F)
C(CO2, c(1, 0, -1)) 2 16.00000 8.000000 7.259259 0.004013532
Residuals 21 23.14286 1.102041
> summary.lm(aov(m ~ C(CO2,c(1,0,-1))))
Coefficients:
Value Std. Error t value Pr(>|t|)
(Intercept) 2.5000 0.2143 11.6667 0.0000
C(CO2, c(1, 0, -1))1 -1.0000 0.2624 -3.8103 0.0010
C(CO2, c(1, 0, -1))2 0.0000 0.3712 0.0000 1.0000
Residual standard error: 1.05 on 21 degrees of freedom
Multiple R-Squared: 0.4088
F-statistic: 7.259 on 2 and 21 degrees of freedom, the p-value is 0.004014
Correlation of Coefficients:
(Intercept) C(CO2, c(1, 0, -1))1
C(CO2, c(1, 0, -1))1 0
C(CO2, c(1, 0, -1))2 0 0
In the meanwhile, I found that the same can be done in R like this:
> library(gregmisc)
> CO2 <- factor(rep(c("A","C","E"),each=8))
> m <- rep(seq(0,3,length=8),3) + rep(c(0,1,2),each=8)
> contrastsCO2 <- rbind("A vs E"=c(1,0,-1))
> a <- aov(m ~ CO2, contrasts=list( CO2=make.contrasts(contrastsCO2)))
> summary(a)
Df Sum Sq Mean Sq F value Pr(>F)
CO2 2 16.000 8.000 7.2593 0.004014 **
Residuals 21 23.143 1.102
> summary.lm(a)
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.500e+00 2.143e-01 11.67 1.22e-10 ***
CO2A vs E -2.000e+00 5.249e-01 -3.81 0.00102 **
CO2C2 9.774e-17 3.712e-01 2.63e-16 1.00000
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
Residual standard error: 1.05 on 21 degrees of freedom
Multiple R-Squared: 0.4088, Adjusted R-squared: 0.3525
F-statistic: 7.259 on 2 and 21 DF, p-value: 0.004014
More information about the R-help
mailing list