[R-sig-ME] Random effect distributions varying by factor level in an nlme() repeated measures model
Brenden, Travis
brenden at anr.msu.edu
Fri Sep 30 22:00:06 CEST 2016
John, Thank you for the response. That is certainly a fall-back option in this particular case. I guess I’m trying to think about this more broadly as to whether/how this type of model could be implemented in nlme(). It is straightforward to implement in something like jags, but perhaps its structure isn’t amenable to nlme() format? I would be fine with this and could then just move on. I just didn’t want to conclude that without checking with others as to whether there is a solution to this that I am just missing. Alternatively, perhaps there isn’t really a “cost” associated with a common variance-covariance for the random effects for all individuals even if one of the factor levels (in this case female sex) is thought to be more variable as it would “subsume” the variability in the other factor level (male sex).
Travis
From: Poe, John [mailto:jdpo223 at g.uky.edu]
Sent: Friday, September 30, 2016 3:28 PM
To: Brenden, Travis <brenden at anr.msu.edu>
Cc: r-sig-mixed-models at r-project.org
Subject: Re: [R-sig-ME] Random effect distributions varying by factor level in an nlme() repeated measures model
If it's only two groups and sample size isn't an issue then why not just estimate a model for males and one for females then predict their random effects distributions and compare them?
On Fri, Sep 30, 2016 at 10:53 AM, Brenden, Travis <brenden at anr.msu.edu<mailto:brenden at anr.msu.edu>> wrote:
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]]
_______________________________________________
R-sig-mixed-models at r-project.org<mailto:R-sig-mixed-models at r-project.org> mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
--
Thanks,
John
John Poe
Doctoral Candidate
Department of Political Science
Research Methodologist
UK Center for Public Health Services & Systems Research
University of Kentucky
111 Washington Avenue, Room 203a
Lexington, KY 40536
www.johndavidpoe.com<http://www.johndavidpoe.com>
[[alternative HTML version deleted]]
More information about the R-sig-mixed-models
mailing list