當(dāng)您可以借用多種模型的信息時(shí)，為什么只選擇一種模型呢？新的 BMA 套件執(zhí)行貝葉斯模型平均，以解釋分析中的模型不確定性。不確定線性回歸模型中應(yīng)包含哪些預(yù)測變量？使用 bmaregress 找出哪些預(yù)測因素很重要。執(zhí)行模型選擇、推理和預(yù)測。使用許多后估計(jì)命令來探索有影響力的模型、模型復(fù)雜性、模型擬合和預(yù)測性能、對模型和預(yù)測變量重要性假設(shè)的敏感性分析等等。

　　▋簡介

　　傳統(tǒng)上，我們選擇一個(gè)模型并根據(jù)該模型進(jìn)行推理和預(yù)測。我們的結(jié)果通常沒有考慮選擇模型的不確定性，因此可能過于樂觀。如果我們選擇的模型與真實(shí)的數(shù)據(jù)生成模型（DGM）有很大不同，它們甚至可能是不正確的。在某些應(yīng)用中，我們可能有關(guān)于 DGM 的強(qiáng)有力的理論或經(jīng)驗(yàn)證據(jù)。在其他應(yīng)用中，通常具有復(fù)雜和不穩(wěn)定的性質(zhì)，例如經(jīng)濟(jì)學(xué)、心理學(xué)和流行病學(xué)中的應(yīng)用，

　　模型平均不是僅依賴一種模型，而是根據(jù)觀察到的數(shù)據(jù)對多個(gè)合理模型的結(jié)果進(jìn)行平均。在 BMA 中，模型的“合理性”由后驗(yàn)?zāi)Ｐ透怕?(PMP) 來描述，該概率是使用基本貝葉斯原理（貝葉斯定理）確定的，并普遍應(yīng)用于所有數(shù)據(jù)分析。

　　在估計(jì)模型參數(shù)和預(yù)測新觀測值時(shí)，BMA 可用于解釋模型不確定性，以避免得出過于樂觀的結(jié)論。它在具有多個(gè)合理模型的應(yīng)用中特別有用，在這些應(yīng)用中沒有明確的理由選擇特定模型而不是其他模型。但即使最終目標(biāo)是選擇單一模型，您也會(huì)發(fā)現(xiàn) BMA 是有益的。它提供了一種原則性的方法來識(shí)別所考慮的模型類別中的重要模型和預(yù)測變量。它的框架允許您了解不同預(yù)測變量之間的相互關(guān)系，即它們一起、單獨(dú)或獨(dú)立出現(xiàn)在模型中的傾向。它可用于評(píng)估最終結(jié)果對不同模型和預(yù)測變量的重要性的各種假設(shè)的敏感性。它提供了對數(shù)分?jǐn)?shù)意義上的最佳預(yù)測。它可用于評(píng)估最終結(jié)果對不同模型和預(yù)測變量的重要性的各種假設(shè)的敏感性。它提供了對數(shù)分?jǐn)?shù)意義上的最佳預(yù)測。它可用于評(píng)估最終結(jié)果對不同模型和預(yù)測變量的重要性的各種假設(shè)的敏感性。它提供了對數(shù)分?jǐn)?shù)意義上的最佳預(yù)測。

　　▋BMA 套件

　　在回歸設(shè)置中，模型不確定性相當(dāng)于模型中應(yīng)包含哪些預(yù)測變量的不確定性。我們可以使用bmaregress來解釋線性回歸中預(yù)測變量的選擇。它要么在可行的情況下使用枚舉選項(xiàng)詳盡地探索模型空間，要么使用帶有采樣的馬爾可夫鏈蒙特卡羅 (MCMC) 模型組合 (MC3)算法選項(xiàng)。它報(bào)告訪問模型的各種摘要，包括預(yù)測變量和模型參數(shù)的后驗(yàn)分布。它允許您指定要從模型中一起包含或排除的預(yù)測變量組以及包含在所有模型中的預(yù)測變量組。它為 mprior()選項(xiàng)中的模型和 gprior() 中的 g 參數(shù)提供各種先驗(yàn)分布，g 參數(shù)控制回歸系數(shù)向零收縮。選項(xiàng)。它還支持因子變量和時(shí)間序列運(yùn)算符，并提供了多種在估計(jì)過程中使用 heredity() 選項(xiàng)處理交互的方法。

　　有許多受支持的后估計(jì)功能，其中還包括一些標(biāo)準(zhǔn)貝葉斯后估計(jì)功能。

　　▋BMA 案例解析

　　下面，我們將使用toy dataset介紹一些BMA功能。該數(shù)據(jù)集包含n=200個(gè)觀測值，p=10個(gè)正交預(yù)測因子，結(jié)果y生成為

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

　　其中，是標(biāo)準(zhǔn)正態(tài)誤差項(xiàng)。

　　. 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è)預(yù)測因子和默認(rèn)固定（基準(zhǔn)）g-prior， bmaregress 使用模型枚舉并探索所有2^10=1,024可能的模型。默認(rèn)模型先驗(yàn)是 β二項(xiàng)式，形狀參數(shù)為 1，在模型大小上是均勻的。先驗(yàn)地，我們的 BMA 模型假設(shè)回歸系數(shù)幾乎沒有收縮到零，因?yàn)?g/ (1+g) = 0.9950接近于1。

　　bmaregress 將 x2和x10確定為非常重要的預(yù)測因子 - 它們的后驗(yàn)包含概率 (PIP) 基本上為 1。它們的回歸系數(shù)的后驗(yàn)均值估計(jì)分別為 1.2（四舍五入）和 5.1，非常接近 1.2 和 1.2 的真實(shí)值5。所有其他預(yù)測變量的 PIP 都低得多，并且系數(shù)估計(jì)值接近于零。我們的 BMA 調(diào)查結(jié)果與真正的 DGM 一致。 

　　讓我們存儲(chǔ)當(dāng)前的估計(jì)結(jié)果以供以后使用。與其他貝葉斯命令一樣，我們必須先保存模擬結(jié)果，然后才能使用 估計(jì)存儲(chǔ)來存儲(chǔ)估計(jì)結(jié)果。

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

　　▋可靠區(qū)間

　　默認(rèn)情況下， bmaregress 不報(bào)告可信區(qū)間 (CrIs)，因?yàn)樗枰獙δＰ蛥?shù)的后驗(yàn)樣本進(jìn)行潛在耗時(shí)的模擬。但我們可以在估計(jì)后計(jì)算它們。

　　我們首先使用 bmacoefsample 從模型參數(shù)的后驗(yàn)分布（包括回歸系數(shù)）生成 MCMC 樣本。然后，我們使用現(xiàn)有的 bayesstats summary 命令來計(jì)算生成的 MCMC 樣本的后驗(yàn)摘要。

　　. bmacoefsample, rseed(18)

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

　　.

　　bayesstats summary

　　Posterior summary statistics                      MCMC sample size =    10,000

　　x2系數(shù)的 95% CrI為 [1.05, 1.34]， x10系數(shù)的 95% CrI為 [4.9, 5.26]，這與我們的 DGM 一致。

　　▋有影響力的模型

　　BMA 系數(shù)估計(jì)值是 1,024 個(gè)模型的平均值。研究哪些模型具有影響力非常重要。我們可以使用 bmastats models來探索 PMPs。

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

　　Variable-inclusion summary

　　Legend: x - estimated

　　毫無疑問，同時(shí)包含 x2 和 x10的模型具有最高的 PMP，約為 63%。事實(shí)上，前兩個(gè)模型對應(yīng)的累積 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

　　該命令默認(rèn)繪制前 100 個(gè)模型，但 如果您想查看更多模型， 可以指定 top() 選項(xiàng)。

　　▋重要預(yù)測因素

　　我們可以通過使用 bmagraph varmap 生成變量包含圖來直觀地探索預(yù)測變量的重要性。

　　. bmagraph varmapComputing model probabilities ...

　　x2 和 x10 均包含在 PMP 排名的所有前 100 名模型中。它們的系數(shù)在所有模型中均為正值。

　　▋模型大小分布

　　我們可以通過使用 bmastats msize 和 bmagraph msize 探索先驗(yàn)和后驗(yàn)?zāi)Ｐ痛笮》植紒硖剿?BMA 模型的復(fù)雜性。

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

. bmagraph msize

　　先驗(yàn)?zāi)Ｐ痛笮》植荚谀Ｐ痛笮∩鲜蔷鶆虻?。后?yàn)?zāi)Ｐ痛笮》植枷蜃髢A斜，眾數(shù)約為 2。先驗(yàn)平均模型大小為 5，后驗(yàn)平均模型大小為 2.48。因此，根據(jù)觀察到的數(shù)據(jù)，與我們的先驗(yàn)假設(shè)相比，BMA 傾向于平均具有大約兩個(gè)預(yù)測變量的較小模型。

　　▋系數(shù)后驗(yàn)分布

　　我們可以使用bmagraph coefdensity來繪制回歸系數(shù)的后驗(yàn)分布。

　　. bmagraph coefdensity {x2} {x3}, combine

　　這些分布是對應(yīng)于預(yù)測變量的不包含概率的零點(diǎn)質(zhì)量和以包含為條件的連續(xù)密度的混合。對于x2的系數(shù)，不包含的概率可以忽略不計(jì)，因此其后驗(yàn)分布本質(zhì)上是一個(gè)以 1.2 為中心的連續(xù)且相當(dāng)對稱的密度。對于x3的系數(shù) ，除了條件連續(xù)密度之外，在零處還有一條紅色垂直線，對應(yīng)于大約 1-.2 = 0.8 的后驗(yàn)不包含概率。

　　▋BMA 預(yù)測

　　bmapredict 產(chǎn)生各種 BMA 預(yù)測。例如，讓我們計(jì)算后驗(yàn)預(yù)測平均值。

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

　　我們還可以計(jì)算預(yù)測 CrIs 來估計(jì)我們預(yù)測的不確定性。（許多傳統(tǒng)模型選擇方法無法實(shí)現(xiàn)這一點(diǎn)。）請注意，這些預(yù)測性 CrIs 還包含模型不確定性。為了計(jì)算預(yù)測 CrIs，我們必須首先保存后驗(yàn) MCMC 模型參數(shù)樣本。

　　. 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 ...

　　我們總結(jié)預(yù)測結(jié)果：

　　. summarize y pmean cri*

　　報(bào)告的預(yù)測平均值觀察平均值以及 95% 預(yù)測 CrI 界限的下限和上限相對于觀察到的結(jié)果 y 看起來是合理的。

　　▋信息模型先驗(yàn)

　　BMA 的優(yōu)勢之一是能夠整合有關(guān)模型和預(yù)測變量的先驗(yàn)信息。這使我們能夠研究結(jié)果對模型和預(yù)測變量重要性的各種假設(shè)的敏感性。假設(shè)我們有可靠的信息，除了x2和x10之外的所有預(yù)測變量與結(jié)果相關(guān)的可能性較小。例如，我們可以為這些預(yù)測變量指定一個(gè)具有低先驗(yàn)包含概率的二項(xiàng)式先驗(yàn)?zāi)Ｐ汀?/p>

　　為了稍后比較 BMA 模型的樣本外性能，我們將樣本隨機(jī)分為兩部分，并使用第一個(gè)子樣本擬合 BMA 模型。我們還保存了這個(gè)信息更豐富的 BMA 模型的結(jié)果。

　　. 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

　　該模型先驗(yàn)的效果是所有不重要預(yù)測變量的 PIP 現(xiàn)在小于 2%。

　　請注意，當(dāng)我們選擇一個(gè)模型時(shí)，我們本質(zhì)上是在用一個(gè)非常強(qiáng)的先驗(yàn)擬合 BMA 模型，即所有選定的預(yù)測變量都必須以 1 的概率包含在內(nèi)。例如，我們可以強(qiáng)制 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 進(jìn)行預(yù)測性能

　　為了比較所考慮的 BMA 模型，我們使用第一個(gè)子樣本重新擬合默認(rèn)的 BMA 模型并存儲(chǔ)估計(jì)結(jié)果。

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

　　我們可以通過使用 bmastats lps 計(jì)算樣本外觀察的對數(shù)預(yù)測得分 (LPS) 來 比較 BMA 模型的樣本外預(yù)測性能。

　　. 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 摘要非常相似。

　　▋清理

　　我們在分析過程中生成了多個(gè)數(shù)據(jù)集。我們不再需要它們，所以我們最后刪除它們。但您可能決定保留它們，特別是當(dāng)它們對應(yīng)于可能耗時(shí)的最終分析時(shí)。

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

　　北京友萬信息科技有限公司作為 Stata 軟件在中國的授權(quán)經(jīng)銷商及合作伙伴，擁有強(qiáng)大的售后服務(wù)團(tuán)隊(duì)，聚合國內(nèi)一線Stata行業(yè)專家為客戶提供優(yōu)質(zhì)的技術(shù)支持服務(wù)。已為國內(nèi)數(shù)十所高等院校及科研院所完成了Stata軟件采購計(jì)劃，用戶覆蓋經(jīng)濟(jì)、管理、生物、歷史、法學(xué)、勞動(dòng)、人口、地理、環(huán)境、教育和心理學(xué)等各個(gè)學(xué)科門類。依托Stata 原廠及自身經(jīng)驗(yàn)豐富的技術(shù)團(tuán)隊(duì)資源，從市場策略、銷售模式和宣傳力度上全面推廣，得到了StataCorp LCC原廠的全面認(rèn)可，并與國內(nèi)重點(diǎn)高校實(shí)驗(yàn)室建立一對一的定點(diǎn)服務(wù)計(jì)劃，樹立了Stata軟件在中國用戶中的良好品牌形象。