一般化線形モデルについて
本記事は、データ解析のための統計モデリング入門(通称緑本)の3,4,5章を参考にしています。 (もし内容に不備等あればご一報ください)
データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学)
R, Pythonコード、使用しているデータはこちら
一般化線形モデルとは
一般化線形モデル(Generalized Linear Model / GLM)とは
一般線形モデルでは目的変数が正規分布に従うことを前提としているが、一般化線形モデルでは目的変数が正規分布に従わなくても適用でき、さらに質的変数であってもよい。また、目的変数と説明変数との関係式は簡単な線形式である必要はなく、以下のように表される。ただしは偏回帰係数を、は説明変数を、は説明変数の数を表す。
$ ここで関数 $g(y) = \beta_0x_0 + \beta_1x_1 + \cdots + \beta_nx_n はリンク関数と言う。ロジスティック回帰分析などに適用できる。一般化線型モデルとも言う。 統計webよりより g(\cdot)
要するに線形モデルのように、どういうデータにおいてもデータが正規分布に従っていると考え、直線で回帰するというのには無理があので、正規分布以外の確率分布も用いることでデータをうまく表現できるようにし、さらに線形ではない関係にも対応できるようにしたのが一般化線形モデル。
以下、例題※を用いて、一般化線形モデルであるポアソン回帰を実施する。
(※「データ解析のための統計モデリング入門」の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)
.dataframe tbody tr th:only-of-type { vertical-align: middle; } <div></div> .dataframe tbody tr th { vertical-align: top; } <div></div> .dataframe thead th { text-align: right; }
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
予測を基にポアソン回帰の統計モデルを立ててみる
モデルの作成
とりあえず、施肥効果
また、ある個体
(比較として種子が体サイズに依存しない一定モデル、種子が体サイズと施肥処理(
各モデルの関数
x モデル
一定モデル
xdモデル
この時それぞれの値は以下のように呼ばれる。
-
(切片)、\beta_1 (傾き):パラメーター\beta_2 -
の部分:線形予測子\beta_1 + \beta_2 x_i -
の\log{\lambda_i} の部分:リンク関数(今回は対数リンク関数と呼ばれる)\log
対数尤度
上記の
となる。 (一定モデル、xdモデルも同様)
対数リンク関数を用いることで「各因子の効果が掛け算で影響する」という効果をうまく反映できている。
最尤推定を実施する
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の正規分布における
の値をとる確率の2倍。-\infty〜0
(この値が大きいほどそのパラメータが0に近いと言える。検定におけるP値のようなもの)
ただし、説明変数をモデルに含めるべきかの判断は、**AIC(Akaike's information criterion)**を使う方が良いとの事。
AIC を使って最良のモデルを選択する
AIC(Akaike's information criterion)
以下の式で表される、良い予測をするモデルが良いモデルである(≠予測精度重視)との考えに基づいた指標。AICが低いモデルほど良いモデル
-
:最大対数尤度\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()
参考
関連書籍
- データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学)
- StanとRでベイズ統計モデリング (Wonderful R)
- Pythonによるベイズ統計モデリング: PyMCでのデータ分析実践ガイド
以上
Discussion