R中Delta法计算灵敏度置信区间的正确性及方法选择咨询

我基于二项式GLM模型计算灵敏度，因变量为测量值，自变量为真实值。已知需用Delta法或Bootstrapping计算置信区间（而非直接用confint函数），以下是我的数据示例、模型及计算代码：

数据与模型代码 #

# create some fake data 
set.seed(1234)
neighb <- sample(1:2, 1846, replace=TRUE, prob=c(0.6, 0.4))
neighb <- factor(neighb)
ind <- sample(1:2, 1846, replace=TRUE, prob=c(0.65, 0.35))
ind <- factor(ind)
df <- cbind.data.frame(ind, neighb)

mq1 <- glm(neighb ~ ind, data=df, family='binomial')
coef_mq1 <- coef(mq1)

灵敏度计算代码 #

m_seq1 <- (1) / (1+ (exp(-(coef_mq1[1]+1*(coef_mq1[2]))))) 
m_seq1

Delta法置信区间计算代码 #

library(msm)
b0 <- coef_mq1[1]
b1 <- coef_mq1[2]
se <- deltamethod(~((1) / (1+ (exp(-(x1+1*(x2)))))) , c(b1,b0), vcov(mq1))
m_seq1 + 1.96*se
m_seq1 - 1.96*se

我的问题：

1. 直接将公式代入deltamethod函数是否正确，是否存在遗漏？
2. 是否有理由优先使用Bootstrapping而非Delta法？

回答 #

关于Delta法代码的正确性 #

你的Delta法代码存在参数顺序匹配错误，需要修正：

• deltamethod的第二个参数是系数向量，顺序必须和公式里的x1、x2对应；同时方差矩阵vcov(mq1)的列顺序与coef(mq1)完全一致（先截距b0，再斜率b1）。
• 你当前公式用x1代表b1、x2代表b0，但传入的系数向量是c(b1,b0)，而vcov(mq1)的顺序是b0在前、b1在后，这会导致方差矩阵与系数不匹配，计算出的标准误差是错误的。

修正后的代码：

# 保持公式中x1对应b0，x2对应b1，与系数、方差矩阵顺序一致
se <- deltamethod(~ 1 / (1 + exp(-(x1 + x2))), coef_mq1, vcov(mq1))
# 计算95%置信区间
ci_upper <- m_seq1 + 1.96 * se
ci_lower <- m_seq1 - 1.96 * se

额外补充：

• 你计算灵敏度的公式是正确的：对于二项式GLM，当自变量ind取第二个水平（对应模型中的x=1，R中因子默认以第一个水平为参照），预测概率就是灵敏度（假设neighb=2代表阳性测量结果，ind=2代表真实阳性）。
• 公式里的1*(coef_mq1[2])可以简化为coef_mq1[2]，不影响结果但更简洁。

优先选择Bootstrapping的场景 #

当满足以下任意一种情况时，优先使用Bootstrapping而非Delta法：

• 模型渐近假设不成立：Delta法依赖GLM的系数渐近正态性假设，如果样本量小、模型存在分离现象或拟合效果差，该假设不成立，Delta法的置信区间会出现偏差。
• 非线性变换复杂度高：灵敏度是logit逆变换（sigmoid函数），Delta法对这种变换表现尚可，但如果衍生统计量是更复杂的非线性函数，Bootstrapping的鲁棒性更强。
• 避免导数推导成本：Delta法需要对目标函数求导，虽然msm包会自动计算，但如果统计量公式复杂，手动验证导数正确性的成本高，Bootstrapping无需推导导数，操作更省心。
• 小样本场景：样本量较小时，Bootstrapping能更准确反映真实抽样分布，而Delta法的渐近假设会导致区间精度不足。

当然，如果样本量足够大、模型拟合良好，Delta法计算速度更快，结果也会和Bootstrapping非常接近。

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