I have tried to apply the MMNL (using the codes in website: MMNL_preference_space) with my data.
There are five attributes: barn, range, oragnic, vitamin, price. The first four are dummy variables, the last one is continuous variable.
I set random heterogeneity for the first four attributes as follows
Code: Select all
randcoeff[["b_barn"]] = -exp(mu_log_b_barn + sigma_log_b_barn * draws_barn)
randcoeff[["b_range"]] = -exp(mu_log_b_range + sigma_log_b_range * draws_range)
randcoeff[["b_organic"]] = -exp(mu_log_b_organic + sigma_log_b_organic * draws_organic)
randcoeff[["b_vitamin"]] = -exp(mu_log_b_vitamin + sigma_log_b_vitamin * draws_vitamin)
Code: Select all
### Clear memory
rm(list = ls())
### Load libraries
library(openxlsx)
library(apollo)
### Initialise code
apollo_initialise()
### Set core controls
apollo_control = list(
modelName = "MMNL_preference_space",
modelDescr = "Mixed logit model choice data, uncorrelated Lognormals in preference space",
indivID = "ID",
nCores = 4,
outputDirectory = "output"
)
# ####################################################### #
#### 2. Data loading ####
# ####################################################### #
dtaApollo <- read.xlsx("dtaApollo.xlsx", sheet = "data", colNames = T)
database <- dtaApollo
# ####################################################### #
#### 3. Parameter definition ####
# ####################################################### #
### Vector of parameters, including any that are kept fixed during estimation
apollo_beta=c(asc_alt3 = 0,
mu_log_b_barn = 0,
sigma_log_b_barn = 0,
mu_log_b_range = 0,
sigma_log_b_range = 0,
mu_log_b_organic = 0,
sigma_log_b_organic = 0,
mu_log_b_vitamin = 0,
sigma_log_b_vitamin = 0,
# b_barn = 0,
# b_range = 0,
b_price = 0)
### Vector with names (in quotes) of parameters to be kept fixed at their starting value in apollo_beta.
### Use apollo_beta_fixed = c() for no fixed parameters.
apollo_fixed = c()
### Set parameters for generating draws
apollo_draws = list(
interDrawsType = "halton",
interNDraws = 500,
interUnifDraws = c(),
interNormDraws = c("draws_barn","draws_range","draws_organic","draws_vitamin"),
intraDrawsType = "halton",
intraNDraws = 0,
intraUnifDraws = c(),
intraNormDraws = c()
)
### Create random parameters
apollo_randCoeff = function(apollo_beta, apollo_inputs){
randcoeff = list()
randcoeff[["b_barn"]] = -exp(mu_log_b_barn + sigma_log_b_barn * draws_barn)
randcoeff[["b_range"]] = -exp(mu_log_b_range + sigma_log_b_range * draws_range)
randcoeff[["b_organic"]] = -exp(mu_log_b_organic + sigma_log_b_organic * draws_organic)
randcoeff[["b_vitamin"]] = -exp(mu_log_b_vitamin + sigma_log_b_vitamin * draws_vitamin)
return(randcoeff)
}
# ####################################################### #
#### 4. Input validation ####
# ####################################################### #
apollo_inputs = apollo_validateInputs()
# ####################################################### #
#### 5. Likelihood definition ####
# ####################################################### #
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"]] = b_barn * (barn_alt1 == 1) + b_range * (range_alt1 == 1) + b_organic * (organic_alt1 == 1) + b_vitamin * (vitamin_alt1 == 1) + b_price * price_alt1
V[["alt2"]] = b_barn * (barn_alt2 == 1) + b_range * (range_alt2 == 1) + b_organic * (organic_alt2 == 1) + b_vitamin * (vitamin_alt2 == 1) + b_price * price_alt2
V[["alt3"]] = asc_alt3 + b_barn * (barn_alt3 == 1) + b_range * (range_alt3 == 1) + b_organic * (organic_alt3 == 1) + b_vitamin * (vitamin_alt3 == 1) + b_price * price_alt3
### Define settings for MNL model component
mnl_settings = list(
alternatives = c(alt1=1, alt2=2, alt3=3),
avail = list(alt1=av_alt1, alt2=av_alt2, alt3=av_alt3),
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)
### 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)
}
# ####################################################### #
#### 6. Model estimation and reporting ####
# ####################################################### #
model = apollo_estimate(apollo_beta, apollo_fixed,
apollo_probabilities,
apollo_inputs)
apollo_modelOutput(model) #to screen
apollo_saveOutput(model) #to file
WARNING: Singular Hessian, cannot calculate s.e.
Could not write hessian to a file.
WARNING: Some eigenvalues of the Hessian are positive, indicating convergence to a saddle point!
For the results, it seems not good as mu_log is very high compared to b_price.
I guess, there is something wrong in my model definition for the mixed logit. Do you have any idea what that might be?Preparing user-defined functions.
Testing likelihood function...
Overview of choices for MNL model component :
alt1 alt2 alt3
Times available 2960.00 2960.00 2960.00
Times chosen 955.00 1307.00 698.00
Percentage chosen overall 32.26 44.16 23.58
Percentage chosen when available 32.26 44.16 23.58
Pre-processing likelihood function...
Creating cluster...
Preparing workers for multithreading...
Testing influence of parameters
Starting main estimation
Initial function value: -4552.701
Initial gradient value:
asc_alt3 mu_log_b_barn sigma_log_b_barn mu_log_b_range sigma_log_b_range
-1.065371e+03 -3.578956e+02 4.327060e-02 -4.779576e+02 -1.088076e-01
mu_log_b_organic sigma_log_b_organic mu_log_b_vitamin sigma_log_b_vitamin b_price
-4.221774e+02 2.523584e-02 -6.998596e+02 -3.929654e-01 1.276747e+03
initial value 4552.700594
iter 2 value 4284.763351
iter 3 value 3064.798184
iter 4 value 3043.632602
iter 5 value 3031.589350
iter 6 value 2855.622805
iter 7 value 2821.452787
iter 8 value 2818.394461
iter 9 value 2809.551204
iter 10 value 2798.110528
iter 11 value 2761.134041
iter 12 value 2751.150277
iter 13 value 2738.469340
iter 14 value 2732.273568
iter 15 value 2726.800654
iter 16 value 2714.178608
iter 17 value 2703.314911
iter 18 value 2691.460221
iter 19 value 2686.312072
iter 20 value 2682.469001
iter 21 value 2679.635570
iter 22 value 2677.734367
iter 23 value 2677.731758
iter 24 value 2677.713236
iter 25 value 2677.534510
iter 26 value 2677.352348
iter 27 value 2676.963312
iter 28 value 2676.203104
iter 29 value 2675.747670
iter 30 value 2675.392379
iter 31 value 2675.183195
iter 32 value 2675.052192
iter 33 value 2674.990024
iter 34 value 2674.967226
iter 35 value 2674.912472
iter 36 value 2674.869342
iter 37 value 2674.840648
iter 38 value 2674.818647
iter 39 value 2674.805682
iter 40 value 2674.798233
iter 41 value 2674.792189
iter 42 value 2674.783433
iter 43 value 2674.764274
iter 44 value 2674.764156
iter 45 value 2674.764008
iter 46 value 2674.760705
iter 47 value 2674.757750
iter 48 value 2674.755915
iter 49 value 2674.754337
iter 50 value 2674.753013
iter 51 value 2674.751115
iter 52 value 2674.750099
iter 53 value 2674.748753
iter 54 value 2674.738955
iter 55 value 2674.711932
iter 56 value 2674.671580
iter 57 value 2674.629577
iter 58 value 2674.281864
iter 59 value 2673.938510
iter 60 value 2673.894603
iter 61 value 2673.859675
iter 62 value 2673.850345
iter 63 value 2673.830013
iter 64 value 2673.819293
iter 65 value 2673.818798
iter 66 value 2673.817939
iter 67 value 2673.813243
iter 68 value 2673.804321
iter 69 value 2673.791031
iter 70 value 2673.775546
iter 70 value 2673.775531
final value 2673.775531
converged
Additional convergence test using scaled estimation. Parameters will be scaled by their
current estimates and additional iterations will be performed.
initial value 2673.775531
iter 2 value 2673.733456
iter 3 value 2673.717491
iter 4 value 2673.712277
iter 5 value 2673.708861
iter 6 value 2673.636402
iter 7 value 2673.610078
iter 8 value 2673.513912
iter 9 value 2673.447813
iter 10 value 2673.398029
iter 11 value 2673.384613
iter 12 value 2673.381160
iter 13 value 2673.377515
iter 13 value 2673.377515
final value 2673.377515
converged
Estimated parameters:
Estimate
asc_alt3 -2.4205
mu_log_b_barn -184.3784
sigma_log_b_barn 205.8515
mu_log_b_range -22.4080
sigma_log_b_range -1.9684
mu_log_b_organic -21.1239
sigma_log_b_organic -2.6224
mu_log_b_vitamin -7.5396
sigma_log_b_vitamin -6.6235
b_price -0.7075
Final LL: -2673.3775
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% (20 NA values)
Computing covariance matrix using numerical methods (numDeriv).
0%....25%....50%....75%....100% (36 NA values)
Computing covariance matrix using numerical methods (maxLik). This may take a while, no
progress bar displayed.
Thanks very much.
--
Cuong