【Stata18 新功能】线性回归的贝叶斯模型平均 (BMA)

　　当您可以借用多种模型的信息时，为什么只选择一种模型呢？新的 BMA 套件执行贝叶斯模型平均，以解释分析中的模型不确定性。不确定线性回归模型中应包含哪些预测变量？使用 bmaregress 找出哪些预测因素很重要。执行模型选择、推理和预测。使用许多后估计命令来探索有影响力的模型、模型复杂性、模型拟合和预测性能、对模型和预测变量重要性假设的敏感性分析等等。

　　▋简介

　　传统上，我们选择一个模型并根据该模型进行推理和预测。我们的结果通常没有考虑选择模型的不确定性，因此可能过于乐观。如果我们选择的模型与真实的数据生成模型（DGM）有很大不同，它们甚至可能是不正确的。在某些应用中，我们可能有关于 DGM 的强有力的理论或经验证据。在其他应用中，通常具有复杂和不稳定的性质，例如经济学、心理学和流行病学中的应用，

　　模型平均不是仅依赖一种模型，而是根据观察到的数据对多个合理模型的结果进行平均。在 BMA 中，模型的“合理性”由后验模型概率 (PMP) 来描述，该概率是使用基本贝叶斯原理（贝叶斯定理）确定的，并普遍应用于所有数据分析。

　　在估计模型参数和预测新观测值时，BMA 可用于解释模型不确定性，以避免得出过于乐观的结论。它在具有多个合理模型的应用中特别有用，在这些应用中没有明确的理由选择特定模型而不是其他模型。但即使最终目标是选择单一模型，您也会发现 BMA 是有益的。它提供了一种原则性的方法来识别所考虑的模型类别中的重要模型和预测变量。它的框架允许您了解不同预测变量之间的相互关系，即它们一起、单独或独立出现在模型中的倾向。它可用于评估最终结果对不同模型和预测变量的重要性的各种假设的敏感性。它提供了对数分数意义上的最佳预测。它可用于评估最终结果对不同模型和预测变量的重要性的各种假设的敏感性。它提供了对数分数意义上的最佳预测。它可用于评估最终结果对不同模型和预测变量的重要性的各种假设的敏感性。它提供了对数分数意义上的最佳预测。

　　▋BMA 套件

　　在回归设置中，模型不确定性相当于模型中应包含哪些预测变量的不确定性。我们可以使用bmaregress来解释线性回归中预测变量的选择。它要么在可行的情况下使用枚举选项详尽地探索模型空间，要么使用带有采样的马尔可夫链蒙特卡罗 (MCMC) 模型组合 (MC3)算法选项。它报告访问模型的各种摘要，包括预测变量和模型参数的后验分布。它允许您指定要从模型中一起包含或排除的预测变量组以及包含在所有模型中的预测变量组。它为 mprior()选项中的模型和 gprior() 中的 g 参数提供各种先验分布，g 参数控制回归系数向零收缩。选项。它还支持因子变量和时间序列运算符，并提供了多种在估计过程中使用 heredity() 选项处理交互的方法。

　　有许多受支持的后估计功能，其中还包括一些标准贝叶斯后估计功能。

　　▋BMA 案例解析

　　下面，我们将使用toy dataset介绍一些BMA功能。该数据集包含n=200个观测值，p=10个正交预测因子，结果y生成为

　　y=0.5+1.2×x2+5×x10+ϵ

　　其中，是标准正态误差项。

　　. webuse bmaintro(Simulated data for BMA example). summarize

　　▋模型枚举

　　我们使用 bmaregress来拟合 y在 x1 到 x10 上的 BMA 线性回归。

. bmaregress y x1-x10Enumerating models ...Computing model probabilities ...Bayesian model averaging                            No. of obs         =    200Linear regression                                   No. of predictors  =     10Model enumeration                                               Groups =     10Always =      0Priors:                                             No. of models      =  1,024Models: Beta-binomial(1, 1)                           For CPMP >= .9 =      9Cons.: Noninformative                            Mean model size    =  2.479Coef.: Zellner's gg: Benchmark, g = 200                        Shrinkage, g/(1+g) = 0.9950sigma2: Noninformative                            Mean sigma2        =  1.272y    Mean   Std. dev.                          Group        PIP    x2    1.198105   .0733478                               2          1    x10    5.08343   .0900953                              10          1    x3    -.0352493   .0773309                               3     .21123    x9    .004321   .0265725                               9    .051516    x1    .0033937   .0232163                               1    .046909    x4    -.0020407   .0188504                               4    .039267    x5    .0005972   .0152443                               5    .033015    x9    -.0005639   .0153214                               8    .032742    x7    -8.23e-06    .015497                               7    .032386    x6    -.0003648   .0143983                               6    .032361    Always    _cons    .5907923   .0804774                               0          1    Note: Coefficient posterior means and std. dev. estimated from 1,024 models.Note: Default priors are used for models and parameter g.

　　通过 10 个预测因子和默认固定（基准）g-prior， bmaregress 使用模型枚举并探索所有2^10=1,024可能的模型。默认模型先验是 β二项式，形状参数为 1，在模型大小上是均匀的。先验地，我们的 BMA 模型假设回归系数几乎没有收缩到零，因为 g/ (1+g) = 0.9950接近于1。

　　bmaregress 将 x2和x10确定为非常重要的预测因子 - 它们的后验包含概率 (PIP) 基本上为 1。它们的回归系数的后验均值估计分别为 1.2（四舍五入）和 5.1，非常接近 1.2 和 1.2 的真实值5。所有其他预测变量的 PIP 都低得多，并且系数估计值接近于零。我们的 BMA 调查结果与真正的 DGM 一致。 

　　让我们存储当前的估计结果以供以后使用。与其他贝叶斯命令一样，我们必须先保存模拟结果，然后才能使用 估计存储来存储估计结果。

. bmaregress, saving(bmareg)note: file bmareg.dta saved.. estimates store bmareg

　　▋可靠区间

　　默认情况下， bmaregress 不报告可信区间 (CrIs)，因为它需要对模型参数的后验样本进行潜在耗时的模拟。但我们可以在估计后计算它们。

　　我们首先使用 bmacoefsample 从模型参数的后验分布（包括回归系数）生成 MCMC 样本。然后，我们使用现有的 bayesstats summary 命令来计算生成的 MCMC 样本的后验摘要。

　　. bmacoefsample, rseed(18)

　　Simulation (10000): ....5000....10000 done

　　.

　　bayesstats summary

　　Posterior summary statistics                      MCMC sample size =    10,000

　　x2系数的 95% CrI为 [1.05, 1.34]， x10系数的 95% CrI为 [4.9, 5.26]，这与我们的 DGM 一致。

　　▋有影响力的模型

　　BMA 系数估计值是 1,024 个模型的平均值。研究哪些模型具有影响力非常重要。我们可以使用 bmastats models来探索 PMPs。

　　. bmastats modelsModel summary         Number of models:Visited = 1,024Reported =     5

　　Variable-inclusion summary

　　Legend: x - estimated

　　毫无疑问，同时包含 x2 和 x10的模型具有最高的 PMP，约为 63%。事实上，前两个模型对应的累积 PMP 约为 77%：

　　. bmastats models, cumulative(0.75)Computing model probabilities ...Model summary           Number of models:Visited = 1,024Reported =     2

　　Variable-inclusion summary

　　Legend:

　　x - estimated

　　我们可以使用 bmagraph pmp 来绘制累积 PMP。

　　. bmagraph pmp, cumulative

　　该命令默认绘制前 100 个模型，但 如果您想查看更多模型， 可以指定 top() 选项。

　　▋重要预测因素

　　我们可以通过使用 bmagraph varmap 生成变量包含图来直观地探索预测变量的重要性。

　　. bmagraph varmapComputing model probabilities ...

　　x2 和 x10 均包含在 PMP 排名的所有前 100 名模型中。它们的系数在所有模型中均为正值。

　　▋模型大小分布

　　我们可以通过使用 bmastats msize 和 bmagraph msize 探索先验和后验模型大小分布来探索 BMA 模型的复杂性。

. Model-size summaryNumber of models = 1,024Model size:Minimum =  0Maximum = 10

. bmagraph msize

　　先验模型大小分布在模型大小上是均匀的。后验模型大小分布向左倾斜，众数约为 2。先验平均模型大小为 5，后验平均模型大小为 2.48。因此，根据观察到的数据，与我们的先验假设相比，BMA 倾向于平均具有大约两个预测变量的较小模型。

　　▋系数后验分布

　　我们可以使用bmagraph coefdensity来绘制回归系数的后验分布。

　　. bmagraph coefdensity {x2} {x3}, combine

　　这些分布是对应于预测变量的不包含概率的零点质量和以包含为条件的连续密度的混合。对于x2的系数，不包含的概率可以忽略不计，因此其后验分布本质上是一个以 1.2 为中心的连续且相当对称的密度。对于x3的系数 ，除了条件连续密度之外，在零处还有一条红色垂直线，对应于大约 1-.2 = 0.8 的后验不包含概率。

　　▋BMA 预测

　　bmapredict 产生各种 BMA 预测。例如，让我们计算后验预测平均值。

　　. bmapredict pmean, meannote: computing analytical posterior predictive means.

　　我们还可以计算预测 CrIs 来估计我们预测的不确定性。（许多传统模型选择方法无法实现这一点。）请注意，这些预测性 CrIs 还包含模型不确定性。为了计算预测 CrIs，我们必须首先保存后验 MCMC 模型参数样本。

　　. bmacoefsample, saving(bmacoef)note: saving existing MCMC simulation results without resampling; specifyoption simulateto force resampling in this case.note: file bmacoef.dtasaved.. bmapredict cri_l cri_u, cri rseed(18)note: computing credible intervals using simulation.Computing predictions ...

　　我们总结预测结果：

　　. summarize y pmean cri*

　　报告的预测平均值观察平均值以及 95% 预测 CrI 界限的下限和上限相对于观察到的结果 y 看起来是合理的。

　　▋信息模型先验

　　BMA 的优势之一是能够整合有关模型和预测变量的先验信息。这使我们能够研究结果对模型和预测变量重要性的各种假设的敏感性。假设我们有可靠的信息，除了x2和x10之外的所有预测变量与结果相关的可能性较小。例如，我们可以为这些预测变量指定一个具有低先验包含概率的二项式先验模型。

　　为了稍后比较 BMA 模型的样本外性能，我们将样本随机分为两部分，并使用第一个子样本拟合 BMA 模型。我们还保存了这个信息更丰富的 BMA 模型的结果。

　　. splitsample, generate(sample) nsplit(2) rseed(18). bmaregress y x1-x10 if sample == 1, mprior(binomial x2 x10 0.5 x1 x3-x9 0.05)

　　saving(bmareg_inf)

　　Enumerating models ...Computing model probabilities ...Bayesian model averaging                          No. of obs         =    100Linear regression                                 No. of predictors  =     10Model enumeration                                             Groups =     10Always =      0Priors:                                           No. of models      =  1,024Models: Binomial, IP varies                         For CPMP >= .9 =      1Cons.: Noninformative                          Mean model size    =  2.072Coef.: Zellner's gg: Benchmark, g = 100                      Shrinkage, g/(1+g) = 0.9901sigma2: Noninformative                          Mean sigma2        =  1.268

　　Note: Coefficient posterior means and std. dev. estimated from 1,024 models.Note: Default prior is used for parameter g.Note: 4 predictors with PIP less than .01 not shown.file bmareg_inf.dtasaved.. estimates store bmareg_inf

　　该模型先验的效果是所有不重要预测变量的 PIP 现在小于 2%。

　　请注意，当我们选择一个模型时，我们本质上是在用一个非常强的先验拟合 BMA 模型，即所有选定的预测变量都必须以 1 的概率包含在内。例如，我们可以强制 bmaregress 包含模型中的所有 变量：

　　. bmaregress y (x1-x10, always) if sample == 1, saving(bmareg_all)Enumerating models ...Computing model probabilities ...Bayesian model averaging                           No. of obs         =    100Linear regression                                  No. of predictors  =     10Model enumeration                                              Groups =      0Always =     10Priors:                                            No. of models      =      1Models: Beta-binomial(1, 1)                          For CPMP >= .9 =      1Cons.: Noninformative                           Mean model size    = 10.000Coef.: Zellner's gg: Benchmark, g = 100                       Shrinkage, g/(1+g) = 0.9901sigma2: Noninformative                           Mean sigma2        =  1.192

　　Note: Coefficient posterior means and std. dev. estimated from 1 model.Note: Default priors are used for models and parameter g.file bmareg_all.dtasaved.. estimates store bmareg_all

　　▋使用 LPS 进行预测性能

　　为了比较所考虑的 BMA 模型，我们使用第一个子样本重新拟合默认的 BMA 模型并存储估计结果。

　　. qui bmaregress y x1-x10 if sample == 1, saving(bmareg, replace). estimates store bmareg

　　我们可以通过使用 bmastats lps 计算样本外观察的对数预测得分 (LPS) 来 比较 BMA 模型的样本外预测性能。

　　. bmastats lps bmareg bmareg_inf bmareg_all if sample == 2, compactLog predictive-score (LPS)Number of observations = 100

　　Notes: Using analytical PMPs.Result bmareg_infhas the smallest mean LPS.

　　信息更丰富的 bmareg_inf 模型的平均 LPS 稍小，但所有模型的 LPS 摘要非常相似。

　　▋清理

　　我们在分析过程中生成了多个数据集。我们不再需要它们，所以我们最后删除它们。但您可能决定保留它们，特别是当它们对应于可能耗时的最终分析时。

　　. erase bmareg.dta. erase bmacoef.dta. erase bmareg_inf.dta. erase bmareg_all.dta

　　北京友万信息科技有限公司作为 Stata 软件在中国的授权经销商及合作伙伴，拥有强大的售后服务团队，聚合国内一线Stata行业专家为客户提供优质的技术支持服务。已为国内数十所高等院校及科研院所完成了Stata软件采购计划，用户覆盖经济、管理、生物、历史、法学、劳动、人口、地理、环境、教育和心理学等各个学科门类。依托Stata 原厂及自身经验丰富的技术团队资源，从市场策略、销售模式和宣传力度上全面推广，得到了StataCorp LCC原厂的全面认可，并与国内重点高校实验室建立一对一的定点服务计划，树立了Stata软件在中国用户中的良好品牌形象。