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
コメント