最尤推定 (一般的なモデル)¶
このチュートリアルでは、statsmodels で新しい最尤推定（MLE）モデルを迅速に実装する方法を説明します。以下の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: Wed, 28 Aug 2024
Time: 10:35:19
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
==============================================================================
Probitの実装を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.46416777e+00 -1.16966613e+00 -1.01173188e-01 -5.94789000e-01]
[-1.16966613e+00 4.81472102e-01 -1.89134582e-02 1.05438217e-01]
[-1.01173188e-01 -1.89134582e-02 7.03758418e-03 2.47189359e-03]
[-5.94789000e-01 1.05438217e-01 2.47189359e-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データセットは、CSV形式でRdatasetsリポジトリにホストされています。データは、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
C:\Users\user\projects\statsmodels\venv\lib\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 "
C:\Users\user\projects\statsmodels\venv\lib\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.06854139 0.06794694 0.01981542]
P-values: [1.22297865e-005 1.76978911e-020 2.01819053e-001 5.97480571e-002
2.15131189e-253 4.62688465e-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: Wed, 28 Aug 2024
Time: 10:35:20
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: Wed, 28 Aug 2024 Pseudo R-squ.: 0.01215
Time: 10:35:20 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桁までしか一致しません。この不一致は、Hessianの数値推定における不正確さに起因しています。現在の文脈では、MASSとstatsmodelsの標準誤差推定の違いは実質的には無関係ですが、非常に精密な推定が必要なユーザーは、数値導関数を使用する際にデフォルト設定に頼らない方が良いことを示しています。そのような場合には、LikelihoodModelクラスの解析的導関数を使用する方が望ましいです。
最終更新日:
2025年01月28日