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]]
