混合效應的可能性和估計 Logistic 回歸

首先讓我們模擬一些具有固定和隨機部分的邏輯回歸數據：

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)

如果我們只想擬合沒有隨機部分的 Logistic 回歸，我們可以使用以下glm函數：

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

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

或者構建我們自己的對數似然函數

在哪裡和 並用於optim()估計使其最大化的參數，如以下示例代碼所示：

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

這當然給出了相同的估計並最大化相同值的對數似然。對於混合效果，我們想要類似的東西

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

但是我們怎樣才能對我們自己的函數做同樣的事情呢？因為可能性是

並且積分沒有封閉形式的表達式，我們需要使用像高斯求積這樣的數值積分。我們可以使用這個包statmod得到一些正交，比如 10

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

**更新：**使用這些正交位置和權重為了(這裡），我們可以近似積分由總和條款與替代和整個術語乘以各自的權重. 因此，我們的似​​然函數現在應該是

另外，我們需要考慮隨機部分的方差，我讀到這可以通過替換在我們的功能與在哪裡, 所以在上面的似然函數中我們實際上替換與的而不是的。

我沒有得到的一個計算問題是如何替換這些術語，因為向量的長度不同。但可能我不明白這一點，因為我在這裡遺漏了一些關鍵的東西，或者誤解了這種方法的工作原理。

我沒有看到“向量的長度不同”，請澄清您的問題。

首先，對於維數小於 4 的積分，像求積這樣的直接數值方法比 MCMC 更有效。我研究了這些問題一段時間，很高興與您討論這個問題。

對於混合效應邏輯回歸，R我找到的唯一顯式代碼來自 Demidenko 教授的《混合模型：理論與應用》一書，您可以通過網頁上的“軟件和數據”欄目下載代碼。logMLEgh()可以在 中找到\mixed_models_data.zip\MixedModels\Chapter07。他沒有使用statmod包來獲取求積，而是編寫了自己的函數gauher()。代碼中有一些小錯誤，我已經和作者討論過了，但是從他的代碼和書開始還是很有幫助的。如果需要，我可以提供更正的版本。

另一個問題是，如果你想得到準確的估計optim()是不夠的，你可能需要使用像 Fisher 評分這樣的方法，如glm().

引用自：https://stats.stackexchange.com/questions/112191