[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 22:15:52 CEST 2016
Using a bayesian framework instead of a standard likelihood based
estimation of random effects should make it much easier to calculate
different random effects by sex. Alternatively, you could *probably* figure
it out in a structural equation modeling framework if you were bound and
determined to stay with MLE. I don't think there's a way to do this
properly in nlme short of just splitting the data.
I think you are right that it's probably not detrimental to your particular
case since you know sex is the major driver of differences in the variance
for the random effect and you are removing the distortion by adding it as a
fixed effect.
On Fri, Sep 30, 2016 at 4:00 PM, Brenden, Travis <brenden at anr.msu.edu>
wrote:
> 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>
> 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
>
--
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