Por que o glmer não atinge a máxima probabilidade (conforme verificado pela aplicação de otimização genérica adicional)?

37

Numericamente derivar a MLE s de GLMM é difícil e, na prática, eu sei, não devemos usar a otimização de força bruta (por exemplo, usando optimem uma maneira simples). Mas, para meu próprio objetivo educacional, quero experimentá-lo para garantir a compreensão correta do modelo (veja o código abaixo). Descobri que sempre obtenho resultados inconsistentes glmer().

Em particular, mesmo que eu use os MLEs glmercomo valores iniciais, de acordo com a função de probabilidade que escrevi ( negloglik), eles não são MLEs ( opt1$valueé menor que opt2). Eu acho que duas razões possíveis são:

  1. negloglik não está bem escrito para que exista muito erro numérico e
  2. a especificação do modelo está incorreta. Para a especificação do modelo, o modelo pretendido é:

L=i=1n(f(yi|N,a,b,ri)g(ri|s)dri)
que é um pmf binomial é um pdf normal. Eu estou tentando estimar , , e . Em particular, quero saber se a especificação do modelo está errada, qual é a especificação correta.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...

Um exemplo mais simples

Para reduzir a possibilidade de ocorrer um grande erro numérico, criei um exemplo mais simples.

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 
queixa
fonte
os MLEs (não as próprias probabilidades logarítmicas) são comparáveis? Ou seja, você está apenas de folga por uma constante?
Ben Bolker 27/05
11
As MLEs estimadas são claramente diferentes ( MLE.glmere MLE.optim) especialmente para o efeito aleatório (veja o novo exemplo), portanto, não é apenas baseado em algum fator constante nos valores de probabilidade, eu acho.
Quibble 27/05
4
@ Ben A definição de um valor elevado de nAGQno glmerfeitas as mles comparável. A precisão padrão de glmernão era muito boa.
Quibble 28/05
5
Vinculando a uma pergunta lme4 similar que @ Steve Walker me ajudou com: stats.stackexchange.com/questions/77313/...
Ben Ogórek
3
Como uma pergunta antiga com muitos votos positivos, isso provavelmente poderia ser adquirido. Não vejo necessidade disso ser fechado.
gung - Restabelece Monica

Respostas:

3

Definir um valor alto de nAGQna glmerchamada fez os MLEs dos dois métodos equivalentes. A precisão padrão de glmernão era muito boa. Isso resolve o problema.

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

Veja a resposta de SteveWalker aqui. Por que não consigo combinar a saída glmer (family = binomial) com a implementação manual do algoritmo de Gauss-Newton? para mais detalhes.

queixa
fonte
11
Mas as probabilidades logísticas estimadas são muito diferentes (presumivelmente por alguma constante); portanto, os diferentes métodos não devem ser misturados.
Quibble 28/05
11
hmm, interessante / surpreendente - obrigado por configurar este exemplo, vou tentar encontrar tempo para analisá-lo.
Ben Bolker 28/05