Blog

Pythonで「データ解析のための統計モデリング入門」:第5章

Pythonで緑本(「データ解析のための統計モデリング入門」)の演習をしていきます。前記事はこちら

第5章は「GLMの尤度比検定と検定の非対称性」ということで、統計モデルの逸脱度の差に注目する尤度比検定についてです。

統計モデルの検定における「帰無仮説・対立仮説」が、AICによるモデル選択では「単純モデル・複雑モデル」となります。一定モデル(単純モデル)とxモデル(複雑モデル)の逸脱度の差が例えば4.5だったとして、尤度比検定では、検定統計量であるこの逸脱度の差が「4.5ぐらいでは改善されていない」と言ってよいのかどうかを調べます。逸脱度の差は以下の式で表せます

\[ \Delta D_{1,2} = -2 \times (\log L_{1} – \log L_{2}) \]

検定の全体の流れは以下のようになります。

(1)まずは帰無仮説である一定モデルが正しいものだと仮定する

(2)観測データに一定モデルをあてはめると、 \(\hat{\beta_{1}} = 2.06 \) となったので、これは真のモデルとほぼ同じと考えよう

(3)この真のモデルからデータを何度も生成し、そのたびに \( \beta_{2} = 0 (k = 1) \) と\( \beta_{2}\neq 0 (k = 2) \) のモデルをあてはめれば、たくさんの \( \Delta D_{1,2} \) が得られるから、 \( \Delta D_{1,2} \) の分布がわかるだろう

(4)そうすれば、一定モデルとxモデルの逸脱度の差が\( \Delta D_{1,2}\ge 4.5 \)となる確率 \( P \) が評価できるだろう

 

5.4.1 方法(1) 汎用性のあるパラメトリックブートストラップ法

パラメトリックブートストラップ法は、前述した検定の流れの(3)における「データを何度も生成し」の仮定を、乱数発生のシミュレーションによって実施する方法です。

この章の例題データのGLMによる推定結果は、一定モデルとxモデルそれぞれfit1とfit2に格納されているとします。

一定モデルとxモデルの逸脱度の差 \( \Delta D_{1,2}\ge 4.5 \) を計算します。

真のモデルから100個体ぶんのデータを新しく生成し、一定モデルとxモデルをこの新データにあてはめてみます。

今回は、逸脱度の差が0.94となりました。この1ステップによって「一定モデルが真のモデルである世界」での逸脱度の差がひとつ得られます。このステップを1000回ほど繰り返すと、「統計検定量の分布」、この例題でいうと「逸脱度の差 \( \Delta D_{1,2} \) の分布」を予測できます。

このPB法を実行するために、自作関数pb()を定義します。

以上の操作によって、逸脱度の差 \( \Delta D_{1,2} \) のサンプルが1000個つくられてdd12に格納されました。実は \( \Delta D_{1,2} \) のサンプル個数は \( 10^3 \) ぐらいでは十分ではなく、\( 10^4 \) かそれ以上にしたほうが良いとのことです。

概要を調べます。

ヒストグラムで図示、そして逸脱度の差が4.5がどのあたりにくるのかを示します。

itsudatsu

合計1000個ある \( \Delta D_{1,2} \) のうちいくつぐらいが、この4.5より右にあるのでしょうか?

1000個中の31個が4.5より大きい、つまり確率は31/1000、すなわち \( P = 0.0031 \) ということになります。ついでに、\( P = 0.05 \) となる逸脱度の差 \( \Delta D_{1,2} \) を調べてみると、

となり、有意水準5%の統計学的検定のわくぐみのもとでは、 \( \Delta D_{1,2}\leq 3.78\) ぐらいまでは「よくある差」とみなされます。

この尤度比検定の結論としては、「逸脱度の差4.5の \( P \) 値は0.031だったので、これは有意水準0.05よりも小さい」ので有意差が有り、「帰無仮説(一定モデル)は棄却され、xモデルが残るのでこれを採択」と判断します。

 

5.4.2 方法(2)  \( \chi^2 \) 分布を使った近似計算法

近似計算法を使うと、もっとお手軽に尤度比検定ができる場合があります。

逸脱度の差 \( \Delta D_{1,2} \) の確率分布は、自由の1の \( \chi^2 \) 分布に近似できる場合があります。Rではanova()関数を使い、anova(fit1, fit2, test=”Chisq”) といったふうにすれば \( P \) 値が得られますが、statsmodelsのanova関数(statsmodels.stats.anova.anova_lm())はGLMの結果に対応していません。

ですが、つまり5.4.1で得た逸脱度の差の確率分布を自由の1の \( \chi^2 \) 分布に近似して作って、その後は同様に計算すればいけるってことですよね・・・?

得られた\( P \) 値は0.0036です。ただ、\( \chi^2 \) 分布近似はサンプルサイズが大きい場合に有効な近似計算なので、\( P \) 値はあまり正確ではない可能性があります。

調査した個体数が多くない小標本のもとでは、PB法を使って逸脱度差の分布をシミュレーションで生成するのがいいでしょう、とのことです。


This Post Has 0 Comments

Leave A Reply