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 (random 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_random.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.877   Mean   :0.1064   Mean   :0.2353  
##  3rd Qu.: 7709   3rd Qu.:4.000   3rd Qu.:0.0000   3rd Qu.:0.0000  
##  Max.   :10278   Max.   :4.000   Max.   :1.0000   Max.   :1.0000  
##       veh2             veh3          wrktodrv      age615todrv    
##  Min.   :0.0000   Min.   :0.000   Min.   :0.000   Min.   :0.0000  
##  1st Qu.:0.0000   1st Qu.:0.000   1st Qu.:0.500   1st Qu.:0.0000  
##  Median :0.0000   Median :0.000   Median :1.000   Median :0.0000  
##  Mean   :0.3333   Mean   :0.325   Mean   :0.687   Mean   :0.1892  
##  3rd Qu.:1.0000   3rd Qu.:1.000   3rd Qu.:1.000   3rd Qu.:0.0000  
##  Max.   :1.0000   Max.   :1.000   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 
## 1094 2418 3426 3340

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_veh0     = 0,
              b_615todrv_veh2     = 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_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.0
## Times chosen                      1094.00  2418.00  3426.00  3340.0
## Percentage chosen overall           10.64    23.53    33.33    32.5
## Percentage chosen when available    10.64    23.53    33.33    32.5
## 
## 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_veh0 
##       856.500000      -151.500000     -1475.500000     -1197.832143 
##  b_615todrv_veh2 b_1824todrv_veh2 b_1824todrv_veh1 b_1824todrv_veh0 
##       277.344643        -3.427381       -25.710714       -84.677381 
##   b_80todrv_veh1     b_inc10_veh2     b_inc10_veh1     b_inc10_veh0 
##        37.787500       -24.750000        37.250000        66.250000 
##   b_inc1030_veh2   b_inc1030_veh1   b_inc1030_veh0   b_inc3060_veh2 
##        -5.250000       164.750000       -32.250000       143.000000 
##   b_inc3060_veh1   b_inc3060_veh0  b_inc60100_veh1  b_inc60100_veh0 
##       135.000000      -314.000000      -175.750000      -465.750000 
## initial  value 14248.333444 
## iter   2 value 13455.116794
## iter   3 value 13254.844582
## iter   4 value 13211.685064
## iter   5 value 13094.669148
## iter   6 value 12999.749602
## iter   7 value 12938.003327
## iter   8 value 12910.151983
## iter   9 value 12882.476386
## iter  10 value 12857.715221
## iter  11 value 12845.390471
## iter  12 value 12834.974511
## iter  13 value 12825.559600
## iter  14 value 12819.301252
## iter  15 value 12809.730226
## iter  16 value 12793.648454
## iter  17 value 12761.585839
## iter  18 value 12745.817958
## iter  19 value 12734.904398
## iter  20 value 12731.921296
## iter  21 value 12731.770058
## iter  22 value 12731.306044
## iter  23 value 12729.699959
## iter  24 value 12728.184005
## iter  25 value 12717.390090
## iter  26 value 12717.015585
## iter  27 value 12716.976055
## iter  28 value 12716.973412
## iter  28 value 12716.973262
## iter  28 value 12716.973260
## final  value 12716.973260 
## converged
## Additional convergence test using scaled estimation. Parameters will be
##   scaled by their current estimates and additional iterations will be
##   performed.
## initial  value 12716.973260 
## iter   1 value 12716.973259
## final  value 12716.973259 
## converged
## Estimated parameters:
##                     Estimate
## asc_veh3             0.00000
## asc_veh2            -0.06776
## asc_veh1            -0.86894
## asc_veh0            -1.68859
## b_wrktodrv_veh0     -0.39984
## b_615todrv_veh2      0.22736
## b_1824todrv_veh2    -1.56179
## b_1824todrv_veh1    -1.57418
## b_1824todrv_veh0    -2.00075
## b_80todrv_veh1       0.32270
## b_inc10_veh2         1.37416
## b_inc10_veh1         2.76353
## b_inc10_veh0         3.91279
## b_inc1030_veh2       0.62826
## b_inc1030_veh1       1.85043
## b_inc1030_veh0       2.39792
## b_inc3060_veh2       0.29759
## b_inc3060_veh1       1.09895
## b_inc3060_veh0       1.13072
## b_inc60100_veh1      0.16799
## b_inc60100_veh0      0.38093
## 
## Final LL: -12716.9733
## 
## 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: -7.280963
## 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 (random model)
## Model run at                                : 2023-11-20 12:38:25
## 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)                : -13467.87
## LL(final)                                   : -12716.97
## Rho-squared vs equal shares                  :  0.1075 
## Adj.Rho-squared vs equal shares              :  0.1061 
## Rho-squared vs observed shares               :  0.0558 
## Adj.Rho-squared vs observed shares           :  0.0543 
## AIC                                         :  25473.95 
## BIC                                         :  25618.7 
## 
## Estimated parameters                        :  20
## Time taken (hh:mm:ss)                       :  00:00:5.9 
##      pre-estimation                         :  00:00:2.06 
##      estimation                             :  00:00:1.21 
##      post-estimation                        :  00:00:2.63 
## Iterations                                  :  32  
## Min abs eigenvalue of Hessian               :  7.280963 
## 
## 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.00000          NA          NA          NA            NA
## asc_veh2            -0.06776     0.03160      -2.144     0.03128        -2.166
## asc_veh1            -0.86894     0.04790     -18.142     0.04782       -18.172
## asc_veh0            -1.68859     0.09668     -17.466     0.10164       -16.614
## b_wrktodrv_veh0     -0.39984     0.08156      -4.902     0.08861        -4.513
## b_615todrv_veh2      0.22736     0.04932       4.610     0.04853         4.685
## b_1824todrv_veh2    -1.56179     0.16388      -9.530     0.16763        -9.317
## b_1824todrv_veh1    -1.57418     0.17376      -9.060     0.19753        -7.969
## b_1824todrv_veh0    -2.00075     0.24952      -8.018     0.30780        -6.500
## b_80todrv_veh1       0.32270     0.08596       3.754     0.08551         3.774
## b_inc10_veh2         1.37416     0.23205       5.922     0.24028         5.719
## b_inc10_veh1         2.76353     0.22322      12.380     0.23175        11.925
## b_inc10_veh0         3.91279     0.23294      16.797     0.24105        16.233
## b_inc1030_veh2       0.62826     0.09280       6.770     0.09480         6.627
## b_inc1030_veh1       1.85043     0.09505      19.468     0.09692        19.092
## b_inc1030_veh0       2.39792     0.12136      19.759     0.12302        19.493
## b_inc3060_veh2       0.29759     0.06402       4.649     0.06431         4.627
## b_inc3060_veh1       1.09895     0.07395      14.860     0.07399        14.852
## b_inc3060_veh0       1.13072     0.10984      10.294     0.11046        10.236
## b_inc60100_veh1      0.16799     0.06800       2.471     0.06786         2.475
## b_inc60100_veh0      0.38093     0.10698       3.561     0.10706         3.558

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_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_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)