Kemungkinan dan perkiraan untuk efek-efek campuran Regresi logistik

8

Pertama mari kita mensimulasikan beberapa data untuk regresi logistik dengan bagian-bagian tetap dan acak:

set.seed(1)
n    <- 100
x    <- runif(n)
z    <- sample(c(0,1), n, replace=TRUE)
b    <- rnorm(2)
beta <- c(0.4, 0.8)
X    <- model.matrix(~x)
Z    <- cbind(z, 1-z)
eta  <- X%*%beta + Z%*%b
pr   <- 1/(1+exp(-eta))
y    <- rbinom(n, 1, pr)

Jika kami hanya ingin menyesuaikan regresi logistik tanpa bagian acak, kami dapat menggunakan glmfungsi:

glm(y~x, family="binomial")

glm(y~x, family="binomial")$coefficients
# (Intercept)           x 
#  -0.2992785   2.1429825 

Atau membangun fungsi kita sendiri dari kemungkinan log

logL(β)=i=1nyilogΛ(ηi)+(1yi)log(1Λ(ηi))

di mana Λ(ηi)=11+exp(ηi) dan ηi=Xiβ dan gunakan optim()untuk memperkirakan parameter yang memaksimalkannya, seperti berikut ini contoh kode:

ll.no.random <- function(theta,X,y){
  beta <- theta[1:ncol(X)]
  eta  <- X%*%beta
  p    <- 1/(1+exp(-eta))
  ll   <- sum( y*log(p) + (1-y)*log(1-p) )
  -ll
}

optim(c(0,1), ll.no.random, X=X, y=y)

optim(c(0,1), ll.no.random, X=X, y=y)$par
# -0.2992456  2.1427484

yang tentu saja memberikan estimasi yang sama dan memaksimalkan log-likelihood untuk nilai yang sama. Untuk efek campuran, kami ingin sesuatu seperti

library(lme4)
glmer(y~x + (1|z), family="binomial")

Tetapi bagaimana kita dapat melakukan hal yang sama dengan fungsi kita sendiri? Karena kemungkinannya adalah

L=j=1JPr(y1j,...,ynjj|x,bj)f(bj)dbj

dan integral tidak memiliki ekspresi bentuk tertutup, kita perlu menggunakan integrasi numerik seperti Gaussian Quadrature. Kita dapat menggunakan paket statmoduntuk mendapatkan beberapa quadratures, katakan 10

library(statmod)    
gq <- gauss.quad(10)
w  <- gq$weights
g  <- gq$nodes

UPDATE: Dengan menggunakan lokasi quadrature ini dan bobot untuk ( sini), kita dapat memperkirakan integral dengan dengan menjumlahkan persyaratan dengan disubstitusi untuk dan keseluruhan istilah dikalikan dengan bobot masing-masing . Jadi, fungsi kemungkinan kita seharusnya sekaranggrwrr=1,...,RR=10bjRgrbjwr

L=j=1Jr=1RPr(y1j,...,ynjj|x,bj=gr)wr

Juga, kita perlu memperhitungkan varians dari bagian acak, saya membaca bahwa ini dapat dicapai dengan mengganti dalam fungsi kami dengan mana , jadi dalam fungsi kemungkinan di atas kita benar-benar mengganti dengan 'dan bukan .bjN(0,σb2)ησjθjθjN(0,1)θgβ

Satu masalah komputasi yang tidak saya dapatkan adalah bagaimana cara mengganti istilah karena vektor tidak akan memiliki panjang yang sama. Tapi mungkin saya tidak mengerti itu, karena saya kehilangan sesuatu yang penting di sini, atau salah paham bagaimana metode ini bekerja.

Steve
sumber
Saya tidak punya waktu untuk melihat dengan baik sekarang, tetapi ini sepertinya digunakan untuk MCMC.
shadowtalker
@ssdecontrol terima kasih, MCMC adalah alternatif yang bagus. Tetapi saya ingin menerapkan pendekatan klasik.
Steve
Apa yang non-klasik tentang MCMC untuk mengevaluasi integral?
shadowtalker
@ssdecontrol dengan baik, saya tidak yakin ... Tapi semua buku yang saya periksa dan lme4, paket R biasa, jangan gunakan MCMC. Jadi, saya ingin tetap dengan itu. Setidaknya di awal.
Steve
1
Sudahkah Anda menanyakan hal ini dalam daftar R-sig-ME ([email protected])? Beberapa orang mungkin dapat membantu Anda lebih lanjut. Selain itu: Saya akan sangat mendesak Anda untuk mempelajari makalah yang Efisien Laplacian dan Adaptive Gaussian Quadrature Algorithms untuk Multilevel Generalized Linear Mixed Models oleh Pinheiro dan Chao. Ini berisi hasil mengenai kinerja AGQ serta aljabar linier di belakangnya. Jika Anda ingin kode mereka .... well, bersiap-siap untuk beberapa dekomposisi matriks serius jarang. : D
usεr11852

Jawaban:

2

Saya tidak melihat bagaimana "vektor-vektor tidak akan memiliki panjang yang sama", mohon klarifikasi pertanyaan Anda.

Pertama-tama, untuk integral dengan dimensi kurang dari 4, metode numerik langsung seperti quadrature lebih efisien daripada MCMC. Saya mempelajari pertanyaan-pertanyaan ini sebentar, dan saya akan senang mendiskusikan masalah ini dengan Anda.

Untuk regresi logistik efek campuran, satu-satunya Rkode eksplisit yang saya temukan adalah dari buku Prof Demidenko, Model Campuran: Teori dan Aplikasi , Anda dapat mengunduh kode tersebut melalui kolom "PERANGKAT LUNAK DAN DATA" di halaman web. The logMLEgh()dapat ditemukan di \mixed_models_data.zip\MixedModels\Chapter07. Dia tidak menggunakan statmodpaket untuk mendapatkan kuadratur, tetapi menulis fungsinya sendiri gauher(). Ada beberapa kesalahan kecil dalam kode dan saya telah membahasnya dengan penulis, tetapi masih sangat membantu untuk memulai dari kode dan bukunya. Saya dapat memberikan versi yang diperbaiki jika diperlukan.

Masalah lain adalah bahwa, jika Anda ingin mendapatkan perkiraan yang akurat, optim()tidak cukup, Anda mungkin perlu menggunakan metode seperti skoring Fisher, seperti pada glm().

Randel
sumber
buku ini sepertinya memiliki banyak info tentang apa yang saya kerjakan. Kode itu sendiri tidak banyak bicara - hanya memesan buku jadi saya harus menunggu untuk itu. Hal tentang vektor adalah bahwa jika kita mengganti syarat, 2 b's dengan 10 node, bagaimana kita akan dapat melipatgandakan matriks Zdan g? Atau saya salah sepenuhnya?
Steve
Saya tahu bahwa saya harus melangkah lebih jauh untuk mendapatkan perkiraan yang akurat, tetapi saya berharap untuk memahami GQ sebagai langkah pertama.
Steve
Anda dapat mempratinjau dua edisi di Google Buku terlebih dahulu. Dalam kode Anda, adalah matriks ? Tapi hanya skalar? Anda dapat mulai dari model intersep acak terlebih dahulu. Jika dimensi efek acak adalah dua, Anda benar-benar membutuhkan 100 node dua dimensi, 10 node pada setiap dimensi. Zn×2bZ = rep(1,n)
Randel
maaf tapi semakin saya memikirkannya semakin membingungkan saya. Jika saya melakukan itu Z=rep(1,n)maka saya akan mendapatkan satu intersepsi acak untuk setiap baris, artinya setiap individu adalah sebuah grup? Dalam contoh saya, saya memiliki dua grup, sehingga kami memiliki dan untuk memberikan kami butuhkan. Tidak? Z%*%b n×22×1n×1
Steve
Oh, saya baru saja memperhatikan bahwa Zdigunakan untuk memisahkan intersep acak untuk setiap cluster, bukan matriks desain untuk efek acak. Maka Anda benar, tetapi Anda harus mengevaluasi integral dan memanfaatkan quadrature secara terpisah untuk setiap cluster. Anda tidak perlu Zlagi, cukup evaluasi integral untuk setiap cluster, dan kemudian jumlahkan bersama-sama. Matriks desain untuk intersep acak hanyalah vektor 1s.
Randel