glmmTMB predicted values and AIC, vs glmer
What allFit indicates to me is that all of the different optimizers
reach very similar log-likelihoods (max difference of 1-4), and if I
recall correctly all of the fixed effects are also very similar, so I
would feel free to disregard the convergence warning as a false positive.
The AIC values are still a lot more different than I would like, and
digging into the gory details (available on request) it seems that they
are actually giving different likelihood calculations -- e.g. for the
parameters reached by glmer, glmmTMB returns a deviance of 607.037 while
glmer returns 609.3362. I suspect something deep in the guts of e.g.
the way the Laplace approximation is computed (not technically a bug in
either code, but some kind of numerical instability).
This is good news if what you're interested in is making sure you get
the correct parameter values. It's bad news if you're interested in
mixing and matching glmmTMB and lme4 fits in a model comparison.
I'm not horribly surprised that things are a bit numerically unstable
- 12 parameters for 95 observations, and the extremely rapid tail
behaviour of cloglog often makes things worse.
I do intend to track down where the infinite prediction is coming
from too, but haven't gotten there yet.
I'm going to include some of my code in case it's useful, but it's
poorly commented, messy *and may be in the wrong order* - I was bouncing
around developing it interactively.
====
Obsglmer1.cll <- glmer(form1, nAGQ=1,
family=binomial(link="cloglog"), data=ffly)
aa <- allFit(Obsglmer1.cll)
summary(aa)$fixef
ll <- -summary(aa)$llik
ll-min(ll)
dotwhisker::dwplot(modList[-2])
library(dotwhisker)
modList <- list(glmmTMB=ObsTMB.cll, nAGQ0=Obsglmer.cll,
nAGQ1=Obsglmer1.cll)
f1 <- ObsTMB.cll$obj$fn
## TMB params, TMB format
p1 <- with(environment(f1),last.par.best[-random])
## TMB params, glmer format
p1B <- c(exp(p1[names(p1)=="theta"]), p1[names(p1)=="beta"])
f3 <- getME(Obsglmer1.cll,"devfun")
## glmer params, glmer format
p3 <- unlist(getME(Obsglmer1.cll,c("theta","beta")))
## glmer params, TMB format
p3B <- c(p3[-(1:2)],log(p3[1:2]))
f3(p1B) ## 609.3559
f3(p3) ## 609.3362
2*f1(p3B) ## 608.037
2*f1(p1) ## 608.0166
Obsglmer2.cll <- update(Obsglmer1.cll,
start=list(theta=p1B[1:2],
fixef=p1B[-(1:2)]))
dwplot(modList)
library(bbmle)
AICtab(modList, logLik=TRUE)
modList2 <- c(modList, list(nAGQ1_newstart=Obsglmer2.cll))
AICtab(modList2, logLik=TRUE)
On 5/21/21 6:24 PM, John Maindonald wrote:
Setting nAGQ=1 does bring the AIC values closer together, albeit with convergence failure when other parameters are left at their defaults. Running allFit() does not provide (to me) any obvious clues on what might work better. ?Or, should one just accept the less strict tolerance? Nor does centering and scale scTime help in this case. llik <- summary(allFit(Obsglmer1.cll))$llik max(lik)-lik ## ? ? ? ? ? ? ? ? ? ? ? bobyqa ? ? ? ? ? ? ? ? ? Nelder_Mead ## ? ? ? ? ? 0.00000000000e+00 ? ? ? ? ? ? 1.31201379190e-04 ## ? ? ? ? ? ? ? ? ? nlminbwrap ? ? ? ? ? ? ? ? ? ? ? ? nmkbw ## ? ? ? ? ? ?1.76015646502e-09 ? ? ? ? ? ? 3.48215053236e-06 ## ? ? ? ? ? ? ?optimx.L-BFGS-B nloptwrap.NLOPT_LN_NELDERMEAD ## ? ? ? ? ? ?1.09246872171e-06 ? ? ? ? ? ? 1.13990722639e-07 ## ? ?nloptwrap.NLOPT_LN_BOBYQA ## ? ? ? ? ? ?5.81495783081e-07 John Maindonaldemail: john.maindonald at anu.edu.au <mailto:john.maindonald at anu.edu.au>
On 22/05/2021, at 01:27, Ben Bolker <bbolker at gmail.com <mailto:bbolker at gmail.com>> wrote: ?I don't know yet, I will dig in and try to see what's going on. ?An infinite predicted value certainly seems like an issue. ?The most obvious difference is that nAGQ=0 is actually doing something slightly different (it's fitting the fixed-effect parameters as part of the PIRLS loop rather than maximizing over them explicitly); I would rather compare the nAGQ=1 case, just to minimize the number of differences. On 5/21/21 6:30 AM, John Maindonald wrote:
The code that follows below demonstrates what I find a very odd
issue with r-sig-mixed-models at r-project.org
<mailto:r-sig-mixed-models at r-project.org><mailto:r-sig-mixed-models at r-project.org
<mailto:r-sig-mixed-models at r-project.org>>, and with the comparison
between
glmmTMB::glmmTMB and lme4::glmer, for models that have the
same model formulae.
1) In the glmmTMB model as fitted to `ffly`, one predicted value is
infinite. ?This, even though all coefficients and standard errors
have finite values. ?This surely indicates that there is something
odd in the way that predictions on the scale of the linear
predictor are calculated.
2) Initially, I'd thought that the `Inf` value explained the
difference between the two AIC values. ?But the difference
pretty much stays the same when the "offending" point is
removed.
I've not checked what happens when there are no observation level
random effects. ?Is there some difference in the way that the
model formulae are interpreted in the two cases?
ffly <-
read.csv('https://github.com/jhmaindonald/dataR/files/6521483/ffly.csv <https://github.com/jhmaindonald/dataR/files/6521483/ffly.csv>')
ffly$obs <- factor(ffly$obs)
form1 <- cbind(Dead,Live)~0+trtGp/TrtTime+(1|obs)+(1|trtGpRep)
library(lme4); library(glmmTMB)
ObsTMB.cll <- glmmTMB(form1,
??????????????????????family=binomial(link="cloglog"), data=ffly)
Obsglmer.cll <- glmer(form1, nAGQ=0,
??????????????????????family=binomial(link="cloglog"), data=ffly)
round(AIC(Obsglmer.cll, ObsTMB.cll), 2)
## ?????????????????????df ???AIC
## Obsglmer.cll 14 639.44
## ObsTMB.cll ?14 636.02
round(c(max(predict(ObsTMB.cll)),max(predict(Obsglmer.cll))),3)
## [1] ??Inf 3.833
range(fixef(ObsTMB.cll))
## [1] -5.28780064553 ?1.15952249272
range(vcov(ObsTMB.cll))
## [1] -0.0546986905338 ?0.2887049942215
## Try also; there are small deviations from the line
plot(predict(Obsglmer.cll)~predict(ObsTMB.cll))
abline(0,1)
## Remove offending point
ffly1 <- ffly[-which.max(predict(ObsTMB.cll)),]
ObsTMB1.cll <- glmmTMB(form1,
??????????????????????family=binomial(link="cloglog"), data=ffly1)
Obsglmer1.cll <- glmer(form1, nAGQ=0,
??????????????????????family=binomial(link="cloglog"), data=ffly1)
cbind(AIC(Obsglmer.cll, ObsTMB.cll), AIC1=AIC(Obsglmer1.cll,
ObsTMB1.cll)[,2])
## ??????????????????????df ???????????????????AIC ?????????????????AIC1
## Obsglmer.cll ?14 639.441888969 639.441889196
## ObsTMB.cll ??14 636.016597730 636.016597723
??## Observe that changes are in the final four decimal places.
round(rbind(glmer=fixef(Obsglmer1.cll),
TMB=unclass(fixef(ObsTMB1.cll)$cond)),3)
??????trtGpspAEgg trtGpspAL2 trtGpspAL3 trtGpspBEgg trtGpspBL2 trtGpspBL3
glmer ??????0.772 ????-2.255 ????-3.299 ?????-1.803 ????-2.632 ????-5.114
TMB ????????0.790 ????-2.367 ????-3.441 ?????-1.889 ????-2.757 ????-5.288
??????trtGpspAEgg:TrtTime trtGpspAL2:TrtTime trtGpspAL3:TrtTime
glmer ??????????????0.231 ?????????????0.372 ?????????????0.563
TMB ????????????????0.288 ?????????????0.398 ?????????????0.602
??????trtGpspBEgg:TrtTime trtGpspBL2:TrtTime trtGpspBL3:TrtTime
glmer ??????????????0.278 ?????????????0.517 ?????????????1.112
TMB ????????????????0.299 ?????????????0.554 ?????????????1.160
John Maindonald email: john.maindonald at anu.edu.au
<mailto:john.maindonald at anu.edu.au><mailto:john.maindonald at anu.edu.au
<mailto:john.maindonald at anu.edu.au>>
[[alternative HTML version deleted]]
_______________________________________________ R-sig-mixed-models at r-project.org <mailto:R-sig-mixed-models at r-project.org> mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models <https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models>
_______________________________________________ R-sig-mixed-models at r-project.org <mailto:R-sig-mixed-models at r-project.org> mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models <https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models>