Pythonで負の二項回帰を行う方法
まず、必要なパッケージをすべてインポートします。
import pandas as pd
from patsy import dmatrices
import numpy as np
import statsmodels.api as sm
import matplotlib.pyplot as plt
次に、カウントデータセットのpandas DataFrameを作成します。
df = pd.read_csv('nyc_bb_bicyclist_counts.csv', header=0, infer_datetime_format=True, parse_dates=, index_col=)
X行列にいくつかの派生回帰変数を追加していきます。
ds = df.index.to_series()
df = ds.dt.month
df = ds.dt.dayofweek
df = ds.dt.day
Date変数は絶対的な日付値を含むので、リグレッサーとしては使用しませんが、Dateはpandas DataFrameのインデックスとしてすでに消費されているので、ドロップするために特別なことをする必要がないのです。
訓練データとテストデータセットを作成しましょう。
mask = np.random.rand(len(df)) < 0.8
df_train = df
df_test = df
print('Training data set length='+str(len(df_train)))
print('Testing data set length='+str(len(df_test)))
STEP 1: これから訓練データセットにポアソン回帰モデルを設定し、フィットさせます。
Patsy記法の回帰式をセットアップします。 BB_COUNTが従属変数であり、回帰変数に依存することをpatsyに伝えている。 DAY, DAY_OF_WEEK, MONTH, HIGH_T, LOW_T and PRECIP.
expr = """BB_COUNT ~ DAY + DAY_OF_WEEK + MONTH + HIGH_T + LOW_T + PRECIP"""
トレーニングとテストデータセットのXとYマトリックスをセットアップします。
y_train, X_train = dmatrices(expr, df_train, return_type='dataframe')
y_test, X_test = dmatrices(expr, df_test, return_type='dataframe')
statsmodelsのGLMクラスを使って、学習データセットでポアソン回帰モデルを学習します。
poisson_training_results = sm.GLM(y_train, X_train, family=sm.families.Poisson()).fit()
これでポアソン回帰モデルの学習は終了です。 学習の結果を見るには、学習サマリーをプリントアウトしてください。
print(poisson_training_results.summary())
これは次のように出力されます：
我々の本当の関心は、学習によって作られた適合率 λ のベクトルにあるのですが、これは、学習した結果ではありません。 このレートベクトルはパラメータpoisson_training_results.muに含まれています。
print(poisson_training_results.mu)
print(len(poisson_training_results.mu))
次の出力は、フィットしたλベクトルの最初の数値と最後の数値です：
これでSTEP1：ポアソン回帰モデルのフィットが完了しました。
ステップ2：データセットに補助的なOLS回帰モデルをあてはめ、あてはめたモデルを使ってαの値を求めます。
apiパッケージをインポートします。
import statsmodels.formula.api as smf
学習データセットのデータフレームに「BB_LAMBDA」という新しい列として、λベクトルを追加してください。 λの次元は(n×1)であることを思い出してください。 この例では、(161 x 1)となります。 また、λベクトルはpoisson_training_results.mu :
df_train = poisson_training_results.mu
で利用可能であることを思い出してください。 この新しい列は、OLS回帰の従属変数の値を格納します。 以下のOLS回帰式の左辺になります:
df_train = df_train.apply(lambda x: ((x - x)**2 - x) / x, axis=1)
上記のコードスニペットで、太字の部分は上記の補助OLSR式の左辺にあたります。
patsyを使ってOLSRのモデル指定を形成してみましょう。 AUX_OLS_DEPが従属変数で、BB_LAMBDA(これはレートベクトルλ)で説明されることをpatsyに伝えたいと思います。 式の最後の’-1’はpatsyのシンタックスで、回帰の切片を使用しないこと、つまり、CameronとTrivediが提案したように、原点を通る直線をフィットさせることを意味しています。
OLSRモデルの設定と適合：
aux_olsr_results = smf.ols(ols_expr, df_train).fit()
回帰パラメータをプリントアウト：
print(aux_olsr_results.params)
単一の回帰変数BB_LAMBDAに対応して以下の単一の係数がプリントアウトされることがわかるでしょう。 この係数は私たちが求めていたαです。
BB_LAMBDA 0.037343
Is α is statistically significant?
ここで私たちは、非常に重要な質問に答える必要があります。 このαの値（0.037343）は統計的に有意なのでしょうか？
なぜこれを調べることがそんなに重要なのでしょうか。 もしαがゼロなら、次の式:
・・・分散＝平均に還元されることを思い出してみてください。 これはポアソン回帰モデルの分散関数です。
αの値が統計的に有意でない場合、負の二項回帰モデルはポアソン回帰モデルよりも学習データセットにフィットする良い仕事ができません。
OLSResultsオブジェクトは回帰係数αのtスコアを含んでいます。
aux_olsr_results.tvalues
これはプリントアウトします：
BB_LAMBDA 4.814096
t値計算機から、99%信頼水準（右側）で、自由度=(172 observations) -( 6 regression variables) – (1 dispersion parameter α)=165 の臨界t値は2.34916であることがわかります。 これはαのt統計量である4.814096より余裕で小さい。
α=0.037343 は統計的に有意であると結論づけられる。
これでSTEP2：αの決定が完了する。
STEP3：STEP2 で求めたαの値をstatsmodels.genmod.families.family.NegativeBinomialクラスに供給し、学習データセットでNB2モデルを学習させる。
これはstatsmodelsのワンステップ操作です:
nb2_training_results = sm.GLM(y_train, X_train,family=sm.families.NegativeBinomial(alpha=aux_olsr_results.params)).fit()
前回と同様に、学習サマリーを出力します:
print(nb2_training_results.summary())
以下のサマリーが出力されます。
ステップ4：トレーニングしたNB2モデルを使って予測をしてみよう。
予測はstatsmodelsでもシングルステップの手順です：
nb2_predictions = nb2_training_results.get_prediction(X_test)
予測を出力してみましょう：
predictions_summary_frame = nb2_predictions.summary_frame()
print(predictions_summary_frame)
以下は出力結果の最初の数行です。
また、テストデータの実際のカウントに対して予測されたカウントをプロットしてみましょう。
predicted_counts=predictions_summary_frameactual_counts = y_testfig = plt.figure()fig.suptitle('Predicted versus actual bicyclist counts on the Brooklyn bridge')predicted, = plt.plot(X_test.index, predicted_counts, 'go-', label='Predicted counts')actual, = plt.plot(X_test.index, actual_counts, 'ro-', label='Actual counts')plt.legend(handles=)plt.show()
出力は以下のとおりです：
悪くありませんね！
このように、NB2モデルによる予測は、自転車に乗っている人の数と、実際の自転車乗車の数を比較することができます。 NB2モデルは自転車数のトレンドにほぼ追従しているようです。 また、ポアソン回帰モデルのパフォーマンスと同様に、その予測は実際の値から大きく外れているケースもあります
。