Dear Tomas
You can produce the results in Montgomery Montgomery with
lme. Please remark that you should indicate the nesting with the
levels in your nested factor. Therefore I recreated your data,
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)
suppli2 15.056 7.528 2.8526 0.07736 .
suppli:batch 9 69.917 7.769 2.9439 0.01667 *
Residuals24 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
--
Credit and Surety PML study: visit our web page www.cs-pml.org
--
Christoph Buser [EMAIL PROTECTED]
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)
suppli2 15.056 7.528 2.8526 0.07736 .
suppli:batch 9 69.917 7.769 2.9439 0.01667 *
Residuals24 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
Residual2.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.417 0.750 1.583
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
312
From VarCorr I obtain the variance component 1.71, but I am not sure if
