[R-sig-ME] Random effect distributions varying by factor level in an nlme() repeated measures model
Poe, John
jdpo223 at g.uky.edu
Fri Sep 30 21:28:06 CEST 2016
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>
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 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
[[alternative HTML version deleted]]
More information about the R-sig-mixed-models
mailing list