R 上泊松回归的预测区间

Posted 2023-03-12

【中文标题】R 上泊松回归的预测区间

【英文标题】：Prediction intervals for poisson regression on R

【发布时间】：2013-07-29 03:50:57

【问题描述】：

这两种方法我都试过了，但我发现两种方法都有困难.. 在我用这两种方法告诉你我的问题之前，我会尝试更好地解释我的问题。

我有数据集“acceptances”，其中我有医院每天的接受数量，其中包含前面描述的独立变量。医院有三个地方可以进行访问。所以在我的数据集中，我每天有 3 行，每个地方都有一行。数据集看起来像：

Date        Place    NumerAccept    weekday month   NoConvention    Rain

2008-01-02  Place1        203       wed     Gen         0             1
2008-01-02  Place2         70       wed     Gen         0             1
2008-01-02  Place3          9       wed     Gen         0             1
2008-01-03  Place1        345       thu     Gen         0             1
2008-01-03  Place2         24       thu     Gen         0             1
2008-01-03  Place3         99       thu     Gen         0             1
2008-01-04  Place1        339       fri     Gen         0             0
2008-01-04  Place2         36       fri     Gen         0             0
2008-01-04  Place3        101       fri     Gen         0             0

....等等...我有数据集直到昨天，所以最后三行是 2013 年 7 月 29 日昨天的接受。 现在我进行泊松回归：

poisson_reg=glm(NumeberAccept ~ 1 + weekday + month + place + NoConvention + Rain, 
                    family = poisson(link = log), data = acceptances)

现在，对于我的预测，我创建了一个新的数据集acceptances_2，我想从中计算未来 2 个月的接受数的预测间隔！！所以第一行是今天的录用数，最后一行是9月29日的录用数。

---

我不知道这个问题是否已经有了答案，但我找不到。我正在尝试在 R 中进行泊松回归，并且我想获得预测区间。我看到lm 的预测函数给它写'interval="prediction"' 但它不适用于predict.glm！

有人知道是否有办法获得这些预测区间吗？有例子可以打代码吗？

所以我要数一数医院每天接受的人数，我有以下代码：

poisson_reg=glm(NumeberAccept ~ 1 + weekday + month + place + NoConvention + Rain, 
                family = poisson(link = log), data = dataset)
summary(poisson_reg)

现在如果我输入 R predict(poisson_reg, newdata, type="responce") 我可以预测每天的接受次数，但我也需要预测间隔！ 我看到，对于预测调用中的 "lm" 类对象，您可以编写：predict(poisson_reg, newdata, interval="prediction")，它给出了 95% 的预测区间。有没有办法用 "glm" 类的对象获得相同的结果？

【问题讨论】：

查看 ?confint。 stat.ethz.ch/R-manual/R-devel/library/MASS/html/confint.html

这没有那么简单。不像在线性模型的情况下，可以简单地将由于参数不确定性引起的方差添加到估计的误差方差中以获得总体预测方差，这里您可能需要从参数的抽样分布（或自举分布）中进行模拟，然后从条件泊松分布中进行模拟，并收集相关的置信包络。

@本：别以为我跟着你..你能举个例子吗？ @Frank：Confint 不行，因为我需要预测区间而不是置信区间..

【参考方案1】：

这可能更像是一个统计问题而不是编程问题，但是：

从上一个问题中窃取示例数据：

ex <- read.table(
  header=TRUE, text='
Number.Accepted  Weekday    Month   Place
  20    6   8   1
  16    7   8   1
  12    4   8   2
  11    7   1   1
  12    1   4   1
  12    7   10  2
  13    5   6   2
')
ex.glm <- glm(Number.Accepted ~ Weekday + Month + Place,
              family = poisson, data = ex)

我们想要预测区间的数据框：

newdata <- data.frame(Weekday=c(5,6),Month=c(9,9),Place=c(1,1))

类似这样的：

bootSimFun <- function(preddata,fit,data) 
    bdat <- data[sample(seq(nrow(data)),size=nrow(data),replace=TRUE),]
    bfit <- update(fit,data=bdat)
    bpred <- predict(bfit,type="response",newdata=preddata)
    rpois(length(bpred),lambda=bpred)

你也可以使用base R中的replicate()，但是plyr::raply()很方便：

library(plyr)
set.seed(101)
simvals <- raply(500,bootSimFun(preddata=newdata,fit=ex.glm,data=ex))
t(apply(simvals,2,quantile,c(0.025,0.975)))
##    2.5% 97.5%
## 1 7.000    40
## 2 7.475    36

【讨论】：

嗨 Ben，我已经尝试过你的方法，但我有一些问题.. 当我在 bdat 中执行 bootSimFun 时，我需要有我的数据集 Acceptances_2，我希望从中获得数字的预置间隔在接下来的两个月中被接受..所以在我的数据集中我没有列 NumberAccepted 因为我想要该变量的预置。所以你的代码给了我一个错误，因为它想要那个变量..我的推理错了吗？

【参考方案2】：

考虑一下 Zelig 包。请参阅此处的泊松插图——http://rss.acs.unt.edu/Rdoc/library/Zelig/doc/poisson.pdf。

Zelig 有一个统一的方法，不仅可以建模（为此，glm() 及其各种链接函数就足够了），还可以提取和绘制感兴趣的数量。特别是，要对预测范围（而不仅仅是预期范围）进行建模，您必须同时模拟您的系数（系统分量）和误差项（随机分量）。简单地对误差项进行平均，我认为这是 predict.glm() 所做的，将为您提供更窄的预期范围。

Zelig 有一个函数 sim()，它可以模拟系统和随机分量，并输出可用于绘制预测和预期范围的内存对象。它还有一个函数 setx()，如果你想模拟你的解释变量给定值的预测不确定性，你可以在 sim() 之前使用它。看这里——http://rss.acs.unt.edu/Rdoc/library/Zelig/html/setx.html。

这一切都始于这篇论文：http://gking.harvard.edu/files/abs/making-abs.shtml。 Zelig 基本上就是 Clarify 长大后的样子。

【讨论】：

嗨，Gabi，我也试过你的方法..但我不明白一些事情..我做了回归和模拟：reg_zelig=zelig(formula=NumeroAccept~1+weekday+month+ place+NoConvention+Rain,data=acceptance,model="poisson") x.out=setx(reg_zelig, data=acceptance_2) 这里我有数据集acceptance_2，我想在其中预测 NumberAccept 所以我没有那个变量..所以出现错误..如果我输入 NumerAccept=0 它可以工作，但是当我执行 s.out

预期值：E(Y|X) 均值 sd 50% 2.5% 97.5% 335.463 3.076 335.401 329.502 341.582 预测值：Y|X 均值 sd 50% 2.5% 97.5% 335.242 18.453 但它只是指我数据集中的第一个观察结果。为什么？也许我把 NumberAccept=0 放在数据集中是错误的，但我没有那个变量，因为我想要它的预测间隔！！我很抱歉，但我不明白这个预测值是否是我的真的需要..你能帮我理解吗？真的谢谢