Secara numerik menurunkan MLE pada GLMM adalah sulit dan, dalam praktiknya, saya tahu, kita tidak boleh menggunakan optimasi brute force (misalnya, menggunakan optim
dengan cara sederhana). Tetapi untuk tujuan pendidikan saya sendiri, saya ingin mencobanya untuk memastikan saya memahami model dengan benar (lihat kode di bawah). Saya menemukan bahwa saya selalu mendapatkan hasil yang tidak konsisten dari glmer()
.
Secara khusus, bahkan jika saya menggunakan MLE dari glmer
sebagai nilai awal, sesuai dengan fungsi likelihood yang saya tulis ( negloglik
), mereka bukan MLE ( opt1$value
lebih kecil dari opt2
). Saya pikir dua alasan potensial adalah:
negloglik
tidak ditulis dengan baik sehingga ada terlalu banyak kesalahan numerik di dalamnya, dan- spesifikasi model salah. Untuk spesifikasi model, model yang dimaksud adalah:
mana f adalah PMF binomial danadalah pdf normal. Saya mencoba memperkirakan,, dan. Secara khusus, saya ingin tahu jika spesifikasi model salah, apa spesifikasi yang benar.
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...
Contoh yang lebih sederhana
Untuk mengurangi kemungkinan memiliki kesalahan numerik yang besar, saya membuat contoh yang lebih sederhana.
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
sumber
MLE.glmer
danMLE.optim
) terutama untuk efek acak (lihat contoh baru), jadi itu bukan hanya berdasarkan pada beberapa faktor konstan dalam nilai-nilai kemungkinan, saya pikir.nAGQ
dalamglmer
membuat MLEs sebanding. Ketepatan standarglmer
tidak terlalu baik.Jawaban:
Menetapkan nilai tinggi
nAGQ
dalamglmer
panggilan membuat MLE dari kedua metode tersebut setara. Ketepatan standarglmer
tidak terlalu baik. Ini menyelesaikan masalah.Lihat jawaban @ SteveWalker di sini. Mengapa saya tidak dapat mencocokkan output glmer (keluarga = binomial) dengan penerapan algoritma Gauss-Newton secara manual? untuk lebih jelasnya.
sumber