
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]]

