I have a question regarding apollo_llCalc in Bayesian estimation. I have a model with two components, both of them are MNL models. Whenever I try to run the function apollo_llCalc I get the following error message:
Code: Select all
Calculating LL of each model component...Not applicable to all components.
[1] NA
So I tried this:It should be noted that when calling this function, the format of apollo_beta needs to be compatible with what is used inside apollo_probabilities. All parameters used inside the model need to be included, with either one value per parameter (classical estimation) or one value per parameter and per observation (Bayesian estimation).
Code: Select all
llCalc.apollo_beta = as.list(apollo_beta)
llCalc.apollo_beta = lapply(llCalc.apollo_beta, rep, times=nrow(database))
apollo_llCalc(apollo_beta=llCalc.apollo_beta,
apollo_probabilities,
apollo_inputs)
Code: Select all
if(!silent) cat("Calculating LL of each model component...")
Pout <- tryCatch(apollo_probabilities(apollo_beta, apollo_inputs, functionality="output"),
error=function(e) return(NA))
if(!anyNA(Pout) && is.list(Pout)){
Code: Select all
apollo_probabilities(apollo_beta=llCalc.apollo_beta, apollo_inputs, functionality = "output")
Is it a bug or am I doing something wrong?
Here is the whole code of the MNL model with two model components:
Code: Select all
# ################################################################# #
#### LOAD LIBRARY AND DEFINE CORE SETTINGS ####
# ################################################################# #
### Clear memory
rm(list = ls())
### Load Apollo library
library(apollo)
### Initialise code
apollo_initialise()
### Set core controls
apollo_control = list(
modelName = "MNL_COV1_B",
modelDescr = "MNL model in utility space (Bayesian estimation)",
indivID ="sys_RespNum",
nCores = 1,
mixing = FALSE,
HB = TRUE,
workinlogs = FALSE
)
# ################################################################# #
#### LOAD DATA AND APPLY ANY TRANSFORMATIONS ####
# ################################################################# #
input.file <- "./input/Survey_SWP_modeldataprep.wide.rds"
database <- readRDS(input.file)
# ################################################################# #
#### DEFINE MODEL PARAMETERS ####
# ################################################################# #
apollo_beta = c(
# ----------------------------------------------------------------- #
#---- Choice Model
# ----------------------------------------------------------------- #
### Alternative specific constants (ASC)
b_asc_1 = 0,
b_asc_2 = 0,
b_asc_3 = 0,
### Mix_Default
b_Mix_Default = 0,
### Mix_Green
b_Mix_Green = 0,
b_Gender_Mix_Green = 0,
b_Age_Mix_Green = 0,
b_Education_Mix_Green = 0,
### Regio_0
b_Regio_0 = 0,
### Regio_50
b_Regio_50 = 0,
b_Gender_Regio_50 = 0,
b_Age_Regio_50 = 0,
b_Education_Regio_50 = 0,
### Regio_100
b_Regio_100 = 0,
b_Gender_Regio_100 = 0,
b_Age_Regio_100 = 0,
b_Education_Regio_100 = 0,
### Label_TUVSud
b_Label_TUVSud = 0,
### Label_OKPower
b_Label_OKPower = 0,
### Label_GrunerStrom
b_Label_GrunerStrom = 0,
### Label_None
b_Label_None = 0,
### Price
b_Price = 0,
b_Income_Price = 0,
b_PriceMonthly_Price = 0,
# ----------------------------------------------------------------- #
#---- Choice Model DR
# ----------------------------------------------------------------- #
### DR_asc_0
b_DR_asc_0 = 0,
### DR_asc_1
b_DR_asc_1 = 0,
b_Age_DR_asc_1 = 0,
b_Education_DR_asc_1 = 0,
b_Income_DR_asc_1 = 0,
b_PriceMonthly_DR_asc_1 = 0
)
### Vector with names (in quotes) of parameters to be kept fixed at their starting value in apollo_beta, use apollo_beta_fixed = c() if none
apollo_fixed = c(
### Choice Model
"b_asc_3",
"b_Mix_Default",
"b_Regio_0",
"b_Label_TUVSud",
### Choice Model DR
"b_DR_asc_0"
)
### Read estimates from other models
# apollo_beta = apollo_readBeta(apollo_beta, apollo_fixed, "flexQgrid_MIXL_FM1_priorHCM_Sobol2000", overwriteFixed=FALSE)
# apollo_beta_new = replace(apollo_beta, names(priors.step2$FC), priors.step2$FC)
# apollo_beta_new = replace(apollo_beta_new, names(priors.step2$svN), priors.step2$svN)
# apollo_beta = apollo_beta_new[names(apollo_beta)]
# ################################################################# #
#### HB settings ####
# ################################################################# #
apollo_HB = list(
hbDist = c(
# ----------------------------------------------------------------- #
#---- Choice Model
# ----------------------------------------------------------------- #
### Alternative specific constants (ASC)
b_asc_1 = "F",
b_asc_2 = "F",
b_asc_3 = "F",
### Mix_Default
b_Mix_Default = "F",
### Mix_Green
b_Mix_Green = "F",
b_Gender_Mix_Green = "F",
b_Age_Mix_Green = "F",
b_Education_Mix_Green = "F",
### Regio_0
b_Regio_0 = "F",
### Regio_50
b_Regio_50 = "F",
b_Gender_Regio_50 = "F",
b_Age_Regio_50 = "F",
b_Education_Regio_50 = "F",
### Regio_100
b_Regio_100 = "F",
b_Gender_Regio_100 = "F",
b_Age_Regio_100 = "F",
b_Education_Regio_100 = "F",
### Label_TUVSud
b_Label_TUVSud = "F",
### Label_OKPower
b_Label_OKPower = "F",
### Label_GrunerStrom
b_Label_GrunerStrom = "F",
### Label_None
b_Label_None = "F",
### Price
b_Price = "F",
b_Income_Price = "F",
b_PriceMonthly_Price = "F",
# ----------------------------------------------------------------- #
#---- Choice Model DR
# ----------------------------------------------------------------- #
### DR_asc_0
b_DR_asc_0 = "F",
### DR_asc_1
b_DR_asc_1 = "F",
b_Age_DR_asc_1 = "F",
b_Education_DR_asc_1 = "F",
b_Income_DR_asc_1 = "F",
b_PriceMonthly_DR_asc_1 = "F"
),
gNCREP = 25000, # burn-in iterations
gNEREP = 25000, # post burn-in iterations
gINFOSKIP = 100,
gFULLCV = FALSE, # indicates whether a full variance-covariance structure should be used for the random parameters. (Defaults to TRUE)
hIW = TRUE, # boolean indicating if a hierarchical Inverted Wishart should be used when sampling in posterior distribution for the covariance matrix
priorVariance = 2, # The prior variance assumed. Ignored if pvMatrix is not NULL. (Defaults to 2).
# Increasing the prior variance tends to place more weight on fitting each individual's data, and places less emphasis on "borrowing" information from the population parameters.
degreesOfFreedom = 5, # Additional degrees of freedom for the prior variance-covariance matrix, not including the number of parameters. (Defaults to 5)
# The higher the value, the greater the influence of the prior variance and more data are needed to change that prior. The scaling for degrees of freedom is relative to the sample size.
gNSKIP = 1, # Number of iterations in between retaining draws for averaging. (Defaults to 1)
targetAcceptanceFixed = 0.3, # The target acceptance rate in the Metropolis-Hastings algorithm for the fixed parameters. (Defaults to 0.3)
targetAcceptanceNormal = 0.3, # The target acceptance rate in the Metropolis-Hastings algorithm for the random parameters. (Defaults to 0.3)
writeModel = FALSE # Indicates whether the model results should be written as a series of CSV files
# fixedA = rep(0,times=8), # random parameters and their values kept fixed during estimation
# fixedD = rep(1,times=8) # random parameters and their values kept fixed during estimation
)
# ################################################################# #
#### GROUP AND VALIDATE INPUTS ####
# ################################################################# #
apollo_inputs = apollo_validateInputs()
# ################################################################# #
#### DEFINE MODEL AND LIKELIHOOD FUNCTION ####
# ################################################################# #
apollo_probabilities = function(apollo_beta, apollo_inputs, functionality) {
### Function initialisation: do not change the following three commands
### 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()
### Create alternative specific constants and coefficients using interactions with covariates
# ----------------------------------------------------------------- #
#---- Choice Model
# ----------------------------------------------------------------- #
### Alternative specific constants (ASC)
b_asc_1_value = b_asc_1
b_asc_1_value2 = b_asc_1_value
b_asc_2_value = b_asc_2
b_asc_2_value2 = b_asc_2_value
b_asc_3_value = b_asc_3
b_asc_3_value2 = b_asc_3_value
### Mix_Default
b_Mix_Default_value = b_Mix_Default
b_Mix_Default_value2 = b_Mix_Default_value
### Mix_Green
b_Mix_Green_value = b_Mix_Green +
b_Gender_Mix_Green * COV_Gender +
b_Age_Mix_Green * COV_Age_centered +
b_Education_Mix_Green * COV_Education_centered
b_Mix_Green_value2 = b_Mix_Green_value
### Regio_0
b_Regio_0_value = b_Regio_0
b_Regio_0_value2 = b_Regio_0_value
### Regio_50
b_Regio_50_value = b_Regio_50 +
b_Gender_Regio_50 * COV_Gender +
b_Age_Regio_50 * COV_Age_centered +
b_Education_Regio_50 * COV_Education_centered
b_Regio_50_value2 = b_Regio_50_value
### Regio_100
b_Regio_100_value = b_Regio_100 +
b_Gender_Regio_100 * COV_Gender +
b_Age_Regio_100 * COV_Age_centered +
b_Education_Regio_100 * COV_Education_centered
b_Regio_100_value2 = b_Regio_100_value
### Label_TUVSud
b_Label_TUVSud_value = b_Label_TUVSud
b_Label_TUVSud_value2 = b_Label_TUVSud_value
### Label_OKPower
b_Label_OKPower_value = b_Label_OKPower
b_Label_OKPower_value2 = b_Label_OKPower_value
### Label_GrunerStrom
b_Label_GrunerStrom_value = b_Label_GrunerStrom
b_Label_GrunerStrom_value2 = b_Label_GrunerStrom_value
### Label_None
b_Label_None_value = b_Label_None
b_Label_None_value2 = b_Label_None_value
### Price
b_Price_value = b_Price +
b_Income_Price * (1-COV_IncomeQuan_miss) * COV_IncomeQuanClassAggregated_centered +
b_PriceMonthly_Price * (COV_PriceMonthly_centered/100)
b_Price_value2 = b_Price_value
# ----------------------------------------------------------------- #
#---- Choice Model DR
# ----------------------------------------------------------------- #
### DR_asc_0
b_DR_asc_0_value = b_DR_asc_0
b_DR_asc_0_value2 = b_DR_asc_0_value
### DR_asc_1
b_DR_asc_1_value = b_DR_asc_1 +
b_Age_DR_asc_1 * COV_Age_centered +
b_Education_DR_asc_1 * COV_Education_centered +
b_Income_DR_asc_1 * (1-COV_IncomeQuan_miss) * COV_IncomeQuanClassAggregated_centered +
b_PriceMonthly_DR_asc_1 * (COV_PriceMonthly_centered/100)
b_DR_asc_1_value2 = b_DR_asc_1_value
### List of utilities: these must use the same names as in mnl_settings, order is irrelevant
### CBC_forced
V = list()
V[['alt1']] = b_asc_1_value2 +
b_Mix_Default_value2 * 0 +
b_Mix_Green_value2 * Mix2.1 +
b_Regio_0_value2 * 0 +
b_Regio_50_value2 * Regio2.1 +
b_Regio_100_value2 * Regio3.1 +
b_Label_TUVSud_value2 * 0 +
b_Label_OKPower_value2 * Label2.1 +
b_Label_GrunerStrom_value2 * Label3.1 +
b_Label_None_value2 * Label4.1 +
b_Price_value2 * Price.1
V[['alt2']] = b_asc_2_value2 +
b_Mix_Default_value2 * 0 +
b_Mix_Green_value2 * Mix2.2 +
b_Regio_0_value2 * 0 +
b_Regio_50_value2 * Regio2.2 +
b_Regio_100_value2 * Regio3.2 +
b_Label_TUVSud_value2 * 0 +
b_Label_OKPower_value2 * Label2.2 +
b_Label_GrunerStrom_value2 * Label3.2 +
b_Label_None_value2 * Label4.2 +
b_Price_value2 * Price.2
V[['alt3']] = b_asc_3_value2 +
b_Mix_Default_value2 * 0 +
b_Mix_Green_value2 * Mix2.3 +
b_Regio_0_value2 * 0 +
b_Regio_50_value2 * Regio2.3 +
b_Regio_100_value2 * Regio3.3 +
b_Label_TUVSud_value2 * 0 +
b_Label_OKPower_value2 * Label2.3 +
b_Label_GrunerStrom_value2 * Label3.3 +
b_Label_None_value2 * Label4.3 +
b_Price_value2 * Price.3
### Define settings for MNL model component
mnl_settings = list(
alternatives = c(alt1=1, alt2=2, alt3=3),
avail = list(alt1=1, alt2=1, alt3=1),
choiceVar = Choice,
V = V,
rows = (CBC_type=="R")
)
# Compute probabilities using MNL model
P[['CBC_forced']] = apollo_mnl(mnl_settings, functionality)
### CBC_free
V = list()
V[['alt1']] = b_DR_asc_1_value2 +
b_Mix_Default_value2 * 0 +
b_Mix_Green_value2 * Mix2.1 +
b_Regio_0_value2 * 0 +
b_Regio_50_value2 * Regio2.1 +
b_Regio_100_value2 * Regio3.1 +
b_Label_TUVSud_value2 * 0 +
b_Label_OKPower_value2 * Label2.1 +
b_Label_GrunerStrom_value2 * Label3.1 +
b_Label_None_value2 * Label4.1 +
b_Price_value2 * Price.1
V[['alt2']] = b_DR_asc_0_value2
### Define settings for MNL model component
mnl_settings = list(
alternatives = c(alt1=1, alt2=2),
avail = list(alt1=1, alt2=1),
choiceVar = Choice,
V = V,
rows = (CBC_type=="RD")
)
# Compute probabilities using MNL model
P[['CBC_free']] = apollo_mnl(mnl_settings, functionality)
### Combined model
P = apollo_combineModels(P, apollo_inputs, 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)
}
# ################################################################# #
#### MODEL ESTIMATION ####
# ################################################################# #
model = apollo_estimate(apollo_beta, apollo_fixed, apollo_probabilities, apollo_inputs)
# ################################################################# #
#### MODEL OUTPUTS ####
# ################################################################# #
# ----------------------------------------------------------------- #
#---- FORMATTED OUTPUT (TO SCREEN) ----
# ----------------------------------------------------------------- #
apollo_modelOutput(model)
# ----------------------------------------------------------------- #
#---- FORMATTED OUTPUT (TO FILE, using model name) ----
# ----------------------------------------------------------------- #
saveOutput_settings = list(
printClassical = TRUE,
printPVal = TRUE,
printDiagnostics = TRUE,
printCovar = TRUE,
printCorr = TRUE,
printOutliers = 50,
printChange = TRUE,
saveEst = TRUE,
saveCov = TRUE,
saveCorr = TRUE,
saveModelObject = TRUE
)
apollo_saveOutput(model, saveOutput_settings)
# ################################################################# #
#### MISC ####
# ################################################################# #
llCalc.apollo_beta = as.list(apollo_beta)
llCalc.apollo_beta = lapply(llCalc.apollo_beta, rep, times=nrow(database))
apollo_llCalc(apollo_beta=llCalc.apollo_beta,
apollo_probabilities,
apollo_inputs)