{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Maximum Likelihood Estimation (Generic models)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This tutorial explains how to quickly implement new maximum likelihood models in `statsmodels`. We give two examples: \n", "\n", "1. Probit model for binary dependent variables\n", "2. Negative binomial model for count data\n", "\n", "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. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Example 1: Probit model" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false, "execution": { "iopub.execute_input": "2026-07-28T10:33:33.246879Z", "iopub.status.busy": "2026-07-28T10:33:33.244058Z", "iopub.status.idle": "2026-07-28T10:33:34.536726Z", "shell.execute_reply": "2026-07-28T10:33:34.535109Z" }, "jupyter": { "outputs_hidden": false } }, "outputs": [], "source": [ "import numpy as np\n", "import statsmodels.api as sm\n", "from scipy import stats\n", "from statsmodels.base.model import GenericLikelihoodModel" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false, "execution": { "iopub.execute_input": "2026-07-28T10:33:34.539547Z", "iopub.status.busy": "2026-07-28T10:33:34.538679Z", "iopub.status.idle": "2026-07-28T10:33:34.553239Z", "shell.execute_reply": "2026-07-28T10:33:34.552687Z" }, "jupyter": { "outputs_hidden": false } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "::\n", "\n", " Number of Observations - 32\n", "\n", " Number of Variables - 4\n", "\n", " Variable name definitions::\n", "\n", " Grade - binary variable indicating whether or not a student's grade\n", " improved. 1 indicates an improvement.\n", " TUCE - Test score on economics test\n", " PSI - participation in program\n", " GPA - Student's grade point average\n", "\n", " GPA TUCE PSI\n", "0 2.66 20.0 0.0\n", "1 2.89 22.0 0.0\n", "2 3.28 24.0 0.0\n", "3 2.92 12.0 0.0\n", "4 4.00 21.0 0.0\n" ] } ], "source": [ "data = sm.datasets.spector.load_pandas()\n", "exog = data.exog\n", "endog = data.endog\n", "print(sm.datasets.spector.NOTE)\n", "print(data.exog.head())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Them, we add a constant to the matrix of regressors:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false, "execution": { "iopub.execute_input": "2026-07-28T10:33:34.613661Z", "iopub.status.busy": "2026-07-28T10:33:34.613399Z", "iopub.status.idle": "2026-07-28T10:33:34.619961Z", "shell.execute_reply": "2026-07-28T10:33:34.618855Z" }, "jupyter": { "outputs_hidden": false } }, "outputs": [], "source": [ "exog = sm.add_constant(exog, prepend=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To create your own Likelihood Model, you simply need to overwrite the loglike method." ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false, "execution": { "iopub.execute_input": "2026-07-28T10:33:34.622269Z", "iopub.status.busy": "2026-07-28T10:33:34.621517Z", "iopub.status.idle": "2026-07-28T10:33:34.626281Z", "shell.execute_reply": "2026-07-28T10:33:34.625205Z" }, "jupyter": { "outputs_hidden": false } }, "outputs": [], "source": [ "class MyProbit(GenericLikelihoodModel):\n", " def loglike(self, params):\n", " exog = self.exog\n", " endog = self.endog\n", " q = 2 * endog - 1\n", " return stats.norm.logcdf(q * np.dot(exog, params)).sum()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Estimate the model and print a summary:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false, "execution": { "iopub.execute_input": "2026-07-28T10:33:34.627908Z", "iopub.status.busy": "2026-07-28T10:33:34.627710Z", "iopub.status.idle": "2026-07-28T10:33:34.682473Z", "shell.execute_reply": "2026-07-28T10:33:34.681777Z" }, "jupyter": { "outputs_hidden": false } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Optimization terminated successfully.\n", " Current function value: 0.400588\n", " Iterations: 292\n", " Function evaluations: 494\n", " MyProbit Results \n", "==============================================================================\n", "Dep. Variable: GRADE Log-Likelihood: -12.819\n", "Model: MyProbit AIC: 33.64\n", "Method: Maximum Likelihood BIC: 39.50\n", "Date: Tue, 28 Jul 2026 \n", "Time: 10:33:34 \n", "No. Observations: 32 \n", "Df Residuals: 28 \n", "Df Model: 3 \n", "==============================================================================\n", " coef std err z P>|z| [0.025 0.975]\n", "------------------------------------------------------------------------------\n", "const -7.4523 2.542 -2.931 0.003 -12.435 -2.469\n", "GPA 1.6258 0.694 2.343 0.019 0.266 2.986\n", "TUCE 0.0517 0.084 0.617 0.537 -0.113 0.216\n", "PSI 1.4263 0.595 2.397 0.017 0.260 2.593\n", "==============================================================================\n" ] } ], "source": [ "sm_probit_manual = MyProbit(endog, exog).fit()\n", "print(sm_probit_manual.summary())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Compare your Probit implementation to ``statsmodels``' \"canned\" implementation:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false, "execution": { "iopub.execute_input": "2026-07-28T10:33:34.684923Z", "iopub.status.busy": "2026-07-28T10:33:34.684715Z", "iopub.status.idle": "2026-07-28T10:33:34.696282Z", "shell.execute_reply": "2026-07-28T10:33:34.695768Z" }, "jupyter": { "outputs_hidden": false } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Optimization terminated successfully.\n", " Current function value: 0.400588\n", " Iterations 6\n" ] } ], "source": [ "sm_probit_canned = sm.Probit(endog, exog).fit()" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false, "execution": { "iopub.execute_input": "2026-07-28T10:33:34.698494Z", "iopub.status.busy": "2026-07-28T10:33:34.698312Z", "iopub.status.idle": "2026-07-28T10:33:34.708431Z", "shell.execute_reply": "2026-07-28T10:33:34.707772Z" }, "jupyter": { "outputs_hidden": false } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "const -7.452320\n", "GPA 1.625810\n", "TUCE 0.051729\n", "PSI 1.426332\n", "dtype: float64\n", "[-7.45233176 1.62580888 0.05172971 1.42631954]\n" ] } ], "source": [ "print(sm_probit_canned.params)\n", "print(sm_probit_manual.params)" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": false, "execution": { "iopub.execute_input": "2026-07-28T10:33:34.710416Z", "iopub.status.busy": "2026-07-28T10:33:34.710230Z", "iopub.status.idle": "2026-07-28T10:33:34.723095Z", "shell.execute_reply": "2026-07-28T10:33:34.722541Z" }, "jupyter": { "outputs_hidden": false } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " const GPA TUCE PSI\n", "const 6.464166 -1.169668 -0.101173 -0.594792\n", "GPA -1.169668 0.481473 -0.018914 0.105439\n", "TUCE -0.101173 -0.018914 0.007038 0.002472\n", "PSI -0.594792 0.105439 0.002472 0.354070\n", "[[ 6.46416807e+00 -1.16966617e+00 -1.01173196e-01 -5.94789032e-01]\n", " [-1.16966617e+00 4.81472120e-01 -1.89134591e-02 1.05438221e-01]\n", " [-1.01173196e-01 -1.89134591e-02 7.03758462e-03 2.47189436e-03]\n", " [-5.94789032e-01 1.05438221e-01 2.47189436e-03 3.54069516e-01]]\n" ] } ], "source": [ "print(sm_probit_canned.cov_params())\n", "print(sm_probit_manual.cov_params())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "## Example 2: Negative Binomial Regression for Count Data\n", "\n", "Consider a negative binomial regression model for count data with\n", "log-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)) + ln \\Gamma (y_i + 1/\\alpha) - ln \\Gamma (y_i+1) - ln \\Gamma (1/\\alpha)\n", "$$\n", "\n", "with a matrix of regressors $X$, a vector of coefficients $\\beta$,\n", "and the negative binomial heterogeneity parameter $\\alpha$. \n", "\n", "Using the ``nbinom`` distribution from ``scipy``, we can write this likelihood\n", "simply as:\n" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": false, "execution": { "iopub.execute_input": "2026-07-28T10:33:34.727488Z", "iopub.status.busy": "2026-07-28T10:33:34.727301Z", "iopub.status.idle": "2026-07-28T10:33:34.735308Z", "shell.execute_reply": "2026-07-28T10:33:34.734772Z" }, "jupyter": { "outputs_hidden": false } }, "outputs": [], "source": [ "import numpy as np\n", "from scipy.stats import nbinom" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": false, "execution": { "iopub.execute_input": "2026-07-28T10:33:34.737657Z", "iopub.status.busy": "2026-07-28T10:33:34.737471Z", "iopub.status.idle": "2026-07-28T10:33:34.745151Z", "shell.execute_reply": "2026-07-28T10:33:34.744433Z" }, "jupyter": { "outputs_hidden": false } }, "outputs": [], "source": [ "def _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" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### New Model Class\n", "\n", "We create a new model class which inherits from ``GenericLikelihoodModel``:" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "collapsed": false, "execution": { "iopub.execute_input": "2026-07-28T10:33:34.747107Z", "iopub.status.busy": "2026-07-28T10:33:34.746932Z", "iopub.status.idle": "2026-07-28T10:33:34.750934Z", "shell.execute_reply": "2026-07-28T10:33:34.749841Z" }, "jupyter": { "outputs_hidden": false } }, "outputs": [], "source": [ "from statsmodels.base.model import GenericLikelihoodModel" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "collapsed": false, "execution": { "iopub.execute_input": "2026-07-28T10:33:34.753333Z", "iopub.status.busy": "2026-07-28T10:33:34.752586Z", "iopub.status.idle": "2026-07-28T10:33:34.759130Z", "shell.execute_reply": "2026-07-28T10:33:34.758089Z" }, "jupyter": { "outputs_hidden": false } }, "outputs": [], "source": [ "class NBin(GenericLikelihoodModel):\n", " def __init__(self, endog, exog, **kwds):\n", " super().__init__(endog, exog, **kwds)\n", "\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", "\n", " def fit(self, start_params=None, maxiter=10000, maxfun=5000, **kwds):\n", " # we have one additional parameter and we need to add it for summary\n", " self.exog_names.append(\"alpha\")\n", " if start_params == None:\n", " # Reasonable starting values\n", " start_params = np.append(np.zeros(self.exog.shape[1]), 0.5)\n", " # intercept\n", " start_params[-2] = np.log(self.endog.mean())\n", " return super().fit(\n", " start_params=start_params, maxiter=maxiter, maxfun=maxfun, **kwds\n", " )" ] }, { "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", " \n", "That's it! You're done!\n", "\n", "### Usage Example\n", "\n", "The [Medpar](https://raw.githubusercontent.com/vincentarelbundock/Rdatasets/doc/COUNT/medpar.html)\n", "dataset is hosted in CSV format at the [Rdatasets repository](https://raw.githubusercontent.com/vincentarelbundock/Rdatasets). We use the ``read_csv``\n", "function from the [Pandas library](https://pandas.pydata.org) to load the data\n", "in memory. We then print the first few columns: \n" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "collapsed": false, "execution": { "iopub.execute_input": "2026-07-28T10:33:34.761090Z", "iopub.status.busy": "2026-07-28T10:33:34.760905Z", "iopub.status.idle": "2026-07-28T10:33:34.764798Z", "shell.execute_reply": "2026-07-28T10:33:34.763718Z" }, "jupyter": { "outputs_hidden": false } }, "outputs": [], "source": [ "import statsmodels.api as sm" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "collapsed": false, "execution": { "iopub.execute_input": "2026-07-28T10:33:34.766565Z", "iopub.status.busy": "2026-07-28T10:33:34.766397Z", "iopub.status.idle": "2026-07-28T10:33:35.369341Z", "shell.execute_reply": "2026-07-28T10:33:35.368166Z" }, "jupyter": { "outputs_hidden": false } }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
loshmowhitediedage80typetype1type2type3provnum
040100110030001
191100110030001
231111110030001
390100110030001
410111110030001
\n", "
" ], "text/plain": [ " los hmo white died age80 type type1 type2 type3 provnum\n", "0 4 0 1 0 0 1 1 0 0 30001\n", "1 9 1 1 0 0 1 1 0 0 30001\n", "2 3 1 1 1 1 1 1 0 0 30001\n", "3 9 0 1 0 0 1 1 0 0 30001\n", "4 1 0 1 1 1 1 1 0 0 30001" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "medpar = sm.datasets.get_rdataset(\"medpar\", \"COUNT\", cache=True).data\n", "\n", "medpar.head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The model we are interested in has a vector of non-negative integers as\n", "dependent variable (``los``), and 5 regressors: ``Intercept``, ``type2``,\n", "``type3``, ``hmo``, ``white``.\n", "\n", "For estimation, we need to create two variables to hold our regressors and the outcome variable. These can be ndarrays or pandas objects." ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "collapsed": false, "execution": { "iopub.execute_input": "2026-07-28T10:33:35.371415Z", "iopub.status.busy": "2026-07-28T10:33:35.371222Z", "iopub.status.idle": "2026-07-28T10:33:35.376818Z", "shell.execute_reply": "2026-07-28T10:33:35.375760Z" }, "jupyter": { "outputs_hidden": false } }, "outputs": [], "source": [ "y = medpar.los\n", "X = medpar[[\"type2\", \"type3\", \"hmo\", \"white\"]].copy()\n", "X[\"constant\"] = 1" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Then, we fit the model and extract some information: " ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "collapsed": false, "execution": { "iopub.execute_input": "2026-07-28T10:33:35.379190Z", "iopub.status.busy": "2026-07-28T10:33:35.378427Z", "iopub.status.idle": "2026-07-28T10:33:35.732180Z", "shell.execute_reply": "2026-07-28T10:33:35.731315Z" }, "jupyter": { "outputs_hidden": false } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Optimization terminated successfully.\n", " Current function value: 3.209014\n", " Iterations: 805\n", " Function evaluations: 1238\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "/opt/hostedtoolcache/Python/3.14.6/x64/lib/python3.14/site-packages/statsmodels/base/model.py:1154: UserWarning: df_model + k_constant + k_extra differs from k_params\n", " genericmlefit = results_class(self, mlefit)\n", "/opt/hostedtoolcache/Python/3.14.6/x64/lib/python3.14/site-packages/statsmodels/base/model.py:1154: UserWarning: df_resid differs from nobs - k_params\n", " genericmlefit = results_class(self, mlefit)\n" ] } ], "source": [ "mod = NBin(y, X)\n", "res = mod.fit()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " Extract parameter estimates, standard errors, p-values, AIC, etc.:" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "collapsed": false, "execution": { "iopub.execute_input": "2026-07-28T10:33:35.740457Z", "iopub.status.busy": "2026-07-28T10:33:35.740277Z", "iopub.status.idle": "2026-07-28T10:33:35.749123Z", "shell.execute_reply": "2026-07-28T10:33:35.748548Z" }, "jupyter": { "outputs_hidden": false } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Parameters: [ 0.2212642 0.70613942 -0.06798155 -0.12903932 2.31026565 0.44575147]\n", "Standard errors: [0.05059259 0.07613047 0.05326096 0.0685414 0.06794696 0.01981542]\n", "P-values: [1.22298084e-005 1.76979047e-020 2.01819053e-001 5.97481232e-002\n", " 2.15207253e-253 4.62685274e-112]\n", "AIC: 9604.95320583016\n" ] } ], "source": [ "print(\"Parameters: \", res.params)\n", "print(\"Standard errors: \", res.bse)\n", "print(\"P-values: \", res.pvalues)\n", "print(\"AIC: \", res.aic)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As usual, you can obtain a full list of available information by typing\n", "``dir(res)``.\n", "We can also look at the summary of the estimation results." ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "collapsed": false, "execution": { "iopub.execute_input": "2026-07-28T10:33:35.753576Z", "iopub.status.busy": "2026-07-28T10:33:35.753185Z", "iopub.status.idle": "2026-07-28T10:33:35.761686Z", "shell.execute_reply": "2026-07-28T10:33:35.761179Z" }, "jupyter": { "outputs_hidden": false } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " NBin Results \n", "==============================================================================\n", "Dep. Variable: los Log-Likelihood: -4797.5\n", "Model: NBin AIC: 9605.\n", "Method: Maximum Likelihood BIC: 9632.\n", "Date: Tue, 28 Jul 2026 \n", "Time: 10:33:35 \n", "No. Observations: 1495 \n", "Df Residuals: 1490 \n", "Df Model: 4 \n", "==============================================================================\n", " coef std err z P>|z| [0.025 0.975]\n", "------------------------------------------------------------------------------\n", "type2 0.2213 0.051 4.373 0.000 0.122 0.320\n", "type3 0.7061 0.076 9.275 0.000 0.557 0.855\n", "hmo -0.0680 0.053 -1.276 0.202 -0.172 0.036\n", "white -0.1290 0.069 -1.883 0.060 -0.263 0.005\n", "constant 2.3103 0.068 34.001 0.000 2.177 2.443\n", "alpha 0.4458 0.020 22.495 0.000 0.407 0.485\n", "==============================================================================\n" ] } ], "source": [ "print(res.summary())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Testing" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can check the results by using the statsmodels implementation of the Negative Binomial model, which uses the analytic score function and Hessian." ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "collapsed": false, "execution": { "iopub.execute_input": "2026-07-28T10:33:35.763579Z", "iopub.status.busy": "2026-07-28T10:33:35.763433Z", "iopub.status.idle": "2026-07-28T10:33:35.813943Z", "shell.execute_reply": "2026-07-28T10:33:35.812795Z" }, "jupyter": { "outputs_hidden": false } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " NegativeBinomial Regression Results \n", "==============================================================================\n", "Dep. Variable: los No. Observations: 1495\n", "Model: NegativeBinomial Df Residuals: 1490\n", "Method: MLE Df Model: 4\n", "Date: Tue, 28 Jul 2026 Pseudo R-squ.: 0.01215\n", "Time: 10:33:35 Log-Likelihood: -4797.5\n", "converged: True LL-Null: -4856.5\n", "Covariance Type: nonrobust LLR p-value: 1.404e-24\n", "==============================================================================\n", " coef std err z P>|z| [0.025 0.975]\n", "------------------------------------------------------------------------------\n", "type2 0.2212 0.051 4.373 0.000 0.122 0.320\n", "type3 0.7062 0.076 9.276 0.000 0.557 0.855\n", "hmo -0.0680 0.053 -1.276 0.202 -0.172 0.036\n", "white -0.1291 0.069 -1.883 0.060 -0.263 0.005\n", "constant 2.3103 0.068 34.001 0.000 2.177 2.443\n", "alpha 0.4457 0.020 22.495 0.000 0.407 0.485\n", "==============================================================================\n" ] } ], "source": [ "res_nbin = sm.NegativeBinomial(y, X).fit(disp=0)\n", "print(res_nbin.summary())" ] }, { "cell_type": "code", "execution_count": 20, "metadata": { "collapsed": false, "execution": { "iopub.execute_input": "2026-07-28T10:33:35.816041Z", "iopub.status.busy": "2026-07-28T10:33:35.815425Z", "iopub.status.idle": "2026-07-28T10:33:35.823141Z", "shell.execute_reply": "2026-07-28T10:33:35.822329Z" }, "jupyter": { "outputs_hidden": false } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "type2 0.221218\n", "type3 0.706173\n", "hmo -0.067987\n", "white -0.129053\n", "constant 2.310279\n", "alpha 0.445748\n", "dtype: float64\n" ] } ], "source": [ "print(res_nbin.params)" ] }, { "cell_type": "code", "execution_count": 21, "metadata": { "collapsed": false, "execution": { "iopub.execute_input": "2026-07-28T10:33:35.825377Z", "iopub.status.busy": "2026-07-28T10:33:35.824693Z", "iopub.status.idle": "2026-07-28T10:33:35.830192Z", "shell.execute_reply": "2026-07-28T10:33:35.829306Z" }, "jupyter": { "outputs_hidden": false } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "type2 0.050592\n", "type3 0.076131\n", "hmo 0.053261\n", "white 0.068541\n", "constant 0.067947\n", "alpha 0.019815\n", "dtype: float64\n" ] } ], "source": [ "print(res_nbin.bse)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Or we could compare them to results obtained using the MASS implementation for R:\n", "\n", " url = 'https://raw.githubusercontent.com/vincentarelbundock/Rdatasets/csv/COUNT/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", "\n", "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." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.14.6" } }, "nbformat": 4, "nbformat_minor": 4 }