statsmodels example in ipython notebook

example_gmle.ipynb

{
 "metadata": {
  "name": "example_gmle"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "# Generic Maximum Likelihood Model\n\nThis 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. \n\n## Example: Negative Binomial Regression for Count Data\n\nConsider a negative binomial regression model for count data with\nlog-likelihood (type NB-2) function expressed as:\n\n$$\n    \\mathcal{L}(\\beta_j; y, \\alpha) = \\sum_{i=1}^n y_i ln \n    \\left ( \\frac{\\alpha exp(X_i'\\beta)}{1+\\alpha exp(X_i'\\beta)} \\right ) -\n    \\frac{1}{\\alpha} ln(1+\\alpha exp(X_i'\\beta)) \\\\\n$$\n$$\n    + ln \\Gamma (y_i + 1/\\alpha) - ln \\Gamma (y_i+1) - ln \\Gamma (1/\\alpha)\n$$\n\nwith a matrix of regressors $X$, a vector of coefficients $\\beta$,\nand the negative binomial heterogeneity parameter $\\alpha$. \n\nUsing the ``nbinom`` distribution from ``scipy``, we can write this likelihood\nsimply as:\n"
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "import numpy as np\nfrom scipy.stats import nbinom\ndef _ll_nb2(y, X, beta, alph):\n    mu = np.exp(np.dot(X, beta))\n    size = 1/alph\n    prob = size/(size+mu)\n    ll = nbinom.logpmf(y, size, prob)\n    return ll",
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "## New Model Class\n\nWe create a new model class which inherits from ``GenericLikelihoodModel``:"
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "from statsmodels.base.model import GenericLikelihoodModel\nclass NBin(GenericLikelihoodModel):\n    def __init__(self, endog, exog, **kwds):\n        super(NBin, self).__init__(endog, exog, **kwds)\n    def nloglikeobs(self, params):\n        alph = params[-1]\n        beta = params[:-1]\n        ll = _ll_nb2(self.endog, self.exog, beta, alph)\n        return -ll \n    def fit(self, start_params=None, maxiter=10000, maxfun=5000, **kwds):\n        if start_params == None:\n            # Reasonable starting values\n            start_params = np.append(np.zeros(self.exog.shape[1]), .5)\n            start_params[0] = np.log(self.endog.mean())\n        return super(NBin, self).fit(start_params=start_params, \n                                     maxiter=maxiter, maxfun=maxfun, \n                                     **kwds) ",
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "Two important things to notice: \n\n+ ``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). \n+ ``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.\n   \nThat's it! You're done!\n\nUsage Example\n-------------\n\nThe [Medpar](http://vincentarelbundock.github.com/Rdatasets/doc/medpar.html)\ndataset is hosted in CSV format at the [Rdatasets repository](http://vincentarelbundock.github.com/Rdatasets). We use the ``read_csv``\nfunction from the [Pandas library](http://pandas.pydata.org) to load the data\nin memory. We then print the first few columns: \n"
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "import pandas as pd\nurl = 'http://vincentarelbundock.github.com/Rdatasets/csv/medpar.csv'\nmedpar = pd.read_csv(url)\nmedpar.head()",
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "html": "<div style=\"max-height:1000px;max-width:1500px;overflow:auto;\">\n<table border=\"1\">\n  <thead>\n    <tr>\n      <th></th>\n      <th>Unnamed: 0</th>\n      <th>los</th>\n      <th>hmo</th>\n      <th>white</th>\n      <th>died</th>\n      <th>age80</th>\n      <th>type</th>\n      <th>type1</th>\n      <th>type2</th>\n      <th>type3</th>\n      <th>provnum</th>\n    </tr>\n  </thead>\n  <tbody>\n    <tr>\n      <td><strong>0</strong></td>\n      <td> 1</td>\n      <td> 4</td>\n      <td> 0</td>\n      <td> 1</td>\n      <td> 0</td>\n      <td> 0</td>\n      <td> 1</td>\n      <td> 1</td>\n      <td> 0</td>\n      <td> 0</td>\n      <td> 30001</td>\n    </tr>\n    <tr>\n      <td><strong>1</strong></td>\n      <td> 2</td>\n      <td> 9</td>\n      <td> 1</td>\n      <td> 1</td>\n      <td> 0</td>\n      <td> 0</td>\n      <td> 1</td>\n      <td> 1</td>\n      <td> 0</td>\n      <td> 0</td>\n      <td> 30001</td>\n    </tr>\n    <tr>\n      <td><strong>2</strong></td>\n      <td> 3</td>\n      <td> 3</td>\n      <td> 1</td>\n      <td> 1</td>\n      <td> 1</td>\n      <td> 1</td>\n      <td> 1</td>\n      <td> 1</td>\n      <td> 0</td>\n      <td> 0</td>\n      <td> 30001</td>\n    </tr>\n    <tr>\n      <td><strong>3</strong></td>\n      <td> 4</td>\n      <td> 9</td>\n      <td> 0</td>\n      <td> 1</td>\n      <td> 0</td>\n      <td> 0</td>\n      <td> 1</td>\n      <td> 1</td>\n      <td> 0</td>\n      <td> 0</td>\n      <td> 30001</td>\n    </tr>\n    <tr>\n      <td><strong>4</strong></td>\n      <td> 5</td>\n      <td> 1</td>\n      <td> 0</td>\n      <td> 1</td>\n      <td> 1</td>\n      <td> 1</td>\n      <td> 1</td>\n      <td> 1</td>\n      <td> 0</td>\n      <td> 0</td>\n      <td> 30001</td>\n    </tr>\n  </tbody>\n</table>\n</div>",
       "output_type": "pyout",
       "prompt_number": 4,
       "text": "   Unnamed: 0  los  hmo  white  died  age80  type  type1  type2  type3  provnum\n0           1    4    0      1     0      0     1      1      0      0    30001\n1           2    9    1      1     0      0     1      1      0      0    30001\n2           3    3    1      1     1      1     1      1      0      0    30001\n3           4    9    0      1     0      0     1      1      0      0    30001\n4           5    1    0      1     1      1     1      1      0      0    30001"
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "The model we are interested in has a vector of non-negative integers as\ndependent variable (``los``), and 5 regressors: ``Intercept``, ``type2``,\n``type3``, ``hmo``, ``white``.\n\nFor estimation, we need to create 2 numpy arrays (pandas DataFrame should also\nwork): a 1d array of length *N* to hold ``los`` values, and a *N* by 5\narray to hold our 5 regressors. These arrays can be constructed manually or\nusing any number of helper functions; the details matter little for our current\npurposes.  Here, we build the arrays we need using the [Patsy](http://patsy.readthedocs.org) package:"
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "import patsy\ny, X = patsy.dmatrices('los~type2+type3+hmo+white', medpar)\nprint y[:5]\nprint X[:5] ",
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": "[[ 4.]\n [ 9.]\n [ 3.]\n [ 9.]\n [ 1.]]\n[[ 1.  0.  0.  0.  1.]\n [ 1.  0.  0.  1.  1.]\n [ 1.  0.  0.  1.  1.]\n [ 1.  0.  0.  0.  1.]\n [ 1.  0.  0.  0.  1.]]\n"
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "Then, we fit the model and extract some information: "
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "mod = NBin(y, X)\nres = mod.fit()",
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": "Optimization terminated successfully.\n         Current function value: 4797.476603\n         Iterations: 805\n         Function evaluations: 1238\n"
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": " Extract parameter estimates, standard errors, p-values, AIC, etc.:"
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "print 'Parameters: ', res.params\nprint 'Standard errors: ', res.bse\nprint 'P-values: ', res.pvalues\nprint 'AIC: ', res.aic",
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": "Parameters:  [ 2.31026565  0.2212642   0.70613942 -0.06798155 -0.12903932  0.44575147]\nStandard errors:  [ 0.06221065  0.05014587  0.07613956  0.05417982  0.06196514  0.01978691]\nP-values:  [ 0.          0.00001022  0.          0.20957345  0.03730137  0.        ]\nAIC:  9606.95320583\n"
      }
     ],
     "prompt_number": 9
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "As usual, you can obtain a full list of available information by typing\n``dir(res)``."
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "## Testing \n \nTo ensure that the above results are sound, we compare them to results\nobtained using the MASS implementation for R:\n\n    url = 'http://vincentarelbundock.github.com/Rdatasets/csv/medpar.csv'\n    medpar = read.csv(url)\n    f = los~factor(type)+hmo+white\n    \n    library(MASS)\n    mod = glm.nb(f, medpar)\n    coef(summary(mod))\n                     Estimate Std. Error   z value      Pr(>|z|)\n    (Intercept)    2.31027893 0.06744676 34.253370 3.885556e-257\n    factor(type)2  0.22124898 0.05045746  4.384861  1.160597e-05\n    factor(type)3  0.70615882 0.07599849  9.291748  1.517751e-20\n    hmo           -0.06795522 0.05321375 -1.277024  2.015939e-01\n    white         -0.12906544 0.06836272 -1.887951  5.903257e-02\n\n## Numerical precision \n\nThe ``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. \n"
    }
   ],
   "metadata": {}
  }
 ]
}