Dlaczego glmer nie osiąga maksymalnego prawdopodobieństwa (potwierdzonego przez dalszą optymalizację ogólną)?

37

Wyprowadzanie liczbowe MLE z GLMM jest trudne i, w praktyce, wiem, nie powinniśmy stosować optymalizacji siły brutalnej (np. Używając optimw prosty sposób). Ale dla własnego celu edukacyjnego chcę go wypróbować, aby upewnić się, że poprawnie rozumiem model (patrz poniższy kod). Odkryłem, że zawsze otrzymuję niespójne wyniki glmer().

W szczególności, nawet jeśli użyję MLE z glmerwartości początkowych, zgodnie z funkcją prawdopodobieństwa, którą napisałem ( negloglik), nie są to MLE ( opt1$valuejest mniejsze niż opt2). Myślę, że dwa potencjalne powody to:

  1. negloglik nie jest dobrze napisany, więc występuje w nim zbyt duży błąd numeryczny, oraz
  2. specyfikacja modelu jest nieprawidłowa. Dla specyfikacji modelu zamierzonym modelem jest:

L=i=1n(f(yi|N,a,b,ri)g(ri|s)dri)
gdzie jest dwumianowym pmf jest normalnym pdf. Próbuję oszacować , i . W szczególności chcę wiedzieć, czy specyfikacja modelu jest nieprawidłowa, jaka jest poprawna specyfikacja.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...

Prostszy przykład

Aby zmniejszyć możliwość wystąpienia dużego błędu numerycznego, stworzyłem prostszy przykład.

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 
spierać się
źródło
czy MLE (a nie same log-prawdopodobieństwo) są porównywalne? To znaczy, czy jesteś po prostu stały?
Ben Bolker
Oszacowane MLE są wyraźnie różne ( MLE.glmeri MLE.optim) szczególnie w przypadku efektu losowego (patrz nowy przykład), więc myślę, że nie opiera się on tylko na pewnym stałym współczynniku wartości prawdopodobieństwa.
spór
4
@Ben Ustawienie wysokiej wartości nAGQin glmersprawiło, że MLE są porównywalne. Domyślna precyzja glmernie była zbyt dobra.
quibble
5
Link do podobnego pytania lme4, które pomogło mi @Steve Walker: stats.stackexchange.com/questions/77313/...
Ben
3
Jako starsze pytanie z dużą liczbą pozytywnych opinii, prawdopodobnie można by to uznać za prawomocne. Nie widzę potrzeby, aby to było zamknięte.
gung - Przywróć Monikę

Odpowiedzi:

3

Ustawienie wysokiej wartości nAGQw glmerwywołaniu sprawiło, że MLE z dwóch metod są równoważne. Domyślna precyzja glmernie była zbyt dobra. To rozwiązuje problem.

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

Zobacz odpowiedź @ SteveWalkera tutaj Dlaczego nie mogę dopasować wyjścia glmer (rodzina = dwumianowa) do ręcznej implementacji algorytmu Gaussa-Newtona? po więcej szczegółów.

spierać się
źródło
1
Jednak szacunkowe prawdopodobieństwa logiczne są bardzo różne (przypuszczalnie pewne stałe), więc nie należy mieszać różnych metod.
quibble
hmm, ciekawe / zaskakujące - dziękuję za skonfigurowanie tego przykładu, postaram się znaleźć czas, aby się temu przyjrzeć.
Ben Bolker