Dear R gurus,
I am trying to work out the problem given in Nested design - Montgomery - Design
of Experiments p.561
I have attached a pdf of the data as well the anova table. It is a mixed model
with Supplier as fixed effect and batches within the supplier as random effects.
I am able to work out the error stratums as below using aov. Which agrees
perfectly with the book example
process=read.table("purity.csv", sep=",", header=TRUE)
proc1=aov(Purity~Supplier+Error(Supplier/Batch), process)
#Results in
print(summary(proc1))
#Error: Supplier -corresponds to fixed effect of Supplier
# Df Sum Sq Mean Sq
#Supplier 2 15.0556 7.5278
#Error: Supplier:Batch -corresponds to random effect of Batch in Supplier
# Df Sum Sq Mean Sq F value Pr(>F)
#Residuals 9 69.917 7.769
#Error: Within --corresponds to error within replications
# Df Sum Sq Mean Sq F value Pr(>F)
#Residuals 24 63.333 2.639
But am not able to get results from the lme model as below. (The model also
seems to be wrong with lot of NaNs produced)
How do I frame the model such that I can get only ~1 | Batch %in% Supplier
Or something else is wrong
Please tell me where I have gone wrong
library(nlme)
proclme=lme(Purity~Supplier,random = ~1|Supplier/Batch,process)
summary(proclme)
VarCorr(proclme)
Linear mixed-effects model fit by REML
Data: process
AIC BIC logLik
154.8440 163.823 -71.42198
Random effects:
Formula: ~1 | Supplier
(Intercept)
StdDev: 1.250514
Formula: ~1 | Batch %in% Supplier
(Intercept) Residual
StdDev: 1.307622 1.624466
Fixed effects: Purity ~ Supplier
Value Std.Error DF t-value p-value
(Intercept) -0.4166667 1.486998 24 -0.2802067 0.7817
SupplierB2 0.7500000 2.102932 0 0.3566449 NaN
SupplierB3 1.5833333 2.102932 0 0.7529169 NaN
Correlation:
(Intr) SpplB2
SupplierB2 -0.707
SupplierB3 -0.707 0.500
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-1.47513431 -0.75010146 0.08125895 0.70609383 1.87201306
Number of Observations: 36
Number of Groups:
Supplier Batch %in% Supplier
3 12
Warning message:
In pt(q, df, lower.tail, log.p) : NaNs produced> VarCorr(proclme)
Variance StdDev
Supplier = pdLogChol(1)
(Intercept) 1.563785 1.250514
Batch = pdLogChol(1)
(Intercept) 1.709877 1.307622
Residual 2.638889 1.624466> intervals(proclme)
Error in intervals.lme(proclme) :
Cannot get confidence intervals on var-cov components: Non-positive definite
approximate variance-covariance
In addition: Warning message:
In qt(p, df, lower.tail, log.p) : NaNs produced
Divaker
Dr. C. Divaker Durairaj, ME, Ph.D
Professor, Farm Machinery
Agricultural Machinery Research Centre
Tamil Nadu Agricultural University
Coimbatore 641003, India
Ph: 91-422-6611204
-------------- next part --------------
A non-text attachment was scrubbed...
Name: montgomry.pdf
Type: application/pdf
Size: 160034 bytes
Desc: not available
Url :
https://stat.ethz.ch/pipermail/r-help/attachments/20071112/ed685b4c/attachment.pdf
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