データ解析のための統計モデリング入門を読んでいる。
その読書メモ。
この章ではさまざまなGLMに触れる。
さまざまな種類のデータで応用できるGLM
GLMの特徴は、確率分布・リンク関数・線形予測子の組み合わせを指定することによってさまざまなタイプのデータを表現できること。
例題：上限のあるカウントデータ
二項分布を使ったGLMについて見ていく。
ポアソン分布は上限のないカウントデータを表現するのに使ったが、二項分布は上限のあるカウントデータに対して使う。
たとえば、「N個体の実験対象に同じ処理をしたら、y個体で反応が陽性、N-y個体で陰性だった」という構造のデータは二項分布で説明を試みる。
ちなみにポアソン分布は、応答変数が「0以上だけど上限がどこにあるのかわからないカウントデータ」を扱う。
二項分布の確率分布は、
$$
p(y | N, q) =
\begin{pmatrix}
N \\
y
\end{pmatrix}
q^{y}
(1-q)^{N-y}
$$
%matplotlib inline
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
df = pd.read_csv("data4a.csv")
df.describe()
|
N |
y |
x |
| count |
100.0 |
100.000000 |
100.000000 |
| mean |
8.0 |
5.080000 |
9.967200 |
| std |
0.0 |
2.743882 |
1.088954 |
| min |
8.0 |
0.000000 |
7.660000 |
| 25% |
8.0 |
3.000000 |
9.337500 |
| 50% |
8.0 |
6.000000 |
9.965000 |
| 75% |
8.0 |
8.000000 |
10.770000 |
| max |
8.0 |
8.000000 |
12.440000 |
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_c['x'], df_t['y'], label='T')
plt.legend()
plt.savefig("survive.png")
二項分布を図示してみる。
ロジスティック回帰とロジットリンク関数
GLMは確率分布・リンク関数・線形予測子を指定する統計モデルであり、
ロジスティック回帰では確率分布は二項分布、そしてリンク関数はロジットリンク関数を指定する。
二項分布では事象が生起する確率をパラメータとして指定する必要がある。 ロジスティック関数の関数形は、
$$
q_i = logistic(z_i) = \frac{1}{1 + exp(-z_i)}
$$
であり変数
は線形予測子
。
ロジスティック関数を変形すると、以下のようになる。
$$
logit(q_i) = log \frac{q_i}{1-q_i}
$$
GLM
import numpy as np
import statsmodels.api as sm
import statsmodels.formula.api as smf
fit = smf.glm(formula='y + I(N -y) ~ x + f', data=df, family=sm.families.Binomial())
res_fit = fit.fit()
res_fit.summary()
Generalized Linear Model Regression Results
| Dep. Variable: | ['y', 'I(N - y)'] | No. Observations: | 100 |
| Model: | GLM | Df Residuals: | 97 |
| Model Family: | Binomial | Df Model: | 2 |
| Link Function: | logit | Scale: | 1.0 |
| Method: | IRLS | Log-Likelihood: | -133.11 |
| Date: | Tue, 14 Mar 2017 | Deviance: | 123.03 |
| Time: | 07:09:30 | Pearson chi2: | 109. |
| No. Iterations: | 8 | | |
| coef | std err | z | P>|z| | [95.0% Conf. Int.] |
| Intercept | -19.5361 | 1.414 | -13.818 | 0.000 | -22.307 -16.765 |
| f[T.T] | 2.0215 | 0.231 | 8.740 | 0.000 | 1.568 2.475 |
| x | 1.9524 | 0.139 | 14.059 | 0.000 | 1.680 2.225 |
交互作用項の入った線形予測子
ここまで使ってきた線形予測子 
を更に複雑化して交互作用項を追加してみる。つまり積の効果を追加する。
交互作用項を加えた線形予測子は
$$
logit(q_i) = \beta_1 + \beta_2 x_i + \beta_3 f_i + \beta_4 x_i f_i
$$
となる。
fit = smf.glm(formula='y + I(N - y) ~ x * f', data=df, family=sm.families.Binomial())
res_fit = fit.fit()
res_fit.summary()
Generalized Linear Model Regression Results
| Dep. Variable: | ['y', 'I(N - y)'] | No. Observations: | 100 |
| Model: | GLM | Df Residuals: | 96 |
| Model Family: | Binomial | Df Model: | 3 |
| Link Function: | logit | Scale: | 1.0 |
| Method: | IRLS | Log-Likelihood: | -132.81 |
| Date: | Tue, 14 Mar 2017 | Deviance: | 122.43 |
| Time: | 07:28:17 | Pearson chi2: | 109. |
| No. Iterations: | 8 | | |
| coef | std err | z | P>|z| | [95.0% Conf. Int.] |
| Intercept | -18.5233 | 1.886 | -9.821 | 0.000 | -22.220 -14.827 |
| f[T.T] | -0.0638 | 2.704 | -0.024 | 0.981 | -5.363 5.235 |
| x | 1.8525 | 0.186 | 9.983 | 0.000 | 1.489 2.216 |
| x:f[T.T] | 0.2163 | 0.280 | 0.772 | 0.440 | -0.333 0.765 |
この例では交互作用項を増やしたモデルでは予測は悪くなっている。
交互作用項を使うときに注意すべきことは、むやみに交互作用項を入れないということ。
説明変数が多い場合には交互作用項の個数が組み合わせ論的爆発で増加してパラメータ推定が困難になる。
割算値の統計モデリングはやめよう
二項分布とロジットリンク関数を使ったロジスティック回帰を使う利点のひとつは、「種子が生存している確率」「処理に応答する確率」といった何かの世紀確率を推定するときに、(観測データ）/（観測データ）の割算値を作り出す必要がなくなること。
観測データどうしの割算によって生じる問題としては、以下のようなものがある。
- 情報が失われる
- 1000打数300安打の打者と10打数3安打の打者が同様に確からしい「三割打者」と扱われてしまう。
- 変換された値の分布はどうなる？
- 分子、分母にそれぞれ誤差の入った数量同士を割算して作られた割算地はどんな確率分布に従う？
割算をしないためのテクニック
例として面積が
である調査地
における人口密度を変数としたくなった際に、
$$
\lambda_i = A_i \times (人口密度） = A_i exp(\beta_1 + \beta_2 x_i)
$$
とモデル化する。
これを変形すると、
$$
\lambda_i = exp(\beta_1 + \beta_2 x_i + log A_i)
$$
となり、 
を線形予測子とする対数リンク関数、ポアソン分布のGLMになる。
確率密度と確率
正規分布など、連続値の確率分布では、確率密度を積分した値が確率と定義されている。
from scipy import stats
y = np.arange(-5, 5, 0.1)
plt.plot(pd.Series(stats.norm.pdf(y, loc=0, scale=1), index=y), label=r'$\mu=0, \sigma=1$')
plt.plot(pd.Series(stats.norm.pdf(y, loc=0, scale=3), index=y), label=r'$\mu=0, \sigma=3$')
plt.plot(pd.Series(stats.norm.pdf(y, loc=2, scale=1), index=y), label=r'$\mu=2, \sigma=1$')
plt.legend()
plt.savefig("Gaussian.png")
まとめ
- GLMでは応答変数のばらつきを表現する確率分布は正規分布だけでなく、ポアソン分布・二項分布・ガンマ分布などが選択できる
- 「N個の観察対象のうちk個で反応が見られた」というタイプのデータに見られるばらつきを表すために二項分布が使える
- 生起確率と線形予測子を結びつけるロジットリンク関数を使ったGLMのあてはめは、ロジスティック回帰と呼ばれる
- 線形予測子の構成要素として、複数の説明変数の積の効果をみる交互作用項が使える
- データ解析でしばしばみられる観測値同士の割算値作成や、応答変数の変数変換の問題点をあげ、ロジスティック回帰やオフセット項の工夫をすれば、情報消失の原因となる「データの加工」は不要になる
- 連続値の確率変数のばらつきを表現する確率分布としては、正規分布・ガンマ分布などがあり、これらを統計モデルの部品として使うときには、離散値の確率分布との違いに注意しなければならない
yusuke-ujitoko.hatenablog.com
