eMDC model with 28 alternatives (Burger ingredients)

[malemu588@gmail.com]

Dear Stephane and David,

I have been working with the eMDC model since May 2025, and I have been able to specify and estimate a model. I read the published paper on eMDC and followed all the recommendations, including those available in this forum. For example, I reduced the number of parameters to be estimated to the lowest possible number.

My concerns, where I need your help, are the following:

1. My model is very sensitive to the starting value of sigma. When I increase the starting value—for example, to 3—the model runs well and a maximum can be found.
2. My model takes a very long time to run; the overall estimation process takes more than a day. Is this something I should be worried about?
3. After obtaining the results, I conducted some predictions, and for some food products that are substitutes (e.g., conventional and organic beef), the cross-price elasticities turn out negative, which does not make sense. Do you think this might relate to the complexity of the eMDC model?

I am sorry for the many questions, but I wish to learn more from this forum in addition to the existing materials.

My R code is attached below.

Kind regards,
Mohammed

### Clear memory
rm(list = ls())
###update apollo
###install.packages("apollo")

### Load Apollo library
library(apollo)

### Set working directory (only works in RStudio)
setwd("C:/Users/zsc245/OneDrive - University of Copenhagen/NOVASOIL/Data analysis/Burger")
getwd()
##apollo_setWorkDir()

### Initialise code
apollo_initialise()

### Set core controls
apollo_control = list(
modelDescr = "Burger",
indivID = "respondent_id",
panelData = FALSE,
workInLogs = TRUE,
outputDirectory = "eMDEOutput")

# ################################################################# #
#### LOAD DATA AND APPLY ANY TRANSFORMATIONS ####
# ################################################################# #

### Load data from within the Apollo package
database = read.csv("burger_analysis_panel.csv",header=TRUE, na.strings=".")
database <-subset(database,!(database$patties==0&buns>0))
database <-subset(database,!(database$patties>0&buns==0))

### Create consumption variables for combined activities
# outside good: time spent at home and travelling

### Randomly split dataset into estimation (70%) and validation (30%)
## set.seed(1)
## database$validation = runif(nrow(database))>0.7
## dbVal = database[ database$validation,] # validation sample
## database = database[!database$validation,] # estimation sample


# ################################################################# #
#### DEFINE MODEL PARAMETERS ####
# ################################################################# #

### Parameters starting values c(name1=value1, name2=value2, ...)
apollo_beta = c(#afemale=0,
sigma = 3,
# Base utility of inside good
b1 = 0,
b2 = 0,
b3 = 0,
b4 = 0,
b5 = 0,
b6 = 0,
b7 = 0,
b8 = 0,
b9 = 0,
b10 = 0,
b11 = 0,
b12 = 0,
b13 = 0,
b14 = 0,
b15 = 0,
b16 = 0,
b17 = 0,
b18 = 0,
b19 = 0,
b20 = 0,
b21 = 0,
b22 = 0,
b23 = 0,
b24 = 0,
b25 = 0,
b26 = 0,
b27 = 0,
b28 = 0,
# Satiation parameters
g1 = 1,
g2 = 1,
g3 = 1,
g4 = 1,
g5 = 1,
g6 = 1,
g7 = 1,
#g8 = 0,
#g9 = 0,
#g10 = 0,
#g11 = 0,
#g12 = 0,
#g13 = 0,
#g14 = 0,
#g15 = 0,
#g16 = 0,
#g17 = 0,
#g18 = 0,
#g19 = 0,
#g20 = 0,
#g21 = 0,
#g22 = 0,
#g23 = 0,
#g24 = 0,
#g25 = 0,
#g26 = 0,
#g27 = 0,
#g28 = 0,
# Complementary/substitution effects
d1_2 = 0,
d1_3 = 0,
d1_4 = 0,
d1_5 = 0,
d1_9 = 0,
d1_13 = 0,
d1_17 = 0,
d1_21 = 0,
d1_25 = 0,
d2_3 = 0,
d2_4 = 0,
d2_5 = 0,
d2_9 = 0,
d2_13 = 0,
d2_17 = 0,
d2_21 = 0,
d2_25 = 0,
d3_4 = 0,
d3_5 = 0,
d3_9 = 0,
d3_13 = 0,
d3_17 = 0,
d3_21 = 0,
d3_25 = 0,
d4_5 = 0,
d4_9 = 0,
d4_13 = 0,
d4_17 = 0,
d4_21 = 0,
d4_25 = 0,
#d5_6 = 0,
#d5_7 = 0,
#d5_8 = 0,
d5_9 = 0,
d5_13 = 0,
d5_17 = 0,
d5_21 = 0,
d5_25 = 0,
#d6_7 = 0,
#d6_8 = 0,
d6_9 = 0,
d6_13 = 0,
d6_17 = 0,
d6_21 = 0,
d6_25 = 0,
#d7_8 = 0,
d7_9 = 0,
d7_13 = 0,
d7_17 = 0,
d7_21 = 0,
d7_25 = 0,
d8_9 = 0,
d8_13 = 0,
d8_17 = 0,
d8_21 = 0,
d8_25 = 0,
#d9_10 = 0,
#d9_11 = 0,
#d9_12 = 0,
d9_13 = 0,
d9_17 = 0,
d9_21 = 0,
d9_25 = 0,
#d10_11 = 0,
#d10_12 = 0,
d10_13 = 0,
d10_17 = 0,
d10_21 = 0,
d10_25 = 0,
#d11_12 = 0,
d11_13 = 0,
d11_17 = 0,
d11_21 = 0,
d11_25 = 0,
d12_13 = 0,
d12_17 = 0,
d12_21 = 0,
d12_25 = 0,
d13_14 = 0,
d13_15 = 0,
d13_16 = 0,
d13_17 = 0,
d13_21 = 0,
d13_25 = 0,
d14_15 = 0,
d14_16 = 0,
d14_17 = 0,
d14_21 = 0,
d14_25 = 0,
d15_16 = 0,
d15_17 = 0,
d15_21 = 0,
d15_25 = 0,
d16_17 = 0,
d16_21 = 0,
d16_25 = 0,
d17_18 = 0,
d17_19 = 0,
d17_20 = 0,
d17_21 = 0,
d17_25 = 0,
d18_19 = 0,
d18_20 = 0,
d18_21 = 0,
d18_25 = 0,
d19_20 = 0,
d19_21 = 0,
d19_25 = 0,
d20_21 = 0,
d20_25 = 0,
d21_22 = 0,
d21_23 = 0,
d21_24 = 0,
d21_25 = 0,
d22_23 = 0,
d22_24 = 0,
d22_25 = 0,
d23_24 = 0,
d23_25 = 0,
d24_25 = 0,
d25_26 = 0,
d25_27 = 0,
d25_28 = 0,
d26_27 = 0,
d26_28 = 0,
d27_28 = 0)
### Names of fixed parameters
apollo_fixed = c('g1', 'g2', 'g3', 'g4', 'g5', 'g6', 'g7', 'd1_5', 'd1_9', 'd1_13', 'd1_17', 'd1_21', 'd1_25', 'd2_5', 'd2_9', 'd2_13',
'd2_17', 'd2_21', 'd2_25', 'd3_5', 'd3_9', 'd3_13', 'd3_17', 'd3_21', 'd3_25', 'd4_5', 'd4_9', 'd4_13', 'd4_17', 'd4_21',
'd4_25', 'd5_9', 'd5_13', 'd5_17', 'd5_21', 'd5_25', 'd6_9', 'd6_13', 'd6_17', 'd6_21', 'd6_25', 'd7_9', 'd7_13', 'd7_17', 'd7_21',
'd7_25', 'd8_9', 'd8_13', 'd8_17', 'd8_21', 'd8_25', 'd9_13', 'd9_17', 'd9_21', 'd9_25', 'd10_13', 'd10_17', 'd10_21', 'd10_25',
'd11_13', 'd11_17', 'd11_21', 'd11_25', 'd12_13', 'd12_17', 'd12_21', 'd12_25', 'd13_17', 'd13_21', 'd13_25', 'd14_17', 'd14_21',
'd14_25', 'd15_17', 'd15_21', 'd15_25', 'd16_17', 'd16_21', 'd16_25', 'd17_21', 'd17_25', 'd18_21', 'd18_25', 'd19_21', 'd19_25',
'd20_21', 'd20_25', 'd21_25', 'd22_25', 'd23_25', 'd24_25')
# ################################################################# #
#### GROUP AND VALIDATE INPUTS ####
# ################################################################# #

apollo_inputs = apollo_validateInputs()

# ################################################################# #
#### DEFINE MODEL AND LIKELIHOOD FUNCTION ####
# ################################################################# #

apollo_probabilities=function(apollo_beta, apollo_inputs,
functionality="estimate"){

### Initialise
apollo_attach(apollo_beta, apollo_inputs)
on.exit(apollo_detach(apollo_beta, apollo_inputs))
P = list()

### Prepare Inputs
alts = c("beefconv", "beefsoilh", "beeforgan", "beefboth", "chickenconv","chickensoilh", "chickenorgan", "chickenboth", "bollerconv", "bollersoilh", "bollerorgan", "bollerboth", "plantconv", "plantsoilh", "plantorgan","plantboth", "salatconv", "salatsoilh", "salatorgan", "salatboth","tomatconv","tomatsoilh","tomatorgan","tomatboth","agurkconv","agurksoilh","agurkorgan","agurkboth" )
nAlt = length(alts)
ones = setNames(as.list(rep(1, nAlt)), alts)
continuousChoice = list(beefconv=quantityItem1,
beefsoilh=quantityItem2,
beeforgan=quantityItem3,
beefboth=quantityItem4,
chickenconv=quantityItem5,
chickensoilh=quantityItem6,
chickenorgan=quantityItem7,
chickenboth=quantityItem8,
plantconv=quantityItem9,
plantsoilh=quantityItem10,
plantorgan=quantityItem11,
plantboth=quantityItem12,
bollerconv=quantityItem13,
bollersoilh=quantityItem14,
bollerorgan=quantityItem15,
bollerboth=quantityItem16,
salatconv=quantityItem17,
salatsoilh=quantityItem18,
salatorgan=quantityItem19,
salatboth=quantityItem20,
tomatconv=quantityItem21,
tomatsoilh=quantityItem22,
tomatorgan=quantityItem23,
tomatboth=quantityItem24,
agurkconv=quantityItem25,
agurksoilh=quantityItem26,
agurkorgan=quantityItem27,
agurkboth=quantityItem28)

cost = list(beefconv=priceItem1,
beefsoilh=priceItem2,
beeforgan=priceItem3,
beefboth=priceItem4,
chickenconv=priceItem5,
chickensoilh=priceItem6,
chickenorgan=priceItem7,
chickenboth=priceItem8,
plantconv=priceItem9,
plantsoilh=priceItem10,
plantorgan=priceItem11,
plantboth=priceItem12,
bollerconv=priceItem13,
bollersoilh=priceItem14,
bollerorgan=priceItem15,
bollerboth=priceItem16,
salatconv=priceItem17,
salatsoilh=priceItem18,
salatorgan=priceItem19,
salatboth=priceItem20,
tomatconv=priceItem21,
tomatsoilh=priceItem22,
tomatorgan=priceItem23,
tomatboth=priceItem24,
agurkconv=priceItem25,
agurksoilh=priceItem26,
agurkorgan=priceItem27,
agurkboth=priceItem28)

utilities = list(
beefconv = b1,
beefsoilh = b2,
beeforgan = b3,
beefboth = b4,
chickenconv = b5,
chickensoilh = b6,
chickenorgan = b7,
chickenboth = b8,
plantconv = b9,
plantsoilh = b10,
plantorgan = b11,
plantboth = b12,
bollerconv = b13,
bollersoilh = b14,
bollerorgan = b15,
bollerboth = b16,
salatconv = b17,
salatsoilh = b18,
salatorgan = b19,
salatboth = b20,
tomatconv = b21,
tomatsoilh = b22,
tomatorgan = b23,
tomatboth = b24,
agurkconv = b25,
agurksoilh = b26,
agurkorgan = b27,
agurkboth = b28
)
gamma = list( beefconv = g1,
beefsoilh = g1,
beeforgan = g1,
beefboth = g1,
chickenconv = g2,
chickensoilh = g2,
chickenorgan = g2,
chickenboth = g2,
plantconv = g3,
plantsoilh = g3,
plantorgan = g3,
plantboth = g3,
bollerconv = g4,
bollersoilh = g4,
bollerorgan = g4,
bollerboth = g4,
salatconv = g5,
salatsoilh = g5,
salatorgan = g5,
salatboth = g5,
tomatconv = g6,
tomatsoilh = g6,
tomatorgan = g6,
tomatboth = g6,
agurkconv = g7,
agurksoilh = g7,
agurkorgan = g7,
agurkboth = g7

)
delta = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
d1_2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
d1_3, d2_3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
d1_4, d2_4, d3_4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
d1_5, d2_5, d3_5, d4_5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
d1_5, d2_5, d3_5, d4_5, d1_2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
d1_5, d2_5, d3_5, d4_5, d1_3, d2_3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
d1_5, d2_5, d3_5, d4_5, d1_4, d2_4, d3_4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
d1_9, d2_9, d3_9, d4_9, d5_9, d6_9, d7_9, d8_9, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
d1_9, d2_9, d3_9, d4_9, d5_9, d6_9, d7_9, d8_9, d1_2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
d1_9, d2_9, d3_9, d4_9, d5_9, d6_9, d7_9, d8_9, d1_3, d2_3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
d1_9, d2_9, d3_9, d4_9, d5_9, d6_9, d7_9, d8_9, d1_4, d2_4, d3_4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
d1_13, d2_13, d3_13, d4_13, d5_13, d6_13, d7_13, d8_13, d9_13, d10_13, d11_13, d12_13, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
d1_13, d2_13, d3_13, d4_13, d5_13, d6_13, d7_13, d8_13, d9_13, d10_13, d11_13, d12_13, d13_14, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
d1_13, d2_13, d3_13, d4_13, d5_13, d6_13, d7_13, d8_13, d9_13, d10_13, d11_13, d12_13, d13_15, d14_15, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
d1_13, d2_13, d3_13, d4_13, d5_13, d6_13, d7_13, d8_13, d9_13, d10_13, d11_13, d12_13, d13_16, d14_16, d15_16, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
d1_17, d2_17, d3_17, d4_17, d5_17, d6_17, d7_17, d8_17, d9_17, d10_17, d11_17, d12_17, d13_17, d14_17, d15_17, d16_17, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
d1_17, d2_17, d3_17, d4_17, d5_17, d6_17, d7_17, d8_17, d9_17, d10_17, d11_17, d12_17, d13_17, d14_17, d15_17, d16_17, d17_18, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
d1_17, d2_17, d3_17, d4_17, d5_17, d6_17, d7_17, d8_17, d9_17, d10_17, d11_17, d12_17, d13_17, d14_17, d15_17, d16_17, d17_19, d18_19, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
d1_17, d2_17, d3_17, d4_17, d5_17, d6_17, d7_17, d8_17, d9_17, d10_17, d11_17, d12_17, d13_17, d14_17, d15_17, d16_17, d17_20, d18_20, d19_20, 0, 0, 0, 0, 0, 0, 0, 0, 0,
d1_21, d2_21, d3_21, d4_21, d5_21, d6_21, d7_21, d8_21, d9_21, d10_21, d11_21, d12_21, d13_21, d14_21, d15_21, d16_21, d17_21, d18_21, d19_21, d20_21, 0, 0, 0, 0, 0, 0, 0, 0,
d1_21, d2_21, d3_21, d4_21, d5_21, d6_21, d7_21, d8_21, d9_21, d10_21, d11_21, d12_21, d13_21, d14_21, d15_21, d16_21, d17_21, d18_21, d19_21, d20_21, d21_22, 0, 0, 0, 0, 0, 0, 0,
d1_21, d2_21, d3_21, d4_21, d5_21, d6_21, d7_21, d8_21, d9_21, d10_21, d11_21, d12_21, d13_21, d14_21, d15_21, d16_21, d17_21, d18_21, d19_21, d20_21, d21_23, d22_23, 0, 0, 0, 0, 0, 0,
d1_21, d2_21, d3_21, d4_21, d5_21, d6_21, d7_21, d8_21, d9_21, d10_21, d11_21, d12_21, d13_21, d14_21, d15_21, d16_21, d17_21, d18_21, d19_21, d20_21, d21_24, d22_24, d23_24, 0, 0, 0, 0, 0,
d1_25, d2_25, d3_25, d4_25, d5_25, d6_25, d7_25, d8_25, d9_25, d10_25, d11_25, d12_25, d13_25, d14_25, d15_25, d16_25, d17_25, d18_25, d19_25, d20_25, d21_25, d22_25, d23_25, d24_25, 0, 0, 0, 0,
d1_25, d2_25, d3_25, d4_25, d5_25, d6_25, d7_25, d8_25, d9_25, d10_25, d11_25, d12_25, d13_25, d14_25, d15_25, d16_25, d17_25, d18_25, d19_25, d20_25, d21_25, d22_25, d23_25, d24_25, d25_26, 0, 0, 0,
d1_25, d2_25, d3_25, d4_25, d5_25, d6_25, d7_25, d8_25, d9_25, d10_25, d11_25, d12_25, d13_25, d14_25, d15_25, d16_25, d17_25, d18_25, d19_25, d20_25, d21_25, d22_25, d23_25, d24_25, d25_27, d26_27, 0, 0,
d1_25, d2_25, d3_25, d4_25, d5_25, d6_25, d7_25, d8_25, d9_25, d10_25, d11_25, d12_25, d13_25, d14_25, d15_25, d16_25, d17_25, d18_25, d19_25, d20_25, d21_25, d22_25, d23_25, d24_25, d25_28, d26_28, d27_28, 0)
delta = matrix(delta, nrow=nAlt, ncol=nAlt, byrow=TRUE)
emdc_settings = list(continuousChoice = continuousChoice,
avail = ones,
utilityOutside = 0,
utilities = utilities,
##budget = real_budget,
sigma = sigma,
gamma = gamma,
delta = delta,
cost = cost)
#maxIterations=100)

P[["model"]] = apollo_emdc(emdc_settings, functionality)

### Comment out as necessary
###P = apollo_panelProd(P, apollo_inputs, functionality)
P = apollo_prepareProb(P, apollo_inputs, functionality)
return(P)
}

# ################################################################# #
#### MODEL ESTIMATION ####
# ################################################################# #

model = apollo_estimate(apollo_beta, apollo_fixed, apollo_probabilities, apollo_inputs, estimate_settings=list(maxIterations=1000, estimationRoutine="BFGS"))

# ################################################################# #
#### MODEL OUTPUTS ####
# ################################################################# #

# ----------------------------------------------------------------- #
#---- FORMATTED OUTPUT (TO SCREEN) ----
# ----------------------------------------------------------------- #

apollo_modelOutput(model)

# ----------------------------------------------------------------- #
#---- FORMATTED OUTPUT (TO FILE, using model name) ----
# ----------------------------------------------------------------- #

apollo_saveOutput(model)

rm(predictions_base, predictions_new, change)

### Use the estimated model to make predictions
predictions_base = apollo_prediction(model, apollo_probabilities, apollo_inputs, prediction_settings=list(runs=3))

###Now imagine the cost for beefsoilh increasses by 1%
database$priceItem12=1.01*database$priceItem12
apollo_inputs=apollo_validateInputs()
predictions_new=apollo_prediction(model,apollo_probabilities,apollo_inputs)

###Return to original data
database$priceItem12=1/1.01*database$priceItem12
apollo_inputs=apollo_validateInputs()

###work with predictions at estimates
predictions_base=predictions_base[["at_estimates"]]

###Compute change in probabilities
##change=(predictions_new-predictions_base)/predictions_base

###Not interested in chosen alternative now,so drop last column
##change=change[,-ncol(change)]

###First two columns (change in ID and task) also not needed
##change=change[,-c(1,2)]

###Summary of changes (possible presence of NAs for unavailable alternatives)
##summary(change)

###Compute own elasticity:
log(sum(predictions_new[,14])/sum(predictions_base[,14]))/log(1.01)

###Compute cross-elasticities for other ingredients
log(sum(predictions_new[,11])/sum(predictions_base[,11]))/log(1.01)

log(sum(predictions_new[,12])/sum(predictions_base[,12]))/log(1.01)

log(sum(predictions_new[,13])/sum(predictions_base[,13]))/log(1.01)

[dpalma]

Hi Mohammed,

The eMDC is a model with several moving parts. In particular, the presence of complementarity/substitution parameters can make estimation tricky, as for certain combinations of values they can lead to decreasing marginal utilities. So it is always a good idea to first find a good model without any complementarity and substitution (delta) parameters, and only after you found that, start introducing them slowly. But let's look at your questions.

1. My model is very sensitive to the starting value of sigma. When I increase the starting value—for example, to 3—the model runs well and a maximum can be found.
This is common, specially for large datasets. It is very common for large models to fail due to the likelihood being too small and hitting the precision limit of the computr (i.e. likelihood ~ 0, so the computer takes it to be 0, and then log(0) produces a failure). Increasing the size of sigma can help with this problem, as also setting apollo_control$panelData = FALSE (to avoid multiplying the likelihood of the observations from teh same individual, if you have multiple observations per individual).

2. My model takes a very long time to run; the overall estimation process takes more than a day. Is this something I should be worried about?
The likelihood of the eMDC involves calculating the determinant of the jacobian of W (see paper). Unfortunately, that does not have any known compact (i.e. easy to write) form. Hence, calculation must be done observation per observation, which greatly slows down the process. I am afraid at the moment there is no way to go around that. But you can always use a bigger computer with more cores, and set apollo_control$nCores accordingly. In most computers apollo_control$nCores=3 should work fine, but if you have a bigger computer, you can set that number higher.

3. After obtaining the results, I conducted some predictions, and for some food products that are substitutes (e.g., conventional and organic beef), the cross-price elasticities turn out negative, which does not make sense. Do you think this might relate to the complexity of the eMDC model?
So you mean that when products A and B are substitute, and A increases price, then demand for B goes down? This is indeed weird. I would guess that the substitution effect is not strong enough to compensate for a third product C becoming comparatively more attractive. Does this happen for all substitutes or just for some? does the magnitude of the deltas have any influence on this?

Best wishes,
David