【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軟件在中國用戶中的良好品牌形象。