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