拓端數據tecdat|R語言貝葉斯線性回歸和多元線性回歸構建工資預測模型

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

 

工資模型

在勞動經濟學領域，收入和工資的研究為從性別歧視到高等教育等問題提供了見解。在本文中，我們將分析橫斷面工資數據，以期在實踐中使用貝葉斯方法，如BIC和貝葉斯模型來構建工資的預測模型。

加載包

在本實驗中，我們將使用dplyr包探索數據，並使用ggplot2包進行數據可視化。我們也可以在其中一個練習中使用MASS包來實現逐步線性回歸。

我們將在實驗室稍后使用此軟件包中使用BAS.LM來實現貝葉斯模型。

數據

本實驗室將使用的數據是在全國935名受訪者中隨機抽取的。

變量	描述
wage	周收入
hours	每周平均工作時間
IQ	智商得分
kww	工作知識分數
educ	受教育年限
exper	工作經驗
tenure	在現任雇主工作多年
age	年齡
married	= 1，如果已婚
black	= 1（如果為黑人）
south	= 1，如果住在南部
urban	= 1，如果居住在都市中
sibs	兄弟姐妹數
brthord	出生順序
meduc	母親的教育程度
feduc	父親的學歷
lwage	工資的自然對數

這是觀察研究還是實驗？

1. 觀察研究

 

探索數據

與任何新數據集一樣，標准的探索性數據分析是一個好的開始。我們將從工資變量開始，因為它是我們模型中的因變量。

1. 關於工資問題，下列哪種說法是錯誤的？
   1. 7名受訪者每周收入低於300元

summary(wage)

1.  
   ## Min. 1st Qu. Median Mean 3rd Qu. Max.
2.  
   ## 115.0 669.0 905.0 957.9 1160.0 3078.0

由於工資是我們的因變量，我們想探討其他變量之間的關系作為預測。
練習：排除工資和工齡，選擇另外兩個你認為可以很好預測工資的變量。使用適當的圖來形象化他們與工資的關系。
受教育程度和工作小時數似乎是工人工資的良好預測因素。

1.  
    
2.  
   ggplot(data = wage, aes(y=wage, x=exper))+geom_point()

ggplot(data = wage, aes(y=wage, x=educ))+geom_point()

簡單的線性回歸

對於我們在數據中看到的工資差異，一個可能的解釋是，更聰明的人賺更多的錢。下圖顯示了周工資和智商得分之間的散點圖。

1.  
   ggplot(data = wage, aes(x = iq, y = wage)) +
2.  
   geom_point()

這個圖是相當雜亂的。雖然智商分數和工資之間可能存在輕微的正線性關系，但智商充其量只是一個粗略的工資預測指標。我們可以通過擬合一個簡單的線性回歸來量化這一點。

m_wage_iq$coefficients

1.  
   ## (Intercept) iq
2.  
   ## 116.991565 8.303064

## [1] 384.7667

回想一下，在模型下

如果使用  和參考先驗    ，然后貝葉斯后驗均值和標准差分別等於頻數估計和標准差。

貝葉斯模型規范假設誤差正態分布且方差為常數。與頻率法一樣，我們通過檢查模型的殘差分布來檢驗這一假設。如果殘差是高度非正態或偏態的，則違反了假設，任何隨后的推斷都是無效的。

1. 檢驗m_wage_iq的殘差。正態分布誤差的假設有效嗎?
   1. 不，因為模型的殘差分布是右偏的。

1.  
   qqnorm(m_wage_iq$residuals)
2.  
   qqline(m_wage_iq$residuals)

練習：重新調整模型，這次使用educ（教育）作為自變量。你對上一個練習的回答有變化嗎？

1.  
   ## (Intercept) educ
2.  
   ## 146.95244 60.21428

summary(m_wage_educ)$sigma

## [1] 382.3203

同樣的結論是，該線性模型的殘差與ϵi∼N（0，σ2）近似正態分布，因此可以在該線性模型的基礎上進行進一步的推斷。

變量轉換

適應數據右偏的一種方法是（自然）對數變換因變量。請注意，這僅在變量嚴格為正時才可能，因為沒有定義負值的對數，並且log（0）=−∞。我們試着用對數工資作為因變量來擬合一個線性模型。問題4將基於這個對數轉換模型。

m_lwage_iq = lm(lwage ~ iq, data = wage)

練習：檢查該模型的殘差。假設正態分布的殘差合理嗎？

基於上述殘差圖，可以假定對數工資線性模型與iq的正態分布。

回想一下，給定σ2的α和β的后驗分布是正態的，但略微遵循一個具有n−p−1自由度的t分布。在這種情況下，p=1，因為智商是我們模型中唯一的對數工資預測因子。因此，α和β的后驗概率都遵循933自由度的t分布，因為df非常大，這些分布實際上是近似正態的。

1. 在參考先驗p（α，β，σ2）∞1/σ2下，給出β的95%后驗置信區間，即IQ系數。
   1. (0.00709, 0.01050)

1.  
    
2.  
   # 從線性模型m_lwage_iq中提取系數值
3.  
    
4.  
   qnorm(c(0.025, 0.975), mean = iq_mean_estimate, sd=iq_sd)

## [1] 0.007103173 0.010511141

練習：智商系數很小，這是意料之中的，因為智商分數提高一分很難對工資產生很高的倍增效應。使系數更易於解釋的一種方法是在將智商放入模型之前將其標准化。從這個新模型來看，智商提高1個標准差（15分）估計工資會增加多少百分比？

智商是用scale函數標准化的，智商提高15分會引起工資的提高

1.  
    
2.  
   coef(summary(m_lwage_scaled_iq))["siq", "Estimate"]*15+coef(summary(m_lwage_scaled_iq))["(Intercept)", "Estimate"]

## [1] 8.767568

多元線性回歸

很明顯，工資可以用很多預測因素來解釋，比如經驗、教育程度和智商。我們可以在回歸模型中包含所有相關的協變量，試圖盡可能多地解釋工資變化。

lm中的.的使用告訴R在模型中包含所有協變量，然后用-wage進一步修改，然后從模型中排除工資變量。
默認情況下，lm函數執行完整的案例分析，因此它會刪除一個或多個預測變量中缺少（NA）值的觀察值。

由於這些缺失的值，我們必須做一個額外的假設，以便我們的推論是有效的。換句話說，我們的數據必須是隨機缺失的。例如，如果所有第一個出生的孩子沒有報告他們的出生順序，數據就不會隨機丟失。在沒有任何額外信息的情況下，我們將假設這是合理的，並使用663個完整的觀測值（與原來的935個相反）來擬合模型。Bayesian和frequentist方法都存在於處理缺失數據的數據集上，但是它們超出了本課程的范圍。

1. 從這個模型來看，誰賺得更多：已婚的黑人還是單身的非黑人？
   1. 已婚黑人

與單一非黑人男子相比，所有其他平等的，已婚的黑人將獲得以下乘數。

1.  
    
2.  
   married_black <- married_coef*1+black_coef*1
3.  
    
4.  
   married_black

## [1] 0.09561888

從線性模型的快速總結中可以看出，自變量的許多系數在統計上並不顯著。您可以根據調整后的R2選擇變量。本文引入了貝葉斯信息准則（BIC），這是一種可用於模型選擇的度量。BIC基於模型擬合，同時根據樣本大小按比例懲罰參數個數。我們可以使用以下命令計算全線性模型的BIC：

BIC(m_lwage_full)

## [1] 586.3732

我們可以比較完整模型和簡化模型的BIC。讓我們試着從模型中刪除出生順序。為了確保觀測值保持不變，可以將數據集指定為na.omit(wage)，它只包含沒有缺失值的觀測值。

m_lwage_nobrthord = lm(lwage ~ . -wage -brthord, data = na.omit(wage))

## [1] 582.4815

如您所見，從回歸中刪除出生順序會減少BIC，我們試圖通過選擇模型來最小化BIC。

1. 從完整模型中消除哪個變量得到最低的BIC?
   1. feduc

1.  
    
2.  
   BIC(m_lwage_sibs)

## [1] 581.4031

BIC(m_lwage_feduc)

## [1] 580.9743

BIC(m_lwage_meduc)

## [1] 582.3722

練習：R有一個函數stepAIC，它將在模型空間中向后運行，刪除變量直到BIC不再降低。它以一個完整的模型和一個懲罰參數k作為輸入。根據BIC（在這種情況下k=log（n）k=log（n））找到最佳模型。

1.  
    
2.  
   #對於AIC，懲罰因子是一個接觸值k。對於step BIC，我們將使用stepAIC函數並將k設置為log（n）
3.  
    
4.  
   step(m_lwage_full1, direction = "backward", k=log(n))
5.  
    

1.  
    
2.  
   ## Residuals:
3.  
   ## Min 1Q Median 3Q Max
4.  
   ## -172.57 -63.43 -35.43 23.39 1065.78
5.  
   ##
6.  
   ## Coefficients:
7.  
   ## Estimate Std. Error t value Pr(>|t|)
8.  
   ## (Intercept) -5546.2382 84.7839 -65.416 < 2e-16 ***
9.  
   ## hours 1.9072 0.6548 2.913 0.0037 **
10.  
    ## tenure -4.1285 0.9372 -4.405 1.23e-05 ***
11.  
    ## lwage 951.0113 11.5041 82.667 < 2e-16 ***
12.  
    ## ---
13.  
    ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
14.  
    ##
15.  
    ## Residual standard error: 120.1 on 659 degrees of freedom
16.  
    ## Multiple R-squared: 0.9131, Adjusted R-squared: 0.9127
17.  
    ## F-statistic: 2307 on 3 and 659 DF, p-value: < 2.2e-16

貝葉斯模型平均

通常，幾個模型都是同樣可信的，只選擇一個模型忽略了選擇模型中包含的變量所涉及的固有不確定性。解決這一問題的一種方法是實現貝葉斯模型平均（Bayesian model averaging，BMA），即對多個模型進行平均，從新數據中獲得系數的后驗值和預測值。我們可以使用它來實現BMA或選擇模型。我們首先將BMA應用於工資數據。

1.  
   bma(lwage ~ . -wage, data = wage_no_na,
2.  
   prior = "BIC",
3.  
   modelprior = uniform())

1.  
    
2.  
   ##
3.  
   ## Marginal Posterior Inclusion Probabilities:
4.  
   ## Intercept hours iq kww educ exper
5.  
   ## 1.00000 0.85540 0.89732 0.34790 0.99887 0.70999
6.  
   ## tenure age married1 black1 south1 urban1
7.  
   ## 0.70389 0.52468 0.99894 0.34636 0.32029 1.00000
8.  
   ## sibs brthord meduc feduc
9.  
   ## 0.04152 0.12241 0.57339 0.23274

summary(bma_lwage)

1.  
   ## P(B != 0 | Y) model 1 model 2 model 3
2.  
   ## Intercept 1.00000000 1.0000 1.0000000 1.0000000
3.  
   ## hours 0.85540453 1.0000 1.0000000 1.0000000
4.  
   ## iq 0.89732383 1.0000 1.0000000 1.0000000
5.  
   ## kww 0.34789688 0.0000 0.0000000 0.0000000
6.  
   ## educ 0.99887165 1.0000 1.0000000 1.0000000
7.  
   ## exper 0.70999255 0.0000 1.0000000 1.0000000
8.  
   ## tenure 0.70388781 1.0000 1.0000000 1.0000000
9.  
   ## age 0.52467710 1.0000 1.0000000 0.0000000
10.  
    ## married1 0.99894488 1.0000 1.0000000 1.0000000
11.  
    ## black1 0.34636467 0.0000 0.0000000 0.0000000
12.  
    ## south1 0.32028825 0.0000 0.0000000 0.0000000
13.  
    ## urban1 0.99999983 1.0000 1.0000000 1.0000000
14.  
    ## sibs 0.04152242 0.0000 0.0000000 0.0000000
15.  
    ## brthord 0.12241286 0.0000 0.0000000 0.0000000
16.  
    ## meduc 0.57339302 1.0000 1.0000000 1.0000000
17.  
    ## feduc 0.23274084 0.0000 0.0000000 0.0000000
18.  
    ## BF NA 1.0000 0.5219483 0.5182769
19.  
    ## PostProbs NA 0.0455 0.0237000 0.0236000
20.  
    ## R2 NA 0.2710 0.2767000 0.2696000
21.  
    ## dim NA 9.0000 10.0000000 9.0000000
22.  
    ## logmarg NA -1490.0530 -1490.7032349 -1490.7102938
23.  
    ## model 4 model 5
24.  
    ## Intercept 1.0000000 1.0000000
25.  
    ## hours 1.0000000 1.0000000
26.  
    ## iq 1.0000000 1.0000000
27.  
    ## kww 1.0000000 0.0000000
28.  
    ## educ 1.0000000 1.0000000
29.  
    ## exper 1.0000000 0.0000000
30.  
    ## tenure 1.0000000 1.0000000
31.  
    ## age 0.0000000 1.0000000
32.  
    ## married1 1.0000000 1.0000000
33.  
    ## black1 0.0000000 1.0000000
34.  
    ## south1 0.0000000 0.0000000
35.  
    ## urban1 1.0000000 1.0000000
36.  
    ## sibs 0.0000000 0.0000000
37.  
    ## brthord 0.0000000 0.0000000
38.  
    ## meduc 1.0000000 1.0000000
39.  
    ## feduc 0.0000000 0.0000000
40.  
    ## BF 0.4414346 0.4126565
41.  
    ## PostProbs 0.0201000 0.0188000
42.  
    ## R2 0.2763000 0.2762000
43.  
    ## dim 10.0000000 10.0000000
44.  
    ## logmarg -1490.8707736 -1490.9381880

輸出model對象和summary命令為我們提供了每個變量的后驗模型包含概率和最可能的模型。例如，模型中包含小時數的后驗概率為0.855。此外，后驗概率為0.0455的最可能模型包括截距、工作時間、智商、教育程度、年齡、婚姻狀況、城市生活狀況和母親教育程度。雖然0.0455的后驗概率聽起來很小，但它比分配給它的統一先驗概率大得多，因為有216個可能的模型。
在模型平均法下，還可以可視化系數的后驗分布。我們將智商系數的后驗分布繪制如下。

1.  
    
2.  
   plot(coef_lwage, subset = c(3,13), ask=FALSE)

我們還可以為這些系數提供95%的置信區間：

1.  
   ## 2.5% 97.5% beta
2.  
   ## Intercept 6.787648e+00 6.841901e+00 6.8142970694
3.  
   ## hours -9.213871e-03 0.000000e+00 -0.0053079979
4.  
   ## iq 0.000000e+00 6.259498e-03 0.0037983313
5.  
   ## kww 0.000000e+00 8.290187e-03 0.0019605787
6.  
   ## educ 2.247901e-02 6.512858e-02 0.0440707549
7.  
   ## exper 0.000000e+00 2.100097e-02 0.0100264057
8.  
   ## tenure 0.000000e+00 1.288251e-02 0.0059357058
9.  
   ## age 0.000000e+00 2.541561e-02 0.0089659753
10.  
    ## married1 1.173088e-01 3.015797e-01 0.2092940731
11.  
    ## black1 -1.900452e-01 0.000000e+00 -0.0441863361
12.  
    ## south1 -1.021332e-01 1.296992e-05 -0.0221757978
13.  
    ## urban1 1.374767e-01 2.621311e-01 0.1981221313
14.  
    ## sibs 0.000000e+00 0.000000e+00 0.0000218455
15.  
    ## brthord -2.072393e-02 0.000000e+00 -0.0019470674
16.  
    ## meduc 0.000000e+00 2.279918e-02 0.0086717156
17.  
    ## feduc -7.996880e-06 1.558576e-02 0.0025125930
18.  
    ## attr(,"Probability")
19.  
    ## [1] 0.95
20.  
    ## attr(,"class")
21.  
    ## [1] "confint.bas"

對於問題7-8，我們將使用簡化的數據集，其中不包括兄弟姐妹數量、出生順序和父母教育。

wage_red = wage %>% dplyr::select(-sibs, -brthord, -meduc, -feduc)

1. 基於這個簡化的數據集，根據貝葉斯模型平均，下列哪一個變量的邊際后驗包含概率最低？
   1. 年齡

1.  
   ##
2.  
   ## Call:
3.  
   ##
4.  
   ##
5.  
   ## Marginal Posterior Inclusion Probabilities:
6.  
   ## Intercept hours iq kww educ exper
7.  
   ## 1.0000 0.8692 0.9172 0.3217 1.0000 0.9335
8.  
   ## tenure age married1 black1 south1 urban1
9.  
   ## 0.9980 0.1786 0.9999 0.9761 0.8149 1.0000

1.  
   ## P(B != 0 | Y) model 1 model 2 model 3
2.  
   ## Intercept 1.0000000 1.0000 1.0000000 1.0000000
3.  
   ## hours 0.8691891 1.0000 1.0000000 1.0000000
4.  
   ## iq 0.9171607 1.0000 1.0000000 1.0000000
5.  
   ## kww 0.3216992 0.0000 1.0000000 0.0000000
6.  
   ## educ 1.0000000 1.0000 1.0000000 1.0000000
7.  
   ## exper 0.9334844 1.0000 1.0000000 1.0000000
8.  
   ## tenure 0.9980015 1.0000 1.0000000 1.0000000
9.  
   ## age 0.1786252 0.0000 0.0000000 0.0000000
10.  
    ## married1 0.9999368 1.0000 1.0000000 1.0000000
11.  
    ## black1 0.9761347 1.0000 1.0000000 1.0000000
12.  
    ## south1 0.8148861 1.0000 1.0000000 0.0000000
13.  
    ## urban1 1.0000000 1.0000 1.0000000 1.0000000
14.  
    ## BF NA 1.0000 0.5089209 0.2629789
15.  
    ## PostProbs NA 0.3311 0.1685000 0.0871000
16.  
    ## R2 NA 0.2708 0.2751000 0.2634000
17.  
    ## dim NA 10.0000 11.0000000 9.0000000
18.  
    ## logmarg NA -2275.4209 -2276.0963811 -2276.7565998
19.  
    ## model 4 model 5
20.  
    ## Intercept 1.0000000 1.0000000
21.  
    ## hours 1.0000000 0.0000000
22.  
    ## iq 1.0000000 1.0000000
23.  
    ## kww 0.0000000 0.0000000
24.  
    ## educ 1.0000000 1.0000000
25.  
    ## exper 1.0000000 1.0000000
26.  
    ## tenure 1.0000000 1.0000000
27.  
    ## age 1.0000000 0.0000000
28.  
    ## married1 1.0000000 1.0000000
29.  
    ## black1 1.0000000 1.0000000
30.  
    ## south1 1.0000000 1.0000000
31.  
    ## urban1 1.0000000 1.0000000
32.  
    ## BF 0.2032217 0.1823138
33.  
    ## PostProbs 0.0673000 0.0604000
34.  
    ## R2 0.2737000 0.2628000
35.  
    ## dim 11.0000000 9.0000000
36.  
    ## logmarg -2277.0143763 -2277.1229445

1. 判斷:包含所有變量的朴素模型的后驗概率大於0.5。(系數使用Zellner-Siow零先驗，模型使用β二項(1,1)先驗)
   1. 真的

bma_lwage_full

1.  
   ##
2.  
   ## Call:
3.  
   ##
4.  
   ##
5.  
   ## Marginal Posterior Inclusion Probabilities:
6.  
   ## Intercept hours iq kww educ exper
7.  
   ## 1.0000 0.9792 0.9505 0.6671 0.9998 0.8951
8.  
   ## tenure age married1 black1 south1 urban1
9.  
   ## 0.9040 0.7093 0.9998 0.7160 0.6904 1.0000
10.  
    ## sibs brthord meduc feduc
11.  
    ## 0.3939 0.5329 0.7575 0.5360

1.  
   ## P(B != 0 | Y) model 1 model 2 model 3 model 4
2.  
   ## Intercept 1.0000000 1.00000000 1.00000000 1.00000000 1.0000
3.  
   ## hours 0.9791805 1.00000000 1.00000000 1.00000000 1.0000
4.  
   ## iq 0.9504649 1.00000000 1.00000000 1.00000000 1.0000
5.  
   ## kww 0.6670580 1.00000000 1.00000000 1.00000000 1.0000
6.  
   ## educ 0.9998424 1.00000000 1.00000000 1.00000000 1.0000
7.  
   ## exper 0.8950911 1.00000000 1.00000000 1.00000000 1.0000
8.  
   ## tenure 0.9040156 1.00000000 1.00000000 1.00000000 1.0000
9.  
   ## age 0.7092839 1.00000000 1.00000000 1.00000000 1.0000
10.  
    ## married1 0.9997881 1.00000000 1.00000000 1.00000000 1.0000
11.  
    ## black1 0.7160065 1.00000000 1.00000000 1.00000000 1.0000
12.  
    ## south1 0.6903763 1.00000000 1.00000000 1.00000000 1.0000
13.  
    ## urban1 1.0000000 1.00000000 1.00000000 1.00000000 1.0000
14.  
    ## sibs 0.3938833 1.00000000 1.00000000 0.00000000 0.0000
15.  
    ## brthord 0.5329258 1.00000000 1.00000000 1.00000000 0.0000
16.  
    ## meduc 0.7575462 1.00000000 1.00000000 1.00000000 1.0000
17.  
    ## feduc 0.5359832 1.00000000 0.00000000 1.00000000 0.0000
18.  
    ## BF NA 0.01282537 0.06040366 0.04899546 1.0000
19.  
    ## PostProbs NA 0.07380000 0.02320000 0.01880000 0.0126
20.  
    ## R2 NA 0.29250000 0.29140000 0.29090000 0.2882
21.  
    ## dim NA 16.00000000 15.00000000 15.00000000 13.0000
22.  
    ## logmarg NA 76.00726935 77.55689364 77.34757164 80.3636
23.  
    ## model 5
24.  
    ## Intercept 1.000000
25.  
    ## hours 1.000000
26.  
    ## iq 1.000000
27.  
    ## kww 1.000000
28.  
    ## educ 1.000000
29.  
    ## exper 1.000000
30.  
    ## tenure 1.000000
31.  
    ## age 1.000000
32.  
    ## married1 1.000000
33.  
    ## black1 1.000000
34.  
    ## south1 1.000000
35.  
    ## urban1 1.000000
36.  
    ## sibs 0.000000
37.  
    ## brthord 1.000000
38.  
    ## meduc 1.000000
39.  
    ## feduc 0.000000
40.  
    ## BF 0.227823
41.  
    ## PostProbs 0.012500
42.  
    ## R2 0.289600
43.  
    ## dim 14.000000
44.  
    ## logmarg 78.884413

練習：用數據集繪制年齡系數的后驗分布圖。

1.  
    
2.  
   plot(coef_bma_wage_red, ask = FALSE)

預測

貝葉斯統計的一個主要優點是預測和預測的概率解釋。大部分貝葉斯預測都是使用模擬技術來完成的。這通常應用於回歸建模中，盡管我們將通過一個僅包含截距項的示例來進行分析。

假設你觀察到y的四個數值觀測值，分別為2、2、0和0，樣本均值y′=1，樣本方差s2=4/3。假設y∼N（μ，σ2），在參考先驗p（μ，σ2）∼1/σ2下，我們的后驗概率變為

以樣本均值為中心

其中a=（n-1）/2和b=s2（n-1）/2=2。

為了得到y5的預測分布，我們可以先從σ2的后驗點模擬，然后再從μ模擬y5。我們對y5年的預測結果將來自一項新的觀測結果的后驗預測分布。下面的示例從y5的后驗預測分布中提取100,000次。

1.  
   set.seed(314)
2.  
   N = 100000
3.  
    
4.  
   y_5 = rnorm(N, mu, sqrt(sigma2))

我們可以通過觀察模擬數據直方圖的平滑版本，查看預測分布的估計值。

新觀測的95%中心置信區間為在這種情況下，L是0.025分位數，U是0.975分位數。我們可以使用分位數函數來獲得這些值，從而找到tracy5的0.025和0.975的樣本分位數。

1. 估計一個新的觀測y595%置信區間
   1. (-3.11, 5.13)

1.  
    
2.  
    
3.  
   quantile(y_5, probs = c(0.025, 0.975))

1.  
   ## 2.5% 97.5%
2.  
   ## -3.109585 5.132511

練習：在上面的簡單例子中，可以使用積分來分析計算后驗預測值。在這種情況下，它是一個具有3個自由度（n−1）的t分布。繪制y的經驗密度和t分布的實際密度。它們之間有什么比較？

1.  
    
2.  
   plot(den_of_y)

BAS預測

在BAS中，用貝葉斯模型平均法構造預測區間是通過仿真實現的，而在模型選擇的情況下，用預測區間進行精確推理往往是可行的。

回到工資數據集，讓我們找到最佳預測模型下的預測值，即預測值最接近BMA和相應的后驗標准差的模型。

predict(bma_lwage, estimator="BPM")

1.  
   ## [1] "Intercept" "hours" "iq" "kww" "educ"
2.  
   ## [6] "exper" "tenure" "age" "married1" "urban1"
3.  
   ## [11] "meduc"

我們可以將其與我們之前發現的最高概率模型和中位概率模型（MPM）進行比較

predict(bma_lwage, estimator="MPM")

1.  
   ## [1] "Intercept" "hours" "iq" "educ" "exper"
2.  
   ## [6] "tenure" "age" "married1" "urban1" "meduc"

MPM除了包含HPM的所有變量外，還包含exper，而BPM除了MPM中的所有變量外，還包含kwh。
練習：使用簡化數據，最佳預測模型、中位概率模型和最高后驗概率模型中包含哪些協變量？
讓我們來看看BPM模型中哪些特征會影響最高工資。

1.  
   ## [,1]
2.  
   ## wage "1586"
3.  
   ## hours "40"
4.  
   ## iq "127"
5.  
   ## kww "48"
6.  
   ## educ "16"
7.  
   ## exper "16"
8.  
   ## tenure "12"
9.  
   ## age "37"
10.  
    ## married "1"
11.  
    ## black "0"
12.  
    ## south "0"
13.  
    ## urban "1"
14.  
    ## sibs "4"
15.  
    ## brthord "4"
16.  
    ## meduc "16"
17.  
    ## feduc "16"
18.  
    ## lwage "7.36897"

可以得到預測對數工資的95%可信區間

1.  
   ## 2.5% 97.5% pred
2.  
   ## 6.661869 8.056452 7.359160

換算成工資，我們可以將區間指數化。

1.  
   ## 2.5% 97.5% pred
2.  
   ## 782.0108 3154.0780 1570.5169

獲得一個95%的預測區間的工資。
如果使用BMA，區間是

1.  
   ## 2.5% 97.5% pred
2.  
   ## 733.3446 2989.2076 1494.9899

練習：使用簡化后的數據，為BPM下預測工資最高的個人構建95%的預測區間。

參考文獻

Wooldridge, Jeffrey. 2000. Introductory Econometrics- A Modern Approach. South-Western College Publishing. 

---

最受歡迎的見解

1.matlab使用貝葉斯優化的深度學習

2.matlab貝葉斯隱馬爾可夫hmm模型實現

3.R語言Gibbs抽樣的貝葉斯簡單線性回歸仿真

4.R語言中的block Gibbs吉布斯采樣貝葉斯多元線性回歸

5.R語言中的Stan概率編程MCMC采樣的貝葉斯模型

6.Python用PyMC3實現貝葉斯線性回歸模型

7.R語言使用貝葉斯 層次模型進行空間數據分析

8.R語言隨機搜索變量選擇SSVS估計貝葉斯向量自回歸（BVAR）模型

9.matlab貝葉斯隱馬爾可夫hmm模型實現