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

Marc Girondot m@rc_grt @end|ng |rom y@hoo@|r
Mon Sep 28 22:25:18 CEST 2020


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)



More information about the R-sig-mixed-models mailing list