ベイズモデル平均化(BMA)

 従来、自分でモデルを選択し、それに基づいて分析を行っていました。この場合得られた分析結果は、選択したモデルに依拠します。妥当性のあるモデルが複数存在する場合、 この分析アプローチは信頼できるものではないかもしれません。モデル平均化では、複数モデルを元にした分析が可能となり、結果におけるモデル設定の不確かさを担保します。 BMAはベイズ理論によって、あらゆるデータ分析方法における、モデルの不確かさを担保します。回帰分析においては、モデルの不確実性は、どの予測変数を回帰モデルに 含める必要があるかについての不確実性を表します。
 新しいコマンドbmaregressは、線形回帰のBMAを実行し、推論、予測、および必要に応じてモデルの選択にも使用できます。次のように行います。

. bmaregress y x1 x2

 これはアウトカムyに対して予測変数x1, x2を含む・除外する、妥当なモデルを、観測データに基づき各モデルがどの程度尤もらしいによって、統合します。 さまざまな事前分布から選択して、結果に対するモデルと予測変数の重要性に関する仮定の影響を調べることができます。
 推定後分析コマンドではモデルの確率、予測変数の重要度、モデルの複雑さの確認、予測平均の取得、予測の評価、回帰係数の推論が可能です。

bmaコマンド群

 回帰分析では、モデルの不確実性は、どの予測子をモデルに含める必要があるかという不確実性に相当します。bmaregressを使用して、線形回帰での予測子の選択を説明できます。 可能な場合は常にenumerationオプションを使用してモデル空間を徹底的に調査するか、samplingオプションを使用して マルコフ連鎖モンテカルロ(MCMC)モデル合成(MC3)アルゴリズムを使用して、モデル空間を探査し、考慮したモデルと含まれる予測変数、およびモデルパラメータの事後分布について様々な概要を報告します。 モデルに含めたり除外したりする予測変数のグループと、すべてのモデルに含まれる予測変数を指定できます。mprior()オプションでモデルのさまざまな事前分布を提供し、 gprior()オプションでゼロに向かう回帰係数の縮小を制御するgパラメータを提供します。また、因子変数と時系列演算子もサポートし、heredity()オプションを使用して 推定中に相互作用を処理する方法をいくつか提供します。
 サポートされている多くの事後推定機能があり、標準のベイジアン事後推定機能もいくつか含まれています。

操作例

 bmaコマンドの操作例を紹介します。使用するデータセットは観測行(n)が200行、列数(p)が10列です。アウトカム y を次式によって生成します。

y = 0.5 + 1.2 x2 + 5 x10 + ε

 誤差項 ε は標準正規分布に従います。

. webuse bmaintro . summarize

モデルの列挙

 bmaregressコマンドでx1からx10までをyに回帰します。

. bmaregress y x1-x10

 10の予測変数とデフォルトのg事前分布では、bmaregressは想定し得る全モデル、210=1,024個を探索します。 デフォルトのモデル事前分布は形状パラメータ1のベータ・二項分布で、モデルサイズに対して一様です。BMAモデルは g/(1+g)=0.9950 が1に近いので、 回帰係数はゼロ方向に若干収縮することを前提としています。
 bmaregressは、x2とx10を非常に重要な予測変数として特定しました。それらの事後包含確率(PIP)は基本的に1です。 それらの回帰係数の事後平均推定値は、それぞれ1.2と5.1であり、真の値である1.2と5です。他のすべての予測因子は、はるかに低いPIPとゼロに近い係数推定値を持っています。 私たちのBMAの調査結果は、真のデータ生成過程と一致しています。
 後で使用するために、現在の推定結果を保存しましょう。 他のベイズコマンドと同様に、estimates storeを使用して推定結果を保存する前に、 最初にシミュレーション結果を保存する必要があります。

. bmaregress, saving(bmareg)
. estimates store bmareg

信用区間

 bmaregressは、モデルパラメータの事後標本のシミュレーションに時間がかかることがあるため、既定では信用区間(CrIs)を報告しません。しかし、推定後にそれらを計算できます。
 最初にbmacoefsampleを使用して、回帰係数を含むモデルパラメータの事後分布からMCMCサンプルを生成します。次に、既存のbayesstats summaryコマンドを使用して、 生成されたMCMCサンプルの事後要約を計算します。

. bmacoefsample, rseed(18)
. bayesstats summary

 x2、x10の95%信用区間はそれぞれ[1.05, 1.34], [4.9, 5.26]でデータ生成過程と一致します。

影響力のあるモデル

 BMA係数の推定値は、1,024モデルにわたって平均化されています。どのモデルが影響力があるかを調査することが重要です。bmastatsモデルを使用してPMPを探索できます。

. bmastats models

 当然のことながら、x2とx10の両方を含むモデルのPMPは約63%と最高です。実際、上位2つのモデルは累積PMPが約77%に相当します。

. bmastats models, cumulative(.75)

bmagrahp pmpで累積PMPをプロットできます。

. bmagraph pmp, cumulative

 このコマンドは、デフォルトで最初の100モデルをプロットしますが、さらに多くのモデルを表示したい場合は、top()オプションを指定できます。

重要な予測変数

 bmagraph varmapを使用して変数包含マップを作成することにより、予測子の重要性を視覚的に調べることができます。

. bmagraph varmap

 x2とx10は、PMP でランク付けされた上位100モデルのすべてに含まれています。それらの係数の符号は、すべてのモデルで正です。

モデルサイズの分布

 bmastats msizeとbmagraph msizeを使用して事前および事後モデルサイズの分布を調べることで、BMAモデルの複雑さを調べることができます。

. bmastats msize
. bmagraph msize

 事前情報におけるモデルサイズの分布は、全体で均一です。事後モデルサイズ分布は、約2の最頻値で左に歪んでいます。事前分布における平均モデルサイズは5で、 事後モデルサイズは2.48です。したがって、観察されたデータに基づいて、BMAは、アプリオリな仮定と比較して、平均して約2つの予測因子を持つより小さなモデルを好む傾向があります。

係数の事後分布

 bmagraph coefdensity を使用して、回帰係数の事後分布をプロットできます。

. bmagraph coefdensity {x2} {x3}, combine

 これらの分布は、予測子の非含有確率に対応するゼロ点の質量と、含有を条件とする連続密度の混合物です。x2の係数の場合、非包含の確率は無視できるため、 その事後分布は基本的に1.2を中心とする連続したかなり対称的な密度になります。x3の係数については、条件付き連続密度に加えて、おおよそ 1 - 0.2 = 0.8 の事後非包含確率に対応するゼロに赤い垂直線があります。

BMAの予測

 bmapredictは、さまざまなBMA予測を生成します。たとえば、事後予測平均を計算してみましょう。

. bmapredict pmean, mean

 また、予測の不確実性を推定するために、予測CrIを計算することもできます。(これは、多くの従来のモデル選択方法では利用できません。)これらの予測CrIには、 モデルの不確実性も組み込まれていることに注意してください。予測CrIを計算するには、最初に事後MCMCモデルパラメータサンプルを保存する必要があります。

. bmacoefsample, saving(bmacoef)
. bmapredict cri_l cri_u, cri rseed(18)

 予測結果を要約します。

. summarize y pmean cri*

 予測平均と95%予測CrIの下限と上限の観察に対する報告された平均は、観察された結果yに比べて妥当に見えます。

モデルの事前分布

 BMAの強みの1つは、モデルと予測因子に関する事前情報を組み込むことができることです。 これにより、モデルと予測変数の重要性に関するさまざまな仮定に対する結果の感度を調べることができます。x2 と x10 を除くすべての予測因子が結果に関連する可能性が 低いという信頼できる情報があるとします。 たとえば、これらの予測子の事前包含確率が低い二項モデル事前指定を指定できます。
 後でBMAモデルのサンプル外のパフォーマンスを比較するために、サンプルをランダムに 2つの部分に分割し、最初のサブサンプルを使用して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)
. estimates store bmareg_inf

 この前のモデルの効果は、すべての重要でない予測変数のPIPが2%未満になったことです。
1つのモデルを選択する場合、基本的に、選択したすべての予測子を確率 1 で含める必要があるという非常に強い事前確率で BMAモデルを適合させることに注意してください。 たとえば、モデルにすべての変数を含めるようにbmaregressを強制できます。

. bmaregress y (x1-x10, always) if sample == 1, saving(bmareg_all)

LPSによる予測パフォーマンス

 考慮されたBMAモデルを比較するために、最初のサブサンプルを使用してデフォルトのBMAモデルを再調整し、推定結果を保存します。

. quietly 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, compact

 より有益なbmareg_infモデルの平均LPSはわずかに小さくなっていますが、すべてのモデルのLPSの要約統計量は非常に類似しています。 相互検証を使用してBMAモデルのパフォーマンスを比較する方法については、英文PDFマニュアルの[BMA] bmastats lps を参照してください。

Stata is a registered trademark of StataCorp LLC, College Station, TX, USA, and the Stata logo is used with the permission of StataCorp.

page_top_icon