¿Por qué glmer no consigue la máxima verosimilitud (como se comprueba al aplicar una optimización genérica adicional)?

Question

Derivar numéricamente el MLE s de GLMM es difícil y, en la práctica, sé que no debemos utilizar la optimización por fuerza bruta (por ejemplo, utilizando optim de forma sencilla). Pero por mi propio propósito educativo, quiero probarlo para asegurarme de que entiendo correctamente el modelo (ver el código de abajo). Encontré que siempre obtengo resultados inconsistentes de glmer() .

En particular, incluso si utilizo los MLEs de glmer como valores iniciales, según la función de probabilidad que escribí ( negloglik ), no son MLEs ( opt1$value es menor que opt2 ). Creo que hay dos posibles razones:

1. negloglik no está bien escrito, por lo que tiene demasiados errores numéricos, y
2. la especificación del modelo es errónea. Para la especificación del modelo, el modelo previsto es:

\begin{equation} L=\prod_{i=1}^{n} \left(\int_{-\infty}^{\infty}f(y_i|N,a,b,r_{i})g(r_{i}|s)dr_{i}\right) \end{equation} donde $f$ es una pmf binomial y $g$ es un pdf normal. Estoy tratando de estimar $a$ , $b$ y $s$ . En particular, quiero saber si la especificación del modelo es incorrecta, cuál es la especificación correcta.

p <- function(x,a,b) exp(a+b*x)/(1+exp(a+b*x))

a <- -4  # fixed effect (intercept)
b <- 1   # fixed effect (slope)
s <- 1.5 # random effect (intercept)
N <- 8
x <- rep(2:6, each=20)
n <- length(x) 
id <- 1:n
r  <- rnorm(n, 0, s) 
y  <- rbinom(n, N, prob=p(x,a+r,b))

negloglik <- function(p, x, y, N){
  a <- p[1]
  b <- p[2]
  s <- p[3]

  Q <- 100  # Inf does not work well
  L_i <- function(r,x,y){
    dbinom(y, size=N, prob=p(x, a+r, b))*dnorm(r, 0, s)
  }

  -sum(log(apply(cbind(y,x), 1, function(x){ 
    integrate(L_i,lower=-Q,upper=Q,x=x[2],y=x[1],rel.tol=1e-14)$value
  })))
}

library(lme4)
(model <- glmer(cbind(y,N-y)~x+(1|id),family=binomial))

opt0 <- optim(c(fixef(model), sqrt(VarCorr(model)$id[1])), negloglik, 
                x=x, y=y, N=N, control=list(reltol=1e-50,maxit=10000)) 
opt1 <- negloglik(c(fixef(model), sqrt(VarCorr(model)$id[1])), x=x, y=y, N=N)
opt0$value  # negative loglikelihood from optim
opt1        # negative loglikelihood using glmer generated parameters
-logLik(model)==opt1 # but these are substantially different...

Un ejemplo más sencillo

Para reducir la posibilidad de tener un gran error numérico, he creado un ejemplo más sencillo.

y  <- c(0, 3)
N  <- c(8, 8)
id <- 1:length(y)

negloglik <- function(p, y, N){
  a <- p[1]
  s <- p[2]
  Q <- 100  # Inf does not work well
  L_i <- function(r,y){
    dbinom(y, size=N, prob=exp(a+r)/(1+exp(a+r)))*dnorm(r,0,s)
  }
  -sum(log(sapply(y, function(x){
    integrate(L_i,lower=-Q, upper=Q, y=x, rel.tol=1e-14)$value
  })))
}

library(lme4)
(model <- glmer(cbind(y,N-y)~1+(1|id), family=binomial))
MLE.glmer <- c(fixef(model), sqrt(VarCorr(model)$id[1]))
opt0 <- optim(MLE.glmer, negloglik, y=y, N=N, control=list(reltol=1e-50,maxit=10000)) 
MLE.optim <- opt0$par
MLE.glmer # MLEs from glmer
MLE.optim # MLEs from optim

L_i <- function(r,y,N,a,s) dbinom(y,size=N,prob=exp(a+r)/(1+exp(a+r)))*dnorm(r,0,s)

L1 <- integrate(L_i,lower=-100, upper=100, y=y[1], N=N[1], a=MLE.glmer[1], 
                s=MLE.glmer[2], rel.tol=1e-10)$value
L2 <- integrate(L_i, lower=-100, upper=100, y=y[2], N=N[2], a=MLE.glmer[1], 
                s=MLE.glmer[2], rel.tol=1e-10)$value

(log(L1)+log(L2)) # loglikelihood (manual computation)
logLik(model)     # loglikelihood from glmer

Answer 1

Fijar un valor alto de nAGQ en el glmer hizo que los MLEs de los dos métodos fueran equivalentes. La precisión por defecto de glmer no era muy bueno. Esto resuelve la cuestión.

glmer(cbind(y,N-y)~1+(1|id),family=binomial,nAGQ=20)

Ver la respuesta de @SteveWalker aquí ¿Por qué no puedo hacer coincidir la salida de glmer (familia=binomial) con la implementación manual del algoritmo de Gauss-Newton? para más detalles.