A quick introdution:
I work with transport studies and there's this software we use, PTV Visum, which provides a variety of choice models to distribute travel demand along paths stochastically (route choice). First we setup the utility of the routes (for example, B0 * time + B1 * distance + B2 * toll) and then we setup which choice model should be used to calculate routes probabilities.
There's a choice model in the software called "Box-cox" which basically applies a Box-Cox Transformation over the Utility of the route before calculating the probabilities. Please note that we must insert just one theta (Box-Cox parameter) when setting up the process, since the transformation is applied over the utility itself, not over the independent variables one by one (time, distance etc.)
The problem:
I have some data from DCE (unlabeled alternatives, route A vs route B) and I'd like to calibrate a model (coefficients and theta) which would be the same one used by Visum in its choice model (to ensure consistency in this process of obtaining and inserting parameters in the software). I was told I could use Apollo to do that, since the package allows free specification of utilities and probabilities.
Right now I'm facing convergence problems. I believe this is due to the fact that the Box-Cox Transformation don't accept/work with negative numbers and perhaps during the convergence process some utilities are returning negative values. Giving more details, since the utilities for route choice are tipically negative (time, cost and distance all should return negative coefficients) I've already built the model as V = - f(V|theta), so that coefficients would take positive values. But there are also categorical variables in the experiment, like road type (twolane vs. multilane), and "multilane" tipically increases the utility (which means, in my "inverted model", that its coefficient is negative). So I believe that, for some few specific choice tasks in my database, some negative utilities may still appear in the middle of the convergence process.
This is my code:
Code: Select all
apollo_initialise()
### Set core controls
apollo_control = list(
modelName = "MNL_SP",
modelDescr = "RS-BNDES Route Choice",
indivID = "id",
outputDirectory = "output",
workInLogs=TRUE
)
### Vector of parameters, including any that are kept fixed in estimation
apollo_beta=c(asc_alt1 = 0,
asc_alt2 = 0,
b_dist = 0.01,
b_tempo = 0.01,
b_custo = 0.01,
b_conc = -0.01,
b_dupla = -0.01,
b_terra = 0.01,
b_concdupla = -0.01,
bc_theta = 0.5
)
### 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_fixed = c("asc_alt2", "bc_theta")
apollo_inputs = apollo_validateInputs()
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[["alt1"]] = -((((asc_alt1 + b_dist * dv1 + b_tempo * tv1 + b_custo * cv1 + b_conc * dummy.conc1 + b_dupla * dummy.dupla1 + b_terra * dummy.terra1 + b_concdupla * dummy.concdupla1)^bc_theta)-1)/bc_theta)
V[["alt2"]] = -((((asc_alt2 + b_dist * dv2 + b_tempo * tv2 + b_custo * cv2 + b_conc * dummy.conc2 + b_dupla * dummy.dupla2 + b_terra * dummy.terra2 + b_concdupla * dummy.concdupla2)^bc_theta)-1)/bc_theta)
### Define settings for MNL model component
mnl_settings = list(
alternatives = c(alt1=0, alt2=1),
avail = 1,
choiceVar = Choice,
utilities = V
)
### Compute probabilities using MNL model
P[["model"]] = apollo_mnl(mnl_settings, functionality)
### Take product across observation for same individual
P = apollo_panelProd(P, apollo_inputs, functionality)
### Prepare and return outputs of function
P = apollo_prepareProb(P, apollo_inputs, functionality)
return(P)
}
estimate_settings = list(
bootstrapSE=0
)
model = apollo_estimate(apollo_beta, apollo_fixed, apollo_probabilities, apollo_inputs, estimate_settings)
modelOutput_settings = list(
printPVal=1
)
apollo_modelOutput(model, modelOutput_settings)
Code: Select all
WARNING: Element bc_theta in 'apollo_fixed' is constrained to a value other than zero or one. This may be intentional. If not, stop this function by pressing the "Escape"
key and adjust the starting values accordingly.
Current process will resume in 5 seconds unless interrupted by the user.....
Preparing user-defined functions.
Testing likelihood function...
Overview of choices for MNL model component :
alt1 alt2
Times available 7488.00 7488.00
Times chosen 3507.00 3981.00
Percentage chosen overall 46.83 53.17
Percentage chosen when available 46.83 53.17
Pre-processing likelihood function...
Testing influence of parameters........
Starting main estimation
BGW is using FD derivatives for model Jacobian. (Caller did not provide derivatives.)
Iterates will be written to:
output/MNL_SP_iterations.csv it nf F RELDF PRELDF RELDX MODEL stppar
0 1 5.623234851e+03
1 4 5.373310894e+03 4.444e-02 4.312e-02 9.29e-02 G 2.46e+00
2 7 5.127290474e+03 4.579e-02 4.327e-02 1.13e-01 G 1.27e+01
3 11 5.051877012e+03 1.471e-02 1.468e-02 2.41e-02 S 1.94e+01
4 14 5.016548615e+03 6.993e-03 6.989e-03 1.22e-02 S 2.94e+01
5 18 5.007878231e+03 1.728e-03 1.725e-03 3.25e-03 S 9.36e+01
6 22 5.005690692e+03 4.368e-04 4.363e-04 9.45e-04 S 2.69e+02
7 25 5.004556737e+03 2.265e-04 2.264e-04 5.00e-04 S 3.94e+02
8 28 5.003956862e+03 1.199e-04 1.198e-04 2.39e-04 S 5.67e+02
9 30 5.003324488e+03 1.264e-04 1.264e-04 2.08e-04 S 4.53e+02
10 33 5.002920567e+03 8.073e-05 8.073e-05 6.78e-05 S 6.15e+02
11 35 5.002491791e+03 8.570e-05 8.571e-05 4.06e-05 S 5.98e+02
12 37 5.001594215e+03 1.794e-04 1.795e-04 4.18e-05 S 9.13e+01
13 41 4.820028613e+03 3.630e-02 3.649e-02 5.29e-03 S 5.79e-01
14 43 4.730828338e+03 1.851e-02 1.564e-02 7.40e-03 S 1.08e-01
15 44 4.702528277e+03 5.982e-03 5.986e-03 5.32e-03 G 1.95e-02
16 47 4.689024603e+03 2.872e-03 2.889e-03 3.52e-02 G 3.12e-03
17 49 4.614805183e+03 1.583e-02 1.525e-02 2.40e-01 G 8.64e-05
18 57 4.611472864e+03 7.221e-04 7.727e-04 1.28e-02 G 8.69e-02
19 59 4.607199623e+03 9.267e-04 9.231e-04 1.67e-02 S 2.81e-02
20 61 4.580590442e+03 5.776e-03 5.823e-03 1.09e-01 S 8.41e-03
21 65 4.565872610e+03 3.213e-03 3.338e-03 8.06e-02 S 1.39e-01
22 68 4.562392852e+03 7.621e-04 1.003e-03 2.23e-02 G 1.30e-01
23 70 4.560316377e+03 4.551e-04 4.252e-04 2.66e-02 S 2.05e-02
24 77 4.560014857e+03 6.612e-05 5.457e-05 1.80e-03 G 2.55e+00
25 85 4.559997754e+03 3.751e-06 3.095e-06 2.31e-05 G 7.90e+01
26 88 4.559984315e+03 2.947e-06 2.355e-06 1.07e-05 G 1.23e+02
27 91 4.559971161e+03 2.885e-06 1.918e-06 4.71e-06 G 1.66e+02
28 97 4.559969062e+03 4.603e-07 3.575e-07 2.32e-07 G 9.85e+02
29 100 4.559967228e+03 4.023e-07 2.625e-07 1.04e-07 G 1.36e+03
30 105 4.559966752e+03 1.043e-07 6.073e-08 1.06e-08 G 5.95e+03
31 108 4.559966360e+03 8.599e-08 3.872e-08 4.66e-09 G 9.35e+03
32 112 4.559966182e+03 3.897e-08 1.264e-08 9.77e-10 G 2.87e+04
33 115 4.559966018e+03 3.596e-08 7.275e-09 4.33e-10 G 4.99e+04
34 120 4.559965977e+03 9.153e-09 1.051e-09 4.72e-11 G 3.45e+05
35 123 4.559965944e+03 7.063e-09 5.449e-10 2.27e-11 G 6.66e+05
36 126 4.559965915e+03 6.369e-09 2.803e-10 1.10e-11 G 1.29e+06
37 130 4.559965895e+03 4.511e-09 7.192e-11 2.67e-12 G 5.05e+06
38 138 4.559965893e+03 3.963e-10 1.145e-12 4.09e-14 G 3.17e+08
39 142 4.559965892e+03 2.531e-10 2.866e-13 1.02e-14 G 1.27e+09
40 148 4.559965891e+03 7.157e-11 1.793e-14 6.60e-16 G 2.02e+10
41 156 4.559965891e+03 4.471e-12 2.803e-16 9.92e-18 G 1.30e+12
42 159 4.559965891e+03 6.106e-12 1.412e-16 4.96e-18 G 2.61e+12
***** Singular convergence *****
Estimated parameters:
Estimate
asc_alt1 -6.6239e-04
asc_alt2 0.000000
b_dist -4.288e-05
b_tempo 0.378654
b_custo 0.002289
b_conc 0.055112
b_dupla -0.123840
b_terra 0.897965
b_concdupla -0.055922
bc_theta 0.500000
Final LL: -4559.9659
WARNING: Estimation failed. No covariance matrix to compute.
Current process will resume in 3 seconds unless interrupted by the user...
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
Your model was estimated using the BGW algorithm. Please acknowledge this by citing Bunch et al. (1993) - DOI 10.1145/151271.151279
Model run by gabriel.souza using Apollo 0.3.1 on R 4.2.1 for Windows.
Please acknowledge the use of Apollo by citing Hess & Palma (2019)
DOI 10.1016/j.jocm.2019.100170
www.ApolloChoiceModelling.com
Model name : MNL_SP
Model description : RS-BNDES Route Choice
Model run at : 2023-12-14 19:54:55
Estimation method : bgw
Model diagnosis : Singular convergence
Number of individuals : 1248
Number of rows in database : 7488
Number of modelled outcomes : 7488
Number of cores used : 1
Model without mixing
LL(start) : -5623.23
LL at equal shares, LL(0) : -5190.29
LL at observed shares, LL(C) : -5175.27
LL(final) : -4559.97
Rho-squared vs equal shares : 0.1214
Adj.Rho-squared vs equal shares : 0.1199
Rho-squared vs observed shares : 0.1189
Adj.Rho-squared vs observed shares : 0.1175
AIC : 9135.93
BIC : 9191.3
Estimated parameters : 8
Time taken (hh:mm:ss) : 00:00:11.26
pre-estimation : 00:00:5.54
estimation : 00:00:5.67
post-estimation : 00:00:0.04
Iterations : 43 (Singular convergence)
Unconstrained optimisation.
Estimates:
Estimate s.e. t.rat.(0) p(1-sided) Rob.s.e. Rob.t.rat.(0) p(1-sided)
asc_alt1 -6.6239e-04 NA NA NA NA NA NA
asc_alt2 0.000000 NA NA NA NA NA NA
b_dist -4.288e-05 NA NA NA NA NA NA
b_tempo 0.378654 NA NA NA NA NA NA
b_custo 0.002289 NA NA NA NA NA NA
b_conc 0.055112 NA NA NA NA NA NA
b_dupla -0.123840 NA NA NA NA NA NA
b_terra 0.897965 NA NA NA NA NA NA
b_concdupla -0.055922 NA NA NA NA NA NA
bc_theta 0.500000 NA NA NA NA NA NA
The most simple thing I could think about to overcome this problem was to insert a "maximum condition" over the utility "pre-transformed" (if it's <0, returns 0). I tried using max() from base R inside the V specification in Apollo but it didn't work.
And just be clear, this "possible negative value in utility" is just a suspicion. I'm not sure that's the cause of the error.
I appreciate very much if anyone could help me on this one.