Wyprowadzanie liczbowe MLE z GLMM jest trudne i, w praktyce, wiem, nie powinniśmy stosować optymalizacji siły brutalnej (np. Używając optim
w 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 glmer
wartości początkowych, zgodnie z funkcją prawdopodobieństwa, którą napisałem ( negloglik
), nie są to MLE ( opt1$value
jest mniejsze niż opt2
). Myślę, że dwa potencjalne powody to:
negloglik
nie jest dobrze napisany, więc występuje w nim zbyt duży błąd numeryczny, oraz- specyfikacja modelu jest nieprawidłowa. Dla specyfikacji modelu zamierzonym modelem jest:
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.
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
r
maximum-likelihood
optimization
lme4-nlme
spierać się
źródło
źródło
MLE.glmer
iMLE.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.nAGQ
inglmer
sprawiło, że MLE są porównywalne. Domyślna precyzjaglmer
nie była zbyt dobra.Odpowiedzi:
Ustawienie wysokiej wartości
nAGQ
wglmer
wywołaniu sprawiło, że MLE z dwóch metod są równoważne. Domyślna precyzjaglmer
nie była zbyt dobra. To rozwiązuje problem.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.
źródło