按年份划分的泊松回归总计的预测区间

Question

我正在使用调查设计进行一项研究，以预测 2017 年至 2030 年数据集中每年的手术总数。我目前正在使用泊松回归来创建这些预测估计值。

  svydesign(
    data = four,
    strata = ~NIS_STRATUMnew,
    ids = ~NISIDnew,
    weights = ~DISCWTnew,
    nest = T
  )

我能够成功地使用以下方法为我们拥有数据的年份（2017 年之前）生成估计值和准确的置信区间：

kneetotal1<- svyby(~kneePJI, ~YEAR, design = mydesign, FUN = svytotal, vartype = "ci")

但是，如果我尝试直接使用我的数据进行预测 (1) 或使用 svyby (2) 结果做出预测：

(1)

AdjPoissonKnee <- svyglm(kneePJI ~ YEAR, family = poisson(), design = mydesign)
years <- data.frame(YEAR = 2018:2030)
predictHip <- predict(AdjPoissonHip, newdata = data.frame(YEAR = 2018:2030), type = "response", se.fit =TRUE, interval = "predict") 

（这似乎只是产生了每个人的程序可能性。我不确定如何为此产生年度累积总和以及置信区间）

(2)

futureyears <- data.table(YEAR = 2018:2030)
AdjPoissonKnee <- glm(kneePJI ~ YEAR, family = poisson(), data=kneetotal1)
kneepredict <- predict(object = AdjPoissonKnee, newdata=futureyears, type = "response")

数据集非常大，因此可能会导致预测间隔非常窄的问题。 Lmk 如果有办法我可以通过添加一些数据的 sn-p 来提供帮助。

示例输出 (2) 与 2002-2017 年的 svyby 相结合：

kneetotal1<- svyby(~kneePJI, ~YEAR, design = mydesign, FUN = svytotal, vartype = "ci")

> kneetotal1 
YEAR    kneePJI 
2002    8205.194 
2003    9648.019 
2004    10362.076 
2005    11771.961 
2006    12099.399 
2007    12993.681 
2008    15438.871 
2009    14303.562 
2010    16051.562 
2011    17158.166 
2012    17055.006 
2013    18064.991 
2014    19080.006 
2015    19429.988 
2016    18070.002 
2017    18399.995

AdjPoissonKnee <- glm(kneePJI ~ YEAR, family = poisson(), data=kneetotal1)
kneese <- predict(object = AdjPoissonKnee, newdata=futureyears, type = "response", interval = "prediction", se.fit =TRUE)

> kneese
$fit
       1        2        3        4        5        6        7        8        9       10       11       12       13 
22107.84 23232.34 24414.04 25655.84 26960.81 28332.15 29773.25 31287.64 32879.07 34551.44 36308.87 38155.70 40096.46 

$se.fit
        1         2         3         4         5         6         7         8         9        10        11        12        13 
 87.14235 100.68424 115.63713 132.05885 150.02228 169.61303 190.92798 214.07449 239.16993 266.34157 295.72666 327.47263 361.73741 

> stderrKnee <- kneese$se.fit
> predfitKnee <- kneese$fit
> reserrsqKnee <- predfitKnee^2*summary(AdjPoissonKnee)$dispersion^2
> prederrsqKnee <- stderrKnee^2 + reserrsqKnee
> lower_pi <- predfitKnee - 1.96 * sqrt(prederrsqKnee)
> upper_pi <- predfitKnee + 1.96 * sqrt(prederrsqKnee)
> lower_pi
        1         2         3         4         5         6         7         8         9        10        11        12        13 
-21223.86 -22303.48 -23438.01 -24630.27 -25883.19 -27199.86 -28583.52 -30037.57 -31565.61 -33171.39 -34858.88 -36632.22 -38495.80
> upper_pi
        1         2         3         4         5         6         7         8         9        10        11        12        13 
 65439.55  68768.16  72266.09  75941.96  79804.81  83864.16  88130.01  92612.86  97323.74 102274.27 107476.62 112943.62 118688.72 

Answer 1

一个问题是predict.svyglm 没有有 interval="predict" 选项（您不会收到错误，因为 R 中的方法需要接受和忽略它们没有的参数理解，以便继承工作）

但是，它不会有太大的不同。平均值为 20000 的 Poisson 变量的标准差为 $\sqrt{20000}$，即大约 141，将 $1.96\times 141$ 添加到区间半角不会有很大的不同。

很可能您的数据相对于泊松分布过度分散：如果您这样做 summary(AdjPoissonKnee)（对于 svyglm），部分输出是估计的分散参数，我猜它会更大大于 1，这意味着预测分布不应该是泊松。

这是为什么predict.svyglm目前不做预测区间的一部分。考虑泊松分布。它是离散的。预测区间应该是多少？

可以创建一个预测区间来估计分散参数，然后通过具有正确方差的正态分布来近似预测分布。 Peter Ellis 为普通 glms 讨论了这种方法here，他的代码可以修改：你会这样做（未经测试，因为我没有你的数据）

AdjPoissonKnee <- svyglm(kneePJI ~ YEAR, 
       family = poisson(), design = mydesign)
years <- data.frame(YEAR = 2018:2030)
predictHip <- predict(AdjPoissonKnee, 
      newdata = data.frame(YEAR = 2018:2030), 
       type = "response", se.fit =TRUE) 

stderr<-SE(predictHip)
predfit<-coef(predictHip)

reserrsq<- predfit^2*summary(AdjPoissonKnee)$dispersion^2

prederrsq <- stderr^2 + reserrsq

lower_pi <- predfit - 1.96 * sqrt(prederrsq)
upper_pi <- predfit + 1.96 * sqrt(prederrsq)