<NME> Prediction with PyLogit.ipynb <BEF> ADDFILE <MSG> Added example of using the predict and panel_predict methods with PyLogit. <DFF> @@ -0,0 +1,1776 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Prediction with PyLogit\n", + "The purpose of this notebook is to demonstrate how one can use PyLogit to predict the probabilities of discrete outcomes. The two main use cases are when one has\n", + "1. already estimated a model using PyLogit and wants to use the estimated parameters to predict probabilities for a different dataset\n", + "2. one or more sets of parameters (obtained from some process other than model estimation with PyLogit) and one wants to use those parameters to predict probabilities.\n", + "\n", + "Both of these uses will be demonstrated below. Essentially, the \"predict\" function of a model object is the place to go to perform both of these tasks. The exception to this is for mixed logit models, where in addition to the predict function, there is also a \"panel predict\" function that allows one to condition on past observations for individuals in one's prediction dataset that were also in one's estimation dataset. For more information on the panel_predict method, see chapter 11 (\"Individual-Level Parameters\") of Discrete Choice Methods with Simulation by Kenneth Train (2009).\n", + "\n", + "The examples below should make all of this clear.\n", + "\n", + "The dataset being used in this notebook is the \"Swissmetro\" dataset used in the Python Biogeme examples. The data can be downloaded at <a href=\"http://biogeme.epfl.ch/examples_swissmetro.html\">http://biogeme.epfl.ch/examples_swissmetro.html</a>, and a detailed explanation of the variables and data-collection procedure can be found at http://www.strc.ch/conferences/2001/bierlaire1.pdf." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/Users/timothyb0912/anaconda/lib/python2.7/site-packages/IPython/html.py:14: ShimWarning: The `IPython.html` package has been deprecated. You should import from `notebook` instead. `IPython.html.widgets` has moved to `ipywidgets`.\n", + " \"`IPython.html.widgets` has moved to `ipywidgets`.\", ShimWarning)\n" + ] + } + ], + "source": [ + "# For recording the model specification \n", + "from collections import OrderedDict\n", + "\n", + "# For making plots pretty\n", + "import seaborn\n", + "# For file input/output\n", + "import pandas as pd\n", + "# For vectorized math operations\n", + "import numpy as np\n", + "# For plotting\n", + "import matplotlib.pyplot as plt\n", + "\n", + "# For model estimation and prediction\n", + "import pylogit as pl\n", + "\n", + "# To display plots inline\n", + "%matplotlib inline " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Load the dataset\n", + "Note that because we want to demonstrate the use of the dataset for prediction, we will randomly remove 2,000 of the response in the dataset. After estimating the various models, we will then predict the probabilities for these held-out responses." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "# Load the dataset\n", + "raw_data = pd.read_csv(\"../data/long_swiss_metro_data.csv\")\n", + "\n", + "# Get a list of the choice situation ids that will be used\n", + "# in estimation and in prediction\n", + "all_situation_ids = np.sort(raw_data[\"custom_id\"].unique())\n", + "\n", + "# Set a random seed for reproducibility\n", + "np.random.seed(61)\n", + "\n", + "# Shuffle and split the ids\n", + "np.random.shuffle(all_situation_ids)\n", + "prediction_ids = all_situation_ids[:2000]\n", + "estimation_ids = all_situation_ids[2000:]\n", + "\n", + "# Create the estimation and prediction datasets\n", + "estimation_df = raw_data.loc[raw_data[\"custom_id\"].isin(estimation_ids)].copy()\n", + "prediction_df = raw_data.loc[raw_data[\"custom_id\"].isin(prediction_ids)].copy()\n" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/html": [ + "<div>\n", + "<table border=\"1\" class=\"dataframe\">\n", + " <thead>\n", + " <tr style=\"text-align: right;\">\n", + " <th></th>\n", + " <th>3</th>\n", + " <th>4</th>\n", + " <th>5</th>\n", + " <th>6</th>\n", + " <th>7</th>\n", + " </tr>\n", + " </thead>\n", + " <tbody>\n", + " <tr>\n", + " <th>custom_id</th>\n", + " <td>2.000000</td>\n", + " <td>2.000000</td>\n", + " <td>2.00</td>\n", + " <td>3.000000</td>\n", + " <td>3.000000</td>\n", + " </tr>\n", + " <tr>\n", + " <th>mode_id</th>\n", + " <td>1.000000</td>\n", + " <td>2.000000</td>\n", + " <td>3.00</td>\n", + " <td>1.000000</td>\n", + " <td>2.000000</td>\n", + " </tr>\n", + " <tr>\n", + " <th>CHOICE</th>\n", + " <td>0.000000</td>\n", + " <td>1.000000</td>\n", + " <td>0.00</td>\n", + " <td>0.000000</td>\n", + " <td>1.000000</td>\n", + " </tr>\n", + " <tr>\n", + " <th>GROUP</th>\n", + " <td>2.000000</td>\n", + " <td>2.000000</td>\n", + " <td>2.00</td>\n", + " <td>2.000000</td>\n", + " <td>2.000000</td>\n", + " </tr>\n", + " <tr>\n", + " <th>SURVEY</th>\n", + " <td>0.000000</td>\n", + " <td>0.000000</td>\n", + " <td>0.00</td>\n", + " <td>0.000000</td>\n", + " <td>0.000000</td>\n", + " </tr>\n", + " <tr>\n", + " <th>SP</th>\n", + " <td>1.000000</td>\n", + " <td>1.000000</td>\n", + " <td>1.00</td>\n", + " <td>1.000000</td>\n", + " <td>1.000000</td>\n", + " </tr>\n", + " <tr>\n", + " <th>ID</th>\n", + " <td>1.000000</td>\n", + " <td>1.000000</td>\n", + " <td>1.00</td>\n", + " <td>1.000000</td>\n", + " <td>1.000000</td>\n", + " </tr>\n", + " <tr>\n", + " <th>PURPOSE</th>\n", + " <td>1.000000</td>\n", + " <td>1.000000</td>\n", + " <td>1.00</td>\n", + " <td>1.000000</td>\n", + " <td>1.000000</td>\n", + " </tr>\n", + " <tr>\n", + " <th>FIRST</th>\n", + " <td>0.000000</td>\n", + " <td>0.000000</td>\n", + " <td>0.00</td>\n", + " <td>0.000000</td>\n", + " <td>0.000000</td>\n", + " </tr>\n", + " <tr>\n", + " <th>TICKET</th>\n", + " <td>1.000000</td>\n", + " <td>1.000000</td>\n", + " <td>1.00</td>\n", + " <td>1.000000</td>\n", + " <td>1.000000</td>\n", + " </tr>\n", + " <tr>\n", + " <th>WHO</th>\n", + " <td>1.000000</td>\n", + " <td>1.000000</td>\n", + " <td>1.00</td>\n", + " <td>1.000000</td>\n", + " <td>1.000000</td>\n", + " </tr>\n", + " <tr>\n", + " <th>LUGGAGE</th>\n", + " <td>0.000000</td>\n", + " <td>0.000000</td>\n", + " <td>0.00</td>\n", + " <td>0.000000</td>\n", + " <td>0.000000</td>\n", + " </tr>\n", + " <tr>\n", + " <th>AGE</th>\n", + " <td>3.000000</td>\n", + " <td>3.000000</td>\n", + " <td>3.00</td>\n", + " <td>3.000000</td>\n", + " <td>3.000000</td>\n", + " </tr>\n", + " <tr>\n", + " <th>MALE</th>\n", + " <td>0.000000</td>\n", + " <td>0.000000</td>\n", + " <td>0.00</td>\n", + " <td>0.000000</td>\n", + " <td>0.000000</td>\n", + " </tr>\n", + " <tr>\n", + " <th>INCOME</th>\n", + " <td>2.000000</td>\n", + " <td>2.000000</td>\n", + " <td>2.00</td>\n", + " <td>2.000000</td>\n", + " <td>2.000000</td>\n", + " </tr>\n", + " <tr>\n", + " <th>GA</th>\n", + " <td>0.000000</td>\n", + " <td>0.000000</td>\n", + " <td>0.00</td>\n", + " <td>0.000000</td>\n", + " <td>0.000000</td>\n", + " </tr>\n", + " <tr>\n", + " <th>ORIGIN</th>\n", + " <td>2.000000</td>\n", + " <td>2.000000</td>\n", + " <td>2.00</td>\n", + " <td>2.000000</td>\n", + " <td>2.000000</td>\n", + " </tr>\n", + " <tr>\n", + " <th>DEST</th>\n", + " <td>1.000000</td>\n", + " <td>1.000000</td>\n", + " <td>1.00</td>\n", + " <td>1.000000</td>\n", + " <td>1.000000</td>\n", + " </tr>\n", + " <tr>\n", + " <th>seat_configuration</th>\n", + " <td>0.000000</td>\n", + " <td>0.000000</td>\n", + " <td>0.00</td>\n", + " <td>0.000000</td>\n", + " <td>0.000000</td>\n", + " </tr>\n", + " <tr>\n", + " <th>travel_time</th>\n", + " <td>103.000000</td>\n", + " <td>60.000000</td>\n", + " <td>117.00</td>\n", + " <td>130.000000</td>\n", + " <td>67.000000</td>\n", + " </tr>\n", + " <tr>\n", + " <th>headway</th>\n", + " <td>30.000000</td>\n", + " <td>10.000000</td>\n", + " <td>0.00</td>\n", + " <td>60.000000</td>\n", + " <td>30.000000</td>\n", + " </tr>\n", + " <tr>\n", + " <th>travel_cost</th>\n", + " <td>48.000000</td>\n", + " <td>49.000000</td>\n", + " <td>84.00</td>\n", + " <td>48.000000</td>\n", + " <td>58.000000</td>\n", + " </tr>\n", + " <tr>\n", + " <th>travel_time_hrs</th>\n", + " <td>1.716667</td>\n", + " <td>1.000000</td>\n", + " <td>1.95</td>\n", + " <td>2.166667</td>\n", + " <td>1.116667</td>\n", + " </tr>\n", + " <tr>\n", + " <th>headway_hrs</th>\n", + " <td>0.500000</td>\n", + " <td>0.166667</td>\n", + " <td>0.00</td>\n", + " <td>1.000000</td>\n", + " <td>0.500000</td>\n", + " </tr>\n", + " <tr>\n", + " <th>free_ticket</th>\n", + " <td>0.000000</td>\n", + " <td>0.000000</td>\n", + " <td>0.00</td>\n", + " <td>0.000000</td>\n", + " <td>0.000000</td>\n", + " </tr>\n", + " <tr>\n", + " <th>travel_cost_hundreth</th>\n", + " <td>0.480000</td>\n", + " <td>0.490000</td>\n", + " <td>0.84</td>\n", + " <td>0.480000</td>\n", + " <td>0.580000</td>\n", + " </tr>\n", + " <tr>\n", + " <th>single_luggage_piece</th>\n", + " <td>0.000000</td>\n", + " <td>0.000000</td>\n", + " <td>0.00</td>\n", + " <td>0.000000</td>\n", + " <td>0.000000</td>\n", + " </tr>\n", + " <tr>\n", + " <th>multiple_luggage_pieces</th>\n", + " <td>0.000000</td>\n", + " <td>0.000000</td>\n", + " <td>0.00</td>\n", + " <td>0.000000</td>\n", + " <td>0.000000</td>\n", + " </tr>\n", + " <tr>\n", + " <th>regular_class</th>\n", + " <td>1.000000</td>\n", + " <td>1.000000</td>\n", + " <td>1.00</td>\n", + " <td>1.000000</td>\n", + " <td>1.000000</td>\n", + " </tr>\n", + " <tr>\n", + " <th>train_survey</th>\n", + " <td>1.000000</td>\n", + " <td>1.000000</td>\n", + " <td>1.00</td>\n", + " <td>1.000000</td>\n", + " <td>1.000000</td>\n", + " </tr>\n", + " </tbody>\n", + "</table>\n", + "</div>" + ], + "text/plain": [ + " 3 4 5 6 7\n", + "custom_id 2.000000 2.000000 2.00 3.000000 3.000000\n", + "mode_id 1.000000 2.000000 3.00 1.000000 2.000000\n", + "CHOICE 0.000000 1.000000 0.00 0.000000 1.000000\n", + "GROUP 2.000000 2.000000 2.00 2.000000 2.000000\n", + "SURVEY 0.000000 0.000000 0.00 0.000000 0.000000\n", + "SP 1.000000 1.000000 1.00 1.000000 1.000000\n", + "ID 1.000000 1.000000 1.00 1.000000 1.000000\n", + "PURPOSE 1.000000 1.000000 1.00 1.000000 1.000000\n", + "FIRST 0.000000 0.000000 0.00 0.000000 0.000000\n", + "TICKET 1.000000 1.000000 1.00 1.000000 1.000000\n", + "WHO 1.000000 1.000000 1.00 1.000000 1.000000\n", + "LUGGAGE 0.000000 0.000000 0.00 0.000000 0.000000\n", + "AGE 3.000000 3.000000 3.00 3.000000 3.000000\n", + "MALE 0.000000 0.000000 0.00 0.000000 0.000000\n", + "INCOME 2.000000 2.000000 2.00 2.000000 2.000000\n", + "GA 0.000000 0.000000 0.00 0.000000 0.000000\n", + "ORIGIN 2.000000 2.000000 2.00 2.000000 2.000000\n", + "DEST 1.000000 1.000000 1.00 1.000000 1.000000\n", + "seat_configuration 0.000000 0.000000 0.00 0.000000 0.000000\n", + "travel_time 103.000000 60.000000 117.00 130.000000 67.000000\n", + "headway 30.000000 10.000000 0.00 60.000000 30.000000\n", + "travel_cost 48.000000 49.000000 84.00 48.000000 58.000000\n", + "travel_time_hrs 1.716667 1.000000 1.95 2.166667 1.116667\n", + "headway_hrs 0.500000 0.166667 0.00 1.000000 0.500000\n", + "free_ticket 0.000000 0.000000 0.00 0.000000 0.000000\n", + "travel_cost_hundreth 0.480000 0.490000 0.84 0.480000 0.580000\n", + "single_luggage_piece 0.000000 0.000000 0.00 0.000000 0.000000\n", + "multiple_luggage_pieces 0.000000 0.000000 0.00 0.000000 0.000000\n", + "regular_class 1.000000 1.000000 1.00 1.000000 1.000000\n", + "train_survey 1.000000 1.000000 1.00 1.000000 1.000000" + ] + }, + "execution_count": 3, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# Look at the first 5 rows of the data\n", + "# Note that for the mode_id, 1-Train, 2-SwissMetro, 3-Car\n", + "estimation_df.head().T" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Create the model specification\n", + "As in previous examples, the model specification being used is\n", + "$$\n", + "\\begin{aligned}\n", + "V_{i, \\textrm{Train}} &= \\textrm{ASC Train} + \\\\\n", + "&\\quad \\beta _{ \\textrm{tt_transit} } \\textrm{Travel Time} _{ \\textrm{Train}} * \\frac{1}{60} + \\\\\n", + "&\\quad \\beta _{ \\textrm{tc_train} } \\textrm{Travel Cost}_{\\textrm{Train}} * \\left( GA == 0 \\right) * 0.01 + \\\\\n", + "&\\quad \\beta _{ \\textrm{headway_train} } \\textrm{Headway} _{\\textrm{Train}} * \\frac{1}{60} + \\\\\n", + "&\\quad \\beta _{ \\textrm{survey} } \\left( \\textrm{Train Survey} == 1 \\right) \\\\\n", + "\\\\\n", + "V_{i, \\textrm{Swissmetro}} &= \\textrm{ASC Swissmetro} + \\\\\n", + "&\\quad \\beta _{ \\textrm{tt_transit} } \\textrm{Travel Time} _{ \\textrm{Swissmetro}} * \\frac{1}{60} + \\\\\n", + "&\\quad \\beta _{ \\textrm{tc_sm} } \\textrm{Travel Cost}_{\\textrm{Swissmetro}} * \\left( GA == 0 \\right) * 0.01 + \\\\\n", + "&\\quad \\beta _{ \\textrm{headway_sm} } \\textrm{Heaway} _{\\textrm{Swissmetro}} * \\frac{1}{60} + \\\\\n", + "&\\quad \\beta _{ \\textrm{seat} } \\left( \\textrm{Seat Configuration} == 1 \\right) \\\\\n", + "&\\quad \\beta _{ \\textrm{survey} } \\left( \\textrm{Train Survey} == 1 \\right) \\\\\n", + "&\\quad \\beta _{ \\textrm{first_class} } \\left( \\textrm{First Class} == 0 \\right) \\\\\n", + "\\\\\n", + "V_{i, \\textrm{Car}} &= \\beta _{ \\textrm{tt_car} } \\textrm{Travel Time} _{ \\textrm{Car}} * \\frac{1}{60} + \\\\\n", + "&\\quad \\beta _{ \\textrm{tc_car}} \\textrm{Travel Cost}_{\\textrm{Car}} * 0.01 + \\\\\n", + "&\\quad \\beta _{\\textrm{luggage}=1} \\left( \\textrm{Luggage} == 1 \\right) + \\\\\n", + "&\\quad \\beta _{\\textrm{luggage}>1} \\left( \\textrm{Luggage} > 1 \\right)\n", + "\\end{aligned}\n", + "$$\n", + "Here, $i$ is used to denote a particular individual." + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "# NOTE: - Specification and variable names must be ordered dictionaries.\n", + "# - Keys should be variables within the long format dataframe.\n", + "# The sole exception to this is the \"intercept\" key.\n", + "# - For the specification dictionary, the values should be lists\n", + "# of integers or or lists of lists of integers. Within a list, \n", + "# or within the inner-most list, the integers should be the \n", + "# alternative ID's of the alternative whose utility specification \n", + "# the explanatory variable is entering. Lists of lists denote \n", + "# alternatives that will share a common coefficient for the variable\n", + "# in question.\n", + "\n", + "basic_specification = OrderedDict()\n", + "basic_names = OrderedDict()\n", + "\n", + "basic_specification[\"intercept\"] = [1, 2]\n", + "basic_names[\"intercept\"] = ['ASC Train',\n", + " 'ASC Swissmetro']\n", + "\n", + "basic_specification[\"travel_time_hrs\"] = [[1, 2,], 3]\n", + "basic_names[\"travel_time_hrs\"] = ['Travel Time, units:hrs (Train and Swissmetro)',\n", + " 'Travel Time, units:hrs (Car)']\n", + "\n", + "basic_specification[\"travel_cost_hundreth\"] = [1, 2, 3]\n", + "basic_names[\"travel_cost_hundreth\"] = ['Travel Cost * (Annual Pass == 0), units: 0.01 CHF (Train)',\n", + " 'Travel Cost * (Annual Pass == 0), units: 0.01 CHF (Swissmetro)',\n", + " 'Travel Cost, units: 0.01 CHF (Car)']\n", + "\n", + "basic_specification[\"headway_hrs\"] = [1, 2]\n", + "basic_names[\"headway_hrs\"] = [\"Headway, units:hrs, (Train)\",\n", + " \"Headway, units:hrs, (Swissmetro)\"]\n", + "\n", + "basic_specification[\"seat_configuration\"] = [2]\n", + "basic_names[\"seat_configuration\"] = ['Airline Seat Configuration, base=No (Swissmetro)']\n", + "\n", + "basic_specification[\"train_survey\"] = [[1, 2]]\n", + "basic_names[\"train_survey\"] = [\"Surveyed on a Train, base=No, (Train and Swissmetro)\"]\n", + "\n", + "basic_specification[\"regular_class\"] = [1]\n", + "basic_names[\"regular_class\"] = [\"First Class == False, (Swissmetro)\"]\n", + "\n", + "basic_specification[\"single_luggage_piece\"] = [3]\n", + "basic_names[\"single_luggage_piece\"] = [\"Number of Luggage Pieces == 1, (Car)\"]\n", + "\n", + "basic_specification[\"multiple_luggage_pieces\"] = [3]\n", + "basic_names[\"multiple_luggage_pieces\"] = [\"Number of Luggage Pieces > 1, (Car)\"]\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Estimate the various models to be used for prediction\n", + "Note that we will first estimate all models that can be estimated with PyLogit. Then we will demonstrate the use of the prediction function for these models." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "# The 'alt_id_column' is the name of a column to be created in the long-format data\n", + "# It will identify the alternative associated with each row.\n", + "alt_id_column = \"mode_id\"\n", + "\n", + "# The \"obs_id_column\" is a custom id column that ignores the fact that this is a \n", + "# panel/repeated-observations dataset. This column denotes each individual choice\n", + "# situation.\n", + "obs_id_column = \"custom_id\"\n", + "\n", + "# The \"choice_column\" records the name of the column that denotes whether or not each\n", + "# individual chose the alternative on a given row.\n", + "choice_column = \"CHOICE\"" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Estimate the MNL model" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Log-likelihood at zero: -4,906.5129\n", + "Initial Log-likelihood: -4,906.5129\n", + "Estimation Time: 0.15 seconds.\n", + "Final log-likelihood: -3,668.2772\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/Users/timothyb0912/anaconda/lib/python2.7/site-packages/scipy/optimize/_minimize.py:382: RuntimeWarning: Method BFGS does not use Hessian information (hess).\n", + " RuntimeWarning)\n" + ] + } + ], + "source": [ + "# Estimate the multinomial logit model (MNL)\n", + "swissmetro_mnl = pl.create_choice_model(data=estimation_df,\n", + " alt_id_col=alt_id_column,\n", + " obs_id_col=obs_id_column,\n", + " choice_col=choice_column,\n", + " specification=basic_specification,\n", + " model_type=\"MNL\",\n", + " names=basic_names)\n", + "\n", + "# Specify the initial values and method for the optimization.\n", + "swissmetro_mnl.fit_mle(np.zeros(14))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Estimate the Asymmetric Logit Model" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Log-likelihood at zero: -4,906.5129\n", + "Initial Log-likelihood: -3,668.2772\n", + "Estimation Time: 0.42 seconds.\n", + "Final log-likelihood: -3,533.2248\n" + ] + } + ], + "source": [ + "# Create the various specification and name dictionaries \n", + "# for the asymmetric logit model. It is the same as for the MNL\n", + "# except that it excludes the intercept from V_{i, alternative}.\n", + "asym_specification = OrderedDict()\n", + "asym_names = OrderedDict()\n", + "\n", + "for col in basic_specification:\n", + " if col != \"intercept\":\n", + " asym_specification[col] = basic_specification[col]\n", + " asym_names[col] = basic_names[col]\n", + "\n", + "# Get the list of intercept names for the asymmetric logit model\n", + "# This is used to tell PyLogit that we will be using 'outside'\n", + "# intercept parameters.\n", + "# See equation 2 of http://arxiv.org/abs/1606.05900 for more\n", + "# details\n", + "asym_intercept_names = basic_names[\"intercept\"]\n", + "# Specify what alternative is not having its intercept estimated\n", + "asym_intercept_ref_pos = 2\n", + "\n", + "# Give names to the shape parameters of the asymmetric logit model\n", + "# Note that all of the shape parameters are not identifiable so we\n", + "# need to restrict one of them. This accounts for why we do not have\n", + "# a \"shape_car\" name.\n", + "asym_shape_names = [\"shape_train\", \"shape_swiss_metro\"]\n", + "\n", + "# Note the index of the alternative whose shape parameter is constrained\n", + "# (i.e. the Car alternative)\n", + "asym_ref = 2\n", + "\n", + "# Create the asymmetric logit's model object\n", + "swissmetro_asym = pl.create_choice_model(data=estimation_df,\n", + " alt_id_col=alt_id_column,\n", + " obs_id_col=obs_id_column,\n", + " choice_col=choice_column,\n", + " specification=asym_specification,\n", + " model_type=\"Asym\",\n", + " names=asym_names,\n", + " shape_names=asym_shape_names,\n", + " intercept_names=asym_intercept_names,\n", + " shape_ref_pos=asym_ref,\n", + " intercept_ref_pos=asym_intercept_ref_pos)\n", + "\n", + "# Note that, below, we use None for initial values and use kwargs to\n", + "# specify our initial values for the optimization. This is not\n", + "# necessary, but it is convenient and an option.\n", + "\n", + "# Note that dividing the index coefficients by log(J), where J is\n", + "# the number of alternatives in the dataset, accounts for the fact\n", + "# that when each shape parameter is 1/J, the value of the asymmetric\n", + "# logit model's estimated index coefficients are equal to the logit\n", + "# model's estimates, divided by log(J).\n", + "# Again, see http://arxiv.org/abs/1606.05900 for more details.\n", + "\n", + "swissmetro_asym.fit_mle(None,\n", + " init_shapes=np.zeros(2),\n", + " init_intercepts=swissmetro_mnl.params.values[:2],\n", + " init_coefs=swissmetro_mnl.params.values[2:] / np.log(3))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Estimate the Uneven Logit Model" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Log-likelihood at zero: -4,906.5129\n", + "Initial Log-likelihood: -3,668.2772\n", + "Estimation Time: 0.31 seconds.\n", + "Final log-likelihood: -3,511.3976\n" + ] + } + ], + "source": [ + "# Create the various specification and name dictionaries \n", + "# for the uneven logit model. It is the same as for the MNL\n", + "# except that it excludes the intercept from V_{i, alternative}.\n", + "uneven_specification = OrderedDict()\n", + "uneven_names = OrderedDict()\n", + "\n", + "for col in basic_specification:\n", + " if col != \"intercept\":\n", + " uneven_specification[col] = basic_specification[col]\n", + " uneven_names[col] = basic_names[col]\n", + "\n", + "# Get the list of intercept names for the uneven logit model\n", + "# This is used to tell PyLogit that we will be using 'outside'\n", + "# intercept parameters.\n", + "# See equation 2 of http://arxiv.org/abs/1606.05900 for more\n", + "# details\n", + "uneven_intercept_names = basic_names[\"intercept\"]\n", + "# Specify what alternative is not having its intercept estimated\n", + "uneven_intercept_ref_pos = 2\n", + "\n", + "# Specify the names of the uneven logit model's shape parameters\n", + "# Note that we include \"shape_car\" because all of the uneven logit\n", + "# model's shape parameters are identifiable.\n", + "uneven_shape_names = [\"shape_train\", \"shape_swiss_metro\", \"shape_car\"]\n", + "\n", + "swissmetro_uneven = pl.create_choice_model(data=estimation_df,\n", + " alt_id_col=alt_id_column,\n", + " obs_id_col=obs_id_column,\n", + " choice_col=choice_column,\n", + " specification=uneven_specification,\n", + " model_type=\"Uneven\",\n", + " names=uneven_names,\n", + " shape_names=uneven_shape_names,\n", + " intercept_names=uneven_intercept_names,\n", + " intercept_ref_pos=uneven_intercept_ref_pos)\n", + "\n", + "\n", + "# Also, note that we use None for initial values and use kwargs to\n", + "# specify our initial values for the optimaztion. This is necessary\n", + "# to use 'outside' intercept parameters with the model.\n", + "swissmetro_uneven.fit_mle(None,\n", + " init_shapes=np.zeros(3),\n", + " init_intercepts=swissmetro_mnl.params.values[:2],\n", + " init_coefs=swissmetro_mnl.params.values[2:])\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Estimate the Scobit Model" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Log-likelihood at zero: -4,906.5129\n", + "Initial Log-likelihood: -3,668.2772\n", + "Estimation Time: 0.50 seconds.\n", + "Final log-likelihood: -3,505.9938\n" + ] + } + ], + "source": [ + "# Create the various specification and name dictionaries \n", + "# for the scobit model. It is the same as for the MNL\n", + "# except that it excludes the intercept from V_{i, alternative}.\n", + "scobit_specification = OrderedDict()\n", + "scobit_names = OrderedDict()\n", + "\n", + "for col in basic_specification:\n", + " if col != \"intercept\":\n", + " scobit_specification[col] = basic_specification[col]\n", + " scobit_names[col] = basic_names[col]\n", + "\n", + "# Get the list of intercept names for the scobit model\n", + "# This is used to tell PyLogit that we will be using 'outside'\n", + "# intercept parameters.\n", + "# See equation 2 of http://arxiv.org/abs/1606.05900 for more\n", + "# details\n", + "scobit_intercept_names = basic_names[\"intercept\"]\n", + "\n", + "# Create the names of the shape parameters that are needed for the scobit model\n", + "scobit_shape_names = [\"shape_train\", \"shape_swiss_metro\", \"shape_car\"]\n", + "\n", + "# Specify which intercept/ASC is not being estimated (namely, the Car intercept)\n", + "scobit_intercept_ref = 2\n", + "\n", + "swissmetro_scobit = pl.create_choice_model(data=estimation_df,\n", + " alt_id_col=alt_id_column,\n", + " obs_id_col=obs_id_column,\n", + " choice_col=choice_column,\n", + " specification=scobit_specification,\n", + " model_type=\"Scobit\",\n", + " names=scobit_names,\n", + " shape_names=scobit_shape_names,\n", + " intercept_ref_pos=scobit_intercept_ref,\n", + " intercept_names=scobit_intercept_names)\n", + "\n", + "# Note that we are using 'outside' intercept parameters for this model\n", + "swissmetro_scobit.fit_mle(None,\n", + " init_shapes=np.zeros(3),\n", + " init_intercepts=swissmetro_mnl.params.values[:2],\n", + " init_coefs=swissmetro_mnl.params.values[2:])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Estimate the Clog-Log Model" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Log-likelihood at zero: -4,906.5129\n", + "Initial Log-likelihood: -4,249.4904\n", + "Estimation Time: 0.41 seconds.\n", + "Final log-likelihood: -3,680.1916\n" + ] + } + ], + "source": [ + "# Create the specification and name dictionaries \n", + "# for the clog-log model.It is the same as for the MNL\n", + "# except that it excludes the intercept from V_{i, alternative}.\n", + "clog_specification = OrderedDict()\n", + "clog_names = OrderedDict()\n", + "\n", + "# Copy the specification dictionary from the logit\n", + "# model, without the intercept parameters so we\n", + "# can place the intercept parameters outside the\n", + "# index\n", + "for col in basic_specification:\n", + " if col != \"intercept\":\n", + " clog_specification[col] = basic_specification[col]\n", + " clog_names[col] = basic_names[col]\n", + "\n", + "# Get the list of intercept names for the clog-log model\n", + "# This is used to tell PyLogit that we will be using 'outside'\n", + "# intercept parameters.\n", + "# See equation 2 of http://arxiv.org/abs/1606.05900 for more\n", + "# details\n", + "clog_intercept_names = basic_names[\"intercept\"]\n", + "\n", + "# Specify which intercept/ASC is not being estimated \n", + "# (i.e, the Car intercept)\n", + "clog_intercept_ref = 2\n", + "\n", + "# Create the Clog-log model object\n", + "swissmetro_clog = pl.create_choice_model(data=estimation_df,\n", + " alt_id_col=alt_id_column,\n", + " obs_id_col=obs_id_column,\n", + " choice_col=choice_column,\n", + " specification=clog_specification,\n", + " model_type=\"Cloglog\",\n", + " intercept_ref_pos=clog_intercept_ref,\n", + " names=clog_names,\n", + " intercept_names=clog_intercept_names)\n", + "\n", + "# Estimate the clog log model. Note we don't pass one single array of\n", + "# initial values but instead pass keyword arguments\n", + "swissmetro_clog.fit_mle(None,\n", + " init_intercepts=swissmetro_mnl.params.values[:2],\n", + " init_coefs=swissmetro_mnl.params.values[2:])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Estimate a Nested Logit Model\n", + "For the purposes of the example we'll simply nest the \"existing\" modes together and leave the future mode (the SwissMetro) in its own nest." + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Log-likelihood at zero: -4,594.2735\n", + "Initial Log-likelihood: -3,666.8367\n", + "Estimation Time: 0.49 seconds.\n", + "Final log-likelihood: -3,575.2522\n" + ] + }, + { + "data": { + "text/html": [ + "<table class=\"simpletable\">\n", + "<caption>Nested Logit Model Regression Results</caption>\n", + "<tr>\n", + " <th>Dep. Variable:</th> <td>CHOICE</td> <th> No. Observations: </th> <td>4,768</td> \n", + "</tr>\n", + "<tr>\n", + " <th>Model:</th> <td>Nested Logit Model</td> <th> Df Residuals: </th> <td>4,752</td> \n", + "</tr>\n", + "<tr>\n", + " <th>Method:</th> <td>MLE</td> <th> Df Model: </th> <td>16</td> \n", + "</tr>\n", + "<tr>\n", + " <th>Date:</th> <td>Thu, 17 Nov 2016</td> <th> Pseudo R-squ.: </th> <td>0.222</td> \n", + "</tr>\n", + "<tr>\n", + " <th>Time:</th> <td>23:52:23</td> <th> Pseudo R-bar-squ.: </th> <td>0.218</td> \n", + "</tr>\n", + "<tr>\n", + " <th>converged:</th> <td>False</td> <th> Log-Likelihood: </th> <td>-3,575.252</td>\n", + "</tr>\n", + "<tr>\n", + " <th> </th> <td> </td> <th> LL-Null: </th> <td>-4,594.273</td>\n", + "</tr>\n", + "</table>\n", + "<table class=\"simpletable\">\n", + "<tr>\n", + " <td></td> <th>coef</th> <th>std err</th> <th>z</th> <th>P>|z|</th> <th>[95.0% Conf. Int.]</th> \n", + "</tr>\n", + "<tr>\n", + " <th>Future Modes</th> <td> 40.0000</td> <td> nan</td> <td> nan</td> <td> nan</td> <td> nan nan</td>\n", + "</tr>\n", + "<tr>\n", + " <th>Existing Modes</th> <td> -0.4513</td> <td> 0.121</td> <td> -3.728</td> <td> 0.000</td> <td> -0.689 -0.214</td>\n", + "</tr>\n", + "<tr>\n", + " <th>ASC Train</th> <td> -0.6825</td> <td> 0.115</td> <td> -5.924</td> <td> 0.000</td> <td> -0.908 -0.457</td>\n", + "</tr>\n", + "<tr>\n", + " <th>ASC Swissmetro</th> <td> -0.1837</td> <td> 0.113</td> <td> -1.619</td> <td> 0.105</td> <td> -0.406 0.039</td>\n", + "</tr>\n", + "<tr>\n", + " <th>Travel Time, units:hrs (Train and Swissmetro)</th> <td> -0.4452</td> <td> 0.026</td> <td> -17.129</td> <td> 0.000</td> <td> -0.496 -0.394</td>\n", + "</tr>\n", + "<tr>\n", + " <th>Travel Time, units:hrs (Car)</th> <td> -0.3540</td> <td> 0.025</td> <td> -14.339</td> <td> 0.000</td> <td> -0.402 -0.306</td>\n", + "</tr>\n", + "<tr>\n", + " <th>Travel Cost * (Annual Pass == 0), units: 0.01 CHF (Train)</th> <td> -0.1932</td> <td> 0.055</td> <td> -3.505</td> <td> 0.000</td> <td> -0.301 -0.085</td>\n", + "</tr>\n", + "<tr>\n", + " <th>Travel Cost * (Annual Pass == 0), units: 0.01 CHF (Swissmetro)</th> <td> -0.1713</td> <td> 0.049</td> <td> -3.474</td> <td> 0.001</td> <td> -0.268 -0.075</td>\n", + "</tr>\n", + "<tr>\n", + " <th>Travel Cost, units: 0.01 CHF (Car)</th> <td> -0.5397</td> <td> 0.062</td> <td> -8.742</td> <td> 0.000</td> <td> -0.661 -0.419</td>\n", + "</tr>\n", + "<tr>\n", + " <th>Headway, units:hrs, (Train)</th> <td> -0.2679</td> <td> 0.050</td> <td> -5.375</td> <td> 0.000</td> <td> -0.366 -0.170</td>\n", + "</tr>\n", + "<tr>\n", + " <th>Headway, units:hrs, (Swissmetro)</th> <td> -0.4575</td> <td> 0.225</td> <td> -2.031</td> <td> 0.042</td> <td> -0.899 -0.016</td>\n", + "</tr>\n", + "<tr>\n", + " <th>Airline Seat Configuration, base=No (Swissmetro)</th> <td> -0.3986</td> <td> 0.103</td> <td> -3.887</td> <td> 0.000</td> <td> -0.600 -0.198</td>\n", + "</tr>\n", + "<tr>\n", + " <th>Surveyed on a Train, base=No, (Train and Swissmetro)</th> <td> 1.4522</td> <td> 0.112</td> <td> 12.994</td> <td> 0.000</td> <td> 1.233 1.671</td>\n", + "</tr>\n", + "<tr>\n", + " <th>First Class == False, (Swissmetro)</th> <td> 0.2518</td> <td> 0.061</td> <td> 4.101</td> <td> 0.000</td> <td> 0.131 0.372</td>\n", + "</tr>\n", + "<tr>\n", + " <th>Number of Luggage Pieces == 1, (Car)</th> <td> 0.1595</td> <td> 0.058</td> <td> 2.759</td> <td> 0.006</td> <td> 0.046 0.273</td>\n", + "</tr>\n", + "<tr>\n", + " <th>Number of Luggage Pieces > 1, (Car)</th> <td> 0.7496</td> <td> 0.167</td> <td> 4.495</td> <td> 0.000</td> <td> 0.423 1.076</td>\n", + "</tr>\n", + "</table>" + ], + "text/plain": [ + "<class 'statsmodels.iolib.summary.Summary'>\n", + "\"\"\"\n", + " Nested Logit Model Regression Results \n", + "==============================================================================\n", + "Dep. Variable: CHOICE No. Observations: 4,768\n", + "Model: Nested Logit Model Df Residuals: 4,752\n", + "Method: MLE Df Model: 16\n", + "Date: Thu, 17 Nov 2016 Pseudo R-squ.: 0.222\n", + "Time: 23:52:23 Pseudo R-bar-squ.: 0.218\n", + "converged: False Log-Likelihood: -3,575.252\n", + " LL-Null: -4,594.273\n", + "==================================================================================================================================\n", + " coef std err z P>|z| [95.0% Conf. Int.]\n", + "----------------------------------------------------------------------------------------------------------------------------------\n", + "Future Modes 40.0000 nan nan nan nan nan\n", + "Existing Modes -0.4513 0.121 -3.728 0.000 -0.689 -0.214\n", + "ASC Train -0.6825 0.115 -5.924 0.000 -0.908 -0.457\n", + "ASC Swissmetro -0.1837 0.113 -1.619 0.105 -0.406 0.039\n", + "Travel Time, units:hrs (Train and Swissmetro) -0.4452 0.026 -17.129 0.000 -0.496 -0.394\n", + "Travel Time, units:hrs (Car) -0.3540 0.025 -14.339 0.000 -0.402 -0.306\n", + "Travel Cost * (Annual Pass == 0), units: 0.01 CHF (Train) -0.1932 0.055 -3.505 0.000 -0.301 -0.085\n", + "Travel Cost * (Annual Pass == 0), units: 0.01 CHF (Swissmetro) -0.1713 0.049 -3.474 0.001 -0.268 -0.075\n", + "Travel Cost, units: 0.01 CHF (Car) -0.5397 0.062 -8.742 0.000 -0.661 -0.419\n", + "Headway, units:hrs, (Train) -0.2679 0.050 -5.375 0.000 -0.366 -0.170\n", + "Headway, units:hrs, (Swissmetro) -0.4575 0.225 -2.031 0.042 -0.899 -0.016\n", + "Airline Seat Configuration, base=No (Swissmetro) -0.3986 0.103 -3.887 0.000 -0.600 -0.198\n", + "Surveyed on a Train, base=No, (Train and Swissmetro) 1.4522 0.112 12.994 0.000 1.233 1.671\n", + "First Class == False, (Swissmetro) 0.2518 0.061 4.101 0.000 0.131 0.372\n", + "Number of Luggage Pieces == 1, (Car) 0.1595 0.058 2.759 0.006 0.046 0.273\n", + "Number of Luggage Pieces > 1, (Car) 0.7496 0.167 4.495 0.000 0.423 1.076\n", + "==================================================================================================================================\n", + "\"\"\"" + ] + }, + "execution_count": 11, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# Specify the nesting specification.\n", + "# I.e. specify the name of the various nests,\n", + "# and which alternatives belong to which nests.\n", + "nest_membership = OrderedDict()\n", + "nest_membership[\"Future Modes\"] = [2]\n", + "nest_membership[\"Existing Modes\"] = [1, 3]\n", + "\n", + "# Create the nested logit model object\n", + "swissmetro_nested = pl.create_choice_model(data=estimation_df,\n", + " alt_id_col=alt_id_column,\n", + " obs_id_col=obs_id_column,\n", + " choice_col=choice_column,\n", + " specification=basic_specification,\n", + " model_type=\"Nested Logit\",\n", + " names=basic_names,\n", + " nest_spec=nest_membership)\n", + "\n", + "# Create the initial values for the nest parameters\n", + "# Note: This should be in terms of the reparameterized values\n", + "# (i.e. the logit of the inverse of the scale parameter) used\n", + "# by PyLogit.\n", + "# Also, note that 40 corresponds to a scale parameter that is\n", + "# essentially equal to 1, and 5 corresponds to a scale\n", + "# parameter just slightly greater than 1.\n", + "init_nests = np.array([40, 5])\n", + "\n", + "# Create the initial values for the estimation\n", + "initial_values = np.concatenate([init_nests,\n", + " swissmetro_mnl.params.values])\n", + "\n", + "# Estimate the nested logit model\n", + "# Note we constrain the nesting parameter for nest 1 since\n", + "# the Future Modes nest only has 1 alternative in it.\n", + "swissmetro_nested.fit_mle(initial_values,\n", + " constrained_pos=[0])\n", + "\n", + "swissmetro_nested.get_statsmodels_summary()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Estimate a Mixed Logit Model\n", + "Here, we'll estimate a relatively simple model, allowing for normally distributed sensitivities to travel time and travel cost. The distribution will be across the various individuals in the dataset(note, not across the choice situations of a given individual). The individuals are denoted by the \"ID\" column.\n", + "\n", + "Given that we're estimating a distribution for sensitivities to time and cost, it makes sense to use a distribution with bounded support like a log-normal distribution. In this way, we could restrict the sensitivities to only being negative. However, Pylogit does not (yet!) support distributions other than the normal distribution for mixed logit." + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Log-likelihood at zero: -4,906.5129\n", + "Initial Log-likelihood: -3,668.2772\n", + "Estimation Time: 3.25 minutes.\n", + "Final log-likelihood: -2,824.8372\n" + ] + }, + { + "data": { + "text/html": [ + "<table class=\"simpletable\">\n", + "<caption>Mixed Logit Model Regression Results</caption>\n", + "<tr>\n", + " <th>Dep. Variable:</th> <td>CHOICE</td> <th> No. Observations: </th> <td>4,768</td> \n", + "</tr>\n", + "<tr>\n", + " <th>Model:</th> <td>Mixed Logit Model</td> <th> Df Residuals: </th> <td>4,749</td> \n", + "</tr>\n", + "<tr>\n", + " <th>Method:</th> <td>MLE</td> <th> Df Model: </th> <td>19</td> \n", + "</tr>\n", + "<tr>\n", + " <th>Date:</th> <td>Thu, 17 Nov 2016</td> <th> Pseudo R-squ.: </th> <td>0.424</td> \n", + "</tr>\n", + "<tr>\n", + " <th>Time:</th> <td>23:55:44</td> <th> Pseudo R-bar-squ.: </th> <td>0.420</td> \n", + "</tr>\n", + "<tr>\n", + " <th>converged:</th> <td>False</td> <th> Log-Likelihood: </th> <td>-2,824.837</td>\n", + "</tr>\n", + "<tr>\n", + " <th> </th> <td> </td> <th> LL-Null: </th> <td>-4,906.513</td>\n", + "</tr>\n", + "</table>\n", + "<table class=\"simpletable\">\n", + "<tr>\n", + " <td></td> <th>coef</th> <th>std err</th> <th>z</th> <th>P>|z|</th> <th>[95.0% Conf. Int.]</th> \n", + "</tr>\n", + "<tr>\n", + " <th>ASC Train</th> <td> -2.8480</td> <td> 0.460</td> <td> -6.190</td> <td> 0.000</td> <td> -3.750 -1.946</td>\n", + "</tr>\n", + "<tr>\n", + " <th>ASC Swissmetro</th> <td> -2.9022</td> <td> 0.398</td> <td> -7.299</td> <td> 0.000</td> <td> -3.682 -2.123</td>\n", + "</tr>\n", + "<tr>\n", + " <th>Travel Time, units:hrs (Train and Swissmetro)</th> <td> -3.0662</td> <td> 0.178</td> <td> -17.262</td> <td> 0.000</td> <td> -3.414 -2.718</td>\n", + "</tr>\n", + "<tr>\n", + " <th>Travel Time, units:hrs (Car)</th> <td> -3.7078</td> <td> 0.169</td> <td> -21.924</td> <td> 0.000</td> <td> -4.039 -3.376</td>\n", + "</tr>\n", + "<tr>\n", + " <th>Travel Cost * (Annual Pass == 0), units: 0.01 CHF (Train)</th> <td> -5.3458</td> <td> 0.634</td> <td> -8.435</td> <td> 0.000</td> <td> -6.588 -4.104</td>\n", + "</tr>\n", + "<tr>\n", + " <th>Travel Cost * (Annual Pass == 0), units: 0.01 CHF (Swissmetro)</th> <td> -2.7045</td> <td> 0.330</td> <td> -8.195</td> <td> 0.000</td> <td> -3.351 -2.058</td>\n", + "</tr>\n", + "<tr>\n", + " <th>Travel Cost, units: 0.01 CHF (Car)</th> <td> -2.8815</td> <td> 0.333</td> <td> -8.664</td> <td> 0.000</td> <td> -3.533 -2.230</td>\n", + "</tr>\n", + "<tr>\n", + " <th>Headway, units:hrs, (Train)</th> <td> -0.5237</td> <td> 0.117</td> <td> -4.488</td> <td> 0.000</td> <td> -0.752 -0.295</td>\n", + "</tr>\n", + "<tr>\n", + " <th>Headway, units:hrs, (Swissmetro)</th> <td> -0.9966</td> <td> 0.356</td> <td> -2.799</td> <td> 0.005</td> <td> -1.694 -0.299</td>\n", + "</tr>\n", + "<tr>\n", + " <th>Airline Seat Configuration, base=No (Swissmetro)</th> <td> -0.5921</td> <td> 0.185</td> <td> -3.195</td> <td> 0.001</td> <td> -0.955 -0.229</td>\n", + "</tr>\n", + "<tr>\n", + " <th>Surveyed on a Train, base=No, (Train and Swissmetro)</th> <td> 5.8038</td> <td> 0.452</td> <td> 12.852</td> <td> 0.000</td> <td> 4.919 6.689</td>\n", + "</tr>\n", + "<tr>\n", + " <th>First Class == False, (Swissmetro)</th> <td> 0.9976</td> <td> 0.183</td> <td> 5.441</td> <td> 0.000</td> <td> 0.638 1.357</td>\n", + "</tr>\n", + "<tr>\n", + " <th>Number of Luggage Pieces == 1, (Car)</th> <td> 1.1681</td> <td> 0.334</td> <td> 3.503</td> <td> 0.000</td> <td> 0.514 1.822</td>\n", + "</tr>\n", + "<tr>\n", + " <th>Number of Luggage Pieces > 1, (Car)</th> <td> 2.7793</td> <td> 1.030</td> <td> 2.699</td> <td> 0.007</td> <td> 0.761 4.798</td>\n", + "</tr>\n", + "<tr>\n", + " <th>Sigma Travel Time, units:hrs (Train and Swissmetro)</th> <td> 2.2175</td> <td> 0.113</td> <td> 19.692</td> <td> 0.000</td> <td> 1.997 2.438</td>\n", + "</tr>\n", + "<tr>\n", + " <th>Sigma Travel Time, units:hrs (Car)</th> <td> 1.4464</td> <td> 0.108</td> <td> 13.426</td> <td> 0.000</td> <td> 1.235 1.658</td>\n", + "</tr>\n", + "<tr>\n", + " <th>Sigma Travel Cost * (Annual Pass == 0), units: 0.01 CHF (Train)</th> <td> 2.4709</td> <td> 0.512</td> <td> 4.826</td> <td> 0.000</td> <td> 1.467 3.474</td>\n", + "</tr>\n", + "<tr>\n", + " <th>Sigma Travel Cost * (Annual Pass == 0), units: 0.01 CHF (Swissmetro)</th> <td> 1.5809</td> <td> 0.276</td> <td> 5.720</td> <td> 0.000</td> <td> 1.039 2.123</td>\n", + "</tr>\n", + "<tr>\n", + " <th>Sigma Travel Cost, units: 0.01 CHF (Car)</th> <td> -1.8658</td> <td> 0.315</td> <td> -5.930</td> <td> 0.000</td> <td> -2.482 -1.249</td>\n", + "</tr>\n", + "</table>" + ], + "text/plain": [ + "<class 'statsmodels.iolib.summary.Summary'>\n", + "\"\"\"\n", + " Mixed Logit Model Regression Results \n", + "==============================================================================\n", + "Dep. Variable: CHOICE No. Observations: 4,768\n", + "Model: Mixed Logit Model Df Residuals: 4,749\n", + "Method: MLE Df Model: 19\n", + "Date: Thu, 17 Nov 2016 Pseudo R-squ.: 0.424\n", + "Time: 23:55:44 Pseudo R-bar-squ.: 0.420\n", + "converged: False Log-Likelihood: -2,824.837\n", + " LL-Null: -4,906.513\n", + "========================================================================================================================================\n", + " coef std err z P>|z| [95.0% Conf. Int.]\n", + "----------------------------------------------------------------------------------------------------------------------------------------\n", + "ASC Train -2.8480 0.460 -6.190 0.000 -3.750 -1.946\n", + "ASC Swissmetro -2.9022 0.398 -7.299 0.000 -3.682 -2.123\n", + "Travel Time, units:hrs (Train and Swissmetro) -3.0662 0.178 -17.262 0.000 -3.414 -2.718\n", + "Travel Time, units:hrs (Car) -3.7078 0.169 -21.924 0.000 -4.039 -3.376\n", + "Travel Cost * (Annual Pass == 0), units: 0.01 CHF (Train) -5.3458 0.634 -8.435 0.000 -6.588 -4.104\n", + "Travel Cost * (Annual Pass == 0), units: 0.01 CHF (Swissmetro) -2.7045 0.330 -8.195 0.000 -3.351 -2.058\n", + "Travel Cost, units: 0.01 CHF (Car) -2.8815 0.333 -8.664 0.000 -3.533 -2.230\n", + "Headway, units:hrs, (Train) -0.5237 0.117 -4.488 0.000 -0.752 -0.295\n", + "Headway, units:hrs, (Swissmetro) -0.9966 0.356 -2.799 0.005 -1.694 -0.299\n", + "Airline Seat Configuration, base=No (Swissmetro) -0.5921 0.185 -3.195 0.001 -0.955 -0.229\n", + "Surveyed on a Train, base=No, (Train and Swissmetro) 5.8038 0.452 12.852 0.000 4.919 6.689\n", + "First Class == False, (Swissmetro) 0.9976 0.183 5.441 0.000 0.638 1.357\n", + "Number of Luggage Pieces == 1, (Car) 1.1681 0.334 3.503 0.000 0.514 1.822\n", + "Number of Luggage Pieces > 1, (Car) 2.7793 1.030 2.699 0.007 0.761 4.798\n", + "Sigma Travel Time, units:hrs (Train and Swissmetro) 2.2175 0.113 19.692 0.000 1.997 2.438\n", + "Sigma Travel Time, units:hrs (Car) 1.4464 0.108 13.426 0.000 1.235 1.658\n", + "Sigma Travel Cost * (Annual Pass == 0), units: 0.01 CHF (Train) 2.4709 0.512 4.826 0.000 1.467 3.474\n", + "Sigma Travel Cost * (Annual Pass == 0), units: 0.01 CHF (Swissmetro) 1.5809 0.276 5.720 0.000 1.039 2.123\n", + "Sigma Travel Cost, units: 0.01 CHF (Car) -1.8658 0.315 -5.930 0.000 -2.482 -1.249\n", + "========================================================================================================================================\n", + "\"\"\"" + ] + }, + "execution_count": 12, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# Specify the parameters, by name, that are being treated\n", + "# as random variables\n", + "mixing_variables = (basic_names[\"travel_time_hrs\"] +\n", + " basic_names[\"travel_cost_hundreth\"])\n", + "\n", + "# Create the mixed logit model object\n", + "swissmetro_mixed = pl.create_choice_model(data=estimation_df,\n", + " alt_id_col=alt_id_column,\n", + " obs_id_col=obs_id_column,\n", + " choice_col=choice_column,\n", + " specification=basic_specification,\n", + " model_type=\"Mixed Logit\",\n", + " names=basic_names,\n", + " mixing_id_col=\"ID\",\n", + " mixing_vars=mixing_variables)\n", + "\n", + "# Create the initial values for estimation.\n", + "# Note that we need one value for each beta\n", + "# coefficient and for each of the standard deviations\n", + "# that will now be estimated.\n", + "initial_values = np.concatenate([swissmetro_mnl.params.values,\n", + " np.zeros(len(mixing_variables))])\n", + "\n", + "# Estimate the mixed logit model\n", + "swissmetro_mixed.fit_mle(initial_values,\n", + " seed=26,\n", + " num_draws=400)\n", + "\n", + "swissmetro_mixed.get_statsmodels_summary()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Predictions using 1 dataset and 1 set of parameters\n", + "Here, we will simply predict the probabilities of each available alternative being chosen in each choice situation from the prediction dataset, using all of the models that we estimated above. Note that we will perform two predictions with the Mixed Logit model. One will use the standard \"predict\" function and the other will use the \"panel_predict\" function." + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "# Create an ordered dictionary mapping shortened model\n", + "# names to estimated model objects.\n", + "model_name_to_obj = OrderedDict()\n", + "model_name_to_obj[\"Clog\"] = swissmetro_clog\n", + "model_name_to_obj[\"MNL\"] = swissmetro_mnl\n", + "model_name_to_obj[\"Nested\"] = swissmetro_nested\n", + "model_name_to_obj[\"Asym\"] = swissmetro_asym\n", + "model_name_to_obj[\"Uneven\"] = swissmetro_uneven\n", + "model_name_to_obj[\"Scobit\"] = swissmetro_scobit\n", + "model_name_to_obj[\"Mixed\"] = swissmetro_mixed\n", + "\n", + "# Initialize a dictionary to store the log-likelihoods\n", + "# on the prediction dataset.\n", + "prediction_log_likelihoods = OrderedDict()\n", + "for key in model_name_to_obj:\n", + " prediction_log_likelihoods[key] = None\n", + "\n", + "# Calculate and store the predictions for all models besides\n", + "# the mixed logit model. The predictions will be stored on\n", + "# prediction_df\n", + "for model_name in model_name_to_obj:\n", + " if model_name == \"Mixed\":\n", + " # Don't make predictions for the Mixed Logit Model\n", + " # using the for-loop\n", + " continue\n", + " else:\n", + " # Get the model object to be used in making predictions\n", + " model_obj = model_name_to_obj[model_name]\n", + " \n", + " # Note that by default, the predict method returns the\n", + " # predicted probabilities for each available alternative\n", + " # for each choice situation.\n", + " prediction_array = model_obj.predict(prediction_df)\n", + " \n", + " # We can use the \"return_long_probs\" and \"choice_col\" \n", + " # keyword arguments to ensure that we only return the \n", + " # probability of the chosen alternative for each\n", + " # choice situation.\n", + " chosen_predictions = model_obj.predict(prediction_df,\n", + " choice_col=choice_column,\n", + " return_long_probs=False)\n", + " \n", + " # Calculate the predicted log-likelihood in two ways and,\n", + " # for demonstration, ensure that they agree.\n", + " log_likelihood_calc_1 = np.log(chosen_predictions).sum()\n", + " \n", + " log_predictions = np.log(prediction_array)\n", + " choice_array = prediction_df[choice_column].values\n", + " log_likelihood_calc_2 = choice_array.dot(log_predictions)\n", + " \n", + " assert np.allclose(log_likelihood_calc_1,\n", + " log_likelihood_calc_2)\n", + " \n", + " # Store the log-likelihood\n", + " prediction_log_likelihoods[model_name] = log_likelihood_calc_1\n", + " \n", + " # Create a column name for the predictions\n", + " prediction_col = model_name + \"_Predictions\"\n", + " \n", + " # Store the predicted probabilities\n", + " prediction_df[prediction_col] = prediction_array" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "# Calculate and store the predictions for the mixed logit \n", + "# model using both the predict and panel_predict functions.\n", + "# Note that the panel_predict function will condition the\n", + "# predictions on the estimation data of each person\n", + "# (if there is any)\n", + "for prediction_type in [\"Mixed\", \"Mixed_Panel\"]:\n", + " if \"Panel\" not in prediction_type:\n", + " # We set the random seed simply for reproducibility\n", + " predictions = swissmetro_mixed.predict(prediction_df,\n", + " seed=2,\n", + " num_draws=400)\n", + " else:\n", + " # The arguments are ['data', 'num_draws']\n", + " args = [prediction_df, 400]\n", + " kwargs = {\"seed\": 2}\n", + "\n", + " # Note that we can change the return type the same\n", + " # way as with the predict function. Pass\n", + " # return_long_probs = False and \n", + " # choice_col = choice_column as keyword arguments\n", + " predictions = swissmetro_mixed.panel_predict(*args,\n", + " **kwargs)\n", + "\n", + " # Create the column name for the predictions\n", + " prediction_col = prediction_type + \"_Predictions\"\n", + " # Store the predicted probabilities\n", + " prediction_df[prediction_col] = predictions\n", + "\n", + " # Calculate the log-likelihood\n", + " actual_choices = prediction_df[choice_column].values\n", + " log_predicted_probs = np.log(predictions)\n", + "\n", + " prediction_log_likelihoods[prediction_type] =\\\n", + " actual_choices.dot(log_predicted_probs)\n", + " " + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "Mixed -1534.343235\n", + "Clog -1516.255638\n", + "MNL -1495.346433\n", + "Asym -1462.023258\n", + "Nested -1461.652262\n", + "Scobit -1451.833670\n", + "Uneven -1447.384942\n", + "Mixed_Panel -947.548620\n", + "Name: Log-likelihood of Predictions, dtype: float64" + ] + }, + "execution_count": 15, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# Look at the predicted log-likelihoods of the various models\n", + "pd.Series(prediction_log_likelihoods,\n", + " name=\"Log-likelihood of Predictions\").sort(inplace=False)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Predictions using 1 dataset and multiple parameters.\n", + "One feature of PyLogit's predict method is that it allows users to pass sets of parameters as a keyword argument (using \"param_list\"). This allows one to make multiple sets of predictions at once. Such a feature is useful when one has, for instance, a sample of MCMC samples from the posterior distribution of one's model parameters. Alternatively, as we shall see below, this feature is also useful when one wants to plot the log-likelihood function.\n", + "\n", + "To demonstrate how one uses PyLogit to predict probabilities when one has 1 dataset and multiple sets of coefficients, we will simply use the Uneven Logit model. Most models in PyLogit support prediction with multiple sets of coefficients, and all models that support this feature make use of the same \"param_list\" keyword argument in the predict method.\n", + "\n", + "The two times where prediction with multiple parameters <strong><u>is not</u></strong> suppored is in the Nested Logit model and in the \"panel_predict\" method of the Mixed Logit model. Note that the \"standard\" predict method of the Mixed Logit model <strong><u>does</u></strong> support prediction with multiple parameters." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Creating multiple sets of parameters\n", + "Here, we will assume that we want to plot the in-sample (i.e. estimation set) log-likelihood with respect to the travel cost parameter for the car.\n", + "\n", + "First, let's look at the summary for the Uneven Logit model." + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/html": [ + "<table class=\"simpletable\">\n", + "<caption>Multinomial Uneven Logit Model Regression Results</caption>\n", + "<tr>\n", + " <th>Dep. Variable:</th> <td>CHOICE</td> <th> No. Observations: </th> <td>4,768</td> \n", + "</tr>\n", + "<tr>\n", + " <th>Model:</th> <td>Multinomial Uneven Logit Model</td> <th> Df Residuals: </th> <td>4,751</td> \n", + "</tr>\n", + "<tr>\n", + " <th>Method:</th> <td>MLE</td> <th> Df Model: </th> <td>17</td> \n", + "</tr>\n", + "<tr>\n", + " <th>Date:</th> <td>Thu, 17 Nov 2016</td> <th> Pseudo R-squ.: </th> <td>0.284</td> \n", + "</tr>\n", + "<tr>\n", + " <th>Time:</th> <td>23:55:54</td> <th> Pseudo R-bar-squ.: </th> <td>0.281</td> \n", + "</tr>\n", + "<tr>\n", + " <th>converged:</th> <td>False</td> <th> Log-Likelihood: </th> <td>-3,511.398</td>\n", + "</tr>\n", + "<tr>\n", + " <th> </th> <td> </td> <th> LL-Null: </th> <td>-4,906.513</td>\n", + "</tr>\n", + "</table>\n", + "<table class=\"simpletable\">\n", + "<tr>\n", + " <td></td> <th>coef</th> <th>std err</th> <th>z</th> <th>P>|z|</th> <th>[95.0% Conf. Int.]</th> \n", + "</tr>\n", + "<tr>\n", + " <th>shape_train</th> <td> 1.0271</td> <td> 0.141</td> <td> 7.274</td> <td> 0.000</td> <td> 0.750 1.304</td>\n", + "</tr>\n", + "<tr>\n", + " <th>shape_swiss_metro</th> <td> 0.1120</td> <td> 0.232</td> <td> 0.482</td> <td> 0.630</td> <td> -0.343 0.567</td>\n", + "</tr>\n", + "<tr>\n", + " <th>shape_car</th> <td> -0.4157</td> <td> 0.561</td> <td> -0.741</td> <td> 0.458</td> <td> -1.515 0.683</td>\n", + "</tr>\n", + "<tr>\n", + " <th>ASC Train</th> <td> -1.1675</td> <td> 0.240</td> <td> -4.873</td> <td> 0.000</td> <td> -1.637 -0.698</td>\n", + "</tr>\n", + "<tr>\n", + " <th>ASC Swissmetro</th> <td> -0.6523</td> <td> 0.234</td> <td> -2.785</td> <td> 0.005</td> <td> -1.111 -0.193</td>\n", + "</tr>\n", + "<tr>\n", + " <th>Travel Time, units:hrs (Train and Swissmetro)</th> <td> -0.3399</td> <td> 0.051</td> <td> -6.621</td> <td> 0.000</td> <td> -0.441 -0.239</td>\n", + "</tr>\n", + "<tr>\n", + " <th>Travel Time, units:hrs (Car)</th> <td> -0.9527</td> <td> 0.591</td> <td> -1.611</td> <td> 0.107</td> <td> -2.112 0.207</td>\n", + "</tr>\n", + "<tr>\n", + " <th>Travel Cost * (Annual Pass == 0), units: 0.01 CHF (Train)</th> <td> -0.2552</td> <td> 0.071</td> <td> -3.573</td> <td> 0.000</td> <td> -0.395 -0.115</td>\n", + "</tr>\n", + "<tr>\n", + " <th>Travel Cost * (Annual Pass == 0), units: 0.01 CHF (Swissmetro)</th> <td> -0.2218</td> <td> 0.061</td> <td> -3.648</td> <td> 0.000</td> <td> -0.341 -0.103</td>\n", + "</tr>\n", + "<tr>\n", + " <th>Travel Cost, units: 0.01 CHF (Car)</th> <td> -0.5969</td> <td> 0.404</td> <td> -1.479</td> <td> 0.139</td> <td> -1.388 0.194</td>\n", + "</tr>\n", + "<tr>\n", + " <th>Headway, units:hrs, (Train)</th> <td> -0.2833</td> <td> 0.056</td> <td> -5.086</td> <td> 0.000</td> <td> -0.392 -0.174</td>\n", + "</tr>\n", + "<tr>\n", + " <th>Headway, units:hrs, (Swissmetro)</th> <td> -0.4631</td> <td> 0.220</td> <td> -2.102</td> <td> 0.036</td> <td> -0.895 -0.031</td>\n", + "</tr>\n", + "<tr>\n", + " <th>Airline Seat Configuration, base=No (Swissmetro)</th> <td> -0.1097</td> <td> 0.108</td> <td> -1.015</td> <td> 0.310</td> <td> -0.321 0.102</td>\n", + "</tr>\n", + "<tr>\n", + " <th>Surveyed on a Train, base=No, (Train and Swissmetro)</th> <td> 2.0947</td> <td> 0.238</td> <td> 8.809</td> <td> 0.000</td> <td> 1.629 2.561</td>\n", + "</tr>\n", + "<tr>\n", + " <th>First Class == False, (Swissmetro)</th> <td> 0.1665</td> <td> 0.064</td> <td> 2.598</td> <td> 0.009</td> <td> 0.041 0.292</td>\n", + "</tr>\n", + "<tr>\n", + " <th>Number of Luggage Pieces == 1, (Car)</th> <td> 0.5712</td> <td> 0.365</td> <td> 1.565</td> <td> 0.118</td> <td> -0.144 1.286</td>\n", + "</tr>\n", + "<tr>\n", + " <th>Number of Luggage Pieces > 1, (Car)</th> <td> 2.1772</td> <td> 0.967</td> <td> 2.251</td> <td> 0.024</td> <td> 0.281 4.073</td>\n", + "</tr>\n", + "</table>" + ], + "text/plain": [ + "<class 'statsmodels.iolib.summary.Summary'>\n", + "\"\"\"\n", + " Multinomial Uneven Logit Model Regression Results \n", + "==========================================================================================\n", + "Dep. Variable: CHOICE No. Observations: 4,768\n", + "Model: Multinomial Uneven Logit Model Df Residuals: 4,751\n", + "Method: MLE Df Model: 17\n", + "Date: Thu, 17 Nov 2016 Pseudo R-squ.: 0.284\n", + "Time: 23:55:54 Pseudo R-bar-squ.: 0.281\n", + "converged: False Log-Likelihood: -3,511.398\n", + " LL-Null: -4,906.513\n", + "==================================================================================================================================\n", + " coef std err z P>|z| [95.0% Conf. Int.]\n", + "----------------------------------------------------------------------------------------------------------------------------------\n", + "shape_train 1.0271 0.141 7.274 0.000 0.750 1.304\n", + "shape_swiss_metro 0.1120 0.232 0.482 0.630 -0.343 0.567\n", + "shape_car -0.4157 0.561 -0.741 0.458 -1.515 0.683\n", + "ASC Train -1.1675 0.240 -4.873 0.000 -1.637 -0.698\n", + "ASC Swissmetro -0.6523 0.234 -2.785 0.005 -1.111 -0.193\n", + "Travel Time, units:hrs (Train and Swissmetro) -0.3399 0.051 -6.621 0.000 -0.441 -0.239\n", + "Travel Time, units:hrs (Car) -0.9527 0.591 -1.611 0.107 -2.112 0.207\n", + "Travel Cost * (Annual Pass == 0), units: 0.01 CHF (Train) -0.2552 0.071 -3.573 0.000 -0.395 -0.115\n", + "Travel Cost * (Annual Pass == 0), units: 0.01 CHF (Swissmetro) -0.2218 0.061 -3.648 0.000 -0.341 -0.103\n", + "Travel Cost, units: 0.01 CHF (Car) -0.5969 0.404 -1.479 0.139 -1.388 0.194\n", + "Headway, units:hrs, (Train) -0.2833 0.056 -5.086 0.000 -0.392 -0.174\n", + "Headway,
Added example of using the predict and panel_predict methods with PyLogit.