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

最近我們被客戶要求撰寫關(guān)于貝葉斯線性回歸的研究報(bào)告，包括一些圖形和統(tǒng)計(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每周平均工作時(shí)間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)？

• 觀察研究

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

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

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

• 7名受訪者每周收入低于300元

summary(wage)##????Min.?1st?Qu.??Median????Mean?3rd?Qu.????Max.?##???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)系。

點(diǎn)擊標(biāo)題查閱往期內(nèi)容

視頻：R語言中的Stan概率編程MCMC采樣的貝葉斯模型

左右滑動查看更多

01

02

03

04

受教育程度和工作小時(shí)數(shù)似乎是工人工資的良好預(yù)測因素。

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

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

簡單的線性回歸

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

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

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

m_wage_iq$coefficients##?(Intercept)??????????iq?##??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è)，任何隨后的推斷都是無效的。

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

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

qqnorm(m_wage_iq$residuals)qqline(m_wage_iq$residuals)

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

##?(Intercept)????????educ?##???146.95244????60.21428summary(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)格為正時(shí)才可能，因?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)的。

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

• (0.00709, 0.01050)

#?從線性模型m_lwage_iq中提取系數(shù)值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分會引起工資的提高

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ī)缺失的。例如，如果所有第一個出生的孩子沒有報(bào)告他們的出生順序，數(shù)據(jù)就不會隨機(jī)丟失。在沒有任何額外信息的情況下，我們將假設(shè)這是合理的，并使用663個完整的觀測值（與原來的935個相反）來擬合模型。Bayesian和frequentist方法都存在于處理缺失數(shù)據(jù)的數(shù)據(jù)集上，但是它們超出了本文的范圍。

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

• 已婚黑人

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

married_black?<-?married_coef*1+black_coef*1 married_black##?[1]?0.09561888

從線性模型的快速總結(jié)中可以看出，自變量的許多系數(shù)在統(tǒng)計(jì)上并不顯著。您可以根據(jù)調(diào)整后的R2選擇變量。本文引入了貝葉斯信息準(zhǔn)則（BIC），這是一種可用于模型選擇的度量。BIC基于模型擬合，同時(shí)根據(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。

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

• feduc

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））找到最佳模型。

#對于AIC，懲罰因子是一個接觸值k。對于step BIC，我們將使用stepAIC函數(shù)并將k設(shè)置為log（n）step(m_lwage_full1,?direction?=?"backward",?k=log(n))##?Residuals:##?????Min??????1Q??Median??????3Q?????Max?##?-172.57??-63.43??-35.43???23.39?1065.78?##?##?Coefficients:##???????????????Estimate?Std.?Error?t?value?Pr(>|t|)????##?(Intercept)?-5546.2382????84.7839?-65.416??<?2e-16?***##?hours???????????1.9072?????0.6548???2.913???0.0037?**?##?tenure?????????-4.1285?????0.9372??-4.405?1.23e-05?***##?lwage?????????951.0113????11.5041??82.667??<?2e-16?***##?---##?Signif.?codes:??0?'***'?0.001?'**'?0.01?'*'?0.05?'.'?0.1?'?'?1##?##?Residual?standard?error:?120.1?on?659?degrees?of?freedom##?Multiple?R-squared:??0.9131,?Adjusted?R-squared:??0.9127?##?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ù)。

bma(lwage?~?.?-wage,?data?=?wage_no_na, ???????????????????prior?=?"BIC",? ???????????????????modelprior?=?uniform())##?##??Marginal?Posterior?Inclusion?Probabilities:?##?Intercept??????hours?????????iq????????kww???????educ??????exper??##???1.00000????0.85540????0.89732????0.34790????0.99887????0.70999??##????tenure????????age???married1?????black1?????south1?????urban1??##???0.70389????0.52468????0.99894????0.34636????0.32029????1.00000??##??????sibs????brthord??????meduc??????feduc??##???0.04152????0.12241????0.57339????0.23274summary(bma_lwage)##???????????P(B?!=?0?|?Y)????model?1???????model?2???????model?3##?Intercept????1.00000000?????1.0000?????1.0000000?????1.0000000##?hours????????0.85540453?????1.0000?????1.0000000?????1.0000000##?iq???????????0.89732383?????1.0000?????1.0000000?????1.0000000##?kww??????????0.34789688?????0.0000?????0.0000000?????0.0000000##?educ?????????0.99887165?????1.0000?????1.0000000?????1.0000000##?exper????????0.70999255?????0.0000?????1.0000000?????1.0000000##?tenure???????0.70388781?????1.0000?????1.0000000?????1.0000000##?age??????????0.52467710?????1.0000?????1.0000000?????0.0000000##?married1?????0.99894488?????1.0000?????1.0000000?????1.0000000##?black1???????0.34636467?????0.0000?????0.0000000?????0.0000000##?south1???????0.32028825?????0.0000?????0.0000000?????0.0000000##?urban1???????0.99999983?????1.0000?????1.0000000?????1.0000000##?sibs?????????0.04152242?????0.0000?????0.0000000?????0.0000000##?brthord??????0.12241286?????0.0000?????0.0000000?????0.0000000##?meduc????????0.57339302?????1.0000?????1.0000000?????1.0000000##?feduc????????0.23274084?????0.0000?????0.0000000?????0.0000000##?BF???????????????????NA?????1.0000?????0.5219483?????0.5182769##?PostProbs????????????NA?????0.0455?????0.0237000?????0.0236000##?R2???????????????????NA?????0.2710?????0.2767000?????0.2696000##?dim??????????????????NA?????9.0000????10.0000000?????9.0000000##?logmarg??????????????NA?-1490.0530?-1490.7032349?-1490.7102938##?????????????????model?4???????model?5##?Intercept?????1.0000000?????1.0000000##?hours?????????1.0000000?????1.0000000##?iq????????????1.0000000?????1.0000000##?kww???????????1.0000000?????0.0000000##?educ??????????1.0000000?????1.0000000##?exper?????????1.0000000?????0.0000000##?tenure????????1.0000000?????1.0000000##?age???????????0.0000000?????1.0000000##?married1??????1.0000000?????1.0000000##?black1????????0.0000000?????1.0000000##?south1????????0.0000000?????0.0000000##?urban1????????1.0000000?????1.0000000##?sibs??????????0.0000000?????0.0000000##?brthord???????0.0000000?????0.0000000##?meduc?????????1.0000000?????1.0000000##?feduc?????????0.0000000?????0.0000000##?BF????????????0.4414346?????0.4126565##?PostProbs?????0.0201000?????0.0188000##?R2????????????0.2763000?????0.2762000##?dim??????????10.0000000????10.0000000##?logmarg???-1490.8707736?-1490.9381880

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

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

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

##????????????????????2.5%????????97.5%??????????beta##?Intercept??6.787648e+00?6.841901e+00??6.8142970694##?hours?????-9.213871e-03?0.000000e+00?-0.0053079979##?iq?????????0.000000e+00?6.259498e-03??0.0037983313##?kww????????0.000000e+00?8.290187e-03??0.0019605787##?educ???????2.247901e-02?6.512858e-02??0.0440707549##?exper??????0.000000e+00?2.100097e-02??0.0100264057##?tenure?????0.000000e+00?1.288251e-02??0.0059357058##?age????????0.000000e+00?2.541561e-02??0.0089659753##?married1???1.173088e-01?3.015797e-01??0.2092940731##?black1????-1.900452e-01?0.000000e+00?-0.0441863361##?south1????-1.021332e-01?1.296992e-05?-0.0221757978##?urban1?????1.374767e-01?2.621311e-01??0.1981221313##?sibs???????0.000000e+00?0.000000e+00??0.0000218455##?brthord???-2.072393e-02?0.000000e+00?-0.0019470674##?meduc??????0.000000e+00?2.279918e-02??0.0086717156##?feduc?????-7.996880e-06?1.558576e-02??0.0025125930##?attr(,"Probability")##?[1]?0.95##?attr(,"class")##?[1]?"confint.bas"

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

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

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

• 年齡

##?##?Call:##?##?##??Marginal?Posterior?Inclusion?Probabilities:?##?Intercept??????hours?????????iq????????kww???????educ??????exper??##????1.0000?????0.8692?????0.9172?????0.3217?????1.0000?????0.9335??##????tenure????????age???married1?????black1?????south1?????urban1??##????0.9980?????0.1786?????0.9999?????0.9761?????0.8149?????1.0000##???????????P(B?!=?0?|?Y)????model?1???????model?2???????model?3##?Intercept?????1.0000000?????1.0000?????1.0000000?????1.0000000##?hours?????????0.8691891?????1.0000?????1.0000000?????1.0000000##?iq????????????0.9171607?????1.0000?????1.0000000?????1.0000000##?kww???????????0.3216992?????0.0000?????1.0000000?????0.0000000##?educ??????????1.0000000?????1.0000?????1.0000000?????1.0000000##?exper?????????0.9334844?????1.0000?????1.0000000?????1.0000000##?tenure????????0.9980015?????1.0000?????1.0000000?????1.0000000##?age???????????0.1786252?????0.0000?????0.0000000?????0.0000000##?married1??????0.9999368?????1.0000?????1.0000000?????1.0000000##?black1????????0.9761347?????1.0000?????1.0000000?????1.0000000##?south1????????0.8148861?????1.0000?????1.0000000?????0.0000000##?urban1????????1.0000000?????1.0000?????1.0000000?????1.0000000##?BF???????????????????NA?????1.0000?????0.5089209?????0.2629789##?PostProbs????????????NA?????0.3311?????0.1685000?????0.0871000##?R2???????????????????NA?????0.2708?????0.2751000?????0.2634000##?dim??????????????????NA????10.0000????11.0000000?????9.0000000##?logmarg??????????????NA?-2275.4209?-2276.0963811?-2276.7565998##?????????????????model?4???????model?5##?Intercept?????1.0000000?????1.0000000##?hours?????????1.0000000?????0.0000000##?iq????????????1.0000000?????1.0000000##?kww???????????0.0000000?????0.0000000##?educ??????????1.0000000?????1.0000000##?exper?????????1.0000000?????1.0000000##?tenure????????1.0000000?????1.0000000##?age???????????1.0000000?????0.0000000##?married1??????1.0000000?????1.0000000##?black1????????1.0000000?????1.0000000##?south1????????1.0000000?????1.0000000##?urban1????????1.0000000?????1.0000000##?BF????????????0.2032217?????0.1823138##?PostProbs?????0.0673000?????0.0604000##?R2????????????0.2737000?????0.2628000##?dim??????????11.0000000?????9.0000000##?logmarg???-2277.0143763?-2277.1229445

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

• 真的

bma_lwage_full##?##?Call:##?##?##??Marginal?Posterior?Inclusion?Probabilities:?##?Intercept??????hours?????????iq????????kww???????educ??????exper??##????1.0000?????0.9792?????0.9505?????0.6671?????0.9998?????0.8951??##????tenure????????age???married1?????black1?????south1?????urban1??##????0.9040?????0.7093?????0.9998?????0.7160?????0.6904?????1.0000??##??????sibs????brthord??????meduc??????feduc??##????0.3939?????0.5329?????0.7575?????0.5360##???????????P(B?!=?0?|?Y)?????model?1?????model?2?????model?3?model?4##?Intercept?????1.0000000??1.00000000??1.00000000??1.00000000??1.0000##?hours?????????0.9791805??1.00000000??1.00000000??1.00000000??1.0000##?iq????????????0.9504649??1.00000000??1.00000000??1.00000000??1.0000##?kww???????????0.6670580??1.00000000??1.00000000??1.00000000??1.0000##?educ??????????0.9998424??1.00000000??1.00000000??1.00000000??1.0000##?exper?????????0.8950911??1.00000000??1.00000000??1.00000000??1.0000##?tenure????????0.9040156??1.00000000??1.00000000??1.00000000??1.0000##?age???????????0.7092839??1.00000000??1.00000000??1.00000000??1.0000##?married1??????0.9997881??1.00000000??1.00000000??1.00000000??1.0000##?black1????????0.7160065??1.00000000??1.00000000??1.00000000??1.0000##?south1????????0.6903763??1.00000000??1.00000000??1.00000000??1.0000##?urban1????????1.0000000??1.00000000??1.00000000??1.00000000??1.0000##?sibs??????????0.3938833??1.00000000??1.00000000??0.00000000??0.0000##?brthord???????0.5329258??1.00000000??1.00000000??1.00000000??0.0000##?meduc?????????0.7575462??1.00000000??1.00000000??1.00000000??1.0000##?feduc?????????0.5359832??1.00000000??0.00000000??1.00000000??0.0000##?BF???????????????????NA??0.01282537??0.06040366??0.04899546??1.0000##?PostProbs????????????NA??0.07380000??0.02320000??0.01880000??0.0126##?R2???????????????????NA??0.29250000??0.29140000??0.29090000??0.2882##?dim??????????????????NA?16.00000000?15.00000000?15.00000000?13.0000##?logmarg??????????????NA?76.00726935?77.55689364?77.34757164?80.3636##?????????????model?5##?Intercept??1.000000##?hours??????1.000000##?iq?????????1.000000##?kww????????1.000000##?educ???????1.000000##?exper??????1.000000##?tenure?????1.000000##?age????????1.000000##?married1???1.000000##?black1?????1.000000##?south1?????1.000000##?urban1?????1.000000##?sibs???????0.000000##?brthord????1.000000##?meduc??????1.000000##?feduc??????0.000000##?BF?????????0.227823##?PostProbs??0.012500##?R2?????????0.289600##?dim???????14.000000##?logmarg???78.884413

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

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次。

set.seed(314)N?=?100000y_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ù)。

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

• (-3.11, 5.13)

quantile(y_5,?probs?=?c(0.025,?0.975))##??????2.5%?????97.5%?##?-3.109585??5.132511

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

?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]?"Intercept"?"hours"?????"iq"????????"kww"???????"educ"????? ##??[6]?"exper"?????"tenure"????"age"???????"married1"??"urban1"??? ##?[11]?"meduc"

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

predict(bma_lwage,?estimator="MPM")##??[1]?"Intercept"?"hours"?????"iq"????????"educ"??????"exper"???? ##??[6]?"tenure"????"age"???????"married1"??"urban1"????"meduc"

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

##?????????[,1]?????##?wage????"1586"???##?hours???"40"?????##?iq??????"127"????##?kww?????"48"?????##?educ????"16"?????##?exper???"16"?????##?tenure??"12"?????##?age?????"37"?????##?married?"1"??????##?black???"0"??????##?south???"0"??????##?urban???"1"??????##?sibs????"4"??????##?brthord?"4"??????##?meduc???"16"?????##?feduc???"16"?????##?lwage???"7.36897"

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

##?????2.5%????97.5%?????pred?##?6.661869?8.056452?7.359160

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

##??????2.5%?????97.5%??????pred?##??782.0108?3154.0780?1570.5169

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

##??????2.5%?????97.5%??????pred?##??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.

本文摘選?《?R語言貝葉斯線性回歸和多元線性回歸構(gòu)建工資預(yù)測模型?》?，點(diǎn)擊“閱讀原文”獲取全文完整資料。

點(diǎn)擊標(biāo)題查閱往期內(nèi)容

Metropolis Hastings采樣和貝葉斯泊松回歸Poisson模型
Matlab用BUGS馬爾可夫區(qū)制轉(zhuǎn)換Markov switching隨機(jī)波動率模型、序列蒙特卡羅SMC、M H采樣分析時(shí)間序列R語言RSTAN MCMC：NUTS采樣算法用LASSO 構(gòu)建貝葉斯線性回歸模型分析職業(yè)聲望數(shù)據(jù)
R語言BUGS序列蒙特卡羅SMC、馬爾可夫轉(zhuǎn)換隨機(jī)波動率SV模型、粒子濾波、Metropolis Hasting采樣時(shí)間序列分析
R語言Metropolis Hastings采樣和貝葉斯泊松回歸Poisson模型
R語言貝葉斯MCMC：用rstan建立線性回歸模型分析汽車數(shù)據(jù)和可視化診斷
R語言貝葉斯MCMC：GLM邏輯回歸、Rstan線性回歸、Metropolis Hastings與Gibbs采樣算法實(shí)例
R語言貝葉斯Poisson泊松-正態(tài)分布模型分析職業(yè)足球比賽進(jìn)球數(shù)
R語言用Rcpp加速M(fèi)etropolis-Hastings抽樣估計(jì)貝葉斯邏輯回歸模型的參數(shù)
R語言邏輯回歸、Naive Bayes貝葉斯、決策樹、隨機(jī)森林算法預(yù)測心臟病
R語言中貝葉斯網(wǎng)絡(luò)（BN）、動態(tài)貝葉斯網(wǎng)絡(luò)、線性模型分析錯頜畸形數(shù)據(jù)
R語言中的block Gibbs吉布斯采樣貝葉斯多元線性回歸
Python貝葉斯回歸分析住房負(fù)擔(dān)能力數(shù)據(jù)集
R語言實(shí)現(xiàn)貝葉斯分位數(shù)回歸、lasso和自適應(yīng)lasso貝葉斯分位數(shù)回歸分析
Python用PyMC3實(shí)現(xiàn)貝葉斯線性回歸模型
R語言用WinBUGS 軟件對學(xué)術(shù)能力測驗(yàn)建立層次（分層）貝葉斯模型
R語言Gibbs抽樣的貝葉斯簡單線性回歸仿真分析
R語言和STAN,JAGS：用RSTAN,RJAG建立貝葉斯多元線性回歸預(yù)測選舉數(shù)據(jù)
R語言基于copula的貝葉斯分層混合模型的診斷準(zhǔn)確性研究
R語言貝葉斯線性回歸和多元線性回歸構(gòu)建工資預(yù)測模型
R語言貝葉斯推斷與MCMC：實(shí)現(xiàn)Metropolis-Hastings 采樣算法示例
R語言stan進(jìn)行基于貝葉斯推斷的回歸模型
R語言中RStan貝葉斯層次模型分析示例
R語言使用Metropolis-Hastings采樣算法自適應(yīng)貝葉斯估計(jì)與可視化
R語言隨機(jī)搜索變量選擇SSVS估計(jì)貝葉斯向量自回歸（BVAR）模型
WinBUGS對多元隨機(jī)波動率模型：貝葉斯估計(jì)與模型比較
R語言實(shí)現(xiàn)MCMC中的Metropolis–Hastings算法與吉布斯采樣
R語言貝葉斯推斷與MCMC：實(shí)現(xiàn)Metropolis-Hastings 采樣算法示例
R語言使用Metropolis-Hastings采樣算法自適應(yīng)貝葉斯估計(jì)與可視化
視頻：R語言中的Stan概率編程MCMC采樣的貝葉斯模型
R語言MCMC:Metropolis-Hastings采樣用于回歸的貝葉斯估計(jì)