【視訊】線性回歸中的貝葉斯推斷與R語言預測勞工工資資料|資料分享

在這個視訊中，我們轉向簡單線性回歸中的​​貝葉斯​​推斷。 我們将使用一個參照先驗分布，它提供了頻率主義解決方案和貝葉斯答案之間的聯系。 然後在R語言中用貝葉斯線性回歸、貝葉斯模型平均 (BMA)來預測勞工工資資料 （檢視文末了解資料擷取方式）。

視訊：線性回歸中的貝葉斯推斷與R語言預測勞工工資資料案例

貝葉斯推斷線性回歸與R語言預測勞工工資資料

，時長09:58

為了說明這些想法，我們将使用一個例子來預測身體脂肪。 獲得準确的體脂測量值是很昂貴的，而且在家裡也不容易做到。 

相反，可以使用現成的測量值（如腹圍）來預測體脂百分比的預測模型是很容易使用的，而且價格低廉。

圖中顯示了從水下稱重得到的體脂百分比和多名男性的腹圍。 

為了預測體脂，我們可以從一個線性回歸開始，在散點圖中加入了回歸線。

這條線的估計斜率為0.63，截距約為-39%。 每增加一厘米，我們預計體脂會增加0.63%。 負的截距過程作為一個實體模型沒有意義，但預測一個腰圍為0厘米的男性也沒有意義。 盡管如此，這個線性回歸可能是一個準确的近似值，用于預測目的。 随着測量結果在這個人群的觀察範圍内。

我們的最佳估計線使用α和β的通常平方估計值來獲得拟合值和預測值。

提供拟合誤差估計的殘差是觀察值和預測值之間的差異，用于診斷以及估計sigma平方。 這是通過平均平方誤差，即平方誤差之和除以自由度。

記住，剩餘自由度是樣本量減去模型中的回歸系數。

模型和先驗

讓我們從基準的角度來看看估計的情況。 我們将從與OLS相同的模型開始，但有一個額外的假設，即誤差是正态分布的常數變體。

我們經常加上這個假設，以獲得具有常數平方的置信區間。 是以，在頻率主義回歸中經常使用的任何診斷圖都可以在這裡用來檢查這個假設。

對于共轭分析，先驗是以西格瑪為條件的回歸系數的雙變量正态分布。Sigma在這裡提供了響應的機關比例。

從邊際上看，α是正态分布，給定西格瑪的平方，以注意不在這個無機關規模參數作為子α的方差控制。 同樣地，β是正态分布，給定西格瑪的平方，其平均值為b0，方差控制由參數S次β。

阿爾法和貝塔之間的協方差，是西格瑪平方乘以一個參數，這個參數描述了我們對阿爾法和貝塔如何共同變化的先驗信念。 如果s alpha beta是0，那麼先驗的信念是，alpha将獨立于beta，條件是sigma平方。

為了完成規範，我們使用共轭先驗sigma平方，其中1動方差有一個自由度為N的gamma分布，sigma平方不是對sigma平方的先驗估計。 這與我們在以前的視訊中使用的先驗相似。 由于共轭的關系，後驗分布将是正态伽馬分布，有簡單的規則來更新參數。 

參照先驗和後驗分布

提供一個貝葉斯分析作為起點是很有用的。 參照先驗分布是作為那些方差參數到無窮大時的正态分布的極限得到的，并且在α和β中是平坦或均勻的。 而共轭伽馬先驗的極限，先驗自由度為0，提供了我們在以前視訊中使用的西格瑪平方的參照先驗。

與之相關的β的參照後驗是一個以OLS估計值為中心的學生t分布，其尺度與OLS标準誤差相同。 自由度是n減2，就像在頻率分析中一樣。

同樣，α的邊際後驗分布也是一個學生t分布，其中心和尺度由OLS估計值給出。

利用這些參數的聯合分布，我們可以得到x的任何值的預期體脂分布，這是一個學生t分布，有n減2的自由度。 估計值由OLS估計值給出。

使用統計軟體，我們可以得到參數估計值、标準差和95%置信區間。 你應該能夠确認，當我們使用這個參照先驗時，貝葉斯估計值與頻率主義估計值相同。

參數估計

參數、拟合值或預測的置信區間都是以後驗均值加或減适當的t-quantile乘以标準差得到的。 主要差別在于對區間的解釋。 例如，根據資料，我們現在認為，腰圍每增加10厘米，體脂将增加5.8%至6.9%的可能性是95%。 當然，這個模型是一個近似值，是以要小心任何因果關系的解釋，然而，它仍然可以對預測有幫助。

預測體脂

對于預測身體脂肪，我們将使用後驗預測分布。

如果我們知道參數，一個新的觀察将隻是通過取平均值來獲得，基于我們的人口回歸方程加上相關的不确定性，描述個體偏離人口平均值和x的程度。

鑒于手頭的資料，我們的後驗預測分布是一個具有n-2個自由度的t分布。 預測新值的最佳估計是人口線的後驗平均值。 這與我們之前計算的拟合值是一樣的，但是基于預測标準差的比例參數。

預測标準差包含了x處的回歸線的後驗不确定性。從最後一個項中，我們可以看到，當我們在資料的平均值附近進行預測時，變異性将是最小的。 随着處于x範圍兩端的x的值的變異性增加。

還有一個額外的估計值sigma squared，它來自于誤差epsilon的不确定性，或者說我們期望觀測值偏離回歸線的程度。 

圖中顯示了資料，預測方程是後驗平均數，用紅色表示，估計總體平均數的95%的點狀區間用灰色虛線表示。 預測體脂的95%可信區間是外側的虛線。 大多數資料都在預測區間内，正如人們所期望的那樣。 然而，橙色圈出的點的體脂率遠遠低于模型的預期。

R語言用貝葉斯線性回歸、貝葉斯模型平均 (BMA)來預測勞工工資案例

下面，貝葉斯資訊準則（BIC）和貝葉斯模型平均法被應用于建構一個簡明的收入預測模型。

這些資料是從 935 名受訪者的随機樣本中收集的。該資料集是_計量經濟學資料集_系列的一部分 。

加載包

資料将首先使用該​

​dplyr​

​​ 包進行探索 ，并使用該​

​ggplot2​

​ 包進行可視化 。稍後，實作逐漸貝葉斯線性回歸和貝葉斯模型平均 (BMA)。

資料

資料集網頁提供了以下變量描述表：

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

探索資料

與任何新資料集一樣，一個好的起點是标準的探索性資料分析。彙總表是簡單的第一步。

1.  # 資料集中所有變量的彙總表--包括連續變量和分類變量
2.  summary(wage)

因變量（工資）的直方圖給出了合理預測應該是什麼樣子的。

1.  #工資資料的簡單柱狀圖
2.  hst(wge$wae, breks = 30)

直方圖還可用于大緻了解哪些地方不太可能出現結果。

1.  # 檢查圖表 "尾部 "的點的數量
2.  sm(wage$ge < 300)
3.  ## \[1\] 6
4.  sm(wae$wge > 2000)
5.  ## \[1\] 20      

簡單線性回歸

由于周工資（'wage'）是該分析中的因變量，我們想探索其他變量作為預測變量的關系。我們在資料中看到的工資變化的一種可能的、簡單的解釋是更聰明的人賺更多的錢。下圖顯示了每周工資和 IQ 分數之間的散點圖。

gplot(wae, es(iq, wge)) + gom\_oint() +gom\_smoth()      

IQ 分數和工資之間似乎存在輕微的正線性關系，但僅靠 IQ 并不能可靠地預測工資。盡管如此，這種關系可以通過拟合一個簡單的線性回歸來量化，它給出：

工資 i = α + β⋅iqi + ϵiwagei = α + β⋅iqi + ϵi
1.  m\_wg\_iq = lm(wge ~ iq, dta = age)
2.  coefients      

工資 i = 116.99 + 8.3 ⋅iqi + ϵiwagei = 116.99 + 8.3 ⋅iqi + ϵi

在轉向貝葉斯改進這個模型之前，請注意貝葉斯模組化假設誤差 (ϵi) 以恒定方差正态分布。通過檢查模型的殘差分布來檢查該假設。如果殘差高度非正态或偏斜，則違反假設并且任何後續推論都無效。要檢查假設，請按如下方式繪制殘差：

1.  # 用散點圖和模型誤差殘差的直方圖來檢查正态性假設
2.  glot(dta = mwag_q, es(x = .ite, y = .rd)) +
3.    gemittr() +

1.  plot(dta = m\_g\_iq, aes(x = .reid)) +
2.    histgm(bnwth = 10)

變量變換

兩個圖都顯示殘差是右偏的。是以，IQ（因為它目前存在于資料集中）不應用作貝葉斯預測模型。但是，對 僅具有正值的偏斜_因_變量使用（自然）對數變換 通常可以解決問題。下面，該模型使用轉換後的工資變量進行了重新拟合。

1.  # 用IQ的自然對數拟合th模型
2.  lm(lage ~ iq, data = wae)

1.  # 殘差sctterplot和轉換後資料的柱狀圖
2.  plt(data = m\_lag\_iq, es(x = .fited, y = .reid)) 
3.    geiter() +

1.  ggpot(dta = m_lwgeiq, as(x = .resd)) +
2.    gostgam(binwth = .1) +

殘差确實大緻呈正态分布。然而，由此産生的 IQ 系數非常小（隻有 0.0088），這是可以預料的，因為 IQ 分數提高 1 分幾乎不會對工資産生太大影響。需要進一步細化。資料集包含更多資訊。

多元線性回歸和 BIC

我們可以首先在​​回歸模型​​中包含所有潛在的解釋變量，來粗略地嘗試解釋盡可能多的工資變化。

1.  # 對資料集中的所有變量運作一個線性模型，使用'.'約定。
2.  full = lm(lwge ~ . - wage, dta = wge)

完整線性模型的上述總結表明，自變量的許多系數在統計上并不顯着（請參閱第 4 個數字列中的 p 值）。選擇模型變量的一種方法是使用貝葉斯資訊準則 (BIC)。BIC 是模型拟合的數值評估，它也會按樣本大小的比例懲罰更多的參數。這是完整線性模型的 BIC：

BIC(full)      

BIC 值越小表示拟合越好。是以，BIC 可以針對各種縮減模型進行計算，然後與完整模型 BIC 進行比較，以找到适合工資預測工作的最佳模型。當然，R 有一個功能可以系統地執行這些 BIC 調整。

1.  # 用step計算模型
2.  pIC(lwge ~ . - wge, dta = na.oi(wge))lg(lgth(na.mit(wge))))

1.  # 顯示逐漸模型的BIC
2.  BIC(se_mol)

調用 step找到産生最低 BIC 的變量組合，并提供它們的系數。很不錯。

貝葉斯模型平均（BMA）

即使BIC處于最低值，我們能有多大把握确定所得到的模型是真正的 "最佳拟合"？答案很可能取決于基礎資料的規模和穩定性。在這些不确定的時候，貝葉斯模型平均化（BMA）是有幫助的。BMA對多個模型進行平均化，獲得系數的後驗值和新資料的預測值。下面，BMA被應用于工資資料（排除NA值後）。

1.  # 不包括NA
2.  a_ona = na.omt(wae)
3.   
4.  # 運作BMA，指定BIC作為判斷結果模型的标準
5.  BMA(wge ~ . -wge, daa= ae\_o\_a,
6.                     pror = "BIC", 
7.                     moepor = ufom())
8.   
9.  # 顯示結果
10.  summary      

結果表顯示了五個最有可能的模型，以及每個系數被包含在真實模型中的機率。我們看到，出生順序和是否有兄弟姐妹是最不可能被包含的變量，而教育和智商變量則被鎖定。BMA模型的排名也可以用圖像圖來顯示，它清楚地顯示哪些變量在所有模型中，哪些變量被排除在所有模型之外，以及那些介于兩者之間的變量。

ge(b_lge, tp.oels)      

我們還可以提供模型系數的95%置信區間。下面的結果支援了關于包括或排除系數的決定。例如，在區間包含零，有大量證據支援排除該變量。

confint(ceflae)      

進行預測

構模組化型後，pediction 隻是插入資料的問題：

1.  # 用一個虛構的勞工的統計資料來預測資料的例子
2.   
3.  # 進行預測
4.  redict = pedct(e_odl, newdt = wrkr,eitr = "BMA")
5.   
6.  # 将結果轉換為元
7.  exp(wk_pedct)      

預計這名化妝從業人員的周薪為 745 元。這到底有多準确？你得問她，但我們對我們的變量選擇很有信心，并對現有的資料盡了最大努力。應用的貝葉斯技術使我們對結果有信心。

資料擷取

在下面公衆号背景回複“工資資料”，可擷取完整資料。