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:
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
MLE.glmer
i 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.
nAGQ
in glmer
sprawiło, że MLE są porównywalne. Domyślna precyzja glmer
nie była zbyt dobra.