Important: Read this before posting to this forum

  1. This forum is for questions related to the use of Apollo. We will answer some general choice modelling questions too, where appropriate, and time permitting. We cannot answer questions about how to estimate choice models with other software packages.
  2. There is a very detailed manual for Apollo available at http://www.ApolloChoiceModelling.com/manual.html. This contains detailed descriptions of the various Apollo functions, and numerous examples are available at http://www.ApolloChoiceModelling.com/examples.html. In addition, help files are available for all functions, using e.g. ?apollo_mnl
  3. Before asking a question on the forum, users are kindly requested to follow these steps:
    1. Check that the same issue has not already been addressed in the forum - there is a search tool.
    2. Ensure that the correct syntax has been used. For any function, detailed instructions are available directly in Apollo, e.g. by using ?apollo_mnl for apollo_mnl
    3. Check the frequently asked questions section on the Apollo website, which discusses some common issues/failures. Please see http://www.apollochoicemodelling.com/faq.html
    4. Make sure that R is using the latest official release of Apollo.
  4. If the above steps do not resolve the issue, then users should follow these steps when posting a question:
    1. provide full details on the issue, including the entire code and output, including any error messages
    2. posts will not immediately appear on the forum, but will be checked by a moderator first. This may take a day or two at busy times. There is no need to submit the post multiple times.

Using apollo_llCalc in Bayesian estimation

Report bugs or highlight issues with Apollo functions. At a minimum, please include the part of the output where you believe the bug is manifested. Ideally, please share your model file.
Post Reply
cybey
Posts: 60
Joined: 26 Apr 2020, 19:38

Using apollo_llCalc in Bayesian estimation

Post by cybey »

Helle Stephane, hello David,

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
A look in the manual revealed that apollo_beta must have as many columns as observations in the database.
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).
So I tried this:

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)
But I still get the same error message. When looking into the source code of apollo_llCalc I saw an if-statement that is always false:

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)){
Pout is not a list when apollo_probabilities is called:

Code: Select all

apollo_probabilities(apollo_beta=llCalc.apollo_beta, apollo_inputs, functionality = "output")
It's a vector.

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)
stephanehess
Site Admin
Posts: 974
Joined: 24 Apr 2020, 16:29

Re: Using apollo_llCalc in Bayesian estimation

Post by stephanehess »

Hi

there was indeed a bug which will be fixed in 0.2.8, so that you would now get a result instead of NA. However, the use of this function doesn't make particular sense with HB estimation, as we will explain in the manual for 0.2.8, as follows:

For classical estimation, no changes are needed, and even with random coefficients, the full distribution with numerical simulation of the log-likelihood would then be used. However, for Bayesian estimation, apollo_beta would need to be passed to apollo\_llCalc as a list, with one element per parameter, where each of these has as many values as there are observations in the data. Of course, this would imply that the log-likelihood is calculated at a point value for each parameter and observation, not recognising the random distributions.

Stephane & David
--------------------------------
Stephane Hess
www.stephanehess.me.uk
cybey
Posts: 60
Joined: 26 Apr 2020, 19:38

Re: Using apollo_llCalc in Bayesian estimation

Post by cybey »

Hi guys,

Thanks for the quick reply. I am aware that with Bayesian estimation the calculation of the log-likelihood function is only point-wise. But this is exactly what I was looking for: If I understood it correctly, to calculate the WAIC, I have to sample from the posteriori distribution and calculate the log-likelihood value with these samples. By the way, this would also be a cool feature for future versions of Apollo: AIC, BIC, DIC, and WAIC for Bayesian estimation.

Nico
stephanehess
Site Admin
Posts: 974
Joined: 24 Apr 2020, 16:29

Re: Using apollo_llCalc in Bayesian estimation

Post by stephanehess »

Nico

by the way, for now, you can already use llcalc by just setting HB=FALSE in apollo_control and still following the guidance below on dimensionality

We will consider adding further functionalities

Stephane
--------------------------------
Stephane Hess
www.stephanehess.me.uk
Post Reply