WTP distribution in preference space

[kkavta]

Dear Prof. Hess,

I have developed a mixed logit model and want to calculate 'Value of Risk Reduction (VRR)' from the two random parameters used in the model. First parameter is travel time (negative log normal distribution) and second is 'Safety information' (positive log normal). I have the following questions:

1) Would the distribution of VRR also be negative log normal? Where, VRR= unconditional(b_safety)/unconditional(b_tt)

2) I have plotted the density for both random variables and VRR using the below code. However, the plot doesn't resemble log normal distribution (Pics attached). What could be the reason for this?

plot(density(unconditionals[["b_tt"]]),col=1,main="Travel time",xlab="beta", ylab="density", xlim = c(-4,0.5), ylim=c(0,17))
legend(-2,4,c("traveltime"),col=c(1),lty=1,cex=1)

plot(density(unconditionals[["b_info"]]),col=2,main="safety info",xlab="beta", ylab="density", xlim = c(0,4), ylim=c(-0.5,15))
legend(2,15,c("safety_info"),col=c(2),lty=1,cex=1)

plot(density(VRR),col=3,main="VRR distribution",xlab="VRR per 10% reduced crashes",ylab="density",xlim = c(-40,3), ylim=c(0,5))
legend(-30,4,c("VRR"),col=c(3),lty=1,cex=1)

3) How can I get summary of VRR values? I have tried this, but have doubts if it is correct. The code used and the result obtained are the following:

code:
VRR = (unconditionals[["b_info"]]/unconditionals[["b_tt"]])
(mean(VRR)); (sd(VRR))

summary((conditionals[["b_info"]]/conditionals[["b_tt"]]))

results:
VRR = (unconditionals[["b_info"]]/unconditionals[["b_tt"]])
> (mean(VRR)); (sd(VRR))
[1] -7.282214
[1] 47.02358
>
> summary((conditionals[["b_info"]]/conditionals[["b_tt"]]))
ID post.mean post.sd
Min. :1 Min. :-12.74364 Min. : 0.04947
1st Qu.:1 1st Qu.: -5.31435 1st Qu.: 0.12425
Median :1 Median : -1.21677 Median : 1.39308
Mean :1 Mean : -2.93442 Mean : 3.75905
3rd Qu.:1 3rd Qu.: -0.20944 3rd Qu.: 7.78408
Max. :1 Max. : -0.06611 Max. :13.80384


Thank you for your support.
K






-----------X------------X-------------X-------------X

Complete output of model:
# ################################################################# #
> #### LOAD LIBRARY AND DEFINE CORE SETTINGS ####
> # ################################################################# #
>
> ### Clear memory
> rm(list = ls())
>
> ### Load Apollo library
> library(apollo)
>
> ### Initialise code
> apollo_initialise()
Apollo ignition sequence completed
>
> ### Set core controls
> apollo_control = list(
+ modelName = "E_info_mix_pref4",
+ modelDescr = "Mixed model (TT&info) for info in preference space",
+ indivID = "ID",
+ nCores = 12,
+ outputDirectory = "out_info_mix_pref"
+ )
>
> # ################################################################# #
> #### LOAD DATA AND APPLY ANY TRANSFORMATIONS ####
> # ################################################################# #
> database = read.csv("E_use_070224.csv",header=TRUE)
>
> # ################################################################# #
> #### DEFINE MODEL PARAMETERS ####
> # ################################################################# #
> ### Vector of parameters, including any that are kept fixed in estimation
> apollo_beta=c(asc_sco = 0,
+ asc_alt = 0,
+ b_ageL34 = 0,
+ b_nonEUcitizen = 0,
+ b_timepres_yes = 0,
+ b_expL1 = 0,
+ mu_log_b_tt = -3,
+ sigma_log_b_tt = -1,
+ mu_log_b_info = 2,
+ sigma_log_b_info = 2)
>
> ### 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_alt")
>
> ### Set parameters for generating draws
> apollo_draws = list(
+ interDrawsType = "halton",
+ interNDraws = 50000,
+ interUnifDraws = c(),
+ interNormDraws = c("eta_tt","eta_info"),
+ intraDrawsType = "halton",
+ intraNDraws = 0,
+ intraUnifDraws = c(),
+ intraNormDraws = c()
+ )
>
> ### Create random parameters
> apollo_randCoeff = function(apollo_beta, apollo_inputs){
+ randcoeff = list()
+
+ randcoeff[["b_tt"]] = -exp( mu_log_b_tt + sigma_log_b_tt * eta_tt )
+ randcoeff[["b_info"]] = exp(mu_log_b_info + sigma_log_b_info * eta_info)
+
+ return(randcoeff)
+ }
>
> # ################################################################# #
> #### GROUP AND VALIDATE INPUTS ####
> # ################################################################# #
> apollo_inputs = apollo_validateInputs()
apollo_draws and apollo_randCoeff were found, so apollo_control$mixing was set to TRUE
Several observations per individual detected based on the value of ID. Setting panelData in
apollo_control set to TRUE.
All checks on apollo_control completed.
All checks on database completed.
Generating inter-individual draws .. Done
>
> # ################################################################# #
> #### 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[["scoober"]] = asc_sco + b_tt * time_sco + b_info * (saf_info/10) +
+ b_nonEUcitizen * (citizen == 3) + b_ageL34 * (age == 1|age==2|age ==3) + b_timepres_yes * (time_pres == 2) +
+ b_expL1 * (job_exp == 1)
+ V[["alt"]] = asc_alt + b_tt * time_alt
+ ### Define settings for MNL model component
+ mnl_settings = list(
+ alternatives = c(scoober=1, alt=2),
+ avail = list(scoober=1, alt=1),
+ choiceVar = choice_info,
+ 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)
+
+ ### Average across inter-individual draws
+ P = apollo_avgInterDraws(P, apollo_inputs, functionality)
+
+ ### Prepare and return outputs of function
+ P = apollo_prepareProb(P, apollo_inputs, functionality)
+ return(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 :
scoober alt
Times available 606.00 606.00
Times chosen 322.00 284.00
Percentage chosen overall 53.14 46.86
Percentage chosen when available 53.14 46.86


Pre-processing likelihood function...
Creating cluster...
Preparing workers for multithreading...

Testing influence of parameters
Starting main estimation

BGW using analytic model derivatives supplied by caller...


Iterates will be written to:
out_info_mix_pref/E_info_mix_pref4_iterations.csv it nf F RELDF PRELDF RELDX MODEL stppar
0 1 6.005727432e+02
1 3 5.198353528e+02 1.344e-01 8.955e-02 4.18e-02 G 8.51e-01
2 4 4.215050586e+02 1.892e-01 1.089e-01 2.57e-01 G 7.92e-03
3 5 3.453171883e+02 1.808e-01 2.496e-01 5.77e-01 S 0.00e+00
4 7 3.328778312e+02 3.602e-02 3.234e-02 1.41e-01 S 4.07e-01
5 8 3.287079915e+02 1.253e-02 1.447e-02 7.90e-02 S 2.13e-01
6 9 3.256073046e+02 9.433e-03 9.400e-03 1.21e-01 G 3.63e-02
7 10 3.250042870e+02 1.852e-03 1.642e-03 7.30e-02 G 0.00e+00
8 11 3.249473422e+02 1.752e-04 1.655e-04 1.84e-02 G 0.00e+00
9 12 3.249370258e+02 3.175e-05 2.642e-05 9.41e-03 G 0.00e+00
10 13 3.249346864e+02 7.199e-06 6.426e-06 4.41e-03 G 0.00e+00
11 14 3.249340743e+02 1.884e-06 1.557e-06 2.03e-03 G 0.00e+00
12 15 3.249339032e+02 5.267e-07 6.396e-07 2.42e-03 S 0.00e+00
13 16 3.249338878e+02 4.746e-08 6.566e-08 3.72e-04 S 0.00e+00
14 17 3.249338840e+02 1.149e-08 1.098e-08 1.18e-04 S 0.00e+00
15 18 3.249338839e+02 4.288e-10 3.356e-10 2.38e-05 S 0.00e+00
16 19 3.249338839e+02 4.533e-11 4.445e-11 1.76e-05 S 0.00e+00

***** Relative function convergence *****
Additional convergence test using scaled estimation. Parameters will be scaled by their
current estimates and additional iterations will be performed.

BGW using analytic model derivatives supplied by caller...


Iterates will be appended to:
out_info_mix_pref/E_info_mix_pref4_iterations.csv it nf F RELDF PRELDF RELDX MODEL stppar
0 1 3.249338839e+02
1 2 3.249338839e+02 1.494e-13 1.117e-13 1.10e-06 G 0.00e+00

***** Relative function convergence *****

Estimated parameters with approximate standard errors from BHHH matrix:
Estimate BHHH se BHH t-ratio (0)
asc_sco 1.601 0.5901 2.712
asc_alt 0.000 NA NA
b_ageL34 -1.529 0.5813 -2.631
b_nonEUcitizen 1.455 0.6047 2.406
b_timepres_yes -1.153 0.4948 -2.331
b_expL1 1.122 0.5730 1.958
mu_log_b_tt -1.410 0.2521 -5.592
sigma_log_b_tt -1.553 0.2963 -5.242
mu_log_b_info -1.371 0.3804 -3.603
sigma_log_b_info 1.217 0.5276 2.308

Final LL: -324.9339

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: -1.239554
Computing score matrix...

Your model was estimated using the BGW algorithm. Please acknowledge this by citing Bunch
et al. (1993) - DOI 10.1145/151271.151279
>
> # ################################################################# #
> #### MODEL OUTPUTS ####
> # ################################################################# #
>
> apollo_modelOutput(model, list(printPVal = 2))
Model run by kkavta using Apollo 0.3.1 on R 4.3.1 for Darwin.
Please acknowledge the use of Apollo by citing Hess & Palma (2019)
DOI 10.1016/j.jocm.2019.100170
www.ApolloChoiceModelling.com

Model name : E_info_mix_pref4
Model description : Mixed model (TT&info) for info in preference space
Model run at : 2024-02-20 18:39:20.756729
Estimation method : bgw
Model diagnosis : Relative function convergence
Optimisation diagnosis : Maximum found
hessian properties : Negative definite
maximum eigenvalue : -1.239554
reciprocal of condition number : 0.0123507
Number of individuals : 202
Number of rows in database : 606
Number of modelled outcomes : 606

Number of cores used : 12
Number of inter-individual draws : 50000 (halton)

LL(start) : -600.57
LL at equal shares, LL(0) : -420.05
LL at observed shares, LL(C) : -418.85
LL(final) : -324.93
Rho-squared vs equal shares : 0.2264
Adj.Rho-squared vs equal shares : 0.205
Rho-squared vs observed shares : 0.2242
Adj.Rho-squared vs observed shares : 0.2051
AIC : 667.87
BIC : 707.53

Estimated parameters : 9
Time taken (hh:mm:ss) : 00:06:58.87
pre-estimation : 00:01:23.36
estimation : 00:01:30.23
initial estimation : 00:01:15.44
estimation after rescaling : 00:00:14.79
post-estimation : 00:04:5.28
Iterations : 17
initial estimation : 16
estimation after rescaling : 1

Unconstrained optimisation.

Estimates:
Estimate s.e. t.rat.(0) p(2-sided) Rob.s.e. Rob.t.rat.(0) p(2-sided)
asc_sco 1.601 0.7504 2.133 0.032931 0.9957 1.608 0.107927
asc_alt 0.000 NA NA NA NA NA NA
b_ageL34 -1.529 0.6597 -2.318 0.020436 0.7981 -1.916 0.055330
b_nonEUcitizen 1.455 0.6152 2.365 0.018022 0.6642 2.191 0.028474
b_timepres_yes -1.153 0.5611 -2.056 0.039815 0.6722 -1.716 0.086223
b_expL1 1.122 0.5960 1.882 0.059821 0.6572 1.707 0.087836
mu_log_b_tt -1.410 0.2623 -5.375 7.667e-08 0.2882 -4.892 9.998e-07
sigma_log_b_tt -1.553 0.2702 -5.747 9.058e-09 0.2755 -5.638 1.721e-08
mu_log_b_info -1.371 0.3833 -3.576 3.4884e-04 0.4005 -3.423 6.2032e-04
sigma_log_b_info 1.217 0.4230 2.878 0.003997 0.3908 3.116 0.001836

> # ----------------------------------------------------------------- #
> #---- FORMATTED OUTPUT (TO FILE, using model name) ----
> # ----------------------------------------------------------------- #
>
> apollo_saveOutput(model, list(printPVal = 2))
Model output saved to out_info_mix_pref/E_info_mix_pref4_output.txt
Estimates saved to out_info_mix_pref/E_info_mix_pref4_estimates.csv
Model object saved to out_info_mix_pref/E_info_mix_pref4.rds
>
> # ################################################################# #
> ##### POST-PROCESSING ####
> # ################################################################# #
>
> ### Print outputs of additional diagnostics to new output file (remember to close file writing when complete)
> apollo_sink()
Writing output to file out_info_mix_pref/E_info_mix_pref4_additional_output.txt. Please run
"apollo_sink()" again after finishing writing results.
>
> # ----------------------------------------------------------------- #
> #---- CONDITIONALS AND UNCONDITIONALS ----
> # ----------------------------------------------------------------- #
> unconditionals <- apollo_unconditionals(model,apollo_probabilities,apollo_inputs)
Unconditional distributions computed
>
> conditionals <- apollo_conditionals(model,apollo_probabilities,apollo_inputs)
Calculating conditionals...
>
> mean(unconditionals[["b_tt"]])
[1] -0.8156153
> sd(unconditionals[["b_tt"]])
[1] 2.592688
> summary(conditionals[["b_tt"]])
ID post.mean post.sd
Min. : 1.00 Min. :-4.05688 Min. :0.06473
1st Qu.: 51.25 1st Qu.:-1.50391 1st Qu.:0.14660
Median :101.50 Median :-0.30396 Median :0.33052
Mean :101.50 Mean :-0.82090 Mean :1.43543
3rd Qu.:151.75 3rd Qu.:-0.13387 3rd Qu.:3.50985
Max. :202.00 Max. :-0.07557 Max. :6.88311
>
> mean(unconditionals[["b_info"]])
[1] 0.5328426
> sd(unconditionals[["b_info"]])
[1] 0.9811726
> summary(conditionals[["b_info"]])
ID post.mean post.sd
Min. : 1.00 Min. :0.1037 Min. :0.08147
1st Qu.: 51.25 1st Qu.:0.2746 1st Qu.:0.32648
Median :101.50 Median :0.4326 Median :0.57204
Mean :101.50 Mean :0.5335 Mean :0.76278
3rd Qu.:151.75 3rd Qu.:0.7318 3rd Qu.:1.27419
Max. :202.00 Max. :1.5186 Max. :2.07384
>
> VRR = (unconditionals[["b_info"]]/unconditionals[["b_tt"]])
> (mean(VRR)); (sd(VRR))
[1] -7.282214
[1] 47.02358
>
> summary((conditionals[["b_info"]]/conditionals[["b_tt"]]))
ID post.mean post.sd
Min. :1 Min. :-12.74364 Min. : 0.04947
1st Qu.:1 1st Qu.: -5.31435 1st Qu.: 0.12425
Median :1 Median : -1.21677 Median : 1.39308
Mean :1 Mean : -2.93442 Mean : 3.75905
3rd Qu.:1 3rd Qu.: -0.20944 3rd Qu.: 7.78408
Max. :1 Max. : -0.06611 Max. :13.80384
>
> plot(density(VRR),col=3,main="VRR distribution",xlab="VRR per 10% reduced crashes",ylab="density",xlim = c(-40,3), ylim=c(0,5))
> legend(-30,4,c("VRR"),col=c(3),lty=1,cex=1)
>
> plot(density(unconditionals[["b_tt"]]),col=1,main="Travel time",xlab="beta", ylab="density", xlim = c(-4,0.5), ylim=c(0,17))
> legend(-2,4,c("traveltime"),col=c(1),lty=1,cex=1)
>
> plot(density(unconditionals[["b_info"]]),col=2,main="safety info",xlab="beta", ylab="density", xlim = c(0,4), ylim=c(-0.5,15))
> legend(2,15,c("safety_info"),col=c(2),lty=1,cex=1)
>
> # ----------------------------------------------------------------- #
> #---- switch off writing to file ----
> # ----------------------------------------------------------------- #
>
> if(sink.number()>0) sink()



Plots:

[stephanehess]

Hi

Yes, a ratio of two lognormals is a lognormals, and if one (and only one) of the two is a negative lognormal, then the ratio is negative lognormal too

Re your plots, are you using a low number of draws? Maybe also try without imposing the limits

Stephane