最大概似估計（通用模型）¶
本教學說明如何在 statsmodels 中快速實作新的最大概似模型。我們將提供兩個範例：
二元應變數的 Probit 模型
計數資料的負二項式模型
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
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：需要提供一維起始值陣列。此陣列的大小決定了將在最佳化中使用的參數數量。
就是這樣！您完成了！
使用範例¶
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。
對於估計，我們需要建立兩個變數來保存迴歸變數和結果變數。這些可以是 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 實作的負二項式模型來檢查結果，該模型使用解析分數函數和 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
數值精度¶
statsmodels 通用 MLE 和 R 參數估計值在小數點後第四位之前一致。但是，標準差僅在小數點後第二位之前一致。此差異是我們的 Hessian 數值估計值不精確所造成的。在目前的情況下，MASS 和 statsmodels 標準差估計值之間的差異在實質上無關緊要，但這突顯了一個事實，即需要非常精確估計值的使用者在使用數值導數時可能不總是希望依賴預設設定。在這種情況下，最好使用 LikelihoodModel 類別的分析導數。
上次更新：2024 年 10 月 03 日