[R-sig-ME] Modelling heterogeneity and crossed random effects
ONKELINX, Thierry
Thierry.ONKELINX at inbo.be
Wed Aug 18 11:17:21 CEST 2010
Dear Amelie,
In my opinion, a correlation structure (e.g. corAR1(~Year) or corExp(~Year)) will do to represent your design. And it will give you information about the difference in variance in a year and among years.
A second option would be to add year as a random slope per individual. Random = ~ factor(Year) - 1|ID
You could even combine both options.
Note that according to Zuur et al. (2009) is random intercept is equivalent to a compound symmetry correlation structure.
lme(Z ~ ..., random = ~ 1|A) is equivalent to gls(Z ~ ..., correlation = corCompSymm(~A))
HTH,
Thierry
----------------------------------------------------------------------------
ir. Thierry Onkelinx
Instituut voor natuur- en bosonderzoek
team Biometrie & Kwaliteitszorg
Gaverstraat 4
9500 Geraardsbergen
Belgium
Research Institute for Nature and Forest
team Biometrics & Quality Assurance
Gaverstraat 4
9500 Geraardsbergen
Belgium
tel. + 32 54/436 185
Thierry.Onkelinx at inbo.be
www.inbo.be
To call in the statistician after the experiment is done may be no more than asking him to perform a post-mortem examination: he may be able to say what the experiment died of.
~ Sir Ronald Aylmer Fisher
The plural of anecdote is not data.
~ Roger Brinner
The combination of some data and an aching desire for an answer does not ensure that a reasonable answer can be extracted from a given body of data.
~ John Tukey
> -----Oorspronkelijk bericht-----
> Van: Amelie Lescroel [mailto:amelie.lescroel at univ-rennes1.fr]
> Verzonden: woensdag 18 augustus 2010 10:41
> Aan: ONKELINX, Thierry; r-sig-mixed-models at r-project.org
> Onderwerp: RE: [R-sig-ME] Modelling heterogeneity and crossed
> random effects
>
> Dear Thierry,
>
> Thanks a lot for your answer. I was hoping that year as a
> random effect would 1) account for the study design (I have
> several points per individual for each year and I wanted to
> quantify the correlation of 2 observations from the same
> individual within a year vs. across years) and 2) capture
> other year effects that would not be accounted for by my
> fixed effects. And indeed, all my models including year as a
> random effect performed better, in terms of AIC, than those
> that did not include year. Otherwise, yes, it would easier to
> model the variance in nlme. In either package though, I'm not
> sure that I found the right structure model that would
> correspond to the study design (longitudinal study with
> replicated points within years) and I would welcome any suggestion.
>
> Best,
>
> Amelie
>
> -----Original Message-----
> From: ONKELINX, Thierry [mailto:Thierry.ONKELINX at inbo.be]
> Sent: Wednesday, August 18, 2010 10:30 AM
> To: Amelie Lescroel; r-sig-mixed-models at r-project.org
> Subject: RE: [R-sig-ME] Modelling heterogeneity and crossed
> random effects
>
> Dear Amelie,
>
> Do you expect a common effect of year on all individuals that
> is not captured by your fixed effects? If not, you do not
> need to add year as a random effect and only a random effect
> of individual will do. Hence you could switch back to nlme
> which has more features in terms of variance and correlation
> structures.
>
> HTH,
>
> Thierry
>
> --------------------------------------------------------------
> --------------
> ir. Thierry Onkelinx
> Instituut voor natuur- en bosonderzoek
> team Biometrie & Kwaliteitszorg
> Gaverstraat 4
> 9500 Geraardsbergen
> Belgium
>
> Research Institute for Nature and Forest team Biometrics &
> Quality Assurance Gaverstraat 4 9500 Geraardsbergen Belgium
>
> tel. + 32 54/436 185
> Thierry.Onkelinx at inbo.be
> www.inbo.be
>
> To call in the statistician after the experiment is done may
> be no more than asking him to perform a post-mortem
> examination: he may be able to say what the experiment died of.
> ~ Sir Ronald Aylmer Fisher
>
> The plural of anecdote is not data.
> ~ Roger Brinner
>
> The combination of some data and an aching desire for an
> answer does not ensure that a reasonable answer can be
> extracted from a given body of data.
> ~ John Tukey
>
>
> > -----Oorspronkelijk bericht-----
> > Van: r-sig-mixed-models-bounces at r-project.org
> > [mailto:r-sig-mixed-models-bounces at r-project.org] Namens Amelie
> > Lescroel
> > Verzonden: woensdag 18 augustus 2010 10:06
> > Aan: r-sig-mixed-models at r-project.org
> > Onderwerp: Re: [R-sig-ME] Modelling heterogeneity and
> crossed random
> > effects
> >
> > Dear all,
> > I did not receive any answer to my questions below. Not that I
> > consider that anybody "owes" me an answer but I would really need
> > advices from people more knowledgeable than I am.
> > Please let me know if I need to reformulate / shorten my
> questions or
> > examples or if they are too "naïve".
> > Best regards,
> > Amelie
> >
> > -----Original Message-----
> > From: r-sig-mixed-models-bounces at r-project.org
> > [mailto:r-sig-mixed-models-bounces at r-project.org] On Behalf
> Of Amelie
> > Lescroel
> > Sent: Tuesday, August 17, 2010 10:16 PM
> > To: r-sig-mixed-models at r-project.org
> > Subject: [R-sig-ME] Modelling heterogeneity and crossed
> random effects
> >
> > Dear all,
> >
> >
> >
> > I am currently trying to model the behavioural response of
> individual
> > seabirds (in terms of foraging efficiency) to the variation
> in sea ice
> > cover
> > (SICdr) of their foraging environment. I have 13 years of
> data, birds
> > are individually marked and followed, I have several records (=
> > foraging efficiency data = CPUEr in my
> > code) per individual (IDr) for each year
> > (YEARr) and individuals are followed across years.
> >
> >
> >
> > I am trying to find the right random effect structure (biologically
> > meaningful and dealing with problems of
> > independence) and to deal with heterogeneity of the
> residual variance
> > at the same time (for all my models, the variance of the residuals
> > increases with increasing fitted values).
> > Regarding the random effect structure, would you say that crossed
> > random effects of the form (1|IDr) + (1|YEARr) would
> correctly reflect
> > the study design? Is there any way to model the variance
> heterogeneity
> > in lmer that would be analogous to the varIdent or varFixed
> functions
> > in nlme? So far, I can model the variance heterogeneity
> with nlme only
> > and the
> > (hopefully) appropriate random effect structure with lmer
> only. Would
> > you have other suggestions for dealing with this heteroscedasticity?
> >
> >
> >
> > Here are a couple of examples regarding the random effect structure
> > with some associated questions:
> >
> >
> >
> > > M1 <- lmer(CPUEr~SEXr+SICdr+(1|IDr))
> >
> > > summary(M1)
> >
> >
> >
> > Linear mixed model fit by REML
> >
> > Formula: CPUEr ~ SEXr + SICdr + (1 | IDr)
> >
> > AIC BIC logLik deviance REMLdev
> >
> > 270.2 297.6 -130.1 234.5 260.2
> >
> > Random effects:
> >
> > Groups Name Variance Std.Dev.
> >
> > IDr (Intercept) 0.010906 0.10443
> >
> > Residual 0.060610 0.24619
> >
> > Number of obs: 1759, groups: IDr, 229
> >
> >
> >
> > Fixed effects:
> >
> > Estimate Std. Error t value
> >
> > (Intercept) 0.3070164 0.0155734 19.714
> >
> > SEXrM 0.0961795 0.0195420 4.922
> >
> > SICdr 0.0026240 0.0008478 3.095
> >
> >
> >
> > Correlation of Fixed Effects:
> >
> > (Intr) SEXrM
> >
> > SEXrM -0.612
> >
> > SICdr -0.478 -0.006
> >
> >
> >
> > Here, the correlation between 2 observations from the same
> individual
> > (irrespective of year) is:
> > 0.010906/(0.010906+0.060610)=0.15
> >
> >
> >
> > > M2 <- lmer(CPUEr~SEXr+SICdr+(1|YEARr))
> >
> > > summary(M2)
> >
> > Linear mixed model fit by REML
> >
> > Formula: CPUEr ~ SEXr + SICdr + (1 | YEARr)
> >
> > AIC BIC logLik deviance REMLdev
> >
> > 117.1 144.5 -53.55 84.8 107.1
> >
> > Random effects:
> >
> > Groups Name Variance Std.Dev.
> >
> > YEARr (Intercept) 0.020395 0.14281
> >
> > Residual 0.059892 0.24473
> >
> > Number of obs: 1759, groups: YEARr, 13
> >
> >
> >
> > Fixed effects:
> >
> > Estimate Std. Error t value
> >
> > (Intercept) 0.36443 0.04367 8.345
> >
> > SEXrM 0.10819 0.01175 9.207
> >
> > SICdr -0.00920 0.00192 -4.793
> >
> >
> >
> > Correlation of Fixed Effects:
> >
> > (Intr) SEXrM
> >
> > SEXrM -0.134
> >
> > SICdr -0.367 0.009
> >
> >
> >
> > Here, the correlation between 2 observations from the same year
> > (irrespective of the bird) is:
> > 0.020395/(0.020395+0.059892)=0.25 How do I get the correlation of 2
> > observations from the same individual within a year? By modeling
> > CPUEr~SEXr+SICdr+(1|YEARr/IDr)?
> >
> >
> >
> > > M3 <- lmer(CPUEr~SEXr+SICdr+(1|YEARr/IDr))
> >
> > > summary(M3)
> >
> > Linear mixed model fit by REML
> >
> > Formula: CPUEr ~ SEXr + SICdr + (1 | YEARr/IDr)
> >
> > AIC BIC logLik deviance REMLdev
> >
> > 51.29 84.12 -19.64 17.21 39.29
> >
> > Random effects:
> >
> > Groups Name Variance Std.Dev.
> >
> > IDr:YEARr (Intercept) 0.0097178 0.09858
> >
> > YEARr (Intercept) 0.0188065 0.13714
> >
> > Residual 0.0500727 0.22377
> >
> > Number of obs: 1759, groups: IDr:YEARr, 543; YEARr, 13
> >
> >
> >
> > Fixed effects:
> >
> > Estimate Std. Error t value
> >
> > (Intercept) 0.357318 0.042408 8.426
> >
> > SEXrM 0.104650 0.014207 7.366
> >
> > SICdr -0.008960 0.001855 -4.831
> >
> >
> >
> > Correlation of Fixed Effects:
> >
> > (Intr) SEXrM
> >
> > SEXrM -0.166
> >
> > SICdr -0.365 0.004
> >
> >
> >
> > Then, would the correlation of 2 observations from the same
> individual
> > within a year be 0.0097178/(0.0097178+0.0500727)=0.16?
> >
> >
> >
> > My best model (in terms of AIC) so far is the following:
> >
> >
> >
> > > M4 <- lmer(CPUEr~SEXr+SICdr+(SICdr|IDr)+(1|YEARr))
> >
> > > summary(M4)
> >
> > Linear mixed model fit by REML
> >
> > Formula: CPUEr ~ SEXr + SICdr + (SICdr | IDr) + (1 | YEARr)
> >
> > AIC BIC logLik deviance REMLdev
> >
> > 12.88 56.66 1.559 -24.55 -3.119
> >
> > Random effects:
> >
> > Groups Name Variance Std.Dev. Corr
> >
> > IDr (Intercept) 8.9314e-03 0.0945058
> >
> > SICdr 2.3781e-05 0.0048766 -0.464
> >
> > YEARr (Intercept) 2.1401e-02 0.1462922
> >
> > Residual 5.0765e-02 0.2253112
> >
> > Number of obs: 1759, groups: IDr, 229; YEARr, 13
> >
> >
> >
> > Fixed effects:
> >
> > Estimate Std. Error t value
> >
> > (Intercept) 0.363366 0.045471 7.991
> >
> > SEXrM 0.100215 0.017188 5.830
> >
> > SICdr -0.009910 0.001974 -5.021
> >
> >
> >
> > Correlation of Fixed Effects:
> >
> > (Intr) SEXrM
> >
> > SEXrM -0.189
> >
> > SICdr -0.357 0.010
> >
> >
> >
> > How should I interpret the random effects?
> >
> >
> >
> > I am using the R package version 0.999375-31 of lme4 and R version
> > 2.9.2.
> >
> >
> >
> > Thanks in advance for your help!
> >
> >
> >
> > Cheers,
> >
> >
> >
> > Amelie
> >
> >
> >
> >
> >
> >
> >
> >
> >
> >
> > [[alternative HTML version deleted]]
> >
> > _______________________________________________
> > R-sig-mixed-models at r-project.org mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> >
> > _______________________________________________
> > R-sig-mixed-models at r-project.org mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> >
>
> Druk dit bericht a.u.b. niet onnodig af.
> Please do not print this message unnecessarily.
>
> Dit bericht en eventuele bijlagen geven enkel de visie van de
> schrijver weer
>
> en binden het INBO onder geen enkel beding, zolang dit
> bericht niet bevestigd is door een geldig ondertekend
> document. The views expressed in this message and any annex
> are purely those of the writer and may not be regarded as
> stating an official position of INBO, as long as the message
> is not confirmed by a duly signed document.
>
>
Druk dit bericht a.u.b. niet onnodig af.
Please do not print this message unnecessarily.
Dit bericht en eventuele bijlagen geven enkel de visie van de schrijver weer
en binden het INBO onder geen enkel beding, zolang dit bericht niet bevestigd is
door een geldig ondertekend document. The views expressed in this message
and any annex are purely those of the writer and may not be regarded as stating
an official position of INBO, as long as the message is not confirmed by a duly
signed document.
More information about the R-sig-mixed-models
mailing list