R語言中回歸模型預(yù)測的不同類型置信區(qū)間應(yīng)用比較分析
原文鏈接:http://tecdat.cn/?p=13913?

?
?
我們討論了使用程序來獲得預(yù)測的置信區(qū)間的方法。我們將討論線性回歸。
?
> plot(cars)
> reg=lm(dist~speed,data=cars)
> abline(reg,col="red")
> n=nrow(cars)
> x=21
> points(x,predict(reg,newdata= data.frame(speed=x)),pch=19,col="red")

?
我們正在這里做出一個預(yù)測。正如在R課堂上(以及在預(yù)測模型的過程中)所回顧的,當(dāng)我們要為預(yù)測提供一個置信區(qū)間時,建議您為預(yù)測器確定置信區(qū)間(這將取決于預(yù)測誤差)參數(shù)的估計)和潛在值的置信區(qū)間(這也取決于模型誤差,即殘差的離散度)。讓我們從預(yù)測的置信區(qū)間開始:
abline(reg,col="light blue")
points(x,predict(reg,newdata=data.frame(speed=x)),pch=19,col="blue")

藍(lán)色值是可能的預(yù)測,可以通過在我們的觀察數(shù)據(jù)庫中重新采樣獲得。殘差(以及因此的斜率和回歸線的常數(shù)的估計值)的正態(tài)性假設(shè)下的置信區(qū)間(為90%)如下
lines(0:30,U[,2],col="red",lwd=2)
lines(0:30,U[,3],col="red",lwd=2)

?
我們可以在這里比較在500個生成的數(shù)據(jù)集上獲得的值的分布,并比較經(jīng)驗分位數(shù)和假設(shè)正態(tài)性下的分位數(shù),
polygon(c(D$x[I],rev(D$x[I])),c(D$y[I],rep(0,length(I))),col="blue",border=NA)

數(shù)量比較
5% ? ? ?95%
58.63689 70.31281
fit ? ? ?lwr ? ? ?upr
65.00149 59.65934 70.34364
現(xiàn)在,讓我們來看另一種類型的置信區(qū)間,即關(guān)注變量的可能值。這次,除了繪制新樣本和計算預(yù)測值之外,我們還將在每次繪制中添加噪聲,我們獲得可能的值。
points(x,Yx[s],pch=19,col="red")

?
同樣,在這里,我們可以比較(以圖形方式開始)通過重采樣獲得的值,以及在正常情況下獲得的值,
polygon(c(D$x[I],rev(D$x[I])),c(D$y[I],rep(0,length(I))),col="blue",border=NA)

?
從數(shù)字上給出以下比較
5% ? ? ?95%
44.43468 96.01357
fit ? ? ?lwr ? ? ?upr
1 67.63136 45.16967 90.09305
這次右側(cè)略有不對稱。顯然,我們不能假設(shè)高斯殘差,因為正值比負(fù)值大??紤]到數(shù)據(jù)的性質(zhì)(距離不能為負(fù)),這是合理的。
然后,我們開始討論使用回歸模型。
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 3209 4372 4411 4428 4435 4456
[2,] 3367 4659 4696 4720 4730 ? NA
[3,] 3871 5345 5398 5420 ? NA ? NA
[4,] 4239 5917 6020 ? NA ? NA ? NA
[5,] 4929 6794 ? NA ? NA ? NA ? NA
[6,] 5217 ? NA ? NA ? NA ? NA ? NA
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 3209 1163 ? 39 ? 17 ? ?7 ? 21
[2,] 3367 1292 ? 37 ? 24 ? 10 ? NA
[3,] 3871 1474 ? 53 ? 22 ? NA ? NA
[4,] 4239 1678 ?103 ? NA ? NA ? NA
[5,] 4929 1865 ? NA ? NA ? NA ? NA
[6,] 5217 ? NA ? NA ? NA ? NA ? NA
然后,我們可以建立一個數(shù)據(jù)。
> head(base,12)
y ? ai bj
1 ?3209 2000 ?0
2 ?3367 2001 ?0
3 ?3871 2002 ?0
4 ?4239 2003 ?0
5 ?4929 2004 ?0
6 ?5217 2005 ?0
7 ?1163 2000 ?1
8 ?1292 2001 ?1
9 ?1474 2002 ?1
10 1678 2003 ?1
11 1865 2004 ?1
12 ? NA 2005 ?1
> tail(base,12)
y ? ai bj
25 ?7 2000 ?4
26 10 2001 ?4
27 NA 2002 ?4
28 NA 2003 ?4
29 NA 2004 ?4
30 NA 2005 ?4
31 21 2000 ?5
32 NA 2001 ?5
33 NA 2002 ?5
34 NA 2003 ?5
35 NA 2004 ?5
36 NA 2005 ?5
然后,我們可以使用基于?Stavros Christofides的對數(shù)增量支付模型的回歸模型,該模型基于對數(shù)正態(tài)模型,該模型最初由Etienne de Vylder于1978年提出。
Residuals:
Min ? ? ? 1Q ? Median ? ? ? 3Q ? ? ?Max
-0.26374 -0.05681 ?0.00000 ?0.04419 ?0.33014
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) ? ? ? ? 7.9471 ? ? 0.1101 ?72.188 6.35e-15 ***
as.factor(ai)2001 ? 0.1604 ? ? 0.1109 ? 1.447 ?0.17849
as.factor(ai)2002 ? 0.2718 ? ? 0.1208 ? 2.250 ?0.04819 *
as.factor(ai)2003 ? 0.5904 ? ? 0.1342 ? 4.399 ?0.00134 **
as.factor(ai)2004 ? 0.5535 ? ? 0.1562 ? 3.543 ?0.00533 **
as.factor(ai)2005 ? 0.6126 ? ? 0.2070 ? 2.959 ?0.01431 *
as.factor(bj)1 ? ? -0.9674 ? ? 0.1109 ?-8.726 5.46e-06 ***
as.factor(bj)2 ? ? -4.2329 ? ? 0.1208 -35.038 8.50e-12 ***
as.factor(bj)3 ? ? -5.0571 ? ? 0.1342 -37.684 4.13e-12 ***
as.factor(bj)4 ? ? -5.9031 ? ? 0.1562 -37.783 4.02e-12 ***
as.factor(bj)5 ? ? -4.9026 ? ? 0.2070 -23.685 4.08e-10 ***
---
Signif. codes: ?0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.1753 on 10 degrees of freedom
(15 observations deleted due to missingness)
Multiple R-squared: 0.9975, Adjusted R-squared: 0.9949
F-statistic: 391.7 on 10 and 10 DF, ?p-value: 1.338e-11
[,1] ? [,2] [,3] [,4] [,5] [,6]
[1,] 2871.2 1091.3 41.7 18.3 ?7.8 21.3
[2,] 3370.8 1281.2 48.9 21.5 ?9.2 25.0
[3,] 3768.0 1432.1 54.7 24.0 10.3 28.0
[4,] 5181.5 1969.4 75.2 33.0 14.2 38.5
[5,] 4994.1 1898.1 72.5 31.8 13.6 37.1
[6,] 5297.8 2013.6 76.9 33.7 14.5 39.3
> sum(base$py[is.na(base$y)])
[1] 2481.857
我們獲得與通過Chain Ladder方法獲得的結(jié)果略有不同。如Hachemeister和Stanard在1975年所建議的,我們還可以嘗試Poisson回歸(具有對數(shù)鏈接),
Deviance Residuals:
Min ? ? ? 1Q ? Median ? ? ? 3Q ? ? ?Max
-2.3426 ?-0.4996 ? 0.0000 ? 0.2770 ? 3.9355
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) ? ? ? ?8.05697 ? ?0.01551 519.426 ?< 2e-16 ***
as.factor(ai)2001 ?0.06440 ? ?0.02090 ? 3.081 ?0.00206 **
as.factor(ai)2002 ?0.20242 ? ?0.02025 ? 9.995 ?< 2e-16 ***
as.factor(ai)2003 ?0.31175 ? ?0.01980 ?15.744 ?< 2e-16 ***
as.factor(ai)2004 ?0.44407 ? ?0.01933 ?22.971 ?< 2e-16 ***
as.factor(ai)2005 ?0.50271 ? ?0.02079 ?24.179 ?< 2e-16 ***
as.factor(bj)1 ? ?-0.96513 ? ?0.01359 -70.994 ?< 2e-16 ***
as.factor(bj)2 ? ?-4.14853 ? ?0.06613 -62.729 ?< 2e-16 ***
as.factor(bj)3 ? ?-5.10499 ? ?0.12632 -40.413 ?< 2e-16 ***
as.factor(bj)4 ? ?-5.94962 ? ?0.24279 -24.505 ?< 2e-16 ***
as.factor(bj)5 ? ?-5.01244 ? ?0.21877 -22.912 ?< 2e-16 ***
---
Signif. codes: ?0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 46695.269 ?on 20 ?degrees of freedom
Residual deviance: ? ?30.214 ?on 10 ?degrees of freedom
(15 observations deleted due to missingness)
AIC: 209.52
Number of Fisher Scoring iterations: 4
> round(matrix(base$py2,n,n),1)
[,1] ? [,2] [,3] [,4] [,5] [,6]
[1,] 3155.7 1202.1 49.8 19.1 ?8.2 21.0
[2,] 3365.6 1282.1 53.1 20.4 ?8.8 22.4
[3,] 3863.7 1471.8 61.0 23.4 10.1 25.7
[4,] 4310.1 1641.9 68.0 26.1 11.2 28.7
[5,] 4919.9 1874.1 77.7 29.8 12.8 32.7
[6,] 5217.0 1987.3 82.4 31.6 13.6 34.7
>
> sum(base$py2[is.na(base$y)])
[1] 2426.985
該預(yù)測與通過鏈梯方法獲得的估計量一致。Klaus Schmidt和AngelaWünsche于1998年在鏈梯,邊際總和和最大似然估計中建立了帶有最小偏差方法的鏈接。
?