Modeling Setup

In this vignette we present the modeling workflow for three examples, two being real world data sets and the third using a simulated data set. The core function described is the setup function.

The setup() Function

For the estimation of a multinomial probit model (MNP model) we need to

The setup function combines the first two of these steps to provide an object of the Rprobit class. We use an R6 class for these objects. For details on the perils of R6 classes see this website.

The input arguments for the setup() function depend on the use case. This will be described for three use cases in the following.

Setup Example for Empirical Choice Data in Wide Format

In the choice data vignette Data and Simulation we discuss setting up the data structure for the Train data set in mlogit (Croissant 2020). For setting up the model including the data set we have to do little more than to set up the choice data structure:

data("Train", package = "mlogit")

train_model <- setup(
  formula = choice ~ price + time + change | 0 | comfort,
  df = Train,
  column_decider_id = "id",
  column_occasion_id = "choiceid"
)
train_model$mod$H_b <- train_model$mod$H_b[, -1]
train_model$mod$f_b[1] <- -1
train_model$mod$H_L <- matrix(c(0, 0, 1), ncol = 1)
train_model$mod$f_L <- c(0, 0, 0)

print(train_model)
#> <Rprobit object>
#> 
#> Deciders: 235
#> Choice occasions: from 5 to 19
#> Formula: choice ~ price + time + change | 0 | comfort
#> Random effects: none
#> 
#> Utility equations:
#> 
#>   U_A = b[1] * price_A + b[2] * time_A + b[3] * change_A + b[4] * comfort_A + eps[1]
#>   U_B = b[1] * price_B + b[2] * time_B + b[3] * change_B + b[5] * comfort_B + eps[2]

This model assumes that the price coefficient is fixed at b[1]=−1b[1] = -1 and includes a free Cholesky diagonal parameter for the error scale of the second alternative. The object train_model now contains the choice data set and the model.

train_model$choice_data
#> <choice_data object>
#> 
#> Data format: wide
#> Deciders: 235
#> Choice occasions: from 5 to 19
#> Alternatives: A, B
#> Alternative-constant covariates: none
#> Alternative-varying covariates: price, time, change, comfort
train_model$mod
#> <mod object>
#> 
#> Formula: choice ~ price + time + change | 0 | comfort
#> Random effects: none
#> Alternative names: A, B
#> Coefficients: 5 fixed, 0 random
#> Free parameters total: 5
#> 
#> Utility equations:
#> 
#>   U_A = b[1] * price_A + b[2] * time_A + b[3] * change_A + b[4] * comfort_A + eps[1]
#>   U_B = b[1] * price_B + b[2] * time_B + b[3] * change_B + b[5] * comfort_B + eps[2]
#> 
#> Parameter restrictions:
#> 
#>   b = H_b theta_b + f_b
#>   length(theta_b): 4
#>   H_b (5 x 4) =
#>                theta_1 theta_2 theta_3 theta_4
#> b[1]:price           0       0       0       0
#> b[2]:time            1       0       0       0
#> b[3]:change          0       1       0       0
#> b[4]:comfort_A       0       0       1       0
#> b[5]:comfort_B       0       0       0       1
#>   f_b (5 x 1) =
#>             
#> price     -1
#> time       0
#> change     0
#> comfort_A  0
#> comfort_B  0
#> 
#> 
#>   No parameter O
#> 
#>   L = vech(chol(Sigma)) = H_L theta_L
#>   length(theta_L): 1
#>   H_L (3 x 1) =
#>       theta_5
#> A_:_A       0
#> B_:_A       0
#> B_:_B       1

Additionally the object contains information initialized to NULL that will be filled during the estimation:

This will be filled after estimation is done. For estimation we need to define an initial estimate and subsequently run the numerical search. Details on this can be found in the vignette Model Fitting. Here we only present the workflow for the case study:

train_model$fit(print.level = 0, cml_type = "probit")
#> Looking for initial values...
#> Starting numerical optimization...
#> Calculating approximated Hessian at estimate...
#> Done.

Now the fields in train_model relating to the estimated model have been filled. Thanks to the R6 classes, there is no need to assign the output of the fitting to a new object. See the vignette Model Fitting for details.

train_model$theta_fit
#> [1]  -19.58247 -222.92894 -647.03749 -664.30029 1155.00549
train_model$ll_fit
#> [1] -1727.507
train_model$theta_vcov
#>           theta_1     theta_2   theta_3    theta_4     theta_5
#> theta_1  2.528401    9.387623   13.4429   14.91971   -1.252003
#> theta_2  9.387623 1546.301778  228.1602  210.88997   98.160764
#> theta_3 13.442896  228.160209 1748.2758 1370.67093 -310.918310
#> theta_4 14.919709  210.889965 1370.6709 1759.67200 -349.156130
#> theta_5 -1.252003   98.160764 -310.9183 -349.15613 2313.558397
train_model$hessian_fit
#>               theta_1       theta_2       theta_3       theta_4       theta_5
#> theta_1  0.4237194903 -2.057622e-03 -9.528337e-04 -2.645602e-03 -2.107169e-04
#> theta_2 -0.0020576223  6.740800e-04 -6.014744e-05 -2.472923e-05 -4.152895e-05
#> theta_3 -0.0009528337 -6.014744e-05  1.479188e-03 -1.130923e-03  3.014801e-05
#> theta_4 -0.0026456017 -2.472923e-05 -1.130923e-03  1.488953e-03  7.234203e-05
#> theta_5 -0.0002107169 -4.152895e-05  3.014801e-05  7.234203e-05  4.488518e-04
train_model$Vhat_fit
#>               theta_1       theta_2       theta_3       theta_4       theta_5
#> theta_1  0.4237194903 -2.057622e-03 -9.528337e-04 -2.645602e-03 -2.107169e-04
#> theta_2 -0.0020576223  6.740800e-04 -6.014744e-05 -2.472923e-05 -4.152895e-05
#> theta_3 -0.0009528337 -6.014744e-05  1.479188e-03 -1.130923e-03  3.014801e-05
#> theta_4 -0.0026456017 -2.472923e-05 -1.130923e-03  1.488953e-03  7.234203e-05
#> theta_5 -0.0002107169 -4.152895e-05  3.014801e-05  7.234203e-05  4.488518e-04

A summary provides more information on the fit:

summary(train_model)
#> 
#> Rprobit model summary
#> 
#> Formula:
#> choice ~ price + time + change | 0 | comfort 
#> 
#> Data:
#>   Deciders: 235 
#>   Choice occasions: 2929 
#>   Choice occasions per decider: 5 to 19 
#>   Alternatives: 2 (A, B)
#> 
#> Response residuals:
#>               A        B
#> Min    -0.98105 -0.97971
#> 1Q     -0.40121 -0.40669
#> Median  0.01398 -0.01398
#> 3Q      0.40669  0.40121
#> Max     0.97971  0.98105
#> 
#> Coefficients:
#>              Estimate Std. Error z value Pr(>|z|)    
#> price           -1.00         NA      NA       NA    
#> time           -19.58       1.59 -12.315  < 2e-16 ***
#> change        -222.93      39.32  -5.669 1.43e-08 ***
#> comfort_A     -647.04      41.81 -15.475  < 2e-16 ***
#> comfort_B     -664.30      41.95 -15.836  < 2e-16 ***
#> Sigma[A,A]       0.00         NA      NA       NA    
#> Sigma[B,A]       0.00         NA      NA       NA    
#> Sigma[B,B] 1334037.67  111110.29  12.006  < 2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Error covariance matrix (Sigma):
#>   A       B
#> A 0       0
#> B 0 1334038
#> 
#> Prediction summary:
#>                                 A      B
#> Observed choice frequency  0.5032 0.4968
#> Mean predicted probability 0.5011 0.4989
#> 
#> Fit statistics:
#>   Log-likelihood: -1728 
#>   AIC:  3465 
#>   BIC:  3495

We can change the normalization and refit:

train_model$mod$H_b <- diag(5)
train_model$mod$f_b[1] <- 0
train_model$mod$H_L <- matrix(0, 3, 0)
train_model$mod$f_L[3] <- 1

train_model$fit(print.level = 0, cml_type = "probit")
#> Looking for initial values...
#> Starting numerical optimization...
#> Calculating approximated Hessian at estimate...
#> Done.
summary(train_model)
#> 
#> Rprobit model summary
#> 
#> Formula:
#> choice ~ price + time + change | 0 | comfort 
#> 
#> Data:
#>   Deciders: 235 
#>   Choice occasions: 2929 
#>   Choice occasions per decider: 5 to 19 
#>   Alternatives: 2 (A, B)
#> 
#> Response residuals:
#>               A        B
#> Min    -0.98110 -0.97748
#> 1Q     -0.40278 -0.40576
#> Median  0.01417 -0.01417
#> 3Q      0.40576  0.40278
#> Max     0.97748  0.98110
#> 
#> Coefficients:
#>              Estimate Std. Error z value Pr(>|z|)    
#> price      -8.576e-04  3.616e-05 -23.715  < 2e-16 ***
#> time       -1.618e-02  1.541e-03 -10.499  < 2e-16 ***
#> change     -1.873e-01  3.540e-02  -5.291 1.22e-07 ***
#> comfort_A  -5.535e-01  3.996e-02 -13.851  < 2e-16 ***
#> comfort_B  -5.747e-01  3.994e-02 -14.389  < 2e-16 ***
#> Sigma[A,A]  0.000e+00         NA      NA       NA    
#> Sigma[B,A]  0.000e+00         NA      NA       NA    
#> Sigma[B,B]  1.000e+00         NA      NA       NA    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Error covariance matrix (Sigma):
#>   A B
#> A 0 0
#> B 0 1
#> 
#> Prediction summary:
#>                                 A      B
#> Observed choice frequency  0.5032 0.4968
#> Mean predicted probability 0.5031 0.4969
#> 
#> Fit statistics:
#>   Log-likelihood: -1728 
#>   AIC:  3465 
#>   BIC:  3495

Setup Example for Empirical Choice Data in Long Format

In the choice data vignette Data and Simulation we discuss setting up the data structure for the ModeCanada data set in mlogit (Croissant 2020). Again setting up the model uses the function setup:

data("ModeCanada", package = "mlogit")

mode_model <- setup(
  formula = choice ~ cost + ivt | income | 0,
  df = ModeCanada,
  format = "long",
  column_decider_id = "case",
  column_alternative = "alt",
  ref_alt = "train"
)
mode_model$fit(cml_type = "FP")
#> Looking for initial values...
#> Starting numerical optimization...
#> Calculating approximated Hessian at estimate...
#> Done.
summary(mode_model)
#> 
#> Rprobit model summary
#> 
#> Formula:
#> choice ~ cost + ivt | income | 0 
#> 
#> Data:
#>   Deciders: 4324 
#>   Choice occasions: 4324 
#>   Choice occasions per decider: 1 
#>   Alternatives: 4 (train, car, bus, air)
#> 
#> Response residuals:
#>           train     car        bus        air
#> Min    -0.46551 -0.8811 -1.464e-01 -9.983e-01
#> 1Q     -0.19965 -0.2449 -4.519e-03 -1.076e-01
#> Median -0.11763  0.1190 -1.703e-04 -3.720e-44
#> 3Q     -0.06957  0.2701 -1.771e-05  1.971e-01
#> Max     0.97157  1.0000  9.997e-01  9.980e-01
#> 
#> Coefficients:
#>                      Estimate Std. Error z value Pr(>|z|)    
#> cost               -0.0001481  0.0003740  -0.396    0.692    
#> ivt                -0.0015422  0.0002333  -6.611 3.81e-11 ***
#> income_car          0.0152307  0.0014536  10.478  < 2e-16 ***
#> income_bus         -0.0310694  0.0664221  -0.468    0.640    
#> income_air          0.0164195  0.0014255  11.519  < 2e-16 ***
#> ASC_car             0.1028614  0.0850548   1.209    0.227    
#> ASC_bus            -0.6537584  1.7105058  -0.382    0.702    
#> ASC_air            -0.2338992  0.1121592  -2.085    0.037 *  
#> Sigma[train,train]  0.0000000         NA      NA       NA    
#> Sigma[car,train]    0.0000000         NA      NA       NA    
#> Sigma[bus,train]    0.0000000         NA      NA       NA    
#> Sigma[air,train]    0.0000000         NA      NA       NA    
#> Sigma[car,car]      1.0000000         NA      NA       NA    
#> Sigma[bus,car]      0.4725130  1.0079150   0.469    0.639    
#> Sigma[air,car]      0.8777184  0.0187349  46.849  < 2e-16 ***
#> Sigma[bus,bus]      1.5359574  4.4947052   0.342    0.733    
#> Sigma[air,bus]      0.5007707  0.9634913   0.520    0.603    
#> Sigma[air,air]      0.7760658  0.0318343  24.378  < 2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Error covariance matrix (Sigma):
#>       train    car    bus    air
#> train     0 0.0000 0.0000 0.0000
#> car       0 1.0000 0.4725 0.8777
#> bus       0 0.4725 1.5360 0.5008
#> air       0 0.8777 0.5008 0.7761
#> 
#> Prediction summary:
#>                             train    car      bus    air
#> Observed choice frequency  0.1441 0.5118 0.003700 0.3404
#> Mean predicted probability 0.1591 0.4987 0.006278 0.3359
#> 
#> Fit statistics:
#>   Log-likelihood: -3046 
#>   AIC:  6118 
#>   BIC:  6201

For estimation the full pairwise composite marginal likelihood is used rather than the probit likelihood, see the vignette Model Fitting for details.

Setup Example for Simulated Choice Data

If no empirical data are supplied, setup() can be used to simulate choice data. In this case, the input argument alt_names defines the number and names of alternatives. The simulation uses centered normal regressor values with a standard deviation set by sd. When no parameter vector is specified it is obtained from independent standard normal draws.

set.seed(1)
sim_model <- setup(
  formula = choice ~ quality + price | income | time,
  re = "quality",
  alt_names = c("bus", "car", "train"),
  ref_alt = "bus",
  corr_re = FALSE,
  error_struc = "diag",
  N = 300,
  Tp = rep(3, 300),
  sd = 2
)

print(sim_model)
#> <Rprobit object>
#> 
#> Deciders: 300
#> Choice occasions: from 3 to 3
#> Formula: choice ~ quality + price | income | time
#> Random effects: quality
#> 
#> Utility equations:
#> 
#>   U_bus = b[1] * price_bus + b[6] * time_bus + b[9] * quality_bus + eps[1]
#>   U_car = b[1] * price_car + b[2] * income + b[4] + b[7] * time_car + b[9] * quality_car + eps[2]
#>   U_train = b[1] * price_train + b[3] * income + b[5] + b[8] * time_train + b[9] * quality_train + eps[3]
head(sim_model$choice_data$df)
#>   id occasion choice quality_bus quality_car quality_train  price_bus
#> 1  1        1    car -1.24248116   -2.937764   -0.17615477  5.3497417
#> 2  1        2    car -4.42939977    0.305373    1.83557786  3.1449302
#> 3  1        3    car  2.24986184    3.547525    0.06287853  1.3568384
#> 4  2        1    car -0.08986722   -1.296142    2.71786116  0.2328252
#> 5  2        2    car -0.03238053   -0.399635    0.22698969 -0.9326263
#> 6  2        3    car  1.88767242    1.378487    0.34861309  3.6924666
#>    price_car price_train      income   time_bus   time_car time_train
#> 1  1.0648922  -0.7934002  2.92707309  0.8736785 -2.4476120   1.405398
#> 2 -2.3520717  -0.2371061 -1.00778735 -0.1318735 -0.4571761  -1.783163
#> 3  4.9673081  -1.9144302  0.05572473  0.6269721  2.5807325  -2.819253
#> 4 -2.6633381   0.8372755 -0.92598589  2.4499124 -3.6550388  -1.907681
#> 5 -2.7215367  -0.2143593 -0.69813561 -2.4666690  1.7468251   1.652745
#> 6 -0.4701849  -0.5978371  2.77431682  0.6281037 -2.7008338   1.621070
sim_model$choice_data$theta_dgp
#>  [1] -0.6264538  0.1836433 -0.8356286  1.5952808  0.3295078 -0.8204684
#>  [7]  0.4874291  0.7383247  0.5757814  0.3053884  1.5117812  0.3898432
sim_model$choice_data$choice_counts
#>   bus   car train 
#>   195   446   259

The simulated data set can be used to estimate a model and compare the results to the true parameter vector:

sim_model$fit(
  theta_init = NULL,
  print.level = 0,
  cml_type = "probit"
)
#> Looking for initial values...
#> Starting numerical optimization...
#> Calculating approximated Hessian at estimate...
#> Done.
summary(sim_model)
#> 
#> Rprobit model summary
#> 
#> Formula:
#> choice ~ quality + price | income | time 
#> 
#> Data:
#>   Deciders: 300 
#>   Choice occasions: 900 
#>   Choice occasions per decider: 3 
#>   Alternatives: 3 (bus, car, train)
#> 
#> Response residuals:
#>               bus        car      train
#> Min    -9.583e-01 -9.436e-01 -0.9392899
#> 1Q     -6.822e-02 -1.419e-01 -0.0479703
#> Median -3.444e-03 -9.248e-06 -0.0002244
#> 3Q     -6.307e-09  1.344e-01  0.0054315
#> Max     9.867e-01  9.859e-01  0.9757626
#> 
#> Coefficients:
#>                          Estimate Std. Error z value Pr(>|z|)    
#> price                  -7.676e-01  1.696e-01  -4.526 6.02e-06 ***
#> income_car              2.286e-01  9.903e-02   2.308   0.0210 *  
#> income_train           -9.104e-01  1.895e-01  -4.805 1.55e-06 ***
#> ASC_car                 1.641e+00  3.462e-01   4.740 2.13e-06 ***
#> ASC_train               1.561e-01  1.575e-01   0.991   0.3217    
#> time_bus               -7.976e-01  1.863e-01  -4.280 1.87e-05 ***
#> time_car                5.455e-01  1.397e-01   3.904 9.47e-05 ***
#> time_train              8.343e-01  1.829e-01   4.562 5.07e-06 ***
#> quality                 5.857e-01  1.329e-01   4.407 1.05e-05 ***
#> Omega[quality,quality]  1.855e-01  1.017e-01   1.824   0.0681 .  
#> Sigma[bus,bus]          1.000e+00         NA      NA       NA    
#> Sigma[car,bus]          0.000e+00         NA      NA       NA    
#> Sigma[train,bus]        0.000e+00         NA      NA       NA    
#> Sigma[car,car]          2.871e+00  1.526e+00   1.882   0.0599 .  
#> Sigma[train,car]        0.000e+00         NA      NA       NA    
#> Sigma[train,train]      3.425e-11  1.170e-02   0.000   1.0000    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Random-effect covariance matrix (Omega):
#>         quality
#> quality  0.1855
#> 
#> Error covariance matrix (Sigma):
#>       bus   car     train
#> bus     1 0.000 0.000e+00
#> car     0 2.871 0.000e+00
#> train   0 0.000 3.425e-11
#> 
#> Prediction summary:
#>                               bus    car  train
#> Observed choice frequency  0.2167 0.4956 0.2878
#> Mean predicted probability 0.2138 0.5002 0.2860
#> 
#> Fit statistics:
#>   Log-likelihood: -399.1 
#>   AIC: 822.2 
#>   BIC: 879.9
rbind(
  "fit" = sim_model$theta_fit,
  "dgp" = sim_model$choice_data$theta_dgp
) |> round(2)
#>      [,1] [,2]  [,3] [,4] [,5]  [,6] [,7] [,8] [,9] [,10] [,11] [,12]
#> fit -0.77 0.23 -0.91 1.64 0.16 -0.80 0.55 0.83 0.59  0.43  1.69  0.00
#> dgp -0.63 0.18 -0.84 1.60 0.33 -0.82 0.49 0.74 0.58  0.31  1.51  0.39

References

Croissant, Y. 2020. “Estimation of Random Utility Models in R: The mlogit Package.” Journal of Statistical Software 95 (11). https://doi.org/10.18637/jss.v095.i11.