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

【データ解析のための統計モデリング入門】6章 GLMの応用範囲をひろげる

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

この章ではさまざまな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")

f:id:yusuke_ujitoko:20170313215716p:plain

二項分布を図示してみる。

f:id:yusuke_ujitoko:20170313215710p:plain

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

GLMは確率分布・リンク関数・線形予測子を指定する統計モデルであり、 ロジスティック回帰では確率分布は二項分布、そしてリンク関数はロジットリンク関数を指定する。

二項分布では事象が生起する確率をパラメータとして指定する必要がある。 ロジスティック関数の関数形は、 {} $$ q_i = logistic(z_i) = \frac{1}{1 + exp(-z_i)} $$

であり変数{z_i}は線形予測子{z_i = \beta_1 + \beta_2 x_i + \cdots}
ロジスティック関数を変形すると、以下のようになる。 {} $$ 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

交互作用項の入った線形予測子

ここまで使ってきた線形予測子 {\beta_1 + \beta_2 x_i + \beta_3 f_i }を更に複雑化して交互作用項を追加してみる。つまり積の効果を追加する。 交互作用項を加えた線形予測子は

{} $$ 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安打の打者が同様に確からしい「三割打者」と扱われてしまう。
  • 変換された値の分布はどうなる?
    • 分子、分母にそれぞれ誤差の入った数量同士を割算して作られた割算地はどんな確率分布に従う?

割算をしないためのテクニック

例として面積が{A_i}である調査地{i}における人口密度を変数としたくなった際に、 {} $$ \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) $$ となり、 {z_i = \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")

f:id:yusuke_ujitoko:20170313230357p:plain

まとめ

  • GLMでは応答変数のばらつきを表現する確率分布は正規分布だけでなく、ポアソン分布・二項分布・ガンマ分布などが選択できる
  • 「N個の観察対象のうちk個で反応が見られた」というタイプのデータに見られるばらつきを表すために二項分布が使える
  • 生起確率と線形予測子を結びつけるロジットリンク関数を使ったGLMのあてはめは、ロジスティック回帰と呼ばれる
  • 線形予測子の構成要素として、複数の説明変数の積の効果をみる交互作用項が使える
  • データ解析でしばしばみられる観測値同士の割算値作成や、応答変数の変数変換の問題点をあげ、ロジスティック回帰やオフセット項の工夫をすれば、情報消失の原因となる「データの加工」は不要になる
  • 連続値の確率変数のばらつきを表現する確率分布としては、正規分布・ガンマ分布などがあり、これらを統計モデルの部品として使うときには、離散値の確率分布との違いに注意しなければならない

yusuke-ujitoko.hatenablog.com