R is a programming language for statistical computing and graphics. To install R, required installation files can be found in the following link:

Download RStudio

Apollo Choice Modeling framework has been employed in running the vehicle ownership model component. We prepared the R notebooks for all models included in the guidebook.

Clear memory

rm(list = ls())

Install and load Apollo library

In first step of data processing in R, we need to call required package

library(apollo)
## Apollo 0.2.8
## http://www.ApolloChoiceModelling.com
## See url for a detailed manual, examples and a user forum.
## Sign up to the user forum to receive updates on new releases.
## 
## Your version of Apollo is more than six months old.
## Using the latest version will ensure you have all
##  current functionality and bug fixes.
## You can update to the latest version by typing:
##  install.packages("apollo")

Initialize code and set core controls

apollo_initialise()
apollo_control = list(
  modelName       = "MNL_Vehicle Ownership",
  modelDescr      = "Simple MNL model on vehicle ownership data (NMO model)",
  indivID         = "id",
  outputDirectory = "output",
  panelData = FALSE
)

Load data and check summary

database = read.csv("D:/Projects/20-102/Project Task/Task 4/R code/VO data_NMO.csv",header=TRUE)
summary(database)
##        id            Choice           veh0             veh1       
##  Min.   :    1   Min.   :1.000   Min.   :0.0000   Min.   :0.0000  
##  1st Qu.: 2570   1st Qu.:2.000   1st Qu.:0.0000   1st Qu.:0.0000  
##  Median : 5140   Median :3.000   Median :0.0000   Median :0.0000  
##  Mean   : 5140   Mean   :2.741   Mean   :0.1159   Mean   :0.2577  
##  3rd Qu.: 7709   3rd Qu.:3.000   3rd Qu.:0.0000   3rd Qu.:1.0000  
##  Max.   :10278   Max.   :4.000   Max.   :1.0000   Max.   :1.0000  
##       veh2             veh3           wrktodrv      age615todrv    
##  Min.   :0.0000   Min.   :0.0000   Min.   :0.000   Min.   :0.0000  
##  1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:0.500   1st Qu.:0.0000  
##  Median :0.0000   Median :0.0000   Median :1.000   Median :0.0000  
##  Mean   :0.3958   Mean   :0.2306   Mean   :0.687   Mean   :0.1892  
##  3rd Qu.:1.0000   3rd Qu.:0.0000   3rd Qu.:1.000   3rd Qu.:0.0000  
##  Max.   :1.0000   Max.   :1.0000   Max.   :4.000   Max.   :4.0000  
##    age80todrv          inc_10           inc_1030         inc_3060     
##  Min.   :0.00000   Min.   :0.00000   Min.   :0.0000   Min.   :0.0000  
##  1st Qu.:0.00000   1st Qu.:0.00000   1st Qu.:0.0000   1st Qu.:0.0000  
##  Median :0.00000   Median :0.00000   Median :0.0000   Median :0.0000  
##  Mean   :0.06482   Mean   :0.04077   Mean   :0.1355   Mean   :0.2109  
##  3rd Qu.:0.00000   3rd Qu.:0.00000   3rd Qu.:0.0000   3rd Qu.:0.0000  
##  Max.   :4.00000   Max.   :1.00000   Max.   :1.0000   Max.   :1.0000  
##    inc_60100           nmo            nmo_wrk         nmo_inc60     
##  Min.   :0.0000   Min.   :0.0000   Min.   :0.0000   Min.   :0.0000  
##  1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:0.0000  
##  Median :0.0000   Median :0.0000   Median :0.0000   Median :0.0000  
##  Mean   :0.2544   Mean   :0.1946   Mean   :0.1615   Mean   :0.1191  
##  3rd Qu.:1.0000   3rd Qu.:0.0000   3rd Qu.:0.0000   3rd Qu.:0.0000  
##  Max.   :1.0000   Max.   :1.0000   Max.   :1.0000   Max.   :1.0000  
##       uno   
##  Min.   :1  
##  1st Qu.:1  
##  Median :1  
##  Mean   :1  
##  3rd Qu.:1  
##  Max.   :1
table(database$Choice)
## 
##    1    2    3    4 
## 1191 2649 4068 2370

Here, Choice = Vehicle ownership level, veh0 = zero vehicle HH, veh1 = HH with one vehicle, veh2= HH with two vehicles, veh3 = HH with three or more vehicles, wrktodrv = Ratio of workers to drivers, age615todrv = Ratio of people age 6-15 to drivers, age80todrv = Ratio of people age 80 or older to drivers, inc_10 = Household income <=$10,000, inc_1030 = $10,000 < Household income <= $30,000, inc_3060 = $30,000 < Household income <= $60,000 inc_60100 = $60,000 < Household income <= $100,000, nmo = NMO availability indicator (urban household indicator), nmo_wrk = NMO availability Presence of worker in HH, nmo_inc60 = NMO availabilityHousehold income >= $60,000

Initialize vector of parameters

Include fixed parameters in estimation process. Vector with names (in quotes) of parameters to be kept fixed at their starting value in apollo_beta, use apollo_beta_fixed = c() if none

apollo_beta=c(asc_veh3            = 0,
              asc_veh2            = 0,
              asc_veh1            = 0,
              asc_veh0            = 0,
              b_wrktodrv_veh1     = 0,
              b_wrktodrv_veh0     = 0,
              b_615todrv_veh2     = 0,
              b_80todrv_veh1      = 0,
              b_inc10_veh2        = 0,
              b_inc10_veh1        = 0,
              b_inc10_veh0        = 0,
              b_inc1030_veh2      = 0,
              b_inc1030_veh1      = 0,
              b_inc1030_veh0      = 0,
              b_inc3060_veh2      = 0,
              b_inc3060_veh1      = 0,
              b_inc3060_veh0      = 0,          
              b_inc60100_veh1     = 0,
              b_nmo_veh2          = 0,
              b_nmo_veh1          = 0,
              b_nmo_veh0          = 0,
              b_nmowrk_veh0       = 0,
              b_nmoinc60_veh1     = 0,
              b_nmoinc60_veh0     = 0)



apollo_fixed = c("asc_veh3")


apollo_inputs = apollo_validateInputs()
## All checks on apollo_control completed.
## All checks on database completed.

Define the utilities of MNL model and compute the likelihood function

# Define model and likelihood function

apollo_probabilities=function(apollo_beta, apollo_inputs, functionality="estimate"){
    
  # Attach inputs and detach after function exit
  apollo_attach(apollo_beta, apollo_inputs)
  on.exit(apollo_detach(apollo_beta, apollo_inputs))

  # Create list of probabilities P
  P = list()
  
  # List of utilities: these must use the same names as in mnl_settings, order is irrelevant
  V = list()
  V[["veh3"]]   =   asc_veh3                                                                                                                                                                
  V[["veh2"]]   =   asc_veh2 +  b_615todrv_veh2 *   age615todrv +   b_inc10_veh2    *   inc_10  +   b_inc1030_veh2  *   inc_1030    +   b_inc3060_veh2  *   inc_3060 +  b_nmo_veh2  *   nmo 
  
  V[["veh1"]]   =   asc_veh1    +   b_wrktodrv_veh1 *   wrktodrv +  b_80todrv_veh1  *   age80todrv  +   b_inc10_veh1    *   inc_10  +   b_inc1030_veh1  *   inc_1030    +   b_inc3060_veh1  *   inc_3060    +   b_inc60100_veh1 *   inc_60100   +   b_nmo_veh1  *   nmo +   b_nmoinc60_veh1 *   nmo_inc60

  V[["veh0"]]   =   asc_veh0    +   b_wrktodrv_veh0 *   wrktodrv +  b_inc10_veh0    *   inc_10  +   b_inc1030_veh0  *   inc_1030    +   b_inc3060_veh0  *   inc_3060 +  b_nmo_veh0  *   nmo +   b_nmowrk_veh0   *   nmo_wrk +   b_nmoinc60_veh0 *   nmo_inc60
  
  # Define settings for MNL model component
  mnl_settings = list(
    alternatives  = c(veh0=1, veh1=2, veh2=3, veh3=4), 
    avail         = list(veh0=uno, veh1=uno, veh2=uno, veh3=uno),
    choiceVar     = Choice,
    utilities     = V
  )
  
  # Compute probabilities using MNL model
  P[["model"]] = apollo_mnl(mnl_settings, functionality)
  
  # Prepare and return outputs of function
  P = apollo_prepareProb(P, apollo_inputs, functionality)
  return(P)
  nrow(p)
}

Model estimation

model = apollo_estimate(apollo_beta, apollo_fixed, apollo_probabilities, apollo_inputs)
## Preparing user-defined functions.
## 
## Testing likelihood function...
## 
## Overview of choices for MNL model component :
##                                      veh0     veh1     veh2     veh3
## Times available                  10278.00 10278.00 10278.00 10278.00
## Times chosen                      1191.00  2649.00  4068.00  2370.00
## Percentage chosen overall           11.59    25.77    39.58    23.06
## Percentage chosen when available    11.59    25.77    39.58    23.06
## 
## Pre-processing likelihood function...
## 
## Testing influence of parameters
## Starting main estimation
## Initial function value: -14248.33 
## Initial gradient value:
##        asc_veh2        asc_veh1        asc_veh0 b_wrktodrv_veh1 b_wrktodrv_veh0 
##      1498.50000        79.50000     -1378.50000       -10.51548     -1127.63214 
## b_615todrv_veh2  b_80todrv_veh1    b_inc10_veh2    b_inc10_veh1    b_inc10_veh0 
##       481.76131        67.15417       -47.75000        60.25000        77.25000 
##  b_inc1030_veh2  b_inc1030_veh1  b_inc1030_veh0  b_inc3060_veh2  b_inc3060_veh1 
##       -23.25000       261.75000       -28.25000       212.00000       219.00000 
##  b_inc3060_veh0 b_inc60100_veh1      b_nmo_veh2      b_nmo_veh1      b_nmo_veh0 
##      -287.00000      -105.75000        -6.00000        67.00000       346.00000 
##   b_nmowrk_veh0 b_nmoinc60_veh1 b_nmoinc60_veh0 
##       245.00000        65.00000        94.00000 
## initial  value 14248.333444 
## iter   2 value 13181.974434
## iter   3 value 12332.803321
## iter   4 value 12235.819677
## iter   5 value 11973.627845
## iter   6 value 11844.393466
## iter   7 value 11772.274982
## iter   8 value 11728.143769
## iter   9 value 11693.267412
## iter  10 value 11685.700346
## iter  11 value 11664.531674
## iter  12 value 11642.004568
## iter  13 value 11629.461109
## iter  14 value 11614.332482
## iter  15 value 11594.985655
## iter  16 value 11576.882833
## iter  17 value 11554.239300
## iter  18 value 11525.510069
## iter  19 value 11499.461155
## iter  20 value 11462.121768
## iter  21 value 11366.116651
## iter  22 value 11359.522761
## iter  23 value 11356.742589
## iter  24 value 11349.094256
## iter  25 value 11268.423244
## iter  26 value 11260.973784
## iter  27 value 11240.682031
## iter  28 value 11233.549023
## iter  29 value 11231.058463
## iter  30 value 11230.743665
## iter  31 value 11230.716202
## iter  32 value 11230.711084
## iter  33 value 11230.710229
## iter  33 value 11230.710158
## iter  33 value 11230.710158
## final  value 11230.710158 
## converged
## Additional convergence test using scaled estimation. Parameters will be
##   scaled by their current estimates and additional iterations will be
##   performed.
## initial  value 11230.710158 
## iter   1 value 11230.710157
## final  value 11230.710157 
## converged
## Estimated parameters:
##                    Estimate
## asc_veh3             0.0000
## asc_veh2             0.2925
## asc_veh1            -1.1952
## asc_veh0            -3.3095
## b_wrktodrv_veh1      0.1789
## b_wrktodrv_veh0     -0.8889
## b_615todrv_veh2      0.4124
## b_80todrv_veh1       0.4388
## b_inc10_veh2         0.9517
## b_inc10_veh1         3.4670
## b_inc10_veh0         5.5465
## b_inc1030_veh2       0.4370
## b_inc1030_veh1       2.4909
## b_inc1030_veh0       3.6446
## b_inc3060_veh2       0.1910
## b_inc3060_veh1       1.5377
## b_inc3060_veh0       1.9189
## b_inc60100_veh1      0.3875
## b_nmo_veh2           1.2466
## b_nmo_veh1           1.6371
## b_nmo_veh0           3.8170
## b_nmowrk_veh0        0.9049
## b_nmoinc60_veh1      0.8111
## b_nmoinc60_veh0      1.0194
## 
## Final LL: -11230.7102
## 
## Calculating log-likelihood at equal shares (LL(0)) for applicable
##   models...
## Calculating log-likelihood at observed shares from estimation data
##   (LL(c)) for applicable models...
## Calculating LL of each model component...
## Calculating other model fit measures
## Computing covariance matrix using analytical gradient.
##  0%....25%....50%....75%...100%
## Negative definite Hessian with maximum eigenvalue: -4.040706
## Computing score matrix...

Model outputs

apollo_modelOutput(model)
## Model run by Surface using Apollo 0.2.8 on R 4.1.3 for Windows.
## www.ApolloChoiceModelling.com
## 
## Model name                                  : MNL_Vehicle Ownership
## Model description                           : Simple MNL model on vehicle ownership data (NMO model)
## Model run at                                : 2023-11-20 12:35:38
## Estimation method                           : bfgs
## Model diagnosis                             : successful convergence 
## Number of individuals                       : 10278
## Number of rows in database                  : 10278
## Number of modelled outcomes                 : 10278
## 
## Number of cores used                        :  1 
## Model without mixing
## 
## LL(start)                                   : -14248.33
## LL at equal shares, LL(0)                   : -14248.33
## LL at observed shares, LL(C)                : -13405.94
## LL(final)                                   : -11230.71
## Rho-squared vs equal shares                  :  0.2118 
## Adj.Rho-squared vs equal shares              :  0.2102 
## Rho-squared vs observed shares               :  0.1623 
## Adj.Rho-squared vs observed shares           :  0.1605 
## AIC                                         :  22507.42 
## BIC                                         :  22673.89 
## 
## Estimated parameters                        :  23
## Time taken (hh:mm:ss)                       :  00:00:6.6 
##      pre-estimation                         :  00:00:2.04 
##      estimation                             :  00:00:1.25 
##      post-estimation                        :  00:00:3.3 
## Iterations                                  :  37  
## Min abs eigenvalue of Hessian               :  4.040706 
## 
## Unconstrained optimisation.
## 
## These outputs have had the scaling used in estimation applied to them.
## Estimates:
##                    Estimate        s.e.   t.rat.(0)    Rob.s.e. Rob.t.rat.(0)
## asc_veh3             0.0000          NA          NA          NA            NA
## asc_veh2             0.2925     0.03271       8.943     0.03283         8.910
## asc_veh1            -1.1952     0.07468     -16.003     0.07732       -15.457
## asc_veh0            -3.3095     0.19085     -17.341     0.21682       -15.264
## b_wrktodrv_veh1      0.1789     0.06382       2.803     0.06816         2.624
## b_wrktodrv_veh0     -0.8889     0.12325      -7.212     0.13324        -6.671
## b_615todrv_veh2      0.4124     0.05005       8.240     0.05041         8.181
## b_80todrv_veh1       0.4388     0.08636       5.082     0.08827         4.972
## b_inc10_veh2         0.9517     0.29209       3.258     0.29260         3.253
## b_inc10_veh1         3.4670     0.27700      12.517     0.27943        12.407
## b_inc10_veh0         5.5465     0.32848      16.885     0.33842        16.389
## b_inc1030_veh2       0.4370     0.10643       4.106     0.10707         4.082
## b_inc1030_veh1       2.4909     0.11128      22.384     0.11224        22.193
## b_inc1030_veh0       3.6446     0.21251      17.150     0.22402        16.269
## b_inc3060_veh2       0.1910     0.06957       2.746     0.07053         2.709
## b_inc3060_veh1       1.5377     0.08392      18.323     0.08333        18.454
## b_inc3060_veh0       1.9189     0.21185       9.058     0.21261         9.025
## b_inc60100_veh1      0.3875     0.06751       5.740     0.06763         5.729
## b_nmo_veh2           1.2466     0.11642      10.708     0.11798        10.567
## b_nmo_veh1           1.6371     0.15649      10.461     0.17325         9.449
## b_nmo_veh0           3.8170     0.19581      19.493     0.20765        18.382
## b_nmowrk_veh0        0.9049     0.16886       5.359     0.17638         5.130
## b_nmoinc60_veh1      0.8111     0.14572       5.566     0.14911         5.439
## b_nmoinc60_veh0      1.0194     0.23485       4.341     0.24566         4.150

Here, asc_veh3 = Alternative specific constant for HH with three or more vehicles, asc_veh2 = Alternative specific constant for HH with two vehicles, asc_veh1 = Alternative specific constant for HH with one vehicle, asc_veh0 = Alternative specific constant for HH with zero vehicle, b_wrktodrv_veh1 = Coefficient for ratio of workers to drivers (HH with one vehicle), b_wrktodrv_veh0 = Coefficient for ratio of workers to drivers (zero vehicle HH), b_615todrv_veh2 = Coefficient for ratio of people age 6-15 to drivers (HH with two or more vehicles), b_80todrv_veh1 = Coefficient for ratio of people age 80 or older to drivers to drivers (HH with one vehicle), b_inc10_veh2 = Coefficient for household income <=$10,000 (HH with two vehicles), b_inc10_veh1 = Coefficient for household income <=$10,000 (HH with one vehicle), b_inc10_veh0 = Coefficient for household income <=$10,000 (zero vehicle HH), b_inc1030_veh2 = Coefficient for $10,000 < household income <= $30,000 (HH with two vehicles), b_inc1030_veh1 = Coefficient for $10,000 < household income <= $30,000 (HH with one vehicle), b_inc1030_veh0 = Coefficient for $10,000 < household income <= $30,000 (zero vehicle HH), b_inc3060_veh2 = Coefficient for $30,000 < household income <= $60,000 (HH with two vehicles), b_inc3060_veh1 = Coefficient for $30,000 < household income <= $60,000 (HH with one vehicle), b_inc3060_veh0 = Coefficient for $30,000 < household income <= $60,000 (zero vehicle HH), b_nmo_veh2 = Coefficient for NMO availability indicator (urban household indicator) (HH with two vehicles), b_nmo_veh1 = Coefficient for NMO availability indicator (urban household indicator) (HH with one vehicle), b_nmo_veh0 = Coefficient for NMO availability indicator (urban household indicator) (zero vehicle HH), b_nmowrk_veh0 = Coefficient for NMO availability Presence of worker in HH (zero vehicle HH), b_nmoinc60_veh1 = Coefficient for NMO availability Household income >= $60,000 (HH with one vehicle), b_nmoinc60_veh0 = Coefficient for NMO availability *Household income >= $60,000 (zero vehicle HH)