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

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

データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学)

作者: 久保拓弥
出版社/メーカー: 岩波書店
発売日: 2012/05/19
メディア: 単行本
購入: 16人クリック: 163回
この商品を含むブログ (29件) を見る

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