[R-sig-ME] Random effect distributions varying by factor level in an nlme() repeated measures model

[Brenden, Travis]

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))
> fixef(fit.01)
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))
> 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.
> head(ranef(fit.02))
            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]]