け日記

最近はPythonでいろいろやってます

Python: statsmodelsでガンマ回帰

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

データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学)

実装はこちらです。

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:plainf:id:ohke:20180317120139p:plain

GLMの設計

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

GLMでは、確率分布・リンク関数・線形予測子を決める必要がありました (1/252/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を使います(参考)。

  1.  \beta_1 : パラメータなし(切片のみ)モデル
  2.  \beta_1+\beta_2bmi : bmiモデル
  3.  \beta_1+\beta_2bmi+\beta_3sex : bmi+sexモデル
  4.  \beta_1 + \beta_2ldl+\beta_3age : ldl+ageモデル
  5.  \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