# [R] aov y lme

From: Tomas Goicoa <tomas.goicoa_at_unavarra.es>
Date: Sat 20 Jan 2007 - 11:42:39 GMT

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@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help