[R] problem with se.contrast()
Jamie Jarabek
jjarabek at stat.ufl.edu
Thu Feb 17 18:56:42 CET 2005
Christoph,
Thank you for your advice. My actual design is indeed more complicated than what
I have indicated here. I was just using this as a toy example illustrate my
particular problem. As suggested by Prof. Ripley I will download R-devel and see
if the fixes included within alleviate my problems.
Jamie Jarabek
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.
>
> -----------------------------------------------------------------------------
>
> 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