[R-sig-ME] nAGQ > 1 in lme4::glmer gives unexpected likelihood

[Vaida, Florin]

Interestingly (and reassuringly), Laplace and nAGQ give consistent results for binomial GLMM.
Could there be a glitch for Poisson GLMM?

#### Likelihood of Binomial GLMM Models

# r-sig-mixed-models GLMM question*
library(lme4)
set.seed(51)

# Simulate some random effect-driven Poisson data*
# random_effects <- rnorm(10, 0, 2)
# random_effects <- rnorm(10, 0, 0.01)
random_effects <- rnorm(10, 0, 1)

group <- rep(1:10, 10)
simulated_data <- data.frame(y = rbinom(n = 100, size = 1, 
    prob = exp(2+random_effects[group])/(1+exp(2+random_effects[group]))), group = group)

# Fit models with Laplace (nAGQ = 1) and nAGQ = 11*
fit_Laplace <- glmer(y ~ (1|group), data = simulated_data, family =
                       binomial())
fit_AGQ <- glmer(y ~ (1|group), data = simulated_data, family = binomial(),
                 nAGQ = 11)
# fit_glm = glm(y ~ group, data=simulated_data, family=binomial)

(logLik(fit_Laplace))
(logLik(fit_AGQ))  # very similar to fit_Laplace
# (logLik(fit_glm)) 


> On Apr 24, 2020, at 7:24 AM, Douglas Bates <bates using stat.wisc.edu> 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> 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>
>> 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 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
>> 
> 
> 	[[alternative HTML version deleted]]
> 
> _______________________________________________
> R-sig-mixed-models using r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models