データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学)を読んだのでメモしておきます。勉強を始めたばかりで間違ってる箇所があるかもしれません。
統計モデリングのコツ
まず最初に以下の記事にある統計モデリングについてのコツを胸に刻んでおきます。
シンプルなモデルから始めるということ
これは、いきなり最終的に実現したい複雑なモデルを構築しようとするのではなく、簡単なモデルから始めてだんだんと拡張していく方がやりやすいということです。
作成したモデルに対していきなり現実のデータを入れるのではなくて、必ずシミュレーションデータへの適用を行うということ
シミュレーションで作成された理想的なデータに対して、モデルがうまくパラメータを推定できるのは当たり前じゃないかと考える人もいると思います。
R で一般化線形モデル: glm() 関数
| 確率分布 | 乱数生成 | パラメーター推定 | |
|---|---|---|---|
| (離散) | ベルヌーイ分布 | rbinom() | glm(family = binomial) |
| 二項分布 | rbinom() | glm(family = binomial) | |
| ポアソン分布 | rpois() | glm(family = poisson) | |
| 負の二項分布 | rnbinom() | glm.nb() in library(MASS) | |
| (連続) | ガンマ分布 | rgamma() | glm(family = gamma) |
| 正規分布 | rnorm() | glm(family = gaussian) |
# glm( モデル式, 確率分布の指定(リンク関数の指定), data.frameの指定)
fit <- glm( y ~ log.x, family = poission(link = "log", data = d)
- モデル式(線形予測子 z): どの説明変数を使うか
- link関数： zと応答変数(y)平均値の関係は？
- family: 度の確率分布を使うか
1. データを理解するために統計モデルを作る
統計モデルの基本部品である「確率分布」を使うとさまざまな「ばらつき」「欠損」を表現できる。また、さまざまな種類の「誤差」もモデル化することができる。
2. 確率分布と統計モデルの最尤推定
確率分布（probability distribution）は統計モデルの本質的な部品であり、データに見られるさまざまな「ばらつき」を表現する。
確率分布とは確率変数の値とそれが出現する確率を対応させたもの。
確率分布にはさまざまなものがあるので、データの特徴に合わせて確率分布を選ばなければならない。
どのような確率分布を使った統計モデルであっても「データに対するあてはまりの良さ」を対数尤度で表すことができ、最尤推定とは対数尤度を最大にするようなパラメーターを探す出すことである。
今得られたデータにあてはまるようなパラメーターを探しだすのが推定、次に得られるデータへのあてはまりを重視するのが予測である。
まず「この統計モデリングでは、このような理由でこの確率分布を使いました」と他人にきちんと説明できるようになることが重要。
ポアソン分布
データが離散値、ゼロ以上の範囲、上限とくになし、平均$\approx$分散
ポアソン分布で指定できるパラメーターはひとつだけであり、それは分布の平均である。
ポアソン分布
p(y|\lambda) = \frac{\lambda^y exp(-\lambda)}{y!}
このp(y|λ)平均がλであるときに、ポアソン分布に従う確率変数がyという値になる確率である。
ポアソン分布の特徴
- すべてのyについて和をとると1になる。
- 率分布の平均はλである。
- 分散と平均は等しい。(λ=平均=分散)
二項分布
データが離散値、ゼロ以上で有限の範囲、分散は平均の関数
正規分布
データが連続値、範囲が$[-∞, +∞]$、分散は平均とは無関係
ガンマ分布
データが連続値、範囲が$[0, +∞]$、分散は平均の関数
最尤推定法
最尤推定法は尤度という「あてはまりの良さ」を表す統計量を最大にするようなパラメーターの値を探そうとするパラメーター推定法である。
尤度の実体は、あるλの値を決めたときにすべての個体iについての$p(y_i | \lambda)$の積である。
尤度はパラメーター関数のため、$L(\lambda)$とすると以下のように定義される。
L(\lambda) = p(y_1 | \lambda) \times p(y_2 | \lambda) \times \dots \times p(y_m | \lambda) \\
= \prod_i p(y_i | \lambda) = \prod_i \frac{\lambda^y exp(-\lambda)}{y_i!}
なぜ積になるのかというと「$y_1$が2である」かつ「$y_2$が２である」かつ ... とi個の事象が同時に真である確率を計算したいからである。
確率分布の選び方
データを見たら次の点に注意してみる。
- 説明した医療は【離散】か【連続】か
- 説明した医療の【範囲】は？
- 説明した医療の標本分散と標本平均の関係は？
3. 一般化線形モデル(GLM)
個体ごとに異なる説明変数(個体の属性)によって平均種子数が変化する統計モデルを一般化線形モデル(GLM)という。
ポアソン回帰とは、観測データに対するポアソン分布を使った統計モデルのあてはめ(fitting)であり、この統計もdるの対数尤度$logL$が最大になるパラメーター$\beta_1, \beta_2$の推定値を決めること。
このモデルでの対数尤度は以下のとおりである。
logL(\beta_1, \beta_2) = \sum_i log \frac{\lambda_i^{y_i} exp(-\lambda_i)}{y_i !}
線形予測子は$log\lambda_i = \beta_1 + \beta_2 x_i$なので$\lambda_i$は$\beta_1$と$\beta_2$の関数となっている。
パラメーター推定ができればどんなモデルでもいいという発想はやめて、数式が現象をどのように表現しているのかという点に注意しながら統計モデルを設計する必要がある。
応答変数yの構造に合わせて適切な確率分布を選ぶことが重要。
GLMは確率分布、リンク関数、線形予測しを指定する統計モデルであり、Rのglm()関数でパラメーターを推定できる。
4. GLMのモデル選択 -AICとモデルの予測の良さ-
あてはまりの悪さ=逸脱度
D = -2logL^*
AICは統計モデルのあてはまりの良さ(goodness of fit)ではなく、予測の良さ(goodness of prediction)を重視するモデル選択基準である。
AIC = -2{ (最大対数尤度) - (最尤推定したパラメーター数) } \\
= -2(logL^* - k) \\
= D + 2k
このAICが「一番小さいモデルが良いモデル」となる。
5. GLMの尤度比検定と検定の非対称性
尤度比検定では下記「逸脱度の差」を検定統計量としてあつかう。
\Delta D_{1,2} = -2 × (logL_1^* - logL_2^*)
Neyman-Pearsonの統計学的検定のわくぐみでは、パラメーター数の少ないモデルを帰無仮説と位置づけ、帰無仮説が棄却できるかどうかの確率評価に専念する。
- 帰無仮説：一定モデル(パラーメーター数 k=1, $\beta_2=0$)
- 対立仮説：xモデル(k=2, $\beta_2\neq0$)
6. GLMの応用範囲を広げる -ロジスティック回帰など-
二項分布の確率分布の定義
p(y|N, q) = \left(\begin{array}{cc} N \\ y \end{array}\right) q^y(1-q)^{N-y} \\
p(y|N, q) ... N個中のy個で事象が発生する確率 \
$ \left(\begin{array}{cc} N \ y \end{array}\right) $ ... N個の観測種子の中からy個の生存種子を選び出す場合の数
ロジスティック関数の定義
q_i = logistic(z_i) = \frac{1}{1+exp(-z_i)} \\
ここで変数z_iは、線形予測子z_i=\beta_1+\beta_2 x_i + \dots
logistic <- function(z) 1 / (1 + exp(-z)) # 関数の定義
z <- seq(-6, 6, 0.1)
plot(z, logistic(z), type = "l")
ロジット関数(logit function)の定義
(ロジット関数はロジスティック関数の逆関数である)
logit(q_i) = log \frac{q_i}{1 - q_i} = z_i
正規分布は連続値のデータを統計モデルで扱うための確率分布である。
正規分布はガウス分布と呼ばれることもある。
正規分布では平均値のパラメーターを$\mu$(ミュー)で、値のばらつきを表す標準偏差をパラメーター$\sigma$(シグマ)で表す。
平均パラメーター$\mu$と標準偏差パラメーター$\sigma$の正規分布の数式表現は以下のようになる。
p(y|\mu, \sigma) = \frac{1}{\sqrt{2\pi\sigma^2}}exp\left(\begin{array} -\frac{(y-\mu)^2}{2\sigma^2} \end{array}\right)
これは正規分布の確率密度関数となる。
連続値の確率分布では、確率密度を積分した量が確率として定義されている。
ガンマ分布は確率変数の取りうる範囲が0以上の連続確率分布であり、確率密度関数は以下のとおりである。
p(y|s,r) = \frac{r^s}{\Gamma(s)}y^{s-1}exp(-ry)
ここでsはshapeパラメーター、rはrateパラメータと呼ばれる。
7. 一般化線形混合モデル(GLMM)
GLMMとは「観測されなかった個体差」を組み込んだGLMである。
現実の観測データでは「説明変数意外は全て均質」といった条件は基本的に満たされないため、平均λは全体共通とはならない。そのため、説明変数が同じであるなら平均も同じというGLMの家庭が成立しないため、期待値は大きくばらつく。
このような「人間が測定できない、測定しなかった個体差」を組み込んだモデルを一般化線形混合モデルという。これはデータのばらつきは「二項分布・ポアソン分布」で、個体のばらつきは「正規分布」であらわすような複数の確率分布を部品とする統計モデルである。
「全個体の生存種子数の分布は、ただひとつの二項分布で説明できる」と期待していたのに失敗してしまった。これを過分散(=ばらつきが大きい)という。
個体iの個体差を表すパラメーターを$r_i$とするとロジット関数は以下のようになる。
logit(q_i) = \beta_1 + \beta_2 x_i + r_i
8. マルコフ連鎖モンテカルロ(MCMC)法とベイズ統計モデル
MCMCアルゴリズムでは、一つのステップの中で前の状態$q$に基づいて新しい状態
$q^{new}$を作り出しているので、マルコフ連鎖(Markov chain)と呼ばれる。
また一般に、乱数を利用した計算アルゴリズムはモンテカルロ法(Monte Carlo methd)と呼ばれている。
MCMCアルゴリズムの目的は、何か特定の値の探索ではなく、ステップ数とともに変化するパラメーターの値の生成である。これをサンプリングと呼びます。
メトロポリス法(Metropolis method)
- パラメーター$q$の初期値を選ぶ。
- $q$を増やすか減らすかをランダムに決める。
- $q^{new}$において尤度が大きくなる(あてはまりが良くなる)なら$q$の値を$q^{new}$に変更する。
- $q^{new}$で尤度が小さくなる(当てはまりが悪くなる)場合であっても、確率$r$で$q$の値を$q^{new}$に変更する。このときの確率$r$は尤度比$r=\frac{L(q^{new})}{L(q)}$とする。
サンプリングされたヒストグラムで示される確率分布をがマルコフ連鎖の定常分布である。
MCMCサンプリングは統計モデルを観測データにあてはめる方法のひとつであり、その結果として与えられたデータとモデルの下でのパラメーターの確率分布が得られる。
しかし、パラメーターの確率分布という考え方は頻度主義のわくぐみでは表現数ることができないのでベイズ統計学で使う統計モデルとして表現する必要がある。
種子の生存確率$q$をベイズの公式で書くと以下のようになる。
p(q|Y) = \frac{p(Y|q) p(q)}{\sum_q p(Y|q) p(q)}
この式の左辺の$p(q|Y)$はデータ$Y$が得られた時に$q$が従う確率分布でベイズ統計学ではこれを事後分布とよぶ。
右辺の$p(Y|q)$は$q$の値が決まっている時にデータ$Y$が観測される確率である。この問題の場合、二項分布の積である尤度$L(q)$がそれに相当するため、$p(Y|q)=L(q)$となる。また、この尤度の後ろに付いている$p(q)$はデータ$Y$がないときの$q$の確率分布で、ベイズ統計モデルでは事前分布とよぶ。
ベイズ統計モデルは以下の様な構造をもつ統計モデルである。
事後分布 = \frac{尤度×事前分布}{データが得られる確率} \propto 尤度×事前分布
### ギブスサンプリング
ギブスサンプリングの利点
- 各MCMCステップにおいてもとの値と更新された値の相関がより小さい。
- MCMCサンプリングの詳細を指定しなくても良い
参考
GLM から始める統計モデリング
glm関数の使い方 ... glm()関数, AIC
実践 統計モデリング入門 【1. 概要・目次】