最大似然估计（通用模型）¶
本教程介绍如何在 statsmodels 中快速实现新的最大似然模型。我们给出两个示例
用于二元因变量的 Probit 模型
用于计数数据的负二项式模型
The GenericLikelihoodModel 类通过提供诸如自动数值微分和统一接口到 scipy 优化函数之类的工具来简化此过程。使用 statsmodels，用户只需“插入”对数似然函数即可拟合新的 MLE 模型。
示例 1：Probit 模型¶
[1]:
import numpy as np
from scipy import stats
import statsmodels.api as sm
from statsmodels.base.model import GenericLikelihoodModel
The 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
==============================================================================
将您的 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.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 类提供自动微分，因此我们无需提供 Hessian 或 Score 函数来计算协方差估计。
示例 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)
需要注意两点
nloglikeobs：此函数应返回每个观测值（即 endog/X 矩阵的行）的负对数似然函数的一个评估值。start_params：需要提供一个一维起始值数组。此数组的大小决定了优化中将使用的参数数量。
就这样！你完成了！
使用示例¶
The 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。
为了进行估计，我们需要创建两个变量来保存我们的回归变量和结果变量。这些变量可以是 ndarrays 或 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 对负二项式模型的实现来检查结果，该实现使用解析得分函数和 Hessian。
[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
数值精度¶
The statsmodels 通用 MLE 和 R 参数估计值在小数点后第四位一致。然而，标准误仅在小数点后第二位一致。这种差异是我们 Hessian 数值估计中精度不足的结果。在当前上下文中，MASS 和 statsmodels 标准误估计值之间的差异在实质上是无关紧要的，但这突出了这样一个事实，即需要非常精确估计的用户可能并不总是想在使用数值导数时依赖默认设置。在这种情况下，最好使用 LikelihoodModel 类中的解析导数。
最后更新：2024 年 10 月 3 日