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