`_ .. raw:: html

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¶

In [ ]:

from __future__ import print_function
   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:

In [ ]:

data = sm.datasets.spector.load_pandas()
   exog = data.exog
   endog = data.endog
   print(sm.datasets.spector.NOTE)
   print(data.exog.head())

Them, we add a constant to the matrix of regressors:

In [ ]:

exog = sm.add_constant(exog, prepend=True)

::
   
       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
   1  2.89    22    0
   2  3.28    24    0
   3  2.92    12    0
   4  4.00    21    0

To create your own Likelihood Model, you simply need to overwrite the loglike method.

In [ ]:

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:

In [ ]:

sm_probit_manual = MyProbit(endog, exog).fit()
   print(sm_probit_manual.summary())

Compare your Probit implementation to statsmodels' "canned" implementation:

In [ ]:

sm_probit_canned = sm.Probit(endog, exog).fit()

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:                Mon, 20 Jul 2015                                         
   Time:                        17:43:29                                         
   No. Observations:                  32                                         
   Df Residuals:                      28                                         
   Df Model:                           3                                         
   ==============================================================================
                    coef    std err          z      P>|z|      [95.0% Conf. Int.]
   ------------------------------------------------------------------------------
   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
   ==============================================================================

In [ ]:

print(sm_probit_canned.params)
   print(sm_probit_manual.params)

Optimization terminated successfully.
            Current function value: 0.400588
            Iterations 6

In [ ]:

print(sm_probit_canned.cov_params())
   print(sm_probit_manual.cov_params())

const   -7.452320
   GPA      1.625810
   TUCE     0.051729
   PSI      1.426332
   dtype: float64
   [-7.45233176  1.62580888  0.05172971  1.42631954]

Notice that the GenericMaximumLikelihood class provides automatic differentiation, so we didn't 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:

$$ \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 [ ]:

import numpy as np
   from scipy.stats import nbinom

          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.46416770e+00  -1.16966617e+00  -1.01173180e-01  -5.94788993e-01]
    [ -1.16966617e+00   4.81472116e-01  -1.89134586e-02   1.05438227e-01]
    [ -1.01173180e-01  -1.89134586e-02   7.03758392e-03   2.47189191e-03]
    [ -5.94788993e-01   1.05438227e-01   2.47189191e-03   3.54069512e-01]]

In [ ]:

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:

In [ ]:

from statsmodels.base.model import GenericLikelihoodModel

In [ ]:

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 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 [ ]:

import statsmodels.api as sm

In [ ]:

medpar = sm.datasets.get_rdataset("medpar", "COUNT", cache=True).data
   
   medpar.head()

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 two variables to hold our regressors and the outcome variable. These can be ndarrays or pandas objects.

In [ ]:

y = medpar.los
   X = medpar[["type2", "type3", "hmo", "white"]]
   X["constant"] = 1

Then, we fit the model and extract some information:

In [ ]:

mod = NBin(y, X)
   res = mod.fit()

/Users/tom.augspurger/Envs/py3/lib/python3.4/site-packages/IPython/kernel/__main__.py:3: SettingWithCopyWarning: 
   A value is trying to be set on a copy of a slice from a DataFrame.
   Try using .loc[row_indexer,col_indexer] = value instead
   
   See the the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy
     app.launch_new_instance()

Extract parameter estimates, standard errors, p-values, AIC, etc.:

In [ ]:

print('Parameters: ', res.params)
   print('Standard errors: ', res.bse)
   print('P-values: ', res.pvalues)
   print('AIC: ', res.aic)

Optimization terminated successfully.
            Current function value: 3.209014
            Iterations: 805
            Function evaluations: 1238

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.

In [ ]:

print(res.summary())

Parameters:  [ 0.2212642   0.70613942 -0.06798155 -0.12903932  2.31026565  0.44575147]
   Standard errors:  [ 0.05059259  0.07613047  0.05326097  0.06854141  0.06794692  0.01981542]
   P-values:  [  1.22297920e-005   1.76978024e-020   2.01819161e-001   5.97481600e-002
      2.15079626e-253   4.62688439e-112]
   AIC:  9604.95320583

Testing¶

We can check the results by using the statsmodels implementation of the Negative Binomial model, which uses the analytic score function and Hessian.

In [ ]:

res_nbin = sm.NegativeBinomial(y, X).fit(disp=0)
   print(res_nbin.summary())

                                 NBin Results                                 
   ==============================================================================
   Dep. Variable:                    los   Log-Likelihood:                -4797.5
   Model:                           NBin   AIC:                             9605.
   Method:            Maximum Likelihood   BIC:                             9632.
   Date:                Mon, 20 Jul 2015                                         
   Time:                        17:43:31                                         
   No. Observations:                1495                                         
   Df Residuals:                    1490                                         
   Df Model:                           4                                         
   ==============================================================================
                    coef    std err          z      P>|z|      [95.0% Conf. Int.]
   ------------------------------------------------------------------------------
   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
   ==============================================================================

In [ ]:

print(res_nbin.params)

                     NegativeBinomial Regression Results                      
   ==============================================================================
   Dep. Variable:                    los   No. Observations:                 1495
   Model:               NegativeBinomial   Df Residuals:                     1490
   Method:                           MLE   Df Model:                            4
   Date:                Mon, 20 Jul 2015   Pseudo R-squ.:                 0.01215
   Time:                        17:43:31   Log-Likelihood:                -4797.5
   converged:                       True   LL-Null:                       -4856.5
                                           LLR p-value:                 1.404e-24
   ==============================================================================
                    coef    std err          z      P>|z|      [95.0% Conf. Int.]
   ------------------------------------------------------------------------------
   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.277      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.4458      0.020     22.495      0.000         0.407     0.485
   ==============================================================================

In [ ]:

print(res_nbin.bse)

type2       0.221231
   type3       0.706175
   hmo        -0.067990
   white      -0.129065
   constant    2.310288
   alpha       0.445758
   dtype: float64

Or we could 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

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 second 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.