➿
二変量プロビットモデルのベイズ推定
二値アウトカムが2つあってその間の相関を考えたモデルを使いたい場合、二変量プロビットモデルが一つの選択肢になります。この記事では、説明変数がなく切片項だけの場合の二変量プロビットモデルをStanで実装する例を紹介します。なお、片方のアウトカムが欠測しているデータも含めてあてはめることができます。
実装するにあたって、二変量正規分布の累積分布関数を計算する必要があります。次のBenjamin Goodrich氏のGitHubの投稿と清水先生のブログ記事の実装を参考にしました。
ここではデータ例として、次のRコードの4～24行目でデータ生成しました。
- 6行目：
N_omは1個目のアウトカムだけが観測され、2個目のアウトカムが欠測の被験者の人数 - 7行目：
N_moは1個目のアウトカムが欠測で、2個目のアウトカムだけが観測された被験者の人数 - 8行目：
N_ooは両方のアウトカムが観測された被験者の人数 - 9～17行目：各観測になる確率の真値の計算。例えば、12行目の
p_om_1x_trueは二変量正規分布を確率楕円として描くと、次の左図の赤いゾーンの積分値に相当します。14行目のp_oo_00_trueは次の右図の赤いゾーンの積分値に相当します。
- 19～24行目：被験者のうちアウトカムが1だった人数。例えば、19行目の
Y_om_1xはN_omのうち1個目のアウトカムが1だった人数を表します。21行目のY_ooは、アウトカムが1だった人数ではなく、両方のアウトカムがそれぞれ(0,0) (1,0) (0,1) (1,1)
Stan実装は次になりました。
- 2～12行目：二変量正規分布の累積分布関数の計算になります。
z1やz2の足し算や掛け算を含むため、片方が無限大の場合にはpositive_infinity()（Stan内の無限大を表す関数）は使えず、代わりに十分大きな数を使用する必要があります。 - 25行目：ここでは切片項だけ推定します。説明変数を使いたい場合は、係数を
parametersブロックに書いて、transformed parametersブロックで説明変数と係数たちからxを作成すればOKです。 - 30～35行目：データ生成したRのプログラムとほぼ対応しています。30行目の
100は前述の部分と関連し、無限大の代わりに十分大きな数を使っています。 - 41行目：Stanの多項分布の引数には合計の被験者数は不要で、データの合計から自動で算出されます。引数の確率は縦ベクトルを渡す必要があるため、
[ ]で行ベクトルを作成して'で転置しています。
推定するには先ほどのRコードの27～30行目を実行します。推定結果の要約は次になりました。
> fit
variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
lp__ -785.31 -784.99 1.20 0.97 -787.71 -783.99 1.00 2045 2250
x[1] 0.28 0.29 0.05 0.05 0.20 0.37 1.00 3473 2473
x[2] 0.15 0.15 0.05 0.05 0.07 0.23 1.00 3203 2945
rho 0.64 0.64 0.06 0.06 0.53 0.73 1.00 3306 2734
p_om_1x 0.39 0.39 0.02 0.02 0.36 0.42 1.00 3473 2473
p_mo_x1 0.44 0.44 0.02 0.02 0.41 0.47 1.00 3203 2945
p_oo_00 0.45 0.45 0.02 0.02 0.41 0.48 1.00 2924 2117
p_oo_10 0.11 0.11 0.02 0.02 0.09 0.14 1.00 3690 2683
p_oo_01 0.16 0.16 0.02 0.02 0.13 0.19 1.00 3922 3165
p_oo_11 0.28 0.28 0.02 0.02 0.25 0.31 1.00 3140 2500
Enjoy!
Discussion