在R中对泊松回归进行预测区间

Question

我已经尝试了这两种方法，但是在使用它们时遇到了一些困难... 在我谈论这两种方法的问题之前，我会更好地解释一下我的问题。

我有一个名为“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中计算出下两个月的接受量的预测区间！因此，第一行将是今天接受量，最后一行将是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="response")，我可以得到每天接受人数的预测结果，但我还需要预测间隔！我发现，在"lm"类对象的预测调用中，您可以写入：predict(poisson_reg, newdata, interval="prediction")，它会给出95%的预测间隔。有没有办法在"glm"类对象中获得相同的结果？

- klair

---

请查看?confint。http://stat.ethz.ch/R-manual/R-devel/library/MASS/html/confint.html - Frank P.

这并不是那么简单。与线性模型情况不同的是，在线性模型中，人们可以将由于参数不确定性引起的方差简单地添加到估计误差方差中以获得总体预测方差，但在这里，您可能需要从参数的抽样分布（或自助法分布）进行模拟，然后从条件泊松分布进行模拟，并收集相关置信区间。 - Ben Bolker

@ Ben：不太明白你的意思，能给个例子吗？ @ Frank：Confint 不行，因为我需要预测区间而不是置信区间。 - klair

Answer 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)
}

你也可以使用 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 Bolker

---

嗨Ben，我尝试了你的方法，但是我有一些问题。当我在bdat中执行bootSimFun时，我需要从我的数据集acceptances_2中获取预测区间，以便获得接下来两个月的接受数量。因此，在我的数据集中，我没有NumberAccepted列，因为我想要该变量的预测值。所以你的代码会给我一个错误，因为它需要那个变量。我的推理错了吗？ - klair

Answer 2

考虑Zelig包。在这里看泊松分布的vignette-- 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

---

嗨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<-sim(reg_zelig, x = x.out)时，它给出的结果是： - klair

预期值：E(Y|X) 平均值 标准差 中位数 2.5%分位数 97.5%分位数 335.463 3.076 335.401 329.502 341.582预测值：Y|X 平均值 标准差 中位数 2.5%分位数 97.5%分位数 335.242 18.45 335 300 373但这只是针对数据集中第一个观测值的情况。为什么呢？也许我在数据集中放置了NumberAccept=0，但实际上我没有该变量，因为我想要针对它进行预测区间！！我很抱歉，但我不确定这些预测值是否真正满足我的需求。您能帮助我理解吗？非常感谢。 - klair