An embedded and charset-unspecified text was scrubbed... Name: not available URL: <https://stat.ethz.ch/pipermail/r-sig-mixed-models/attachments/20100817/e514c91c/attachment.pl>
Modelling heterogeneity and crossed random effects
5 messages · Amelie LESCROEL, ONKELINX, Thierry
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
_______________________________________________
R-sig-mixed-models at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
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.
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.
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.