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 glm
fungsi:
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
di mana dan
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
dan integral tidak memiliki ekspresi bentuk tertutup, kita perlu menggunakan integrasi numerik seperti Gaussian Quadrature. Kita dapat menggunakan paket statmod
untuk 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 sekarang
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 .
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.
Jawaban:
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
R
kode 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. ThelogMLEgh()
dapat ditemukan di\mixed_models_data.zip\MixedModels\Chapter07
. Dia tidak menggunakanstatmod
paket untuk mendapatkan kuadratur, tetapi menulis fungsinya sendirigauher()
. 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 padaglm()
.sumber
b
's dengan 10 node, bagaimana kita akan dapat melipatgandakan matriksZ
dang
? Atau saya salah sepenuhnya?Z = rep(1,n)
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
Z
digunakan 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 perluZ
lagi, cukup evaluasi integral untuk setiap cluster, dan kemudian jumlahkan bersama-sama. Matriks desain untuk intersep acak hanyalah vektor 1s.