R语言负二项回归模型的预测区间求解及新数据报错问题

问题根源 #

ciTools::add_pi()处理glm.nb模型时，内部模拟逻辑依赖原始训练数据集的样本量（这里是50）。当你在原数据后追加3行新数据，总样本量变为53，模拟生成的随机数长度与新数据长度不匹配，从而触发报错。

方案1：单独构建新预测数据集（最简解决） #

不要在原cars数据中追加新行，而是单独创建仅包含待预测点的数据框，再传入add_pi()：

# 训练负二项回归模型
fit <- MASS::glm.nb(dist ~ speed, data = cars)

# 创建新的预测数据
new_data <- tibble::tibble(speed = c(58, 63, 121))

# 生成预测区间
ciTools::add_pi(new_data, fit)

如果需要合并原数据与新数据的预测结果，可以分别处理后再拼接：

# 原数据的预测结果
cars_preds <- ciTools::add_pi(cars, fit)

# 新数据的预测结果
new_preds <- ciTools::add_pi(new_data, fit)

# 合并两个结果
combined_preds <- dplyr::bind_rows(cars_preds, new_preds)

# 查看前5行原数据结果和后3行新数据结果
head(combined_preds, 5)
tail(combined_preds, 3)

方案2：手动模拟预测区间（不依赖ciTools） #

基于负二项分布的统计性质，手动模拟生成预测区间，避免第三方包的内部限制：

# 训练模型
fit <- MASS::glm.nb(dist ~ speed, data = cars)

# 定义新预测数据
new_data <- tibble::tibble(speed = c(58, 63, 121))

# 预测每个新点的均值响应
mu <- predict(fit, newdata = new_data, type = "response")

# 提取模型的离散参数theta
theta <- fit$theta

# 模拟10000次负二项分布的响应值（固定种子保证可复现）
set.seed(123)
simulated_vals <- replicate(10000, rnbinom(nrow(new_data), size = theta, mu = mu))

# 计算95%预测区间的分位数
pred_intervals <- t(apply(simulated_vals, 1, quantile, probs = c(0.025, 0.975)))

# 整理最终结果
final_result <- cbind(
  new_data,
  pred = round(mu, 4),
  LPB0.025 = pred_intervals[, 1],
  UPB0.975 = pred_intervals[, 2]
)

final_result

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