[R] R - lme

[S Ellison]

Divaker,

Thanks for the data.

For me,
> summary(aov(Purity~Supplier/Batch, process))

gives exactly the same results for mean squares as
> aov(Purity~Supplier+Error(Supplier/Batch), process)
except that the latter gives no p-values (because Supplier appears as
both error term and fixed effect, there isn't anything left to test
supplier against)
For compactness, I'll use the mean squares in

summary(aov(Purity~Supplier/Batch,data=process))
               Df Sum Sq Mean Sq F value  Pr(>F)  
Supplier        2 15.056   7.528  2.8526 0.07736 .
Supplier:Batch  9 69.917   7.769  2.9439 0.01667 *
Residuals      24 63.333   2.639  

The component of variance for batch is
(7.769-2.639)/3 = 1.71(1.307 as SD)

and for supplier, the between-supplier component of variance comes out
negative because it's 7.53-7.77 etc... and if you were testing the
supplier as a fixed effect, testing against the batch MS would imply
that supplier is not a significant effect.


For a mixed-effects model, I normally expected lme's syntax to be

lme(Purity~Supplier, random=~1|Batch, data=P),

but in this case this does not treat the batch as nested, because the
same batch levels appear in all Suppliers. I get the classical result
from 
> process$BS<-factor(process$Supplier:process$Batch)
> lme(Purity~Supplier, random=~1|BS, data=process)
Linear mixed-effects model fit by REML
  Data: P 
  Log-restricted-likelihood: -71.42198
  Fixed: Purity ~ Supplier 
(Intercept)  SupplierT2  SupplierT3 
 -0.4166667   0.7500000   1.5833333 

Random effects:
 Formula: ~1 | BS
        (Intercept) Residual
StdDev:    1.307622 1.624466

Which returns the expected 1.307 between-batch SD. 

Lesson, subject to correction from the gods of lme: lme's random
effects are not treated as nested within fixed effects groups unless the
level identifiers for the random effects are unique to the fixed effects
group they are in. 

And I obviously need to read Pinheiro and Bates again, or I'd have
known that without thinking about it.

Steve E