[R] problem with se.contrast()
Christoph Buser
buser at stat.math.ethz.ch
Thu Feb 17 17:31:10 CET 2005
Dear Jamie
As Prof. Ripley explained your analysis is equivalent to the fixed effect
models for the means, so you can calculate it by (if this is your design):
> Lab <- factor(rep(c("1","2","3"),each=12))
> Material <- factor(rep(c("A","B","C","D"),each=3,times=3))
> Measurement <- c(12.20,12.28,12.16,15.51,15.02,15.29,18.14,18.08,18.21,
> 18.54,18.36,18.45,12.59,12.30,12.67,14.98,15.46,15.22,
> 18.54,18.31,18.60,19.21,18.77,18.69,12.72,12.78,12.66,
> 15.33,15.19,15.24,18.00,18.15,17.93,18.88,18.12,18.03)
> testdata <- data.frame(Lab,Material,Measurement)
> rm(list=c("Lab","Material","Measurement"))
You can aggregate your data
> dat.mean <- aggregate(testdata$Measurement,
> by = list(Material=testdata$Material,Lab=testdata$Lab),
> FUN = mean)
> names(dat.mean)[3] <- "Measurement"
> test.red.aov1 <- aov(Measurement ~ Lab + Material, data = dat.mean)
> se.contrast(test.red.aov1,
> list(Material=="A",Material=="B",Material=="C",Material=="D"),
> coef=c(0.5,0.5,-0.5,-0.5),dat.mean)
> [1] 0.1220339
By aggregating the data you bypass the problem in se.contrast and you do
not need R-devel.
-----------------------------------------------------------------------------
The second way to get the same is to set your contrast for the factor
"Material" and calculate you model with this contrast and use summary.lm:
> dat.mean$Material <- C(dat.mean$Material, c(-0.5,-0.5,0.5,0.5),
> how.many = 3)
> test.red.aov2 <- aov(Measurement ~ Lab + Material, data = dat.mean)
> summary.lm(test.red.aov2)
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 16.02000 0.10568 151.583 5.56e-12 ***
Lab2 0.25833 0.14946 1.728 0.135
Lab3 0.06583 0.14946 0.440 0.675
Material1 4.52278 0.12203 37.062 2.58e-08 ***
Material2 1.21056 0.12203 9.920 6.06e-05 ***
Material3 1.55389 0.12203 12.733 1.44e-05 ***
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
Material1 is now the contrast you are interested in and you get beside the
Std. Error the Estimate, too. Material2 and Material3 are just orthogonal
contrasts to complete.
-----------------------------------------------------------------------------
The third way (which I prefer) is using lme from the package nlme
> library(nlme)
Here I use the origianl data set (not aggregated) and set the desired
contrast, too.
> testdata$Material <- C(testdata$Material, c(-0.5,-0.5,0.5,0.5),
> how.many = 3)
> test.lme <- lme(fixed = Measurement ~ Material, data = testdata,
> random = ~1|Lab/Material)
With anova you get the F-Test for the fixed factor "Material"
> anova(test.lme)
> numDF denDF F-value p-value
(Intercept) 1 24 43301.14 <.0001
Material 3 6 544.71 <.0001
and with the summary you have your contrast with standard error:
> summary(test.lme)
Fixed effects: Measurement ~ Material
Value Std.Error DF t-value p-value
(Intercept) 16.128056 0.07750547 24 208.08925 0e+00
Material1 4.522778 0.12203325 6 37.06185 0e+00
Material2 1.210556 0.12203325 6 9.91988 1e-04
Material3 1.553889 0.12203325 6 12.73332 0e+00
-----------------------------------------------------------------------------
Last but not least I tried it with R-devel and the original data frame:
First I reset the contrast on the default value:
testdata$Material <- C(testdata$Material, "contr.treatment", how.many = 3)
I used your syntax and se.contrast():
> test.aov <- with(testdata,aov(Measurement ~ Material + Error(Lab/Material)))
> se.contrast(test.aov,
> list(Material=="A",Material=="B",Material=="C",Material=="D"),
> coef=c(0.5,0.5,-0.5,-0.5),data=testdata)
> [1] 0.1432572
I got a different result and I have admit that I didn't understand why
there is a differnce between the lme model and this one. There are some
comments in the help pages but I'm not sure if this is the answer.
-----------------------------------------------------------------------------
I hope some of the code above can help to analyse your data. Maybe Prof.
Ripley was right and you have another design. Then you just can ignore this
and your life is much more easier :-)
Best regards,
Christoph Buser
--------------------------------------------------------------
Christoph Buser <buser at stat.math.ethz.ch>
Seminar fuer Statistik, LEO C11
ETH (Federal Inst. Technology) 8092 Zurich SWITZERLAND
phone: x-41-1-632-5414 fax: 632-1228
http://stat.ethz.ch/~buser/
--------------------------------------------------------------
Jamie Jarabek writes:
> I am having trouble getting standard errors for contrasts using se.contrast() in
> what appears to be a simple case to me. The following test example illustrates
> my problem:
>
> Lab <- factor(rep(c("1","2","3"),each=12))
> Material <- factor(rep(c("A","B","C","D"),each=3,times=3))
> Measurement <- c(12.20,12.28,12.16,15.51,15.02,15.29,18.14,18.08,18.21,18.54,18.36
> ,18.45,12.59,12.30,12.67,14.98,15.46,15.22,18.54,18.31,18.60,19.21,18.77
> ,18.69,12.72,12.78,12.66,15.33,15.19,15.24,18.00,18.15,17.93,18.88,18.12,18.03)
>
> testdata <- data.frame(Lab,Material,Measurement)
> rm(list=c("Lab","Material","Measurement"))
>
> test.aov <- with(testdata,aov(Measurement ~ Material + Error(Lab/Material)))
>
> This gives me the desired ANOVA table. I next want to get the standard
> errors for certain contrasts and following the help page for
> se.contrast() I tried the following but I get an error:
>
> >se.contrast(test.aov,list(Material=="A",Material=="B",Material=="C",Material=="D"),coef=c(1,1,-1,-1),data=testdata)
> Error in matrix(0, length(asgn), ncol(effects), dimnames = list(nm[1 + :
> length of dimnames [1] not equal to array extent
>
> I have tested this on R 2.0.1 on Windows XP and Solaris and get the same
> error on both systems. I am unsure as to what I am doing wrong here. Thanks for
> any help.
>
> Jamie Jarabek
>
> ______________________________________________
> 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
More information about the R-help
mailing list