R is a programming language for statistical computing and graphics. To install R, required installation files can be found in the following link:
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.
rm(list = ls())
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")
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
)
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
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 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 = 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...
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)