.ipynb .pdf

repository open issue

1: MTC MNL Mode Choice

import pandas as pd
import larch.numba as lx

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

d = lx.examples.MTC(format='dataset')
d

<xarray.Dataset>
Dimensions:    (caseid: 5029, altid: 6)
Coordinates:
  * caseid     (caseid) int64 1 2 3 4 5 6 7 ... 5024 5025 5026 5027 5028 5029
  * altid      (altid) int64 1 2 3 4 5 6
    alt_names  (altid) <U7 'DA' 'SR2' 'SR3+' 'Transit' 'Bike' 'Walk'
Data variables: (12/38)
    hhid       (caseid) int64 2 3 5 6 8 8 12 ... 9429 9430 9433 9434 9436 9438
    perid      (caseid) int64 1 1 1 1 2 3 2 2 1 1 2 2 ... 2 1 1 2 1 1 1 1 1 1 1
    numalts    (caseid) int64 2 2 2 2 2 2 2 2 2 2 2 2 ... 2 2 2 2 2 2 2 2 2 2 2
    dist       (caseid) float64 7.69 11.62 4.1 14.58 ... 4.34 12.06 2.95 0.73
    wkzone     (caseid) int64 664 738 696 665 679 ... 1002 1012 1009 1030 1021
    hmzone     (caseid) int64 726 9 667 8 704 665 ... 978 1023 1019 1030 1021
    ...         ...
    ivtt       (caseid, altid) float64 13.38 18.38 20.38 25.9 ... 1.59 6.55 0.0
    ovtt       (caseid, altid) float64 2.0 2.0 2.0 15.2 2.0 ... 4.5 16.0 4.5 0.0
    tottime    (caseid, altid) float64 15.38 20.38 22.38 ... 17.59 11.05 19.1
    totcost    (caseid, altid) float64 70.63 35.32 20.18 115.6 ... 75.0 0.0 0.0
    altnum     (caseid, altid) int64 1 2 3 4 5 0 1 2 3 4 ... 3 4 0 6 1 2 3 4 5 6
    avail      (caseid, altid) int8 1 1 1 1 1 0 1 1 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])

  • 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)
  • hhid
    (caseid)
    int64
    2 3 5 6 8 ... 9433 9434 9436 9438

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

  • 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])

  • 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])

  • 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])

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

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

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

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

  • 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])

  • 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])

  • 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])

  • 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])

  • 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])

  • 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])

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

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

  • 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])

  • 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])

  • 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])

  • 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])

  • 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])

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

    array([ 7,  7, 12, ..., 99, 99, 99])

  • 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])

  • 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])

  • 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])

  • 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])

  • 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])

  • 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])

  • 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])

  • 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])

  • 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])

  • 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.  ])

  • 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])

  • 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])

  • chose
    (caseid, altid)
    int64
    1 0 0 0 0 0 0 0 ... 0 0 1 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]])

  • 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.  ]])

  • 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.  ]])

  • 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 ]])

  • 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.  ]])

  • altnum
    (caseid, altid)
    int64
    1 2 3 4 5 0 1 2 ... 0 6 1 2 3 4 5 6

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

  • 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]], dtype=int8)

• Attributes: (2)

  _caseid_ :

  caseid

  _altid_ :

  altid

m = lx.Model(d)

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 P, X, PX
m.utility_co[2] = P("ASC_SR2")  + P("hhinc#2") * X("hhinc")
m.utility_co[3] = P("ASC_SR3P") + P("hhinc#3") * X("hhinc")
m.utility_co[4] = P("ASC_TRAN") + P("hhinc#4") * X("hhinc")
m.utility_co[5] = P("ASC_BIKE") + P("hhinc#5") * X("hhinc")
m.utility_co[6] = P("ASC_WALK") + P("hhinc#6") * X("hhinc")

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("tottime") + PX("totcost")

Lastly, we need to identify :ref:idca 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_var = 'avail'
m.choice_ca_var = 'chose'

And let’s give our model a descriptive title.

m.title = "MTC Example 1 (Simple MNL)"

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
altid
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

Having created this model, we can then estimate it:

m.maximize_loglike()

key	value
loglike	-3626.18625551293
x	0 ASC_BIKE -2.376e+00 ASC_SR2 -2.178e+00 ASC_SR3P -3.725e+00 ASC_TRAN -6.709e-01 ASC_WALK -2.068e-01 hhinc#2 -2.170e-03 hhinc#3 3.577e-04 hhinc#4 -5.286e-03 hhinc#5 -1.281e-02 hhinc#6 -9.686e-03 totcost -4.920e-03 tottime -5.134e-02
tolerance	2.0071751679019565e-06
steps	array([1., 1., 1., 1., 1., 1., 1., 1., 1., 1.])
message	'Optimization terminated successfully.'
elapsed_time	0:00:00.230719
method	'bhhh'
n_cases	5029
iteration_number	0
logloss	0.7210551313408093
__verbose_repr__	True

m.calculate_parameter_covariance()

m.parameter_summary()

	Value	Std Err	t Stat	Signif	Null Value
ASC_BIKE	-2.38	0.305	-7.80	***	0.00
ASC_SR2	-2.18	0.105	-20.81	***	0.00
ASC_SR3P	-3.73	0.178	-20.96	***	0.00
ASC_TRAN	-0.671	0.133	-5.06	***	0.00
ASC_WALK	-0.207	0.194	-1.07		0.00
hhinc#2	-0.00217	0.00155	-1.40		0.00
hhinc#3	0.000358	0.00254	0.14		0.00
hhinc#4	-0.00529	0.00183	-2.89	**	0.00
hhinc#5	-0.0128	0.00532	-2.41	*	0.00
hhinc#6	-0.00969	0.00303	-3.19	**	0.00
totcost	-0.00492	0.000239	-20.60	***	0.00
tottime	-0.0513	0.00310	-16.57	***	0.00

It is a little tough to read this report because the parameters can show up in pretty much any order, as they are not sorted when they are automatically discovered by Larch. We can use the reorder method to fix this:

m.ordering = (
    ("LOS", "totcost", "tottime", ),
    ("ASCs", "ASC.*", ),
    ("Income", "hhinc.*", ),
)

m.parameter_summary()

		Value	Std Err	t Stat	Signif	Null Value
Category	Parameter
LOS	totcost	-0.00492	0.000239	-20.60	***	0.00
	tottime	-0.0513	0.00310	-16.57	***	0.00
ASCs	ASC_BIKE	-2.38	0.305	-7.80	***	0.00
	ASC_SR2	-2.18	0.105	-20.81	***	0.00
	ASC_SR3P	-3.73	0.178	-20.96	***	0.00
	ASC_TRAN	-0.671	0.133	-5.06	***	0.00
	ASC_WALK	-0.207	0.194	-1.07		0.00
Income	hhinc#2	-0.00217	0.00155	-1.40		0.00
	hhinc#3	0.000358	0.00254	0.14		0.00
	hhinc#4	-0.00529	0.00183	-2.89	**	0.00
	hhinc#5	-0.0128	0.00532	-2.41	*	0.00
	hhinc#6	-0.00969	0.00303	-3.19	**	0.00

m.estimation_statistics()

Statistic	Aggregate	Per Case
Number of Cases	5029
Log Likelihood at Convergence	-3626.19	-0.72
Log Likelihood at Null Parameters	-7309.60	-1.45
Rho Squared w.r.t. Null Parameters	0.504

previous

MTC Work Mode Choice Data

next

2: MTC MNL Mode Choice, Motorized