Stata分析機能例題を利用してもっと知ろう
ベイズ分析

サンプルデータをダウンロードし、データを確認します。

データの確認

　　　　　

Kuehl (2000, 551)の酸素の取り込みと運動の研究例を考えます。
研究の目的は、2種類の運動プログラム：12週間のステップエアロビクスまたは12週間の平地ランニングを行った際の酸素取り込み量の変化を分析する事。
例題データセット、oxygen.dtaをダウンロードしましょう。
例題データセットのサンプルサイズは12件です。

データセットには次の変数があります。

• change ： 酸素取り込み量の変化(l/m)
• group ： 運動グループ(0: ランニング、1:エアロビクス)
• age ： 年齢
• ageXgr ： 年齢とグループの交差項

Kuehl (2000)では、このデータを分析するために共分散分析を使用しますが、ここでは以下の線形回帰モデルを考えます。

例題1 ： OLSでのフィット

　　　　　

まずは、最小二乗法モデルにフィットした場合を確認しましょう！

サンプルサイズが12件での結果をそのまま信頼していいのでしょうか？

帰無仮説H_o: β_group=0は、p値が0.014なので、有意水準5％で棄却され、groupとageはどちらもアウトカムに対して有意であることがわかります。
しかし、p値は、使用したデータにおいてどの程度、帰無仮説が正しいのかという疑問に答えるものであり、帰無仮説そのものの確率ではありません。

groupの95％信頼区間は[1.38 9.51]であり、間に0を含みません、このためgroupはchangeの有意な予測因子であると言えます。
ただし、groupの真の係数が95％の確率で、1.38と9.51の間に存在すると結論付けることは出来ません。
パラメータの確率的な区間の解釈は、ベイズ統計でのみ可能なものです。

例題2 ： 事前情報のないベイズ正規線形回帰

事前情報のあるベイズ線形回帰で、正確性を高めてみましょう。

先程の例題1の頻度統計では、パラメータは未知の固定量であるとして、パラメータの確率的要約を得ることはできないと述べました。
ベイズ統計はパラメータの事後分布を様々な観点から推定することに焦点を置き観測されたデータの持つ情報で事前分布を更新します。
つまり事後分布はパラメータの事前分布とパラメータを与えるデータの尤度関数で表されます。

bayesmhコマンドで、oxygen.dtaをベイズ線形回帰でフィットしましょう。

ここのベイズ線形モデルでは4つのパラメータ（3つの回帰係数とデータの分散）を利用します。
アウトカム、changeは正規分布に従うと仮定し、分散はJeffreysの無情報事前分布で開始します。
ジェフリーズの事前分布では、係数の同時事前分布と分散が分散の逆数に従います。
Xはデザイン行列、β=(β_0,β_group,β_age)'は係数のベクトルです。

コマンドの後に従属変数名、さらに共変量と続きます。
アウトカムの分布はlikelihood()オプションで指定し、事前分布はprior()オプションで指定します。

分散パラメータを指定しなければなりませんので、{var}として定義しています。
3つの回帰係数{change:group}、{change:age}、{change:_cons}はbayesmhによって自動的に作成されます。

分布をlikelihood()オプション内でnormal({var})と設定し、尤度関数の分散パラメータを{var}としています。
この設定と回帰の設定を合わせて、アウトカムchangeの尤度モデルが定義されます。
事前分布にはprior({change:}, flat)で全ての回帰係数に密度1の一様分布を適用します。
{change:}はモデル名changeの全てのパラメータの省略表現です。
最後に事前分布jeffreysで分散パラメータ{var}に密度1/σ^2を指定します。

乱数シードを設定して、分析結果を再現可能にします。
以降の例題ではすべて乱数シード値に14を使用します。

bayesmhはモデルの要約のヘッダの右側では、MCMCの反復回数が12,500であり、これはMCMCサンプルから除外されるデフォルト2,500回のバーンインとデフォルトが10,000回のMCMCサンプルまたはMCMCサンプルサイズの合計です。

受容率は推定されたパラメータ値をアルゴリズムが採択する割合です。
この例の受容率0.14は10,000のパラメータの14％がアルゴリズムで採択されることを示します。
この数値は通常30％に満たないことが大半です。受容率が低い(10％未満)場合は、収束に問題が生じていることがあります。

最後に、bayesmhは推定結果の要約表を報告します。
Meanは事後平均の推定結果、パラメータの周辺事後分布の平均を報告しています。事後平均推定は例題1のOLS推定の結果に非常に近くなります。
この結果は、無情報事前分布を使用しているためです。

Std. Dev.は、推定された事後標準誤差を報告します、これは事後分布の標準誤差です。
この値はパラメータの事後分布の変動性を示し、OLSの標準誤差に相当します。

Medianは事後分布の中央値を報告し、分布の対称性を判断することができます。
事後平均と中央値は回帰係数に近い値を取っていますので、事後分布が対称であることが考えられます。

最後の2列は、パラメータの信用区間を表します。信頼区間とは異なり、これらの区間は直感的に確率で解釈できます。
例えば、groupの係数が1.16と9.25の間にある確率はおよそ0.95です。
区間の下限は0より大きいので、運動プログラムは酸素取り込みの変化に影響を与えていると結論付けることができます。

例題3 ： 事前情報のあるベイズ線形回帰

単純化するためにすべての係数が独立で平均0、分散σ^2の正規分布に従い、分散パラメータは上記のように逆ガンマ分布に従うものとします。

このモデルをbayesmhコマンドでフィットします。
normal(0, {var})で係数の事前分布を、igamma(2.5, 2.5)で分散の事前分布を指定します。

例題4 ： 多変量事前情報のベイズ正規線形回帰

正規線形回帰の回帰係数でよく利用されるZellnerのg事前分布(Zellner 1986)を使用します。
gは事前標本サイズ、ν_0逆ガンマ分布の事前自由度、σ_0^2は逆ガンマ分布の事前分散です。
この事前情報は係数間で依存関係があります。ここでは、Hoff (2009)の値に近いパラメータ：g=12,ν_0=1,σ_0^2=8を使用します。

bayesmhコマンドはzellnersg0()を使用して上記の事前情報を設定します。
はじめに、分散の次元を記述します。この例では3とします。
次に事前分布の自由度として12を使用します。
最後に分散のパラメータを記述します。この例では{var}です。

例題5 ： 収束の確認

bayesgraphコマンドでパラメータ推定におけるMCMCの収束を視覚的に確認できます。
診断で良く用いられるグラフは、bayesgraph diagnosticsでまとめて表示できます。

表示される診断結果には、トレースプロット、自己相関プロット、ヒストグラム、MCMCサンプルの全体、前半、後半の密度とカーネル密度推定のオーバーレイです。
トレースプロットと自己相関プロットは高い自己相関を示しています。
ヒストグラムは単峰分布ではありません。
これらから、収束に問題があることがわかります。

有効なサンプルサイズと収束に関する統計量をbayesstats essコマンドで表示します。

MCMCサンプルの自己相関が低いほど、ESS推定量はMCMCサンプルサイズに近くなり、推定したパラメータはより正確なものになります。
このような場合、分散パラメータを独立したブロックに分割して、回帰係数とは別にパラメータのアップデート行うことができます。
bayesmhコマンドでは、block()オプションで指定します。

estimates store agegroupを使用して、推定結果をagegroupに保存します。
bayesmhの後にestimates storeを使用するには、saving()オプションでbayesmhのシミュレーション結果をStataのデータセットとして保存する必要があります。

例題6 ： 推定後要約

モデルパラメータとパラメータの関数について、bayesstats summaryコマンドで推定後に要約して表示することができます。

例題7 ： ベイズ予測

ベイズ予測は、モデルの適合度を確認し、観測値の将来の値を予測するのに適しています。
bayespredictを使って、アウトカム変数changeの反復標本を作成し、新しいデータセットchange_pred.dtaとして保存します。

bayesstats summaryコマンドで予測された観測値の事後要約を計算させます。
推定を行います。

1列目は、事後平均、事後予測分布による期待されるアウトカムのMCMC推定です。

複製標本と観測標本を比較します。これら2つの差異は、bayesstats ppvaluesコマンドで事後予測p値として計測されます。

推定された事後予測p値はすべて0.05と0.95の間(_ysim1_2を除いて)に収まっています。この結果から各観測値が適切にフィットされていることがわかります。

新しく2つの観測値を追加し、bayespredictで標本外予測を行います。

例題8 ： モデルの比較

ここまで推定モデルと誤差項を含む完全モデルを比較してみましょう。
Zellnerのg事前分布と分散の逆ガンマ事前分布をもう一度使用します。
事前パラメータの値は例題4のものを使用します。

モデルを比較するには、bayesstats icコマンドを使用します。

例題9 ： 仮説検定

例題8を使用して、モデルの実際の確率を計算します。これには、bayestest modelコマンドを使用します。

どちらのモデルの事前確率が等しいという仮定の下では、交差項を含まないモデル、agegroupの確率は、0.78で、fullモデルでは0.22です。
fullモデルより優れていることをより強く裏付けるには、agegroupの確率がより大きく(0.9以上)なければなりません。

パラメータ区間の仮設検定を行うこともできます。
交差項の無いモデルにおいて、groupの係数が4と8の間にある確率を計算しましょう。

例題10 ： シミュレーションデータセットの削除

目的の分析ができたら、bayesmhで作成したシミュレーション用のデータセットは不要なので削除するようにしましょう。