[R-sig-ME] Convergence problem example

Paul Johnson pauljohn32 at gmail.com
Wed Apr 25 20:18:16 CEST 2018


Thanks. I confirm your suggestion.

There are still some moving pieces I don't understand. Mostly, I
wonder if the convergence warning that I saw is a meaningful thing,
whether users should be notified about it at all.  I'm not a
specialist on this.

The estimates differ at the 6th decimal point. In the fixed estimates
(m1 is default, m1a is bobyqa):

> fixef(m1) - fixef(m1a)
      (Intercept)             lncfs              aiai heiferprimiparous
     2.489050e-05     -5.759644e-06      2.962373e-06     -1.161505e-06

Maybe the convergence warning does find a real issue.

The anova() comparison thinks the models are identical:

> anova(m1a, m1)
Data: dairy
Models:
m1a: fscr ~ lncfs + ai + heifer + (1 | cow)
m1: fscr ~ lncfs + ai + heifer + (1 | cow)
    Df    AIC    BIC  logLik deviance Chisq Chi Df Pr(>Chisq)
m1a  5 4017.8 4047.9 -2003.9   4007.8
m1   5 4017.8 4047.9 -2003.9   4007.8     0      0          1

You have to turn up the digits to find a difference:

> print(anova(m1a, m1), digits = 20)
Data: dairy
Models:
m1a: fscr ~ lncfs + ai + heifer + (1 | cow)
m1: fscr ~ lncfs + ai + heifer + (1 | cow)
        Df                   AIC                     BIC
                  logLik
m1a  5 4017.8128979117764175 4047.8895344568850305 -2003.9064489558882087
m1   5 4017.8128979273369623 4047.8895344724455754 -2003.9064489636684812
                 deviance Chisq Chi Df Pr(>Chisq)
m1a 4007.8128979117764175
m1   4007.8128979273369623     0      0          1
>



pj


On Wed, Apr 25, 2018 at 10:12 AM, Phillip Alday <phillip.alday at mpi.nl> wrote:
> Changing the optimizer seems to address the convergence warnings and has
> no impact on the other aspects of the fit:
>
>> m1 <- update(m1, control=glmerControl(optimizer="bobyqa"))
>> summary(m1)
> Generalized linear mixed model fit by maximum likelihood (Adaptive
>   Gauss-Hermite Quadrature, nAGQ = 30) [glmerMod]
>  Family: binomial  ( logit )
> Formula: fscr ~ lncfs + ai + heifer + (1 | cow)
>    Data: dairy
> Control: glmerControl(optimizer = "bobyqa")
>
>      AIC      BIC   logLik deviance df.resid
>   4017.8   4047.9  -2003.9   4007.8     3022
>
> Scaled residuals:
>     Min      1Q  Median      3Q     Max
> -1.8136 -0.7652 -0.6479  1.0420  1.6213
>
> Random effects:
>  Groups Name        Variance Std.Dev.
>  cow    (Intercept) 0.3419   0.5847
> Number of obs: 3027, groups:  cow, 1575
>
> Fixed effects:
>                   Estimate Std. Error z value Pr(>|z|)
> (Intercept)       -1.55696    0.43481  -3.581 0.000343 ***
> lncfs              0.52222    0.10007   5.218  1.8e-07 ***
> aiai              -1.09560    0.12341  -8.878  < 2e-16 ***
> heiferprimiparous -0.08259    0.09718  -0.850 0.395418
> ---
> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>
> Correlation of Fixed Effects:
>             (Intr) lncfs  aiai
> lncfs       -0.965
> aiai        -0.193 -0.046
> heifrprmprs -0.054  0.016 -0.051
>
> Phillip
>
>
>
> On 04/25/2018 04:56 PM, Paul Johnson wrote:
>> In the book Multilevel and Longitudinal Modeling using Stata,
>> Rabe-Hesketh and Skrondal have a lot of exercises and over the years
>> I've been trying to write Stata and R code to demonstrate
>> similarities/differences.
>>
>> I've run into an example where Stata (either with builtin methods or
>> the addon gllamm) seems to think it gets estimates but lme4
>> diagnostics say there is a convergence failure.  I want to impose on
>> you, get your advice about it. We see convergence warnings quite often
>> with lme4, but they are usually the "rescale your variables" errors,
>> not as blunt as this.
>>
>> The first part retrieves "dairy.dta" from the book website
>>
>> library(foreign)
>> library(lme4)
>>
>> fn <- "dairy"
>> if (!file.exists(paste0(fn, ".dta12"))) {
>>     download.file(paste0("http://www.stata-press.com/data/mlmus3/", fn, ".dta"),
>>                   destfile = paste0(fn, ".dta12"))
>> }
>>
>> dairy <- read.dta(paste0(fn, ".dta12"))
>>
>> #1
>> m1 <- glmer(fscr ~ lncfs + ai + heifer + (1 | cow),
>>             family = binomial(link = "logit"), nAGQ = 30,
>>             data = dairy)
>> summary(m1)
>>
>> From that, I see:
>>> m1 <- glmer(fscr ~ lncfs + ai + heifer + (1 | cow),
>> +             family = binomial(link = "logit"), nAGQ = 30,
>> +             data = dairy)
>> Warning messages:
>> 1: 'rBind' is deprecated.
>>  Since R version 3.2.0, base's rbind() should work fine with S4 objects
>> 2: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
>>   Model failed to converge with max|grad| = 0.00396932 (tol = 0.001,
>> component 1)
>>> summary(m1)
>> Generalized linear mixed model fit by maximum likelihood (Adaptive
>>   Gauss-Hermite Quadrature, nAGQ = 30) [glmerMod]
>>  Family: binomial  ( logit )
>> Formula: fscr ~ lncfs + ai + heifer + (1 | cow)
>>    Data: dairy
>>
>>      AIC      BIC   logLik deviance df.resid
>>   4017.8   4047.9  -2003.9   4007.8     3022
>>
>> Scaled residuals:
>>     Min      1Q  Median      3Q     Max
>> -1.8136 -0.7652 -0.6479  1.0420  1.6213
>>
>> Random effects:
>>  Groups Name        Variance Std.Dev.
>>  cow    (Intercept) 0.3419   0.5847
>> Number of obs: 3027, groups:  cow, 1575
>>
>> Fixed effects:
>>                   Estimate Std. Error z value Pr(>|z|)
>> (Intercept)       -1.55693    0.43481  -3.581 0.000343 ***
>> lncfs              0.52221    0.10007   5.218 1.81e-07 ***
>> aiai              -1.09560    0.12341  -8.878  < 2e-16 ***
>> heiferprimiparous -0.08259    0.09718  -0.850 0.395410
>> ---
>> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>>
>> Correlation of Fixed Effects:
>>             (Intr) lncfs  aiai
>> lncfs       -0.965
>> aiai        -0.193 -0.046
>> heifrprmprs -0.054  0.016 -0.051
>> convergence code: 0
>> Model failed to converge with max|grad| = 0.00396932 (tol = 0.001, component 1)
>>
>>
>> Stata output for comparison purposes, using gllamm:
>>
>> . gllamm fscr lncfs ai heifer, i(cow) link(logit) fam(binom) adapt
>>
>> Running adaptive quadrature
>> Iteration 0:    log likelihood = -2004.6011
>> Iteration 1:    log likelihood = -2003.9085
>> Iteration 2:    log likelihood = -2003.9069
>>
>>
>> Adaptive quadrature has converged, running Newton-Raphson
>> Iteration 0:   log likelihood = -2003.9069
>> Iteration 1:   log likelihood = -2003.9069
>> Iteration 2:   log likelihood = -2003.9064
>> Iteration 3:   log likelihood = -2003.9064
>>
>> number of level 1 units = 3027
>> number of level 2 units = 1575
>>
>> Condition Number = 47.123877
>>
>> gllamm model
>>
>> log likelihood = -2003.9064
>>
>> ------------------------------------------------------------------------------
>>         fscr |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
>> -------------+----------------------------------------------------------------
>>        lncfs |   .5222157   .1000693     5.22   0.000     .3260834     .718348
>>           ai |  -1.095598   .1234099    -8.88   0.000    -1.337477   -.8537193
>>       heifer |  -.0825878   .0971811    -0.85   0.395    -.2730593    .1078837
>>        _cons |  -1.556961   .4348008    -3.58   0.000    -2.409155   -.7047673
>> ------------------------------------------------------------------------------
>>
>>
>> Variances and covariances of random effects
>> ------------------------------------------------------------------------------
>>
>>
>> ***level 2 (cow)
>>
>>     var(1): .34188062 (.1263136)
>> ------------------------------------------------------------------------------
>>
>>
>> Using the newer meglm that is provided with Stata
>>
>> . meglm fscr lncfs ai heifer || cow: , family(binom) link(logit)
>>
>> Fitting fixed-effects model:
>>
>> Iteration 0:   log likelihood = -2011.8253
>> Iteration 1:   log likelihood = -2009.1421
>> Iteration 2:   log likelihood = -2009.1412
>> Iteration 3:   log likelihood = -2009.1412
>>
>> Refining starting values:
>>
>> Grid node 0:   log likelihood = -2015.6021
>>
>> Fitting full model:
>>
>> Iteration 0:   log likelihood = -2015.6021
>> Iteration 1:   log likelihood =  -2006.709
>> Iteration 2:   log likelihood = -2003.9174
>> Iteration 3:   log likelihood = -2003.9065
>> Iteration 4:   log likelihood = -2003.9065
>>
>> Mixed-effects GLM                               Number of obs     =      3,027
>> Family:                binomial
>> Link:                     logit
>> Group variable:             cow                 Number of groups  =      1,575
>>
>>                                                 Obs per group:
>>                                                               min =          1
>>                                                               avg =        1.9
>>                                                               max =          5
>>
>> Integration method: mvaghermite                 Integration pts.  =          7
>>
>>                                                 Wald chi2(3)      =     103.86
>> Log likelihood = -2003.9065                     Prob > chi2       =     0.0000
>> ------------------------------------------------------------------------------
>>         fscr |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
>> -------------+----------------------------------------------------------------
>>        lncfs |   .5222154   .1000694     5.22   0.000      .326083    .7183479
>>           ai |  -1.095597   .1234095    -8.88   0.000    -1.337476   -.8537191
>>       heifer |  -.0825878   .0971811    -0.85   0.395    -.2730593    .1078836
>>        _cons |  -1.556961   .4348012    -3.58   0.000    -2.409155   -.7047659
>> -------------+----------------------------------------------------------------
>> cow          |
>>    var(_cons)|   .3418776   .1263105                      .1657237    .7052721
>> ------------------------------------------------------------------------------
>> LR test vs. logistic model: chibar2(01) = 10.47       Prob >= chibar2 = 0.0006
>>
>>
>> As far as I can see, parameter estimates are the same. I did not
>> compare conditional modes.
>>
>> If lme4's converge warning is valid, then we have a concrete case
>> where Stata is reporting the wrong thing.  I can see some upside there
>> :)
>>
>> R details
>>
>>> sessionInfo()
>> R version 3.4.4 (2018-03-15)
>> Platform: x86_64-pc-linux-gnu (64-bit)
>> Running under: Ubuntu 17.10
>>
>> Matrix products: default
>> BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.7.1
>> LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.7.1
>>
>> locale:
>>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C
>>  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8
>>  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8
>>  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C
>>  [9] LC_ADDRESS=C               LC_TELEPHONE=C
>> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
>>
>> attached base packages:
>> [1] stats     graphics  grDevices utils     datasets  methods   base
>>
>> other attached packages:
>> [1] lme4_1.1-15    Matrix_1.2-14  foreign_0.8-69
>>
>> loaded via a namespace (and not attached):
>>  [1] minqa_1.2.4     MASS_7.3-49     compiler_3.4.4  tools_3.4.4
>>  [5] Rcpp_0.12.15    splines_3.4.4   nlme_3.1-137    grid_3.4.4
>>  [9] nloptr_1.0.4    lattice_0.20-35
>>
>>



-- 
Paul E. Johnson   http://pj.freefaculty.org
Director, Center for Research Methods and Data Analysis http://crmda.ku.edu

To write to me directly, please address me at pauljohn at ku.edu.



More information about the R-sig-mixed-models mailing list