泊松回归回归线剖面95%置信区间计算结果过大问题排查

核心错误：响应均值剖面置信区间的计算逻辑错误 #

你当前代码中得到的橙色区间根本不是预测值的95%剖面置信区间，错误出在这一段逻辑：

(pCI <- confint(glm.pois))
nice.xs <- sort(data$x)
pred.pCI   <- apply(t(pCI), 1, FUN=function(x){exp( x[1]+x[2]*nice.xs)})

confint(glm.pois)返回的是两个回归系数各自的边际剖面置信区间，输出矩阵的行对应截距项、x的系数，列对应2.5%、97.5%分位数。你直接转置后按行提取系数组合，相当于强制同时取两个系数的2.5%分位值计算下限、同时取97.5%分位值计算上限，完全忽略了两个回归系数估计值的协相关性，相当于把两个独立极端事件叠加，得到的区间置信水平远高于95%，所以宽度会远大于标准误法的结果。

次要问题：标准误法的分位数近似误差 #

你计算标准误法置信区间时用了2 * se_link近似95%正态分位数，严格来说应该使用qnorm(0.975)（约1.96），不过这个差异很小，不是导致两个区间差距过大的核心原因。

修正方案 #

如果需要计算响应均值的正确置信区间，可通过系数协方差矩阵计算线性组合的标准误，再反变换到响应尺度，参考修正代码如下：

# 计算基于系数协方差矩阵的正确联合置信区间
vcov_mat <- vcov(glm.pois)
pred_pCI <- t(sapply(nice.xs, function(x_val){
  lin_com <- c(1, x_val)
  fit_lin <- sum(lin_com * coef(glm.pois))
  se_lin <- sqrt(t(lin_com) %*% vcov_mat %*% lin_com)
  ilink(fit_lin + c(-1.96, 1.96)*se_lin)
}))

绘图时把原来的pred.pCI[,1]和pred.pCI[,2]替换为pred_pCI[,1]和pred_pCI[,2]，两个区间就会基本一致。

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