[R-sig-ME] How to get estimates of time (Age) as well as Lots in a mixed model
Robert A. LaBudde
ral at lcfltd.com
Tue Jun 29 18:44:14 CEST 2010
I have a small dataset that represents measurements across Lots of a
processed food, where each lot has a particular Age since manufacturer.
The object of the study is twofold: 1) Determine the Lot-Lot random
effect size, and 2) determine the degradation slope with Age.
Unfortunately, there is only one unique Lot for each particular
unique Age. (A Lot is all manufacture from a particular date.) So the
Lot effect is confounded to some extent with the Age effect.
My question is: What is the best way to get estimate of the two
effects from this dataset? Should I use the Lot effect from a random
effects model with no Age variable (so the Lot effect will include
the Age degradation effect), or should I fit a mixed model anyway,
and let the software sort them out (which it apparently does by
dropping a level of Lots)?
I happen to be using the "nlme" package, so here is the example (I
get a little larger Lot standard deviation when I include Age in the model):
> blots<- read.table('pseudolots.txt', header=TRUE)
> blots$ID<- factor(blots$ID)
> blots$Unit<- factor(blots$Unit)
> blots
ID Lot Unit Age Conc
1 1 A 1 3 1.44
2 2 A 2 3 1.56
3 3 B 3 41 1.03
4 4 B 4 41 1.57
5 5 C 5 229 1.49
6 6 C 6 229 1.66
7 7 D 7 238 0.88
8 8 D 8 238 0.93
9 9 E 9 349 1.43
10 10 E 10 349 1.22
11 11 F 11 391 1.42
12 12 F 12 391 1.46
> #Modeling
> require('nlme')
Loading required package: nlme
> fit1<- lme(Conc ~ Age, data=blots, random=~1|Lot)
> summary(fit1)
Linear mixed-effects model fit by REML
Data: blots
AIC BIC logLik
22.74858 23.95893 -7.374292
Random effects:
Formula: ~1 | Lot
(Intercept) Residual
StdDev: 0.2320592 0.1786757
Fixed effects: Conc ~ Age
Value Std.Error DF t-value p-value
(Intercept) 1.3710359 0.18970706 6 7.227121 0.0004
Age -0.0001449 0.00074846 4 -0.193538 0.8560
Correlation:
(Intr)
Age -0.823
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-1.59441923 -0.45471612 -0.06071212 0.52440223 1.42781630
Number of Observations: 12
Number of Groups: 6
> anova(fit1)
numDF denDF F-value p-value
(Intercept) 1 6 154.51072 <.0001
Age 1 4 0.03746 0.856
>
> fit2<- lme(Conc ~ 1, data=blots, random=~1|Lot)
> summary(fit2)
Linear mixed-effects model fit by REML
Data: blots
AIC BIC logLik
8.122374 9.31606 -1.061187
Random effects:
Formula: ~1 | Lot
(Intercept) Residual
StdDev: 0.2010265 0.1786757
Fixed effects: Conc ~ 1
Value Std.Error DF t-value p-value
(Intercept) 1.340833 0.09693139 6 13.83281 0
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-1.57582754 -0.56625710 -0.01917496 0.56893404 1.44640790
Number of Observations: 12
Number of Groups: 6
> anova(fit2)
numDF denDF F-value p-value
(Intercept) 1 6 191.3466 <.0001
================================================================
Robert A. LaBudde, PhD, PAS, Dpl. ACAFS e-mail: ral at lcfltd.com
Least Cost Formulations, Ltd. URL: http://lcfltd.com/
824 Timberlake Drive Tel: 757-467-0954
Virginia Beach, VA 23464-3239 Fax: 757-467-2947
"Vere scire est per causas scire"
More information about the R-sig-mixed-models
mailing list