如何提取MASS包polr有序logistic回归预测值的标准误与变异

polr默认的predict()方法不返回预测概率的标准误，这类由模型参数非线性转换得到的统计量，标准误可以通过Delta方法（速度快，大样本下表现好）或参数Bootstrap方法（稳健，小样本下更可靠）计算，对应示例代码的可复现实现如下：

方法1：Delta方法（推荐日常使用） #

Delta方法通过参数的方差协方差矩阵，结合非线性函数的梯度推导统计量的抽样方差，不需要大量重复抽样，计算速度快。

# 提取模型核心组件：注意polr的参数顺序是先阈值项(zeta)，后自变量系数(coefficients)
coef_all <- c(house.plr$zeta, house.plr$coefficients)
vcov_all <- vcov(house.plr)
n_response <- nlevels(housing$Sat) # 响应变量的类别数，示例中为3

# 自定义预测概率计算函数：输入参数、新数据集，返回和原生predict()一致的展平概率向量
calc_predprob <- function(params, newdata, n_level) {
  tau <- params[1:(n_level-1)] # 前n_level-1个为阈值
  beta <- params[n_level:length(params)] # 剩余为自变量回归系数
  # 构造和polr内部编码一致的模型矩阵（无截距）
  X <- model.matrix(~ Infl + Type + Cont - 1, data = newdata)
  linpred <- as.vector(X %*% beta)
  # 计算累积概率
  cum_prob <- matrix(1, nrow = nrow(newdata), ncol = n_level)
  for (k in 1:(n_level-1)) {
    cum_prob[,k] <- plogis(tau[k] - linpred)
  }
  # 累积概率转各分类的边际概率
  marg_prob <- cbind(cum_prob[,1], cum_prob[,2:n_level] - cum_prob[,1:(n_level-1)])
  colnames(marg_prob) <- levels(housing$Sat)
  return(c(t(marg_prob))) # 按行展平，和后续输出格式对齐
}

# 校验：自定义函数输出是否和原生predict结果完全一致
stopifnot(all.equal(
  as.vector(t(predict(house.plr, test_dat, type = "p"))),
  calc_predprob(coef_all, test_dat, n_response)
))

# 用Delta方法计算标准误，调用msm包的deltamethod函数实现数值梯度计算
# 若未安装msm包，先运行 install.packages("msm")
library(msm)
pred_se <- deltamethod(
  g = calc_predprob,
  mean = coef_all,
  cov = vcov_all,
  newdata = test_dat,
  n_level = n_response
)

# 整理结果输出
pred_result <- data.frame(
  test_id = rep(1:nrow(test_dat), each = n_response),
  Type = rep(test_dat$Type, each = n_response),
  sat_level = rep(levels(housing$Sat), nrow(test_dat)),
  pred_prob = as.vector(t(predict(house.plr, test_dat, type = "p"))),
  pred_se = pred_se
)
print(pred_result, digits = 3)

运行后输出的pred_se列就是对应预测概率的标准误，示例数据下第一组观测（Type为Tower）Low满意度的预测概率0.378对应的标准误约为0.042。

方法2：参数Bootstrap方法 #

如果不信任大样本渐近假设，可以用参数Bootstrap：从模型系数的渐近多元正态分布中重复抽样，每次抽样都计算对应的预测概率，最终取预测概率抽样分布的标准差作为标准误。

set.seed(123) # 固定随机种子保证结果可复现
n_boot <- 1000 # 抽样次数，1000次即可达到足够精度
# 从系数的渐近分布抽样
coef_boot_sample <- MASS::mvrnorm(n = n_boot, mu = coef_all, Sigma = vcov_all)
# 预分配矩阵存储每次抽样的预测结果
boot_pred_mat <- matrix(NA, nrow = n_boot, ncol = nrow(test_dat)*n_response)
for (i in 1:n_boot) {
  boot_pred_mat[i,] <- calc_predprob(coef_boot_sample[i,], test_dat, n_response)
}
# 计算各预测概率的抽样标准差，即标准误
pred_se_boot <- apply(boot_pred_mat, 2, sd)
# 结果和Delta方法的输出差异通常在1e-3量级，一致性很高

踩坑提示 #

• 不要用重采样原始观测的非参数Bootstrap计算该标准误，非参数Bootstrap得到的是包含观测抽样变异的预测区间标准差，不是模型参数估计带来的预测均值标准误，和需要的结果不匹配
• 自定义预测概率函数时，必须保证因子编码、变量顺序、参数顺序和polr内部逻辑完全一致，代码里加的stopifnot校验步骤不要省略，避免编码错误导致结果偏差
• 如果模型中包含偏移项（offset）、权重项，需要同步在calc_predprob函数中加入对应计算逻辑，否则结果会出错

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