拓端數據tecdat|R語言基於Bootstrap的線性回歸預測置信區間估計方法

原文鏈接：http://tecdat.cn/?p=21625 

 

我們知道參數的置信區間的計算，這些都服從一定的分布(t分布、正態分布），因此在標准誤前乘以相應的t分值或Z分值。但如果我們找不到合適的分布時，就無法計算置信區間了嗎？幸運的是，有一種方法幾乎可以用於計算各種參數的置信區間，這就是Bootstrap 法。

本文使用BOOTSTRAP來獲得預測的置信區間。我們將在線性回歸基礎上討論。

1.  
    
2.  
   > reg=lm(dist~speed,data=cars)
3.  
   > points(x,predict(reg,newdata= data.frame(speed=x)))

 

這是一個單點預測。當我們想給預測一個置信區間時，預測的置信區間取決於參數估計誤差。

預測置信區間

讓我們從預測的置信區間開始

1.  
    
2.  
   > for(s in 1:500){
3.  
   + indice=sample(1:n,size=n,
4.  
   + replace=TRUE)
5.  
   + points(x,predict(reg,newdata=data.frame(speed=x)),pch=19,col="blue")
6.  
    
7.  
    

藍色值是通過在我們的觀測數據庫中重新取樣獲得的可能預測值。值得注意的是，在殘差正態性假設下（回歸線的斜率和常數估計值），置信區間（90%）如下所示：

predict(reg,interval ="confidence",

 

在這里，我們可以比較500個生成數據集上的值分布，並將經驗分位數與正態假設下的分位數進行比較，

1.  
   > hist(Yx,proba=TRUE
2.  
   > boxplot(Yx,horizontal=TRUE
3.  
   > polygon(c( x ,rev(x I]))))

 

可以看出，經驗分位數與正態假設下的分位數是可以比較的。

1.  
   > quantile(Yx,c(.05,.95))
2.  
   5% 95%
3.  
   58.63689 70.31281
4.  
   + level=.9,newdata=data.frame(speed=x))
5.  
   fit lwr upr
6.  
   1 65.00149 59.65934 70.34364

感興趣變量的可能值

現在讓我們看看另一種類型的置信區間，關於感興趣變量的可能值。這一次，除了提取新樣本和計算預測外，我們還將在每次繪制時添加噪聲，以獲得可能的值。

1.  
   > for(s in 1:500){
2.  
   + indice=sample(1:n,size=n,
3.  
   + base=cars[indice,]
4.  
   + erreur=residuals(reg)
5.  
   + predict(reg,newdata=data.frame(speed=x))+E

 

在這里，我們可以（首先以圖形方式）比較通過重新取樣獲得的值和在正態假設下獲得的值，

1.  
   > hist(Yx,proba=TRUE)
2.  
   > boxplot(Yx) abline(v=U[2:3)
3.  
   > polygon(c(D$x[I,rev(D$x[I])

 

數值上給出了以下比較

1.  
   > quantile(Yx,c(.05,.95))
2.  
   5% 95%
3.  
   44.43468 96.01357
4.  
   U=predict(reg,interval ="prediction"
5.  
   fit lwr upr
6.  
   1 67.63136 45.16967 90.09305

這一次，右側有輕微的不對稱。顯然，我們不能假設高斯殘差，因為有更大的正值，而不是負值。考慮到數據的性質，這是有意義的（制動距離不能是負數）。

然后開始討論在供應中使用回歸模型。為了獲得具有獨立性，有人認為必須使用增量付款的數據，而不是累計付款。

可以創建一個數據庫，解釋變量是行和列。

1.  
   > base=data.frame(
2.  
   + y
3.  
    
4.  
   > head(base,12)
5.  
   y ai bj
6.  
   1 3209 2000 0
7.  
   2 3367 2001 0
8.  
   3 3871 2002 0
9.  
   4 4239 2003 0
10.  
    5 4929 2004 0
11.  
    6 5217 2005 0
12.  
    7 1163 2000 1
13.  
    8 1292 2001 1
14.  
    9 1474 2002 1
15.  
    10 1678 2003 1
16.  
    11 1865 2004 1
17.  
    12 NA 2005 1

然后，我們可以從基於對數增量付款數據的回歸模型開始，該模型基於對數正態模型

1.  
    
2.  
    
3.  
   Coefficients:
4.  
   Estimate Std. Error t value Pr(>|t|)
5.  
   (Intercept) 7.9471 0.1101 72.188 6.35e-15 ***
6.  
   as.factor(ai)2001 0.1604 0.1109 1.447 0.17849
7.  
   as.factor(ai)2002 0.2718 0.1208 2.250 0.04819 *
8.  
   as.factor(ai)2003 0.5904 0.1342 4.399 0.00134 **
9.  
   as.factor(ai)2004 0.5535 0.1562 3.543 0.00533 **
10.  
    as.factor(ai)2005 0.6126 0.2070 2.959 0.01431 *
11.  
    as.factor(bj)1 -0.9674 0.1109 -8.726 5.46e-06 ***
12.  
    as.factor(bj)2 -4.2329 0.1208 -35.038 8.50e-12 ***
13.  
    as.factor(bj)3 -5.0571 0.1342 -37.684 4.13e-12 ***
14.  
    as.factor(bj)4 -5.9031 0.1562 -37.783 4.02e-12 ***
15.  
    as.factor(bj)5 -4.9026 0.2070 -23.685 4.08e-10 ***
16.  
    ---
17.  
    Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
18.  
     
19.  
    Residual standard error: 0.1753 on 10 degrees of freedom
20.  
    (15 observations deleted due to missingness)
21.  
    Multiple R-squared: 0.9975, Adjusted R-squared: 0.9949
22.  
    F-statistic: 391.7 on 10 and 10 DF, p-value: 1.338e-11
23.  
     
24.  
    >
25.  
    exp(predict(reg1,
26.  
    + newdata=base)+summary(reg1)$sigma^2/2)
27.  
     
28.  
    [,1] [,2] [,3] [,4] [,5] [,6]
29.  
    [1,] 2871.2 1091.3 41.7 18.3 7.8 21.3
30.  
    [2,] 3370.8 1281.2 48.9 21.5 9.2 25.0
31.  
    [3,] 3768.0 1432.1 54.7 24.0 10.3 28.0
32.  
    [4,] 5181.5 1969.4 75.2 33.0 14.2 38.5
33.  
    [5,] 4994.1 1898.1 72.5 31.8 13.6 37.1
34.  
    [6,] 5297.8 2013.6 76.9 33.7 14.5 39.3
35.  
     
36.  
    > sum(py[is.na(y)])
37.  
    [1] 2481.857

這與鏈式梯度法的結果略有不同，但仍然具有可比性。我們也可以嘗試泊松回歸（用對數鏈接）

1.  
   glm(y~
2.  
   + as.factor(ai)+
3.  
   + as.factor(bj),data=base,
4.  
   + family=poisson)
5.  
    
6.  
    
7.  
   Coefficients:
8.  
   Estimate Std. Error z value Pr(>|z|)
9.  
   (Intercept) 8.05697 0.01551 519.426 < 2e-16 ***
10.  
    as.factor(ai)2001 0.06440 0.02090 3.081 0.00206 **
11.  
    as.factor(ai)2002 0.20242 0.02025 9.995 < 2e-16 ***
12.  
    as.factor(ai)2003 0.31175 0.01980 15.744 < 2e-16 ***
13.  
    as.factor(ai)2004 0.44407 0.01933 22.971 < 2e-16 ***
14.  
    as.factor(ai)2005 0.50271 0.02079 24.179 < 2e-16 ***
15.  
    as.factor(bj)1 -0.96513 0.01359 -70.994 < 2e-16 ***
16.  
    as.factor(bj)2 -4.14853 0.06613 -62.729 < 2e-16 ***
17.  
    as.factor(bj)3 -5.10499 0.12632 -40.413 < 2e-16 ***
18.  
    as.factor(bj)4 -5.94962 0.24279 -24.505 < 2e-16 ***
19.  
    as.factor(bj)5 -5.01244 0.21877 -22.912 < 2e-16 ***
20.  
    ---
21.  
    Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
22.  
     
23.  
    (Dispersion parameter for poisson family taken to be 1)
24.  
     
25.  
    Null deviance: 46695.269 on 20 degrees of freedom
26.  
    Residual deviance: 30.214 on 10 degrees of freedom
27.  
    (15 observations deleted due to missingness)
28.  
    AIC: 209.52
29.  
     
30.  
    Number of Fisher Scoring iterations: 4
31.  
     
32.  
    > predict(reg2,
33.  
    newdata=base,type="response")
34.  
     
35.  
    > sum(py2[is.na(y)])
36.  
    [1] 2426.985

預測結果與鏈式梯度法得到的估計值吻合。克勞斯·施密特（Klaus Schmidt）和安吉拉·溫什（Angela Wünsche）於1998年在鏈式梯度法、邊際和最大似然估計中建立了與最小偏差方法的聯系。

---

最受歡迎的見解

1.R語言多元Logistic邏輯回歸 應用案例

2.面板平滑轉移回歸(PSTR)分析案例實現

3.matlab中的偏最小二乘回歸（PLSR）和主成分回歸（PCR）

4.R語言泊松Poisson回歸模型分析案例

5.R語言回歸中的Hosmer-Lemeshow擬合優度檢驗

6.r語言中對LASSO回歸，Ridge嶺回歸和Elastic Net模型實現

7.在R語言中實現Logistic邏輯回歸

8.python用線性回歸預測股票價格

9.R語言如何在生存分析與Cox回歸中計算IDI，NRI指標