Note
Go to the end to download the full example code.
16. Discrete mixture with panel dataΒΆ
Example of a discrete mixture of logit models, also called latent class model. The class membership model includes socio-economic variables. The datafile is organized as panel data.
Michel Bierlaire, EPFL Mon Jun 23 2025, 16:29:45
from IPython.core.display_functions import display
See the data processing script: Panel data preparation for Swissmetro.
from swissmetro_panel import (
CAR_AV_SP,
CAR_CO_SCALED,
CAR_TT_SCALED,
CHOICE,
INCOME,
SM_AV,
SM_COST_SCALED,
SM_TT_SCALED,
TRAIN_AV_SP,
TRAIN_COST_SCALED,
TRAIN_TT_SCALED,
database,
)
import biogeme.biogeme_logging as blog
from biogeme.biogeme import BIOGEME
from biogeme.expressions import (
Beta,
Draws,
ExpressionOrNumeric,
MonteCarlo,
PanelLikelihoodTrajectory,
log,
)
from biogeme.models import logit
from biogeme.results_processing import (
EstimationResults,
get_pandas_estimated_parameters,
)
logger = blog.get_screen_logger(level=blog.INFO)
logger.info('Example b16_panel_discrete_socio_eco.py')
Example b16_panel_discrete_socio_eco.py
Parameters to be estimated. One version for each latent class.
NUMBER_OF_CLASSES = 2
b_cost = [Beta(f'b_cost_class{i}', 0, None, None, 0) for i in range(NUMBER_OF_CLASSES)]
Define a random parameter, normally distributed across individuals, designed to be used for Monte-Carlo simulation.
b_time = [Beta(f'b_time_class{i}', 0, None, None, 0) for i in range(NUMBER_OF_CLASSES)]
It is advised not to use 0 as starting value for the following parameter.
b_time_s = [
Beta(f'b_time_s_class{i}', 1, None, None, 0) for i in range(NUMBER_OF_CLASSES)
]
b_time_rnd: list[ExpressionOrNumeric] = [
b_time[i] + b_time_s[i] * Draws(f'b_time_rnd_class{i}', 'NORMAL_ANTI')
for i in range(NUMBER_OF_CLASSES)
]
We do the same for the constants, to address serial correlation.
asc_car = [
Beta(f'asc_car_class{i}', 0, None, None, 0) for i in range(NUMBER_OF_CLASSES)
]
asc_car_s = [
Beta(f'asc_car_s_class{i}', 1, None, None, 0) for i in range(NUMBER_OF_CLASSES)
]
asc_car_rnd = [
asc_car[i] + asc_car_s[i] * Draws(f'asc_car_rnd_class{i}', 'NORMAL_ANTI')
for i in range(NUMBER_OF_CLASSES)
]
asc_train = [
Beta(f'asc_train_class{i}', 0, None, None, 0) for i in range(NUMBER_OF_CLASSES)
]
asc_train_s = [
Beta(f'asc_train_s_class{i}', 1, None, None, 0) for i in range(NUMBER_OF_CLASSES)
]
asc_train_rnd = [
asc_train[i] + asc_train_s[i] * Draws(f'asc_train_rnd_class{i}', 'NORMAL_ANTI')
for i in range(NUMBER_OF_CLASSES)
]
asc_sm = [Beta(f'asc_sm_class{i}', 0, None, None, 1) for i in range(NUMBER_OF_CLASSES)]
asc_sm_s = [
Beta(f'asc_sm_s_class{i}', 1, None, None, 0) for i in range(NUMBER_OF_CLASSES)
]
asc_sm_rnd = [
asc_sm[i] + asc_sm_s[i] * Draws(f'asc_sm_rnd_class{i}', 'NORMAL_ANTI')
for i in range(NUMBER_OF_CLASSES)
]
Parameter groups used in the generated reports.
Each group contains the estimated parameters associated with one latent class.
Fixed parameters and parameters that are replaced by identification constraints
are not listed here, as they do not appear in the estimation results. The
parameters of the class membership model are not included in either group, and
will therefore be reported in the automatically generated Other parameters
section.
PARAMETER_GROUPS = {
'Class 0': [
'b_cost_class0',
'asc_car_class0',
'asc_car_s_class0',
'asc_train_class0',
'asc_train_s_class0',
'asc_sm_s_class0',
],
'Class 1': [
'b_cost_class1',
'b_time_class1',
'b_time_s_class1',
'asc_car_class1',
'asc_car_s_class1',
'asc_train_class1',
'asc_train_s_class1',
'asc_sm_s_class1',
],
}
Parameters for the class membership model.
class_cte = Beta('class_cte', 0, None, None, 0)
class_inc = Beta('class_inc', 0, None, None, 0)
In class 0, it is assumed that the time coefficient is zero
b_time_rnd[0] = 0
Utility functions.
v_train_per_class = [
asc_train_rnd[i] + b_time_rnd[i] * TRAIN_TT_SCALED + b_cost[i] * TRAIN_COST_SCALED
for i in range(NUMBER_OF_CLASSES)
]
v_swissmetro_per_class = [
asc_sm_rnd[i] + b_time_rnd[i] * SM_TT_SCALED + b_cost[i] * SM_COST_SCALED
for i in range(NUMBER_OF_CLASSES)
]
v_car_per_class = [
asc_car_rnd[i] + b_time_rnd[i] * CAR_TT_SCALED + b_cost[i] * CAR_CO_SCALED
for i in range(NUMBER_OF_CLASSES)
]
v = [
{1: v_train_per_class[i], 2: v_swissmetro_per_class[i], 3: v_car_per_class[i]}
for i in range(NUMBER_OF_CLASSES)
]
Associate the availability conditions with the alternatives
av = {1: TRAIN_AV_SP, 2: SM_AV, 3: CAR_AV_SP}
The choice model is a discrete mixture of logit, with availability conditions We calculate the conditional probability for each class.
choice_probability_per_class = [
PanelLikelihoodTrajectory(logit(v[i], av, CHOICE)) for i in range(NUMBER_OF_CLASSES)
]
Class membership model.
score_class_0 = class_cte + class_inc * INCOME
prob_class0 = logit({0: score_class_0, 1: 0}, None, 0)
prob_class1 = logit({0: score_class_0, 1: 0}, None, 1)
Conditional on the random variables, likelihood for the individual.
conditional_choice_probability = (
prob_class0 * choice_probability_per_class[0]
+ prob_class1 * choice_probability_per_class[1]
)
We integrate over the random variables using Monte-Carlo
log_probability = log(MonteCarlo(conditional_choice_probability))
The model is complex, and there are numerical issues when calculating the second derivatives. Therefore, we instruct Biogeme not to evaluate the second derivatives. As a consequence, the statistics reported after estimation are based on the BHHH matrix instead of the Rao-Cramer bound.
the_biogeme = BIOGEME(
database,
log_probability,
number_of_draws=5_000,
seed=1223,
calculating_second_derivatives='never',
group_of_parameters=PARAMETER_GROUPS,
)
the_biogeme.model_name = 'b16_panel_discrete_socio_eco'
Biogeme parameters read from biogeme.toml.
Estimate the parameters.
try:
results = EstimationResults.from_yaml_file(
filename=f'saved_results/{the_biogeme.model_name}.yaml'
)
except FileNotFoundError:
results = the_biogeme.estimate()
print(results.short_summary())
Results for model b16_panel_discrete_socio_eco
Nbr of parameters: 16
Sample size: 752
Observations: 6768
Excluded data: 0
Final log likelihood: -3524.399
Akaike Information Criterion: 7080.798
Bayesian Information Criterion: 7154.762
pandas_results = get_pandas_estimated_parameters(
estimation_results=results,
group_of_parameters=PARAMETER_GROUPS,
)
for group_name, pandas_table in pandas_results.items():
display(group_name if group_name else 'Estimated parameters')
display(pandas_table)
Class 0
Name Value BHHH std err. BHHH t-stat. BHHH p-value
2 asc_train_class0 -0.957442 0.677489 -1.413220 0.157591
3 asc_train_s_class0 3.041669 0.640831 4.746447 0.000002
4 b_cost_class0 -1.174028 0.549083 -2.138161 0.032504
5 asc_sm_s_class0 0.071794 6.145112 0.011683 0.990678
6 asc_car_class0 -4.869044 1.646061 -2.957998 0.003096
7 asc_car_s_class0 6.155319 1.883396 3.268203 0.001082
Class 1
Name Value BHHH std err. BHHH t-stat. BHHH p-value
8 asc_train_class1 -0.351554 0.287497 -1.222810 2.214016e-01
9 asc_train_s_class1 1.754333 0.458415 3.826956 1.297378e-04
10 b_time_class1 -6.953075 0.365774 -19.009212 0.000000e+00
11 b_time_s_class1 3.277191 0.368806 8.885960 0.000000e+00
12 b_cost_class1 -4.747135 0.252359 -18.811054 0.000000e+00
13 asc_sm_s_class1 1.707598 0.332173 5.140696 2.737226e-07
14 asc_car_class1 1.013010 0.216171 4.686158 2.783819e-06
15 asc_car_s_class1 2.874361 0.261181 11.005234 0.000000e+00
Other parameters
Name Value BHHH std err. BHHH t-stat. BHHH p-value
0 class_cte -1.276545 0.483859 -2.638259 0.008333
1 class_inc -0.201005 0.179155 -1.121966 0.261877
Total running time of the script: (0 minutes 0.149 seconds)