From: Tomas Goicoa <tomas.goicoa_at_unavarra.es>

Date: Fri 19 Jan 2007 - 19:32:30 GMT

Date: Fri 19 Jan 2007 - 19:32:30 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 asReceived on Sat Jan 20 06:49:25 2007

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

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 Fri 19 Jan 2007 - 21:30:28 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.
*