本記事は、データ解析のための統計モデリング入門（通称緑本）の3,4,5章を参考にしています。
（もし内容に不備等あればご一報ください）
R, Pythonコード、使用しているデータはこちら
一般化線形モデルとは
一般化線形モデル（Generalized Linear Model / GLM）とは
一般線形モデルでは目的変数が正規分布に従うことを前提としているが、一般化線形モデルでは目的変数が正規分布に従わなくても適用でき、さらに質的変数であってもよい。また、目的変数と説明変数との関係式は簡単な線形式である必要はなく、以下のように表される。ただしは偏回帰係数を、は説明変数を、は説明変数の数を表す。
$$g(y) = \beta_0x_0 + \beta_1x_1 + \cdots + \beta_nx_n$$
ここで関数$g(\cdot)$はリンク関数と言う。ロジスティック回帰分析などに適用できる。一般化線型モデルとも言う。
統計webよりより
要するに線形モデルのように、どういうデータにおいてもデータが正規分布に従っていると考え、直線で回帰するというのには無理があので、正規分布以外の確率分布も用いることでデータをうまく表現できるようにし、さらに線形ではない関係にも対応できるようにしたのが一般化線形モデル。
以下、例題※を用いて、一般化線形モデルであるポアソン回帰を実施する。
（※「データ解析のための統計モデリング入門」の3章より）
例題
架空の植物100個体からなる集団を調査していて、各個体の種子数を数えたものがデータだとする。
また、各個体の体サイズや各個体への施肥処理の有無もデータとして記録しているものとする。
条件は以下（解析データファイルは data3a.csv）
- 植物個体：$i$
- 各個体における種子の数：$y_i$
- 各個体の体サイズ：$x_i$（正の実数）
- 各個体への施肥処理の有無：$f_i$（$T$：処理あり, $C$：処理なし）
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
%matplotlib inline
# データの読み込み
data = pd.read_csv("data3a.csv",dtype={"y": int,"x":float, "f":"category"})
data.head(5)| y | x | f | |
|---|---|---|---|
| 0 | 6 | 8.31 | C |
| 1 | 6 | 9.44 | C |
| 2 | 6 | 9.50 | C |
| 3 | 12 | 9.07 | C |
| 4 | 10 | 10.16 | C |
データの観察
# データ全体を見る
data_c = data.loc[data['f']=='C']
data_t = data.loc[data['f']=='T']
plt.scatter(data_c["x"], data_c["y"], color='red', label='C：施肥処理なし')
plt.scatter(data_t["x"], data_t["y"], color='blue', label='T：施肥処理あり')
plt.xlabel('体サイズ')
plt.ylabel('種子の数')
plt.legend(loc='upper left')
plt.show()
作成した散布図から以下の予想が立てられる。
- 体サイズ $x$ の増加に伴って、種子数 $y$ が増えているように見える
- 施肥処理の効果 $f$ は無さそう
予測を基にポアソン回帰の統計モデルを立ててみる
モデルの作成
とりあえず、施肥効果 $f$ は無視して、ある個体 $i$ の種子数が $y_i$ である確率は、ポアソン分布に従っていると仮定する。
$$P(y_i|\lambda_i) = \frac{\lambda_i^{y_i}\exp(-\lambda_i)}{y_i !}$$
また、ある個体 $i$ の平均種子数 $\lambda_{i}$ が、体サイズ $x_i$ に関与しているという見立てから次の関数（$x$モデル）を定義する 。
（比較として種子が体サイズに依存しない一定モデル、種子が体サイズと施肥処理($d_i$)の両方に依存するxdモデルも用意する）
各モデルの関数
x モデル
$$ \lambda_i = \exp(\beta_1 + \beta_2 x_i) $$
$$ \log{\lambda_i} = \beta_1 + \beta_2 x_i $$
一定モデル
$$ \lambda_i = \exp(\beta_1) $$
$$ \log{\lambda_i} = \beta_1 $$
xdモデル
$$ \lambda_i = \exp(\beta_1 + \beta_2 x_i + \beta_3 d_i) $$
$$ \log{\lambda_i} = \beta_1 + \beta_2 x_i + \beta_3 d_i$$
$$d_i = \left\{
\begin{array}{ll}
1 & (f_i=C) \\
0 & (f_i=T)
\end{array}
\right.
$$
この時それぞれの値は以下のように呼ばれる。
- $\beta_1$（切片）、$\beta_2$（傾き）：パラメーター
- $\beta_1 + \beta_2 x_i$ の部分：線形予測子
- $\log{\lambda_i}$ の $\log$ の部分：リンク関数（今回は対数リンク関数と呼ばれる）
対数尤度
上記の $\lambda_i$ の定義を用いると、この 統計モデルの対数尤度は下記のようになる。
$$ \log{L(\beta_1,\beta_2)} = \sum_i{\log{\frac{\lambda_i ^{y_i}\exp(-\lambda_i)}{y_i !}}} $$
$$(\lambda_i = \exp(\beta_1+\beta_2 x_i))$$
となる。
（一定モデル、xdモデルも同様）
対数リンク関数を用いることで「各因子の効果が掛け算で影響する」という効果をうまく反映できている。
$$ \lambda_i = \exp(\beta_1 + \beta_2 x_i) $$
最尤推定を実施する
statsmodelを使用
import statsmodels.api as sm
import statsmodels.formula.api as smf
# ポアソン回帰の推定量の導出
model = smf.glm('y ~ x', data=data, family=sm.families.Poisson())
result = model.fit()
# 詳細な結果の確認
print(result.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.0000
Method: IRLS Log-Likelihood: -235.39
Date: Sat, 12 Feb 2022 Deviance: 84.993
Time: 15:31:47 Pearson chi2: 83.8
No. Iterations: 4
Covariance Type: nonrobust
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
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
==============================================================================以上の結果から最尤推定値は以下のように求まった。
- $\hat{\beta_1} = 1.29$
- $\hat{\beta_2} = 0.0757$
各パラメータの説明
- std err：標準誤差（最尤推定値がどれくらいバラつくかの指標）
- z（value）：最尤推定値をSEで割った値。この値からWald統計量を算出して、推定値がゼロから十分離れているか確認できる。
- P>|z|：平均|z|、標準偏差1の正規分布における $-\infty〜0$ の値をとる確率の2倍。
（この値が大きいほどそのパラメータが0に近いと言える。検定におけるP値のようなもの）
ただし、説明変数をモデルに含めるべきかの判断は、AIC（Akaike's information criterion）を使う方が良いとの事。
AIC を使って最良のモデルを選択する
AIC（Akaike's information criterion）
以下の式で表される、良い予測をするモデルが良いモデルである（≠予測精度重視）との考えに基づいた指標。AICが低いモデルほど良いモデル
$$AIC = -2(\log{L^\ast}-k)$$
- $\log{L^\ast}$：最大対数尤度
- $k$：予測に用いたパラメーターの数
パラメーター数が多い統計モデルほどデータへの当てはまりが予測精度、最大対数尤度）がよくなるが、これはたまたま得られた観測データへの当てはまりの良さを示しており、真のモデルに対する当てはまりの良さではない。
AICではパラメーター数の情報を指標に組み込む事で、観測データに比べて真のモデルに対する当てはまりが悪くなってしまうことを防いでいる。
それぞれのモデルのAICを比較してみる
# AIC の計算
# 一定モデル
model_const = smf.glm('y ~ 1', data=data, family=sm.families.Poisson())
result_const = model_const.fit()
print('AIC(一定モデル):',round(result_const.aic,2))
# xモデル
model_x = smf.glm('y ~ x', data=data, family=sm.families.Poisson())
result_x = model_x.fit()
print('AIC(xモデル):',round(result_x.aic,2))
# xdモデル
model_xd = smf.glm('y ~ x + f', data=data, family=sm.families.Poisson())
result_xd = model_xd.fit()
print('AIC(xdモデル):',round(result_xd.aic,2))結果
- 一定モデル: 477.29
- xモデル: 474.77
- xdモデル: 476.59
（xdモデル：体サイズと施肥効果の両方を考慮したモデル)
結果は以上のようになり、xモデル（体サイズ）のみを考慮したモデルが最も良いモデルだと結論づけられる。
結果のプロット
# ポアソン回帰の推定結果を使って、回帰直線を図示
x = np.linspace(7, 13, 100)
y = np.exp(result.params["Intercept"] + x*result.params["x"])
data_c = data.loc[data['f']=='C']
data_t = data.loc[data['f']=='T']
plt.scatter(data_c["x"], data_c["y"], color='red', label='C：施肥処理なし')
plt.scatter(data_t["x"], data_t["y"], color='blue', label='T：施肥処理あり')
plt.plot(x,y,"r--")
plt.xlabel('体サイズ')
plt.ylabel('種子の数')
plt.legend(loc='upper left')
plt.show()
参考
関連書籍
