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

[Phillip Alday]

Potentially. In my experience, the best-tested portions of GLMM software
(including lme4 and MixedModels) are related to Binomial models, then
Poisson, then any model family with a dispersion parameter.

That said, families without a dispersion parameter (Binomial and
Poisson) tend to work equally well in terms of the code, though Poisson
might be less well-behaved numerically.

On 25/04/2020 16:40, Vaida, Florin wrote:
> 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
> _______________________________________________
> R-sig-mixed-models using r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models