はじめに

これはRICORA Advent Calendar 12 日目の記事になります。

みなさん初めまして。

さぶろくです。

今回は、データ解析のための統計モデリング入門という本を5章まで読んだので、4、5章をとりあえずまとめてみました。

統計学の基礎的な理論については、説明していません。

また、統計モデリングに今まで触れたことがなかったので、間違いが含まれている可能性があると思いますがご了承ください。

実行環境は、R version 4.1.0 です。

では、始めましょう。

本の概要

対象読者と章ごとの説明の流れを紹介します。

対象読者

作者曰く、この本の読者として対象としているのは、「数理モデルで現象を表現・説明する」基礎訓練を受けていない人達のようです。

大学で、数理モデリングや解析プログラミングの講義・演習を受けないまま、研究者になった後に複雑なデータ解析に取り組まなければならない人が増えています。

つまり、「今後データ解析をするだろうけど、やったことがない」と言った方にお勧めできる本のようです。

説明の流れ

説明の流れを本に載っていたようなフローチャート(まえがき p.vii)にしてみました。

この表を見てみると、ちょうど第5章のあたりまでがキリが良かったので、今回はここまでの内容をまとめていきます。

検定とモデル選択について

本の概要を話したところで本題に移りましょう。

これから話すのは、どのような統計モデルを選べば良いかと言った話です。
統計モデルとは、手に入れたデータを要約するようなものです。

全てのデータを持って来れば正しく要約できるので、どのモデルを選べば良いかなどと考える必要はないと思う方もいるでしょう。

ですが、ある事象の全てのデータを持ってくるというのは、データの欠損やそもそも個体数が多すぎて調べきれないなど困難です。

したがって、私たちは実際のデータから全てのデータをほぼ完璧に説明できるものが作れるように目指します。それが作れれば、調べられていないデータを大方予測できたり、それを見ることでデータの特徴を掴むことができたりします。

では、我々が作ったモデルが良いモデルであるかどうかを判断する方法は何でしょう？

その際に使われる方法が検定とモデル選択です。

今回、検定ではNeyman-Pearsonによる検定の枠組みをモデル選択ではAIC(Akaike's information criterion)を使います。

2つの違いについては、下図(p.96)のようです。
見てみると、統計モデルのパラメータの最尤推定計算までは同じのようです。

AIC

AICとは何であるか説明します。
その前に、逸脱度というものを定義します。
逸脱度とはあてはまりの良さを表す最大対数尤度に-2をかけたものであり、あてはまりの悪さを表す指標です。

このとき、最大対数尤度を$\log L^{*} $, 逸脱度を$D$とすると、

D = -2\log L^{*}

のように表せます。

ここで、AICとは予測の良さを重視する基準であり、AICが低いほど予測の良いモデルとなります。

先ほど定義した逸脱度$D$と最尤推定したパラメータ数$k$を用いて、

\begin{align}
AIC &= -2 \{(最大対数尤度) - (最尤推定したパラメータ数)\}\\
&= -2(\log L^{*} - k)\\
&= D + 2k
\end{align}

式を見る感じ、あてはまりの悪さの指標である逸脱度にパラメータ数を足して補正して予測の悪さを表現していることがわかります。

なぜ、この式が予測の悪さを表すかについて、数値実験的な説明はこの本の4章5節、数理統計学的な証明については坂本・石黒・北川『情報量統計学』あるいは小西・北川『情報量基準』を参照してください。

Neyman-Pearsonの検定の枠組み

Neyman-Pearsonの検定の枠組みの流れを追ってみましょう。

統計モデルにおける検定では、「帰無仮説は正しい」という命題が否定できるかどうか、その点だけを調べます。

まず、モデルのあてはまりの良さなどを表す検定統計量を作ります。

次に、帰無仮説が真のモデルであると仮定し、検定統計量の値がとりうる確率分布を調べ、「ありがちな範囲」を指定します。
この「ありがちな範囲」が95%である時、「有意水準を5%に設定した」と言います。

最後に、対立仮説のモデルで得られた統計量が「ありがちな範囲」の内部にあるか、外部にあるかを確認します。
ここで、外部にあれば、帰無仮説は棄却され、対立仮説が支持された言い、内部にあれば、帰無仮説は棄却されなかったと結論づけます。

以上がNeyman-Pearsonの検定の流れです。

今回は、Neyman-Pearsonの検定の枠組みの中でも、尤度比検定を扱います。
名前にある通り、この検定では尤度比といったものを使います。

尤度比とは、帰無仮説でのパラメータを代入した最大尤度を対立仮説でのパラメータの最大尤度で割ったものと定義されます。

ここで、

帰無仮説での最大尤度を$ {L_{1}}^{*}$

対立仮説での最大尤度を$ {L_{2}}^{*}$

とすると、それぞれの最大対数尤度を使って下のように表せます。

\frac{{L_{1}}^{*}}{{L_{2}}^{*}} = \frac{\exp({\log {L_{1}}^{*}})}{\exp(\log {L_{2}}^{*})}

ここで、上式の両辺の対数をとって、-2をかける。

つまり、逸脱度の差をとったもの

\begin{align}
\Delta D_{1,2} &= -2(\log {L_{1}}^{*} - \log {L_{2}}^{*})\\
&= D_{1} - D_{2}

\end{align}

を検定統計量として使います。

( $D_{1}$ と $D_{2}$ はそれぞれ帰無仮説、対立仮説の逸脱度とする。)

そして、この検定統計量がモデルを改善するか否かを判断する際に、先ほどの検定の流れに沿っていきます。

次の章で、具体的な例題に取り組みながら、検定のやり方について説明をします。

実際にやってみた

今回扱うデータは、架空植物100個体の種子数を集めたものです。

そして、作る統計モデルは、

$\lambda_{i} = \exp(\beta_{1} + \beta_{2}x_i)$を平均とする**ポアソン分布のGLM(一般線形モデル)**です。

この時、$\lambda_{i}$を種子数の平均, $x_{i}$を植物の体サイズとします。

(当ブログでは、GLMについて詳しく説明しません。詳しくは、3章の一般線形モデル(GLM)を参照してください。)

そして、比較するモデルは次の2つです。

一定モデル: $\lambda_{i}$が定数、$x_{i}$に依存しないモデル(傾き$\beta_{2}=0$, パラメータ数$k = 1$)
xモデル: $\lambda_{i}$が$x_{i}$に依存するモデル(傾き$\beta_{2}\not=0$, パラメータ数$k = 2$)

AICによるモデル選択

まず、データの読み込みをします。

d <- read.csv("data3a.csv")

次に、統計モデルを作ります。

一定モデルは、

fit1 <- glm(y.rnd ~ 1, data = d, family = poisson)

となり、

xモデルは、

fit2 <- glm(y.rnd ~ x, data = d, family = poisson)

となります。

ここで、それぞれのモデルの詳細を見てみよう。(fit1, fit2と実行すれば表示される。)

まず、一定モデルは下のようになります。

$\beta_{1}$の最尤推定量は2.058, AICが477.3であることがわかりました。

次に、xモデルは下のようになります。

$\beta_{1}$, $\beta_{2}$ の最尤推定量は1.29172, 0.07566、AICは474.8となります。

2つのモデルのAICを比較してみると、xモデルのAICの方が値が小さいので、

xモデルは一定モデルよりも予測の良いモデルと言えます。

以上で、AICによるモデル選択は終わりです。

統計モデリングの検定

2つの統計モデルを作るまでは、AICと同じである。
ただ、検定の枠組みを使うので、

帰無仮説:一定モデル(パラメータ数$k=1$, $\beta_{2}=0$)
対立仮説:xモデル(k=2, $\beta_{2}=0$)

のように設定する。

それでは、この２つのモデルの逸脱度の差を見てみよう。

fit1$deviance - fit2$deviance

検定の枠組みより、帰無仮説(一定モデル)が正しいという命題を否定できるかどうか確かめる。

つまり、逸脱度が4.5ぐらいではモデルが改善されていないと言って良いのか調べてみましょう。

ここで、検定における２種類の過誤について説明しましょう。

下の表のようになります。

↓帰無仮説は	帰無仮説は棄却	棄却できる
真のモデルである	第1種の過誤	(問題なし)
真のモデルではない	(問題なし)	第2種の過誤

表を見れば、わかりますがこの2つの誤り両方を避けることはできません。

そこで、Neyman-Pearsonの検定の枠組みでは第1種の過誤の検討にだけ専念します。

それでは、全体の流れを追っていきましょう。

帰無仮説である一定モデルが正しいものだと仮定する
観測データからパラメータ(今回だと $\hat{\beta_{1}}=2.06$ )が得られたので、これを真のモデルとほぼ同じであると仮定する。
このモデルからデータを何度も生成し、そのたびに$\beta_{2}=0(k=1)$と$\beta_{2}\not=0(k=2)$のモデルをあてはめれば、多くの$\Delta D_{1,2}$が得られるので$\Delta D_{1, 2}$の分布がわかる。
$\Delta D_{1,2} \geq 4.5$となる確率$P$($P$値)を評価します。

ここで、P値の扱いについて説明しましょう。

$P$値が「大きい」:$\Delta D_{1,2}=4.5$はよくあること→帰無仮説は棄却できない
$P$値が「小さい」:$\Delta D_{1,2}=4.5$はとても珍しいこと→帰無仮説を棄却しよう。残ったxモデルが正しいモデルだ！

そして、前の章の検定の流れより、$P$値を有意水準より大きいか、小さいかを判断するわけです。

この時、適当に有意水準$\alpha$を決めて、

$P\geq \alpha$:帰無仮説は棄却できない
$P<\alpha$:帰無仮説は棄却できる

これで、モデルの評価方法がわかりました。
今回は有意水準を5%として検定をします。

それでは、$P$値を計算するための、パラメトリックブートストラップ法(PB法)とカイ二乗分布を使った近似計算法を紹介しましょう。

1.パラメトリックブートストラップ法(PB法)

PB法とは、データを何度も生成する過程を、乱数発生を使ったシュミレーションによって実装されます。

それでは、PB法を使って、$\Delta D_{1,2}$の分布を作りましょう。

まず、Rの自作関数pb()を定義しましょう。

#PBの自作関数
get.dd <- function(d) #データの生成と逸脱度差の評価
{
    n.sample <- nrow(d) #データ数
    y.mean <- mean(d$y) #標本平均
    d$y.rnd <- rpois(n.sample, lambda = y.mean)
    fit1 <- glm(y.rnd ~ 1, data = d, family = poisson)
    fit2 <- glm(y.rnd ~ x, data = d, family = poisson)
    fit1$deviance - fit2$deviance #逸脱度の差を返す
}
pb <- function(d, n.bootstrap)
{
    replicate(n.bootstrap, get.dd(d))
}

この関数を呼び出してみましょう。

source("pb.R") #pb.Rを読み込み
dd12 <- pb(d, n.bootstrap = 10000)

R上での操作によって、$\Delta D_{1,2}$のサンプルが10000個を作り、dd12に格納されました。

それでは、$\Delta D_{1,2}$の分布をヒストグラムで見てみましょう。

hist(dd12, 100)
abline(v = 4.5, lty = 2)

点線は$\Delta D_{1,2}=4.5$となる位置です。
この線より、右側にある$\Delta D_{1,2}$はいくつあるでしょうか。

数えてみましょう。

sum(dd12 >= 4.5)

結果は324個と出ました。

つまり、$P$値は$324/10000=0.0324$となります。
これは、有意水準5%より小さいので帰無仮説は棄却されます。

したがって、残ったxモデルが採択されることになります。

2.カイ二乗分布を使った近似計算法

逸脱度の差$\Delta D_{1,2}$の確率分布というのは自由度1のカイ二乗分布に近似できる場合があります。
(詳しくはHoelの『統計理論入門』を参照すると良いらしい。)

そのような近似を行うためにRのanova関数を使います。

anova(fit1, fit2, test = "Chisq") 
#anovaは逸脱度を調べる関数
#Chisqでカイ二乗分布に近似

実行結果は以下のようになり、$P$値は0.034であるから、有意水準5%より小さいので帰無仮説は棄却されます。

したがって、この手法でもxモデルが採択されることがわかりました。

この方法はサンプルサイズが大きい場合に有効な近似計算ができるので、データ数が100くらいだと、PB法の方が優れていると言えるでしょう。

以上で、統計モデルの検定は終わりです。

結局どっちを使えば良いの？

ここまで、AICによるモデル選択と、統計モデルの検定を紹介しました。

では、どちらを使えば良いでしょう。

結論から言うと、この2つの比較方法は使う目的が違うので、どちらが良いかとは一概に言えません。

AICによるモデル選択は、「良い予測をするモデル」を選ぶと言う目的を持ちます。

Neyman-Pearsonの枠組みのもとでの統計的検定では、帰無仮説の安全な棄却を目的とします。
ここで注意しなければいけないのは、採択された対立仮説がどのような意味で良いかは明確ではありません。

どちらの比較方法を使った方が良いかについては、目的によって変わりそうですね。
これらの手法を使う際は、推定した統計モデルの意味合いを把握しておく必要がありそうです。

ちなみに、筆者は「推定された統計モデルの解釈は、それぞれの研究ごとに固有なものであり、分野ごとに異なる自然現象の捉え方に依存しているので、その文脈の中で検討すべき問題です。」と述べています。

まとめ

今回のブログのまとめは以下です。

未知のデータを予測したり、データの特徴を掴むために統計モデルが作られる。
作った統計モデルが"良いモデル"であるかどうかを確かめるために、AICによるモデル選択やNeyman-Pearsonの検定の枠組みに基づいたモデルの検定を使う。
モデル選択と統計モデルの検定はどちらが良いかは一概に言えない。推定された統計モデルの解釈を文脈の中で検討すべきである。

以上で、当ブログは終わりです。
統計モデリングの良い勉強になったので、6章以降もブログにしていこうと思います。
読んでいただきありがとうございました！

参考文献