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

原文出處：拓端數(shù)據(jù)部落公眾號

工資模型

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

加載包

在本實(shí)驗(yàn)中，我們將使用dplyr包探索數(shù)據(jù)，并使用ggplot2包進(jìn)行數(shù)據(jù)可視化。我們也可以在其中一個練習(xí)中使用MASS包來實(shí)現(xiàn)逐步線性回歸。

我們將在實(shí)驗(yàn)室稍后使用此軟件包中使用BAS.LM來實(shí)現(xiàn)貝葉斯模型。

數(shù)據(jù)

本實(shí)驗(yàn)室將使用的數(shù)據(jù)是在全國935名受訪者中隨機(jī)抽取的。

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

這是觀察研究還是實(shí)驗(yàn)？

1. 觀察研究

探索數(shù)據(jù)

與任何新數(shù)據(jù)集一樣，標(biāo)準(zhǔn)的探索性數(shù)據(jù)分析是一個好的開始。我們將從工資變量開始，因?yàn)樗俏覀兡Ｐ椭械囊蜃兞俊?/p>

1. 關(guān)于工資問題，下列哪種說法是錯誤的？

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

由于工資是我們的因變量，我們想探討其他變量之間的關(guān)系作為預(yù)測。
練習(xí)：排除工資和工齡，選擇另外兩個你認(rèn)為可以很好預(yù)測工資的變量。使用適當(dāng)?shù)膱D來形象化他們與工資的關(guān)系。
受教育程度和工作小時數(shù)似乎是工人工資的良好預(yù)測因素。

1. 


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

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

簡單的線性回歸

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

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

2. geom_point()

這個圖是相當(dāng)雜亂的。雖然智商分?jǐn)?shù)和工資之間可能存在輕微的正線性關(guān)系，但智商充其量只是一個粗略的工資預(yù)測指標(biāo)。我們可以通過擬合一個簡單的線性回歸來量化這一點(diǎn)。

m_wage_iq$coefficients

1. ## (Intercept) ? ? ? ? ?iq

2. ## ?116.991565 ? ?8.303064

## [1] 384.7667

回想一下，在模型下

如果使用?

?和參考先驗(yàn)?

? ?，然后貝葉斯后驗(yàn)均值和標(biāo)準(zhǔn)差分別等于頻數(shù)估計(jì)和標(biāo)準(zhǔn)差。

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

1. 檢驗(yàn)m_wage_iq的殘差。正態(tài)分布誤差的假設(shè)有效嗎?

1. 不，因?yàn)槟Ｐ偷臍埐罘植际怯移摹?/p>

1. qqnorm(m_wage_iq$residuals)

2. qqline(m_wage_iq$residuals)

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

1. ## (Intercept) ? ? ? ?educ

2. ## ? 146.95244 ? ?60.21428

summary(m_wage_educ)$sigma## [1] 382.3203

同樣的結(jié)論是，該線性模型的殘差與?i～N（0，σ2）近似正態(tài)分布，因此可以在該線性模型的基礎(chǔ)上進(jìn)行進(jìn)一步的推斷。

變量轉(zhuǎn)換

適應(yīng)數(shù)據(jù)右偏的一種方法是（自然）對數(shù)變換因變量。請注意，這僅在變量嚴(yán)格為正時才可能，因?yàn)闆]有定義負(fù)值的對數(shù)，并且log（0）=?∞。我們試著用對數(shù)工資作為因變量來擬合一個線性模型。問題4將基于這個對數(shù)轉(zhuǎn)換模型。

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

練習(xí)：檢查該模型的殘差。假設(shè)正態(tài)分布的殘差合理嗎？

基于上述殘差圖，可以假定對數(shù)工資線性模型與iq的正態(tài)分布。

回想一下，給定σ2的α和β的后驗(yàn)分布是正態(tài)的，但略微遵循一個具有n?p?1自由度的t分布。在這種情況下，p=1，因?yàn)橹巧淌俏覀兡Ｐ椭形ㄒ坏膶?shù)工資預(yù)測因子。因此，α和β的后驗(yàn)概率都遵循933自由度的t分布，因?yàn)閐f非常大，這些分布實(shí)際上是近似正態(tài)的。

1. 在參考先驗(yàn)p（α，β，σ2）∞1/σ2下，給出β的95%后驗(yàn)置信區(qū)間，即IQ系數(shù)。

1. (0.00709, 0.01050)

1. 


2. # 從線性模型m_lwage_iq中提取系數(shù)值

3. 


4. qnorm(c(0.025, 0.975), mean = iq_mean_estimate, sd=iq_sd)

## [1] 0.007103173 0.010511141

練習(xí)：智商系數(shù)很小，這是意料之中的，因?yàn)橹巧谭謹(jǐn)?shù)提高一分很難對工資產(chǎn)生很高的倍增效應(yīng)。使系數(shù)更易于解釋的一種方法是在將智商放入模型之前將其標(biāo)準(zhǔn)化。從這個新模型來看，智商提高1個標(biāo)準(zhǔn)差（15分）估計(jì)工資會增加多少百分比？

智商是用scale函數(shù)標(biāo)準(zhǔn)化的，智商提高15分會引起工資的提高

1. 


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

## [1] 8.767568

多元線性回歸

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

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

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

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

1. 已婚黑人

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

1. 


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

3. 


4. married_black

## [1] 0.09561888

從線性模型的快速總結(jié)中可以看出，自變量的許多系數(shù)在統(tǒng)計(jì)上并不顯著。您可以根據(jù)調(diào)整后的R2選擇變量。本文引入了貝葉斯信息準(zhǔn)則（BIC），這是一種可用于模型選擇的度量。BIC基于模型擬合，同時根據(jù)樣本大小按比例懲罰參數(shù)個數(shù)。我們可以使用以下命令計(jì)算全線性模型的BIC：

BIC(m_lwage_full)## [1] 586.3732

我們可以比較完整模型和簡化模型的BIC。讓我們試著從模型中刪除出生順序。為了確保觀測值保持不變，可以將數(shù)據(jù)集指定為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.4031BIC(m_lwage_feduc)## [1] 580.9743BIC(m_lwage_meduc)## [1] 582.3722

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

1. 


2. #對于AIC，懲罰因子是一個接觸值k。對于step BIC，我們將使用stepAIC函數(shù)并將k設(shè)置為log（n）

3. 


4. step(m_lwage_full1, direction = "backward", k=log(n))

5. 


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

貝葉斯模型平均

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

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命令為我們提供了每個變量的后驗(yàn)?zāi)Ｐ桶怕屎妥羁赡艿哪Ｐ汀＠?，模型中包含小時數(shù)的后驗(yàn)概率為0.855。此外，后驗(yàn)概率為0.0455的最可能模型包括截距、工作時間、智商、教育程度、年齡、婚姻狀況、城市生活狀況和母親教育程度。雖然0.0455的后驗(yàn)概率聽起來很小，但它比分配給它的統(tǒng)一先驗(yàn)概率大得多，因?yàn)橛?16個可能的模型。
在模型平均法下，還可以可視化系數(shù)的后驗(yàn)分布。我們將智商系數(shù)的后驗(yàn)分布繪制如下。

1. 


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

我們還可以為這些系數(shù)提供95%的置信區(qū)間：

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，我們將使用簡化的數(shù)據(jù)集，其中不包括兄弟姐妹數(shù)量、出生順序和父母教育。

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

1. 基于這個簡化的數(shù)據(jù)集，根據(jù)貝葉斯模型平均，下列哪一個變量的邊際后驗(yàn)包含概率最低？

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. 判斷:包含所有變量的樸素模型的后驗(yàn)概率大于0.5。(系數(shù)使用Zellner-Siow零先驗(yàn)，模型使用β二項(xiàng)(1,1)先驗(yàn))

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

練習(xí)：用數(shù)據(jù)集繪制年齡系數(shù)的后驗(yàn)分布圖。

1. 


2. plot(coef_bma_wage_red, ask = FALSE)

預(yù)測

貝葉斯統(tǒng)計(jì)的一個主要優(yōu)點(diǎn)是預(yù)測和預(yù)測的概率解釋。大部分貝葉斯預(yù)測都是使用模擬技術(shù)來完成的。這通常應(yīng)用于回歸建模中，盡管我們將通過一個僅包含截距項(xiàng)的示例來進(jìn)行分析。

假設(shè)你觀察到y(tǒng)的四個數(shù)值觀測值，分別為2、2、0和0，樣本均值y′=1，樣本方差s2=4/3。假設(shè)y～N（μ，σ2），在參考先驗(yàn)p（μ，σ2）～1/σ2下，我們的后驗(yàn)概率變?yōu)?/p>

以樣本均值為中心

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

為了得到y(tǒng)5的預(yù)測分布，我們可以先從σ2的后驗(yàn)點(diǎn)模擬，然后再從μ模擬y5。我們對y5年的預(yù)測結(jié)果將來自一項(xiàng)新的觀測結(jié)果的后驗(yàn)預(yù)測分布。下面的示例從y5的后驗(yàn)預(yù)測分布中提取100,000次。

1. set.seed(314)

2. N = 100000

3. 


4. y_5 = rnorm(N, mu, sqrt(sigma2))

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

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

1. 估計(jì)一個新的觀測y595%置信區(qū)間

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

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

1. 


2. plot(den_of_y)

BAS預(yù)測

在BAS中，用貝葉斯模型平均法構(gòu)造預(yù)測區(qū)間是通過仿真實(shí)現(xiàn)的，而在模型選擇的情況下，用預(yù)測區(qū)間進(jìn)行精確推理往往是可行的。

回到工資數(shù)據(jù)集，讓我們找到最佳預(yù)測模型下的預(yù)測值，即預(yù)測值最接近BMA和相應(yīng)的后驗(yàn)標(biāo)準(zhǔn)差的模型。

predict(bma_lwage, estimator="BPM")

1. ## ?[1] "Intercept" "hours" ? ? "iq" ? ? ? ?"kww" ? ? ? "educ"

2. ## ?[6] "exper" ? ? "tenure" ? ?"age" ? ? ? "married1" ?"urban1"

3. ## [11] "meduc"

我們可以將其與我們之前發(fā)現(xiàn)的最高概率模型和中位概率模型（MPM）進(jìn)行比較

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。
練習(xí)：使用簡化數(shù)據(jù)，最佳預(yù)測模型、中位概率模型和最高后驗(yàn)概率模型中包含哪些協(xié)變量？
讓我們來看看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"

可以得到預(yù)測對數(shù)工資的95%可信區(qū)間

1. ## ? ? 2.5% ? ?97.5% ? ? pred

2. ## 6.661869 8.056452 7.359160

換算成工資，我們可以將區(qū)間指數(shù)化。

1. ## ? ? ?2.5% ? ? 97.5% ? ? ?pred

2. ## ?782.0108 3154.0780 1570.5169

獲得一個95%的預(yù)測區(qū)間的工資。
如果使用BMA，區(qū)間是

1. ## ? ? ?2.5% ? ? 97.5% ? ? ?pred

2. ## ?733.3446 2989.2076 1494.9899

練習(xí)：使用簡化后的數(shù)據(jù)，為BPM下預(yù)測工資最高的個人構(gòu)建95%的預(yù)測區(qū)間。

參考文獻(xiàn)

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

最受歡迎的見解

1.matlab使用貝葉斯優(yōu)化的深度學(xué)習(xí)

2.matlab貝葉斯隱馬爾可夫hmm模型實(shí)現(xiàn)

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

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

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

6.Python用PyMC3實(shí)現(xiàn)貝葉斯線性回歸模型

7.R語言使用貝葉斯 層次模型進(jìn)行空間數(shù)據(jù)分析

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

9.matlab貝葉斯隱馬爾可夫hmm模型實(shí)現(xiàn)