I would like to estimate a mixed-logit from SP and RP data. The data comes from different RP/SP individuals, the SP data has 4 choice tasks for each individual, the RP data is just a single observation. I do not have the socio-demographics for the RP observations so they should only be used to scale the model.
The model is estimating on a single core, but this takes obviously very long with the complexity. When specifying multiple cores, I am encounting issues. Exactly the same has happened with a standard MNL, but there the computation speed was not such an issue.
First approach, error: "checkForRemoteErrors(lapply(cl, recvResult)) :
10 nodes produced errors; first error: SPECIFICATION ISSUE - You attempted to call the function apollo_panelProd at a stage where the probabilities are already at the individual level!"
Code: Select all
##############################################################
#### 1. LOAD LIBRARY & CORE SETTINGS ####
##############################################################
library(apollo)
apollo_initialise()
apollo_control = list(
modelName = "MMNL_RP_SP",
modelDescr = "RP–SP mixed logit",
indivID = "RID",
nCores = 4, # adjust to your machine
outputDirectory = "output"
)
##############################################################
#### 2. LOAD DATA ####
##############################################################
database <- data
##############################################################
#### 3. DEFINE FIXED & MEAN/S.D. PARAMETERS ####
##############################################################
apollo_beta = c(
## ─── Random coefficients: means & st. devs ─────────────────
mu_asc_opt1 = 0, sigma_asc_opt1 = 1,
mu_asc_opt2 = 0, sigma_asc_opt2 = 1,
mu_asc_opt3 = 0, sigma_asc_opt3 = 1,
mu_log_b_price = -3, sigma_log_b_price = .5, # log-normal, keeps sign < 0
## ─── Fixed taste parameters ────────────────────────────────
b_opt3 = 0,
## SP interactions: gender
b_gender_opt1 = 0, b_gender_opt2 = 0, b_gender_opt3 = 0,
## SP interactions: city centre
b_cityCenter_opt1 = 0, b_cityCenter_opt2 = 0, b_cityCenter_opt3 = 0,
## SP interactions: age categories
b_age18_29_opt1 = 0, b_age18_29_opt2 = 0, b_age18_29_opt3 = 0,
b_age30_44_opt1 = 0, b_age30_44_opt2 = 0, b_age30_44_opt3 = 0,
b_age45_64_opt1 = 0, b_age45_64_opt2 = 0, b_age45_64_opt3 = 0,
## SP interactions: category segments 1–5
b_cat1_opt1 = 0, b_cat1_opt2 = 0, b_cat1_opt3 = 0,
b_cat2_opt1 = 0, b_cat2_opt2 = 0, b_cat2_opt3 = 0,
b_cat3_opt1 = 0, b_cat3_opt2 = 0, b_cat3_opt3 = 0,
b_cat4_opt1 = 0, b_cat4_opt2 = 0, b_cat4_opt3 = 0,
b_cat5_opt1 = 0, b_cat5_opt2 = 0, b_cat5_opt3 = 0,
## SP interactions: PT subscription dummy
b_ptsub_opt1 = 0, b_ptsub_opt2 = 0, b_ptsub_opt3 = 0,
## Scale parameters
mu_RP = 1, # fixed
mu_SP = 1 # to be estimated
)
apollo_fixed = c("mu_RP") # keep RP scale normalised
##############################################################
#### 4. RANDOM COEFFICIENTS & DRAWS ####
##############################################################
apollo_draws = list(
interDrawsType = "halton",
interNDraws = 500,
interNormDraws = c("d_asc_opt1","d_asc_opt2","d_asc_opt3","d_b_price"),
intraDrawsType = "none", intraNDraws = 0
)
apollo_randCoeff = function(apollo_beta, apollo_inputs){
rc = list()
## Normally-distributed ASCs
rc[["asc_opt1"]] = mu_asc_opt1 + sigma_asc_opt1 * d_asc_opt1
rc[["asc_opt2"]] = mu_asc_opt2 + sigma_asc_opt2 * d_asc_opt2
rc[["asc_opt3"]] = mu_asc_opt3 + sigma_asc_opt3 * d_asc_opt3
## Log-normal (negative) price coefficient
rc[["b_price"]] = -exp(mu_log_b_price + sigma_log_b_price * d_b_price)
return(rc)
}
##############################################################
#### 5. GROUP & VALIDATE INPUTS ####
##############################################################
apollo_inputs = apollo_validateInputs()
##############################################################
#### 6. LIKELIHOOD FUNCTION ####
##############################################################
apollo_probabilities = function(apollo_beta, apollo_inputs, functionality="estimate"){
apollo_attach(apollo_beta, apollo_inputs)
on.exit(apollo_detach(apollo_beta, apollo_inputs))
P = list() # container
## ── (1) SP-only ASC terms with interactions ─────────────────
asc_opt1_val = asc_opt1
asc_opt2_val = asc_opt2
asc_opt3_val = asc_opt3
if(functionality=="estimate"){ # only needed when estimating
asc_opt1_val = asc_opt1_val +
ifelse(rp==0, b_gender_opt1 * gender, 0) +
ifelse(rp==0, b_cityCenter_opt1* city_center, 0) +
ifelse(rp==0, b_age18_29_opt1 * age_cat_18_29, 0) +
ifelse(rp==0, b_age30_44_opt1 * age_cat_30_44, 0) +
ifelse(rp==0, b_age45_64_opt1 * age_cat_45_64, 0) +
ifelse(rp==0, b_cat1_opt1 * category_1, 0) +
ifelse(rp==0, b_cat2_opt1 * category_2, 0) +
ifelse(rp==0, b_cat3_opt1 * category_3, 0) +
ifelse(rp==0, b_cat4_opt1 * category_4, 0) +
ifelse(rp==0, b_cat5_opt1 * category_5, 0) +
ifelse(rp==0, b_ptsub_opt1 * ptsub, 0)
asc_opt2_val = asc_opt2_val +
ifelse(rp==0, b_gender_opt2 * gender, 0) +
ifelse(rp==0, b_cityCenter_opt2* city_center, 0) +
ifelse(rp==0, b_age18_29_opt2 * age_cat_18_29, 0) +
ifelse(rp==0, b_age30_44_opt2 * age_cat_30_44, 0) +
ifelse(rp==0, b_age45_64_opt2 * age_cat_45_64, 0) +
ifelse(rp==0, b_cat1_opt2 * category_1, 0) +
ifelse(rp==0, b_cat2_opt2 * category_2, 0) +
ifelse(rp==0, b_cat3_opt2 * category_3, 0) +
ifelse(rp==0, b_cat4_opt2 * category_4, 0) +
ifelse(rp==0, b_cat5_opt2 * category_5, 0) +
ifelse(rp==0, b_ptsub_opt2 * ptsub, 0)
asc_opt3_val = asc_opt3_val +
ifelse(rp==0, b_gender_opt3 * gender, 0) +
ifelse(rp==0, b_cityCenter_opt3* city_center, 0) +
ifelse(rp==0, b_age18_29_opt3 * age_cat_18_29, 0) +
ifelse(rp==0, b_age30_44_opt3 * age_cat_30_44, 0) +
ifelse(rp==0, b_age45_64_opt3 * age_cat_45_64, 0) +
ifelse(rp==0, b_cat1_opt3 * category_1, 0) +
ifelse(rp==0, b_cat2_opt3 * category_2, 0) +
ifelse(rp==0, b_cat3_opt3 * category_3, 0) +
ifelse(rp==0, b_cat4_opt3 * category_4, 0) +
ifelse(rp==0, b_cat5_opt3 * category_5, 0) +
ifelse(rp==0, b_ptsub_opt3 * ptsub, 0)
}
## ── (2) Utilities ───────────────────────────────────────────
V = list(
opt1 = asc_opt1_val + b_price * price_opt1,
opt2 = asc_opt2_val + b_price * price_opt2,
opt3 = asc_opt3_val + b_price * price_opt3 + b_opt3 * opt3,
opt4 = 0
)
## ── (3) RP model block ──────────────────────────────────────
mnl_RP = list(
alternatives = c(opt1=1, opt2=2, opt3=3, opt4=4),
avail = list(opt1=avail_opt1, opt2=avail_opt2, opt3=avail_opt3, opt4=1),
choiceVar = pref1,
utilities = lapply(V, function(x) mu_RP * x),
rows = (rp==1)
)
P[["RP"]] = apollo_mnl(mnl_RP, functionality)
## ── (4) SP model block ──────────────────────────────────────
mnl_SP = list(
alternatives = c(opt1=1, opt2=2, opt3=3, opt4=4),
avail = list(opt1=avail_opt1, opt2=avail_opt2, opt3=avail_opt3, opt4=1),
choiceVar = pref1,
utilities = lapply(V, function(x) mu_SP * x),
rows = (rp==0)
)
P[["SP"]] = apollo_mnl(mnl_SP, functionality)
## ── (5) Combine, panel product, average draws, return ───────
P = apollo_combineModels(P, apollo_inputs, functionality)
P = apollo_panelProd( P, apollo_inputs, functionality)
P = apollo_avgInterDraws( P, apollo_inputs, functionality)
P = apollo_prepareProb( P, apollo_inputs, functionality)
return(P)
}
##############################################################
#### 7. ESTIMATE ####
##############################################################
model = apollo_estimate(
apollo_beta, apollo_fixed,
apollo_probabilities, apollo_inputs
)
##############################################################
#### 8. OUTPUT ####
##############################################################
apollo_modelOutput(model)
apollo_saveOutput(model)SPECIFICATION ISSUE - Observations from the same individual must be combined (i.e. multiplied) before averaging over inter-individual draws."
Code: Select all
##############################################################
#### 1. LIBRARIES & CORE CONTROLS ####
##############################################################
library(apollo)
apollo_initialise()
apollo_control = list(
modelName = "RP–SP mixed logit",
modelDescr = "RP–SP mixed logit",
indivID = "RID",
nCores = 4,
outputDirectory = "output"
)
##############################################################
#### 2. LOAD DATA ####
##############################################################
database <- data
##############################################################
#### 3. STARTING VALUES ####
##############################################################
apollo_beta = c(
## random-coeff means & s.d.'s
mu_asc_opt1 = 0, sigma_asc_opt1 = 1,
mu_asc_opt2 = 0, sigma_asc_opt2 = 1,
mu_asc_opt3 = 0, sigma_asc_opt3 = 1,
mu_log_b_price = -3, sigma_log_b_price = .5,
## fixed taste parameters
b_opt3 = 0,
b_gender_opt1 = 0, b_gender_opt2 = 0, b_gender_opt3 = 0,
b_cityCenter_opt1 = 0, b_cityCenter_opt2 = 0, b_cityCenter_opt3 = 0,
b_age18_29_opt1 = 0, b_age18_29_opt2 = 0, b_age18_29_opt3 = 0,
b_age30_44_opt1 = 0, b_age30_44_opt2 = 0, b_age30_44_opt3 = 0,
b_age45_64_opt1 = 0, b_age45_64_opt2 = 0, b_age45_64_opt3 = 0,
b_cat1_opt1 = 0, b_cat1_opt2 = 0, b_cat1_opt3 = 0,
b_cat2_opt1 = 0, b_cat2_opt2 = 0, b_cat2_opt3 = 0,
b_cat3_opt1 = 0, b_cat3_opt2 = 0, b_cat3_opt3 = 0,
b_cat4_opt1 = 0, b_cat4_opt2 = 0, b_cat4_opt3 = 0,
b_cat5_opt1 = 0, b_cat5_opt2 = 0, b_cat5_opt3 = 0,
b_ptsub_opt1 = 0, b_ptsub_opt2 = 0, b_ptsub_opt3 = 0,
## scales
mu_RP = 1,
mu_SP = 1
)
apollo_fixed = c("mu_RP") # keep RP scale normalised
##############################################################
#### 4. DRAWS & RANDOM COEFFICIENTS ####
##############################################################
apollo_draws = list(
interDrawsType = "halton",
interNDraws = 500,
interNormDraws = c("d_asc_opt1","d_asc_opt2","d_asc_opt3","d_b_price"),
intraDrawsType = "none",
intraNDraws = 0
)
apollo_randCoeff = function(apollo_beta, apollo_inputs){
rc = list()
rc[["asc_opt1"]] = mu_asc_opt1 + sigma_asc_opt1 * d_asc_opt1
rc[["asc_opt2"]] = mu_asc_opt2 + sigma_asc_opt2 * d_asc_opt2
rc[["asc_opt3"]] = mu_asc_opt3 + sigma_asc_opt3 * d_asc_opt3
rc[["b_price"]] = -exp(mu_log_b_price + sigma_log_b_price * d_b_price) # log-normal <0
return(rc)
}
##############################################################
#### 5. VALIDATE INPUTS ####
##############################################################
apollo_inputs = apollo_validateInputs()
##############################################################
#### 6. LIKELIHOOD ####
##############################################################
apollo_probabilities = function(apollo_beta, apollo_inputs,
functionality = "estimate"){
apollo_attach(apollo_beta, apollo_inputs)
on.exit(apollo_detach(apollo_beta, apollo_inputs))
P = list()
## --- 6.1 SP-only ASC adjustments ---------------------------
asc_opt1_val = asc_opt1
asc_opt2_val = asc_opt2
asc_opt3_val = asc_opt3
if(functionality == "estimate"){
asc_opt1_val = asc_opt1_val +
ifelse(rp==0, b_gender_opt1 * gender, 0) +
ifelse(rp==0, b_cityCenter_opt1* city_center, 0) +
ifelse(rp==0, b_age18_29_opt1 * age_cat_18_29,0) +
ifelse(rp==0, b_age30_44_opt1 * age_cat_30_44,0) +
ifelse(rp==0, b_age45_64_opt1 * age_cat_45_64,0) +
ifelse(rp==0, b_cat1_opt1 * category_1, 0) +
ifelse(rp==0, b_cat2_opt1 * category_2, 0) +
ifelse(rp==0, b_cat3_opt1 * category_3, 0) +
ifelse(rp==0, b_cat4_opt1 * category_4, 0) +
ifelse(rp==0, b_cat5_opt1 * category_5, 0) +
ifelse(rp==0, b_ptsub_opt1 * ptsub, 0)
asc_opt2_val = asc_opt2_val +
ifelse(rp==0, b_gender_opt2 * gender, 0) +
ifelse(rp==0, b_cityCenter_opt2* city_center, 0) +
ifelse(rp==0, b_age18_29_opt2 * age_cat_18_29,0) +
ifelse(rp==0, b_age30_44_opt2 * age_cat_30_44,0) +
ifelse(rp==0, b_age45_64_opt2 * age_cat_45_64,0) +
ifelse(rp==0, b_cat1_opt2 * category_1, 0) +
ifelse(rp==0, b_cat2_opt2 * category_2, 0) +
ifelse(rp==0, b_cat3_opt2 * category_3, 0) +
ifelse(rp==0, b_cat4_opt2 * category_4, 0) +
ifelse(rp==0, b_cat5_opt2 * category_5, 0) +
ifelse(rp==0, b_ptsub_opt2 * ptsub, 0)
asc_opt3_val = asc_opt3_val +
ifelse(rp==0, b_gender_opt3 * gender, 0) +
ifelse(rp==0, b_cityCenter_opt3* city_center, 0) +
ifelse(rp==0, b_age18_29_opt3 * age_cat_18_29,0) +
ifelse(rp==0, b_age30_44_opt3 * age_cat_30_44,0) +
ifelse(rp==0, b_age45_64_opt3 * age_cat_45_64,0) +
ifelse(rp==0, b_cat1_opt3 * category_1, 0) +
ifelse(rp==0, b_cat2_opt3 * category_2, 0) +
ifelse(rp==0, b_cat3_opt3 * category_3, 0) +
ifelse(rp==0, b_cat4_opt3 * category_4, 0) +
ifelse(rp==0, b_cat5_opt3 * category_5, 0) +
ifelse(rp==0, b_ptsub_opt3 * ptsub, 0)
}
## --- 6.2 Utilities -----------------------------------------
V = list(
opt1 = asc_opt1_val + b_price * price_opt1,
opt2 = asc_opt2_val + b_price * price_opt2,
opt3 = asc_opt3_val + b_price * price_opt3 + b_opt3 * opt3,
opt4 = 0
)
## --- 6.3 RP likelihood (one row per person) ----------------
mnl_RP = list(
alternatives = c(opt1=1, opt2=2, opt3=3, opt4=4),
avail = list(opt1=avail_opt1, opt2=avail_opt2, opt3=avail_opt3, opt4=1),
choiceVar = pref1,
utilities = lapply(V, function(x) mu_RP * x),
rows = (rp == 1)
)
P[["RP"]] = apollo_mnl(mnl_RP, functionality)
## --- 6.4 SP likelihood (panel) -----------------------------
mnl_SP = list(
alternatives = c(opt1=1, opt2=2, opt3=3, opt4=4),
avail = list(opt1=avail_opt1, opt2=avail_opt2, opt3=avail_opt3, opt4=1),
choiceVar = pref1,
utilities = lapply(V, function(x) mu_SP * x),
rows = (rp == 0)
)
P[["SP"]] = apollo_mnl(mnl_SP, functionality)
## --- 6.5 Combine RP & SP (still obs × draw) ----------------
P = apollo_combineModels(P, apollo_inputs, functionality)
## --- 6.6 Collapse across panel observations ----------------
if(any(apollo_inputs$rowsPerIndiv > 1)){
P = apollo_panelProd(P, apollo_inputs, functionality)
}
## --- 6.7 Average draws & tidy ------------------------------
P = apollo_avgInterDraws(P, apollo_inputs, functionality)
P = apollo_prepareProb(P, apollo_inputs, functionality)
return(P)
}
##############################################################
#### 7. ESTIMATE ####
##############################################################
model = apollo_estimate(
apollo_beta, apollo_fixed,
apollo_probabilities, apollo_inputs
)
##############################################################
#### 8. OUTPUT ####
##############################################################
apollo_modelOutput(model)
apollo_saveOutput(model)