[R-sig-ME] nAGQ > 1 in lme4::glmer gives unexpected likelihood
Phillip Alday
ph||||p@@|d@y @end|ng |rom mp|@n|
Sun Apr 26 17:32:58 CEST 2020
The nAGQ setting doesn't affect the constant in question, but rather how
accurately the deviance is evaluated. *Very coarsely*, think of nAGQ as
impacting how much "rounding" you do at intermediate steps. (It's more
complicated than that, but the general intuition is not completely off.)
So for a lot of models, changing that setting won't have much of an
impact, but there could be some situations where that setting matters more.
The "off by a constant" issue is for the deviance and not the log
likelihood per se. While we often think of the deviance as just being -2
log likelihood, that's a very special case for the "easy" models you
learn in intro to stats (like the classical fixed-effects linear
regression). The term "deviance" actually comes the fact that it
describes the "deviation" in fit (as measured by log likelihood) from a
fully saturated model, which is relatively easy to define for classical LM.
Like so many obvious and easy things from classical regression, there is
no obvious and easy equivalent in the GLMM world. (GLM breaks things
from LM, and LMM breaks things from LM, GLMM just breaks soooo much.)
However, because we usually care about differences of deviance, this
doesn't matter too much -- the additive constant just cancels out.
Likewise, it doesn't matter for treating deviance as the objective in an
optimization problem.
This is then the problem with different types of "deviance". Like the
log likelihood, the deviance is uniquely defined for a given probability
model (if we only know what the saturated model is). So we can talk
about the deviance of an unconditional or conditional model, but it's a
little bit problematic to talk about the "unconditional deviance",
although I understand the temptation to use that as a convenient shorthand.
I've simplified a bit and I'm hoping that I've cleared up more
inaccuracies than I've created, lest the stats deities smite me ....
Phillip
On 25/04/2020 16:41, Vaida, Florin wrote:
> Looks like you guys thought about deviance issues a lot:
> https://github.com/lme4/lme4/issues/375
>
> Speaking about the principle of "least surprise", it is surprising to get completely different behaviors in deviance for nAGQ=1 vs nAGQ=2 (say), in the same model - (what Ben Goldstein pointed out).
>
> A couple further questions:
> 1. Is the "off by" constant the same, for a given model, if I only change nAGQ? E.g., comparing nAGQ=5 and nAGQ=20?
> 2. If so, can you calculate the "off by" constant by comparing with formula-based Laplace approximation with a bona fide Adaptive Gaussian Quadrature using existing algorithm but n=1?
>
> Florin
>
> On Apr 24, 2020, at 7:59 AM, Ben Bolker <bbolker using gmail.com<mailto:bbolker using gmail.com>> wrote:
>
>
> I found the note I was looking for. In ?deviance.glmerMod it says:
>
> If adaptive Gauss-Hermite quadrature is used, then
> ‘logLik(object)’ is currently only proportional to the
> absolute-unconditional log-likelihood.
>
> (see the discussion on this page for more context); see also https://github.com/lme4/lme4/blob/master/misc/logLikGLMM/logLikGLMM.R
>
>
> On 4/24/20 10:24 AM, Douglas Bates wrote:
> Having said that, I do see that the fits in the MixedModels package for
> Julia produce similar values of the deviance with the Laplace approximation
> and nAGQ = 7
>
> julia> m1 = fit(MixedModel, @formula(y ~ 1 + (1|group)), dd, Poisson())
> Generalized Linear Mixed Model fit by maximum likelihood (nAGQ = 1)
> y ~ 1 + (1 | group)
> Distribution: Poisson{Float64}
> Link: LogLink()
>
> Deviance: 193.5587
>
> Variance components:
> Column Variance Std.Dev.
> group (Intercept) 3.9577026 1.9893975
>
> Number of obs: 100; levels of grouping factors: 10
>
> Fixed-effects parameters:
> ──────────────────────────────────────────────────
> Estimate Std.Error z value P(>|z|)
> ──────────────────────────────────────────────────
> (Intercept) 2.65175 0.632317 4.19 <1e-4
> ──────────────────────────────────────────────────
>
> julia> m1 = fit(MixedModel, @formula(y ~ 1 + (1|group)), dd, Poisson(),
> nAGQ=7)
> Generalized Linear Mixed Model fit by maximum likelihood (nAGQ = 7)
> y ~ 1 + (1 | group)
> Distribution: Poisson{Float64}
> Link: LogLink()
>
> Deviance: 193.5104
>
> Variance components:
> Column Variance Std.Dev.
> group (Intercept) 3.9577026 1.9893975
>
> Number of obs: 100; levels of grouping factors: 10
>
> Fixed-effects parameters:
> ──────────────────────────────────────────────────
> Estimate Std.Error z value P(>|z|)
> ──────────────────────────────────────────────────
> (Intercept) 2.65175 0.632317 4.19 <1e-4
> ──────────────────────────────────────────────────
>
> As the person who wrote the first version of the nAGQ code in R I would not
> be surprised if there was a constant dropped somewhere. It is difficult
> code.
>
> And the results here in the Julia package make me uncomfortable because the
> values of the parameter estimates are identical in the two fits. I would
> expect them to be close but not identical.
>
> Isn't it good to know that there is still room for research in this area?
> :-)
>
> On Fri, Apr 24, 2020 at 9:05 AM Douglas Bates <bates using stat.wisc.edu<mailto:bates using stat.wisc.edu>> wrote:
>
> There's a lot of variability in your lambdas
>
> exp(3 + random_effects)
> [1] 91.5919358 6.9678749 4.1841478 78.0771666 890.6931394
> 20.8558107
> [7] 3.0037864 0.3049416 2.1675995 40.6209684
>
> Do you really expect that some groups will have a mean count of nearly 900
> whereas others will have a mean count less than 1?
>
>
> On Wed, Apr 22, 2020 at 5:58 PM Ben Goldstein <ben.goldstein using berkeley.edu<mailto:ben.goldstein using berkeley.edu>>
> wrote:
>
> Hi all,
>
> I'm using lme4::glmer to estimate Poisson mixed models in a very simple
> context (single random effect). I'm interested in the model
> likelihood/AIC
> across many simulated datasets.
>
> To investigate whether the Laplace approximation was appropriate for my
> data context, I explored using the argument nAGQ to improve the accuracy
> of
> the likelihood estimation. When I changed nAGQ to a value > 1, I saw an
> unexpectedly huge change in the likelihood; log-likelihoods tended to be
> off by ~200. Other statistics packages (e.g. GLMMadaptive) yield
> estimates
> that agree with lme4's Laplace approximation, as did a manual likelihood
> estimate, and not with the nAGQ > 2 estimate.
>
> The following code reproduces the problem I'm encountering.
>
> *# r-sig-mixed-models GLMM question*
> library(lme4)
> set.seed(51)
>
> *# Simulate some random effect-driven Poisson data*
> random_effects <- rnorm(10, 0, 2)
> group <- rep(1:10, 10)
> simulated_data <- data.frame(y = rpois(n = 100, lambda = exp(3 +
> random_effects[group])),
> group = group)
>
> *# Fit models with Laplace (nAGQ = 1) and nAGQ = 11*
> fit_Laplace <- glmer(y ~ (1|group), data = simulated_data, family =
> poisson())
> fit_AGQ <- glmer(y ~ (1|group), data = simulated_data, family =
> poisson(),
> nAGQ = 11)
>
> logLik(fit_Laplace)
> logLik(fit_AGQ)
> logLik(fit_Laplace) - logLik(fit_AGQ) *# Huge difference!*
>
> When I execute the above code, I see a difference in likelihood of
> -218.8894. I've tested across many simulations and on 2 different
> machines
> (Mac and Linux). My version of lme4 is up to date.
>
> Has anyone run into this issue before? Am I using the glmer function
> wrong,
> or is it possible there's something going on under the hood?
>
> Thanks,
> Ben
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> 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
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> 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
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> 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
>
> _______________________________________________
> 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
>
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-mixed-models using r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
More information about the R-sig-mixed-models
mailing list