Manipulasi model regresi logistik

Saya ingin memahami apa yang dilakukan kode berikut. Orang yang menulis kode tidak lagi berfungsi di sini dan hampir tidak berdokumen. Saya diminta untuk menyelidikinya oleh seseorang yang berpikir " ini model regresi logistik bayesian "

bglm <- function(Y,X) {
    # Y is a vector of binary responses
    # X is a design matrix

    fit <- glm.fit(X,Y, family = binomial(link = logit))
    beta <- coef(fit)
    fs <- summary.glm(fit)
    M <- t(chol(fs$cov.unscaled))
    betastar <- beta + M %*% rnorm(ncol(M))
    p <- 1/(1 + exp(-(X %*% betastar)))
    return(runif(length(p)) <= p)
}

Saya bisa melihat bahwa itu cocok dengan model logistik, mengambil transpos factorisation Cholseky dari matriks kovarian yang diperkirakan, mengalikannya dengan vektor gambar dari dan kemudian ditambahkan ke perkiraan model. Ini kemudian ditakdirkan sebelumnya oleh matriks desain, logit terbalik ini diambil, dibandingkan dengan vektor gambar dari dan vektor biner yang dihasilkan dikembalikan. Tapi apa semua ini berarti secara statistik?$N(0,1)$$U(0,1)$

Apa fungsi ini:
Pada dasarnya, fungsi menghasilkan data respon pseudorandom baru (yaitu, ) dari model data Anda. Model yang digunakan adalah model frequentist standar. Seperti biasa, ini mengasumsikan bahwa data * Anda adalah konstanta yang dikenal - mereka tidak dijadikan sampel dengan cara apa pun. Apa yang saya lihat sebagai fitur penting dari fungsi ini adalah menggabungkan ketidakpastian tentang parameter yang diestimasi. $Y$$X$

* Perhatikan bahwa Anda harus secara manual menambahkan vektor sebagai kolom paling kiri dari matriks Anda sebelum memasukkannya ke fungsi, kecuali jika Anda ingin menekan intersep (yang umumnya bukan ide yang baik).$1$$X$

Apa gunanya fungsi ini:
Jujur saya tidak tahu. Itu bisa menjadi bagian dari rutinitas MCMC Bayesian, tetapi itu hanya akan menjadi satu bagian - Anda akan membutuhkan lebih banyak kode di tempat lain untuk benar-benar menjalankan analisis Bayesian. Saya tidak merasa cukup ahli dalam metode Bayesian untuk mengomentari hal ini secara pasti, tetapi fungsinya tidak 'terasa' bagi saya seperti apa yang biasanya digunakan.

Itu juga bisa digunakan dalam analisis daya berbasis simulasi. (Lihat jawaban saya di sini: Simulasi analisis daya regresi logistik - eksperimen yang dirancang , untuk informasi tentang hal semacam ini.) Perlu dicatat bahwa analisis daya berdasarkan data sebelumnya yang tidak memperhitungkan ketidakpastian estimasi parameter ke dalam perhitungan sering optimis. (Saya membahas poin itu di sini: Ukuran efek yang diinginkan vs ukuran efek yang diharapkan .)

Jika Anda ingin menggunakan fungsi ini:
Seperti yang dicatat @whuber dalam komentar, fungsi ini tidak efisien. Jika Anda ingin menggunakan ini untuk melakukan (misalnya) analisis daya, saya akan membagi fungsi menjadi dua fungsi baru. Yang pertama akan membaca data Anda dan menampilkan parameter dan ketidakpastian. Fungsi baru kedua akan menghasilkan data pseudorandom baru . Berikut ini adalah contoh (walaupun dimungkinkan untuk memperbaikinya lebih lanjut): $Y$

simulationParameters <- function(Y,X) {
                        # Y is a vector of binary responses
                        # X is a design matrix, you don't have to add a vector of 1's 
                        #   for the intercept

                        X    <- cbind(1, X)  # this adds the intercept for you
                        fit  <- glm.fit(X,Y, family = binomial(link = logit))
                        beta <- coef(fit)
                        fs   <- summary.glm(fit)
                        M    <- t(chol(fs$cov.unscaled))

                        return(list(betas=beta, uncertainties=M))
}

simulateY <- function(X, betas, uncertainties, ncolM, N){

             # X      <- cbind(1, X)  # it will be slightly faster if you input w/ 1's
             # ncolM  <- ncol(uncertainties) # faster if you input this
             betastar <- betas + uncertainties %*% rnorm(ncolM)
             p        <- 1/(1 + exp(-(X %*% betastar)))

             return(rbinom(N, size=1, prob=p))
}