[R] aov y lme

[Christoph Buser]

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.