13: Destination in CBD

The models considered to this point include variables that describe the attributes of alternatives, modes, and the characteristics of decision-makers (the work commuters). The mode choice decision also is influenced by variables that describe the context in which the trip is made.

For example, a work trip to the regional central business district (CBD) is more likely to be made by transit than an otherwise similar trip to a suburban work place because the CBD is generally well-served by transit, has more opportunities to make additional stops by walking and is less auto friendly due to congestion and limited and expensive parking. This suggests that the model specification can be enhanced by including variables related to the context of the trip, such as destination zone location.

Model 13 adds the alternative specific CBD dummy variable which indicates whether the destination zone (workplace) is located in the CBD. (pp. 122)

import larch

larch.__version__

'6.0.37.dev32+gdab7641'

This example is a mode choice model built using the MTC example dataset. First we create the DB and Model objects:

d = larch.examples.MTC(format="dataset")
d

<xarray.Dataset> Size: 2MB
Dimensions:    (caseid: 5029, altid: 6)
Coordinates:
  * caseid     (caseid) int64 40kB 1 2 3 4 5 6 ... 5024 5025 5026 5027 5028 5029
  * altid      (altid) int64 48B 1 2 3 4 5 6
    alt_names  (altid) <U7 168B 'DA' 'SR2' 'SR3+' 'Transit' 'Bike' 'Walk'
Data variables: (12/38)
    chose      (caseid, altid) float32 121kB 1.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0
    ivtt       (caseid, altid) float64 241kB 13.38 18.38 20.38 ... 1.59 6.55 0.0
    ovtt       (caseid, altid) float64 241kB 2.0 2.0 2.0 15.2 ... 16.0 4.5 0.0
    tottime    (caseid, altid) float64 241kB 15.38 20.38 22.38 ... 11.05 19.1
    totcost    (caseid, altid) float64 241kB 70.63 35.32 20.18 ... 75.0 0.0 0.0
    hhid       (caseid) int64 40kB 2 3 5 6 8 8 ... 9429 9430 9433 9434 9436 9438
    ...         ...
    corredis   (caseid) int64 40kB 0 1 0 0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 0 0
    vehbywrk   (caseid) float64 40kB 4.0 1.0 0.33 1.0 0.0 ... 2.0 2.0 3.0 3.0
    vocc       (caseid) int64 40kB 1 0 1 0 2 0 1 1 1 1 0 ... 1 2 1 1 0 1 2 1 1 1
    wgt        (caseid) int64 40kB 1 1 1 1 1 1 1 1 1 1 1 ... 1 1 1 1 1 1 1 1 1 1
    _avail_    (caseid, altid) int8 30kB 1 1 1 1 1 0 1 1 1 ... 1 0 1 1 1 1 1 1 1
    avail      (caseid, altid) int8 30kB 1 1 1 1 1 0 1 1 1 ... 1 0 1 1 1 1 1 1 1
Attributes:
    _caseid_:  caseid
    _altid_:   altid

xarray.Dataset

• Dimensions:
  • caseid: 5029
  • altid: 6
• Coordinates: (3)
  • caseid
    (caseid)
    int64
    1 2 3 4 5 ... 5026 5027 5028 5029

    array([   1,    2,    3, ..., 5027, 5028, 5029], shape=(5029,))

  • altid
    (altid)
    int64
    1 2 3 4 5 6

    array([1, 2, 3, 4, 5, 6])

  • alt_names
    (altid)
    <U7
    'DA' 'SR2' 'SR3+' ... 'Bike' 'Walk'

    array(['DA', 'SR2', 'SR3+', 'Transit', 'Bike', 'Walk'], dtype='<U7')

• Data variables: (38)
  • chose
    (caseid, altid)
    float32
    1.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0

    array([[1., 0., 0., 0., 0., 0.],
           [0., 0., 0., 1., 0., 0.],
           [1., 0., 0., 0., 0., 0.],
           ...,
           [1., 0., 0., 0., 0., 0.],
           [1., 0., 0., 0., 0., 0.],
           [1., 0., 0., 0., 0., 0.]], shape=(5029, 6), dtype=float32)

  • ivtt
    (caseid, altid)
    float64
    13.38 18.38 20.38 ... 1.59 6.55 0.0

    array([[13.38, 18.38, 20.38, 25.9 , 40.5 ,  0.  ],
           [29.92, 34.92, 21.92, 22.96, 58.95,  0.  ],
           [ 8.6 , 13.6 , 15.6 , 16.87,  0.  ,  0.  ],
           ...,
           [17.26, 22.26, 24.26,  0.  ,  0.  ,  0.  ],
           [17.2 , 22.2 , 24.2 ,  8.87,  0.  ,  0.  ],
           [ 6.54, 11.54, 13.54,  1.59,  6.55,  0.  ]], shape=(5029, 6))

  • ovtt
    (caseid, altid)
    float64
    2.0 2.0 2.0 15.2 ... 16.0 4.5 0.0

    array([[ 2.  ,  2.  ,  2.  , 15.2 ,  2.  ,  0.  ],
           [10.  , 10.  , 10.  , 14.2 , 10.  ,  0.  ],
           [ 6.  ,  6.  ,  6.  , 21.4 ,  0.  ,  0.  ],
           ...,
           [ 3.  ,  3.  ,  3.  ,  0.  ,  0.  ,  0.  ],
           [ 2.  ,  2.  ,  2.  , 19.02,  0.  ,  0.  ],
           [ 4.5 ,  4.5 ,  4.5 , 16.  ,  4.5 ,  0.  ]], shape=(5029, 6))

  • tottime
    (caseid, altid)
    float64
    15.38 20.38 22.38 ... 11.05 19.1

    array([[15.38, 20.38, 22.38, 41.1 , 42.5 ,  0.  ],
           [39.92, 44.92, 31.92, 37.16, 68.95,  0.  ],
           [14.6 , 19.6 , 21.6 , 38.27,  0.  ,  0.  ],
           ...,
           [20.26, 25.26, 27.26,  0.  ,  0.  ,  0.  ],
           [19.2 , 24.2 , 26.2 , 27.89,  0.  , 61.  ],
           [11.04, 16.04, 18.04, 17.59, 11.05, 19.1 ]], shape=(5029, 6))

  • totcost
    (caseid, altid)
    float64
    70.63 35.32 20.18 ... 75.0 0.0 0.0

    array([[ 70.63,  35.32,  20.18, 115.64,   0.  ,   0.  ],
           [390.81, 195.4 ,  97.97, 185.  ,   0.  ,   0.  ],
           [ 37.76,  18.88,  10.79, 105.  ,   0.  ,   0.  ],
           ...,
           [112.05,  56.03,  32.01,   0.  ,   0.  ,   0.  ],
           [ 41.42,  20.71,  11.83,  75.  ,   0.  ,   0.  ],
           [ 11.42,   5.71,   3.26,  75.  ,   0.  ,   0.  ]], shape=(5029, 6))

  • hhid
    (caseid)
    int64
    2 3 5 6 8 ... 9433 9434 9436 9438

    array([   2,    3,    5, ..., 9434, 9436, 9438], shape=(5029,))

  • perid
    (caseid)
    int64
    1 1 1 1 2 3 2 2 ... 2 1 1 1 1 1 1 1

    array([1, 1, 1, ..., 1, 1, 1], shape=(5029,))

  • numalts
    (caseid)
    int64
    2 2 2 2 2 2 2 2 ... 2 2 2 2 2 2 2 2

    array([2, 2, 2, ..., 2, 2, 2], shape=(5029,))

  • dist
    (caseid)
    float64
    7.69 11.62 4.1 ... 12.06 2.95 0.73

    array([ 7.69, 11.62,  4.1 , ..., 12.06,  2.95,  0.73], shape=(5029,))

  • wkzone
    (caseid)
    int64
    664 738 696 665 ... 1009 1030 1021

    array([ 664,  738,  696, ..., 1009, 1030, 1021], shape=(5029,))

  • hmzone
    (caseid)
    int64
    726 9 667 8 ... 1023 1019 1030 1021

    array([ 726,    9,  667, ..., 1019, 1030, 1021], shape=(5029,))

  • rspopden
    (caseid)
    float64
    15.52 35.81 85.35 ... 2.98 14.77

    array([15.52, 35.81, 85.35, ..., 18.33,  2.98, 14.77], shape=(5029,))

  • rsempden
    (caseid)
    float64
    9.96 53.33 105.9 ... 19.06 40.51

    array([  9.96,  53.33, 105.93, ...,  39.67,  19.06,  40.51], shape=(5029,))

  • wkpopden
    (caseid)
    float64
    37.26 32.91 29.45 ... 2.98 14.77

    array([37.26, 32.91, 29.45, ...,  3.82,  2.98, 14.77], shape=(5029,))

  • wkempden
    (caseid)
    float64
    3.48 764.2 22.42 ... 19.06 40.51

    array([  3.48, 764.19,  22.42, ...,  19.55,  19.06,  40.51], shape=(5029,))

  • vehavdum
    (caseid)
    int64
    1 1 1 1 0 0 1 1 ... 1 1 1 1 1 1 1 1

    array([1, 1, 1, ..., 1, 1, 1], shape=(5029,))

  • femdum
    (caseid)
    int64
    0 0 1 1 1 1 0 0 ... 0 0 1 0 1 0 0 1

    array([0, 0, 1, ..., 0, 0, 1], shape=(5029,))

  • age
    (caseid)
    int64
    35 40 28 34 43 ... 58 33 34 35 37

    array([35, 40, 28, ..., 34, 35, 37], shape=(5029,))

  • drlicdum
    (caseid)
    int64
    1 1 1 1 0 0 1 1 ... 1 1 1 1 1 1 1 1

    array([1, 1, 1, ..., 1, 1, 1], shape=(5029,))

  • noncadum
    (caseid)
    int64
    0 0 1 1 1 1 1 1 ... 1 0 1 0 0 1 0 0

    array([0, 0, 1, ..., 1, 0, 0], shape=(5029,))

  • numveh
    (caseid)
    int64
    4 1 1 2 0 0 2 2 ... 2 2 2 2 2 2 3 3

    array([4, 1, 1, ..., 2, 3, 3], shape=(5029,))

  • hhsize
    (caseid)
    int64
    1 1 5 3 3 3 2 4 ... 4 3 2 4 3 4 3 5

    array([1, 1, 5, ..., 4, 3, 5], shape=(5029,))

  • hhinc
    (caseid)
    float64
    42.5 17.5 12.5 ... 47.5 42.5 47.5

    array([42.5, 17.5, 12.5, ..., 47.5, 42.5, 47.5], shape=(5029,))

  • famtype
    (caseid)
    int64
    7 7 12 3 6 6 8 ... 5 5 99 99 99 99

    array([ 7,  7, 12, ..., 99, 99, 99], shape=(5029,))

  • hhowndum
    (caseid)
    int64
    0 0 0 1 1 1 0 1 ... 1 1 1 0 1 1 1 0

    array([0, 0, 0, ..., 1, 1, 0], shape=(5029,))

  • numemphh
    (caseid)
    int64
    1 1 3 2 3 3 2 2 ... 2 1 1 1 1 1 1 1

    array([1, 1, 3, ..., 1, 1, 1], shape=(5029,))

  • numadlt
    (caseid)
    int64
    1 1 3 2 3 3 2 2 ... 2 3 1 1 1 1 1 1

    array([1, 1, 3, ..., 1, 1, 1], shape=(5029,))

  • nmlt5
    (caseid)
    int64
    0 0 2 1 0 0 0 0 ... 2 0 0 1 1 0 0 2

    array([0, 0, 2, ..., 0, 0, 2], shape=(5029,))

  • nm5to11
    (caseid)
    int64
    0 0 0 0 0 0 0 1 ... 0 0 1 0 1 3 2 2

    array([0, 0, 0, ..., 3, 2, 2], shape=(5029,))

  • nm12to16
    (caseid)
    int64
    0 0 0 0 0 0 0 1 ... 0 0 0 2 0 0 0 0

    array([0, 0, 0, ..., 0, 0, 0], shape=(5029,))

  • wkccbd
    (caseid)
    int64
    0 1 0 1 0 0 0 0 ... 0 0 0 0 0 0 0 0

    array([0, 1, 0, ..., 0, 0, 0], shape=(5029,))

  • wknccbd
    (caseid)
    int64
    0 0 1 0 1 1 1 0 ... 0 0 0 0 0 0 0 0

    array([0, 0, 1, ..., 0, 0, 0], shape=(5029,))

  • corredis
    (caseid)
    int64
    0 1 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0

    array([0, 1, 0, ..., 0, 0, 0], shape=(5029,))

  • vehbywrk
    (caseid)
    float64
    4.0 1.0 0.33 1.0 ... 2.0 3.0 3.0

    array([4.  , 1.  , 0.33, ..., 2.  , 3.  , 3.  ], shape=(5029,))

  • vocc
    (caseid)
    int64
    1 0 1 0 2 0 1 1 ... 1 1 0 1 2 1 1 1

    array([1, 0, 1, ..., 1, 1, 1], shape=(5029,))

  • wgt
    (caseid)
    int64
    1 1 1 1 1 1 1 1 ... 1 1 1 1 1 1 1 1

    array([1, 1, 1, ..., 1, 1, 1], shape=(5029,))

  • _avail_
    (caseid, altid)
    int8
    1 1 1 1 1 0 1 1 ... 0 1 1 1 1 1 1 1

    array([[1, 1, 1, 1, 1, 0],
           [1, 1, 1, 1, 1, 0],
           [1, 1, 1, 1, 0, 0],
           ...,
           [1, 1, 1, 0, 0, 0],
           [1, 1, 1, 1, 0, 1],
           [1, 1, 1, 1, 1, 1]], shape=(5029, 6), dtype=int8)

  • avail
    (caseid, altid)
    int8
    1 1 1 1 1 0 1 1 ... 0 1 1 1 1 1 1 1

    array([[1, 1, 1, 1, 1, 0],
           [1, 1, 1, 1, 1, 0],
           [1, 1, 1, 1, 0, 0],
           ...,
           [1, 1, 1, 0, 0, 0],
           [1, 1, 1, 1, 0, 1],
           [1, 1, 1, 1, 1, 1]], shape=(5029, 6), dtype=int8)

• Indexes: (2)
  • caseid
    PandasIndex

    PandasIndex(Index([   1,    2,    3,    4,    5,    6,    7,    8,    9,   10,
           ...
           5020, 5021, 5022, 5023, 5024, 5025, 5026, 5027, 5028, 5029],
          dtype='int64', name='caseid', length=5029))

  • altid
    PandasIndex

    PandasIndex(Index([1, 2, 3, 4, 5, 6], dtype='int64', name='altid'))

• Attributes: (2)

  _caseid_ :

  caseid

  _altid_ :

  altid

m = larch.Model(d, compute_engine="numba")

Then we can build up the utility function. We’ll use some :ref:idco data first, using the Model.utility.co attribute. This attribute is a dict-like object, to which we can assign :class:LinearFunction objects for each alternative code.

from larch import PX, P, X

for a in [2, 3]:
    m.utility_co[a] = +P("hhinc#2,3") * X("hhinc")
for a in [4, 5, 6]:
    m.utility_co[a] = +P(f"hhinc#{a}") * X("hhinc")

Sometimes we may want to define a part of the utility function that is common across all (or almost all) of the alternatives. We can access a dictionary (more generically called a “mapping”) of alternative codes to alternative names, which can be found via the Dataset.dc.alts_mapping attribute:

d.dc.alts_mapping

{np.int64(1): np.str_('DA'),
 np.int64(2): np.str_('SR2'),
 np.int64(3): np.str_('SR3+'),
 np.int64(4): np.str_('Transit'),
 np.int64(5): np.str_('Bike'),
 np.int64(6): np.str_('Walk')}

Using this like a standard Python dictionary, we can iterate over all the alternatives, skipping 1, and setting alternative specific constants (ASC’s) for the rest.

for a, name in d.dc.alts_mapping.items():
    if a == 1:
        continue
    m.utility_co[a] += (
        +P("ASC_" + name)
        + P("vehbywrk_" + name) * X("vehbywrk")
        + P("wkcbd_" + name) * X("wkccbd + wknccbd")
    )

Next we’ll use some idca data, with the utility_ca attribute. This attribute is only a single :class:LinearFunction that is applied across all alternatives using :ref:idca data. Because the data is structured to vary across alternatives, the parameters (and thus the structure of the :class:LinearFunction) does not need to vary across alternatives.

m.utility_ca = (
    +PX("totcost")
    + P("motorized_time") * X("(altid <= 4) * tottime")
    + P("nonmotorized_time") * X("(altid > 4) * tottime")
    + P("motorized_ovtbydist") * X("(altid <= 4) * ovtt/dist")
)

Lastly, we need to identify idca Format data that gives the availability for each alternative, as well as the number of times each alternative is chosen. (In traditional discrete choice analysis, this is often 0 or 1, but it need not be binary, or even integral.)

m.availability_ca_var = "avail"
m.choice_ca_var = "chose"

And let’s give our model a descriptive title.

m.title = "MTC Example 13, CBD"

We can view a summary of the choices and alternative availabilities to make sure the model is set up correctly.

m.choice_avail_summary()

	name	chosen	available
1	DA	3637	4755
2	SR2	517	5029
3	SR3+	161	5029
4	Transit	498	4003
5	Bike	50	1738
6	Walk	166	1479
< Total All Alternatives >		5029	<NA>

We’ll set a parameter cap (bound) at +/- 20, which helps improve the numerical stability of the optimization algorithm used in estimation.

m.set_cap(20)

Having created this model, we can then estimate it:

assert m.compute_engine == "numba"

result = m.maximize_loglike(stderr=True, method="bhhh")
m.calculate_parameter_covariance()
m.loglike()

Iteration 009 [Optimization terminated successfully]

Best LL = -3440.643668408318

	value	best	initvalue	minimum	maximum	nullvalue	holdfast
param_name
ASC_Bike	-1.650870	-1.650870	0.0	-20.0	20.0	0.0	0
ASC_SR2	-1.634346	-1.634346	0.0	-20.0	20.0	0.0	0
ASC_SR3+	-3.536895	-3.536895	0.0	-20.0	20.0	0.0	0
ASC_Transit	-0.201823	-0.201823	0.0	-20.0	20.0	0.0	0
ASC_Walk	0.083785	0.083785	0.0	-20.0	20.0	0.0	0
hhinc#2,3	-0.001735	-0.001735	0.0	-20.0	20.0	0.0	0
hhinc#4	-0.006149	-0.006149	0.0	-20.0	20.0	0.0	0
hhinc#5	-0.011116	-0.011116	0.0	-20.0	20.0	0.0	0
hhinc#6	-0.007836	-0.007836	0.0	-20.0	20.0	0.0	0
motorized_ovtbydist	-0.150112	-0.150112	0.0	-20.0	20.0	0.0	0
motorized_time	-0.028598	-0.028598	0.0	-20.0	20.0	0.0	0
nonmotorized_time	-0.046414	-0.046414	0.0	-20.0	20.0	0.0	0
totcost	-0.003286	-0.003286	0.0	-20.0	20.0	0.0	0
vehbywrk_Bike	-0.697974	-0.697974	0.0	-20.0	20.0	0.0	0
vehbywrk_SR2	-0.415452	-0.415452	0.0	-20.0	20.0	0.0	0
vehbywrk_SR3+	-0.212059	-0.212059	0.0	-20.0	20.0	0.0	0
vehbywrk_Transit	-0.910910	-0.910910	0.0	-20.0	20.0	0.0	0
vehbywrk_Walk	-0.719448	-0.719448	0.0	-20.0	20.0	0.0	0
wkcbd_Bike	0.375960	0.375960	0.0	-20.0	20.0	0.0	0
wkcbd_SR2	0.255890	0.255890	0.0	-20.0	20.0	0.0	0
wkcbd_SR3+	1.057301	1.057301	0.0	-20.0	20.0	0.0	0
wkcbd_Transit	1.356216	1.356216	0.0	-20.0	20.0	0.0	0
wkcbd_Walk	0.174672	0.174672	0.0	-20.0	20.0	0.0	0

np.float64(-3440.643668408318)

m.parameter_summary()

	Value	Std Err	t Stat	Signif	Null Value
Parameter
ASC_Bike	-1.65	0.428	-3.85	***	0.00
ASC_SR2	-1.63	0.136	-11.99	***	0.00
ASC_SR3+	-3.54	0.204	-17.35	***	0.00
ASC_Transit	-0.202	0.243	-0.83		0.00
ASC_Walk	0.0838	0.348	0.24		0.00
hhinc#2,3	-0.00174	0.00139	-1.25		0.00
hhinc#4	-0.00615	0.00199	-3.09	**	0.00
hhinc#5	-0.0111	0.00524	-2.12	*	0.00
hhinc#6	-0.00784	0.00319	-2.45	*	0.00
motorized_ovtbydist	-0.150	0.0197	-7.61	***	0.00
motorized_time	-0.0286	0.00377	-7.59	***	0.00
nonmotorized_time	-0.0464	0.00572	-8.11	***	0.00
totcost	-0.00329	0.000252	-13.03	***	0.00
vehbywrk_Bike	-0.698	0.256	-2.73	**	0.00
vehbywrk_SR2	-0.415	0.0770	-5.40	***	0.00
vehbywrk_SR3+	-0.212	0.111	-1.92		0.00
vehbywrk_Transit	-0.911	0.115	-7.93	***	0.00
vehbywrk_Walk	-0.719	0.168	-4.28	***	0.00
wkcbd_Bike	0.376	0.321	1.17		0.00
wkcbd_SR2	0.256	0.110	2.33	*	0.00
wkcbd_SR3+	1.06	0.172	6.14	***	0.00
wkcbd_Transit	1.36	0.161	8.41	***	0.00
wkcbd_Walk	0.175	0.225	0.78		0.00

It is a little tough to read this report because the parameters show up in alphabetical order. We can use the reorder method to fix this and group them systematically:

m.ordering = (
    (
        "LOS",
        ".*cost.*",
        ".*time.*",
        ".*dist.*",
    ),
    (
        "Income",
        "hhinc.*",
    ),
    ("Ownership", "vehbywrk.*"),
    ("Zonal", "wkcbd.*"),
    (
        "ASCs",
        "ASC.*",
    ),
)

m.parameter_summary()

		Value	Std Err	t Stat	Signif	Null Value
Category	Parameter
LOS	totcost	-0.00329	0.000252	-13.03	***	0.00
	motorized_time	-0.0286	0.00377	-7.59	***	0.00
	nonmotorized_time	-0.0464	0.00572	-8.11	***	0.00
	motorized_ovtbydist	-0.150	0.0197	-7.61	***	0.00
Income	hhinc#2,3	-0.00174	0.00139	-1.25		0.00
	hhinc#4	-0.00615	0.00199	-3.09	**	0.00
	hhinc#5	-0.0111	0.00524	-2.12	*	0.00
	hhinc#6	-0.00784	0.00319	-2.45	*	0.00
Ownership	vehbywrk_Bike	-0.698	0.256	-2.73	**	0.00
	vehbywrk_SR2	-0.415	0.0770	-5.40	***	0.00
	vehbywrk_SR3+	-0.212	0.111	-1.92		0.00
	vehbywrk_Transit	-0.911	0.115	-7.93	***	0.00
	vehbywrk_Walk	-0.719	0.168	-4.28	***	0.00
Zonal	wkcbd_Bike	0.376	0.321	1.17		0.00
	wkcbd_SR2	0.256	0.110	2.33	*	0.00
	wkcbd_SR3+	1.06	0.172	6.14	***	0.00
	wkcbd_Transit	1.36	0.161	8.41	***	0.00
	wkcbd_Walk	0.175	0.225	0.78		0.00
ASCs	ASC_Bike	-1.65	0.428	-3.85	***	0.00
	ASC_SR2	-1.63	0.136	-11.99	***	0.00
	ASC_SR3+	-3.54	0.204	-17.35	***	0.00
	ASC_Transit	-0.202	0.243	-0.83		0.00
	ASC_Walk	0.0838	0.348	0.24		0.00

Finally, let’s print model statistics. Note that if you want LL at constants you need to run a separate model.

m.estimation_statistics()

Statistic	Aggregate	Per Case
Number of Cases	5029
Log Likelihood at Convergence	-3440.64	-0.68
Log Likelihood at Null Parameters	-7309.60	-1.45
Rho Squared w.r.t. Null Parameters	0.529