On Thu, 17 Feb 2005, Christoph Buser wrote:
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.
Note that this is not the same contrast, and you need to double the value for comparability with the original problem.
-----------------------------------------------------------------------------
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.
It is. You used the default `constrasts', which are not actually contrasts. With
options(contrasts=c("contr.helmert", "contr.poly"))it gives the same answer as the other two. Because you used non-contrasts there was an efficiency loss (to the Intercept stratum).
-- Brian D. Ripley, [EMAIL PROTECTED] Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UK Fax: +44 1865 272595
______________________________________________ R-help@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
