医療職からデータサイエンティストへ医療職からデータサイエンティストへ統計学、機械学習に関する記事をまとめています。stan前回はベイズ推定を使って単回帰分析を行いました。今回はさらにレベルを上げて階層性のあるデータ解析に挑戦していきます。medi-data.hatenablog.comwww.medi-08-data-06.work階層性のあるデータの概略から、通常の単回帰、線型混合モデルそして、階層ベイズモデルを用いたモデリングを行います。それではいきましょう。階層性のあるデータとはそもそも階層性のあるデータとはどのようなデータでしょうか?前回と同じようにテストの点数を勉強時間から推定する場合を考えてみましょう。今回のデータは生徒数50人のテストの点数と勉強時間を仮定した擬似データセットです。今回は4校から生徒をランダムにサンプリングしたと仮定しています。以下を実行してst_df(student_data_frame)というデータセットを作成してみましょう。勉強時間X、学校ID(1~4)、点数Yが並ぶ50行のデータセットです。学校ごとにサンプリングされた人数も異なっています。)まずは、勉強時間とテストの点数の散布図を書いてみます。school))どうやら勉強時間とテストの点数には相関がありそうです。通常の回帰分析を行います。626.5傾き15、切片-65の回帰直線になりましたね。グラフに追加してみましょう。"lm",se=F)ここまでは通常の回帰分析です。これで勉強時間が1時間増えるほど点数15点増えると言って良いのでしょうか?よくみると学校3に属する生徒に関しては傾きが緩いようです。そのことはモデルには反映されていません。これは、学校によらず勉強時間と点数との関係は一定であるという前提に基づいた解析を行ったからです。現実世界では基礎学力が高い生徒が集まっていたり、勉強をすればするほど点数が伸びる生徒が集まっていたり学校によって特徴があると考えられます。学校ごとに切片や傾きが異なるのが普通です。いった個人データの中に属性があるようなデータを階層性のあるデータと呼びます。今回は混合モデルと階層ベイズモデルの二つを使って階層性のあるデータを解析していきます。混合モデルまずは学校ごとに回帰直線をグラフに書いてみます。点線が通常の回帰分析の結果を表しています。予想通り、学校ごとで回帰直線の切片と傾きが大きく異なっているようです。ここで、混合モデルの登場です。通常の回帰分析は、各個人の値それぞれが独立して回帰直線の周りに散らばっているとして推定しています。イメージとしてはこんな感じです。式で表すともしくは、です。正規分布に従う個人差を足した値として表されます。混合モデルでは、この回帰直線から個人の値までのバラツキを個人差と学校差に分けて考えます。式で表すとつまり、回帰直線の切片と傾きにも正規分布に従う学校差によるバラツキを反映させるというのが、混合モデルになります。イメージとしてはこんな感じです。Rで混合モデルを実行するにはいくつかパケージがあるのですが、今回は{lme4}のlmer関数をを使います。属性)は属性のバラツキを説明変数の切片と傾きに反映させろという意味です。実際にやってみましょう。effectsとあるのが、切片と傾きのバラツキを表しています。今回の場合は切片の標準偏差が57ぐらい、傾きの標準偏差が5ぐらいなので、勉強時間0時間の基礎学力は学校によってだいたい57点のバラツキがあり、勉強時間と点数の関係は5点ぐらいはバラツキがあるということになります。effectsが全体の回帰式の切片と傾きになります。通常の回帰式と今回の結果を比べてみましょう。lmerで作成されたモデルは特殊な形式になっているので、値を取り出すにはsummaryを変数に代入するか、{lme4}のfixef関数を使います。fixef(mix_model)[2],col="black",size=1)混合モデルの結果では、切片が5.8、傾きが11.64と通常の回帰式よりも緩やかになりました。これは、学校差を考慮したことで、個人差の範囲では反映されなかった学校3の値がモデルに反映されたからでしょう。以上が混合モデルによる解析です。混合モデルを使った解析はまだまだ奥が深いので、そのうち記事にします。ベイズ推定さて、次はみなさんお待ちかね階層ベイズモデルです。概略から説明していきます。求めることを行いました。事前分布を仮定します。よって確率的に変動する値となります。この事前分布の事前分布(ややこしい...,)は超事前分布と呼び、分布が階層構造になっているため、階層ベイズモデルと呼ばれます。これを式で表すとこれは、先ほどの混合モデルの式と意味は同じですね。階層ベイズモデルは、切片と傾きの学校差を確率分布のバラツキとして表現したモデルになります。上記をstanで実装していきましょう。}}上記をstratified.stanとして保存します。Rスクリプトに戻り、以下を実行します。problemsという警告が出ましたが、収束に問題はないので無視してもいいみたいです。結果を確認してみましょう。切片が約8、傾きが約12、学校差による切片の標準偏差が110ぐらい、傾きの標準偏差が11ぐらいとなりました。混合モデルよりも標準偏差が大きくなっていますが、切片と傾きの推定値はだいたい同じになりました。グラフに加えてみましょう。赤色がベイズ推定の結果になります。回帰係数の切片と傾きのベイズ95%信用区間と混合モデルの95%信頼区間を比較してみます。混合モデルの信頼区間は{bloom}パッケージのtidyを使うと表示できます。17.241931階層ベイズモデルの方が幅が広く推定されているみたいですね。切片と傾きのバラツキ(sig_a,sig_b)も確率分布として表現できるのも階層ベイズの凄さです。傾きb0の95%信用区間に負の値が含まれているため、必ずしも勉強時間が増えたからと言って、点数が伸びる訳ではなさそうです。sig_a,sig_bの超事前分布には何も指定しなかったので、無情報分布として推定されましたが、他の分布を仮定することもできます。切片、つまり勉強時間が0時間の基礎学力のバラツキが2.9~429ととても広く推定されています。そこで、学校による基礎学力差のバラツキは100点ぐらいのはずだという仮定をモデルに組み込んでみます。するかというと、超事前分布に標準偏差100、自由度3or4の半t分布を使用します。(半t分布は標準偏差の事前分布として推奨されている)半t分布とはt分布の中で0以上の値のみを仮定する分布なのですが、sig_aで<lower=0>を仮定しているので、必然的に半t分布になります。(正規分布の場合はさらに仮定が強くなる)イメージとしては以下の図のような形で、0~100の間に多くの値が集まっていてて、遠くに離れるほどその確率密度が小さくなるような分布です。stanコードのmodelブロックに以下のように追加し、stratifed2.stanとして保存します。}}Rで実行してみます。33.796667sig_aのベイズ信用区間がだいぶ狭くなりました。切片平均a0の幅も小さくなっていますね。このように前提条件を仮定したモデルを作ることもできます。まとめ今回は階層性のあるデータに対して、線形混合モデルと階層ベイズモデルを適応して、結果を比較してみました。単縦な階層性であれば、線形混合モデルと階層ベイズでは同じような結果になることが確認できました。解析したいデータと仮定する前提条件などを考慮してモデル選択したいですね。※本記事は筆者が個人的に学んだことをまとめた記事なります。数学の記法や詳細な理論、筆者の勘違い等で誤りがあった際はご指摘頂けると幸いです。参考こちらの内容をアレンジして記事を書きました。これほどわかりやすくベイズモデリングを学べる書籍はほとんどないでしょう。見る有名なみどりぼんです。一般化線形モデルから、混合モデル、階層ベイズまでの繋がりを理解するのは、この本が一番だと思います。見る混合モデルを詳しく学びたい方は、こちらがおすすめです。混合モデルの概念から、結果の解釈、注意点などが網羅的に書かれています。学んだこと、日々の所感をまとめていきます。※本ブログは筆者が個人的に学んだこと感じたことをまとめたブログになります。所属する組織の意見・見解とは無関係です。medi_data0826引用をストックしましたストック一覧を見る閉じる引用するにはまずログインしてくださいログイン閉じる引用をストックできませんでした。再度お試しください閉じる限定公開記事のため引用できません。