データ解析のための統計モデリング入門を読んでいる。
その読書メモ。
データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・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列ごとの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()
これらからxとyの関係やfの効果は読み取れない。
ポアソン回帰の統計モデル
前章と異なり個体ごとに、ポアソン分布のパラメータが異なるとする。
まず最初に個体の体サイズ x だけに依存する統計モデルを考える。
その個体ごとに異なるパラメータを下記のようにおく。
$$
\lambda_i = exp(\beta_1 + \beta_2 x_i)
$$
対数を取ると、
$$
\log \lambda_i = \beta_1 + \beta_2 x_i
$$
となる。
この右辺は線形予測子と呼ばれる。
また左辺はリンク関数と呼ばれる。
β1とβ2を変化させた時のλは以下のようになる。
ポアソン回帰の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()
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 |
この出力を見ると、は1.29、は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))
説明変数が因子型の統計モデル
因子型の説明変数を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()
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()
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ととなり、少しだけよくなっている。