I've checked the archives for similar posts to my question below and found
similar questions but never a satisfying answer. Links to these previous
questions are here https://stat.ethz.ch/pipermail/r-sig-mixed-models/
2013q3/020655.html
I'm trying to fit a nonlinear mixed-effects growth model using nlme() to
some length and age data for a fish population. The data are repeated
measures from nearly 400 individuals of both male and female fish. It is
easy enough to fit a nonlinear model incorporating sex as a fixed effect
for each of the model parameters and then include a random effect for
individuals for each of the parameters.
library(FSAdata)
library(nlme)
#load the data and subset to the year that will be analyzed
data(WalleyeML)
wml.01 <- filterD(WalleyeML,Year==2001)
#Define the nonlinear growth function to be estimated
gomplog <- function(T,logLinf,tstar,gstar) {
+ exp(logLinf)*exp(-exp(-gstar*(T-tstar)))
+ }
#Fit a model with sex-specific population average parameter estimates
#and individual-specific deviations assuming a constant distribution
#across sexes for the deviations
fit.01 <- nlme(BC.Len~gomplog(BC.Age,logLinf,tstar,gstar),data=wml.01,
+ fixed=logLinf+tstar+gstar~Sex,
+ random=logLinf+tstar+gstar~1|ID,
+ start=list(fixed=c(6.0,0.0,2.0,0.0,0.5,0.0)),
+ weights=varPower(form=~BC.Age))
logLinf.(Intercept) logLinf.SexM tstar.(Intercept)
tstar.SexM gstar.(Intercept)
6.36279069 -0.14087373 1.92551175
-0.23136899 0.54801947
gstar.SexM
0.05194933 head(ranef(fit.01))
logLinf.(Intercept) tstar.(Intercept) gstar.(Intercept)
2001.2392.F 0.015789108 -0.05055326 -0.036448467
2001.2470.M 0.114106528 0.10101465 -0.110922370
2001.2889.F 0.024976125 -0.17716804 -0.036499377
2001.3138.F 0.040140160 -0.06361609 0.048073911
2001.3139.M -0.006431748 0.31297354 -0.009468046
2001.3147.M 0.039636492 0.48169423 -0.256690199
(the first column above is the identifier for an individual fish where the
last letter indicates if that individual is male (M) or female (F).
However, I'm interested in having the underlying probability distributions
for the random effects be different for males and females. In other words,
I would like to have the random effects for males generated from a
MVN(0,Sigma) distribution and the random effects for females generated
from a MVN(0,Sigma1) distribution where Sigma and Sigma1 are different.
Males and females can have very different growth patterns, with female fish
often getting quite a bit larger than males. Therefore it seems reasonable
that there perhaps would be greater variability in random effects for
females than for males. I can't say this for sure but would like to at
least see. I originally thought the following nlme() call would work.
fit.02 <- nlme(BC.Len~gomplog(BC.Age,logLinf,tstar,gstar),data=wml.01,
+ fixed=logLinf+tstar+gstar~Sex,
+ random=logLinf+tstar+gstar~Sex|ID,
+ start=list(fixed=c(6.0,0.0,2.0,0.0,0.5,0.0)),
+ weights=varPower(form=~BC.Age))
But if you look at the random effect estimates from this fitted model, it
appears that 2 sets of random effects are generated for each parameter and
that they apply to every individual regardless of their sex, which is not
what I am looking for although perhaps I am mis-interpreting this output.
If I had coded this correctly, I would expect to see the .(Intercept)
random effects functionally be 0 for males and the .SexM random effects
functionally be 0 for females.
logLinf.(Intercept) logLinf.SexM tstar.(Intercept)
tstar.SexM gstar.(Intercept)
2001.2392.F 0.01492824 -0.0021143909 -0.05017410
0.03296103 -0.033336005
2001.2470.M 0.13741234 -0.0251483452 0.19731797
-0.10774971 -0.113328192
2001.2889.F 0.02000375 -0.0025014450 -0.19655379
0.08813044 -0.026224542
2001.3138.F 0.04430510 -0.0112679068 -0.05814585
-0.01257797 0.042181315
2001.3139.M -0.01154091 -0.0001006231 0.18061890
0.12099219 -0.003301668
2001.3147.M 0.05160132 -0.0134577925 0.40190910
0.06730332 -0.205363908
gstar.SexM
2001.2392.F -3.205562e-03
2001.2470.M 8.861980e-05
2001.2889.F -9.362266e-05
2001.3138.F 1.386004e-02
2001.3139.M -4.230159e-03
2001.3147.M -5.382077e-02
It doesn't make sense to treat Sex as a random effect because there are
only two levels. If I could treat Sex as a random effect then the
following model call would work with the random effect of ID nested within
the random effect of Sex. But I'm rather interested in having the random
effect of ID nested within the fixed effect of Sex. There shouldn't be
any issue of sample size as the dataset are there are about 188 individual
males (1010 total observations) and 199 individual females (1138 total
observations) in the dataset.
fit.03 <- nlme(BC.Len~gomplog(BC.Age,logLinf,tstar,gstar),data=wml.01,
+ fixed=logLinf+tstar+gstar~1,
+ random=list(Sex=logLinf+tstar+gstar~1,ID=logLinf+tstar+
gstar~1
+ start=list(fixed=c(6.0,2.0,0.5)),
+ weights=varPower(form=~BC.Age))
I've tried other model calls but with no success. Again, I realize
questions similar to this have posted but I have not seen an answer as to
how this could be implement in nlme() if at all. Any suggestions would be
appreciated. Thank you.
Travis Brenden
Department of Fisheries and Wildlife
Michigan State University
[[alternative HTML version deleted]]