Closed lukeCLarter closed 1 week ago
bbolker
commented
1 year ago Sorry not to respond sooner. Can you post or send a link to the data? The coefficient estimates don't look that different, but the standard errors from the "new R" version definitely look wonky.
lukeCLarter
commented
1 year ago No worries! Here is the csv and the rmarkdown file I used to run the code. Thanks again for your help!
On Wed, Apr 12, 2023 at 8:43 AM Ben Bolker @.***> wrote:
Sorry not to respond sooner. Can you post or send a link to the data? The coefficient estimates don't look that different, but the standard errors from the "new R" version definitely look wonky.
— Reply to this email directly, view it on GitHub https://github.com/lme4/lme4/issues/720#issuecomment-1505303221, or unsubscribe https://github.com/notifications/unsubscribe-auth/APFYELKC4NPHJ5MRRXCB3KTXA2WORANCNFSM6AAAAAAWEJNP64 . You are receiving this because you authored the thread.Message ID: @.***>
lukeCLarter
commented
1 year ago Apologies for rushing you, but the deadline for me to submit my manuscript following receiving the reviews is approaching, so I would very much appreciate a response by the end of the weekend if possible. Thanks again for your time!
On Wed, Apr 12, 2023 at 4:44 PM Luke C Larter @.***> wrote:
No worries! Here is the csv and the rmarkdown file I used to run the code. Thanks again for your help!
On Wed, Apr 12, 2023 at 8:43 AM Ben Bolker @.***> wrote:
Sorry not to respond sooner. Can you post or send a link to the data? The coefficient estimates don't look that different, but the standard errors from the "new R" version definitely look wonky.
— Reply to this email directly, view it on GitHub https://github.com/lme4/lme4/issues/720#issuecomment-1505303221, or unsubscribe https://github.com/notifications/unsubscribe-auth/APFYELKC4NPHJ5MRRXCB3KTXA2WORANCNFSM6AAAAAAWEJNP64 . You are receiving this because you authored the thread.Message ID: @.***>
bbolker
commented
1 year ago Oops. You responded by e-mail to github, so attachments were stripped. Please re-send by e-mail to maintainer("lme4") ?
bbolker
commented
1 year ago So far I can't reproduce this with new R (r-devel) and new (development) lme4. Since you said the problem seems to be with R version rather than lme4 version, I haven't tried alternative versions of lme4. I do wonder whether it might be an issue with versions of the Matrix package, which has also been under development ... ?? It could also be OS-specific. Can you please run the following minimal code (including summary()/sessionInfo()) on your system and post results ... ?
(assumes availability of data file "free_chucks_model_data.csv")
library(lme4)
data <- read.csv('free_chucks_model_data.csv')
data$chucks_free_rest_binary=ifelse(data$chucks_free_rest>0, 1, 0)
#Full model with interaction and random slope
full <- glmer(chucks_free_rest_binary ~ call_complexity*chorus_size + (1|chorus_name) +
(1+call_complexity|chorus_name:toe_clip_number),
family = binomial(link="logit"), data=data, control = glmerControl(optimizer ="bobyqa"))
summary(full)
sessionInfo()Generalized linear mixed model fit by maximum likelihood (Laplace
Approximation) [glmerMod]
Family: binomial ( logit )
Formula: chucks_free_rest_binary ~ call_complexity * chorus_size + (1 |
chorus_name) + (1 + call_complexity | chorus_name:toe_clip_number)
Data: data
Control: glmerControl(optimizer = "bobyqa")
AIC BIC logLik deviance df.resid
1692.6 1738.0 -838.3 1676.6 2132
Scaled residuals:
Min 1Q Median 3Q Max
-5.6471 0.1582 0.2877 0.4378 1.4995
Random effects:
Groups Name Variance Std.Dev. Corr
chorus_name:toe_clip_number (Intercept) 2.16683 1.472
call_complexity 0.74820 0.865 -0.90
chorus_name (Intercept) 0.03762 0.194
Number of obs: 2140, groups: chorus_name:toe_clip_number, 72; chorus_name, 16
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 5.4742 1.1325 4.834 1.34e-06 ***
call_complexity -0.4734 0.6853 -0.691 0.4896
chorus_size -1.0391 0.2219 -4.683 2.82e-06 ***
call_complexity:chorus_size 0.2634 0.1327 1.985 0.0471 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) cll_cm chrs_s
cll_cmplxty -0.938
chorus_size -0.962 0.900
cll_cmplx:_ 0.910 -0.965 -0.939I am particularly interested in seeing platform, Matrix products, version of lme4 and Matrix
> R Under development (unstable) (2023-04-20 r84291)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Pop!_OS 22.04 LTS
Matrix products: default
BLAS/LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so; LAPACK version 3.10.0
[locale info removed]
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] lme4_1.1-33 Matrix_1.5-4
loaded via a namespace (and not attached):
[1] minqa_1.2.5 MASS_7.3-58.4 compiler_4.4.0 Rcpp_1.0.10 splines_4.4.0
[6] nlme_3.1-162 grid_4.4.0 nloptr_2.0.3 boot_1.3-28.1 lattice_0.21-8Here is the output in ver 4.1.2:
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
Family: binomial ( logit )
Formula: chucks_free_rest_binary ~ call_complexity * chorus_size + (1 | chorus_name) + (1 + call_complexity | chorus_name:toe_clip_number)
Data: data
Control: glmerControl(optimizer = "bobyqa")
AIC BIC logLik deviance df.resid
1692.6 1738.0 -838.3 1676.6 2132
Scaled residuals:
Min 1Q Median 3Q Max
-5.6471 0.1582 0.2877 0.4378 1.4995
Random effects:
Groups Name Variance Std.Dev. Corr
chorus_name:toe_clip_number (Intercept) 2.16683 1.472
call_complexity 0.74820 0.865 -0.90
chorus_name (Intercept) 0.03762 0.194
Number of obs: 2140, groups: chorus_name:toe_clip_number, 72; chorus_name, 16
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 5.4742 1.1330 4.832 1.35e-06 ***
call_complexity -0.4734 0.6855 -0.691 0.4898
chorus_size -1.0391 0.2220 -4.681 2.85e-06 ***
call_complexity:chorus_size 0.2634 0.1327 1.984 0.0472 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) cll_cm chrs_s
cll_cmplxty -0.938
chorus_size -0.962 0.900
cll_cmplx:_ 0.910 -0.965 -0.939
R version 4.1.2 (2021-11-01)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 22621)
Matrix products: default
locale:
[1] LC_COLLATE=English_United States.1252 LC_CTYPE=English_United States.1252 LC_MONETARY=English_United States.1252 LC_NUMERIC=C
[5] LC_TIME=English_United States.1252
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] lmerTest_3.1-3 lme4_1.1-27.1 Matrix_1.3-4
loaded via a namespace (and not attached):
[1] tidyselect_1.1.1 sjlabelled_1.1.8 xfun_0.27 sjPlot_2.8.10 performance_0.8.0 purrr_0.3.4 splines_4.1.2 lattice_0.20-45
[9] parameters_0.15.0 colorspace_2.0-2 vctrs_0.5.1 generics_0.1.1 utf8_1.2.2 rlang_1.0.6 nloptr_1.2.2.3 pillar_1.8.1
[17] glue_1.4.2 DBI_1.1.2 effectsize_0.5 modelr_0.1.8 emmeans_1.8.0 lifecycle_1.0.3 sjmisc_2.8.7 munsell_0.5.0
[25] gtable_0.3.0 bayestestR_0.11.5 mvtnorm_1.1-3 coda_0.19-4 knitr_1.37 datawizard_0.2.1 fansi_0.5.0 broom_0.7.10
[33] Rcpp_1.0.7 xtable_1.8-4 DHARMa_0.4.6 scales_1.1.1 backports_1.3.0 ggeffects_1.1.3 ggplot2_3.3.5 insight_0.18.5
[41] dplyr_1.0.10 numDeriv_2016.8-1.1 grid_4.1.2 cli_3.1.0 tools_4.1.2 sjstats_0.18.1 magrittr_2.0.1 tibble_3.1.5
[49] tidyr_1.1.4 pkgconfig_2.0.3 MASS_7.3-54 ellipsis_0.3.2 estimability_1.4.1 assertthat_0.2.1 minqa_1.2.4 rstudioapi_0.13
[57] R6_2.5.1 boot_1.3-28 nlme_3.1-153 compiler_4.1.2And here it is for ver 4.2.2:
Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
Model failed to converge with max|grad| = 0.0739749 (tol = 0.002, component 1)
Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
Model is nearly unidentifiable: very large eigenvalue
- Rescale variables?
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
Family: binomial ( logit )
Formula: chucks_free_rest_binary ~ call_complexity * chorus_size + (1 | chorus_name) + (1 + call_complexity | chorus_name:toe_clip_number)
Data: data
Control: glmerControl(optimizer = "bobyqa")
AIC BIC logLik deviance df.resid
1692.6 1738.0 -838.3 1676.6 2132
Scaled residuals:
Min 1Q Median 3Q Max
-5.6810 0.1610 0.2856 0.4378 1.4950
Random effects:
Groups Name Variance Std.Dev. Corr
chorus_name:toe_clip_number (Intercept) 2.17984 1.4764
call_complexity 0.75669 0.8699 -0.90
chorus_name (Intercept) 0.03472 0.1863
Number of obs: 2140, groups: chorus_name:toe_clip_number, 72; chorus_name, 16
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 5.3246761 0.0009618 5536.0 <2e-16 ***
call_complexity -0.3864900 0.0009624 -401.6 <2e-16 ***
chorus_size -1.0107972 0.0009624 -1050.2 <2e-16 ***
call_complexity:chorus_size 0.2478955 0.0009634 257.3 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) cll_cm chrs_s
cll_cmplxty 0.000
chorus_size 0.000 0.001
cll_cmplx:_ 0.000 0.000 -0.001
optimizer (bobyqa) convergence code: 0 (OK)
Model failed to converge with max|grad| = 0.0739749 (tol = 0.002, component 1)
Model is nearly unidentifiable: very large eigenvalue
- Rescale variables?
R version 4.2.2 (2022-10-31 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 22621)
Matrix products: default
locale:
[1] LC_COLLATE=English_United States.1252 LC_CTYPE=English_United States.1252 LC_MONETARY=English_United States.1252 LC_NUMERIC=C
[5] LC_TIME=English_United States.1252
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] lmerTest_3.1-3 lme4_1.1-31 Matrix_1.5-1
loaded via a namespace (and not attached):
[1] tidyselect_1.2.0 sjlabelled_1.2.0 xfun_0.37 sjPlot_2.8.12 performance_0.10.2 purrr_1.0.1 splines_4.2.2 lattice_0.20-45
[9] colorspace_2.1-0 vctrs_0.5.2 generics_0.1.3 htmltools_0.5.4 yaml_2.3.7 utf8_1.2.3 rlang_1.0.6 nloptr_2.0.3
[17] pillar_1.8.1 glue_1.6.2 modelr_0.1.10 emmeans_1.8.4-1 lifecycle_1.0.3 sjmisc_2.8.9 munsell_0.5.0 gtable_0.3.1
[25] bayestestR_0.13.0 mvtnorm_1.1-3 evaluate_0.20 coda_0.19-4 knitr_1.42 fastmap_1.1.1 datawizard_0.6.5 fansi_1.0.4
[33] broom_1.0.3 Rcpp_1.0.10 xtable_1.8-4 scales_1.2.1 backports_1.4.1 DHARMa_0.4.6 ggeffects_1.2.0 digest_0.6.31
[41] ggplot2_3.4.1 insight_0.19.0 dplyr_1.1.0 numDeriv_2016.8-1.1 grid_4.2.2 cli_3.6.0 tools_4.2.2 sjstats_0.18.2
[49] magrittr_2.0.3 tibble_3.1.8 tidyr_1.3.0 pkgconfig_2.0.3 MASS_7.3-58.1 estimability_1.4.1 minqa_1.2.5 rmarkdown_2.20
[57] rstudioapi_0.14 R6_2.5.1 boot_1.3-28 nlme_3.1-160 compiler_4.2.2 In case I haven't said so explicitly yet, the only thing that seems problematic (and really bad) are the standard error estimates under the 'new' (R 4.2.2) configuration. The coefficients are not very different (the call complexity coefficient seems to be the most different, -0.39 vs -0.47, and according to the more plausible standard error estimates, the SE of this coefficient is about 0.6, so this is a difference of about 1/7 of a standard error, so doesn't seem like a big deal).
The crazy SE estimates are definitely concerning me; I'd like to figure out what's going on.
I just ran this under an R Docker ('rocker') container using R version 4.2.2 with the specified package versions and got the correct/sensible standard error estimates (details below). Short of installing a Windows virtual machine and then installing the relevant R/package versions, it's not clear how much more debugging I can do on this end ...
Can you try the following code with the 'bad' configuration (if necessary installing the dfoptim and optimx packages first) and see what the results are? (In particular, I want to know if the bad SD estimates occur with all optimizers, or just with a subset of them ...)
aa <- allFit(full)
summary(aa)$sdcor(If nothing else, if it works for some optimizers with the newer version of R then you could switch ...)
install.packages(c("remotes", "dfoptim", "optimx"))
## good: R 4.1.2, lme4 1.1-27.1, Matrix 1.3-4
## bad : R 4.2.2, lme4 1.1-31, Matrix 1.5-1
remotes::install_version("lme4", "1.1-31")
remotes::install_version("Matrix", "1.5-1")docker run --rm -v ~/R/misc:/home/r/ -ti rocker/r-ver:4.2.2 ## from glmer
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 5.4742 1.1325 4.834 1.34e-06 ***
call_complexity -0.4734 0.6853 -0.691 0.4896
chorus_size -1.0391 0.2219 -4.683 2.82e-06 ***
call_complexity:chorus_size 0.2634 0.1327 1.985 0.0471 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
## from sessionInfo():
R version 4.2.2 (2022-10-31)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 22.04.2 LTS
Matrix products: default
BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] lme4_1.1-31 Matrix_1.5-1
loaded via a namespace (and not attached):
[1] minqa_1.2.5 MASS_7.3-58.1 compiler_4.2.2 tools_4.2.2
[5] Rcpp_1.0.10 remotes_2.4.2 splines_4.2.2 nlme_3.1-160
[9] grid_4.2.2 nloptr_2.0.3 boot_1.3-28 lattice_0.20-45Ok weird that you got the sensible results when running it on the docker with ver 4.2.2...
I ran allfit and the results are below. All of the optimizers give some kind of warning (see second chunk).
Here's jusr the $sdcor from all fit (full allfit results in next chunk):
$sdcor
chorus_name:toe_clip_number.(Intercept) chorus_name:toe_clip_number.call_complexity.(Intercept) chorus_name:toe_clip_number.call_complexity
bobyqa 1.476429 0.8698804 -0.8996158
Nelder_Mead 1.479415 0.8694965 -0.8997842
nlminbwrap 1.466426 0.8731048 -0.8986772
optimx.L-BFGS-B 1.451177 0.8588845 -0.8975563
nloptwrap.NLOPT_LN_NELDERMEAD 1.468504 0.8665632 -0.8984417
nloptwrap.NLOPT_LN_BOBYQA 1.471985 0.8649975 -0.8997654
chorus_name.(Intercept)
bobyqa 0.1863221
Nelder_Mead 0.1863085
nlminbwrap 0.1924563
optimx.L-BFGS-B 0.1877743
nloptwrap.NLOPT_LN_NELDERMEAD 0.1954401
nloptwrap.NLOPT_LN_BOBYQA 0.1933269Here are the full allfit results. All optimizers give some kind of warning, but also give very similar results:
$which.OK
bobyqa Nelder_Mead nlminbwrap nmkbw optimx.L-BFGS-B nloptwrap.NLOPT_LN_NELDERMEAD
TRUE TRUE TRUE FALSE TRUE TRUE
nloptwrap.NLOPT_LN_BOBYQA
TRUE
$msgs
$msgs$bobyqa
$msgs$bobyqa[[1]]
[1] "Model failed to converge with max|grad| = 0.0739749 (tol = 0.002, component 1)"
$msgs$bobyqa[[2]]
[1] "Model is nearly unidentifiable: very large eigenvalue\n - Rescale variables?"
$msgs$Nelder_Mead
$msgs$Nelder_Mead[[1]]
[1] "Model failed to converge with max|grad| = 0.0739521 (tol = 0.002, component 1)"
$msgs$Nelder_Mead[[2]]
[1] "Model is nearly unidentifiable: very large eigenvalue\n - Rescale variables?"
$msgs$nlminbwrap
$msgs$nlminbwrap[[1]]
[1] "Model failed to converge with max|grad| = 0.0755479 (tol = 0.002, component 1)"
$msgs$nlminbwrap[[2]]
[1] "Model is nearly unidentifiable: very large eigenvalue\n - Rescale variables?"
$msgs$`optimx.L-BFGS-B`
[1] "Model failed to converge with max|grad| = 0.466605 (tol = 0.002, component 1)"
$msgs$nloptwrap.NLOPT_LN_NELDERMEAD
$msgs$nloptwrap.NLOPT_LN_NELDERMEAD[[1]]
[1] "Model failed to converge with max|grad| = 0.0748281 (tol = 0.002, component 1)"
$msgs$nloptwrap.NLOPT_LN_NELDERMEAD[[2]]
[1] "Model is nearly unidentifiable: very large eigenvalue\n - Rescale variables?"
$msgs$nloptwrap.NLOPT_LN_BOBYQA
[1] "Model failed to converge with max|grad| = 0.0048812 (tol = 0.002, component 1)"
$fixef
(Intercept) call_complexity chorus_size call_complexity:chorus_size
bobyqa 5.324676 -0.3864900 -1.010797 0.2478955
Nelder_Mead 5.388769 -0.4141549 -1.022838 0.2525267
nlminbwrap 5.424477 -0.4523884 -1.033613 0.2621747
optimx.L-BFGS-B 5.277750 -0.3623151 -1.005134 0.2442498
nloptwrap.NLOPT_LN_NELDERMEAD 5.436028 -0.4438554 -1.032669 0.2589604
nloptwrap.NLOPT_LN_BOBYQA 5.471997 -0.4722368 -1.038680 0.2631687
$llik
bobyqa Nelder_Mead nlminbwrap optimx.L-BFGS-B nloptwrap.NLOPT_LN_NELDERMEAD nloptwrap.NLOPT_LN_BOBYQA
-838.3173 -838.3108 -838.3149 -838.3346 -838.3084 -838.3166
$sdcor
chorus_name:toe_clip_number.(Intercept) chorus_name:toe_clip_number.call_complexity.(Intercept) chorus_name:toe_clip_number.call_complexity
bobyqa 1.476429 0.8698804 -0.8996158
Nelder_Mead 1.479415 0.8694965 -0.8997842
nlminbwrap 1.466426 0.8731048 -0.8986772
optimx.L-BFGS-B 1.451177 0.8588845 -0.8975563
nloptwrap.NLOPT_LN_NELDERMEAD 1.468504 0.8665632 -0.8984417
nloptwrap.NLOPT_LN_BOBYQA 1.471985 0.8649975 -0.8997654
chorus_name.(Intercept)
bobyqa 0.1863221
Nelder_Mead 0.1863085
nlminbwrap 0.1924563
optimx.L-BFGS-B 0.1877743
nloptwrap.NLOPT_LN_NELDERMEAD 0.1954401
nloptwrap.NLOPT_LN_BOBYQA 0.1933269
$theta
chorus_name:toe_clip_number.(Intercept) chorus_name:toe_clip_number.call_complexity.(Intercept) chorus_name:toe_clip_number.call_complexity
bobyqa 1.476429 -0.7825582 0.3798613
Nelder_Mead 1.479415 -0.7823592 0.3793919
nlminbwrap 1.466426 -0.7846394 0.3829530
optimx.L-BFGS-B 1.451177 -0.7708972 0.3786821
nloptwrap.NLOPT_LN_NELDERMEAD 1.468504 -0.7785565 0.3805016
nloptwrap.NLOPT_LN_BOBYQA 1.471985 -0.7782949 0.3774623
chorus_name.(Intercept)
bobyqa 0.1863221
Nelder_Mead 0.1863085
nlminbwrap 0.1924563
optimx.L-BFGS-B 0.1877743
nloptwrap.NLOPT_LN_NELDERMEAD 0.1954401
nloptwrap.NLOPT_LN_BOBYQA 0.1933269
$times
user.self sys.self elapsed user.child sys.child
bobyqa 1.54 0.00 3.29 NA NA
Nelder_Mead 3.45 0.01 7.25 NA NA
nlminbwrap 2.50 0.03 4.55 NA NA
optimx.L-BFGS-B 1.28 0.02 4.21 NA NA
nloptwrap.NLOPT_LN_NELDERMEAD 1.06 0.00 3.16 NA NA
nloptwrap.NLOPT_LN_BOBYQA 1.27 0.00 2.83 NA NA
$feval
bobyqa Nelder_Mead nlminbwrap optimx.L-BFGS-B nloptwrap.NLOPT_LN_NELDERMEAD nloptwrap.NLOPT_LN_BOBYQA
430 826 NA 62 695 831
attr(,"class")
[1] "summary.allFit"If it's useful, when specifying different optimizers in my actual model, the available ones all give the weird SEs: Bobyqa:
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
Family: binomial ( logit )
Formula: chucks_free_rest_binary ~ call_complexity * chorus_size + (1 | chorus_name) + (1 + call_complexity | chorus_name:toe_clip_number)
Data: data
Control: glmerControl(optimizer = "bobyqa")
AIC BIC logLik deviance df.resid
1692.6 1738.0 -838.3 1676.6 2132
Scaled residuals:
Min 1Q Median 3Q Max
-5.6810 0.1610 0.2856 0.4378 1.4950
Random effects:
Groups Name Variance Std.Dev. Corr
chorus_name:toe_clip_number (Intercept) 2.17984 1.4764
call_complexity 0.75669 0.8699 -0.90
chorus_name (Intercept) 0.03472 0.1863
Number of obs: 2140, groups: chorus_name:toe_clip_number, 72; chorus_name, 16
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 5.3246761 0.0009618 5536.0 <2e-16 ***
call_complexity -0.3864900 0.0009624 -401.6 <2e-16 ***
chorus_size -1.0107972 0.0009624 -1050.2 <2e-16 ***
call_complexity:chorus_size 0.2478955 0.0009634 257.3 <2e-16 ***Nelder_Mead:
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
Family: binomial ( logit )
Formula: chucks_free_rest_binary ~ call_complexity * chorus_size + (1 | chorus_name) + (1 + call_complexity | chorus_name:toe_clip_number)
Data: data
Control: glmerControl(optimizer = "Nelder_Mead")
AIC BIC logLik deviance df.resid
1692.6 1738.0 -838.3 1676.6 2132
Scaled residuals:
Min 1Q Median 3Q Max
-5.6804 0.1594 0.2865 0.4384 1.4998
Random effects:
Groups Name Variance Std.Dev. Corr
chorus_name:toe_clip_number (Intercept) 2.18867 1.4794
call_complexity 0.75602 0.8695 -0.90
chorus_name (Intercept) 0.03471 0.1863
Number of obs: 2140, groups: chorus_name:toe_clip_number, 72; chorus_name, 16
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 5.3887689 0.0009608 5608.4 <0.0000000000000002 ***
call_complexity -0.4141549 0.0009614 -430.8 <0.0000000000000002 ***
chorus_size -1.0228383 0.0009613 -1064.0 <0.0000000000000002 ***
call_complexity:chorus_size 0.2525267 0.0009621 262.5 <0.0000000000000002 ***nlminbwrap:
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
Family: binomial ( logit )
Formula: chucks_free_rest_binary ~ call_complexity * chorus_size + (1 | chorus_name) + (1 + call_complexity | chorus_name:toe_clip_number)
Data: data
Control: glmerControl(optimizer = "nlminbwrap")
AIC BIC logLik deviance df.resid
1692.6 1738.0 -838.3 1676.6 2132
Scaled residuals:
Min 1Q Median 3Q Max
-5.6760 0.1599 0.2848 0.4364 1.4949
Random effects:
Groups Name Variance Std.Dev. Corr
chorus_name:toe_clip_number (Intercept) 2.15041 1.4664
call_complexity 0.76231 0.8731 -0.90
chorus_name (Intercept) 0.03704 0.1925
Number of obs: 2140, groups: chorus_name:toe_clip_number, 72; chorus_name, 16
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 5.4244771 0.0009398 5772.2 <0.0000000000000002 ***
call_complexity -0.4523884 0.0009403 -481.1 <0.0000000000000002 ***
chorus_size -1.0336132 0.0009404 -1099.2 <0.0000000000000002 ***
call_complexity:chorus_size 0.2621747 0.0009415 278.5 <0.0000000000000002 ***OK, something weird is going on here: I don't understand why/how allFit would be giving different answers from the individual fits themselves. Are you willing to save the results of allFit(), along with the individual fits, in a single .rds file and send them by e-mail ... ? (save(list = c("obj1", "obj2", "obj3"), file = "lme4_GH720_results.rds"))
bbolker
commented
1 year ago I've localized the problem somewhat, although I'm still not clear where it's failing.
## note: renamed 'data' to 'dat'
## fit with glmmTMB for comparison
library(glmmTMB)
gt <- glmmTMB(chucks_free_rest_binary ~ call_complexity*chorus_size + (1|chorus_name) + (1+call_complexity|chorus_name:toe_clip_number),
family = binomial(link="logit"), data=dat)
## helper functions
pp <- function(x) {
unname(signif(sqrt(diag(x)),3))
}
## copy of internal vcov-from-hessian function
ntheta <- 4
calc.vcov.hess <- function(h) {
## invert 2*Hessian, catching errors and forcing symmetric result
## ~= forceSymmetric(solve(h/2)[i,i]) : solve(h/2) = 2*solve(h)
h <- tryCatch(solve(h),
error=function(e) matrix(NA_real_,nrow=nrow(h),ncol=ncol(h)))
i <- -seq_len(ntheta) ## drop var-cov parameters
h <- h[i,i]
forceSymmetric(h + t(h))
}
## get deviance function
fn <- getME(obj1, "devfun")
pars <- unlist(getME(obj1, c("theta", "beta")))
## compute hessian (2d-deriv matrix) 3 ways
H1 <- optimHess(pars, fn, control = list(ndeps = rep(1e-4, length(pars))))
H2 <- lme4:::deriv12(fn, pars)$Hessian
H3 <- numDeriv::hessian(fn, pars)
cbind(
optimHess = pp(calc.vcov.hess(H1)), ## from base R
deriv12 = pp(calc.vcov.hess(H2)), ## what's actually used in lme4
numDeriv = pp(calc.vcov.hess(H3)), ## should (??) be more accurate, but runs into trouble?
old = pp(calc.vcov.hess(obj1@optinfo$derivs$Hessian)), ## stored value
RX = pp(sigma(obj1)^2*obj1@pp$unsc()), ## internal comp (see lmer vignette eq. 55)
glmmTMB = pp(vcov(gt)$cond)) ## glmmTMB, for comparisonResults
optimHess deriv12 numDeriv old RX glmmTMB
[1,] 1.130 1.130 NaN 0.000962 1.150 1.160
[2,] 0.686 0.686 NaN 0.000962 0.704 0.702
[3,] 0.221 0.221 NaN 0.000962 0.225 0.226
[4,] 0.133 0.133 NaN 0.000963 0.136 0.136The same thing happens for most but not all of the optimizers tried in allFit() (I was mistakenly looking at the $sdcor component before, which is wrong! Those are RE SDs and correlations, not SEs of fixed-effect estimates ...)
ss <- sapply(allfit[-4], function(x) sqrt(diag(vcov(x))))
rownames(ss) <- NULL
print(ss, width = 1000)
bobyqa Nelder_Mead nlminbwrap optimx.L-BFGS-B nloptwrap.NLOPT_LN_NELDERMEAD nloptwrap.NLOPT_LN_BOBYQA
[1,] 0.0009618270 0.0009608303 0.0009397562 1.1216788 0.0009463433 1.1332905
[2,] 0.0009623890 0.0009613592 0.0009402971 0.6829854 0.0009468645 0.6857589
[3,] 0.0009624384 0.0009613403 0.0009403515 0.2199460 0.0009468828 0.2220086
[4,] 0.0009633823 0.0009621015 0.0009414895 0.1322074 0.0009477788 0.1327727summary(obj1, use.hessian = FALSE) you get reasonable answers (but a warning telling you that the finite-differences Hessian and the RX hessian differ significantly ...)lme4:::deriv12 on Windows ...numDeriv::hessian() (which should be better than finite differences, as it uses Richardson extrapolation) goes off the rails (the computed Hessian is far from positive definite)The only other thing I've noticed about this problem is that the fixed-effect parameters are strongly correlated (abs(corr) ranges from 0.91 to 0.97, not including self-correlations ...). The condition number of the 'good' Hessians appears to be about 4000 - is that enough that we should get crazy answers?
Note that these problems only apply to the SE estimates; the point estimates of the coefficients all seem fine.
aa <- t(sapply(allfit[-4], fixef))
rownames(aa) <- paste("AF", rownames(aa), sep = ".")
print(width=1000, rbind(
bobyqa = fixef(obj1),
NM = fixef(obj2),
nlminb = fixef(obj3),
aa))
(Intercept) call_complexity chorus_size call_complexity:chorus_size
bobyqa 5.324676 -0.3864900 -1.010797 0.2478955
NM 5.388769 -0.4141549 -1.022838 0.2525267
nlminb 5.424477 -0.4523884 -1.033613 0.2621747
AF.bobyqa 5.324676 -0.3864900 -1.010797 0.2478955
AF.Nelder_Mead 5.388769 -0.4141549 -1.022838 0.2525267
AF.nlminbwrap 5.424477 -0.4523884 -1.033613 0.2621747
AF.optimx.L-BFGS-B 5.277750 -0.3623151 -1.005134 0.2442498
AF.nloptwrap.NLOPT_LN_NELDERMEAD 5.436028 -0.4438554 -1.032669 0.2589604
AF.nloptwrap.NLOPT_LN_BOBYQA 5.471997 -0.4722368 -1.038680 0.2631687Ok, seems like a complicated problem. Thanks so much for all of your sleuthing!
As I need to submit my manuscript revisions right away, including my code and data, I think I will include some comments in the code stating something along the lines of "there appear to be compatability issues between lme4 and R 4.2.2 which are currently being worked out (link to this thread), so run this model on R 4.1.2, or use the GLMMTMB code below to run it if using R 4.2.2". Does this seem reasonable to you? As the person looking into the issue, do you have any preferred wording regarding the note about the compatibility issues?
bbolker
commented
1 year ago Does it solve/work around your problems if you use use.hessian = FALSE when you call summary() (or vcov()) ?
lukeCLarter
commented
1 year ago Yes, it does, so I will go that route. I would still like to include a note regarding why that is needed for this model in r ver 4.2.2 and link to this thread in a comment. Thats ok with you?
bbolker
commented
1 year ago Sure.
lukeCLarter
commented
1 year ago OK, great. Thanks again for all of your help!
bbolker
commented
1 year ago may I include a copy of your data in the GH repo so the example is publicly accessible if we want to work on it?
lukeCLarter
commented
1 year ago Yeah no problem, the data is publicly available elsewhere too
On Thu, Aug 17, 2023 at 1:18 PM Ben Bolker @.***> wrote:
may I include a copy of your data in the GH repo so the example is publicly accessible if we want to work on it?
— Reply to this email directly, view it on GitHub https://github.com/lme4/lme4/issues/720#issuecomment-1682753454, or unsubscribe https://github.com/notifications/unsubscribe-auth/APFYELM5WMFWGROSHJK5XXDXVZN65ANCNFSM6AAAAAAWEJNP64 . You are receiving this because you authored the thread.Message ID: @.***>
cemalley
commented
9 months ago Hi @bbolker just wanted to add to the thread for posterity that I have the same issue today (same style of glmer model, optimizer type, warnings, and discrepancy with different lme4 versions). I think it may be complicated by the update to Matrix this month. I'm sure I'm not the only one.
bbolker
commented
9 months ago Did the suggestion to specify use.hessian = FALSE in summary()/vcov() help?
cemalley
commented
9 months ago use.hessian=FALSE does not seem to solve my issue, which is "wrong" odds ratios and confidence intervals. They are really, really narrow CIs. I have a machine with older versions of lme4 and Matrix that produces "correct" ORs and CIs. I am able to run glmmTMB() with the same model formula and get correct ORs, much faster too. I saw that package mentioned in this thread. thanks.
bbolker
commented
9 months ago Can you provide data (privately/by e-mail if necessary)? While I am very unlikely to dig into the differences between old/new lme4 versions (you can see from the example how much work that can be), I do want to make sure that the newest version of lme4 gives reasonable answers — or if not, understand why not.
cemalley
commented
9 months ago Yes, I may be able to email a few examples of the data anonymized, and my code. I'll send to bbolker@gmail.com.
cemalley
commented
9 months ago Hi again, after creating a small reproducible example I can't get the issue to reproduce! It only happens with my original code which has data merging steps. It's most likely not a lme4 problem. thanks for your response regardless.
Hello, I have a manuscript under review that fit the following model for analysis:
When I initially ran the analysis, I had an older R version (ver 4.1.2), and used lme4 version 1.1-27.1, and the model with the random slope converged fine with no warnings (and indeed it converges when I use groundhog to load this older lme4 version and run it on the older version of R). However, I just updated R to the most recent version a couple of days ago and reinstalled a bunch of my packages, and with the new version of R and lme4 the model fails to converge and produces the following warnings:
These same warnings occur if I run the model with the old version of lme4 running in the new R version. The variables in question are count variables and they do not vary in scale (call complexity goes from 1-4, chorus size from 2-6). Additionally, it is the random slope which is causing the problems; when this is removed, the model converges fine in all package/version combinations.
Thus, it seesm to be the change in R versions causing the issues; the summry output for the models with old and new lme4 versions running on the new R version are the same (see below) and the same warnings are given. Do you know why these different R versions might cause this discrepancy? Do you have any advice regarding which output is more trustworthy (Coefficient estimates seem different enough between the diff version combinations to be worrisome (see below), and the model that doesnt converge has just weird values all over the place)? I'd be very grateful for any light you can shed on this. Thank you for your time.
Model summary info for the different package/version combinations are presented below:
New R, old lme4 (new R, new lme4 produces same summary):
Old R, old lme4: