[R-sig-ME] lme4.0 vs. lme4 version 1.1-7: newer version returns inflated SE using glmer() with bobyqa optimizer
Ben Bolker
bbolker at gmail.com
Wed Jul 30 15:57:36 CEST 2014
On 14-07-30 09:12 AM, Chiara Gambi wrote:
> Dear all,
>
> I apologize if this issue has already been covered in the list.
>
> I am analyzing accuracy data on a picture naming task. Correct responses
> are coded as 1, errors as 0. I have 2 fixed-effects predictors:
>
> 1. Condition, 3 levels, within participants and items
> 2. Freq_group: 2 levels, within participants but between items
>
> Here are the contrasts:
>
> contrasts(joint$Condition)<-cbind(c(-2/3,1/3,1/3), c(0,-1/2,1/2))
> contrasts(joint$Freq_group)<-cbind(c(-0.5,0.5))
>
> In addition, I have a between-participants continuous predictor
> (E_proficiency_SP: English Spoken Language Proficiency) that has been
> centered, but not scaled.
>
> I analyzed these data some time ago using a version of lme4 <1. When I ran
> the analyses in R 3.0.3 using lme4.0, I obtained the same output as in my
> old analyses. See below.
>
> ----------------------------------------------------------------------
> Generalized linear mixed model fit by the Laplace approximation
> Formula: Resp.bin ~ 1 + Freq_group * Condition * E_proficiency_SP + (1
> + Freq_group * Condition | Subject) + (1 + Condition | Item_numb)
> Data: joint.NS
> AIC BIC logLik deviance
> 2094 2330 -1008 2016
> Random effects:
> Groups Name Variance Std.Dev.
> Corr
> Item_numb (Intercept) 0.9656165
> 0.982658
> Condition1 1.6699788 1.292277
> -0.035
> Condition2 0.0077218 0.087874 -0.052
> 1.000
> Subject (Intercept) 0.1608360
> 0.401044
> Freq_group1 0.1649723 0.406168
> -0.417
> Condition1 0.4448590 0.666978 -0.678
> 0.269
> Condition2 0.1776654 0.421504 0.274 -0.739
> 0.319
> Freq_group1:Condition1 0.7216097 0.849476 -0.940 0.701 0.647
> -0.478
> Freq_group1:Condition2 0.7144539 0.845254 -0.345 -0.683 0.079
> 0.375 0.010
> Number of obs: 3107, groups: Item_numb, 131; Subject, 24
>
> Fixed effects:
> Estimate Std. Error z value
> Pr(>|z|)
> (Intercept) 2.70027 0.14160 19.070 <
> 2e-16 ***
> Freq_group1 -1.04832 0.24404 -4.296
> 1.74e-05 ***
> Condition1 0.68127 0.22876 2.978
> 0.00290 **
> Condition2 -0.56823 0.20482 -2.774
> 0.00553 **
> E_proficiency_SP 0.04677 0.08473 0.552
> 0.58096
> Freq_group1:Condition1 -0.32467 0.40460 -0.802
> 0.42229
> Freq_group1:Condition2 -0.10520 0.40897 -0.257
> 0.79700
> Freq_group1:E_proficiency_SP -0.20632 0.12643 -1.632
> 0.10268
> Condition1:E_proficiency_SP -0.13266 0.15364 -0.863
> 0.38791
> Condition2:E_proficiency_SP -0.39620 0.16135 -2.456
> 0.01407 *
> Freq_group1:Condition1:E_proficiency_SP -0.26202 0.25468 -1.029
> 0.30357
> Freq_group1:Condition2:E_proficiency_SP 0.14582 0.31966 0.456
> 0.64827
> ---
> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>
> Correlation of Fixed Effects:
> (Intr) Frq_g1 Cndtn1 Cndtn2 E_p_SP Fr_1:C1 Fr_1:C2 F_1:E_
> C1:E__ C2:E__ F_1:C1:
> Freq_group1
> -0.202
>
> Condition1 -0.155
> -0.008
> Condition2 -0.032 -0.048
> 0.010
> E_prfcnc_SP 0.033 -0.052 0.007
> -0.094
> Frq_grp1:C1 -0.299 0.214 0.032 -0.038
> -0.037
> Frq_grp1:C2 -0.036 -0.214 0.064 -0.252 0.077
> -0.075
> Frq_1:E__SP -0.061 0.052 -0.043 0.101 -0.323 0.012
> -0.126
> Cndt1:E__SP 0.006 -0.032 0.044 -0.077 -0.281 -0.045 0.063
> 0.035
> Cndt2:E__SP -0.079 0.071 -0.073 0.068 -0.037 0.065 -0.119 -0.119
> -0.017
> F_1:C1:E__S -0.034 0.008 -0.050 0.074 -0.443 0.060 -0.093 0.343
> 0.103 -0.078
> F_1:C2:E__S 0.063 -0.094 0.059 -0.117 -0.079 -0.087 0.066 -0.317
> 0.039 -0.186 -0.130
> ----------------------------------------------------------------------------------------------------------------
>
> However, having read recently that lmer version 1.1-7 obtains better fits
> than lme4.0, I fitted the same model using glmer with the bobyqa optimizer
> and maxfun=16000. Here is the output. There were no convergence issues or
> any other warnings (when I increased the number of iterations as suggested).
>
> -------------------------------------------------------------------------------------------------------------------
> Generalized linear mixed model fit by maximum likelihood (Laplace
> Approximation) ['glmerMod']
> Family: binomial ( logit )
> Formula: Resp.bin ~ 1 + Freq_group * Condition * E_proficiency_SP + (1
> + Freq_group * Condition | Subject) + (1 + Condition | Item_numb)
> Data: joint.NS
> Control: glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 16000))
>
> AIC BIC logLik deviance df.resid
> 2094.4 2330.0 -1008.2 2016.4 3068
>
> Scaled residuals:
> Min 1Q Median 3Q Max
> -6.7090 0.1461 0.2251 0.3482 1.6999
>
> Random effects:
> Groups Name Variance Std.Dev.
> Corr
> Item_numb (Intercept) 0.965582
> 0.98264
> Condition1 1.669926 1.29226
> -0.04
> Condition2 0.007722 0.08788 -0.05
> 1.00
> Subject (Intercept) 0.160826
> 0.40103
> Freq_group1 0.164960 0.40615
> -0.42
> Condition1 0.444855 0.66697 -0.68
> 0.27
> Condition2 0.177659 0.42150 0.27 -0.74
> 0.32
> Freq_group1:Condition1 0.721550 0.84944 -0.94 0.70 0.65
> -0.48
> Freq_group1:Condition2 0.714430 0.84524 -0.35 -0.68 0.08
> 0.37 0.01
> Number of obs: 3107, groups: Item_numb, 131; Subject, 24
>
> Fixed effects:
> Estimate Std. Error z value
> Pr(>|z|)
> (Intercept) 2.70023 0.16124 16.747 <
> 2e-16 ***
> Freq_group1 -1.04829 0.25693 -4.080
> 4.5e-05 ***
> Condition1 0.68123 0.27441 2.483
> 0.0130 *
> Condition2 -0.56822 0.29971 -1.896
> 0.0580 .
> E_proficiency_SP 0.04675 0.08527 0.548
> 0.5835
> Freq_group1:Condition1 -0.32464 0.43299 -0.750
> 0.4534
> Freq_group1:Condition2 -0.10517 0.44365 -0.237
> 0.8126
> Freq_group1:E_proficiency_SP -0.20633 0.12692 -1.626
> 0.1040
> Condition1:E_proficiency_SP -0.13266 0.15570 -0.852
> 0.3942
> Condition2:E_proficiency_SP -0.39620 0.16214 -2.443
> 0.0145 *
> Freq_group1:Condition1:E_proficiency_SP -0.26201 0.25735 -1.018
> 0.3086
> Freq_group1:Condition2:E_proficiency_SP 0.14583 0.31962 0.456
> 0.6482
> ---
> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>
> Correlation of Fixed Effects:
> (Intr) Frq_g1 Cndtn1 Cndtn2 E_p_SP Fr_1:C1 Fr_1:C2 F_1:E_
> C1:E__ C2:E__ F_1:C1:
> Freq_group1
> -0.249
>
> Condition1 -0.032
> -0.056
> Condition2 -0.085 -0.059
> -0.056
> E_prfcnc_SP 0.046 -0.055 0.001
> -0.094
> Frq_grp1:C1 -0.289 0.211 -0.065 -0.045
> -0.028
> Frq_grp1:C2 -0.042 -0.174 0.036 -0.328 0.079
> -0.053
> Frq_1:E__SP -0.060 0.047 -0.044 0.080 -0.313 0.019
> -0.117
> Cndt1:E__SP 0.001 -0.024 0.061 -0.070 -0.291 -0.048 0.056
> 0.031
> Cndt2:E__SP -0.092 0.071 -0.074 0.069 -0.032 0.060 -0.122 -0.122
> -0.012
> F_1:C1:E__S -0.038 0.015 -0.047 0.059 -0.443 0.044 -0.083 0.321
> 0.116 -0.077
> F_1:C2:E__S 0.061 -0.081 0.050 -0.115 -0.077 -0.071 0.076 -0.305
> 0.035 -0.177 -0.124
> ---------------------------------------------------------------------------------------------------------------
>
> The two models are very similar, with almost identical estimates and the
> same AIC, BIC, and logLik. However, the model fitted with the newer glmer
> has generally higher standard errors.
>
> Any idea as to what might be behind this? And is there any reason to prefer
> one method to the other a priori?
>
> Many thanks,
>
> Chiara Gambi
My guess is that it's related to this issue:
https://github.com/lme4/lme4/issues/47
>From https://github.com/lme4/lme4/blob/master/inst/NEWS.Rd
\item Standard errors of fixed effects are now computed
from the approximate Hessian by default (see the
\code{use.hessian} argument in \code{vcov.merMod}); this
gives better (correct) answers when the estimates of
the random- and fixed-effect parameters are correlated
(Github #47)
}
If I'm correct, you should be able to see the difference in
standard errors by comparing
sqrt(diag(vcov(fitted_model,use.hessian=TRUE))) ## new default
sqrt(diag(vcov(fitted_model,use.hessian=FALSE))) ## old default
As far as we know, the new standard errors are better/more correct.
Ben Bolker
More information about the R-sig-mixed-models
mailing list