「データ解析のための統計モデリング入門」はRで書かれた本ですが、これをpythonで書き直しながら読み進めていきます。
内容の理解と読み進めることが目的ですので、ベストなコードではないことをご理解ください。
たくさんの人がやっているような内容だとは思いますが、ご了承ください。

前回

確率分布と統計モデリングの最尤推定

2. 一般化線形モデル(GLM)

2.1 例題：個体ごとに平均種子数が異なる場合

第1章の例題と似ているが、個体の大きさがさまざまであったり、あるいは植物たちを二つのグループに分けて、
それぞれで異なる実験処理を施していたりする点で異なる。
架空植物100個体の調査をして得られた、個体ごとの種子数のデータがあるとする。
植物個体$i$の種子数は$y_i$個であり、体サイズ$x_i$が観測されている。
さらに、全個体のうち50個体($i\in\{1,2,\cdots,50\}$)は何も処理していない(処理C)、
残り50体($i\in\{51,52,\cdots,100\}$)には肥料を与える処理(施肥処理、処理T)を施す。

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

import pandas as pd

d = pd.read_csv('kubobook_2012/data/data3a.csv')
d

d.dtypes

y      int64
x    float64
f     object
dtype: object

d.describe()

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

視覚的に把握するためにも、いろいろな図にしてよく見ることは大事である。

import matplotlib.pyplot as plt

for s in ['C', 'T']:
    d_select = d[d['f']==s]
    plt.scatter(d_select['x'], d_select['y'], label=s)
plt.legend(loc='upper left')
plt.xlabel('x')
plt.ylabel('y')

import seaborn as sns

sns.boxplot(x='f', y='y', data=d);

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

最初に個体$i$の体サイス$x_i$だけに依存する統計モデルについて考える。
説明変数は$x_i$で、応答変数は種子数$y_i$である。$f_i$はとりあえず無視する。
ある個体$i$において種子数$y_i$である確率$p(y_i|\lambda_i)$はポアソン分布に従っていて、

$$
p(y_i|\lambda_i)=\frac{\lambda_i^{y_i}\exp{(-\lambda_i)}}{y_i!}
$$

と仮定する。

個体ごとに異なる平均$\lambda_i$を説明変数$x_i$の関数として定義する。
ここでは、ある個体$i$の平均種子数$\lambda_i$が、
$$
\lambda_i=\exp{(\beta_1+\beta_2x_i)}
$$
であるとする。
$\beta_1$や$\beta_2$をパラメータと呼び、とりあえず$\beta_1$を切片、$\beta_2$を傾きと呼ぶことにする。
ここで$\lambda_i=\exp{(\beta_1+\beta_2x_i)}$を図示すると。

import numpy as np

def lmd(x, b1, b2):
    return np.exp(b1+b2*x)

x = np.arange(-3,6,0.1)

plt.plot(x, lmd(x, -2, -0.8), label='{b1,b2}={-2,-0.8}');
plt.plot(x, lmd(x, -1,  0.4), label='{b1,b2}={-1,0.4}');
plt.legend();
plt.xlabel('x')
plt.ylabel('lambda')

このモデルの平均種子数$\lambda_i$の式は、
$$
\log{\lambda_i}=\beta_1+\beta_2x_i
$$
と変形できる。
このとき右辺の$\beta_1+\beta_2x_i$は線形予測子と呼ばれる。($\beta_1$と$\beta_2$の線形結合)
上の式は、$\log{\lambda_i}=(線形予測子)$となっているが、左辺の関数はリンク関数と呼ばれる。
今回の場合は、対数関数が指定されているため、対数リンク関数と呼ばれる。
ポアソン回帰をする場合、たいていは対数リンク関数を使用する。
また、数学的によい性質があるので、ポアソン分布・二項分布の正準リンク関数と呼ぶ。

ポアソン回帰とは、観測データに対するポアソン分布を使った統計モデルのあてはめであり、
この統計モデル対数尤度$\log{L}$が最大になるパラメータ$\hat{\beta}_1$と$\hat{\beta}_2$の推定値を決めることである。
データ$\boldsymbol{Y}$のもとでの、このモデルの対数尤度は、
$$
\log{L}(\beta_1,\beta_2)=\sum_i\log{\frac{\lambda_i^{y_i}\exp{(-\lambda_i)}}{y_i!}}
$$
となり、$\lambda_i$は$\beta_1$と$\beta_2$の関数であることに注意する。

import statsmodels.api as sm

model = sm.GLM(d.y, sm.add_constant(d.x), family=sm.families.Poisson())
#model.fit_regularized(method='elastic_net', alpha=1.0, L1_wt=0.0) # 正則化
result = model.fit()

result.summary()

# sklearnでも用意されている
from sklearn.linear_model import PoissonRegressor

glm = PoissonRegressor(
                        alpha=0,  # ペナルティ項
                        fit_intercept=True,  # 切片
                        max_iter=100,  # ソルバーの試行回数
                      )
glm.fit(d.x.values.reshape(-1, 1), d.y.values.reshape(-1, 1))

print(glm.coef_ )
print(glm.intercept_  )

[0.07566189]
1.2917210225621196

constは切片$\beta_1$、説明変数$x$の係数は傾き$\beta_2$に対応している。
std errorはパラメータの標準誤差の推定値で、 標準誤差(SE) とは推定値$\beta_1$と$\beta_2$のばらつきを標準偏差で表したものである。
簡単に「同じ調査方法で同数の別データをとりなおしたときのばらつき具合」ということにする。
zはz値と呼ばれる統計量であり、最尤推定値をSEで除した値である。
これにより、Wald信頼区間というものが構成できて、推定値がゼロから十分に離れているかの荒い目安となる。
P>|z|は、この値(確率)が大きいほどz値がゼロに近くなり、$\hat{\beta_1}$と$\hat{\beta_2}$がゼロに近いことを表現する。

ある説明変数をモデルに含めるべきか否かといった判断は、Wald信頼区間ではなく、4章で説明するモデル選択を使用するほうが良いかもしれない。
モデル選択はより良い選択をする統計モデルを探し出すもので、当てはまりの改善ではなく、予測の改善を目的としているからである。

ここで、最大対数尤度を当てはまりの良さと呼ぶ。
このモデルの最大対数尤度は-235.4ということがわかる。

ここで、ポアソン回帰の結果を図示する。

xx = np.arange(np.min(d.x), np.max(d.x), 0.5)
yy = result.predict(sm.add_constant(xx))

for s in ['C', 'T']:
    d_select = d[d['f']==s]
    plt.scatter(d_select['x'], d_select['y'], label=s)
plt.plot(xx, yy, '-', color='red')
plt.legend(loc='upper left')
plt.xlabel('x')
plt.ylabel('y')

2.5 説明変数がカテゴリ変数の統計モデル

カテゴリ変数はダミー変数に置き換えてあつかう。
施肥効果$f_i$だけが影響するモデルの平均値を、

$$
\lambda_i=\exp{(\beta_1+\beta_3d_i)}
$$

として、$\beta_1$は切片、$\beta_3$は施肥の効果を表す。
ここで、$d_i$は以下のような値をとるダミー変数である。

$$
d_i=\left\{
\begin{array}{11}
0\ (f_i=C)\\
1\ (f_i=T)
\end{array}
\right.
$$

個体$i$が肥料なし($f_i=C$)の場合、
$$
\lambda_i=\exp{(\beta_1)}
$$

個体$i$が施肥処理をした($f_i=T$)、
$$
\lambda_i=\exp{(\beta_1+\beta_3)}
$$
となる。

d.loc[d['f']=='C', 'f_dummy'] = 0
d.loc[d['f']=='T', 'f_dummy'] = 1

model = sm.GLM(d.y, sm.add_constant(d.f_dummy), family=sm.families.Poisson())
result = model.fit()

result.summary()

個体$i$が肥料なし($f_i=C$)の場合、
$$
\lambda_i=\exp{(2.05+0)}=\exp{(2.05)}=7.77
$$

個体$i$が施肥処理をした($f_i=T$)、
$$
\lambda_i=\exp{(2.05+0.0128)}=\exp{(2.0628)}=7.87
$$
となり、肥料をやると平均種子数は少しだけ増えると予測される。
最大対数尤度は-237.63と$x_i$だけのモデルより当てはまりが悪くなっている。

2.6 説明変数が数量型+カテゴリ型の統計モデル

体サイズ$x_i$と施肥効果$f_i$の複数の説明変数を同時に組み込んだ統計モデルを作る。
$$
\log{\lambda_i}=\beta_1+\beta_2x_i+\beta_3d_i
$$

model1 = sm.GLM(d.y, sm.add_constant(d[['x', 'f_dummy']]), family=sm.families.Poisson())
result1 = model1.fit()

result1.summary()

$f_i$だけのモデルでは効果(係数)がプラスであったが、このモデルではマイナスとなっている。
最大対数尤度に関しては、$x_i$だけのモデルと比較すると、少しだけ当てはまりがよくなっている。

「対数kリンク関数では要因の効果が積で表される」という性質を持つ。
個体$i$が肥料なし($f_i=C$)の場合、
$$
\lambda_i=\exp{(1.26+0.08x_i)}
$$

個体$i$が施肥処理をした($f_i=T$)、
$$
\lambda_i=\exp{(1.26+0.08x_i-0.032)}
$$
となり、以下のように分解できる。
$$
\begin{align}
\lambda_i&=\exp{(1.26)}\exp{(0.08x_i)}\exp{(-0.032)}\
=&(定数)×(サイズの効果)×(施肥処理の効果)
\end{align}
$$

サイズ$x_i$について「説明変数$x_i$が1増加すると、$\lambda_i$は$\exp{(0.08×1)}=1.08$倍に増える」と予測される。
施肥$f_i$については「説明変数$d_i$が1の場合、$\lambda_i$は$\exp{(-0.032×1)}=0.969$倍になる」と予測される。

ここで、平均が線形予測子に等しい恒等リンク関数を使ってポアソン回帰を行ってみる。

model2 = sm.GLM(d.y, sm.add_constant(d[['x', 'f_dummy']]), family=sm.families.Poisson(link=sm.genmod.families.links.identity))
result2 = model2.fit()

result2.summary()

x = np.arange(5, 20, 0.1)

def pred(x, d, result):
    xx = sm.add_constant(x)
    dd = np.array([d]*len(xx))
    X = pd.DataFrame(np.concatenate([xx.reshape(-1,2),dd.reshape(-1,1)], axis=1))
    y_pred = result.predict(X)
    return y_pred

fig = plt.figure(figsize=(12,4))

ax1 = fig.add_subplot(1, 2, 1)
ax1.plot(x, pred(x, 1, result1), label='T');
ax1.plot(x, pred(x, 0, result1), label='C');
ax1.legend()
ax1.set_title('log link')

ax2 = fig.add_subplot(1, 2, 2)
ax2.plot(x, pred(x, 1, result2), label='T');
ax2.plot(x, pred(x, 0, result2), label='C');
ax2.legend()
ax2.set_title('identity link')

2.7 「何でも正規分布」「何でも直線」には無理がある

GLMにおいて、恒等リンク関数を指定すると、一般化ではない線形モデルあるいは一般線形モデルと呼ばれる。
直線回帰はLMのあてはめのひとつである。
直線回帰はGLMの一部であり、統計モデルの特徴は以下のようになる。

観測値{$x_1,x_2,\cdots,x_n$}と{$y_1,y_2,\cdots,y_n$}のペアがあり、$\boldsymbol{X}=\{ x_i\}$を説明変数、$\boldsymbol{Y}=\{ y_i\}$を応答変数と呼ぶ。
$\boldsymbol{Y}$は平均$\mu_i$で標準偏差$\sigma$の正規分布に従う
あるデータ点$i$において平均値が$\mu_i=\beta_1+\beta_2x_i$となる

今回の例のように応答変数がカウントデータだったとすると、以下のような疑問点が生じる。

正規分布は連続的な値を扱う確率分布だったはずでは？
カウントデータなのに平均値の予測がマイナスになる理由は？
図で見ると「ばらつき一定」ではなさそうなのに、分散一定を仮定？

このように直線回帰の統計モデルは「現実ばなれ」している。
これに対してポアソン分布を使ったGLMでは、上に挙げた問題点は、

ポアソン分布を使っているのでカウントデータに正しく対応
対数リンク関数を使えば平均値は常に非負
$y$のばらつきは平均とともに増大する

のようにうまく解決できるように見える。
また、応答変数を対数変換して直線回帰することと、ポアソン回帰はまったく別物であることに注意する。

最後に応答変数を対数変換して直線回帰したものと、ポアソン回帰の結果の比較をする。

model_poisson = sm.GLM(d.y, sm.add_constant(d[['x', 'f_dummy']]), family=sm.families.Poisson())
result_poisson = model_poisson.fit()

result_poisson.summary()

model_linear = sm.OLS(np.log(d.y), sm.add_constant(d[['x', 'f_dummy']]))
result_linear = model_linear.fit()

result_linear.summary()

x = np.arange(6, 15, 0.1)

fig = plt.figure(figsize=(12,4))

ax1 = fig.add_subplot(1, 2, 1)
for s in ['C', 'T']:
    d_select = d[d['f']==s]
    ax1.scatter(d_select['x'], d_select['y'], label=s, alpha=.5)
ax1.plot(x, pred(x, 1, result_poisson), label='T');
ax1.plot(x, pred(x, 0, result_poisson), label='C');
ax1.set_title('GLM')
ax1.legend()

ax2 = fig.add_subplot(1, 2, 2)
for s in ['C', 'T']:
    d_select = d[d['f']==s]
    ax2.scatter(d_select['x'], d_select['y'], label=s, alpha=.5)
ax2.plot(x, np.exp(pred(x, 1, result_linear)), label='T');
ax2.plot(x, np.exp(pred(x, 0, result_linear)), label='C');
ax2.set_title('OLS')
ax2.legend()

同じような結果が得られたが、今回の場合は上記の理由によりポアソン回帰の方が適当と考えられる。

次回

GLMのモデルの選択

参考書