[R-sig-ME] correlated random effects in lmer and false convergence
Douglas Bates
bates at stat.wisc.edu
Tue Jun 22 00:04:08 CEST 2010
A quick response as I must catch a bus soon.
On Mon, Jun 21, 2010 at 4:18 PM, Ailene Kane <ailene at u.washington.edu> wrote:
> Dear Dr. Bates and/or other expert lmer users:
>
> We are having problems fitting a linear mixed effects model using the lmer()
> function in the lme4 package, and are hoping that you will be able (and
> willing!) to answer some of our questions. We have read several chapters in
> your book on lmer, and other online resources which have helped resolve many
> of our questions, but our data set is pretty complex and we wanted to make
> very sure that we are interpreting the output (and some error messages that
> we get) correctly. Thank you so much for your time!
>
> As background, we are exploring the influence of climative factors (e.g.
> snow, temperature, etc) on annual tree growth (“rwi”= ring width index, a
> detrended index for tree growth that ranges from ~-1 to +2). We are
> evaluating the fit of 29 different models, including various climate
> variables and combinations of climate variables as fixed effects (e.g. in
> the model below, GST=Growing season Temperature, and GPT=growing season
> precipitation). Note: the climate variables have been standardized by
> substracting the mean and dividing by the standard deviation. There are 2
> random effects: “ind” (the individual tree; 20 trees were cored at each
> location and they differ in sensitivity to climate) and “yr” (the year of
> tree growth - this is sort of like a block effect, because you could have
> years with similar temperature values). A subset of our dataset is attached
> as a .csv, and some sample code and output is listed below (and in the
> attached file).
>
> We have the following questions:
>
> 1) We originally thought that test13 and test13b were the same thing:
>
> test13<-lmer(rwi~GST*GPT+(0+GST*GPT|ind)+(1|yrs),control=list(maxIter=5000))
> test13b<-lmer(rwi~GST*GPT+(0+GST|ind)+(0+GPT|ind)+(0+GST:GPT|ind)+(1|yrs),control=list(maxIter=5000))
No, they are not the same as you have seen in your later
investigation. The first model generates a vector-valued random
effect for each level of ind with possible correlation in the
variance-covariance matrix. In the second model the components of the
random effects are independent.
It seems to me that there is very little variation associated with the
levels of ind. I enclose a transcript of fitting a series of models
ending up with one that doesn't have any random effects for ind and it
seems to compare well to the others.
I must go to the bus. I'll explain about the model-building strategy later.
> However, we recently noticed that they result in differences in the
> correlations between random effects, and different AIC values (output is
> below and can be generated by the code below, if you are interested). In
> test13, the random effects on ind for GST,GPT, and their interaction are
> perfectly correlated (1.00). In test13b the random effects are NOT
> PERFECTLY correlated (r= 0.6349112). Does this mean that test13 ASSUMES
> correlated random effects (because of the way that it is coded)? We
> definitely don't want this and wanted to make sure that we understood
> correctly how to code the model appropriately. We are comparing the 29
> potential models to a null model and using the delta AIC (in part) to
> evaluate the best-fit model. The null model is coded as follows:
>
> test0<-lmer(rwi~1+(0+1|ind)+(1|yrs),control=list(maxIter=5000)
>
> 2) We have been getting warning message for some models (for example,
> test13): “Warning message: In mer_finalize(ans) : singular convergence (7).”
> We have explored our data and get this error message only for complex
> models (e.g. 2 climate variables with an interaction), and mainly at
> locations where tree growth (of individual trees) is not correlated to any
> of our climate variables. Do you think this interpretation makes sense, in
> which case we can conclude that the model fits poorly, indicating that other
> factors besides climate influence growth? Or, is there something else that
> you think we should do to avoid singular convergence? Again, the output
> generated by the code (test13b) is below
>
> and attached, in case you are interested.
> Thank you so much for your help!
>
> Best wishes,
>
> Ailene
>
> R OUTPUT FROM ATTACHED CODE:
>> summary(test13)
> Linear mixed model fit by REML
> Formula: rwi ~ GST * GPT + (0 + GST * GPT | ind) + (1 | yrs)
> AIC BIC logLik deviance REMLdev
> -528.4 -463 276.2 -577 -552.4
> Random effects:
> Groups Name Variance Std.Dev. Corr
> yrs (Intercept) 0.02764111 0.166256
> ind GST 0.00033501 0.018303
> GPT 0.00029068 0.017049 1.000
> GST:GPT 0.00023243 0.015246 1.000 1.000
> Residual 0.03587474 0.189406
> Number of obs: 1727, groups: yrs, 94; ind, 20
> Fixed effects:
> Estimate Std. Error t value
> (Intercept) 0.999112 0.018777 53.21
> GST 0.071152 0.019373 3.67
> GPT 0.036372 0.019266 1.89
> GST:GPT -0.004302 0.019494 -0.22
> Correlation of Fixed Effects:
> (Intr) GST GPT
> GST 0.033
> GPT -0.024 0.340
> GST:GPT 0.320 0.140 -0.032
>> summary(test13b)
> Linear mixed model fit by REML
> Formula: rwi ~ GST * GPT + (0 + GST | ind) + (0 + GPT | ind) + (0 + GST:GPT
> | ind) + (1 | yrs)
> AIC BIC logLik deviance REMLdev
> -522.7 -473.6 270.4 -565.3 -540.7
> Random effects:
> Groups Name Variance Std.Dev.
> yrs (Intercept) 2.7626e-02 0.1662108
> ind GST:GPT 8.7606e-05 0.0093598
> ind GPT 3.2627e-05 0.0057120
> ind GST 2.3199e-05 0.0048165
> Residual 3.6397e-02 0.1907813
> Number of obs: 1727, groups: yrs, 94; ind, 20
> Fixed effects:
> Estimate Std. Error t value
> (Intercept) 0.999101 0.018778 53.21
> GST 0.071329 0.018967 3.76
> GPT 0.036529 0.018928 1.93
> GST:GPT -0.004107 0.019310 -0.21
> Correlation of Fixed Effects:
> (Intr) GST GPT
> GST 0.034
> GPT -0.024 0.310
> GST:GPT 0.323 0.106 -0.069
> Ailene Kane Ettinger
> PhD Candidate
> Biology Department
> University of Washington
> Box 351800
> Seattle, Washington 98195-1800
> ailene at u.washington.edu
>
>
>
>
-------------- next part --------------
R version 2.11.1 (2010-05-31)
Copyright (C) 2010 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
Natural language support but running in an English locale
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> library(lme4a)
Loading required package: Matrix
Loading required package: lattice
Attaching package: 'Matrix'
The following object(s) are masked from 'package:base':
det
Loading required package: minqa
Loading required package: Rcpp
Attaching package: 'lme4a'
The following object(s) are masked from 'package:stats':
AIC
> hdat <- within(read.csv("~/tests/abamparadat.csv", header=TRUE),
+ {
+ ind <- factor(inds)
+ yrs <- factor(yrs)
+ })
> xyplot(rwi ~ SWE | inds, hdat, pch=16, type =c("p","r"),
+ main=paste("Abam-snow"), lwd=2, ylim=c(0.5,1.5))
>
> (test13b <- lmer(rwi ~ GST*GPT+(1|inds)+(0+GST|inds)+(0+GPT|inds)+(0+GST:GPT|inds)+(1|yrs),
+ hdat, REML=FALSE))
Linear mixed model fit by maximum likelihood ['merMod']
Formula: rwi ~ GST * GPT + (1 | inds) + (0 + GST | inds) + (0 + GPT | inds) + (0 + GST:GPT | inds) + (1 | yrs)
Data: hdat
AIC BIC logLik deviance
-545.3849 -490.8435 282.6925 -565.3849
Random effects:
Groups Name Variance Std.Dev.
yrs (Intercept) 2.637e-02 0.162397
inds GST:GPT 8.551e-05 0.009247
inds GPT 3.040e-05 0.005514
inds GST 2.080e-05 0.004560
inds (Intercept) 0.000e+00 0.000000
Residual 3.640e-02 0.190796
Number of obs: 1727, groups: yrs, 94; inds, 20
Fixed effects:
Estimate Std. Error t value
(Intercept) 0.999094 0.018377 54.37
GST 0.071329 0.018559 3.84
GPT 0.036532 0.018522 1.97
GST:GPT -0.004108 0.018898 -0.22
Correlation of Fixed Effects:
(Intr) GST GPT
GST 0.034
GPT -0.025 0.310
GST:GPT 0.323 0.106 -0.069
> (test13c <- lmer(rwi ~ GST*GPT+(0+GST|inds)+(0+GPT|inds)+(0+GST:GPT|inds)+(1|yrs),
+ hdat, REML=FALSE))
Linear mixed model fit by maximum likelihood ['merMod']
Formula: rwi ~ GST * GPT + (0 + GST | inds) + (0 + GPT | inds) + (0 + GST:GPT | inds) + (1 | yrs)
Data: hdat
AIC BIC logLik deviance
-547.3849 -498.2977 282.6925 -565.3849
Random effects:
Groups Name Variance Std.Dev.
yrs (Intercept) 2.637e-02 0.162397
inds GST:GPT 8.551e-05 0.009247
inds GPT 3.040e-05 0.005514
inds GST 2.080e-05 0.004560
Residual 3.640e-02 0.190796
Number of obs: 1727, groups: yrs, 94; inds, 20
Fixed effects:
Estimate Std. Error t value
(Intercept) 0.999094 0.018377 54.37
GST 0.071329 0.018559 3.84
GPT 0.036532 0.018522 1.97
GST:GPT -0.004108 0.018898 -0.22
Correlation of Fixed Effects:
(Intr) GST GPT
GST 0.034
GPT -0.025 0.310
GST:GPT 0.323 0.106 -0.069
> (test13d <- lmer(rwi ~ GST*GPT+(0+GST|inds)+(0+GPT|inds)+(1|yrs),
+ hdat, REML=FALSE))
Linear mixed model fit by maximum likelihood ['merMod']
Formula: rwi ~ GST * GPT + (0 + GST | inds) + (0 + GPT | inds) + (1 | yrs)
Data: hdat
AIC BIC logLik deviance
-549.0160 -505.3829 282.5080 -565.0160
Random effects:
Groups Name Variance Std.Dev.
yrs (Intercept) 2.637e-02 0.162376
inds GPT 3.646e-05 0.006038
inds GST 1.146e-05 0.003385
Residual 3.649e-02 0.191035
Number of obs: 1727, groups: yrs, 94; inds, 20
Fixed effects:
Estimate Std. Error t value
(Intercept) 0.999095 0.018377 54.37
GST 0.071351 0.018545 3.85
GPT 0.036539 0.018530 1.97
GST:GPT -0.004121 0.018782 -0.22
Correlation of Fixed Effects:
(Intr) GST GPT
GST 0.034
GPT -0.025 0.310
GST:GPT 0.325 0.107 -0.069
> (test13e <- lmer(rwi ~ GST + GPT+(0+GST|inds)+(0+GPT|inds)+(1|yrs),
+ hdat, REML=FALSE))
Linear mixed model fit by maximum likelihood ['merMod']
Formula: rwi ~ GST + GPT + (0 + GST | inds) + (0 + GPT | inds) + (1 | yrs)
Data: hdat
AIC BIC logLik deviance
-550.9679 -512.7889 282.4840 -564.9679
Random effects:
Groups Name Variance Std.Dev.
yrs (Intercept) 2.638e-02 0.162421
inds GPT 3.656e-05 0.006046
inds GST 1.150e-05 0.003391
Residual 3.649e-02 0.191035
Number of obs: 1727, groups: yrs, 94; inds, 20
Fixed effects:
Estimate Std. Error t value
(Intercept) 1.00041 0.01738 57.56
GST 0.07179 0.01844 3.89
GPT 0.03626 0.01849 1.96
Correlation of Fixed Effects:
(Intr) GST
GST -0.001
GPT -0.002 0.320
> (test13f <- lmer(rwi ~ GST + GPT +(1|yrs), hdat, REML=FALSE))
Linear mixed model fit by maximum likelihood ['merMod']
Formula: rwi ~ GST + GPT + (1 | yrs)
Data: hdat
AIC BIC logLik deviance
-554.8997 -527.6289 282.4498 -564.8997
Random effects:
Groups Name Variance Std.Dev.
yrs (Intercept) 0.02637 0.1624
Residual 0.03654 0.1912
Number of obs: 1727, groups: yrs, 94
Fixed effects:
Estimate Std. Error t value
(Intercept) 1.00040 0.01738 57.56
GST 0.07180 0.01843 3.90
GPT 0.03631 0.01844 1.97
Correlation of Fixed Effects:
(Intr) GST
GST -0.001
GPT -0.002 0.321
>
> anova(test13f, test13e, test13d, test13c, test13b)
Data: hdat
Models:
test13f: rwi ~ GST + GPT + (1 | yrs)
test13e: rwi ~ GST + GPT + (0 + GST | inds) + (0 + GPT | inds) + (1 |
test13e: yrs)
test13d: rwi ~ GST * GPT + (0 + GST | inds) + (0 + GPT | inds) + (1 |
test13d: yrs)
test13c: rwi ~ GST * GPT + (0 + GST | inds) + (0 + GPT | inds) + (0 +
test13c: GST:GPT | inds) + (1 | yrs)
test13b: rwi ~ GST * GPT + (1 | inds) + (0 + GST | inds) + (0 + GPT |
test13b: inds) + (0 + GST:GPT | inds) + (1 | yrs)
Df AIC BIC logLik Chisq Chi Df Pr(>Chisq)
test13f 5 -554.90 -527.63 282.45
test13e 7 -550.97 -512.79 282.48 0.0683 2 0.9664
test13d 8 -549.02 -505.38 282.51 0.0481 1 0.8263
test13c 9 -547.38 -498.30 282.69 0.3689 1 0.5436
test13b 10 -545.38 -490.84 282.69 0.0000 1 1.0000
> qqmath(ranef(test13f, postVar = TRUE))
$yrs
>
> proc.time()
user system elapsed
10.980 0.210 11.414
-------------- next part --------------
A non-text attachment was scrubbed...
Name: Rplots.pdf
Type: application/pdf
Size: 142396 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-mixed-models/attachments/20100621/3a106705/attachment.pdf>
More information about the R-sig-mixed-models
mailing list