Loglik of model with maxfun = 0 is not the one expected
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 at 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 at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models