Multinomial logit model#

This guide presents an introduction to a simple MNL model example using PyCMTensor.

Include all necessary packages, including pandas and numpy:

1import numpy as np
2import pandas as pd
3import pycmtensor as cmt
4from pycmtensor.expressions import Beta
5from pycmtensor.models import MNL
6from pycmtensor.statistics import elasticities

Load the swissmetro.dat dataset and define the Data object. The arguments for Data are as follows:

  • df: the pandas DataFrame object

  • choice: column label in the DataFrame object which corresponds to the dependent choice variable

  • drop: list of conditional statement(s) which is used to drop certain rows

  • autoscale: boolean variable which indicates whether to scale the variables to between -10 and 10

  • autoscale_except: list of column label(s) which to ignore when autoscale is True

  • split: split ratio for training:validation datasets (1.0 indicates all samples are used for training and no validation samples)

1sm = pd.read_csv("swissmetro.dat", sep="\t")
2db = cmt.Data(
3    df=sm,
4    choice="CHOICE",
5    drop=[sm["CHOICE"] == 0],
6    autoscale=True,
7    autoscale_except=["ID", "ORIGIN", "DEST", "CHOICE"],
8    split=0.8,
9)

Model specification#

Describe the model by defining the \(\beta\) variables Beta() and specify the utility as a list of equations. Optionally, include the availability conditions if required.

Variables are referenced by using db["<name_of_column>"]

 1b_cost = Beta("b_cost", 0.0, None, None, 0)
 2b_time = Beta("b_time", 0.0, None, None, 0)
 3asc_train = Beta("asc_train", 0.0, None, None, 0)
 4asc_car = Beta("asc_car", 0.0, None, None, 0)
 5asc_sm = Beta("asc_sm", 0.0, None, None, 1)
 6
 7U_1 = b_cost * db["TRAIN_CO"] + b_time * db["TRAIN_TT"] + asc_train
 8U_2 = b_cost * db["SM_CO"] + b_time * db["SM_TT"] + asc_sm
 9U_3 = b_cost * db["CAR_CO"] + b_time * db["CAR_TT"] + asc_car
10
11# specify the utility function and the availability conditions
12U = [U_1, U_2, U_3]  # utility
13AV = [db["TRAIN_AV"], db["SM_AV"], db["CAR_AV"]]  # availability

Define the model#

Define the MNL model MNL(db=,params=,utility=,av=,**kwargs). Argument options explanation:

  • db: the Data object

  • params: dict of params where the keys are the name of the parameter and the values are the Beta objects. For example, params={"asc_car":asc_car, "asc_train":asc_train}. A simple shortcut is to use locals() to load all locally defined Python variables.

keyword arguments are the model hyperparameters. For example, optimizer=SGD sets the optimizer algorithm to stochastic gradient descent.

1mymodel = MNL(db=db, params=locals(), utility=U, av=AV)

Train the model#

Train the model by calling the function .train(db, **kwargs).

Keyword arguments are the training hyperparameters. Possible options are max_steps:int, patience:int, lr_scheduler:scheduler.Scheduler, batch_size:int. For more information and other possible options, see config.Config.hyperparameters

1mymodel.train(db, max_steps=100)

11:01:26 [INFO] Start (n=8575)
11:01:35 [INFO] End (t=00:00:10, VE=40.224%, LL=-7529.683422816989, S=23)

Output results and statistics#

Beta parameters

1print(mymodel.results.beta_statistics())

              value   std err     t-test   p-value rob. std err rob. t-test  \
asc_car   -0.853707  0.047113 -18.120447       0.0     0.001773 -481.494189   
asc_sm          0.0         -          -         -            -           -   
asc_train -1.845829   0.05049 -36.558574       0.0     0.010834 -170.368325   
b_cost     0.025854  0.019678   1.313796  0.188915      0.00523    4.943254   
b_time      -0.5217  0.053767  -9.702955       0.0     0.003695 -141.184022   

          rob. p-value  
asc_car            0.0  
asc_sm               -  
asc_train          0.0  
b_cost        0.000001  
b_time             0.0

Model statistics

1print(mymodel.results.model_statistics())

                                          value
Number of training samples used          8575.0
Number of validation samples used        2143.0
Init. log likelihood               -8885.386433
Final log likelihood               -7529.683423
Accuracy                                 59.78%
Likelihood ratio test                2711.40602
Rho square                             0.152577
Rho square bar                         0.152014
Akaike Information Criterion       15069.366846
Bayesian Information Criterion     15104.649877
Final gradient norm                    0.003011

Correlation matrix

1print(mymodel.results.model_correlation_matrix())

             b_cost    b_time  asc_train   asc_car
b_cost     1.000000  0.215873   0.234180 -0.019954
b_time     0.215873  1.000000   0.750515  0.814288
asc_train  0.234180  0.750515   1.000000  0.685223
asc_car   -0.019954  0.814288   0.685223  1.000000

Benchmark

1print(mymodel.results.benchmark())

                       value
Seed                     980
Model build time    00:00:13
Model train time    00:00:10
iterations per sec  345.14/s

Predictions

1# predictions
2print(mymodel.predict(db, return_choices=False))
3print(np.unique(mymodel.predict(db), return_counts=True))

[[0.12138601 0.55492529 0.3236887 ]
 [0.10724806 0.5969908  0.29576114]
 [0.10717803 0.60330066 0.28952132]
 ...
 [0.11440961 0.55831256 0.32727783]
 [0.12859771 0.5241943  0.34720798]
 [0.10920566 0.58416641 0.30662793]]
(array([1, 2], dtype=int64), array([2068,   75], dtype=int64))

Elasticities

1print(elasticities(mymodel, db, 0, "TRAIN_TT"))

[-0.06027963 -0.15127238 -0.20175932 ... -0.02805637 -0.05351249
 -0.04710554]