Loading [MathJax]/jax/output/HTML-CSS/jax.js
MENU
Edit
In Place Editing Advanced Editing

Maximum Likelihood Estimation (Generic models)

Maximum Likelihood Estimation (Generic models)

Link to Notebook GitHub

This tutorial explains how to quickly implement new maximum likelihood models in statsmodels. We give two examples:

  1. Probit model for binary dependent variables
  2. 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 [1]:
1
2
3
4
5
<span class="kn">from</span> <span class="nn">__future__</span> <span class="kn">import</span> <span class="n">print_function</span>
<span class="kn">import</span> <span class="nn">numpy</span> <span class="kn">as</span> <span class="nn">np</span>
<span class="kn">from</span> <span class="nn">scipy</span> <span class="kn">import</span> <span class="n">stats</span>
<span class="kn">import</span> <span class="nn">statsmodels.api</span> <span class="kn">as</span> <span class="nn">sm</span>
<span class="kn">from</span> <span class="nn">statsmodels.base.model</span> <span class="kn">import</span> <span class="n">GenericLikelihoodModel</span>

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 [2]:
1
2
3
4
5
<span class="n">data</span> <span class="o">=</span> <span class="n">sm</span><span class="o">.</span><span class="n">datasets</span><span class="o">.</span><span class="n">spector</span><span class="o">.</span><span class="n">load_pandas</span><span class="p">()</span>
<span class="n">exog</span> <span class="o">=</span> <span class="n">data</span><span class="o">.</span><span class="n">exog</span>
<span class="n">endog</span> <span class="o">=</span> <span class="n">data</span><span class="o">.</span><span class="n">endog</span>
<span class="k">print</span><span class="p">(</span><span class="n">sm</span><span class="o">.</span><span class="n">datasets</span><span class="o">.</span><span class="n">spector</span><span class="o">.</span><span class="n">NOTE</span><span class="p">)</span>
<span class="k">print</span><span class="p">(</span><span class="n">data</span><span class="o">.</span><span class="n">exog</span><span class="o">.</span><span class="n">head</span><span class="p">())</span>
::

    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

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

In [3]:
1
<span class="n">exog</span> <span class="o">=</span> <span class="n">sm</span><span class="o">.</span><span class="n">add_constant</span><span class="p">(</span><span class="n">exog</span><span class="p">,</span> <span class="n">prepend</span><span class="o">=</span><span class="bp">True</span><span class="p">)</span>

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

In [4]:
1
2
3
4
5
6
<span class="k">class</span> <span class="nc">MyProbit</span><span class="p">(</span><span class="n">GenericLikelihoodModel</span><span class="p">):</span>
    <span class="k">def</span> <span class="nf">loglike</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">params</span><span class="p">):</span>
        <span class="n">exog</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">exog</span>
        <span class="n">endog</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">endog</span>
        <span class="n">q</span> <span class="o">=</span> <span class="mi">2</span> <span class="o">*</span> <span class="n">endog</span> <span class="o">-</span> <span class="mi">1</span>
        <span class="k">return</span> <span class="n">stats</span><span class="o">.</span><span class="n">norm</span><span class="o">.</span><span class="n">logcdf</span><span class="p">(</span><span class="n">q</span><span class="o">*</span><span class="n">np</span><span class="o">.</span><span class="n">dot</span><span class="p">(</span><span class="n">exog</span><span class="p">,</span> <span class="n">params</span><span class="p">))</span><span class="o">.</span><span class="n">sum</span><span class="p">()</span>

Estimate the model and print a summary:

In [5]:
1
2
<span class="n">sm_probit_manual</span> <span class="o">=</span> <span class="n">MyProbit</span><span class="p">(</span><span class="n">endog</span><span class="p">,</span> <span class="n">exog</span><span class="p">)</span><span class="o">.</span><span class="n">fit</span><span class="p">()</span>
<span class="k">print</span><span class="p">(</span><span class="n">sm_probit_manual</span><span class="o">.</span><span class="n">summary</span><span class="p">())</span>
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, 02 Dec 2014
Time:                        12:54:20
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
==============================================================================

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

In [6]:
1
<span class="n">sm_probit_canned</span> <span class="o">=</span> <span class="n">sm</span><span class="o">.</span><span class="n">Probit</span><span class="p">(</span><span class="n">endog</span><span class="p">,</span> <span class="n">exog</span><span class="p">)</span><span class="o">.</span><span class="n">fit</span><span class="p">()</span>
Optimization terminated successfully.
         Current function value: 0.400588
         Iterations 6
In [7]:
1
2
<span class="k">print</span><span class="p">(</span><span class="n">sm_probit_canned</span><span class="o">.</span><span class="n">params</span><span class="p">)</span>
<span class="k">print</span><span class="p">(</span><span class="n">sm_probit_manual</span><span class="o">.</span><span class="n">params</span><span class="p">)</span>
const   -7.452320
GPA      1.625810
TUCE     0.051729
PSI      1.426332
dtype: float64
[-7.4523  1.6258  0.0517  1.4263]
In [8]:
1
2
<span class="k">print</span><span class="p">(</span><span class="n">sm_probit_canned</span><span class="o">.</span><span class="n">cov_params</span><span class="p">())</span>
<span class="k">print</span><span class="p">(</span><span class="n">sm_probit_manual</span><span class="o">.</span><span class="n">cov_params</span><span class="p">())</span>
          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.4642 -1.1697 -0.1012 -0.5948]
 [-1.1697  0.4815 -0.0189  0.1054]
 [-0.1012 -0.0189  0.007   0.0025]
 [-0.5948  0.1054  0.0025  0.3541]]

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:

L(βj;y,α)=n∑i=1yiln(αexp(X′iβ)1+αexp(X′iβ))−1αln(1+αexp(X′iβ))+lnΓ(yi+1/α)−lnΓ(yi+1)−lnΓ(1/α)

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 [9]:
1
2
<span class="kn">import</span> <span class="nn">numpy</span> <span class="kn">as</span> <span class="nn">np</span>
<span class="kn">from</span> <span class="nn">scipy.stats</span> <span class="kn">import</span> <span class="n">nbinom</span>
In [10]:
1
2
3
4
5
6
<span class="k">def</span> <span class="nf">_ll_nb2</span><span class="p">(</span><span class="n">y</span><span class="p">,</span> <span class="n">X</span><span class="p">,</span> <span class="n">beta</span><span class="p">,</span> <span class="n">alph</span><span class="p">):</span>
    <span class="n">mu</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">exp</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">dot</span><span class="p">(</span><span class="n">X</span><span class="p">,</span> <span class="n">beta</span><span class="p">))</span>
    <span class="n">size</span> <span class="o">=</span> <span class="mi">1</span><span class="o">/</span><span class="n">alph</span>
    <span class="n">prob</span> <span class="o">=</span> <span class="n">size</span><span class="o">/</span><span class="p">(</span><span class="n">size</span><span class="o">+</span><span class="n">mu</span><span class="p">)</span>
    <span class="n">ll</span> <span class="o">=</span> <span class="n">nbinom</span><span class="o">.</span><span class="n">logpmf</span><span class="p">(</span><span class="n">y</span><span class="p">,</span> <span class="n">size</span><span class="p">,</span> <span class="n">prob</span><span class="p">)</span>
    <span class="k">return</span> <span class="n">ll</span>

New Model Class

We create a new model class which inherits from GenericLikelihoodModel:

In [11]:
1
<span class="kn">from</span> <span class="nn">statsmodels.base.model</span> <span class="kn">import</span> <span class="n">GenericLikelihoodModel</span>
In [12]:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
<span class="k">class</span> <span class="nc">NBin</span><span class="p">(</span><span class="n">GenericLikelihoodModel</span><span class="p">):</span>
    <span class="k">def</span> <span class="nf">__init__</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">endog</span><span class="p">,</span> <span class="n">exog</span><span class="p">,</span> <span class="o">**</span><span class="n">kwds</span><span class="p">):</span>
        <span class="nb">super</span><span class="p">(</span><span class="n">NBin</span><span class="p">,</span> <span class="bp">self</span><span class="p">)</span><span class="o">.</span><span class="n">__init__</span><span class="p">(</span><span class="n">endog</span><span class="p">,</span> <span class="n">exog</span><span class="p">,</span> <span class="o">**</span><span class="n">kwds</span><span class="p">)</span>
 
    <span class="k">def</span> <span class="nf">nloglikeobs</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">params</span><span class="p">):</span>
        <span class="n">alph</span> <span class="o">=</span> <span class="n">params</span><span class="p">[</span><span class="o">-</span><span class="mi">1</span><span class="p">]</span>
        <span class="n">beta</span> <span class="o">=</span> <span class="n">params</span><span class="p">[:</span><span class="o">-</span><span class="mi">1</span><span class="p">]</span>
        <span class="n">ll</span> <span class="o">=</span> <span class="n">_ll_nb2</span><span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">endog</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">exog</span><span class="p">,</span> <span class="n">beta</span><span class="p">,</span> <span class="n">alph</span><span class="p">)</span>
        <span class="k">return</span> <span class="o">-</span><span class="n">ll</span>
 
    <span class="k">def</span> <span class="nf">fit</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">start_params</span><span class="o">=</span><span class="bp">None</span><span class="p">,</span> <span class="n">maxiter</span><span class="o">=</span><span class="mi">10000</span><span class="p">,</span> <span class="n">maxfun</span><span class="o">=</span><span class="mi">5000</span><span class="p">,</span> <span class="o">**</span><span class="n">kwds</span><span class="p">):</span>
        <span class="c"># we have one additional parameter and we need to add it for summary</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">exog_names</span><span class="o">.</span><span class="n">append</span><span class="p">(</span><span class="s">'alpha'</span><span class="p">)</span>
        <span class="k">if</span> <span class="n">start_params</span> <span class="o">==</span> <span class="bp">None</span><span class="p">:</span>
            <span class="c"># Reasonable starting values</span>
            <span class="n">start_params</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">append</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">exog</span><span class="o">.</span><span class="n">shape</span><span class="p">[</span><span class="mi">1</span><span class="p">]),</span> <span class="o">.</span><span class="mi">5</span><span class="p">)</span>
            <span class="c"># intercept</span>
            <span class="n">start_params</span><span class="p">[</span><span class="o">-</span><span class="mi">2</span><span class="p">]</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">log</span><span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">endog</span><span class="o">.</span><span class="n">mean</span><span class="p">())</span>
        <span class="k">return</span> <span class="nb">super</span><span class="p">(</span><span class="n">NBin</span><span class="p">,</span> <span class="bp">self</span><span class="p">)</span><span class="o">.</span><span class="n">fit</span><span class="p">(</span><span class="n">start_params</span><span class="o">=</span><span class="n">start_params</span><span class="p">,</span>
                                     <span class="n">maxiter</span><span class="o">=</span><span class="n">maxiter</span><span class="p">,</span> <span class="n">maxfun</span><span class="o">=</span><span class="n">maxfun</span><span class="p">,</span>
                                     <span class="o">**</span><span class="n">kwds</span><span class="p">)</span>

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 [13]:
1
<span class="kn">import</span> <span class="nn">statsmodels.api</span> <span class="kn">as</span> <span class="nn">sm</span>
In [14]:
1
2
3
<span class="n">medpar</span> <span class="o">=</span> <span class="n">sm</span><span class="o">.</span><span class="n">datasets</span><span class="o">.</span><span class="n">get_rdataset</span><span class="p">(</span><span class="s">"medpar"</span><span class="p">,</span> <span class="s">"COUNT"</span><span class="p">,</span> <span class="n">cache</span><span class="o">=</span><span class="bp">True</span><span class="p">)</span><span class="o">.</span><span class="n">data</span>
 
<span class="n">medpar</span><span class="o">.</span><span class="n">head</span><span class="p">()</span>
Out[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

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 [15]:
1
2
3
<span class="n">y</span> <span class="o">=</span> <span class="n">medpar</span><span class="o">.</span><span class="n">los</span>
<span class="n">X</span> <span class="o">=</span> <span class="n">medpar</span><span class="p">[[</span><span class="s">"type2"</span><span class="p">,</span> <span class="s">"type3"</span><span class="p">,</span> <span class="s">"hmo"</span><span class="p">,</span> <span class="s">"white"</span><span class="p">]]</span>
<span class="n">X</span><span class="p">[</span><span class="s">"constant"</span><span class="p">]</span> <span class="o">=</span> <span class="mi">1</span>

Then, we fit the model and extract some information:

In [16]:
1
2
<span class="n">mod</span> <span class="o">=</span> <span class="n">NBin</span><span class="p">(</span><span class="n">y</span><span class="p">,</span> <span class="n">X</span><span class="p">)</span>
<span class="n">res</span> <span class="o">=</span> <span class="n">mod</span><span class="o">.</span><span class="n">fit</span><span class="p">()</span>
Optimization terminated successfully.
         Current function value: 3.209014
         Iterations: 805
         Function evaluations: 1238

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

In [17]:
1
2
3
4
<span class="k">print</span><span class="p">(</span><span class="s">'Parameters: '</span><span class="p">,</span> <span class="n">res</span><span class="o">.</span><span class="n">params</span><span class="p">)</span>
<span class="k">print</span><span class="p">(</span><span class="s">'Standard errors: '</span><span class="p">,</span> <span class="n">res</span><span class="o">.</span><span class="n">bse</span><span class="p">)</span>
<span class="k">print</span><span class="p">(</span><span class="s">'P-values: '</span><span class="p">,</span> <span class="n">res</span><span class="o">.</span><span class="n">pvalues</span><span class="p">)</span>
<span class="k">print</span><span class="p">(</span><span class="s">'AIC: '</span><span class="p">,</span> <span class="n">res</span><span class="o">.</span><span class="n">aic</span><span class="p">)</span>
Parameters:  [ 0.2213  0.7061 -0.068  -0.129   2.3103  0.4458]
Standard errors:  [ 0.0506  0.0761  0.0533  0.0685  0.0679  0.0198]
P-values:  [ 0.      0.      0.2018  0.0597  0.      0.    ]
AIC:  9604.95320583

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 [18]:
1
<span class="k">print</span><span class="p">(</span><span class="n">res</span><span class="o">.</span><span class="n">summary</span><span class="p">())</span>
                                 NBin Results
==============================================================================
Dep. Variable:                    los   Log-Likelihood:                -4797.5
Model:                           NBin   AIC:                             9605.
Method:            Maximum Likelihood   BIC:                             9632.
Date:                Tue, 02 Dec 2014
Time:                        12:54:24
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
==============================================================================

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 [19]:
1
2
<span class="n">res_nbin</span> <span class="o">=</span> <span class="n">sm</span><span class="o">.</span><span class="n">NegativeBinomial</span><span class="p">(</span><span class="n">y</span><span class="p">,</span> <span class="n">X</span><span class="p">)</span><span class="o">.</span><span class="n">fit</span><span class="p">(</span><span class="n">disp</span><span class="o">=</span><span class="mi">0</span><span class="p">)</span>
<span class="k">print</span><span class="p">(</span><span class="n">res_nbin</span><span class="o">.</span><span class="n">summary</span><span class="p">())</span>
                     NegativeBinomial Regression Results
==============================================================================
Dep. Variable:                    los   No. Observations:                 1495
Model:               NegativeBinomial   Df Residuals:                     1490
Method:                           MLE   Df Model:                            4
Date:                Tue, 02 Dec 2014   Pseudo R-squ.:                 0.01215
Time:                        12:54:24   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 [20]:
1
<span class="k">print</span><span class="p">(</span><span class="n">res_nbin</span><span class="o">.</span><span class="n">params</span><span class="p">)</span>
type2       0.221231
type3       0.706175
hmo        -0.067990
white      -0.129065
constant    2.310288
alpha       0.445758
dtype: float64
In [21]:
1
<span class="k">print</span><span class="p">(</span><span class="n">res_nbin</span><span class="o">.</span><span class="n">bse</span><span class="p">)</span>
type2       0.050592
type3       0.076132
hmo         0.053261
white       0.068542
constant    0.067947
alpha       0.019816
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.

doc_statsmodels
doc_statsmodels
2025-01-10 15:47:30
Comments
Leave a Comment

Please login to continue.