# Discrete Choice Models Overview¶

[1]:

import numpy as np
import statsmodels.api as sm


## Data¶

Load data from Spector and Mazzeo (1980). Examples follow Greene’s Econometric Analysis Ch. 21 (5th Edition).

[2]:

spector_data = sm.datasets.spector.load()


Inspect the data:

[3]:

print(spector_data.exog.head())

    GPA  TUCE  PSI  const
0  2.66  20.0  0.0    1.0
1  2.89  22.0  0.0    1.0
2  3.28  24.0  0.0    1.0
3  2.92  12.0  0.0    1.0
4  4.00  21.0  0.0    1.0
0    0.0
1    0.0
2    0.0
3    0.0
4    1.0


## Linear Probability Model (OLS)¶

[4]:

lpm_mod = sm.OLS(spector_data.endog, spector_data.exog)
lpm_res = lpm_mod.fit()
print("Parameters: ", lpm_res.params[:-1])

Parameters:  GPA     0.463852
TUCE    0.010495
PSI     0.378555
dtype: float64


## Logit Model¶

[5]:

logit_mod = sm.Logit(spector_data.endog, spector_data.exog)
logit_res = logit_mod.fit(disp=0)
print("Parameters: ", logit_res.params)

Parameters:  GPA       2.826113
TUCE      0.095158
PSI       2.378688
const   -13.021347
dtype: float64


Marginal Effects

[6]:

margeff = logit_res.get_margeff()
print(margeff.summary())

        Logit Marginal Effects
=====================================
Method:                          dydx
At:                           overall
==============================================================================
dy/dx    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
GPA            0.3626      0.109      3.313      0.001       0.148       0.577
TUCE           0.0122      0.018      0.686      0.493      -0.023       0.047
PSI            0.3052      0.092      3.304      0.001       0.124       0.486
==============================================================================


As in all the discrete data models presented below, we can print a nice summary of results:

[7]:

print(logit_res.summary())

                           Logit Regression Results
==============================================================================
Dep. Variable:                  GRADE   No. Observations:                   32
Model:                          Logit   Df Residuals:                       28
Method:                           MLE   Df Model:                            3
Date:                Sat, 25 May 2024   Pseudo R-squ.:                  0.3740
Time:                        15:57:58   Log-Likelihood:                -12.890
converged:                       True   LL-Null:                       -20.592
Covariance Type:            nonrobust   LLR p-value:                  0.001502
==============================================================================
coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
GPA            2.8261      1.263      2.238      0.025       0.351       5.301
TUCE           0.0952      0.142      0.672      0.501      -0.182       0.373
PSI            2.3787      1.065      2.234      0.025       0.292       4.465
const        -13.0213      4.931     -2.641      0.008     -22.687      -3.356
==============================================================================


## Probit Model¶

[8]:

probit_mod = sm.Probit(spector_data.endog, spector_data.exog)
probit_res = probit_mod.fit()
probit_margeff = probit_res.get_margeff()
print("Parameters: ", probit_res.params)
print("Marginal effects: ")
print(probit_margeff.summary())

Optimization terminated successfully.
Current function value: 0.400588
Iterations 6
Parameters:  GPA      1.625810
TUCE     0.051729
PSI      1.426332
const   -7.452320
dtype: float64
Marginal effects:
Probit Marginal Effects
=====================================
Method:                          dydx
At:                           overall
==============================================================================
dy/dx    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
GPA            0.3608      0.113      3.182      0.001       0.139       0.583
TUCE           0.0115      0.018      0.624      0.533      -0.025       0.048
PSI            0.3165      0.090      3.508      0.000       0.140       0.493
==============================================================================


## Multinomial Logit¶

Load data from the American National Election Studies:

[9]:

anes_data = sm.datasets.anes96.load()
anes_exog = anes_data.exog


Inspect the data:

[10]:

print(anes_data.exog.head())

   logpopul  selfLR   age  educ  income
0 -2.302585     7.0  36.0   3.0     1.0
1  5.247550     3.0  20.0   4.0     1.0
2  3.437208     2.0  24.0   6.0     1.0
3  4.420045     3.0  28.0   6.0     1.0
4  6.461624     5.0  68.0   6.0     1.0
0    6.0
1    1.0
2    1.0
3    1.0
4    0.0
Name: PID, dtype: float64


Fit MNL model:

[11]:

mlogit_mod = sm.MNLogit(anes_data.endog, anes_exog)
mlogit_res = mlogit_mod.fit()
print(mlogit_res.params)

Optimization terminated successfully.
Current function value: 1.548647
Iterations 7
0         1         2         3         4          5
const    -0.373402 -2.250913 -3.665584 -7.613843 -7.060478 -12.105751
logpopul -0.011536 -0.088751 -0.105967 -0.091557 -0.093285  -0.140881
selfLR    0.297714  0.391669  0.573451  1.278772  1.346962   2.070080
age      -0.024945 -0.022898 -0.014851 -0.008681 -0.017904  -0.009433
educ      0.082491  0.181043 -0.007152  0.199828  0.216939   0.321926
income    0.005197  0.047874  0.057575  0.084498  0.080958   0.108894


## Poisson¶

Load the Rand data. Note that this example is similar to Cameron and Trivedi’s Microeconometrics Table 20.5, but it is slightly different because of minor changes in the data.

[12]:

rand_data = sm.datasets.randhie.load()
rand_exog = rand_data.exog


Fit Poisson model:

[13]:

poisson_mod = sm.Poisson(rand_data.endog, rand_exog)
poisson_res = poisson_mod.fit(method="newton")
print(poisson_res.summary())

Optimization terminated successfully.
Current function value: 3.091609
Iterations 6
Poisson Regression Results
==============================================================================
Dep. Variable:                  mdvis   No. Observations:                20190
Model:                        Poisson   Df Residuals:                    20180
Method:                           MLE   Df Model:                            9
Date:                Sat, 25 May 2024   Pseudo R-squ.:                 0.06343
Time:                        15:57:58   Log-Likelihood:                -62420.
converged:                       True   LL-Null:                       -66647.
Covariance Type:            nonrobust   LLR p-value:                     0.000
==============================================================================
coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
lncoins       -0.0525      0.003    -18.216      0.000      -0.058      -0.047
idp           -0.2471      0.011    -23.272      0.000      -0.268      -0.226
lpi            0.0353      0.002     19.302      0.000       0.032       0.039
fmde          -0.0346      0.002    -21.439      0.000      -0.038      -0.031
physlm         0.2717      0.012     22.200      0.000       0.248       0.296
disea          0.0339      0.001     60.098      0.000       0.033       0.035
hlthg         -0.0126      0.009     -1.366      0.172      -0.031       0.005
hlthf          0.0541      0.015      3.531      0.000       0.024       0.084
hlthp          0.2061      0.026      7.843      0.000       0.155       0.258
const          0.7004      0.011     62.741      0.000       0.678       0.722
==============================================================================


## Negative Binomial¶

The negative binomial model gives slightly different results.

[14]:

mod_nbin = sm.NegativeBinomial(rand_data.endog, rand_exog)
res_nbin = mod_nbin.fit(disp=False)
print(res_nbin.summary())

/opt/hostedtoolcache/Python/3.10.14/x64/lib/python3.10/site-packages/statsmodels/base/model.py:607: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals
warnings.warn("Maximum Likelihood optimization failed to "

                     NegativeBinomial Regression Results
==============================================================================
Dep. Variable:                  mdvis   No. Observations:                20190
Model:               NegativeBinomial   Df Residuals:                    20180
Method:                           MLE   Df Model:                            9
Date:                Sat, 25 May 2024   Pseudo R-squ.:                 0.01845
Time:                        15:58:00   Log-Likelihood:                -43384.
converged:                      False   LL-Null:                       -44199.
Covariance Type:            nonrobust   LLR p-value:                     0.000
==============================================================================
coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
lncoins       -0.0579      0.006     -9.515      0.000      -0.070      -0.046
idp           -0.2678      0.023    -11.802      0.000      -0.312      -0.223
lpi            0.0412      0.004      9.938      0.000       0.033       0.049
fmde          -0.0381      0.003    -11.216      0.000      -0.045      -0.031
physlm         0.2691      0.030      8.985      0.000       0.210       0.328
disea          0.0382      0.001     26.080      0.000       0.035       0.041
hlthg         -0.0441      0.020     -2.201      0.028      -0.083      -0.005
hlthf          0.0173      0.036      0.478      0.632      -0.054       0.088
hlthp          0.1782      0.074      2.399      0.016       0.033       0.324
const          0.6635      0.025     26.786      0.000       0.615       0.712
alpha          1.2930      0.019     69.477      0.000       1.256       1.329
==============================================================================


## Alternative solvers¶

The default method for fitting discrete data MLE models is Newton-Raphson. You can use other solvers by using the method argument:

[15]:

mlogit_res = mlogit_mod.fit(method="bfgs", maxiter=250)
print(mlogit_res.summary())

Optimization terminated successfully.
Current function value: 1.548647
Iterations: 111
Function evaluations: 117
MNLogit Regression Results
==============================================================================
Dep. Variable:                    PID   No. Observations:                  944
Model:                        MNLogit   Df Residuals:                      908
Method:                           MLE   Df Model:                           30
Date:                Sat, 25 May 2024   Pseudo R-squ.:                  0.1648
Time:                        15:58:02   Log-Likelihood:                -1461.9
converged:                       True   LL-Null:                       -1750.3
Covariance Type:            nonrobust   LLR p-value:                1.822e-102
==============================================================================
PID=1       coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
const         -0.3734      0.630     -0.593      0.553      -1.608       0.861
logpopul      -0.0115      0.034     -0.337      0.736      -0.079       0.056
selfLR         0.2977      0.094      3.180      0.001       0.114       0.481
age           -0.0249      0.007     -3.823      0.000      -0.038      -0.012
educ           0.0825      0.074      1.121      0.262      -0.062       0.227
income         0.0052      0.018      0.295      0.768      -0.029       0.040
------------------------------------------------------------------------------
PID=2       coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
const         -2.2509      0.763     -2.949      0.003      -3.747      -0.755
logpopul      -0.0888      0.039     -2.266      0.023      -0.166      -0.012
selfLR         0.3917      0.108      3.619      0.000       0.180       0.604
age           -0.0229      0.008     -2.893      0.004      -0.038      -0.007
educ           0.1810      0.085      2.123      0.034       0.014       0.348
income         0.0479      0.022      2.149      0.032       0.004       0.092
------------------------------------------------------------------------------
PID=3       coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
const         -3.6656      1.157     -3.169      0.002      -5.932      -1.399
logpopul      -0.1060      0.057     -1.858      0.063      -0.218       0.006
selfLR         0.5734      0.159      3.617      0.000       0.263       0.884
age           -0.0149      0.011     -1.311      0.190      -0.037       0.007
educ          -0.0072      0.126     -0.057      0.955      -0.255       0.240
income         0.0576      0.034      1.713      0.087      -0.008       0.123
------------------------------------------------------------------------------
PID=4       coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
const         -7.6139      0.958     -7.951      0.000      -9.491      -5.737
logpopul      -0.0916      0.044     -2.091      0.037      -0.177      -0.006
selfLR         1.2788      0.129      9.921      0.000       1.026       1.531
age           -0.0087      0.008     -1.031      0.302      -0.025       0.008
educ           0.1998      0.094      2.123      0.034       0.015       0.384
income         0.0845      0.026      3.226      0.001       0.033       0.136
------------------------------------------------------------------------------
PID=5       coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
const         -7.0605      0.844     -8.362      0.000      -8.715      -5.406
logpopul      -0.0933      0.039     -2.371      0.018      -0.170      -0.016
selfLR         1.3470      0.117     11.494      0.000       1.117       1.577
age           -0.0179      0.008     -2.352      0.019      -0.033      -0.003
educ           0.2169      0.085      2.552      0.011       0.050       0.384
income         0.0810      0.023      3.524      0.000       0.036       0.126
------------------------------------------------------------------------------
PID=6       coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
const        -12.1058      1.060    -11.421      0.000     -14.183     -10.028
logpopul      -0.1409      0.042     -3.343      0.001      -0.223      -0.058
selfLR         2.0701      0.143     14.435      0.000       1.789       2.351
age           -0.0094      0.008     -1.160      0.246      -0.025       0.007
educ           0.3219      0.091      3.534      0.000       0.143       0.500
income         0.1089      0.025      4.304      0.000       0.059       0.158
==============================================================================


Last update: May 25, 2024