Pourquoi glmer n'atteint-il pas le maximum de vraisemblance (comme le prouve l'application d'une optimisation générique supplémentaire)?

37

Dériver numériquement les MLE de GLMM est difficile et, dans la pratique, je sais que nous ne devrions pas utiliser l'optimisation de la force brute (par exemple, en utilisant optimune méthode simple). Mais pour mon propre but éducatif, je veux l'essayer pour m'assurer de bien comprendre le modèle (voir le code ci-dessous). J'ai trouvé que je reçois toujours des résultats incohérents glmer().

En particulier, même si j'utilise les MLE glmercomme valeurs initiales, selon la fonction de vraisemblance que j'ai écrite ( negloglik), ce ne sont pas des MLE ( opt1$valueest plus petit que opt2). Je pense que deux raisons potentielles sont:

  1. negloglik n'est pas bien écrit pour qu'il y ait trop d'erreur numérique, et
  2. la spécification du modèle est fausse. Pour la spécification du modèle, le modèle prévu est:

L=i=1n(f(yi|N,a,b,ri)g(ri|s)dri)
où est un binôme pmf et est un pdf normal. J'essaie d'estimer , et . En particulier, je veux savoir si la spécification du modèle est erronée, quelle est la spécification correcte.fgabs
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 exemple plus simple

Pour réduire le risque d'erreur d'erreur importante, j'ai créé un exemple plus simple.

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 
ergoter
la source
Les MLE (et non les log-vraisemblances elles-mêmes) sont-ils comparables? C'est-à-dire, êtes-vous en train de perdre une constante?
Ben Bolker
1
Les MLE estimés sont clairement différents ( MLE.glmeret MLE.optim) en particulier pour l'effet aléatoire (voir le nouvel exemple), de sorte que ce n'est pas uniquement basé sur un facteur constant dans les valeurs de vraisemblance, je pense.
chicaner
4
@Ben Définir une valeur élevée de nAGQdans glmerrend les MLE comparables. La précision par défaut de glmern'était pas très bonne.
chicaner
5
Lien vers une question similaire lme4 que @Steve Walker m'a aidée à résoudre
Ben Ogorek
3
Comme une question plus ancienne avec beaucoup de votes positifs, cela pourrait probablement être préservé. Je ne vois pas la nécessité de fermer cela.
Gay - Rétablir Monica

Réponses:

3

Définir une valeur élevée nAGQdans l' glmerappel rendait les MLE des deux méthodes équivalentes. La précision par défaut de glmern'était pas très bonne. Cela règle la question.

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

Voir la réponse de @ SteveWalker ici. Pourquoi ne puis-je pas associer la sortie glmer (famille = binomial) à la mise en œuvre manuelle de l'algorithme de Gauss-Newton? pour plus de détails.

ergoter
la source
1
Mais les probabilités loglik estimées sont très différentes (probablement par une constante), de sorte que les différentes méthodes ne doivent pas être mélangées.
chicaner
1
hmm, intéressant / surprenant - merci d’avoir mis en place cet exemple, je vais essayer de trouver le temps de le regarder.
Ben Bolker