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 (base 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_base.csv",header=TRUE)
summary(database)
##        id            Choice           veh0              veh1       
##  Min.   :    1   Min.   :1.000   Min.   :0.00000   Min.   :0.0000  
##  1st Qu.: 2570   1st Qu.:2.000   1st Qu.:0.00000   1st Qu.:0.0000  
##  Median : 5140   Median :3.000   Median :0.00000   Median :0.0000  
##  Mean   : 5140   Mean   :2.942   Mean   :0.04145   Mean   :0.2525  
##  3rd Qu.: 7709   3rd Qu.:4.000   3rd Qu.:0.00000   3rd Qu.:1.0000  
##  Max.   :10278   Max.   :4.000   Max.   :1.00000   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.4284   Mean   :0.2777   Mean   :0.687   Mean   :0.1892  
##  3rd Qu.:1.0000   3rd Qu.:1.0000   3rd Qu.:1.000   3rd Qu.:0.0000  
##  Max.   :1.0000   Max.   :1.0000   Max.   :4.000   Max.   :4.0000  
##   age1824todrv       age80todrv          inc_10           inc_1030     
##  Min.   :0.00000   Min.   :0.00000   Min.   :0.00000   Min.   :0.0000  
##  1st Qu.:0.00000   1st Qu.:0.00000   1st Qu.:0.00000   1st Qu.:0.0000  
##  Median :0.00000   Median :0.00000   Median :0.00000   Median :0.0000  
##  Mean   :0.04756   Mean   :0.06482   Mean   :0.04077   Mean   :0.1355  
##  3rd Qu.:0.00000   3rd Qu.:0.00000   3rd Qu.:0.00000   3rd Qu.:0.0000  
##  Max.   :2.00000   Max.   :4.00000   Max.   :1.00000   Max.   :1.0000  
##     inc_3060        inc_60100           uno   
##  Min.   :0.0000   Min.   :0.0000   Min.   :1  
##  1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:1  
##  Median :0.0000   Median :0.0000   Median :1  
##  Mean   :0.2109   Mean   :0.2544   Mean   :1  
##  3rd Qu.:0.0000   3rd Qu.:1.0000   3rd Qu.:1  
##  Max.   :1.0000   Max.   :1.0000   Max.   :1
table(database$Choice)
## 
##    1    2    3    4 
##  426 2595 4403 2854

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, age1824todrv = Ratio of people age 18-24 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_60100 = $60,000 < Household income <= $100,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_615todrv_veh1     = 0,
              b_1824todrv_veh2    = 0,
              b_1824todrv_veh1    = 0,
              b_1824todrv_veh0    = 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_inc60100_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_1824todrv_veh2    *   age1824todrv    +   b_inc10_veh2    *   inc_10 +    b_inc1030_veh2  *   inc_1030 + b_inc3060_veh2   *   inc_3060
  
  V[["veh1"]]   =   asc_veh1 +  b_wrktodrv_veh1 *   wrktodrv + b_615todrv_veh1  *   age615todrv +   b_1824todrv_veh1    *   age1824todrv +  b_80todrv_veh1  *   age80todrv +    b_inc10_veh1    *   inc_10  +   b_inc1030_veh1  *   inc_1030 + b_inc3060_veh1   *   inc_3060 +  b_inc60100_veh1 *   inc_60100
  
  V[["veh0"]]   =   asc_veh0 +  b_wrktodrv_veh0 *   wrktodrv +  b_1824todrv_veh0    *   age1824todrv    +   b_inc10_veh0    *   inc_10  +   b_inc1030_veh0  *   inc_1030    +   b_inc3060_veh0  *   inc_3060    +   b_inc60100_veh0 *   inc_60100
  
  # 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                       426.00  2595.00  4403.00  2854.00
## Percentage chosen overall            4.14    25.25    42.84    27.77
## Percentage chosen when available     4.14    25.25    42.84    27.77
## 
## 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 
##       1833.50000         25.50000      -2143.50000        -31.81548 
##  b_wrktodrv_veh0  b_615todrv_veh2  b_615todrv_veh1 b_1824todrv_veh2 
##      -1700.64881        569.09464       -124.37202         18.55595 
## b_1824todrv_veh1 b_1824todrv_veh0   b_80todrv_veh1     b_inc10_veh2 
##        -29.96071       -101.21071         67.95417        -45.75000 
##     b_inc10_veh1     b_inc10_veh0   b_inc1030_veh2   b_inc1030_veh1 
##         92.25000         42.25000         11.75000        330.75000 
##   b_inc1030_veh0   b_inc3060_veh2   b_inc3060_veh1   b_inc3060_veh0 
##       -150.25000        332.00000        261.00000       -502.00000 
##  b_inc60100_veh1  b_inc60100_veh0 
##       -177.75000       -643.75000 
## initial  value 14248.333444 
## iter   2 value 12227.941889
## iter   3 value 12058.141183
## iter   4 value 11973.668443
## iter   5 value 11787.824685
## iter   6 value 11627.504735
## iter   7 value 11466.245417
## iter   8 value 11434.282848
## iter   9 value 11400.840325
## iter  10 value 11375.511280
## iter  11 value 11346.398967
## iter  12 value 11316.210466
## iter  13 value 11300.979497
## iter  14 value 11257.279339
## iter  15 value 11216.255733
## iter  16 value 11183.256269
## iter  17 value 11073.182585
## iter  18 value 10990.647270
## iter  19 value 10968.864338
## iter  20 value 10954.186065
## iter  21 value 10928.740821
## iter  22 value 10918.270067
## iter  23 value 10899.532190
## iter  24 value 10821.465758
## iter  25 value 10819.919001
## iter  26 value 10789.512055
## iter  27 value 10773.261402
## iter  28 value 10765.802585
## iter  29 value 10765.000379
## iter  30 value 10764.949772
## iter  31 value 10764.942332
## iter  32 value 10764.936238
## iter  33 value 10764.934516
## iter  34 value 10764.933895
## iter  35 value 10764.933702
## iter  35 value 10764.933642
## iter  35 value 10764.933640
## final  value 10764.933640 
## converged
## Additional convergence test using scaled estimation. Parameters will be
##   scaled by their current estimates and additional iterations will be
##   performed.
## initial  value 10764.933640 
## iter   1 value 10764.933639
## final  value 10764.933639 
## converged
## Estimated parameters:
##                     Estimate
## asc_veh3              0.0000
## asc_veh2              0.3378
## asc_veh1             -1.3170
## asc_veh0             -2.3433
## b_wrktodrv_veh1       0.3571
## b_wrktodrv_veh0      -2.1606
## b_615todrv_veh2       0.3817
## b_615todrv_veh1      -0.1599
## b_1824todrv_veh2     -2.0673
## b_1824todrv_veh1     -2.4203
## b_1824todrv_veh0     -1.7042
## b_80todrv_veh1        0.4839
## b_inc10_veh2          1.2573
## b_inc10_veh1          4.0763
## b_inc10_veh0          5.2541
## b_inc1030_veh2        0.6774
## b_inc1030_veh1        2.8862
## b_inc1030_veh0        3.3790
## b_inc3060_veh2        0.4068
## b_inc3060_veh1        1.7835
## b_inc3060_veh0        1.0693
## b_inc60100_veh1       0.5179
## b_inc60100_veh0      -0.7604
## 
## Final LL: -10764.9336
## 
## 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: -3.72853
## 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 (base model)
## Model run at                                : 2023-11-20 12:28:00
## 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)                : -12317.2
## LL(final)                                   : -10764.93
## Rho-squared vs equal shares                  :  0.2445 
## Adj.Rho-squared vs equal shares              :  0.2429 
## Rho-squared vs observed shares               :  0.126 
## Adj.Rho-squared vs observed shares           :  0.1242 
## AIC                                         :  21573.87 
## BIC                                         :  21733.1 
## 
## Estimated parameters                        :  22
## Time taken (hh:mm:ss)                       :  00:00:5.28 
##      pre-estimation                         :  00:00:1.74 
##      estimation                             :  00:00:1.15 
##      post-estimation                        :  00:00:2.4 
## Iterations                                  :  39  
## Min abs eigenvalue of Hessian               :  3.72853 
## 
## 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.3378     0.03128      10.799     0.03061        11.033
## asc_veh1             -1.3170     0.07571     -17.395     0.08059       -16.342
## asc_veh0             -2.3433     0.20022     -11.704     0.23622        -9.920
## b_wrktodrv_veh1       0.3571     0.06427       5.556     0.07028         5.080
## b_wrktodrv_veh0      -2.1606     0.17838     -12.112     0.23294        -9.275
## b_615todrv_veh2       0.3817     0.05882       6.490     0.05686         6.713
## b_615todrv_veh1      -0.1599     0.07630      -2.095     0.08405        -1.902
## b_1824todrv_veh2     -2.0673     0.16023     -12.902     0.17564       -11.770
## b_1824todrv_veh1     -2.4203     0.18646     -12.981     0.23465       -10.314
## b_1824todrv_veh0     -1.7042     0.32380      -5.263     0.41855        -4.072
## b_80todrv_veh1        0.4839     0.08787       5.507     0.09094         5.321
## b_inc10_veh2          1.2573     0.29383       4.279     0.31558         3.984
## b_inc10_veh1          4.0763     0.27892      14.615     0.30222        13.488
## b_inc10_veh0          5.2541     0.33837      15.528     0.36886        14.244
## b_inc1030_veh2        0.6774     0.10383       6.524     0.10831         6.254
## b_inc1030_veh1        2.8862     0.11019      26.193     0.11533        25.026
## b_inc1030_veh0        3.3790     0.22197      15.223     0.24499        13.793
## b_inc3060_veh2        0.4068     0.06558       6.204     0.06637         6.129
## b_inc3060_veh1        1.7835     0.08106      22.003     0.08172        21.823
## b_inc3060_veh0        1.0693     0.24871       4.299     0.26031         4.108
## b_inc60100_veh1       0.5179     0.07229       7.164     0.07193         7.200
## b_inc60100_veh0      -0.7604     0.36574      -2.079     0.36535        -2.081

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_615todrv_veh1 = Coefficient for ratio of people age 6-15 to drivers (HH with one vehicle), b_1824todrv_veh2 = Coefficient for ratio of people age 18-24 to drivers (HH with two or more vehicles), b_1824todrv_veh1 = Coefficient for ratio of people age 18-24 to drivers to drivers (HH with one vehicle), b_1824todrv_veh0 = Coefficient for ratio of people age 18-24 to drivers to drivers (HH with zero vehicle), 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_inc60100_veh1 = Coefficient for $60,000 < household income <= $100,000 (HH with one vehicle), b_inc60100_veh0 = Coefficient for $60,000 < household income <= $100,000 (zero vehicle HH)