医療職からデータサイエンティストへ医療職からデータサイエンティストへ統計学、機械学習に関する記事をまとめています。ベイズ前回は階層性のあるデータに対して、線形混合モデルと階層ベイズを用いて解析を行いました。medi-data.hatenablog.com今回は、より複雑な階層ベイズモデルに挑戦していきます。前回は扱わなった、新たなデータに対しての予測も行っていきましょう。前回は、生徒50人のテストの点数と勉強時間のデータを用いて、4つの学校ごとに切片と傾きが異なる(点数と勉強時間の関係性が異なる)データを解析しました。今回は、3つの都道府県から18の学校を抽出し、そこから生徒450人をサンプリングした仮想データを用います。以下はサンプル抽出のイメージ図です。都道県ごと、そして学校ごとでも切片と傾きが異なることを想定しています。以下を実行してデータセット(st_df)を作成し、グラフ化してみましょう。"lm",se=F)上が都道府県ごとのグラフ、下が都道府県1の学校ごとのグラフです。ご覧の通り、都道府県ごと、学校ごとで切片と傾きが異なっているようです。このようにデータの関係性が複雑になってしまうと、前回の線形混合モデルでも扱うことはできなくなります。こんなデータでも解析することができます。階層ベイズの概略まずは、いつも通り概略から行きましょう。2層の階層ベイズでは、学校ごとに異なる切片と傾きは、全体に共通の回帰直線の周りで確率的に変動する値として推定しました。今回は、学校ごとの回帰直線を都道府県の回帰直線の周りで確率的に変動する値として、都道府県ごとの回帰直線を全体に共通の回帰直線の周りで変動する値として推定します。(ややこしい....)イメージとしては、まず都道府県ごとの回帰直線が、全体に共通の回帰直線の周りで変動する値として推定されます。学校ごと回帰直線が都道府県ごとの回帰直線の周りで確率的に変動する値として推定されるという順番です。ベイズ式とモデル式は以下のようになります。上記をstanでモデリングしていきます。}}少し複雑に見えますが、一つずつ見ていけばそんなに難しくないかと思います。これをstratified3.stanとして保存します。^;)Rのスクリプトファイルに戻り、stanを実行しましょう。dataを渡すときにprefだけは、学校数分の繰り返しになるので少し工夫が必要です。data,seed=123)無事に収束しました。それぞれのパラメーターの95%ベイズ信用区間を見てみます。31.09704a0、b0が回帰直線の全体平均、sigがバラツキになります。これを見ると学校ごとのバラツキ(sig_a,sig_b)が小さく、都道府県ごとのバラツキ(sig_ap,sig_bp)が大きいことが分かりますね。今回の結果から、都道府県によって勉強時間と点数の関係が大きく異なり、同じ県内の学校間ではそれほど大きな違いはないことが推測されます。全体傾きの信用区間に負の値が含まれていることから、必ずしも勉強時間と点数は正の相関ではない可能性があることも分かります。このように階層ベイズでは、パラメーターの確率分布から現象の解釈に幅を持たせることができます。次に全体平均の回帰直線をグラフに書いてみましょう。intercept,slope=slope))この直線の周りに都道府県ごとの回帰直線がバラツキつくことになります。予測さて、ここからは新たなデータに対しての予測を行ってみましょう。都道府県ごと、学校ごとの回帰直線全てに対して、予測ベイズ信用区間を求めることができるのですが、全部は大変なので都道府県2に属する学校15だけ抜き出して予測信用区間を求めてみます。学校15の勉強時間Xの最小値、最大値から新たなデータX_newを作ります。max(st_df$X[st_df$school=="15"]))次にX_newと推定された学校15のパラメーターから予測値をサンプリングします。学校ごとの切片や傾きは、行にサンプリング値、列にそれぞれの値が格納されたマトリックス形式になっているので、学校15は15列目に格納されています。#それぞれの値がマトリックス形式になっている。matrix(nrow=nrow(stf3_res$a),ncol=length(X_new))#学校15の予測値をX_newの値ごとにサンプリングする。)}最後にそれぞれのX_newごとの95%信用区間を求めて、グラフに反映します。alpha=0.2,fill="red")いい感じで予測ができているようですね!まとめ今回は複雑な関係性にあるデータを、階層ベイズを用いて解析しました。階層ベイズは予測もできますが、データ解釈の幅が広がるのが大きなメリットでしょう。前回のように超事前分布に特定の分布を指定すれば、前提とする仮定を組み込むこともできます。自由なモデリングができるのが、線形モデルや線形混合モデルではできないところなので、是非活用していきたいです!NGモデルこの記事を書くにあたってつまったところがあるので、参考までに紹介しておきます。stanでサンプリングした時のことです。30.93676んん!?b0とsig_bpの桁が明らかにおかしい....。データセットを作り直してみたり、seed値をかえてみたりで、3時間ほど格闘した末に犯人を見つけました。傾きb[k]のところがbpではなく、apになってます。皆さんも収束がおかしい場合は、stanファイルを確認してみることをお勧めします。※本記事は筆者が個人的に学んだことをまとめた記事なります。数学の記法や詳細な理論、筆者の勘違い等で誤りがあった際はご指摘頂けると幸いです。www.medi-08-data-06.work参考こちらの内容をアレンジして記事を書きました。これほどわかりやすくベイズモデリングを学べる書籍はほとんどないでしょう。見る有名なみどりぼんです。一般化線形モデルから、混合モデル、階層ベイズまでの繋がりを理解するのは、この本が一番だと思います。前回はベイズ推定を使って単回帰分析を行いました。学んだこと、日々の所感をまとめていきます。※本ブログは筆者が個人的に学んだこと感じたことをまとめたブログになります。所属する組織の意見・見解とは無関係です。medi_data0826引用をストックしましたストック一覧を見る閉じる引用するにはまずログインしてくださいログイン閉じる引用をストックできませんでした。再度お試しください閉じる限定公開記事のため引用できません。