最尤推定（汎用モデル）¶
このチュートリアルでは、`statsmodels` で新しい最尤モデルを迅速に実装する方法について説明します。2つの例を示します。
二値従属変数のためのプロビットモデル
カウントデータのための負の二項モデル
`GenericLikelihoodModel` クラスは、自動数値微分や `scipy` 最適化関数への統一インターフェースなどのツールを提供することにより、プロセスを容易にします。 `statsmodels` を使用すると、ユーザーは対数尤度関数を「プラグイン」するだけで、新しいMLEモデルを適合させることができます。
例1：プロビットモデル¶
[1]:
import numpy as np
from scipy import stats
import statsmodels.api as sm
from statsmodels.base.model import GenericLikelihoodModel
`Spector` データセットは `statsmodels` と 함께 배포됩니다. 従属変数（`endog`）の値のベクトルと回帰変数の行列（`exog`）には、次のようにアクセスできます。
[2]:
data = sm.datasets.spector.load_pandas()
exog = data.exog
endog = data.endog
print(sm.datasets.spector.NOTE)
print(data.exog.head())
::
Number of Observations - 32
Number of Variables - 4
Variable name definitions::
Grade - binary variable indicating whether or not a student's grade
improved. 1 indicates an improvement.
TUCE - Test score on economics test
PSI - participation in program
GPA - Student's grade point average
GPA TUCE PSI
0 2.66 20.0 0.0
1 2.89 22.0 0.0
2 3.28 24.0 0.0
3 2.92 12.0 0.0
4 4.00 21.0 0.0
次に、回帰変数の行列に定数項を追加します。
[3]:
exog = sm.add_constant(exog, prepend=True)
独自の尤度モデルを作成するには、`loglike` メソッドをオーバーライドするだけです。
[4]:
class MyProbit(GenericLikelihoodModel):
def loglike(self, params):
exog = self.exog
endog = self.endog
q = 2 * endog - 1
return stats.norm.logcdf(q*np.dot(exog, params)).sum()
モデルを推定し、サマリーを出力します。
[5]:
sm_probit_manual = MyProbit(endog, exog).fit()
print(sm_probit_manual.summary())
Optimization terminated successfully.
Current function value: 0.400588
Iterations: 292
Function evaluations: 494
MyProbit Results
==============================================================================
Dep. Variable: GRADE Log-Likelihood: -12.819
Model: MyProbit AIC: 33.64
Method: Maximum Likelihood BIC: 39.50
Date: Thu, 03 Oct 2024
Time: 15:50:26
No. Observations: 32
Df Residuals: 28
Df Model: 3
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
const -7.4523 2.542 -2.931 0.003 -12.435 -2.469
GPA 1.6258 0.694 2.343 0.019 0.266 2.986
TUCE 0.0517 0.084 0.617 0.537 -0.113 0.216
PSI 1.4263 0.595 2.397 0.017 0.260 2.593
==============================================================================
独自に実装したプロビットモデルと `statsmodels` の「既製」の実装を比較します。
[6]:
sm_probit_canned = sm.Probit(endog, exog).fit()
Optimization terminated successfully.
Current function value: 0.400588
Iterations 6
[7]:
print(sm_probit_canned.params)
print(sm_probit_manual.params)
const -7.452320
GPA 1.625810
TUCE 0.051729
PSI 1.426332
dtype: float64
[-7.45233176 1.62580888 0.05172971 1.42631954]
[8]:
print(sm_probit_canned.cov_params())
print(sm_probit_manual.cov_params())
const GPA TUCE PSI
const 6.464166 -1.169668 -0.101173 -0.594792
GPA -1.169668 0.481473 -0.018914 0.105439
TUCE -0.101173 -0.018914 0.007038 0.002472
PSI -0.594792 0.105439 0.002472 0.354070
[[ 6.46416776e+00 -1.16966614e+00 -1.01173187e-01 -5.94788999e-01]
[-1.16966614e+00 4.81472101e-01 -1.89134577e-02 1.05438217e-01]
[-1.01173187e-01 -1.89134577e-02 7.03758407e-03 2.47189354e-03]
[-5.94788999e-01 1.05438217e-01 2.47189354e-03 3.54069513e-01]]
`GenericMaximumLikelihood` クラスは自動微分を提供するため、共分散推定値を計算するためにヘッセ行列やスコア関数を指定する必要はありませんでした。
例2：カウントデータの負の二項回帰¶
対数尤度（タイプNB-2）関数が次のように表されるカウントデータの負の二項回帰モデルを考えます。
\[\mathcal{L}(\beta_j; y, \alpha) = \sum_{i=1}^n y_i ln \left ( \frac{\alpha exp(X_i'\beta)}{1+\alpha exp(X_i'\beta)} \right ) - \frac{1}{\alpha} ln(1+\alpha exp(X_i'\beta)) + ln \Gamma (y_i + 1/\alpha) - ln \Gamma (y_i+1) - ln \Gamma (1/\alpha)\]
ここで、\(X\) は回帰変数の行列、\(\beta\) は係数のベクトル、\(\alpha\) は負の二項不均一性パラメータです。
`scipy` の `nbinom` 分布を使用すると、この尤度は次のように簡単に記述できます。
[9]:
import numpy as np
from scipy.stats import nbinom
[10]:
def _ll_nb2(y, X, beta, alph):
mu = np.exp(np.dot(X, beta))
size = 1/alph
prob = size/(size+mu)
ll = nbinom.logpmf(y, size, prob)
return ll
新しいモデルクラス¶
`GenericLikelihoodModel` を継承する新しいモデルクラスを作成します。
[11]:
from statsmodels.base.model import GenericLikelihoodModel
[12]:
class NBin(GenericLikelihoodModel):
def __init__(self, endog, exog, **kwds):
super(NBin, self).__init__(endog, exog, **kwds)
def nloglikeobs(self, params):
alph = params[-1]
beta = params[:-1]
ll = _ll_nb2(self.endog, self.exog, beta, alph)
return -ll
def fit(self, start_params=None, maxiter=10000, maxfun=5000, **kwds):
# we have one additional parameter and we need to add it for summary
self.exog_names.append('alpha')
if start_params == None:
# Reasonable starting values
start_params = np.append(np.zeros(self.exog.shape[1]), .5)
# intercept
start_params[-2] = np.log(self.endog.mean())
return super(NBin, self).fit(start_params=start_params,
maxiter=maxiter, maxfun=maxfun,
**kwds)
注意すべき重要な点が2つあります。
`nloglikeobs`：この関数は、データセット内の観測値ごとに（つまり、endog / X行列の行ごとに）負の対数尤度関数の1つの評価値を返す必要があります。
`start_params`：1次元配列の開始値を指定する必要があります。この配列のサイズによって、最適化で使用されるパラメータの数が決まります。
これで完了です！
使用例¶
Medpar データセットは、Rdatasets リポジトリ でCSV形式でホストされています。 Pandas ライブラリ の `read_csv` 関数を使用して、データをメモリにロードします。次に、最初の数列を表示します。
[13]:
import statsmodels.api as sm
[14]:
medpar = sm.datasets.get_rdataset("medpar", "COUNT", cache=True).data
medpar.head()
[14]:
| los | hmo | white | died | age80 | type | type1 | type2 | type3 | provnum | |
|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 4 | 0 | 1 | 0 | 0 | 1 | 1 | 0 | 0 | 30001 |
| 1 | 9 | 1 | 1 | 0 | 0 | 1 | 1 | 0 | 0 | 30001 |
| 2 | 3 | 1 | 1 | 1 | 1 | 1 | 1 | 0 | 0 | 30001 |
| 3 | 9 | 0 | 1 | 0 | 0 | 1 | 1 | 0 | 0 | 30001 |
| 4 | 1 | 0 | 1 | 1 | 1 | 1 | 1 | 0 | 0 | 30001 |
対象のモデルは、従属変数（`los`）として非負整数のベクトルを持ち、5つの回帰変数：`Intercept`、`type2`、`type3`、`hmo`、`white` を持ちます。
推定のため、回帰変数と結果変数を保持する2つの変数を作成する必要があります。これらはndarrayまたはpandasオブジェクトです。
[15]:
y = medpar.los
X = medpar[["type2", "type3", "hmo", "white"]].copy()
X["constant"] = 1
次に、モデルを適合させ、いくつかの情報を抽出します。
[16]:
mod = NBin(y, X)
res = mod.fit()
Optimization terminated successfully.
Current function value: 3.209014
Iterations: 805
Function evaluations: 1238
/opt/hostedtoolcache/Python/3.10.15/x64/lib/python3.10/site-packages/statsmodels/base/model.py:2748: UserWarning: df_model + k_constant + k_extra differs from k_params
warnings.warn("df_model + k_constant + k_extra "
/opt/hostedtoolcache/Python/3.10.15/x64/lib/python3.10/site-packages/statsmodels/base/model.py:2752: UserWarning: df_resid differs from nobs - k_params
warnings.warn("df_resid differs from nobs - k_params")
パラメータ推定値、標準誤差、p値、AICなどを抽出します。
[17]:
print('Parameters: ', res.params)
print('Standard errors: ', res.bse)
print('P-values: ', res.pvalues)
print('AIC: ', res.aic)
Parameters: [ 0.2212642 0.70613942 -0.06798155 -0.12903932 2.31026565 0.44575147]
Standard errors: [0.05059259 0.07613047 0.05326096 0.0685414 0.06794696 0.01981542]
P-values: [1.22298084e-005 1.76979047e-020 2.01819053e-001 5.97481232e-002
2.15207253e-253 4.62688811e-112]
AIC: 9604.95320583016
通常どおり、`dir(res)` と入力することで、利用可能な情報の完全なリストを取得できます。推定結果のサマリーを見ることもできます。
[18]:
print(res.summary())
NBin Results
==============================================================================
Dep. Variable: los Log-Likelihood: -4797.5
Model: NBin AIC: 9605.
Method: Maximum Likelihood BIC: 9632.
Date: Thu, 03 Oct 2024
Time: 15:50:28
No. Observations: 1495
Df Residuals: 1490
Df Model: 4
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
type2 0.2213 0.051 4.373 0.000 0.122 0.320
type3 0.7061 0.076 9.275 0.000 0.557 0.855
hmo -0.0680 0.053 -1.276 0.202 -0.172 0.036
white -0.1290 0.069 -1.883 0.060 -0.263 0.005
constant 2.3103 0.068 34.001 0.000 2.177 2.443
alpha 0.4458 0.020 22.495 0.000 0.407 0.485
==============================================================================
テスト¶
解析スコア関数とヘッセ行列を使用するstatsmodelsの負の二項モデルの実装を使用して、結果を確認できます。
[19]:
res_nbin = sm.NegativeBinomial(y, X).fit(disp=0)
print(res_nbin.summary())
NegativeBinomial Regression Results
==============================================================================
Dep. Variable: los No. Observations: 1495
Model: NegativeBinomial Df Residuals: 1490
Method: MLE Df Model: 4
Date: Thu, 03 Oct 2024 Pseudo R-squ.: 0.01215
Time: 15:50:28 Log-Likelihood: -4797.5
converged: True LL-Null: -4856.5
Covariance Type: nonrobust LLR p-value: 1.404e-24
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
type2 0.2212 0.051 4.373 0.000 0.122 0.320
type3 0.7062 0.076 9.276 0.000 0.557 0.855
hmo -0.0680 0.053 -1.276 0.202 -0.172 0.036
white -0.1291 0.069 -1.883 0.060 -0.263 0.005
constant 2.3103 0.068 34.001 0.000 2.177 2.443
alpha 0.4457 0.020 22.495 0.000 0.407 0.485
==============================================================================
[20]:
print(res_nbin.params)
type2 0.221218
type3 0.706173
hmo -0.067987
white -0.129053
constant 2.310279
alpha 0.445748
dtype: float64
[21]:
print(res_nbin.bse)
type2 0.050592
type3 0.076131
hmo 0.053261
white 0.068541
constant 0.067947
alpha 0.019815
dtype: float64
または、RのMASS実装を使用して得られた結果と比較することもできます。
url = 'https://raw.githubusercontent.com/vincentarelbundock/Rdatasets/csv/COUNT/medpar.csv'
medpar = read.csv(url)
f = los~factor(type)+hmo+white
library(MASS)
mod = glm.nb(f, medpar)
coef(summary(mod))
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.31027893 0.06744676 34.253370 3.885556e-257
factor(type)2 0.22124898 0.05045746 4.384861 1.160597e-05
factor(type)3 0.70615882 0.07599849 9.291748 1.517751e-20
hmo -0.06795522 0.05321375 -1.277024 2.015939e-01
white -0.12906544 0.06836272 -1.887951 5.903257e-02
数値精度¶
`statsmodels` の汎用MLEと `R` のパラメータ推定値は、小数点以下4桁まで一致します。ただし、標準誤差は小数点以下2桁までしか一致しません。この食い違いは、ヘッセ行列の数値推定値の不正確さに起因します。現在のコンテキストでは、`MASS` と `statsmodels` の標準誤差推定値の差は実質的に無関係ですが、非常に正確な推定値を必要とするユーザーは、数値微分を使用する場合にデフォルト設定に依存したくない場合があるという事実を浮き彫りにしています。このような場合は、`LikelihoodModel` クラスで解析微分を使用することをお勧めします。
最終更新日：2024年10月3日