<jupyter_start><jupyter_text>Multinomial Logit This is a step-by-step guide on how to estimate Multinomial Logit models using the `xlogit` package. You can interactively execute the code in this guide by opening it Google Colab using the following link:[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/arteagac/xlogit/blob/master/examples/multinomial_model.ipynb)  Install `xlogit` package Install `xlogit` using `pip` as shown below.<jupyter_code>!pip install xlogit<jupyter_output>Collecting xlogit
  Downloading https://files.pythonhosted.org/packages/60/5f/9bc576d180c366af77bc04e268536e9e34be23c52a520918aa0cb56b438e/xlogit-0.1.3-py3-none-any.whl
Requirement already satisfied: numpy>=1.13.1 in /usr/local/lib/python3.7/dist-packages (from xlogit) (1.19.5)
Requirement already satisfied: scipy>=1.0.0 in /usr/local/lib/python3.7/dist-packages (from xlogit) (1.4.1)
Installing collected packages: xlogit
Successfully installed xlogit-0.1.3<jupyter_text>Route Choice Dataset This dataset contains choices of 151 commuters among three Home-to-work route alternatives. The three alternatives are arterial, rural, and freeway roads. This dataset was taken from Example 13.1 of the book "Statistical and econometric methods for transportation data analysis" [(Washintong et. al., 2011) ](https://engineering.purdue.edu/~flm/StatEconBook.htm).  Read data We start by importing the data using pandas and renaming the columns of interest (choice, distance, male, and vehicle model). In addition, we create a column with the name of the alternatives and a column that uniquely identifies every observation in the dataset. Note that this dataset is long format.<jupyter_code>import pandas as pd
import numpy as np
df = pd.read_csv("https://engineering.purdue.edu/~flm/StatEcon-Files/Ex13-1.txt",
                 sep="\t", header=None, prefix="x")
df.rename(columns={'x0': 'choice', 'x6': 'dist', 'x10': 'male', 'x14': 'vehmodel'},
          inplace=True)  # Rename columns of interest
df['alt'] = np.tile(['arterial', 'rural', 'freeway'], len(df)//3)  # Add column with alternatives
df['ids'] = np.repeat(np.arange(len(df)//3), 3)  # Add column with unique ids
df['vehage'] = 86 - df['vehmodel']
df<jupyter_output><empty_output><jupyter_text>Create model specification The following code creates the model specification by including additional variables in the dataframe to accomodate the specification needs. The newly added variables are simply the product of existing variables with dummy variables for the different alternatives. Note that the specification in the code below corresponds to the following utility maximization formulation: \begin{equation}\begin{array}{}V_{arterial} & = & \quad \beta_3DIST_{arterial} \\V_{rural} & = \beta_1ASC_{rural} & +  \beta_4DIST_{rural} & + \beta_6VEHAGE_{rural}  \\V_{freeway} & = \beta_2ASC_{freew} & +  \beta_5DIST_{freeway} & +  \beta_7{VEHAGE_{freeway}} & + \beta_8MALE_{freeway} \end{array}\end{equation}<jupyter_code># Alternative specific constants
df['asc_rural'] = np.ones(len(df)) * (df['alt'] == 'rural')
df['asc_freeway'] = np.ones(len(df)) * (df['alt'] == 'freeway')

# Distance
df['dist_arterial'] = df['dist'] * (df['alt'] == 'arterial')
df['dist_rural'] = df['dist'] * (df['alt'] == 'rural')
df['dist_freeway'] = df['dist'] * (df['alt'] == 'freeway')

# Vehicle age
df['vehage_rural'] = df['vehage'] * (df['alt'] == 'rural')
df['vehage_freeway'] = df['vehage'] * (df['alt'] == 'freeway')

# Male driver
df['male_freeway'] = df['male'] * (df['alt'] == 'freeway')<jupyter_output><empty_output><jupyter_text>Estimate model After creating the model specification, we can use `xlogit` to estimate the model as follows:<jupyter_code>from xlogit import MultinomialLogit
varnames=['asc_rural', 'asc_freeway', 'dist_arter', 'dist_rural',
          'dist_freew', 'vehage_rural', 'vehage_freew', 'male_freew']
model = MultinomialLogit()
model.fit(X=df[varnames], y=df['choice'], varnames=varnames,
          ids=df['ids'], alts=df['alt'])
model.summary()<jupyter_output>Estimation time= 0.0 seconds
---------------------------------------------------------------------------
Coefficient              Estimate      Std.Err.         z-val         P>|z|
---------------------------------------------------------------------------
asc_rural               2.8131275     1.1504189     2.4453072        0.0416 *  
asc_freeway            -2.6868817     2.2817329    -1.1775619         0.398    
dist_arterial          -0.1229102     0.0240053    -5.1201358      4.14e-06 ***
dist_rural             -0.1773579     0.0279818    -6.3383224       1.3e-08 ***
dist_freeway           -0.0956391     0.0369681    -2.5870738        0.0295 *  
vehage_rural            0.1236721     0.0535597     2.3090535         0.057 .  
vehage_freeway          0.2268642     0.0755401     3.0032267       0.00969 ** 
male_freeway            0.5990000     0.6202114     0.9657996         0.499    
---------------------------------------------------------------------------
Significance:  0 '***' 0.00[...]<jupyter_text>Swissmetro Dataset In this example, we will estimate a Multinomial Logit where each alternative is defined with a different utility specification. The swissmetro dataset is an SP/RP survey dataset popularly used in Biogeme and Pylogit examples. The dataset is available at http://transp-or.epfl.ch/data/swissmetro.dat and [Bierlaire et. al., (2001)](https://transp-or.epfl.ch/documents/proceedings/BierAxhaAbay01.pdf) provides a detailed discussion of the data as wells as its context and collection process. . Note that the dataset is available in wide format; therefore, we need to convert it to long format for `xlogit`.  Read data The dataset is imported to the Python environment using `pandas`. Then, two types of samples, ones with a trip purpose different to commute or business and ones with an unknown choice, are filtered out. The original dataset contains 10,729 records, but after filtering, 6,768 records remain for following analysis. Finally, a new column that uniquely identifies each sample is added to the dataframe and the `CHOICE` column, which originally contains a numerical coding of the choices, is mapped to a description that is consistent with the alternatives in the column names.<jupyter_code>import pandas as pd
import numpy as np

df_wide = pd.read_table("http://transp-or.epfl.ch/data/swissmetro.dat", sep='\t')

# Keep only observations for commute and business purposes that contain known choices
df_wide = df_wide[(df_wide['PURPOSE'].isin([1, 3]) & (df_wide['CHOICE'] != 0))]

df_wide['custom_id'] = np.arange(len(df_wide))  # Add unique identifier
df_wide['CHOICE'] = df_wide['CHOICE'].map({1: 'TRAIN', 2:'SM', 3: 'CAR'})
df_wide<jupyter_output><empty_output><jupyter_text>Reshape data Given that `xlogit` requires the dataset to be provided in the long format, we reshape the dataset using the `wide_to_long` utility provided by `xlogit`. This function takes as input the column that uniquely identifies each sample (`id_col`), the name of column to save the alternatives (`alt_name`), the list of alternatives (`alt_list`), the columns that vary across alternatives (`varying`), and whether the alternative names are prefix in the column names (`alt_is_prefix`).The `wide_to_long` method fills with `NaN` the columns that do not have certain alternatives (e.g. `SEATS` and `HE`). Depending on your specification needs, you can ignore the `NaN` or replace them with zeros. In this case we replaced them with zeros using the `empty_val` parameter. Additional details about the `wide_to_long` function can be found in the [xlogit's documentation](https://xlogit.readthedocs.io/en/latest/api/utils.html).<jupyter_code>from xlogit.utils import wide_to_long

df = wide_to_long(df_wide, id_col='custom_id', alt_name='alt', sep='_',
                  alt_list=['TRAIN', 'SM', 'CAR'], empty_val=0,
                  varying=['TT', 'CO', 'HE', 'AV', 'SEATS'], alt_is_prefix=True)
df<jupyter_output><empty_output><jupyter_text>We then scale some variables as per the [examples in Biogeme and Pylogit](https://https://github.com/timothyb0912/pylogit/blob/master/examples/notebooks/Main%20PyLogit%20Example.ipynb). The time and headway variables are converted to hours, the price is scaled, and the new variables `single_luggage`, `free ticket`, `multiple_luggage`, `regular_class` and `train_survey` are created to accommodate the model specification as shown below:<jupyter_code># Scale travel time and headway variables to hours
df['time'] = df['TT'] / 60.0 
df['headway'] = df['HE'] / 60.0 
df['cost'] = df['CO'] / 100.0

# We set the cost as zero for individuals with an annual pass paid by employer
annual_pass = (df['alt'].isin(['TRAIN', 'SM'])) & ((df["GA"] == 1) | (df["WHO"] == 2))
df["cost"] = df["cost"] * (~annual_pass)

#Travellers carrying only single luggage
df["single_luggage"] = (df["LUGGAGE"] == 1).astype(int)

#Travellers carrying more than one luggage
df["multiple_luggage"] = (df["LUGGAGE"] == 3).astype(int)

# Travellers travelling in classes other than First class
df["regular_class"] = 1 - df["FIRST"]

# Travellers who responded to the survey while on a train
df["train_survey"] = 1 - df["SURVEY"]
df<jupyter_output><empty_output><jupyter_text>Create model specification By operating the dataframe columns, highly-flexible utility specifications can be modeled in `xlogit`. As shown below, alternative specific constants or coefficients can be included in the specification by strategically creating new columns and setting their values depending on the alternative. This flexibility allows even the specification of one or multiple coefficients per alternative or group of alternatives.<jupyter_code># Create model specification
# Alternative Specific Constants
df['asc_train'] = np.ones(len(df))*(df['alt'] == 'TRAIN')
df['asc_sm'] = np.ones(len(df))*(df['alt'] == 'SM')

# Travel cost (One coefficient per alternative)
df['cost_train'] = df['cost']*(df['alt'] == 'TRAIN')
df['cost_sm'] = df['cost']*(df['alt'] == 'SM')
df['cost_car'] = df['cost']*(df['alt'] == 'CAR')

# Travel time (One coefficient for train and sm and other for car)
df['time_train_sm'] = df['time']*((df['alt'] == 'TRAIN') | (df['alt'] == 'SM'))
df['time_car'] = df['time']*(df['alt'] == 'CAR')

# Headway (One coefficient per alternative, except for car)
df['headway_train'] = df['headway']*(df['alt'] == 'TRAIN')
df['headway_sm'] = df['headway']*(df['alt'] == 'SM')

# Seat config (Coefficient only for swissmetro)
df['seatconf_sm'] = df['SEATS']*(df['alt'] == 'SM')

# Train Survey (Coefficient only for swissmetro)
df['survey_train_sm'] = df['train_survey']* ((df['alt'] == 'TRAIN') | (df['alt'] == 'SM'))

# Regular class (Coefficient only for swissmetro)
df['regular_class_sm'] = df['regular_class']*(df['alt'] == 'TRAIN')

# Luggage (Coefficient only for car)
df['single_lug_car'] = df['single_luggage']*(df['alt'] == 'CAR')
df['multip_lug_car'] = df['multiple_luggage']*(df['alt'] == 'CAR')<jupyter_output><empty_output><jupyter_text>Estimate model The swissmetro dataset contains unbalanced choice situations across individuals (i.e., some individuals do not have observations for all alternatives). The `avail` option enables estimation for such datasets. `avail` takes the values that indicate the availability of each alternative across individuals.Once the model specification is complete, the model is estimated as follows:<jupyter_code>from xlogit import MultinomialLogit

varnames=['asc_train', 'asc_sm', 'time_train_sm', 'time_car', 'cost_train',
          'cost_sm', 'cost_car', 'headway_train', 'headway_sm', 'seatconf_sm',
          'survey_train_sm', 'regular_class_sm', 'single_lug_car',
          'multip_lug_car']
model = MultinomialLogit()
model.fit(X=df[varnames],
          y=df['CHOICE'],
          varnames=varnames,
          alts=df['alt'],
          ids=df['custom_id'],
          avail=df['AV'])
model.summary()<jupyter_output>Estimation time= 0.1 seconds
---------------------------------------------------------------------------
Coefficient              Estimate      Std.Err.         z-val         P>|z|
---------------------------------------------------------------------------
asc_train              -1.2929512     0.1237556   -10.4476139      2.43e-24 ***
asc_sm                 -0.5026152     0.1032927    -4.8659312      5.87e-06 ***
time_train_sm          -0.6990098     0.0396510   -17.6290608      8.13e-67 ***
time_car               -0.7229887     0.0442625   -16.3340968      1.18e-57 ***
cost_train             -0.5618773     0.0807075    -6.9618968      2.59e-11 ***
cost_sm                -0.2816843     0.0417373    -6.7489902       1.1e-10 ***
cost_car               -0.5139009     0.0970406    -5.2957304      6.66e-07 ***
headway_train          -0.3143519     0.0505955    -6.2130370      3.49e-09 ***
headway_sm             -0.3773753     0.1652542    -2.2836046        0.0589 .  
seatconf_sm            [...]<jupyter_text>The estimates are identical to those provided by PyLogit (Brathwaite and Walker 2018) as shown in [this link](https://github.com/timothyb0912/pylogit/blob/master/examples/notebooks/Main%20PyLogit%20Example.ipynb)   Fishing Dataset The following example illustrates the estimation of a Multinomial Logit model for choices of 1,182 individuals for sport fishing modes using `xlogit`. The goal is to analyze the market shares of four alternatives (i.e., beach, pier, boat, and charter) based on their cost and fish catch rate. [Cameron (2005)](http://cameron.econ.ucdavis.edu/mmabook/mma.html) provides additional details about this dataset. The following code illustrates how to use `xlogit` to estimate the model parameters.  Read data The data to be analyzed can be imported to Python using any preferred method. In this example, the data in CSV format was imported using the popular `pandas` Python package. However, it is worth highlighting that `xlogit` does not depend on the `pandas` package, as `xlogit` can take any array-like structure as input. This represents an additional advantage because `xlogit` can be used with any preferred dataframe library, and not only with `pandas`.<jupyter_code>import pandas as pd
df = pd.read_csv("https://raw.github.com/arteagac/xlogit/master/examples/data/fishing_long.csv")
df<jupyter_output><empty_output><jupyter_text>Estimate the model Once the data is in the `Python` environment, `xlogit` can be used to fit the model, as shown below. The `MultinomialLogit` class is imported from `xlogit`, and its constructor is used to initialize a new model. The `fit` method estimates the model using the input data and estimation criteria provided as arguments to the method's call. The arguments of the `fit` methods are described in [`xlogit`'s documentation](https://https://xlogit.readthedocs.io/en/latest/api/).<jupyter_code>from xlogit import MultinomialLogit

varnames = ['income','price', 'catch']
model = MultinomialLogit()
model.fit(X=df[varnames],
          y=df['choice'],
          varnames=varnames,
          isvars=['income'],
          ids=df['id'],
          alts=df['alt'],
          fit_intercept=True)
model.summary()<jupyter_output>Estimation time= 0.0 seconds
---------------------------------------------------------------------------
Coefficient              Estimate      Std.Err.         z-val         P>|z|
---------------------------------------------------------------------------
_intercept.boat         0.5273413     0.2017396     2.6139700        0.0264 *  
_intercept.charter      1.6943827     0.2096186     8.0831712       1.2e-14 ***
_intercept.pier         0.7779899     0.2051062     3.7931076      0.000622 ***
income.boat             0.0000894     0.0000473     1.8906643         0.134    
income.charter         -0.0000333     0.0000480    -0.6936462         0.627    
income.pier            -0.0001276     0.0000467    -2.7335425        0.0192 *  
price                  -0.0251161     0.0015240   -16.4800784      5.89e-54 ***
catch                   0.3578374     0.0985344     3.6315992       0.00113 ** 
---------------------------------------------------------------------------
Significance:  0 '***' 0.00[...]<jupyter_text>Heating Dataset For this example, we use the Heating dataset from R's mlogit package (Croissant 2020), which contains choice of heating systems in California house. The dataset contains 90 observations, with 8 explanatory variables and more information can be found in https://cran.r-project.org/web/packages/mlogit/vignettes/e1mlogit.html.  Read data<jupyter_code>import pandas as pd
import numpy as np

df_wide = pd.read_csv("https://raw.github.com/arteagac/xlogit/master/examples/data/heating_wide.csv")
df_wide<jupyter_output><empty_output><jupyter_text>Reshape data The dataset is available in wide format. Since `xlogit` requires the data in long format, we convert it as shown below:<jupyter_code>from xlogit.utils import wide_to_long
df = wide_to_long(df_wide, id_col='idcase', alt_name='alt', varying=['ic', 'oc'],
                  alt_list=['ec', 'er', 'gc', 'gr', 'hp'], sep='.', empty_val=0)
df<jupyter_output><empty_output><jupyter_text>Estimate the model We now import `MultinomialLogit` from xlogit and estimate the model as shown below:<jupyter_code>from xlogit import MultinomialLogit

varnames = ['ic', 'oc']
model = MultinomialLogit()
model.fit(X=df[varnames],
          y=df['depvar'],
          varnames=varnames,
          alts=df['alt'],
          ids=df['idcase'])
model.summary()<jupyter_output>Estimation time= 0.0 seconds
---------------------------------------------------------------------------
Coefficient              Estimate      Std.Err.         z-val         P>|z|
---------------------------------------------------------------------------
ic                     -0.0062318     0.0003516   -17.7222802      2.16e-59 ***
oc                     -0.0045800     0.0003208   -14.2767999      9.14e-41 ***
---------------------------------------------------------------------------
Significance:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Log-Likelihood= -1095.237
AIC= 2194.474
BIC= 2204.079
{'fraction_non_alphanumeric': 0.16405001928055968, 'fraction_numerical': 0.055638186525643146, 'mean_word_length': 3.1599450045829514, 'pattern_counts': {'":': 0, '<': 52, '<?xml version=': 0, '>': 58, 'https://': 15, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '14', 'dir': 'jupyter-structured-clean-dedup', 'id': '368759', 'license': 'MIT', 'n_tokens_mistral': 6248, 'n_tokens_neox': 5565, 'n_words': 1655, 'path': 'examples/multinomial_model.ipynb', 'repo_name': 'arteagac/pymlogit', 'size': '18153'}