[R-sig-ME] glmmTMB predicted values and AIC, vs glmer

John Maindonald john@m@|ndon@|d @end|ng |rom @nu@edu@@u
Sat May 22 05:52:39 CEST 2021 

Thanks for the informative responses.

Well, it does look as though dealing with the `Inf` is as simple as
calculating eta_predict directly from the model for the linear predictor
and returning that.  Why go the roundabout way and introduce
unnecessary inaccuracy that will in some cases manifest as a lurgy?

I’ve thought if it would be useful to collect a list of alternatives to default
control settings that users of these models have found to assist
convergence and/or avoid complaints of singularity in particular cases.
It does seem that some complaints are a consequence of the route
taken to get to convergence.  Also, there is a dearth of published work
that comments on how the various choices of link function and error
model pan out in practical data analysis contexts, and (especially in
the 99% lethal time estimates for insects on which I have worked in
a plant quarantine context in the past)response , large elements of unrealism
about the confidence that can be placed in any upper tail limits and
confidence intervals.

Incidentally, I have been experimenting with observation level random
effects because they add to the binomial variance on the scale of the
linear predictor, rather than (as for betabinomial and quasibinomial)
multiplying on the scale of the response.  For the dataset with which I
started, the variance then has to be ‘corrected’ by specifying a degree 2
or 3 polynomial or spline for the glmmTMB dispformula.


John Maindonald             email: john.maindonald using anu.edu.au<mailto:john.maindonald using anu.edu.au>


On 22/05/2021, at 13:29, Ben Bolker <bbolker using gmail.com<mailto:bbolker using gmail.com>> wrote:

 I figured out the infinite-prediction thing (didn't solve it, but figured it out):

https://github.com/glmmTMB/glmmTMB/issues/696

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 using anu.edu.au<mailto:john.maindonald using anu.edu.au> <mailto:john.maindonald using anu.edu.au>
On 22/05/2021, at 01:27, Ben Bolker <bbolker using gmail.com<mailto:bbolker using gmail.com> <mailto:bbolker using 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 using r-project.org<mailto:r-sig-mixed-models using r-project.org> <mailto:r-sig-mixed-models using r-project.org><mailto:r-sig-mixed-models using r-project.org <mailto:r-sig-mixed-models using 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 using anu.edu.au<mailto:john.maindonald using anu.edu.au> <mailto:john.maindonald using anu.edu.au><mailto:john.maindonald using anu.edu.au <mailto:john.maindonald using anu.edu.au>>
[[alternative HTML version deleted]]
_______________________________________________
R-sig-mixed-models using r-project.org<mailto:R-sig-mixed-models using r-project.org> <mailto:R-sig-mixed-models using 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 using r-project.org<mailto:R-sig-mixed-models using r-project.org> <mailto:R-sig-mixed-models using 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>


	[[alternative HTML version deleted]]

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