Page 1 of 1

CALCULATION ISSUE - Log-likelihood calculation fails at values close to the starting values!

Posted: 20 Jul 2023, 15:38
by thoscha
Hi there,

I am estimating a mode choice model with data from an SP experiment on mobility coins. First step was to estimate a simple MNL in WTP space with two components.
Model run by ge75mum using Apollo 0.2.9 on R 4.3.1 for Windows.
www.ApolloChoiceModelling.com

Model name : 4.1_mnl_wtp_space.R
Model description : mnl in wtp, No income individuals removed. Preliminary model. Includes shifts (multipliers) for purposes
Model run at : 2023-07-12 00:53:30.562108
Estimation method : bgw
Model diagnosis : Relative function convergence
Optimisation diagnosis : Maximum found
hessian properties : Negative definitive
maximum eigenvalue : -1.61998
Number of individuals : 1053
Number of rows in database : 11172
Number of modelled outcomes : 11172
mocono : 5448
mocosi : 5724

Number of cores used : 1
Model without mixing

LL(start) : -9359.71
LL (whole model) at equal shares, LL(0) : -12822.54
LL (whole model) at observed shares, LL(C) : -12146.8
LL(final, whole model) : -9349.34
Rho-squared vs equal shares : 0.2709
Adj.Rho-squared vs equal shares : 0.2664
Rho-squared vs observed shares : 0.2303
Adj.Rho-squared vs observed shares : 0.2261
AIC : 18812.68
BIC : 19189.05

LL(0,mocono) : -6283.17
LL(final,mocono) : -4602.63
LL(0,mocosi) : -6539.37
LL(final,mocosi) : -4746.71

Estimated parameters : 57
Time taken (hh:mm:ss) : 00:03:2.82
pre-estimation : 00:00:43.12
estimation : 00:00:8.31
initial estimation : 00:00:6.57
estimation after rescaling : 00:00:1.73
post-estimation : 00:02:11.39
Iterations : 11
initial estimation : 10
estimation after rescaling : 1

Unconstrained optimisation.

Estimates:
Estimate s.e. t.rat.(0) Rob.s.e. Rob.t.rat.(0)
asc_pt_base 0.000000 NA NA NA NA
asc_c_base -1.254459 0.374202 -3.3524 0.678908 -1.84776
asc_c_work 0.000000 NA NA NA NA
asc_c_errand 0.532274 0.194014 2.7435 0.325964 1.63293
asc_c_leisure -0.453211 0.197837 -2.2908 0.347810 -1.30304
asc_c_female 0.000000 NA NA NA NA
asc_c_male 0.299004 0.061825 4.8363 0.108138 2.76502
asc_c_divers 1.618516 0.632749 2.5579 0.811620 1.99418
asc_c_age 0.040300 0.012786 3.1519 0.023521 1.71337
asc_c_age_2 -4.4839e-04 1.3464e-04 -3.3302 2.4721e-04 -1.81377
asc_b_base -3.850031 0.319654 -12.0444 0.535656 -7.18751
asc_b_work 0.000000 NA NA NA NA
asc_b_errand 0.203590 0.149451 1.3623 0.261974 0.77714
asc_b_leisure 0.106287 0.148497 0.7157 0.255373 0.41620
asc_b_female 0.000000 NA NA NA NA
asc_b_male 0.175134 0.055085 3.1793 0.090286 1.93976
asc_b_divers 0.146116 0.439456 0.3325 1.130482 0.12925
asc_b_age 0.058708 0.010796 5.4381 0.017788 3.30043
asc_b_age_2 -5.9628e-04 1.1541e-04 -5.1666 1.9361e-04 -3.07982
asc_w_base -2.587003 0.654354 -3.9535 1.193279 -2.16798
asc_w_work 0.000000 NA NA NA NA
asc_w_errand 0.299264 0.412286 0.7259 0.854242 0.35033
asc_w_leisure 0.072829 0.421915 0.1726 0.857031 0.08498
asc_w_female 0.000000 NA NA NA NA
asc_w_male -0.177903 0.103643 -1.7165 0.184262 -0.96549
asc_w_divers 2.212046 0.449991 4.9158 0.711920 3.10715
asc_w_age 0.052144 0.020034 2.6027 0.032396 1.60957
asc_w_age_2 -3.2780e-04 2.0570e-04 -1.5936 3.4243e-04 -0.95727
s_mocono 1.000000 NA NA NA NA
s_mocosi 0.972736 0.029385 33.1034 0.032661 29.78236
b_cost_base -0.150005 0.011513 -13.0297 0.018597 -8.06591
elast_cost_dist -0.717950 0.053069 -13.5287 0.107998 -6.64782
b_mcost_rembudget75days15 -0.036301 0.009344 -3.8850 0.010686 -3.39720
b_mcost_rembudget75days20 -0.013583 0.008767 -1.5493 0.009071 -1.49732
b_mcost_rembudget75days25 -0.008614 0.009316 -0.9246 0.011301 -0.76224
b_mcost_rembudget50days15 -0.080285 0.014135 -5.6800 0.017842 -4.49975
b_mcost_rembudget50days20 -0.011427 0.012140 -0.9413 0.012592 -0.90748
b_mcost_rembudget50days25 -0.013329 0.012619 -1.0563 0.014097 -0.94551
b_mcost_rembudget25days15 -0.080614 0.011172 -7.2158 0.013679 -5.89342
b_mcost_rembudget25days20 -0.054174 0.010026 -5.4033 0.012696 -4.26700
b_mcost_rembudget25days25 -0.058214 0.011242 -5.1783 0.014191 -4.10215
elast_mcost_dist -0.782566 0.086440 -9.0533 0.101831 -7.68494
elast_mcost_inc -0.255097 0.094963 -2.6863 0.116134 -2.19657
v_tt_pt_base 0.400504 0.050325 7.9584 0.078902 5.07597
v_tt_c_base 0.289455 0.066391 4.3599 0.102945 2.81174
v_tt_b_base 0.407471 0.047868 8.5124 0.076805 5.30525
v_tt_w_base 0.535308 0.079314 6.7492 0.156265 3.42565
elast_vtt_dist 0.546130 0.069166 7.8960 0.134232 4.06855
mult_vtt_pt_leisure 0.910207 0.126314 7.2059 0.183174 4.96910
mult_vtt_pt_shopping 0.627084 0.107013 5.8599 0.159394 3.93419
mult_vtt_c_leisure 0.429602 0.188720 2.2764 0.310528 1.38346
mult_vtt_c_shopping 0.458857 0.191549 2.3955 0.292903 1.56659
mult_vtt_w_leisure 0.740412 0.111332 6.6505 0.211563 3.49972
mult_vtt_w_shopping 0.757512 0.114156 6.6358 0.235887 3.21134
mult_vtt_b_leisure 0.865357 0.087197 9.9241 0.129897 6.66187
mult_vtt_b_shopping 0.693743 0.079515 8.7247 0.125469 5.52918
v_aet_pt 0.116845 0.031697 3.6863 0.064021 1.82512
v_freq_pt 0.061717 0.030121 2.0490 0.046273 1.33377
v_trans_pt 0.705802 0.156862 4.4995 0.205693 3.43134
b_mgain_b_base 0.208319 0.086838 2.3989 0.099356 2.09670
b_qual_b_bad 0.000000 NA NA NA NA
b_qual_b_med 0.762972 0.061527 12.4007 0.069484 10.98061
b_qual_b_vgood 0.942095 0.062026 15.1886 0.067958 13.86289
b_weather_c_rainy 0.000000 NA NA NA NA
b_weather_c_sunny -0.373262 0.068852 -5.4212 0.083909 -4.44843
b_weather_b_rainy 0.000000 NA NA NA NA
b_weather_b_sunny 2.830420 0.073280 38.6249 0.096785 29.24438
b_weather_w_rainy 0.000000 NA NA NA NA
b_weather_w_sunny 2.104389 0.111653 18.8475 0.165460 12.71845

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 cnl_settings, order is irrelevant
V = list()

# ascs with shifts
asc_pt = asc_pt_base
asc_c = asc_c_base +
asc_c_work * (purpose == 1) + asc_c_errand * (purpose == 2) + asc_c_leisure * (purpose == 3) +
asc_c_female * (p_gender == "weiblich") + asc_c_male * (p_gender == "männlich") + asc_c_divers * (p_gender == "divers") +
asc_c_age * p_age + asc_c_age_2 * p_age**2
asc_b = asc_b_base +
asc_b_work * (purpose == 1) + asc_b_errand * (purpose == 2) + asc_b_leisure * (purpose == 3) +
asc_b_female * (p_gender == "weiblich") + asc_b_male * (p_gender == "männlich") + asc_b_divers * (p_gender == "divers") +
asc_b_age * p_age + asc_b_age_2 * p_age**2
asc_w = asc_w_base +
asc_w_work * (purpose == 1) + asc_w_errand * (purpose == 2) + asc_w_leisure * (purpose == 3) +
asc_w_female * (p_gender == "weiblich") + asc_w_male * (p_gender == "männlich") + asc_w_divers * (p_gender == "divers") +
asc_w_age * p_age + asc_w_age_2 * p_age**2

# b_cost_base = cost sensitivity at mean distance, elast_cost_dist < 0, decreasing cost elasticity in relation to distance
b_cost = (b_cost_base) *
((trip_dist/mean(trip_dist)) ^ elast_cost_dist)

#b_mcost = b_mcost_base

b_mcost = (b_mcost_rembudget75days15 * (cbudget == 75) * (days_into_month == 15) +
b_mcost_rembudget75days20 * (cbudget == 75) * (days_into_month == 20) +
b_mcost_rembudget75days25 * (cbudget == 75) * (days_into_month == 25) +
b_mcost_rembudget50days15 * (cbudget == 50) * (days_into_month == 15) +
b_mcost_rembudget50days20 * (cbudget == 50) * (days_into_month == 20) +
b_mcost_rembudget50days25 * (cbudget == 50) * (days_into_month == 25) +
b_mcost_rembudget25days15 * (cbudget == 25) * (days_into_month == 15) +
b_mcost_rembudget25days20 * (cbudget == 25) * (days_into_month == 20) +
b_mcost_rembudget25days25 * (cbudget == 25) * (days_into_month == 25)) *
((trip_dist/mean(trip_dist)) ^ elast_mcost_dist) *
((p_inc_lin/mean(p_inc_lin)) ^ elast_mcost_inc)

b_mgain_b = (b_mgain_b_base)

v_tt_pt = v_tt_pt_base *
((trip_dist/mean(trip_dist)) ^ elast_vtt_dist) *
(1*(purpose == 1) + mult_vtt_pt_shopping*(purpose == 2) + mult_vtt_pt_leisure * (purpose == 3))

v_tt_c = v_tt_c_base *
((trip_dist/mean(trip_dist)) ^ elast_vtt_dist) *
(1*(purpose == 1) + mult_vtt_c_shopping*(purpose == 2) + mult_vtt_c_leisure * (purpose == 3))

v_tt_b = v_tt_b_base *
((trip_dist/mean(trip_dist)) ^ elast_vtt_dist) *
(1*(purpose == 1) + mult_vtt_b_shopping*(purpose == 2) + mult_vtt_b_leisure * (purpose == 3))

v_tt_w = v_tt_w_base *
((trip_dist/mean(trip_dist)) ^ elast_vtt_dist) *
(1*(purpose == 1) + mult_vtt_w_shopping*(purpose == 2) + mult_vtt_w_leisure * (purpose == 3))

# mocono
# reference alternative
V[["pt_mocono"]] = asc_pt +
b_cost * (cost_pt + v_tt_pt * tt_pt + v_aet_pt * aett_pt + v_freq_pt * freq_pt + v_trans_pt * trans_pt)
# other alternatives
V[["c_mocono"]] = asc_c +
b_cost * (cost_c + v_tt_c * tt_c) +
b_weather_c_rainy * (weather == "rainy") + b_weather_c_sunny * (weather == "sunny")
V[["b_mocono"]] = asc_b +
b_cost * (v_tt_b * tt_b) +
b_qual_b_bad * (quality_b == "bad") + b_qual_b_med * (quality_b == "medium") + b_qual_b_vgood * (quality_b == "very good") +
b_weather_b_rainy * (weather == "rainy") + b_weather_b_sunny * (weather == "sunny")
V[["w_mocono"]] = asc_w +
b_cost * (v_tt_w * tt_w) +
b_weather_w_rainy * (weather == "rainy") + b_weather_w_sunny * (weather == "sunny")

# mocosi
# reference alternative
V[["pt_mocosi"]] = asc_pt +
b_cost * (cost_pt + v_tt_pt * tt_pt + v_aet_pt * aett_pt + v_freq_pt * freq_pt + v_trans_pt * trans_pt) +
b_mcost * ccost_pt
# other alternatives
V[["c_mocosi"]] = asc_c +
b_cost * (cost_c + v_tt_c * tt_c) +
b_mcost * ccost_c +
b_weather_c_rainy * (weather == "rainy") + b_weather_c_sunny * (weather == "sunny")
V[["b_mocosi"]] = asc_b +
b_cost * (v_tt_b * tt_b) +
b_mgain_b * cgain_b +
b_qual_b_bad * (quality_b == "bad") + b_qual_b_med * (quality_b == "medium") + b_qual_b_vgood * (quality_b == "very good") +
b_weather_b_rainy * (weather == "rainy") + b_weather_b_sunny * (weather == "sunny")
V[["w_mocosi"]] = asc_w +
b_cost * (v_tt_w * tt_w) +
b_weather_w_rainy * (weather == "rainy") + b_weather_w_sunny * (weather == "sunny")

# Compute probabilities for the mocono mode choice part of the data using MNL model
mnl_settings_mocono = list(
alternatives = c(pt_mocono = 1, c_mocono = 2, b_mocono = 3, w_mocono = 4),
avail = list(pt_mocono = pt_availability, c_mocono = car_availability, b_mocono = bicycle_availability, w_mocono = walk_availability),
choiceVar = CHOICE,
utilities = list(pt_mocono = s_mocono * V[["pt_mocono"]],
c_mocono = s_mocono * V[["c_mocono"]],
b_mocono = s_mocono * V[["b_mocono"]],
w_mocono = s_mocono * V[["w_mocono"]]),
rows = (exp_nr == 0)
)

P[["mocono"]] = apollo_mnl(mnl_settings_mocono, functionality)

# Compute probabilities for the RP mode choice part of the data using MNL model
mnl_settings_mocosi = list(
alternatives = c(pt_mocosi = 5, c_mocosi = 6, b_mocosi = 7, w_mocosi = 8),
avail = list(pt_mocosi = pt_availability, c_mocosi = car_availability, b_mocosi = bicycle_availability, w_mocosi = walk_availability),
choiceVar = CHOICE,
utilities = list(pt_mocosi = s_mocosi * V[["pt_mocosi"]],
c_mocosi = s_mocosi * V[["c_mocosi"]],
b_mocosi = s_mocosi * V[["b_mocosi"]],
w_mocosi = s_mocosi * V[["w_mocosi"]]),
rows = (exp_nr == 1)
)

P[["mocosi"]] = apollo_mnl(mnl_settings_mocosi, functionality)

# Combined model
P = apollo_combineModels(P, apollo_inputs, functionality, components = c("mocono","mocosi"))

# 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)
}
All good so far. Now I would like to estimate an MMNL with random parameters for 4 different travel times variables (walk, bike, car and PT) as well as one for the cost parameter. see here:

Code: Select all

### Set parameters for generating draws
apollo_draws = list(
  interDrawsType = "sobol",
  interNDraws    = 1000,
  interUnifDraws = c(),
  interNormDraws = c("draws_tt_pt", "draws_tt_c", "draws_tt_b","draws_tt_w","draws_tc")
)

### Create random parameters
apollo_randCoeff = function(apollo_beta, apollo_inputs){
  randcoeff = list()
  
  # Log normal
  randcoeff[["v_tt_c"]]  = exp( mu_log_v_tt_c_base  + sigma_log_v_tt_c_base  * draws_tt_c  ) *
    ((trip_dist/mean(trip_dist)) ^ elast_vtt_dist)  *
    (1*(purpose == 1) + mult_vtt_c_shopping*(purpose == 2) + mult_vtt_c_leisure * (purpose == 3))
  
  randcoeff[["v_tt_w"]]  = exp( mu_log_v_tt_w_base  + sigma_log_v_tt_w_base  * draws_tt_w  ) *
    ((trip_dist/mean(trip_dist)) ^ elast_vtt_dist)  *
    (1*(purpose == 1) + mult_vtt_w_shopping*(purpose == 2) + mult_vtt_w_leisure * (purpose == 3))
  
  randcoeff[["v_tt_b"]]  = exp( mu_log_v_tt_b_base  + sigma_log_v_tt_b_base  * draws_tt_b  ) *
    ((trip_dist/mean(trip_dist)) ^ elast_vtt_dist)  *
    (1*(purpose == 1) + mult_vtt_b_shopping*(purpose == 2) + mult_vtt_b_leisure * (purpose == 3))
  
  randcoeff[["v_tt_pt"]] = exp( mu_log_v_tt_pt_base + sigma_log_v_tt_pt_base * draws_tt_pt ) *
    ((trip_dist/mean(trip_dist)) ^ elast_vtt_dist)  *
    (1*(purpose == 1) + mult_vtt_pt_shopping*(purpose == 2) + mult_vtt_pt_leisure * (purpose == 3))
  
  # Negative log normal
  randcoeff[["b_cost"]]  = -exp( mu_b_cost_base +  sigma_b_cost_base  * draws_tc ) *
    ((trip_dist/mean(trip_dist)) ^ elast_cost_dist)
  
  return(randcoeff)
}
For the vtt's we assume a lognormal whereas for the cost coefficient we assume a negative lognormal. These are the starting values:

Code: Select all

mu_b_cost_base = -1,
sigma_b_cost_base = 0,
elast_cost_dist = 0,
mu_log_v_tt_pt_base = 2,
  sigma_log_v_tt_pt_base = 0,
  mu_log_v_tt_c_base = 2,
  sigma_log_v_tt_c_base = 0,
  mu_log_v_tt_b_base = 2,
  sigma_log_v_tt_b_base = 0,
  mu_log_v_tt_w_base = 2,
  sigma_log_v_tt_w_base = 0,
  elast_vtt_dist = 0,
  mult_vtt_pt_leisure = 1,
  mult_vtt_pt_shopping = 1,
  mult_vtt_c_leisure = 1,
  mult_vtt_c_shopping = 1,
  mult_vtt_w_leisure = 1,
  mult_vtt_w_shopping = 1,
  mult_vtt_b_leisure = 1,
  mult_vtt_b_shopping = 1,
For the other parameters, I read in the MNL's estimates as starting values. However, I run into troubles =/.

Would appreciate any help here... I assume it might be related to the starting values... but I am stuck at the moment.

Thomas

Re: CALCULATION ISSUE - Log-likelihood calculation fails at values close to the starting values!

Posted: 20 Jul 2023, 15:44
by stephanehess
Hi

the issue is your starting values. With lognormals, you should not start at positive values for the mean of the log of the coefficient - see http://apollochoicemodelling.com/faq.html

Stephane