statsmodelsでNegative Binomial Regressionをかけようと思ったら，想像以上の落とし穴があったのでここに記録を示す．
Negative binomial RegressionはPoisson Regressionの平均と分散が同じ制約を緩めたregressionと捉える．実際に，Possion Regressionのモデルの当てはまりを向上させるためには3つのregressionを考える．Negative Binomial Regressionはその方法の一つ，かつ，0が多いデータに対してもモデルの当てはまりが良い．0が多いデータでも当てられるregressionとしては，Zero inflated poisson regressionが上げられる，しかし，0になる確率pと平均のパラメターλの両方に説明変数を効かせると解釈がしにくい．最後に，過分散を踏まえて標準偏差を計算する方法としてquasi-possion regressionがある．これは点推定の結果は一緒だが，標準偏差が通常のPoisson regressionより広くなる．データに対して頑健な結果を出してくれる．
さて，以上を踏まえてNegative Binomial Regressionをstatsmodelsで実装していく．
Negative Binomial Regressionの実装
実装自体は簡単．
from patsy import dmatrix
import statsmodels.api as sm
data = sm.datasets.get_rdataset("warpbreaks", "datasets").data
X = dmatrix("~ wool + C(tension, Treatment('L'))", data=data, return_type="dataframe")
y = data["breaks"]
model = sm.NegativeBinomial(y, X).fit()
model.summary()
それか，formulaを使うなら，
import statsmodels.formula.api as smf
import statsmodels.api as sm
data = sm.datasets.get_rdataset("warpbreaks", "datasets").data
model = smf.negativebinomial("breaks ~ wool + C(tension, Treatment('L'))", data=data, ).fit()
model.summary()
んで結果が，
NegativeBinomial Regression Results
==============================================================================
Dep. Variable: breaks No. Observations: 54
Model: NegativeBinomial Df Residuals: 50
Method: MLE Df Model: 3
Date: Fri, 10 Jun 2022 Pseudo R-squ.: 0.04391
Time: 07:36:08 Log-Likelihood: -199.38
converged: True LL-Null: -208.54
Covariance Type: nonrobust LLR p-value: 0.0003792
===================================================================================================
coef std err z P>|z| [0.025 0.975]
---------------------------------------------------------------------------------------------------
Intercept 3.6734 0.095 38.516 0.000 3.486 3.860
wool[T.B] -0.1862 0.101 -1.835 0.066 -0.385 0.013
C(tension, Treatment('L'))[T.H] -0.5114 0.124 -4.134 0.000 -0.754 -0.269
C(tension, Treatment('L'))[T.M] -0.2992 0.122 -2.446 0.014 -0.539 -0.059
alpha 0.1006 0.026 3.882 0.000 0.050 0.151
===================================================================================================
となる．alphaがRにおけるθの逆数でdispersionを表す．
今回は，MLEの推定にIterationsが11回で終わっているので良いが，ときどきconvergenceがFailuerするときがある．このようなときは，
model = smf.negativebinomial("breaks ~ wool + C(tension, Treatment('L'))", data=data, ).fit(max_iter=100)
のように，iterationの上限を35から更に上にしておく．
また，Quassi-NegativeBinomial Regressionを行うときは，
model = smf.negativebinomial("breaks ~ wool + C(tension, Treatment('L'))", data=data, ).fit(cov_type="HC0")
のように，cov_typeをHC0, HC1, HC2, HC3のいずれかにする．詳細は，statsmodels.regression.linear_model.RegressionResults に書かれているHeteroscedasticity robust covariance matrix，おそらく，White’s (1980)，MacKinnon and White’s (1985) を参照のこと．
落とし穴
Negative Binomial Regressionを用いるときの落とし穴としては，当初以下のように実装するものだと思っていた．
import statsmodels.api as sm
import statsmodels.formula.api as smf
data = sm.datasets.get_rdataset("warpbreaks", "datasets").data
model = smf.glm("breaks ~ wool + C(tension, Treatment('L'))", data=data, family=sm.families.NegativeBinomial() ).fit()
display(model.summary())
一般的なglmの実装と同じだが，現在，こちらの連続量として取り扱うNegativeBinomialは挙動が怪しくて，RのMASSパッケージのglm.nbと結果が異なってくる．開発者もdiscreteなNegativeBinomialを用いることを推奨している．また，discreteな場合を用いないとdispersion parameterの推定も行ってくれない．ここが一番，documentsなどから読み取れなかったことだった．
その2．さらにNegativeBinomialPを用いることも推奨されたが，これは挙動が不安定．NegativeBinomialPはNegative Binomial Regressionを行う際の平均と分散の関係を指定するが，その次数を指定出来る．しかしながら，使えない．
まだまだ開発をまたなければいけないよう．
ちなみに，MLEの結果を調べたければ，model.mle_settingsかmodel.mle_retvalsの中にある．特にconvergedの結果に関しては，model.mle_retvals[“converged”]で取得可能．v0.14の段階で記法の改善を期待．
参考文献
・FAQ: NegativeBinomial regression’s default uses quassi method? Reason for differences with R’s glm.nb.
・Python Negative Binomial Regression – Results Don’t Match those from R
・BUG/Review GLM scale type negative binomial
・statsmodels.regression.linear_model.RegressionResults
・statsmodels.discrete.discrete_model.NegativeBinomial
・statsmodels.discrete.discrete_model.NegativeBinomialP
・William Greene (2008) Functional forms for the negative binomial model for count data
コメント