Python: statsmodelsでガンマ回帰

これまで読み進めてきた緑本を改めて復習するために、statsmodelsで一般化線形モデル(GLM)を作り、糖尿病に関するデータの回帰分析を行います。

実装はこちらです。

blog/statsmodels_gamma.ipynb at master · ohke/blog · GitHub

データセット

今回はscikit-learn付属の糖尿病データセットを使います。

import pandas as pd
from sklearn.datasets import load_diabetes

diabetes = load_diabetes()

data = pd.DataFrame(diabetes.data, columns=("age", "sex", "bmi", "map", "tc", "ldl", "hdl", "tch", "ltg", "glu"))
data['y'] = diabetes.target

データは442レコード有り、各レコードは11個の変数を持ちます。
1年後の糖尿病の進行状況が目的変数 $y$ (値が大きいほど進行していることを表す)となっており、それ以外は説明変数です。説明変数は全て平均0で正規化されています。この10個の説明変数から1年後の進行状況 $y$ を予測する、回帰分析のタスクとなります。

f:id:ohke:20180317113506p:plain

各変数の相関係数を見ます。
$y$ との相関係数では、 $bmi$ (BMI)が最も高く、 $ltg$ (おそらく中性脂肪)が続いています。
また、 $age$ (年齢)は各変数と緩やかに相関があります。加齢とともに $tc$ (総コレステロール)や $ltg$ が増加している傾向は見て取れます。

f:id:ohke:20180317114039p:plain

$y$ の分布も確認しておきます。平均152、標準偏差77となっています。

f:id:ohke:20180317120129p:plain f:id:ohke:20180317120139p:plain

GLMの設計

それでは $y$ を予測するモデルをGLMで設計します。

GLMでは、確率分布・リンク関数・線形予測子を決める必要がありました (1/25と2/12の投稿も参考にしてください)。

0以上で上限の無い連続値のようですので、確率分布にはガンマ分布を使います
ガンマ分布ですので、リンク関数には対数リンク関数を使います

ガンマ分布は、以下の確率密度関数で表される0以上の連続確率分布です。
$\Gamma(s)$ はガンマ関数で、 $s$ が正の整数なら $\Gamma(s)=(s-1)!$ となります。ガンマ分布の平均は $s\mu$ 、分散は $s\mu^2$ です。

$$ p(y \mid s,\mu)=\frac{\exp(-\frac{y}{\mu})}{\Gamma (s)\mu^s}y^{s-1} $$

線形予測子で指定するパラメータは以下の5パターンで比較します。比較基準にはAICを使います(参考)。

$\beta_1$ : パラメータなし(切片のみ)モデル
$\beta_1+\beta_2bmi$ : bmiモデル
$\beta_1+\beta_2bmi+\beta_3sex$ : bmi+sexモデル
$\beta_1 + \beta_2ldl+\beta_3age$ : ldl+ageモデル
$\beta_1 + \beta_2ldl+\beta_3age+\beta_4 \times ldl \times age$ : ldl+age+ldl×ageモデル(交互作用項有り)

まずはパラメータなしモデルで、ベースラインとなるAICを確認します。

AICは5041
このモデルの最大尤度は-2519.6 (残差逸脱度は126.8)

# pip install statsmodels
import statsmodels.api as sm
import statsmodels.formula.api as smf

# y~1とすることで、切片のみモデルを作る
model0 = smf.glm('y ~ 1', data=data, family=sm.families.Gamma(link=sm.families.links.log))
result0 = model0.fit(disp=0)

print('AIC:', result0.aic)
# AIC: 5041.186062839457

result0.summary()

f:id:ohke:20180317123522p:plain

次に、bmiモデルで最尤推定させたところ、AICは4884となりました。パラメータ無しのモデルと比較すると、-157改善されています。

model1 = smf.glm('y ~ bmi', data=data, family=sm.families.Gamma(link=sm.families.links.log))
result1 = model1.fit(disp=0)

print('AIC:', result1.aic)
# AIC: 4883.66124171387

もう1つ、bmi+sexモデルでも最尤推定します。AICは4885となり、bmiモデルよりも大きいため、複雑度が増すだけのモデルと言えます。したがって、bmiモデルの方が良い予測を行えるモデルです。

model2 = smf.glm('y ~ bmi+sex', data=data, family=sm.families.Gamma(link=sm.families.links.log))
result2 = model2.fit(disp=0)

print('AIC:', result2.aic)
# AIC: 4885.365291972465

ここまでの3つのモデルの予測値を比較してみます。

パラメータ無し(model0)は、単一の平均値(152)で推定されています
bmiモデル(model1)とbmi+sexモデル(model2)はほぼ同じ予測値となっています

f:id:ohke:20180317130045p:plain

次に、ldl+ageモデルとldl+age+ldl×ageモデルから、交互作用項の有無による差を見ていきます。
一般的に、ldl(LDLコレステロール)は加齢に伴って上がる検査値です(相関係数からも緩い正相関が見られます)。つまり、2つは完全に独立したパラメータではありません。そこで交互作用項ldl×ageを導入します。加齢によるLDLコレステロール値の上昇を打ち消す(つまり負の係数となる)ことを期待しています。

2つモデルのAICの差は5.5で、交互作用項の導入による予測精度の向上が期待できます。

model3 = smf.glm('y ~ ldl+age', data=data, family=sm.families.Gamma(link=sm.families.links.log))
result3 = model3.fit(disp=0)
print('AIC:', result3.aic)
# AIC: 5019.939140159227

model4 = smf.glm('y ~ ldl+age+ldl:age', data=data, family=sm.families.Gamma(link=sm.families.links.log))
result4 = model4.fit(disp=0)
print('AIC:', result4.aic)
# AIC: 5014.413561464431

ldl+age+ldl×ageの各パラメータを見ると、 ldlとageは正の係数(それぞれ1.49、1.64)となっていますが、 ldl×ageは負の係数(-28.92)となっており、期待通りの結果となっていることがわかります。

f:id:ohke:20180318105207p:plain