Blog

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

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

今回は第6章。この章では、ポアソン分布以外の確率分布の統計モデル、そしてGLMのオフセット項の使いかたについて書かれています。

※書籍の内容全てを網羅している訳ではありません、自分なりの解釈でまとめております。

 

6.2 例題:上限のあるカウントデータ

二項分布は上限のあるカウントデータをあらわすために使います。

この例題では、ある架空植物の個体 ( i ) それぞれにおいて「 ( N_{i} ) 個の観察種子のうち生きていて発芽能力のあるものは ( y_{i} ) 、死んだ種子は ( N_{i} – y_{i} ) 個」といった観察データが得られたとします。ある個体( i ) から得られた1個の種子が生きている確率を、生存確率 ( q_{i} ) とします。

説明変数としては、個体の大きさをあらわす体サイズ ( x_{i} ) 、施肥処理の有無 ( f_{i} ) となります。ここで調べたいことは、ある個体の生存確率( q_{i} ) が ( x_{i} ) や ( f_{i} ) といった説明変数によってどう変化するのかという点です。

例題データはこちらからダウンロードしてください。

まず、各列のまとめを表示します。

散布図をプロットします。

scatter

体サイズ ( x_{i} ) が大きくなると生存種子数 ( y_{i} ) が多くなるらしい。肥料をやると (( f_{i} = T )) 生存種子数 ( y_{i} ) が多くなるらしい。といったことが分かります。

 

6.3 二項分布で表現する「あり・なし」カウントデータ

この例題のように、「 ( N ) 個のうち ( y ) 個が生存していた」といった構造のカウントデータを統計モデルで表現するときには、二項分布がよく使われています。

 \(\left( \begin{array}{rr} N \\ y \\ \end{array} \right) \) は「場合の数」で、ここでは「\( N \) 個の観察種子の中から \( y \) 個の生存種子を選び出す場合の数」となります。

binom

6.4 ロジスティック回帰とロジットリンク関数

二項分布を使ったGLMのひとつであるロジスティック回帰の統計モデルについてです。

ロジスティック回帰では確率分布は二項分布、そしてリンク関数はロジットリンク関数を指定します。

ロジスティック関数の関数形は

\[ q_{i} = \rm{logistic} (z_{i}) = \frac{1}{1 + \exp (-z_{i})} \]

であり、変数 \( z_{i} \) は線形予測子 \( z_{i} = \beta_{1} + \beta_{2} x_{i} + \) …です。

作図してみます。

logistic

ロジスティック関数の逆関数であるロジット関数は次のようになります。

\[ \rm{logit} (q_{i}) = \log \frac{q_{i}}{1 – q_{i}} \]

この統計モデルをデータにあてはめてパラメーターを推定します。

Rでは目的変数にcbind(y, N-y)を取るようにしますが( \( y \) は生存していた種子数、 \( N-y \) は死んだ種子数)、statsmodelsでは、’y + I(N-y) ~ 説明変数’ のように記述することで実行できます。

おおよそ ( { \hat{\beta_{1}},\hat{\beta_{2}}, \hat{\beta_{3}} } = {-19.5, 1.95, 2.02 } ) となっています。これらの推定値を使って、予測の曲線を書きます。ほぼほぼRのpredict()関数と同じように使えますね。

predict

第4章で説明したAICによるモデル選択によって、これらのモデルの中から良い予測をするものを選びます。RであればMASS packageのstepAIC()関数というものを使えば、ネストしているモデルたちのAICを自動的に比較しながら、AIC最小のモデルを選択できるようですが、statsmodelsでは無いので、自力でやるしかなさそうです。

次に、交互作用を追加してみます。

predict

6.6 割算値の統計モデリングはやめよう

割算値、つまり変形した観測値を統計モデルの応答変数にするのは、不必要であるばかりでなく、まちがった結果を導きかねないものです。そこで、統計モデルを作る際にオフセット項わざを使います。

本では架空データを使ってポアソン回帰におけるオフセット項わざを説明していますが、ここでは省略してstatsmodelsの中のオフセット項の指定の方法だけ記述します。

 

6.7 正規分布とその尤度

正規分布は連続値のデータを統計モデルであつかうための確率分布です。Pythonで正規分布の密度関数を図示してみます。

Gaussian

6.8 ガンマ分布のGLM

ガンマ分布は確率変数のとりうる範囲が0以上の連続確率分布です。

gamma

ガンマ関数を使ったGLMです。

graph

This Post Has 0 Comments

Leave A Reply