kla*_*air 4 regression r poisson prediction intervals
我已经尝试了这两种方法,但我发现这两种方法都有困难..在我告诉你我对这两种方法的问题之前,我试图更好地解释我的问题。
我有一个数据集“接受”,其中我有医院每天接受的数量,其中包含前面描述的独立变量。医院有三个地方供我们就诊……所以在我的数据集中,我每天有 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
Run Code Online (Sandbox Code Playgroud)
.... 等等...我有直到昨天的数据集,所以最后三行是昨天 2013 年 7 月 29 日的接受。现在我做我的泊松回归:
poisson_reg=glm(NumeberAccept ~ 1 + weekday + month + place + NoConvention + Rain,
family = poisson(link = log), data = acceptances)
Run Code Online (Sandbox Code Playgroud)
现在,对于我的预测,我创建了一个新的数据集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)
Run Code Online (Sandbox Code Playgroud)
现在,如果我输入 R,predict(poisson_reg, newdata, type="responce")我就可以预测每天的接受次数,但我也需要预测间隔!我"lm"在 predict 调用中看到了一个类的对象,你可以写:predict(poisson_reg, newdata, interval="prediction")它给出了 95% 的预测间隔。有没有办法用 class 的对象获得相同的东西"glm"?
这可能更像是一个统计问题而不是一个编程问题,但是:
从上一个问题中窃取示例数据:
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)
Run Code Online (Sandbox Code Playgroud)
我们想要预测区间的数据帧:
newdata <- data.frame(Weekday=c(5,6),Month=c(9,9),Place=c(1,1))
Run Code Online (Sandbox Code Playgroud)
像这样的东西:
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)
}
Run Code Online (Sandbox Code Playgroud)
您也可以replicate()从基础 R使用,但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
Run Code Online (Sandbox Code Playgroud)