Mengapa glmer tidak mencapai kemungkinan maksimum (seperti yang diverifikasi dengan menerapkan optimasi generik lebih lanjut)?

37

Secara numerik menurunkan MLE pada GLMM adalah sulit dan, dalam praktiknya, saya tahu, kita tidak boleh menggunakan optimasi brute force (misalnya, menggunakan optimdengan 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 glmersebagai nilai awal, sesuai dengan fungsi likelihood yang saya tulis ( negloglik), mereka bukan MLE ( opt1$valuelebih kecil dari opt2). Saya pikir dua alasan potensial adalah:

  1. negloglik tidak ditulis dengan baik sehingga ada terlalu banyak kesalahan numerik di dalamnya, dan
  2. 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.

L.=saya=1n(-f(ysaya|N,Sebuah,b,rsaya)g(rsaya|s)drsaya)
fa b sgSebuahbs
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 
berdalih
sumber
apakah MLE (bukan log-kemungkinan itu sendiri) sebanding? Yaitu, apakah Anda hanya pergi dengan konstan?
Ben Bolker
1
Perkiraan MLE jelas berbeda ( MLE.glmerdan MLE.optim) terutama untuk efek acak (lihat contoh baru), jadi itu bukan hanya berdasarkan pada beberapa faktor konstan dalam nilai-nilai kemungkinan, saya pikir.
berdalih
4
@Ben Menetapkan nilai tinggi nAGQdalam glmermembuat MLEs sebanding. Ketepatan standar glmertidak terlalu baik.
berdalih
5
Menautkan ke pertanyaan lme4 serupa yang @Steve Walker membantu saya dengan: stats.stackexchange.com/questions/77313/…
Ben Ogorek
3
Sebagai pertanyaan yang lebih tua dengan banyak peningkatan, ini mungkin bisa dilakukan kakek. Saya tidak melihat perlunya ditutup.
gung - Reinstate Monica

Jawaban:

3

Menetapkan nilai tinggi nAGQdalam glmerpanggilan membuat MLE dari kedua metode tersebut setara. Ketepatan standar glmertidak terlalu baik. Ini menyelesaikan masalah.

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

Lihat jawaban @ SteveWalker di sini. Mengapa saya tidak dapat mencocokkan output glmer (keluarga = binomial) dengan penerapan algoritma Gauss-Newton secara manual? untuk lebih jelasnya.

berdalih
sumber
1
Tetapi perkiraan kemungkinan log sangat berbeda (mungkin oleh beberapa konstanta), sehingga metode yang berbeda tidak boleh dicampur.
berdalih
1
hmm, menarik / mengejutkan - terima kasih telah menyiapkan contoh ini, saya akan mencoba mencari waktu untuk melihatnya.
Ben Bolker