Maximum Likelihood Estimation (Generic models)¶
This tutorial explains how to quickly implement new maximum likelihood models in statsmodels
. We give two examples:
Probit model for binary dependent variables
Negative binomial model for count data
The GenericLikelihoodModel
class eases the process by providing tools such as automatic numeric differentiation and a unified interface to scipy
optimization functions. Using statsmodels
, users can fit new MLE models simply by “plugging-in” a log-likelihood function.
Example 1: Probit model¶
[1]:
import numpy as np
from scipy import stats
import statsmodels.api as sm
from statsmodels.base.model import GenericLikelihoodModel
The Spector
dataset is distributed with statsmodels
. You can access a vector of values for the dependent variable (endog
) and a matrix of regressors (exog
) like this:
[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
Them, we add a constant to the matrix of regressors:
[3]:
exog = sm.add_constant(exog, prepend=True)
To create your own Likelihood Model, you simply need to overwrite the loglike method.
[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()
Estimate the model and print a summary:
[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: Tue, 14 Feb 2023
Time: 22:28:59
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
==============================================================================
Compare your Probit implementation to statsmodels
’ “canned” implementation:
[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]]
Notice that the GenericMaximumLikelihood
class provides automatic differentiation, so we did not have to provide Hessian or Score functions in order to calculate the covariance estimates.
Example 2: Negative Binomial Regression for Count Data¶
Consider a negative binomial regression model for count data with log-likelihood (type NB-2) function expressed as:
with a matrix of regressors \(X\), a vector of coefficients \(\beta\), and the negative binomial heterogeneity parameter \(\alpha\).
Using the nbinom
distribution from scipy
, we can write this likelihood simply as:
[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
New Model Class¶
We create a new model class which inherits from 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)
Two important things to notice:
nloglikeobs
: This function should return one evaluation of the negative log-likelihood function per observation in your dataset (i.e. rows of the endog/X matrix).start_params
: A one-dimensional array of starting values needs to be provided. The size of this array determines the number of parameters that will be used in optimization.
That’s it! You’re done!
Usage Example¶
The epilepsy dataset is hosted in CSV format at the Rdatasets repository. We use the read_csv
function from the Pandas library to load the data in memory. We then print the first few entries:
[13]:
import statsmodels.api as sm
[14]:
epilepsy = sm.datasets.get_rdataset("epilepsy", "robustbase", cache=True).data
epilepsy.head()
[14]:
ID | Y1 | Y2 | Y3 | Y4 | Base | Age | Trt | Ysum | Age10 | Base4 | |
---|---|---|---|---|---|---|---|---|---|---|---|
0 | 104 | 5 | 3 | 3 | 3 | 11 | 31 | placebo | 14 | 3.1 | 2.75 |
1 | 106 | 3 | 5 | 3 | 3 | 11 | 30 | placebo | 14 | 3.0 | 2.75 |
2 | 107 | 2 | 4 | 0 | 5 | 6 | 25 | placebo | 11 | 2.5 | 1.50 |
3 | 114 | 4 | 4 | 1 | 4 | 8 | 36 | placebo | 13 | 3.6 | 2.00 |
4 | 116 | 7 | 18 | 9 | 21 | 66 | 22 | placebo | 55 | 2.2 | 16.50 |
The model we are interested in has a vector of non-negative integers as dependent variable (Ysum
), and 3 regressors: Intercept
, Base
, Trt
.
For estimation, we need to create two variables to hold our regressors and the outcome variable. These can be ndarrays or pandas objects.
[15]:
y = epilepsy.Ysum
epilepsy["Trtnum"]=epilepsy["Trt"].map({"placebo": 0, "progabide": 1})
X = epilepsy[["Base", "Trtnum"]].copy()
X["constant"] = 1
Then, we fit the model and extract some information:
[16]:
mod = NBin(y, X)
res = mod.fit()
Optimization terminated successfully.
Current function value: 3.967477
Iterations: 249
Function evaluations: 412
/usr/lib/python3/dist-packages/statsmodels/base/model.py:2694: UserWarning: df_model + k_constant differs from nparams
warnings.warn("df_model + k_constant differs from nparams")
/usr/lib/python3/dist-packages/statsmodels/base/model.py:2696: UserWarning: df_resid differs from nobs - nparams
warnings.warn("df_resid differs from nobs - nparams")
Extract parameter estimates, standard errors, p-values, AIC, etc.:
[17]:
print('Parameters: ', res.params)
print('Standard errors: ', res.bse)
print('P-values: ', res.pvalues)
print('AIC: ', res.aic)
Parameters: [ 0.02733076 -0.21885176 2.4392306 0.30877386]
Standard errors: [0.00338564 0.1553704 0.1573023 0.06454799]
P-values: [6.88375123e-16 1.58959157e-01 3.12810776e-54 1.72155318e-06]
AIC: 474.1623351041066
As usual, you can obtain a full list of available information by typing dir(res)
. We can also look at the summary of the estimation results.
[18]:
print(res.summary())
NBin Results
==============================================================================
Dep. Variable: Ysum Log-Likelihood: -234.08
Model: NBin AIC: 474.2
Method: Maximum Likelihood BIC: 480.4
Date: Tue, 14 Feb 2023
Time: 22:28:59
No. Observations: 59
Df Residuals: 56
Df Model: 2
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
Base 0.0273 0.003 8.073 0.000 0.021 0.034
Trtnum -0.2189 0.155 -1.409 0.159 -0.523 0.086
constant 2.4392 0.157 15.507 0.000 2.131 2.748
alpha 0.3088 0.065 4.784 0.000 0.182 0.435
==============================================================================
Testing¶
We can check the results by using the statsmodels implementation of the Negative Binomial model, which uses the analytic score function and Hessian.
[19]:
res_nbin = sm.NegativeBinomial(y, X).fit(disp=0)
print(res_nbin.summary())
NegativeBinomial Regression Results
==============================================================================
Dep. Variable: Ysum No. Observations: 59
Model: NegativeBinomial Df Residuals: 56
Method: MLE Df Model: 2
Date: Tue, 14 Feb 2023 Pseudo R-squ.: 0.1203
Time: 22:28:59 Log-Likelihood: -234.08
converged: True LL-Null: -266.10
Covariance Type: nonrobust LLR p-value: 1.247e-14
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
Base 0.0273 0.003 8.073 0.000 0.021 0.034
Trtnum -0.2188 0.155 -1.408 0.159 -0.523 0.086
constant 2.4392 0.157 15.507 0.000 2.131 2.748
alpha 0.3088 0.065 4.784 0.000 0.182 0.435
==============================================================================
[20]:
print(res_nbin.params)
Base 0.027330
Trtnum -0.218803
constant 2.439247
alpha 0.308767
dtype: float64
[21]:
print(res_nbin.bse)
Base 0.003386
Trtnum 0.155369
constant 0.157299
alpha 0.064546
dtype: float64
Or we could compare them to results obtained using the MASS implementation for R:
[22]:
%load_ext rpy2.ipython
[23]:
%R f = Ysum~factor(Trt)+Base
%R data(epilepsy,package='robustbase')
%R library(MASS)
%R mod = glm.nb(f, epilepsy)
%R print(coef(summary(mod)))
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.43924704 0.144237213 16.911357 3.710627e-64
factor(Trt)progabide -0.21880303 0.155745198 -1.404878 1.600575e-01
Base 0.02732961 0.002807957 9.732917 2.182424e-22
[23]:
array([[ 2.43924704e+00, 1.44237213e-01, 1.69113572e+01,
3.71062703e-64],
[-2.18803032e-01, 1.55745198e-01, -1.40487820e+00,
1.60057501e-01],
[ 2.73296140e-02, 2.80795707e-03, 9.73291730e+00,
2.18242392e-22]])
Numerical precision¶
The statsmodels
generic MLE and R
parameter estimates agree up to the fourth decimal. The standard errors, however, agree only up to the first decimal. This discrepancy is the result of imprecision in our Hessian numerical estimates. In the current context, the difference between MASS
and statsmodels
standard error estimates is substantively irrelevant, but it highlights the fact that users who need very precise estimates may not always want to rely on default settings when using
numerical derivatives. In such cases, it is better to use analytical derivatives with the LikelihoodModel
class.