带区间的预测图答案

【问题标题】：prediction plot with intervals带区间的预测图
【发布时间】：2017-11-10 05:03:41
【问题描述】：

我试图创建一个预测区间，其中包括我的模型的预测（负二项式）。型号为：

Call:
glm.nb(formula = TOT.N ~ D.PARK + OPEN.L + L.WAT.C + sqrt(L.P.ROAD), 
    init.theta = 4.979895131, link = log)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-3.08218  -0.70494  -0.09268   0.55575   1.67860  

Coefficients:
                 Estimate Std. Error z value Pr(>|z|)    
(Intercept)     4.032e+00  3.363e-01  11.989  < 2e-16 ***
D.PARK         -1.154e-04  1.061e-05 -10.878  < 2e-16 ***
OPEN.L         -1.085e-02  3.122e-03  -3.475  0.00051 ***
L.WAT.C         1.597e-01  7.852e-02   2.034  0.04195 *  
sqrt(L.P.ROAD)  4.924e-01  3.101e-01   1.588  0.11231    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for Negative Binomial(4.9799) family taken to be 1)

    Null deviance: 197.574  on 51  degrees of freedom
Residual deviance:  51.329  on 47  degrees of freedom
AIC: 383.54

Number of Fisher Scoring iterations: 1

我首先创建了一个观察到的与预测的图（如下）。从这个图中我会说预测数据似乎与实际数据非常吻合。

然后我试图用置信区间制作预测图。为此，我决定让变量OPEN.L 变化，而其他变量保持不变 - 均值。我使用的代码如下：

  varying OPEN.L
  minOPEN.L <- min(OPEN.L)
  maxOPEN.L <- max(OPEN.L)
  grid <- seq(minOPEN.L, maxOPEN.L, 1)
  mean.D.PARK <- mean(D.PARK)

  new <- data.frame(D.PARK = mean.D.PARK, OPEN.L = grid, L.WAT.C = mean.L.WAT.C, L.P.ROAD = mean.L.P.ROAD)

  confidece.kills <- predict(final.model, new, se = T, interval = "confidence")
  predict.kills <- predict(final.model, new, se = T, interval = "prediction")

  par(mfrow=c(1, 2), pty="m")
  matplot(grid, predict.kills$fit ,lty=c(1,2,2),type="l",lwd=3,
          xlab="OPEN.L",ylab="TOT.N",
          cex.lab=1.5,cex.axis=1.3)

剧情上没什么可看的（下）：

dput(head(road.data, 55))

    dput(head(road.data, 55))
structure(list(TOT.N = c(22L, 14L, 65L, 55L, 88L, 104L, 49L, 
66L, 26L, 47L, 35L, 55L, 44L, 30L, 33L, 29L, 34L, 64L, 76L, 32L, 
34L, 32L, 35L, 22L, 34L, 25L, 18L, 14L, 14L, 7L, 7L, 17L, 10L, 
3L, 6L, 5L, 2L, 3L, 2L, 2L, 7L, 3L, 5L, 4L, 7L, 12L, 7L, 14L, 
10L, 4L, 11L, 3L), OPEN.L = c(22.684, 24.657, 30.121, 50.277, 
43.609, 31.385, 24.81, 56.228, 48.735, 15.633, 9.999, 39.942, 
10.382, 2.507, 0.738, 15.725, 43.866, 45.102, 39.46, 19.988, 
13.369, 6.848, 2.946, 3.219, 3.218, 34.168, 22.839, 7.258, 8.513, 
23.394, 26.945, 71.436, 62.203, 82.391, 97.574, 94.947, 89.294, 
68.779, 62.173, 67.834, 67.618, 83.357, 70.684, 30.907, 26.687, 
9.571, 26.687, 16.478, 26.365, 39.609, 33.511, 24.438), MONT.S = c(0, 
0, 0.258, 1.783, 2.431, 0, 0, 0, 1.108, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 5.235, 3.658, 5.049, 0.224, 9.426, 0, 0, 0, 0, 0, 0.763, 
7.134, 0, 0, 1.039, 4.326, 0, 0, 0, 0, 0, 0, 0, 1.455, 0, 0, 
0, 4.347, 0, 1.376, 4.347, 1.796, 0, 0.259), POLIC = c(4.811, 
2.224, 1.946, 0.625, 0.791, 0.054, 0.022, 11.263, 1.238, 0.119, 
0.024, 0, 0.038, 0, 0, 0, 0.06, 0.125, 1.7, 0, 0.308, 0.364, 
0.013, 0, 0, 0, 0.529, 0.313, 0.063, 0.202, 0, 0, 0, 0, 0.206, 
0.259, 0.278, 0.812, 0.03, 0.018, 0.206, 0.375, 0.086, 0.05, 
0.06, 0, 0.06, 0, 0, 0.044, 1.861, 0.151), D.PARK = c(250.214, 
741.179, 1240.08, 1739.885, 2232.13, 2724.089, 3215.511, 3709.401, 
4206.477, 4704.176, 5202.328, 5700.669, 6199.342, 6698.151, 7187.762, 
7668.833, 8152.155, 8633.224, 9101.411, 9573.578, 10047.63, 10523.939, 
11002.496, 11482.896, 11976.232, 12470.968, 12968.285, 13465.914, 
13961.321, 14432.954, 14904.995, 15377.983, 15854.389, 16335.936, 
16810.109, 17235.045, 17673.064, 18167.269, 18656.949, 19149.507, 
19645.717, 20141.987, 20640.729, 21138.903, 21631.542, 22119.102, 
22613.647, 23113.45, 23606.088, 24046.886, 24444.874, 24884.803
), SHRUB = c(0.406, 0.735, 0.474, 0.607, 0.173, 0.325, 0.055, 
0.092, 1.744, 0, 0.67, 0.783, 0, 0.178, 0, 0, 0.094, 0.107, 0.702, 
0.827, 1.025, 0, 0.01, 0.012, 0.088, 0.02, 0.087, 0.116, 0.062, 
0, 0, 0.033, 0.133, 0.047, 0.077, 0.182, 0.067, 0.208, 0.063, 
0.122, 0.038, 0.095, 0, 0.02, 0.064, 0.137, 0.064, 0.214, 0.14, 
0.622, 0, 0.18), WAT.RES = c(0.043, 0.182, 0.453, 0.026, 0, 0.039, 
0.114, 0.224, 0.177, 0, 6.309, 2.26, 0.137, 0, 0, 0.402, 0.077, 
0.042, 0, 0.479, 0.36, 0, 0.078, 0, 0, 0, 0.188, 0, 0, 0, 0.213, 
2.452, 0.061, 0, 0, 0, 0.284, 0.579, 0.215, 0, 0, 0, 0.127, 0, 
0.198, 0.473, 0.198, 0, 0, 0, 0, 0.319), L.WAT.C = c(583, 1419, 
2005, 1924, 2167, 2391, 1165, 2428, 2416, 211, 292, 650, 1896, 
2194, 1375, 0, 1655, 1702, 2721, 1694, 1192, 589, 476, 345, 1621, 
1023, 357, 0, 0, 7, 878, 883, 1921, 1479, 1237, 1898, 3951, 1931, 
1365, 591, 868, 1198, 2334, 3525, 3087, 2444, 3087, 3934, 2214, 
2122, 1290, 2471), L.P.ROAD = c(1975, 1761, 1250, 666, 653, 1309, 
685, 677, 664, 654, 696, 678, 652, 665, 655, 627, 1159, 2201, 
2290, 1617, 866, 640, 620, 645, 853, 1370, 631, 603, 609, 605, 
1374, 685, 594, 1075, 595, 676, 684, 733, 1739, 891, 730, 652, 
668, 645, 602, 571, 602, 953, 765, 1578, 2960, 1407), D.WAT.COUR = c(735, 
134.052, 269.029, 48.751, 126.102, 344.444, 95.133, 243.23, 187.084, 
236.004, 15.184, 118.865, 332.257, 28.498, 168.818, 560, 104.839, 
204.943, 256.812, 566.152, 689.823, 694, 300, 132.934, 253.305, 
34.119, 515.233, 825, 1165, 1025, 754.938, 585, 137.112, 80.916, 
35.426, 43, 31.01, 290.029, 68.496, 405, 785, 257, 118.579, 237.041, 
45.832, 44.744, 120.855, 24.313, 178.837, 21.336, 111.764, 225.514
)), .Names = c("TOT.N", "OPEN.L", "MONT.S", "POLIC", "D.PARK", 
"SHRUB", "WAT.RES", "L.WAT.C", "L.P.ROAD", "D.WAT.COUR"), row.names = c(NA, 
-52L), class = c("tbl_df", "tbl", "data.frame"))

如何获得具有相应预测区间的预测图？

任何输入 - cmets、建设性批评、tipps - 都会得到认可。谢谢

【问题讨论】：

标签： r plot prediction

【解决方案1】：

您可以像这样自己估算间隔：

attach(road.data)
library("MASS")
final.model <- glm.nb(formula = TOT.N ~ D.PARK + OPEN.L + L.WAT.C + sqrt(L.P.ROAD),
       init.theta = 4.979895131, link = log)

minOPEN.L <- min(OPEN.L)
maxOPEN.L <- max(OPEN.L)
grid <- seq(minOPEN.L, maxOPEN.L, 1)
mean.D.PARK <- mean(D.PARK)
mean.L.WAT.C <- mean(L.WAT.C)
mean.L.P.ROAD <- mean(L.P.ROAD)

new <- data.frame(D.PARK = mean.D.PARK, OPEN.L = grid, L.WAT.C = mean.L.WAT.C, L.P.ROAD = mean.L.P.ROAD)

predict.kills <- predict(final.model, new, se = T, interval = "prediction")

alpha <- 0.05 ## you want a 95% interval
z <- qnorm(1 - alpha / 2)
lower <- predict.kills$fit - z * predict.kills$se.fit
upper <- predict.kills$fit + z * predict.kills$se.fit

par(mfrow=c(1, 2), pty="m")
matplot(grid, cbind(predict.kills$fit, lower, upper) ,lty=c(1,2,2),type="l",lwd=3,
        xlab="OPEN.L",ylab="log(TOT.N)",
        cex.lab=1.5,cex.axis=1.3)

matplot(grid, exp(cbind(predict.kills$fit, lower, upper)) ,lty=c(1,2,2),type="l",lwd=3,
        xlab="OPEN.L",ylab="TOT.N",
        cex.lab=1.5,cex.axis=1.3)

你会得到这样的结果：

【讨论】：

@一致性，感谢您的帮助。我之前已经计算了置信区间（在我的情况下，它们应该比你的更窄，因为我正在拟合负二项式模型）：`估计 2.5 % 97.5 %（截距）4.0319970271 3.4094551079 4.655064e+00 D.PARK - 0.000.0.0.0170224260 - 0.192080e + 00 + 00/12101155517031-0.011.0155513031-0.087656504我怎样才能插入lower \ upper？
init.theta = 4.979895131 代表过度分散常数对吗？它到底是做什么的？
@Danka 那么您已经计算了所有系数估计的 95% 区间，并且您想使用这些区间来获得 95% 的区间来进行响应吗？据我所知，没有简单的方法可以做到这一点。您可以将这些 95% 区间的结果结合起来得到最终的 95% 区间，但这太宽了；要么你需要知道这些预测变量之间的关系。
@Danka 我认为theta是过度分散参数，如果不是，则与过度分散参数密切相关。过度离散是，如果你使用glm模型，均值和方差之间会有一些关系（在泊松情况下，均值等于方差），但这可能不是我们想要的，所以我们使用过度离散参数来处理有了这个，再回到泊松的例子，我们可以让方差 = theta * mean，theta 将是过度离散参数。
那你为什么选择 theta 为 4.979895131？