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)
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
--------------------------------------------------------------
Credit and Surety PML study: visit our web page www.cs-pml.org
--------------------------------------------------------------
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
> PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.