[Python] 一般化線形混合モデル(glmm)をGPBoostで実装する

Pythonはよく,Rで当たり前のように出来ることが実装されていないと言われる(実際そう感じる).その状況を徐々に改善していこうとPythonista達が新しいパッケージを作成していってくれている.

今回は,Rで出来るけどPythonでなかなか出来ないシリーズの一つ,glmmの実装に関して説明していく.statsmodelsでは,Generalized Linear Mixed Effects Models は推定そのものがエラーでたり推定値がおかしかったりと,使い物になりません.また,MCMC packagesのPyMC3やemceeで実装出来ますが,MCMCなので当たり前ですがやたら時間がかかります.そこで,glmmをlikelihoodをベースに近似で計算してくれるパッケージが欲しいと切実に思う.Rでは,lm4がガウス・エルミート求積法を用いて行っていますが,C++で書かれたLaGaBoost algorithmを用いた実装がPythonでも使うえるパッケージ,GPBoostが出てきた.このパッケージの妥当性を検証していきたい.

検証

株式会社Nospareさんの一般化線形モデルと一般化線形混合モデル が非常に分かりやすく例として適応しやすかったので,これを参考にする.Rのlm4の結果は,こちらのサイトを参照すること.

*注 : version 0.7.8以上で再現が可能になっている.glmmではLaplace approximationを用いていて,non-linear fixed effectsモデルを取り扱うときは,LaGaBoost algorithmを用いているよう.
https://github.com/fabsig/GPBoost/issues/68

さて,Pythonで同様の内容を実装するには以下のようにする.

import gpboost as gpb
from patsy import dmatrix
import statsmodels.api as sm

data = sm.datasets.get_rdataset("grouseticks", "lme4").data
group = data["LOCATION"]
y = data["TICKS"]
X = dmatrix("C(YEAR)", data=data)
model = gpb.GPModel(group_data=group, likelihood="poisson")
model.fit(y=y, X=X, params={"std_dev":True})
model.summary()

余談だが,statsmodelsのformulaを使えないときは,patsyのdmatrixを用いると,formulaを用いてデータを生成することが可能.つまるところ,patsyのdmatrixを使えば,pd.get_dummiesとかちまちま使う必要が消える.

さて,僕のPCではすぐに結果が出る.

[GPBoost] ===================================================
Model summary:
 Log-lik     AIC    BIC
-1149.15 2306.31 2322.3
Nb. observations: 403
Nb. groups: 63 (LOCATION)
---------------------------------------------------
Covariance parameters (random effects):
          Param.
LOCATION  1.6367
---------------------------------------------------
Linear regression coefficients (fixed effects):
               Param.  Std. dev.  z value  P(>|z|)
Intercept      0.6730     0.1822   3.6944   0.0002
C(YEAR)[T.96]  0.9438     0.0852  11.0781   0.0000
C(YEAR)[T.97] -1.4351     0.1233 -11.6347   0.0000
===================================================

さらに,パラメーターを引き出したいときは,

model.get_coef()
model.get_cov_pars()

一瞬で計算出来るglmmパッケージがPythonにも出来たのはとても喜ばしいと考える.

———-雑感(`・ω・´)———-
ちなみに僕は一般的なBayesian modelをMCMCで実装しなければいけないときはPyMC3を.カスタムなlikelihoodをMCMCで推定するときはemceeを使うという風に使い分けている.

丁度,アウトプットの見た目がよくなったcommitがこの部分のbasic.py

https://github.com/fabsig/GPBoost/commit/28492b0be84a34e13292661722cdd8844139cd29#diff-c8a283f07121de87b9e1aa922970e8aa586f3dc3819affd50c5db42dfee69ec2

コメント

タイトルとURLをコピーしました