[R] aov y lme

From: Tomas Goicoa <tomas.goicoa_at_unavarra.es>
Date: Sat 20 Jan 2007 - 11:42:39 GMT

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.
Received on Sat Jan 20 22:46:45 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.