Skip to content

Mixed logit

This example shows a multinomial logit model with one random Beta variable for time. We use a "random coefficients" method of defining the model and a simulation based estimation: 

Beta("b_time") + Beta("s_time") * RandomDraws("rnd_time", "lognormal", 100)

The arguments for RandomDraws(name, draw_type, draws) are:

  • name (str) the name of the draw
  • draw type (str) the distribution of the draws
  • draws (int) the number of random draws to sample from
#import packages
import pandas as pd
from pycmtensor.dataset import Dataset
from pycmtensor.expressions import Beta, RandomDraws
from pycmtensor.models import MNL


# load data
lpmc = pd.read_csv("../../data/lpmc.dat", sep="\t")
ds = Dataset(df=lpmc, choice="travel_mode")
ds.split(0.8)


# Beta parameters
asc_walk = Beta("asc_walk", 0.0, None, None, 1)
asc_cycle = Beta("asc_cycle", 0.0, None, None, 0)
asc_pt = Beta("asc_pt", 0.0, None, None, 0)
asc_drive = Beta("asc_drive", 0.0, None, None, 0)
b_cost = Beta("b_cost", 0.0, None, None, 0)
b_time = Beta("b_time", 0.0, None, None, 0)
s_time = Beta("s_time", 0.5, None, None, 0)
b_purpose = Beta("b_purpose", 0.0, None, None, 0)
b_licence = Beta("b_licence", 0.0, None, None, 0)
b_car_own = Beta("b_car_own", 0.0, None, None, 0)


# additional parameter for variance of b_time
s_time = Beta("s_time", 0.5, None, None, 0)


# define the sampling variable
rnd_time = RandomDraws('rnd_time', 'lognormal', 100)


# composite variables
dur_pt = ds["dur_pt_rail"] + ds["dur_pt_bus"] + ds["dur_pt_int"]
cost_drive = ds["cost_driving_fuel"] + ds["cost_driving_ccharge"]


# define the random beta time coefficient
b_time_rnd = b_time + s_time * rnd_time


# utility functions
U_walk = asc_walk + b_time_rnd * ds["dur_walking"]
U_cycle = asc_cycle + b_time_rnd * ds["dur_cycling"] 
U_pt = asc_pt + b_time_rnd * dur_pt + b_cost * ds["cost_transit"]
U_drive = (asc_drive + b_time_rnd * ds["dur_driving"] + b_cost * cost_drive 
        + b_licence * ds["driving_license"] + b_purpose * ds["purpose"] 
        + b_car_own * ds["car_ownership"])

U = [U_walk, U_cycle, U_pt, U_drive]


# the choice model is MNL
mymodel = MNL(ds, locals(), U)


# estimate the model
from pycmtensor.models import train
train(mymodel, ds)

The following code displays the results of the model estimation, model correlation matrices, predicitons, and elasticities

print(mymodel.results.beta_statistics())
print(mymodel.results.model_statistics())
print(mymodel.results.benchmark())

# correlation matrix
print(mymodel.results.model_correlation_matrix())
print(mymodel.results.model_robust_correlation_matrix())

# probability prediction
prob = mymodel.predict(ds)
pd.DataFrame(prob)

# choice prediction
choice = mymodel.predict(ds, =False)
pd.DataFrame(choice)

# elasticity w.r.t. choice 1
elas = mymodel.elasticities(ds, wrt_choice=1)
pd.DataFrame(elas)