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 optim
une 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 glmer
comme valeurs initiales, selon la fonction de vraisemblance que j'ai écrite ( negloglik
), ce ne sont pas des MLE ( opt1$value
est plus petit que opt2
). Je pense que deux raisons potentielles sont:
negloglik
n'est pas bien écrit pour qu'il y ait trop d'erreur numérique, et- la spécification du modèle est fausse. Pour la spécification du modèle, le modèle prévu est:
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.
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
MLE.glmer
etMLE.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.nAGQ
dansglmer
rend les MLE comparables. La précision par défaut deglmer
n'était pas très bonne.Réponses:
Définir une valeur élevée
nAGQ
dans l'glmer
appel rendait les MLE des deux méthodes équivalentes. La précision par défaut deglmer
n'était pas très bonne. Cela règle la question.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.
la source