Re: [R] aov y lme

From: Christoph Buser <buser_at_stat.math.ethz.ch>
Date: Mon 22 Jan 2007 - 09:08:54 GMT

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@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@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.

______________________________________________ R-help@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.

Received on Mon Jan 22 20:13:42 2007

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.1.8, at Mon 22 Jan 2007 - 09:30:32 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.