Численно получить MLE из GLMM сложно, и на практике, я знаю, мы не должны использовать оптимизацию методом грубой силы (например, используя optim
простой способ). Но для моих собственных образовательных целей я хочу попробовать, чтобы убедиться, что я правильно понимаю модель (см. Код ниже). Я обнаружил, что всегда получаю противоречивые результаты glmer()
.
В частности, даже если я использую MLE в glmer
качестве начальных значений, в соответствии с функцией правдоподобия, которую я написал ( negloglik
), они не являются MLE ( opt1$value
меньше, чем opt2
). Я думаю, что две возможные причины:
negloglik
плохо написано, так что в нем слишком много числовых ошибок, и- спецификация модели неверна. Для спецификации модели предполагаемая модель:
где - биномиальный коэффициент pmf, а - нормальный pdf. Я пытаюсь оценить , и . В частности, я хочу знать, является ли спецификация модели неправильной, какова правильная спецификация.
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...
Более простой пример
Чтобы уменьшить вероятность получения большой числовой ошибки, я создал более простой пример.
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
каламбур
источник
источник
MLE.glmer
иMLE.optim
), особенно для случайного эффекта (см. Новый пример), поэтому я думаю, что он основан не только на некотором постоянном факторе значений вероятности.nAGQ
inglmer
сделало MLE сопоставимыми. Точность по умолчаниюglmer
была не очень хорошей.Ответы:
Установка высокого значения
nAGQ
вglmer
вызове сделала MLEs из двух методов эквивалентными. Точность по умолчаниюglmer
была не очень хорошей. Это решает проблему.См. Ответ @ SteveWalker здесь. Почему я не могу сопоставить вывод glmer (family = binomial) с ручной реализацией алгоритма Гаусса-Ньютона? Больше подробностей.
источник