# [R] aov y lme

Christoph Buser buser at stat.math.ethz.ch
Mon Jan 22 10:08:54 CET 2007

```Dear Tomas

You can produce the results in Montgomery Montgomery with
lme. Please remark that you should indicate the nesting with the
but used 1,...,12 for the levels of batch instead of 1,...,4.

purity<-c(1,-2,-2,1,-1,-3,0,4, 0,-4, 1,0, 1,0,-1,0,-2,4,0,3,
-3,2,-2,2,2,-2,1,3,4,0,-1,2,0,2,2,1)
suppli<-factor(c(rep(1,12),rep(2,12),rep(3,12)))
batch<-factor(c(rep(1:4,3), rep(5:8,3), rep(9:12,3)))
material<-data.frame(purity,suppli,batch)

As you remarked you can use aov

summary(material.aov<-aov(purity~suppli+suppli:batch,data=material))
Df Sum Sq Mean Sq F value  Pr(>F)
suppli        2 15.056   7.528  2.8526 0.07736 .
suppli:batch  9 69.917   7.769  2.9439 0.01667 *
Residuals    24 63.333   2.639
---
Signif. codes:  0 \$,1rx(B***\$,1ry(B 0.001 \$,1rx(B**\$,1ry(B 0.01 \$,1rx(B*\$,1ry(B 0.05 \$,1rx(B.\$,1ry(B 0.1 \$,1rx(B \$,1ry(B 1

Remark that the test of "suppli" is not the same as in
Montgomery. Here it is wrongly tested against the Residuals and
you should perform the calculate the test with:

(Fsuppi <- summary(material.aov)[[1]][1,"Mean Sq"]/
summary(material.aov)[[1]][2,"Mean Sq"])
pf(Fsuppi, df1 = 2, df2 = 9)

To use lme the correct level numbering is now important to
indicate the nesting. You should specify your random component
as

random=~1|batch

If you use "random=~1|suppli/batch" instead, random components
for batch AND suppli are included in the model and you specify a
model that incorporates suppli as random and fixed. Therefore
the correct syntax is

library(nlme)
material.lme<-lme(purity~suppli,random=~1|batch,data=material)
## You obtain the F-test for suppli using "anova"
anova(material.lme)
summary(material.lme)
## Remark that in the summary output, the random effects are
## standard deviations and not variance components and you
## should square them to compare them with Montgomery
## 1.307622^2 = 1.71             1.624466^2 = 2.64

## Or you can use
VarCorr(material.lme)

I hope this helps you.

Best regards,

Christoph Buser

--------------------------------------------------------------

--------------------------------------------------------------
Christoph Buser <buser at stat.math.ethz.ch>
Seminar fuer Statistik, LEO C13
ETH Zurich	8092 Zurich	 SWITZERLAND
phone: x-41-44-632-4673		fax: 632-1228
http://stat.ethz.ch/~buser/
--------------------------------------------------------------

Tomas Goicoa writes:
>
>
>
> Dear R user,
>
> I am trying to reproduce the results in Montgomery D.C (2001, chap 13,
> example 13-1).
>
> Briefly, there are three suppliers, four batches nested within suppliers
> and three determinations of purity (response variable) on each batch. It is
> a two stage nested design, where suppliers are fixed and batches are random.
>
> y_ijk=mu+tau_i+beta_j(nested in tau_i)+epsilon_ijk
>
> Here are the data,
>
> purity<-c(1,-2,-2,1,
>           -1,-3, 0,4,
>            0,-4, 1, 0,
>            1,0,-1,0,
>            -2,4,0,3,
>            -3,2,-2,2,
>            2,-2,1,3,
>            4,0,-1,2,
>            0,2,2,1)
>
> suppli<-factor(c(rep(1,12),rep(2,12),rep(3,12)))
> batch<-factor(rep(c(1,2,3,4),9))
>
> material<-data.frame(purity,suppli,batch)
>
> If I use the function aov, I get
>
> material.aov<-aov(purity~suppli+suppli:batch,data=material)
> summary(material.aov)
>               Df Sum Sq Mean Sq F value  Pr(>F)
> suppli        2 15.056   7.528  2.8526 0.07736 .
> suppli:batch  9 69.917   7.769  2.9439 0.01667 *
> Residuals    24 63.333   2.639
> ---
> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>
> and I can estimate the variance component for the batches as
>
> (7.769- 2.639)/3=1.71
>
> which is the way it is done in Montgomery, D.
>
> I want to use the function lme because I would like to make a diagnosis of
> the model, and I think it is more appropriate.
>
> Looking at Pinheiro and Bates, I have tried the following,
>
> library(nlme)
> material.lme<-lme(purity~suppli,random=~1|suppli/batch,data=material)
> VarCorr(material.lme)
>
>              Variance     StdDev
> suppli =    pdLogChol(1)
> (Intercept) 1.563785     1.250514
> batch =     pdLogChol(1)
> (Intercept) 1.709877     1.307622
> Residual    2.638889     1.624466
>
> material.lme
>
> Linear mixed-effects model fit by REML
>    Data: material
>    Log-restricted-likelihood: -71.42198
>    Fixed: purity ~ suppli
> (Intercept)     suppli2     suppli3
>   -0.4166667   0.7500000   1.5833333
>
> Random effects:
>   Formula: ~1 | suppli
>          (Intercept)
> StdDev:    1.250514
>
>   Formula: ~1 | batch %in% suppli
>          (Intercept) Residual
> StdDev:    1.307622 1.624466
>
> Number of Observations: 36
> Number of Groups:
>             suppli batch %in% suppli
>                  3                12
>
>  From VarCorr I obtain the variance component 1.71, but I am not sure if
> this is the way to fit the model for the nested design. Here, I also have a
> variance component for suppli and this is a fixed factor. Can anyone give
> me a clue?
> 	[[alternative HTML version deleted]]
>
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help