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