eMDC Starting Values

[Jiayu Sun]

Hello,

I am estimating a model with 21 alternatives, which correspond to 21 dishes on the menu. I would like to estimate the eMDC model with full substitution and complementarity effects. The code runs well when I set all the interaction terms between the 21 dishes to a fixed value of 0. However, after I remove these interactions from the fixed value list, the model returns the following error:

Error in apollo_estimate(apollo_beta, apollo_fixed, apollo_probabilities, :
CALCULATION ISSUE - Log-likelihood calculation fails at values close to the starting values!
In addition: Warning messages:
1: In log(Gi) : NaNs produced
2: In log(Gi) : NaNs produced
3: In log(Gi) : NaNs produced

I have tried multiple ways to find the starting values for the interactions, but they all did work. Could you please provide any suggestions on this?

the code is here:

# ################################################################# #
#### LOAD LIBRARY AND DEFINE CORE SETTINGS ####
# ################################################################# #

### Clear memory
rm(list = ls())

### Load Apollo library
library(apollo)

### Set working directory (only works in RStudio)
##apollo_setWorkDir()

### Initialise code
apollo_initialise()

### Set core controls
apollo_control = list(
modelDescr = "Fast",
indivID = "Q_ResponseId",
outputDirectory = "output")

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


### Load data from within the Apollo package
database = Fast

### 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(aincentive=0.066,
sigma = 2.469,
atax=0.219,
# Satiation
g1 = 0.942,
g2 = 1.033,
g3 = 0.995,
g4 = 0.980,
g5 = 1.053,
g6 = 1.016,
g7 = 0.999,
g8 = 1.083,
g9 = 1.163,
g10 = 0.778,
g11 = 0.779,
g12 = 0.732,
g13 = 0.788,
g14 = 0.625,
g15 = 0.412,
g16 = 0.719,
g17 = 0.559,
g18 = 0.574,
g19 = 0.546,
g20 = 0.630,
g21 = 0.683,
# Base utility
b1 = -2.794,
b2 = -3.298,
b3 = -2.772,
b4 = -3.912,
b5 = -3.172,
b6 = -4.422,
b7 = -4.007,
b8 = -4.962,
b9 = -3.398,
b10 = -3.221,
b11 = -2.954,
b12 = -3.666,
b13 = -3.413,
b14 = -2.528,
b15 = -0.773,
b16 = -5.021,
b17 = -2.508,
b18 = -3.479,
b19 = -2.745,
b20 = -3.807,
b21 = -4.485,
btax1 =-0.296,
btax2 =-0.396,
btax3 =-0.388,
btax4 =-0.368,
btax5 =-0.233,
btax6 =-0.075,
btax7 =-0.124,
btax8 =0.099,
btax9 =-0.090,
btax10 =-0.408,
btax11 =0.007,
btax12 =-0.239,
btax13 =-0.122,
btax14 =-0.317,
btax15 =-0.599,
btax16 =0.148,
btax17 =-0.297,
btax18 =-0.348,
btax19 =-0.595,
btax20 =-0.004,
btax21 =-0.172,
# Compl/subst
d1_2=0,
d1_3=0,
d1_4=0,
d1_5=0,
d1_6=0,
d1_7=0,
d1_8=0,
d1_9=0,
d1_10=0,
d1_11=0,
d1_12=0,
d1_13=0,
d1_14=0,
d1_15=0,
d1_16=0,
d1_17=0,
d1_18=0,
d1_19=0,
d1_20=0,
d1_21=0,
d2_3=0,
d2_4=0,
d2_5=0,
d2_6=0,
d2_7=0,
d2_8=0,
d2_9=0,
d2_10=0,
d2_11=0,
d2_12=0,
d2_13=0,
d2_14=0,
d2_15=0,
d2_16=0,
d2_17=0,
d2_18=0,
d2_19=0,
d2_20=0,
d2_21=0,
d3_4=-0,
d3_5=0,
d3_6=-0,
d3_7=-0,
d3_8=-0,
d3_9=0,
d3_10=0,
d3_11=0,
d3_12=0,
d3_13=0,
d3_14=0,
d3_15=0,
d3_16=-0,
d3_17=0,
d3_18=0,
d3_19=0,
d3_20=0,
d3_21=0,
d4_5=-0,
d4_6=0,
d4_7=0,
d4_8=0,
d4_9=-0,
d4_10=0,
d4_11=0,
d4_12=0,
d4_13=0,
d4_14=0,
d4_15=0,
d4_16=0,
d4_17=0,
d4_18=0,
d4_19=0,
d4_20=0,
d4_21=0,
d5_6=-0,
d5_7=-0,
d5_8=-0,
d5_9=-0,
d5_10=0,
d5_11=0,
d5_12=0,
d5_13=0,
d5_14=0,
d5_15=0,
d5_16=0,
d5_17=0,
d5_18=0,
d5_19=0,
d5_20=0,
d5_21=0,
d6_7=0,
d6_8=0,
d6_9=-0,
d6_10=0,
d6_11=0,
d6_12=0,
d6_13=0,
d6_14=0,
d6_15=0,
d6_16=0,
d6_17=0,
d6_18=0,
d6_19=0,
d6_20=0,
d6_21=0,
d7_8=0,
d7_9=-0,
d7_10=0,
d7_11=0,
d7_12=0,
d7_13=0,
d7_14=0,
d7_15=0,
d7_16=-0,
d7_17=0,
d7_18=0,
d7_19=0,
d7_20=0,
d7_21=0,
d8_9=-0,
d8_10=0,
d8_11=0,
d8_12=0,
d8_13=0,
d8_14=0,
d8_15=0,
d8_16=0,
d8_17=-0,
d8_18=0,
d8_19=0,
d8_20=0,
d8_21=0,
d9_10=0,
d9_11=0,
d9_12=0,
d9_13=0,
d9_14=0,
d9_15=-0,
d9_16=0,
d9_17=0,
d9_18=0,
d9_19=0,
d9_20=0,
d9_21=0,
d10_11=-0,
d10_12=-0,
d10_13=-0,
d10_14=0,
d10_15=0,
d10_16=0,
d10_17=0,
d10_18=0,
d10_19=0,
d10_20=0,
d10_21=0,
d11_12=0,
d11_13=0,
d11_14=0,
d11_15=0,
d11_16=0,
d11_17=0,
d11_18=0,
d11_19=0,
d11_20=0,
d11_21=0,
d12_13=0,
d12_14=0,
d12_15=0,
d12_16=0,
d12_17=0,
d12_18=0,
d12_19=0,
d12_20=0,
d12_21=0,
d13_14=0,
d13_15=0,
d13_16=0,
d13_17=0,
d13_18=0,
d13_19=0,
d13_20=0,
d13_21=0,
d14_15=0,
d14_16=0,
d14_17=0,
d14_18=0,
d14_19=0,
d14_20=0,
d14_21=0,
d15_16=0,
d15_17=0,
d15_18=0,
d15_19=0,
d15_20=0,
d15_21=0,
d16_17=0,
d16_18=0,
d16_19=0,
d16_20=0,
d16_21=0,
d17_18=0,
d17_19=0,
d17_20=0,
d17_21=0,
d18_19=0,
d18_20=0,
d18_21=0,
d19_20=0,
d19_21=0,
d20_21=0)
### Names of fixed parameters
apollo_fixed = c()
# ################################################################# #
#### 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("item1", "item2", "item3", "item4", "item5","item6", "item7", "item8", "item9", "item10", "item11", "item12", "item13", "item14", "item15","item16", "item17", "item18", "item19", "item20","item21" )
nAlt = length(alts)
ones = setNames(as.list(rep(1, nAlt)), alts)
continuousChoice = list(item1=item1/1,
item2=item2/1,
item3=item3/1,
item4=item4/1,
item5=item5/1,
item6=item6/1,
item7=item7/1,
item8=item8/1,
item9=item9/1,
item10=item10/1,
item11=item11/1,
item12=item12/1,
item13=item13/1,
item14=item14/1,
item15=item15/1,
item16=item16/1,
item17=item17/1,
item18=item18/1,
item19=item19/1,
item20=item20/1,
item21=item21/1)

cost = list(item1=p1/1,
item2=p2/1,
item3=p3/1,
item4=p4/1,
item5=p5/1,
item6=p6/1,
item7=p7/1,
item8=p8/1,
item9=p9/1,
item10=p10/1,
item11=p11/1,
item12=p12/1,
item13=p13/1,
item14=p14/1,
item15=p15/1,
item16=p16/1,
item17=p17/1,
item18=p18/1,
item19=p19/1,
item20=p20/1,
item21=p21/1)

utilities = list(
item1 = b1 + btax1*tax,
item2 = b2 + btax2*tax,
item3 = b3 + btax3*tax,
item4 = b4 + btax4*tax,
item5 = b5 + btax5*tax,
item6 = b6 + btax6*tax,
item7 = b7 + btax7*tax,
item8 = b8 + btax8*tax,
item9 = b9 + btax9*tax,
item10 = b10 + btax10*tax,
item11 = b11 + btax11*tax,
item12 = b12 + btax12*tax,
item13 = b13 + btax13*tax,
item14 = b14 + btax14*tax,
item15 = b15 + btax15*tax,
item16 = b16 + btax16*tax,
item17 = b17 + btax17*tax,
item18 = b18 + btax18*tax,
item19 = b19 + btax19*tax,
item20 = b20 + btax20*tax,
item21 = b21 + btax21*tax
)
gamma = list( item1 = g1,
item2 = g2,
item3 = g3,
item4 = g4,
item5 = g5,
item6 = g6,
item7 = g7,
item8 = g8,
item9 = g9,
item10 = g10,
item11 = g11,
item12 = g12,
item13 = g13,
item14 = g14,
item15 = g15,
item16 = g16,
item17 = g17,
item18 = g18,
item19 = g19,
item20 = g20,
item21 = g21

)
delta = c(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,
d1_3, d2_3, 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,
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,
d1_6, d2_6, d3_6, d4_6, d5_6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
d1_7, d2_7, d3_7, d4_7, d5_7, d6_7, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
d1_8, d2_8, d3_8, d4_8, d5_8, d6_8, d7_8, 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,
d1_10, d2_10, d3_10, d4_10, d5_10, d6_10, d7_10, d8_10, d9_10, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
d1_11, d2_11, d3_11, d4_11, d5_11, d6_11, d7_11, d8_11, d9_11, d10_11, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
d1_12, d2_12, d3_12, d4_12, d5_12, d6_12, d7_12, d8_12, d9_12, d10_12, d11_12, 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,
d1_14, d2_14, d3_14, d4_14, d5_14, d6_14, d7_14, d8_14, d9_14, d10_14, d11_14, d12_14, d13_14, 0, 0, 0, 0, 0, 0, 0, 0,
d1_15, d2_15, d3_15, d4_15, d5_15, d6_15, d7_15, d8_15, d9_15, d10_15, d11_15, d12_15, d13_15, d14_15, 0, 0, 0, 0, 0, 0, 0,
d1_16, d2_16, d3_16, d4_16, d5_16, d6_16, d7_16, d8_16, d9_16, d10_16, d11_16, d12_16, d13_16, d14_16, d15_16, 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,
d1_18, d2_18, d3_18, d4_18, d5_18, d6_18, d7_18, d8_18, d9_18, d10_18, d11_18, d12_18, d13_18, d14_18, d15_18, d16_18, d17_18, 0, 0, 0, 0,
d1_19, d2_19, d3_19, d4_19, d5_19, d6_19, d7_19, d8_19, d9_19, d10_19, d11_19, d12_19, d13_19, d14_19, d15_19, d16_19, d17_19, d18_19, 0, 0, 0,
d1_20, d2_20, d3_20, d4_20, d5_20, d6_20, d7_20, d8_20, d9_20, d10_20, d11_20, d12_20, d13_20, d14_20, d15_20, d16_20, d17_20, d18_20, d19_20, 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)
delta = matrix(delta, nrow=nAlt, ncol=nAlt, byrow=TRUE)
emdc_settings = list(continuousChoice = continuousChoice,
avail = ones,
utilityOutside = aincentive*tenper+atax*tax,
utilities = utilities,
budget = real_budget,
sigma = sigma,
gamma = gamma,
delta = delta,
cost = cost,
maxIterations=20)

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)

[dpalma]

Hi,

Finding good starting values for eMDC si not easy. The model is quite sensitive to the values of delta (complementarity/substitution) parameters, so it cat easily fail for values that are not close to the optimum. The failure is caused because before starting estimation, Apollo adds a small value to each parameter and tries to evaluate apollo_probabilities. The problem is that even that small value is too large for the model, and it makes it fail.

I would recommend the following:

• Start with a simplified delta matrix first. You could, for example, assume all delta parameters to be equal across dishes of the same type. Let's say there are 7 meat (M) dishes, 7 fish (F) dishes, and 7 vegetarian (V) dishes, then you could assume that there are just six delta parameters: dMM, dFM, dFF, dVM, dVF, and dVV. Once you manage to estimate this model, you can try making the delta matrix more complicated, as then you will have better starting values (the estimates of the simple model) to use.
• If you have more than one observation per individual, I recommend you set apollo_control$panelData = FALSE. This can help reduce the number of cases when the likelihood becomes zero.
• You can also try setting apollo_control$workInLogs = TRUE, but then you will also have to set estimate_settings$estimationRoutine="BFGS" when calling apollo_estimate.

Finally, I see that you are using the eMDC with observed budget. You may also want to try the version with unobserved budget. To use this version of the model, you simply need to remove the definition of the budget from emdc_settings. The benefit of the unobserved budget is that you do not need to know the budget to forecast. This is useful if the main objective of your modelling is forecasting. If instead you just want to understand how people choose, then using the version with observed budget is fine.

Hope this helps.

Bets wishes,
David

[Jiayu Sun]

Hello David,

Thank you for your responses. They’ve been very helpful.

I have a follow-up question related to the model. After estimating it, I used the validation data to make predictions with the following code. The output I received includes the following sets of results:

1) outside, item1–item21

2) outside.1, item1.1–item21.1

3) outside.2, item1.2–item21.2

4) outside.3, item1.3–item21.3

Could you please provide any insights on what these repeated groups (with .1, .2, .3) represent?

Also, is it possible for me get the standard errors of the predictions? Thank you!

Code:
model = apollo_loadModel(apollo_inputs$apollo_control$modelName)
apollo_inputs = apollo_validateInputs(database=dbVal)
database=dbVal
apollo_inputs$apollo_control$nCores = 4
pred = apollo_prediction(model, apollo_probabilities, apollo_inputs)
XObs = cbind(item1=dbVal$item1/1,
item2=dbVal$item2/1,
item3=dbVal$item3/1,
item4=dbVal$item4/1,
item5=dbVal$item5/1,
item6=dbVal$item6/1,
item7=dbVal$item7/1,
item8=dbVal$item8/1,
item9=dbVal$item9/1,
item10=dbVal$item10/1,
item11=dbVal$item11/1,
item12=dbVal$item12/1,
item13=dbVal$item13/1,
item14=dbVal$item14/1,
item15=dbVal$item15/1,
item16=dbVal$item16/1,
item17=dbVal$item17/1,
item18=dbVal$item18/1,
item19=dbVal$item19/1,
item20=dbVal$item20/1,
item21=dbVal$item21/1)


XPre = pred[,4:24]

round(sqrt(colMeans((XObs - XPre)^2)),2)
round(sqrt(mean((colSums(XObs) - colSums(XPre))^2)),2)

Output:

outside item1 item2 item3 item4 item5 item6 item7 item8 item9 item10 item11 item12 item13
Aggregate 70887.49 2298.95 1286.49 1561.79 1590.02 1437.71 2112.23 1694.04 2312.61 1018.54 1689.79 2061.29 3489.10 1653.89
Average 15.27 0.50 0.28 0.34 0.34 0.31 0.46 0.37 0.50 0.22 0.36 0.44 0.75 0.36
item14 item15 item16 item17 item18 item19 item20 item21 outside.1 item1.1 item2.1 item3.1 item4.1 item5.1
Aggregate 3215.66 9975.14 4459.89 3260.42 4442.20 6887.68 4219.61 2969.56 79610.21 7280.63 4708.33 5052.36 6249.99 4982.75
Average 0.69 2.15 0.96 0.70 0.96 1.48 0.91 0.64 17.15 1.57 1.01 1.09 1.35 1.07
item6.1 item7.1 item8.1 item9.1 item10.1 item11.1 item12.1 item13.1 item14.1 item15.1 item16.1 item17.1 item18.1
Aggregate 8597.18 6620.33 10095.31 3849.72 6051.69 6723.84 12487.58 6099.66 9467.60 20505.08 18879.84 9690.27 15392.43
Average 1.85 1.43 2.18 0.83 1.30 1.45 2.69 1.31 2.04 4.42 4.07 2.09 3.32
item19.1 item20.1 item21.1 outside.2 item1.2 item2.2 item3.2 item4.2 item5.2 item6.2 item7.2 item8.2 item9.2 item10.2
Aggregate 20808.25 15167.96 12316.56 4641 4640.32 4637.12 4634.06 4640.98 4636.1 4641 4640.9 4641 4632.98 4640.06
Average 4.48 3.27 2.65 1 1.00 1.00 1.00 1.00 1.0 1 1.0 1 1.00 1.00
item11.2 item12.2 item13.2 item14.2 item15.2 item16.2 item17.2 item18.2 item19.2 item20.2 item21.2 outside.3 item1.3
Aggregate 4640.16 4641 4640.54 4640.9 4641 4641 4640.86 4641 4641 4641 4640.98 0 4.99
Average 1.00 1 1.00 1.0 1 1 1.00 1 1 1 1.00 0 0.00
item2.3 item3.3 item4.3 item5.3 item6.3 item7.3 item8.3 item9.3 item10.3 item11.3 item12.3 item13.3 item14.3 item15.3
Aggregate 23.41 37.92 0.28 26.92 0 1.02 0 44.29 6.69 6.31 0 3.66 0.94 0
Average 0.01 0.01 0.00 0.01 0 0.00 0 0.01 0.00 0.00 0 0.00 0.00 0
item16.3 item17.3 item18.3 item19.3 item20.3 item21.3
Aggregate 0 1.36 0 0 0 0.28
Average 0 0.00 0 0 0 0.00

[dpalma]

Hi,

The forecast contains the following columns (in the following order):

• Expected consumption of outside, then all inside goods.
• The s.d. of the previous values
• Expected probability of consuming (i.e. non-zero consumption) of each good
• The s.d. of the previous values.

Thanks for noticing the names were not informative, we will use better names for the columns in the next version of Apollo.

Best wishes,
David