111: Swissmetro Cross-Nested Logit Mode Choice

This example is a mode choice model built using the Swissmetro example dataset. First we create the DB and Model objects. When we create the DB object, we will redefine the weight value:

d = larch.DB.Example('SWISSMETRO')
m = larch.Model(d)

We can attach a title to the model. The title does not affect the calculations as all; it is merely used in various output report styles.

m.title = "swissmetro example 11 (cross nested logit)"

The swissmetro dataset, as with all Biogeme data, is only in co format.

from larch.roles import P,X
m.utility[1] = ( P.ASC_TRAIN
               + P.B_TIME * X.TRAIN_TT
               + P.B_COST * X("TRAIN_CO*(GA==0)") )
m.utility[2] = ( P.B_TIME * X.SM_TT
               + P.B_COST * X("SM_CO*(GA==0)") )
m.utility[3] = ( P.ASC_CAR
               + P.B_TIME * X.CAR_TT
               + P.B_COST * X("CAR_CO") )

To create a new nest, we can use the new_nest command, although we’ll need to know what the alternative codes are for the alternatives in our dataset. To find out, we can do:

>>> m.df.alternatives()
[(1, 'Train'), (2, 'SM'), (3, 'Car')]

For this example, we want to nest together the Train and Car modes into a “existing” modes nest, and we want to nest Train and SM together into a “public” modes nest. This creates a structure different from a traditional nested logit model, because the Train mode is “cross-nested”: it appears in more than one nest. The desired nesting structure then looks like this:

digraph Desired_Nesting_Structure {
bgcolor="transparent"
root [label="Root", shape="oval"]
train [label="Train", shape=box, style="rounded", penwidth=2]
sm [label="SM", shape=box, style="rounded", penwidth=2]
car [label="Car", shape=box, style="rounded", penwidth=2]
public [label="Public", shape=oval]
exist [label="Existing", shape=oval]
exist -> train
public -> train
public -> sm
exist -> car
root -> exist
root -> public
}

To do this we can use the new_nest command like this:

existing_id = m.new_nest("existing", parent=m.root_id, children=[1,3])
public_id = m.new_nest("public", parent=m.root_id, children=[1,2])

The new_nest() method returns an identifying code number for the newly created nest. We’ll use that code number below.

For a cross-nested model, we need to assign an allocation level to each graph link for all entering links of any node that has more than one predecessor. In this case, that applies only to the “Train” node.

Larch employs a logit-type function to manage this allocation, instead of specifying the allocation directly as a parameter. So, the allocation on the link Public->Train (PT) is given by

\[\alpha_{PT} = \frac{\exp ( \phi_{PT} Z )}{\exp ( \phi_{PT} Z ) + \exp ( \phi_{ET} Z )},\]

where \(\phi_{PT}\) is a vector of zero or more parameters associated with the link PT, \(\phi_{ET}\) is a similar vector of parameters for the link Public->Existing (ET) and \(Z\) is a vector of idco Format variables.

If we give our model no other commands, the length of these parameter and data vectors will be zero, and the allocation parameters for links PT and ET will default to 0.5 each. But, we can relax this constrained default value by putting some defined parameters into the formula:

m.link[existing_id, 1](param="phi_et", data="1")

Note that, as for alternative specific constants in the utility function, we leave out one of the possibilities (here, we skip ET) to avoid overspecifying the model.

Constructing the allocation parameters like this has a few benefits. It automatically ensures that the total allocation for each node to it’s incoming links is always 1. It also allows the parameter phi_pt to be estimated as unconstrained, as it can take on any value without bound and still produce a properly normalized model. Lastly, it allows the allocation to be a function of the idco Format data, and not merely a fixed value. One drawback is that there is no statistical test to use on the value of the estimated phi parameters to determine if the true value of the alloction should be 1 or 0, because those allocations are associated with infinite values for phi.

Larch will find all the parameters in the model, but we’d like to output them in a particular order, so we want to reorder the parameters. We can use the reorder method to fix this:

m.reorder_parameters("ASC", "B_", "existing", "public", "phi")

We can estimate the models and check the results match up with those given by Biogeme:

>>> result = m.maximize_loglike()
>>> print(result.message)
Optimization terminated successfully...

>>> print(m.report('txt', sigfigs=3))
=========================================================================================...
swissmetro example 11 (cross nested logit)
=========================================================================================...
Model Parameter Estimates
-----------------------------------------------------------------------------------------...
Parameter       InitValue       FinalValue      StdError        t-Stat          NullValue
ASC_TRAIN        0.0             0.0983          0.0563          1.74            0.0
ASC_CAR          0.0            -0.24            0.0384         -6.26            0.0
B_TIME           0.0            -0.00777         0.000558       -13.9            0.0
B_COST           0.0            -0.00819         0.000446       -18.4            0.0
existing         1.0             0.398           0.0276         -21.8            1.0
public           1.0             0.243           0.0336         -22.5            1.0
phi_et           0.0            -0.0197          0.116          -0.17            0.0
=========================================================================================...
Model Estimation Statistics
-----------------------------------------------------------------------------------------...
Log Likelihood at Convergence           -5214.05
Log Likelihood at Null Parameters       -6964.66
-----------------------------------------------------------------------------------------...
Rho Squared w.r.t. Null Parameters      0.251
=========================================================================================...

Tip

If you want access to the model in this example without worrying about assembling all the code blocks together on your own, you can load a read-to-estimate copy like this:

m = larch.Model.Example(111)