[R-sig-ME] Loglik of model with maxfun = 0 is not the one expected

[Ben Bolker]

   This is a perfectly reasonable, and an unexpectedly difficult, 
problem.  The glmer-fitting machinery has a lot of internal state that 
is hard to set exactly right.

   In this particular case, adding 
nAGQ0initStep=FALSE to glmerControl() works [it disables the first-stage 
optimization step; presumably at some point during this step, or the 
pre- and post-processing surrounding it, the GLMM likelihood evaluation 
is called and messes up the internal state)

Here is another way to reproduce the log-likelihood:

m3 <- update(m1, devFunOnly=TRUE)
-0.5*m3(c(theta,beta))

The development version of lme4 will also allow

-0.5*getME(m1,"devfun")(c(theta,beta))

(the internal code is more efficient and self-contained than the 
update(., "devFunOnly") version)

On 9/28/20 4:25 PM, Marc Girondot via R-sig-mixed-models wrote:
> Thanks a lot to Thierry Onkelinx <thierry.onkelinx using inbo.be> to answer to 
> my previous demand. It works now !
> 
> Another point that I don't understand.
> 
> I fit a model. I get -log likelihood
> I extract the fitted parameters.
> I inject these fitted parameters to estimate the likelihood of the model 
> and I obtain a different -log likelihood.
> 
> So I suppose that I don't extract all the fitted parameters... but which 
> one ?
> 
> Thanks if you can help me to understand what's happened here.
> 
> Marc
> 
> ___________________
> 
> library(lme4)
> 
> # Ei is simply data that will be used to fit the model
> 
> Ei <- structure(list(Ma = c(15, 14, 28, 5, 26, 7, 25, 4, 9, 25,
>                           8, 11, 2, 3, 1, 1, 0, 0, 10, 5, 0, 0, 9, 1, 2, 
> 0, 10, 0, 0, 0,
>                           8, 1, 0, 0, 12, 3, 29, 10, 26, 9, 18, 5, 14, 
> 6, 1, 2, 0, 0),
>                 Fa = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 5, 0, 4, 5, 18, 2,
>                             11, 2, 10, 0, 3, 7, 3, 0, 9, 8, 8, 0, 10, 9, 
> 9, 0, 8, 5,
>                             4, 0, 0, 0, 0, 5, 1, 2, 0, 8, 0, 8, 3, 10, 3),
>                 IT = c(26.4,
>                                            26.9, 27.9, 27.9, 28.4, 28.4, 
> 29, 29, 29, 29.5, 29.5, 29.5,
>                                            30, 30, 30.3, 30.3, 30.8, 
> 30.8, 28, 29.5, 31, 32.5, 28, 29.5,
>                                            31, 32.5, 28, 29.5, 31, 32.5, 
> 28, 29.5, 31, 32.5, 28.15,
>                                            28.15, 28.65, 28.65, 29.15, 
> 29.15, 29.55, 29.55, 29.75, 29.75,
>                                            30.05, 30.05, 30.65, 30.65),
>                 RMU = structure(c(2L, 2L, 2L,
>                                   2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
> 2L, 2L, 2L, 2L, 2L, 2L,
>                                   3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
> 3L, 3L, 3L, 3L, 3L, 3L,
>                                   3L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
> 1L, 1L, 1L, 1L, 1L, 1L),
>                                 .Label = c("ASW", "AW",
>                                            "PSW"), class = "factor"),
>                 C = structure(c(4L,
>                                      4L, 4L, 3L, 4L, 3L, 4L, 7L, 8L, 4L, 
> 7L, 8L, 7L, 8L, 7L, 8L,
>                                      7L, 8L, 5L, 5L, 5L, 5L, 6L, 6L, 6L, 
> 6L, 10L, 10L, 10L, 10L,
>                                      9L, 9L, 9L, 9L, 1L, 2L, 1L, 2L, 1L, 
> 2L, 1L, 2L, 1L, 2L, 1L,
>                                      2L, 1L, 2L),
>                                    .Label = c("A", "B", "C", "D",
>                                               "E", "F", "G", "H", "I", 
> "J"), class = "factor")),
>            class = "data.frame", row.names = 302:349)
> 
> # m1 if the fitted model
> 
> m1 <- glmer(formula = cbind(Ma, Fa) ~ 0 + IT + (1 | RMU / C) + (0 + IT | 
> RMU),
>              data=Ei,
>              control = glmerControl(optimizer="bobyqa", optCtrl = 
> list(maxfun = 1000000)),
>              family = binomial(link = "logit"))
> 
> # I check that all is ok.
> logLik(m1)
> ranef(m1)
> fixef(m1)
> 
> theta <- getME(m1, "theta")
> beta <- getME(m1, "beta")
> 
> # use the fitted values of m1 that I introduce in m2. Note that the 
> model is not fitted because maxfun = 0. Log likelihood is just calculated
> m2 <- glmer(formula = cbind(Ma, Fa) ~ 0 + IT + (1 | RMU / C) + (0 + IT | 
> RMU),
>              data=Ei,
>              control = glmerControl(optimizer="bobyqa", optCtrl = 
> list(maxfun = 0)),
>              start = list(theta = theta,
>                           fixef = beta),
>              family = binomial(link = "logit"))
> 
> # Both likelihood should be identical... they are not
> logLik(m1)
> logLik(m2)
> 
> _______________________________________________
> R-sig-mixed-models using r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models