[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.7以上で再現が可能になっています.
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, "optimizer_cov": "nelder_mead", 
                            "delta_rel_conv": 1e-12})
model.summary()

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

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

[GPBoost] [Warning] Standard deviations of linear regression coefficients for non-Gaussian data can be "very approximative". 
===================================================
Model summary:
    Log-lik        AIC         BIC
-1149.14448 2306.28896 2322.284706
Nb. observations: 403
Nb. groups:  LOCATION
       63
---------------------------------------------------
Covariance parameters (random effects):
          Param.
LOCATION  1.6464
---------------------------------------------------
Linear regression coefficients (fixed effects):
             Param.  Std. dev.  z value  P(>|z|)
Covariate_1  0.6488     0.1827   3.5519   0.0004
Covariate_2  0.9437     0.0852  11.0758   0.0000
Covariate_3 -1.4365     0.1234 -11.6391   0.0000
===================================================

少し結果が見にくいので,design_infoを用いて,

coef = model.get_coef()
coef.columns = X.design_info.column_names
coef

と置き換えておくと,coefに元々指定したパラメーターの内容で結果を見ることが可能だ.

今回の例は,likelihoodがかなり平坦なために計算精度を上げないとRの結果と同じものを得られなかった.が,一瞬で計算出来るglmmパッケージがPythonにも出来たのはとても喜ばしいと考える.

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

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

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

コメント

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