This tutorial explains how to quickly implement new maximum likelihood models in statsmodels. 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.
Consider a negative binomial regression model for count data with log-likelihood (type NB-2) function expressed as:
\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)
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:
In [1]: import numpy as np
In [2]: from scipy.stats import nbinom
In [3]: 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
...:
We create a new model class which inherits from GenericLikelihoodModel:
In [4]: from statsmodels.base.model import GenericLikelihoodModel
In [5]: class NBin(GenericLikelihoodModel):
...: def __init__(self, endog, exog, **kwds):
...: super(NBin, self).__init__(endog, exog, **kwds)
...:
In [6]: def nloglikeobs(self, params):
...: alph = params[-1]
...: beta = params[:-1]
...: ll = _ll_nb2(self.endog, self.exog, beta, alph)
...: return -ll
...:
# Reasonable starting values
Two important things to notice:
That’s it! You’re done!
The Medpar 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 columns:
In [7]: import pandas as pd
In [8]: url = 'http://vincentarelbundock.github.com/Rdatasets/csv/COUNT/medpar.csv'
In [9]: medpar = pd.read_csv(url)
In [10]: medpar.head()
Out[10]:
Unnamed: 0 los hmo white died age80 type type1 type2 type3 \
0 1 4 0 1 0 0 1 1 0 0
1 2 9 1 1 0 0 1 1 0 0
2 3 3 1 1 1 1 1 1 0 0
3 4 9 0 1 0 0 1 1 0 0
4 5 1 0 1 1 1 1 1 0 0
provnum
0 30001
1 30001
2 30001
3 30001
4 30001
The model we are interested in has a vector of non-negative integers as dependent variable (los), and 5 regressors: Intercept, type2, type3, hmo, white.
For estimation, we need to create 2 numpy arrays (pandas DataFrame should also work): a 1d array of length N to hold los values, and a N by 5 array to hold our 5 regressors. These arrays can be constructed manually or using any number of helper functions; the details matter little for our current purposes. Here, we build the arrays we need using the Patsy package:
In [11]: import patsy
In [12]: y, X = patsy.dmatrices('los~type2+type3+hmo+white', medpar)
In [13]: print y[:5]
[[ 4.]
[ 9.]
[ 3.]
[ 9.]
[ 1.]]
In [14]: print X[:5]
[[ 1. 0. 0. 0. 1.]
[ 1. 0. 0. 1. 1.]
[ 1. 0. 0. 1. 1.]
[ 1. 0. 0. 0. 1.]
[ 1. 0. 0. 0. 1.]]
Then, we fit the model and extract some information:
In [15]: mod = NBin(y, X)
In [16]: res = mod.fit()
---------------------------------------------------------------------------
AttributeError Traceback (most recent call last)
<ipython-input-16-deef2687e692> in <module>()
----> 1 res = mod.fit()
/build/statsmodels-iqIBhp/statsmodels-0.5.0+git13-g8e07d34/debian/python-statsmodels/usr/lib/python2.7/dist-packages/statsmodels/base/model.pyc in fit(self, start_params, method, maxiter, full_output, disp, callback, retall, **kwargs)
813 method=method, maxiter=maxiter,
814 full_output=full_output,
--> 815 disp=disp, callback=callback, **kwargs)
816 genericmlefit = GenericLikelihoodModelResults(self, mlefit)
817
/build/statsmodels-iqIBhp/statsmodels-0.5.0+git13-g8e07d34/debian/python-statsmodels/usr/lib/python2.7/dist-packages/statsmodels/base/model.pyc in fit(self, start_params, method, maxiter, full_output, disp, fargs, callback, retall, **kwargs)
355 disp=disp, maxiter=maxiter, callback=callback,
356 retall=retall, full_output=full_output,
--> 357 hess=hess)
358
359 if not full_output: # xopt should be None and retvals is argmin
/build/statsmodels-iqIBhp/statsmodels-0.5.0+git13-g8e07d34/debian/python-statsmodels/usr/lib/python2.7/dist-packages/statsmodels/base/model.pyc in _fit_mle_nm(f, score, start_params, fargs, kwargs, disp, maxiter, callback, retall, full_output, hess)
481 ftol=ftol, maxiter=maxiter, maxfun=maxfun,
482 full_output=full_output, disp=disp, retall=retall,
--> 483 callback=callback)
484 if full_output:
485 if not retall:
/usr/lib/python2.7/dist-packages/scipy/optimize/optimize.pyc in fmin(func, x0, args, xtol, ftol, maxiter, maxfun, full_output, disp, retall, callback)
371 'return_all': retall}
372
--> 373 res = _minimize_neldermead(func, x0, args, callback=callback, **opts)
374 if full_output:
375 retlist = res['x'], res['fun'], res['nit'], res['nfev'], res['status']
/usr/lib/python2.7/dist-packages/scipy/optimize/optimize.pyc in _minimize_neldermead(func, x0, args, callback, xtol, ftol, maxiter, maxfev, disp, return_all, **unknown_options)
436 if retall:
437 allvecs = [sim[0]]
--> 438 fsim[0] = func(x0)
439 nonzdelt = 0.05
440 zdelt = 0.00025
/usr/lib/python2.7/dist-packages/scipy/optimize/optimize.pyc in function_wrapper(*wrapper_args)
279 def function_wrapper(*wrapper_args):
280 ncalls[0] += 1
--> 281 return function(*(wrapper_args + args))
282
283 return ncalls, function_wrapper
/build/statsmodels-iqIBhp/statsmodels-0.5.0+git13-g8e07d34/debian/python-statsmodels/usr/lib/python2.7/dist-packages/statsmodels/base/model.pyc in <lambda>(params, *args)
327
328 nobs = self.endog.shape[0]
--> 329 f = lambda params, *args: -self.loglike(params, *args) / nobs
330 score = lambda params: -self.score(params) / nobs
331 try:
/build/statsmodels-iqIBhp/statsmodels-0.5.0+git13-g8e07d34/debian/python-statsmodels/usr/lib/python2.7/dist-packages/statsmodels/base/model.pyc in loglike(self, params)
762
763 def loglike(self, params):
--> 764 return self.loglikeobs(params).sum(0)
765
766 def nloglike(self, params):
/build/statsmodels-iqIBhp/statsmodels-0.5.0+git13-g8e07d34/debian/python-statsmodels/usr/lib/python2.7/dist-packages/statsmodels/base/model.pyc in loglikeobs(self, params)
768
769 def loglikeobs(self, params):
--> 770 return -self.nloglikeobs(params)
771
772 def score(self, params):
AttributeError: 'NBin' object has no attribute 'nloglikeobs'
Extract parameter estimates, standard errors, p-values, AIC, etc.:
In [17]: res.params
Out[17]:
array([-0.0168, 0.0099, -0.0187, -0.0142, 0.2545, 0.2407, 0.0804,
-1.9522, -0.3341, -0.169 , 0.0049, -0.0036, -0.0141, -0.004 ,
-0.0039, 0.0917, 0.049 , 0.008 , 0.0002, -0.0022, 2.9589])
In [18]: res.bse
Out[18]:
array([ 0.0004, 0.0006, 0.0007, 0.0004, 0.0299, 0.0571, 0.0139,
0.3168, 0.0613, 0.0327, 0.0013, 0.0002, 0.0019, 0.0005,
0.001 , 0.0145, 0.0075, 0.0015, 0. , 0.0003, 1.5467])
In [19]: res.pvalues
Out[19]:
array([ 0. , 0. , 0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0.0001, 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. , 0. , 0.0557])
In [20]: res.aic
Out[20]: 6039.2251179879531
As usual, you can obtain a full list of available information by typing dir(res).
To ensure that the above results are sound, we compare them to results obtained using the MASS implementation for R:
url = 'http://vincentarelbundock.github.com/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 and R parameter estimates agree up to the fourth decimal. The standard errors, however, agree only up to the second decimal. This discrepancy may be 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 may be better to use analytical derivatives with the LikelihoodModel class.