如何從二元邏輯回歸模型中獲得兩個概率之間差異的置信區間？
我在這裡看到了這個問題，但沒有明確的答案。
假設我有一個簡單的二元邏輯回歸模型 $ x: $
$$ \log(p/(1-p)) = \beta_0 + \beta_1x $$
然後我知道：
$$ p/(1-p) = e^{\beta_0 + \beta_1x} $$
和
$$ p = e^{\beta_0 + \beta_1x}/(1 + e^{\beta_0 + \beta_1x}) $$
因此，如果 $ x=0, $ 那麼模型變為：
$$ \log(p/(1-p)) = \beta_0 $$
而如果 $ x=1, $ 那麼模型變為：
$$ \log(p/(1-p)) = \beta_0 + \beta_1 $$
獲得置信區間 $ p $ 什麼時候 $ x=0, $ 我這樣做了：
model <- glm(y~., family=binomial(), data) #For x=0 Bigma = vcov(model) sig = sqrt(Bigma[1,1]) logit_p = coef(model)[1][[1]] #The intercept theta_L = logit_p - 1.96*sig theta_U = logit_p + 1.96*sig p_L = plogis(theta_L) p_U = plogis(theta_U)置信區間為 (p_L, p_U)。
那麼當 x=1 時 p 的置信區間：
sig_x1 = sqrt(Bigma[1,1] + Bigma[2,2] + 2*Bigma[1,2]) logit_p_x1 = coef(model)[1][[1]] + coef(model)[2][[1]] #beta_0 + beta_1 theta_L_x1 = logit_p_x1 - 1.96*sig_x1 theta_U_x1 = logit_p_x1 + 1.96*sig_x1 p_L_x1 = plogis(theta_L_x1) p_U_x1 = plogis(theta_U_x1)置信區間為 (p_L_x1, p_U_x1)。
現在我想要一個成功概率差異的置信區間 $ x=0 $ 和 $ x=1. $
我可以獲得差異的點估計：
$$ p_{x_1} - p_{x_0} = [e^{\beta_0 + \beta_1} - e^{\beta_0}]/[(1+e^{\beta_0 + \beta_1})(1+e^{\beta_0})] $$
我知道下一步是計算差異的標準誤差，但我不知道該怎麼做。
問題：當兩個概率的差異時，置信區間的公式和R代碼是什麼？ $ x=0 $ 什麼時候 $ x=1? $
delta 方法狀態
$$ \operatorname{Var}(g(X)) = [g'(X)]^2 \operatorname{Var}(X) $$
因為這個問題涉及到兩個參數，我們可以把它擴展到多元delta方法
$$ =\nabla g^T , \Sigma , \nabla g $$
這裡，
$$ g = \left[e^{\beta_{0}+\beta_{1}}-e^{\beta_{0}}\right] /\left[\left(1+e^{\beta_{0}+\beta_{1}}\right)\left(1+e^{\beta_{0}}\right)\right] $$
和 $ \Sigma $ 是模型的方差協方差矩陣。 $ \nabla g $ 是……噁心。我不打算用手來做，計算機代數雖然很快會產生一堆符號。但是，您可以使用自微分計算梯度。一旦你計算出方差，那麼它只是你對概率差的估計加/減 1.96 乘以標準偏差（方差的根）。注意，這種方法將產生低於 0 或高於 1 的答案。
我們可以通過以下方式在 R 中執行此操作（注意您需要安裝
autodiffr包）。library(autodiffr) g = function(b) (exp(b[1] + b[2]) - exp(b[1])) / ((1+ exp(b[1] + b[2]))*(1+exp(b[1]))) x = rbinom(100, 1, 0.5) eta = -0.8 + 0.2*x p = plogis(eta) y = rbinom(100, 1, p) model = glm(y~x, family=binomial()) Bigma = vcov(model) grad_g = makeGradFunc(g) nabla_g = grad_g(coef(model)) se = as.numeric(sqrt(nabla_g %*% Bigma %*% nabla_g)) estimate = diff(predict(model, newdata=list(x=c(0, 1)), type='response')) estimate + c(-1, 1)*1.96*se對這個適度的例子重複這個過程表明得到的置信區間接近名義覆蓋率，這是一件好事，但我想隨著概率接近 0 或 1，情況會變得更糟。