共计 15879 个字符，预计需要花费 40 分钟才能阅读完成。
原文链接：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
工资的自然对数
这是察看钻研还是试验？
- 察看钻研
摸索数据
与任何新数据集一样，规范的探索性数据分析是一个好的开始。咱们将从工资变量开始，因为它是咱们模型中的因变量。
- 对于工资问题，下列哪种说法是谬误的？
- 7 名受访者每周支出低于 300 元
summary(wage)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 115.0 669.0 905.0 957.9 1160.0 3078.0
因为工资是咱们的因变量，咱们想探讨其余变量之间的关系作为预测。
练习：排除工资和工龄，抉择另外两个你认为能够很好预测工资的变量。应用适当的图来形象化他们与工资的关系。
受教育水平和工作小时数仿佛是工人工资的良好预测因素。
ggplot(data = wage, aes(y=wage, x=exper))+geom_point()
ggplot(data = wage, aes(y=wage, x=educ))+geom_point()
简略的线性回归
对于咱们在数据中看到的工资差别，一个可能的解释是，更聪慧的人赚更多的钱。下图显示了周工资和智商得分之间的散点图。
ggplot(data = wage, aes(x = iq, y = wage)) +
geom_point()
这个图是相当芜杂的。尽管智商分数和工资之间可能存在轻微的正线性关系，但智商充其量只是一个粗略的工资预测指标。咱们能够通过拟合一个简略的线性回归来量化这一点。
m\_wage\_iq$coefficients## (Intercept) iq
## 116.991565 8.303064## \[1\] 384.7667回忆一下，在模型下
如果应用 和参考先验 ，而后贝叶斯后验均值和标准差别离等于频数预计和标准差。
贝叶斯模型标准假如误差正态分布且方差为常数。与频率法一样，咱们通过查看模型的残差散布来测验这一假如。如果残差是高度非正态或偏态的，则违反了假如，任何随后的推断都是有效的。
- 测验 m\_wage\_iq 的残差。正态分布误差的假如无效吗?
- 不，因为模型的残差散布是右偏的。
qqnorm(m\_wage\_iq$residuals)
qqline(m\_wage\_iq$residuals)
练习：从新调整模型，这次应用 educ（教育）作为自变量。你对上一个练习的答复有变动吗？
## (Intercept) educ
## 146.95244 60.21428summary(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 十分大，这些散布实际上是近似正态的。
- 在参考先验 p（α，β，σ2）∞1/σ2 下，给出 β 的 95% 后验置信区间，即 IQ 系数。
- (0.00709, 0.01050)
# 从线性模型 m\_lwage\_iq 中提取系数值
qnorm(c(0.025, 0.975), mean = iq\_mean\_estimate, sd=iq_sd)## \[1\] 0.007103173 0.010511141练习：智商系数很小，这是意料之中的，因为智商分数进步一分很难对工资产生很高的倍增效应。使系数更易于解释的一种办法是在将智商放入模型之前将其标准化。从这个新模型来看，智商进步 1 个标准差（15 分）预计工资会减少多少百分比？
智商是用 scale 函数标准化的，智商进步 15 分会引起工资的进步
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 办法都存在于解决缺失数据的数据集上，然而它们超出了本课程的范畴。
- 从这个模型来看，谁赚得更多：已婚的黑人还是独身的非黑人？
- 已婚黑人
与繁多非黑人女子相比，所有其余平等的，已婚的黑人将取得以下乘数。
married\_black <- married\_coef\*1+black_coef\*1
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。
- 从残缺模型中打消哪个变量失去最低的 BIC?
feduc
BIC(m\_lwage\_sibs)## \[1\] 581.4031BIC(m\_lwage\_feduc)## \[1\] 580.9743BIC(m\_lwage\_meduc)## \[1\] 582.3722练习：R 有一个函数 stepAIC，它将在模型空间中向后运行，删除变量直到 BIC 不再升高。它以一个残缺的模型和一个惩办参数 k 作为输出。依据 BIC（在这种状况下 k =log（n）k=log（n））找到最佳模型。
# 对于 AIC，惩办因子是一个接触值 k。对于 step BIC，咱们将应用 stepAIC 函数并将 k 设置为 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贝叶斯模型均匀
通常，几个模型都是同样可信的，只抉择一个模型疏忽了抉择模型中蕴含的变量所波及的固有不确定性。解决这一问题的一种办法是实现贝叶斯模型均匀（Bayesian model averaging，BMA），即对多个模型进行均匀，从新数据中取得系数的后验值和预测值。咱们能够应用它来实现 BMA 或抉择模型。咱们首先将 BMA 利用于工资数据。
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 命令为咱们提供了每个变量的后验模型蕴含概率和最可能的模型。例如，模型中蕴含小时数的后验概率为 0.855。此外，后验概率为 0.0455 的最可能模型包含截距、工作工夫、智商、教育水平、年龄、婚姻状况、城市生活状况和母亲教育水平。尽管 0.0455 的后验概率听起来很小，但它比调配给它的对立先验概率大得多，因为有 216 个可能的模型。
在模型平均法下，还能够可视化系数的后验散布。咱们将智商系数的后验散布绘制如下。
plot(coef_lwage, subset = c(3,13), ask=FALSE)
咱们还能够为这些系数提供 95% 的置信区间：
## 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，咱们将应用简化的数据集，其中不包含兄弟姐妹数量、出世程序和父母教育。
wage_red = wage %>% dplyr::select(-sibs, -brthord, -meduc, -feduc)- 基于这个简化的数据集，依据贝叶斯模型均匀，下列哪一个变量的边际后验蕴含概率最低？
- 年龄
##
## 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- 判断: 蕴含所有变量的奢侈模型的后验概率大于 0.5。(系数应用 Zellner-Siow 零先验，模型应用 β 二项 (1,1) 先验)
- 真的
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练习：用数据集绘制年龄系数的后验分布图。
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 次。
set.seed(314)
N = 100000
y_5 = rnorm(N, mu, sqrt(sigma2))咱们能够通过观察模仿数据直方图的平滑版本，查看预测散布的估计值。
新观测的 95% 核心置信区间为在这种状况下，L 是 0.025 分位数，U 是 0.975 分位数。咱们能够应用分位数函数来取得这些值，从而找到 tracy5 的 0.025 和 0.975 的样本分位数。
- 预计一个新的观测 y595% 置信区间
- (-3.11, 5.13)
quantile(y_5, probs = c(0.025, 0.975))## 2.5% 97.5%
## -3.109585 5.132511练习：在下面的简略例子中，能够应用积分来剖析计算后验预测值。在这种状况下，它是一个具备 3 个自由度（n−1）的 t 散布。绘制 y 的教训密度和 t 散布的理论密度。它们之间有什么比拟？
plot(den\_of\_y)
BAS 预测
在 BAS 中，用贝叶斯模型平均法结构预测区间是通过仿真实现的，而在模型抉择的状况下，用预测区间进行准确推理往往是可行的。
回到工资数据集，让咱们找到最佳预测模型下的预测值，即预测值最靠近 BMA 和相应的后验标准差的模型。
predict(bma_lwage, estimator="BPM")## \[1\] "Intercept" "hours" "iq" "kww" "educ"
## \[6\] "exper" "tenure" "age" "married1" "urban1"
## \[11\] "meduc"咱们能够将其与咱们之前发现的最高概率模型和中位概率模型（MPM）进行比拟
predict(bma_lwage, estimator="MPM")## \[1\] "Intercept" "hours" "iq" "educ" "exper"
## \[6\] "tenure" "age" "married1" "urban1" "meduc"MPM 除了蕴含 HPM 的所有变量外，还蕴含 exper，而 BPM 除了 MPM 中的所有变量外，还蕴含 kwh。
练习：应用简化数据，最佳预测模型、中位概率模型和最高后验概率模型中蕴含哪些协变量？
让咱们来看看 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"能够失去预测对数工资的 95% 可信区间
## 2.5% 97.5% pred
## 6.661869 8.056452 7.359160换算成工资，咱们能够将区间指数化。
## 2.5% 97.5% pred
## 782.0108 3154.0780 1570.5169取得一个 95% 的预测区间的工资。
如果应用 BMA，区间是
## 2.5% 97.5% pred
## 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 模型实现
