[R-sig-ME] How to get estimates of time (Age) as well as Lots in a mixed model

[Jesus Frias]

Dear Robert,

I think that your situation is similar to many other mixed effect 
problems (e.g. a clinical trial in phase II where there are a lot of 
subjects but there is very few information on each individual patient). 

In my personal opinion (and I have used R for something very similar to 
what you are trying to do) the nlme library is well able to solve your 
problem. However you probably will need a lot more data than the 
example you have there to have an estimate of the lot random with  some 
confidence.

In my opinion

fit1<- lme(Conc ~ Age, data=blots, random=~1|Lot)

will do what you want. You have at least two replicates from each lot 
with the same Age. 

Be mindful however that the lot-to-lot random effect  you are 
estimating is only at the intercept. All your lots are behaving with 
the same time slope. If you want to study if the lot is affecting how 
fast the food product degrades you might need a model like:

fit1<- lme(Conc ~ Age, data=blots, random=~Age|Lot)

And that is where your confounding might give you trouble.

Is nice to see that other people in the area of food tech are using 
this fantastic tool!

regards,

Jesus

Robert A. LaBudde wrote:


>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"
>
>_______________________________________________
>R-sig-mixed-models at r-project.org mailing list
>https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>

--------------------------------------------------
Jesús María Frías Celayeta
Head of the Dept of Food Science 
School of Food Sci. and Env. Health 
DIT 
t + 353 1 402 4459 
f + 353 1 402 4495
http://fseh.dit.ie/o4/StaffListing/JesusFrias.html

This message has been scanned for content and viruses by the DIT Information Services E-Mail Scanning Service, and is believed to be clean. http://www.dit.ie