[R] problem with se.contrast()

[Prof Brian Ripley]

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,                  ripley at stats.ox.ac.uk
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