【読後メモ】データ解析のための統計モデリング入門（その2）

今更感ありますが、改めて「データ解析のための統計モデリング入門」を読みました(一部飛ばしたけど)。本記事は読後のメモです。
その2として一般化線形混合モデル(第7章)からMCMC法(第8章)までをまとめます

データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学)

作者: 久保拓弥
出版社/メーカー: 岩波書店
発売日: 2012/05/19
メディア: 単行本
購入: 16人クリック: 163回
この商品を含むブログ (29件) を見る

【目次】

はじめに
第7章一般化線形混合モデル（GLMM）
第8章マルコフ連鎖モンテカルロ（MCMC）法とベイズ統計モデル
ここまでのまとめ
参考文献

はじめに

「データ解析のための統計モデリング入門」を読んだので、どのようなことが書かれているのかといったメモを残しておきます。

この書籍では、統計モデリングとして「一般化線形モデル」に着目し、その基礎と発展について書かれています。よくある統計学の入門書では、確率の基礎や確率分布の紹介が多くの紙面を割いていることがあります。しかし本書では、一般化線形モデルに焦点を絞ることで、確率・統計の基礎についてはあまり触れられてはいません。これだけ書くと難しいように感じるかもしれませんが、逆に焦点を絞ることで本書を通して一貫した話題になっており、一般化線形モデルについて理解しやすいのではないかと思います。また、例題をベースにしているのも理解が進みやすいポイントかなと思います。

本書籍の読後メモは全3回でまとめる予定です。

本記事は第2回です。第2回では、第7章一般化線形混合モデル（GLMM）と第8章 MCMC アルゴリズムについてまとめていきます。

【トップに戻る】

第1回 : learning-with-machine.hatenablog.com

第3回 : Coming Soon

第7章一般化線形混合モデル（GLMM）

第6章までGLM（一般化線形モデル）を扱ってきました。 GLMは応用範囲が広く、多くの実例があります（特にロジスティック回帰はよく見る気がします）。しかし現実の問題では、GLMでは当てはまらない問題も多々あります。

そこで第7章では、また一つモデルを複雑にした「一般化線形混合モデル（GLMM）」を扱っています。

一般化線形混合モデル（GLMM）とは

GLMMとは、GLMの線形予測子に「ランダム効果」と書籍内で呼ばれている効果を付加したモデルです。

ランダム効果とは何かということを見ていくために、ロジスティック回帰を例にします。 GLMでのロジスティック回帰は以下の図のように表現されます（参照：【読後メモ】データ解析のための統計モデリング入門（その1） - 機械と学習する）。 f:id:hippy-hikky:20190220004228p:plain

これは、特徴量 $x$ が同じであれば、二項分布のパラメータである「事象が生起する確率 $q$ 」は同じであることを仮定しています。
しかしながら、現実には計測した特徴量が同じでもバラツキが生じることはあります。これを本書籍では「個体差」または「ランダム効果」と呼んでいます。計測していない、または、計測できていない特徴量を一纏めにしたものと見ることができます*1。

このランダム効果をモデルに組み込むと以下のようなモデルになります。

f:id:hippy-hikky:20190226002921p:plain

ここで、ランダム効果は $r_i$ と表記しているものです。 $r_i$ は個体 $i$ によって変わるものですが（「個体差」なので）、個体毎に予測するわけにはいきません。そこで、 $r_i$ は何らかの分布に従って生じると考えるのです。上記の例では、正規分布 $N(0, s^{2})$ に従うと仮定しています。

このように、ランダム効果を付加したGLMをGLMM（一般化線形混合モデル）と呼びます。
（何が「混合」なのかは次の節で）

GLMMのパラメータ推定

2.1に記載のモデルから、推定するべきパラメータは $\beta, r_i$ となりそうですが、 2.1でも述べたように、 $r_i$ を全て推定するのはモデルとしての意味がありません。そのため、ランダム効果は $r_i$ は何らかの分布に従って生じると考え、 $r_i$ について積分することで、予測対象から $r_i$ を消してしまいます。

${ \displaystyle L_i = \int_{-\infty}^{\infty} p(y_i | \beta, r_i) p(r_i | s) dr_i }$

これは、無限の $r_i$ を持った分布( $p(y_i | \beta, r_i)$ )を $r_i$ の分布( $p(r_i | s)$ )で重み付けして混ぜ合わせていると考えることができます。これで、予測するパラメータは、［ $\beta, r_i$ ］から［ $\beta, s$ ］になります。 2.1びGLMMの図にあるように、階層的なモデルになるわけです。

実際のパラメータの推定については、本書籍ではRを使って計算する方法が紹介されています。数秒検索したところ、Pythonで計算できるライブラリは無さそう。。。
Pythonで計算するなら、階層ベイズモデルと考えてMCMC アルゴリズム（後述）で計算するのが良いかもしれません。（この計算については、後日まとめたいと思います）

現実のデータ解析にはGLMMが必要

書籍では、過分散か否かでGLMMが必要かどうかを考えることもあるが、そもそものデータの取得方法でGLMで良いのかGLMMが必要なのかを見極める必要があると述べています。

GLMMは個体差や場所差などデータ取得の条件が異なる（擬似反復）ようにデータが構成される可能性があるなら GLMMで分析が必要とのことです。地域毎のアンケート調査のように、地域差が出る可能性がある場合には、場所差を考慮したモデルが必要ということなのかなと思います。

【トップに戻る】

第8章マルコフ連鎖モンテカルロ（MCMC）法とベイズ統計モデル

GLMMでは、計測仕切れない要因を考慮したモデルとして個体差や場所差を含んだモデルを考えてきました。個体差のようなものを個別にパラメータ推論するのは現実的ではないので、 2.2節に式を載せたように積分してその要因の期待値をとるような方法を学びました。個体差だけでなく場所差やもっと他の要因が絡むような場合、モデルはより複雑になり、 2.2節の積分が多重化していきます。すると、計算が困難となります。

このような複雑なモデルを解析的ではなく、数値的に近似計算するアルゴリズムとして「MCMC アルゴリズム(Markov Chain Monte Carlo method)」があります。

MCMC アルゴリズムとは

MCMC アルゴリズムとは、「多変量の確率分布からの乱数発生方法」と紹介されています。が、多変量でなくても良くて（書籍内の例でも最初は1変数）、ある確率分布から乱数を得るためのアルゴリズムと理解できます。

なぜ、「確率分布から乱数を得る」必要があるのでしょうか？
実際、乱数を得ることにはそこまで重要な意味はなく、 ”ある確率分布”というのが解析的には不明な場合にも、他の分布とMCMC アルゴリズムを利用することで不明な”ある確率分布”から乱数を得ることができるのです。乱数がたくさん得られれば、その平均が期待値になりますし、乱数の分布から”ある確率分布”の形も推定できます。

MCMC アルゴリズムの計算

MCMC アルゴリズムの大まかな流れは次の手順です。

解析対象のデータのモデルを決める
- 書籍では、20個体の植物それぞれから8個の種子を採取し、その生死のモデルを例にしています
- 8個中何個の種子が生きているかをモデル化するので二項分布になりますね
- 二項分布のパラメータ $q$ が不明で、この $q$ を求めたいわけです
モデルに基づいて尤度を決める
- モデルが決まれば、自ずと尤度が決まります
- 未知パラメータ $q$ の関数です( $L(q)$ , 対数をとった対数尤度 $\log{L(q)}$ が良く使われる)
パラメータ $q$ の初期値を適当に決める
$q$ を増やすか減らすかをランダムに決める
新しい $q^{new}$ で計算した $L(q^{new})$ と $L(q^{old})$ の値の大小に合わせて $q^{new}$ を採択するか $q^{old}$ を採択するかを決める
4に戻って繰り返す

これを多数回繰り返すことで、尤度 $L(q)$ の高くなる $q$ に近づいていき、 $q$ の最尤付近でランダムウォークするはずです。この $q$ の系列をMCMCサンプルと呼び、 $q$ の分布 $p(q|Y)$ （Yはデータ）が推定できます。

$p(q|Y)$ は尤度 $L(q)$ に応じてサンプルした $q$ から構成しており、尤度 $L(q)$ に比例する確率分布です。 $p(q|Y)$ はデータ $Y$ が与えられた際の $q$ がとる確率分布なので、この平均をとれば $q$ の期待値を推定できます。

MCMC アルゴリズムによってなぜ $p(q|Y)$ からのサンプリングが得られるのかについての詳細は別途ブログにまとめて行こうと思います。

MCMCサンプリングとベイズ統計モデル

ベイズ統計の枠組みで $p(q|Y)$ を捉えると、ベイズの定理から以下の式が求められます。

${ \displaystyle p(q|Y) = \eta p(Y|q)p(q) }$

ここで、 $\eta$ は規格化定数で右辺を確率の定義に当てはまるようにした定数です（積分が1になる）。右辺第1項はあるパラメタ $q$ が与えられた条件での $Y$ が生成される確率なので、尤度 $L(q)$ のことです。 $p(q)$ は事前分布と呼ばれていて、データがない時の $q$ の分布を表します。

ベイズ統計では、上記の式のように、知りたい未知のパラメータについての確率分布（これを「事後分布」と呼ぶ）を出力として得ます。事後分布の平均値、中央値、パーセンタイルなどを利用して予測値を活用していきます。パラメータには真の値があり、その推定値をある決まった値で算出するという頻度主義の考え方とは異なります。

先に示したベイズ統計の枠組みの式を上記の尤度などを使って書き換えると、次のようになります。

${ \displaystyle p(q|Y) \propto L(q) p(q) }$

事前分布 $p(q)$ が定数であれば3.2で示したように、 $p(q|Y)$ は尤度 $L(q)$ に比例するという式そのものです。ベイズ統計では、このように事前分布を加えることで、最尤法での計算に加えて尤度 $q$ の取りやすさや値の範囲などに追加の情報を使っているのです。事前分布 $p(q)$ を定数（一様分布）を利用すれば、最尤法と実質的に同じものになります。

【トップに戻る】

ここまでのまとめ

本記事では、「データ分析のための統計モデリング」の7,8章のメモをまとめてきました。

第7章では、一般化線形混合モデル(GLMM)と呼ばれる、未知の個体差を組み込んだモデルを扱いました。データとしては計測できていない（できない）要因がある場合に、 GLMMを利用して未知の要因を含めたモデリングが必要になってくるということでした。

第8章では、GLMMで扱う未知の要因が多いなど複雑なモデルを考える必要がある場合の計算方法として MCMC アルゴリズムを紹介しました。 MCMC アルゴリズムによると、未知のパラメータの事後分布からのサンプルが得られるということでした。 GLMMの枠組みで具体的な計算はここでは紹介されていませんが、第10章で階層ベイズモデルとして扱われます。

以上、今回は二つの章を扱いましたが、いよいよ現実的なモデルを扱うようになってきました。次回は、MCMC法を活用して、GLMのベイズモデル化、またGLMMのベイズモデル化（階層ベイズモデル）を扱います。