[R-sig-ME] Loglik of model with maxfun = 0 is not the one expected
Ben Bolker
bbo|ker @end|ng |rom gm@||@com
Mon Sep 28 22:54:29 CEST 2020
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
More information about the R-sig-mixed-models
mailing list