読者です 読者をやめる 読者になる 読者になる

【データ解析のための統計モデリング入門】3章 一般化線形モデル(GLM)

データ解析のための統計モデリング入門を読んでいる。
その読書メモ。

3章では、説明変数を組み込んだ統計モデルを扱う。

本書で提供されているcsvファイルをいじっていく。
なおデータをいじる際にpandasを使った。
チュートリアルとして下記が有用だった。

10 Minutes to pandas — pandas 0.17.0 documentation GitHub - ajcr/100-pandas-puzzles: 100 short data puzzles for pandas, ranging from very easy to super tricky (55% complete)

観測されたデータの概要を調べる

import pandas as pd

df = pd.read_csv("data3a.csv")

# 冒頭5行を取得
df.head(5)

# 末尾5行を取得
df.tail(5)

# データの型の確認
df.info()

# 統計量を確認
df.describe()

# データの型の確認
type(df)

# 列のデータ型の確認
type(df['x'])

# x列の取得
df['x']

# y列の取得
df['y']

# f列の取得
df['f']

統計モデリングの前にデータを図示する

横軸にx列、縦軸にy列をとった散布図をかく

# fの値でdfを分ける
df_c = df[df.f == 'C']
df_t = df[df.f == 'T']

plt.scatter(df_c['x'], df_c['y'], label="C")
plt.scatter(df_t['x'], df_t['y'], label="T")

plt.legend()
plt.show()

f:id:yusuke_ujitoko:20170311205856p:plain

f列ごとのy列の箱ひげ図をかく

# fの値でdfを分ける
df_c = df[df.f == 'C']
df_t = df[df.f == 'T']

plt.boxplot([df_c['y'], df_t['y']])
plt.show()

f:id:yusuke_ujitoko:20170311210617p:plain

これらからxとyの関係やfの効果は読み取れない。

ポアソン回帰の統計モデル

前章と異なり個体ごとに、ポアソン分布のパラメータが異なるとする。
まず最初に個体の体サイズ x だけに依存する統計モデルを考える。
その個体ごとに異なるパラメータを下記のようにおく。
{} $$ \lambda_i = exp(\beta_1 + \beta_2 x_i) $$ 対数を取ると、 {} $$ \log \lambda_i = \beta_1 + \beta_2 x_i $$ となる。
この右辺は線形予測子と呼ばれる。
また左辺はリンク関数と呼ばれる。

β1とβ2を変化させた時のλは以下のようになる。 f:id:yusuke_ujitoko:20170311235148p:plain

ポアソン回帰のGLMで対数リンク関数を用いるのは、「推定計算に都合が良く」かつ「わかりやすい」ため。
推定計算に都合がよい点は、らい数リンクが常に非負である点。
わかりやすさは、要因の効果が積で表されるため。

GLMのあてはめをしてみる

import statsmodels.api as sm
import statsmodels.formula.api as smf

fit = smf.glm(formula='y ~ x', data=df, family=sm.families.Poisson())
fit.fit().summary()
Generalized Linear Model Regression Results
Dep. Variable: y No. Observations: 100
Model: GLM Df Residuals: 98
Model Family: Poisson Df Model: 1
Link Function: log Scale: 1.0
Method: IRLS Log-Likelihood: -235.39
Date: Sun, 12 Mar 2017 Deviance: 84.993
Time: 14:16:42 Pearson chi2: 83.8
No. Iterations: 7


coef std err z P>|z| [95.0% Conf. Int.]
Intercept 1.2917 0.364 3.552 0.000 0.579 2.005
x 0.0757 0.036 2.125 0.034 0.006 0.145

この出力を見ると、{\beta_1}は1.29、{\beta_2}は0.0757であることがわかる。
std errは標準誤差であり、推定値のバラつきを標準偏差で表したもの。
zはz値と呼ばれる統計量であり、最尤推定値をSEで割ったもの。
これによってWald信頼区間というものを構成でき、推定値がゼロから十分離れているかどうかの粗い目安になる。
P|z|は平均がz値の絶対値、標準偏差1の正規分布における、マイナス無限大からゼロまでの値を取る確率の2倍。
この確率が大きいほどz値がゼロに近くなる。

Log-Likelyhoodは最大対数尤度であり、あてはまりの良さを表す。
この値が大きいほど、あてはまりが良いとする。

ポアソン回帰モデルによる予測

# fの値でdfを分ける
df_c = df[df.f == 'C']
df_t = df[df.f == 'T']

plt.scatter(df_c['x'], df_c['y'], label="C")
plt.scatter(df_t['x'], df_t['y'], label="T")

# 予測
x = np.arange(np.min(df['x']), np.max(df['y']))
plt.plot(x, np.exp(1.29 + 0.0757 * x))

f:id:yusuke_ujitoko:20170312143149p:plain

説明変数が因子型の統計モデル

因子型の説明変数をGLMであつかうために、因子型の説明変数をダミー変数に置き換えて立式する。

import statsmodels.api as sm
import statsmodels.formula.api as smf

fit = smf.glm(formula='y ~ f', data=df, family=sm.families.Poisson())
fit.fit().summary()
Generalized Linear Model Regression Results
Dep. Variable: y No. Observations: 100
Model: GLM Df Residuals: 98
Model Family: Poisson Df Model: 1
Link Function: log Scale: 1.0
Method: IRLS Log-Likelihood: -237.63
Date: Sun, 12 Mar 2017 Deviance: 89.475
Time: 14:34:10 Pearson chi2: 87.1
No. Iterations: 7


coef std err z P>|z| [95.0% Conf. Int.]
Intercept 2.0516 0.051 40.463 0.000 1.952 2.151
f[T.T] 0.0128 0.071 0.179 0.858 -0.127 0.153

最大対数尤度は -237.63になっていて、あてはまりはさっきより悪くなっている。

説明変数が数量型+因子型の統計モデル

import statsmodels.api as sm
import statsmodels.formula.api as smf

fit = smf.glm(formula='y ~ x + f', data=df, family=sm.families.Poisson())
fit.fit().summary()
Generalized Linear Model Regression Results
Dep. Variable: y No. Observations: 100
Model: GLM Df Residuals: 97
Model Family: Poisson Df Model: 2
Link Function: log Scale: 1.0
Method: IRLS Log-Likelihood: -235.29
Date: Sun, 12 Mar 2017 Deviance: 84.808
Time: 14:36:20 Pearson chi2: 83.8
No. Iterations: 7


coef std err z P>|z| [95.0% Conf. Int.]
Intercept 1.2631 0.370 3.417 0.001 0.539 1.988
f[T.T] -0.0320 0.074 -0.430 0.667 -0.178 0.114
x 0.0801 0.037 2.162 0.031 0.007 0.153

最大対数尤度は-235.29ととなり、少しだけよくなっている。

まとめ

  • 一般化線形モデル(GLM)はポアソン回帰やロジスティック回帰など、いくつかの制約を満たしている統計モデルたちの総称
  • 統計モデルを作るにはデータの図示がとても大切
  • GLMは確率分布・リンク関数・線形予測子を指定する統計モデル
  • 統計モデルの因子型の説明変数は、ダミー変数という考え方で理解できる
  • GLMでは数量型・因子型の両タイプの説明変数を同時に組み込んでよく、またそのときに対数リンク関数を使っていると説明変数の効果が、それぞれの積として表現できるので理解しやすい。
  • GLMの設計では、データを上手く表現できる確率分布を選ぶという発想なので、「何でも正規分布」といった考え方から脱却できる。

yusuke-ujitoko.hatenablog.com