R语言中Logistic回归：特定结果y对应预测变量x的置信区间

核心结论 #

不能直接用模型参数的置信区间来推导对应概率的年龄置信区间，因为Logistic模型的逆变换（从概率求年龄）是非线性的，参数的线性置信区间经过非线性变换后不再是准确的置信区间，必须用delta方法或者bootstrap抽样来计算。

方法一：Delta方法（推荐，基于渐近正态性） #

Delta方法通过对逆变换函数求导，结合参数的协方差矩阵计算年龄值的标准误，进而构造置信区间。

步骤与代码实现 #

1. 明确逆变换函数：对于失败概率$p$，对应的年龄$x = \frac{\text{logit}(p) - b_0}{b_1}$，其中$\text{logit}(p) = \ln\left(\frac{p}{1-p}\right)$
2. 对$x$关于模型参数$b_0$、$b_1$求偏导：
   • $\frac{\partial x}{\partial b_0} = -\frac{1}{b_1}$
   • $\frac{\partial x}{\partial b_1} = -\frac{x}{b_1}$
3. 结合参数协方差矩阵计算$x$的方差与标准误，最终得到95%置信区间（$x \pm 1.96 \times SE$）

# 假设logit2是已拟合完成的Logistic回归模型
library(MASS)

# 定义逆变换函数
fun2 <- function(prob, b0, b1){
  (log(prob/(1-prob)) - b0)/b1
}

# 以p(失败)=0.5（对应p(成功)=0.5）为例计算
target_prob <- 0.5
b0 <- coef(logit2)[1]
b1 <- coef(logit2)[2]

# 年龄点估计
x_hat <- fun2(target_prob, b0, b1)

# 获取参数协方差矩阵
vcov_mat <- vcov(logit2)
var_b0 <- vcov_mat[1,1]
var_b1 <- vcov_mat[2,2]
cov_b0b1 <- vcov_mat[1,2]

# 计算偏导数
db0 <- -1/b1
db1 <- -x_hat/b1

# 计算年龄的方差与标准误
var_x <- db0^2 * var_b0 + db1^2 * var_b1 + 2 * db0 * db1 * cov_b0b1
se_x <- sqrt(var_x)

# 生成95%置信区间
ci_x <- c(x_hat - 1.96*se_x, x_hat + 1.96*se_x)
cat("p(失败)=0.5对应的年龄95%置信区间：", ci_x, "\n")

# 批量计算多个概率的置信区间（排除p=1.0，因logit(p)无意义）
prob_vec <- c(0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9)
compute_ci <- function(prob){
  x_hat <- fun2(prob, b0, b1)
  db0 <- -1/b1
  db1 <- -x_hat/b1
  var_x <- db0^2 * var_b0 + db1^2 * var_b1 + 2 * db0 * db1 * cov_b0b1
  se_x <- sqrt(var_x)
  c(lower=x_hat-1.96*se_x, upper=x_hat+1.96*se_x)
}

ci_results <- t(sapply(prob_vec, compute_ci))
cbind(prob=prob_vec, age_hat=sapply(prob_vec, fun2, b0, b1), ci_results)

方法二：Bootstrap抽样（更稳健，适合小样本） #

通过重复抽样原数据并拟合模型，得到多组参数估计值，再计算每组参数对应的年龄值，最后取分位数作为置信区间。

代码实现 #

set.seed(123) # 固定随机种子保证结果可重复
n_boot <- 1000 # 抽样次数
boot_x <- numeric(n_boot)

for(i in 1:n_boot){
  # 有放回抽样原数据
  boot_data <- logit2$data[sample(nrow(logit2$data), replace=TRUE),]
  # 重新拟合Logistic模型
  boot_model <- glm(y ~ age, data=boot_data, family=binomial)
  # 计算目标概率对应的年龄
  boot_x[i] <- fun2(target_prob, coef(boot_model)[1], coef(boot_model)[2])
}

# 取2.5%和97.5%分位数作为95%置信区间
boot_ci <- quantile(boot_x, c(0.025, 0.975))
cat("Bootstrap法得到的p(失败)=0.5对应的年龄95%置信区间：", boot_ci, "\n")

常见误区解析 #

直接将$b_0$、$b_1$的置信区间代入逆变换函数计算年龄区间是错误的——逆变换是非线性操作，参数的线性置信区间经过变换后会偏离真实的年龄置信区间，无法保证准确性。

内容的提问来源于stack exchange，提问作者laura