ApolloChoiceModelling forum

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       = "LC_with_covariates",
  modelDescr      = "LC model with covariates in class allocation model on Swiss route choice data",
  indivID         = "id",
  nCores          = 3, 
  outputDirectory = "output"
)

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

### Loading data from package
### if data is to be loaded from a file (e.g. called data.csv), 
### the code would be: database = read.csv("data.csv",header=TRUE)
library(readxl)
database = read_excel("C:/Users/marah/Australian National University/ECA application - Documents/Estimation/Apollo BW/apollobwexcel.xlsx")
#database = read_excel("apollobwexcel.xlsx")
#database = read_excel("apollobwexcel.xlsx")
#View(database)
#attach(database)
### for data dictionary, use ?apollo_drugChoiceData

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

### Vector of parameters, including any that are kept fixed in estimation
apollo_beta=c(bafford_a     = 0,
              bpayback_a      = 0,
              bmargins_a      = 0,
              brentalmarket_a     = 0,
              bsavings_a     = 0,
              brent_a    = 0,
              bresale_a     = 0,
              bnoinstall_a     = 0,
              btime_a       = 0,
              bnoconsider_a  = 0,
              bgov_a       = 0,
              bafford_b     = 0,
              bpayback_b      = 0,
              bmargins_b      = 0,
              brentalmarket_b     = 0,
              bsavings_b     = 0,
              brent_b    = 0,
              bresale_b     = 0,
              bnoinstall_b     = 0,
              btime_b       = 0,
              bnoconsider_b  = 0,
              bgov_b       = 0,
              asc_alt1 = 0,
              asc_alt2 = 0,
              asc_alt3 = 0,
              asc_alt4 = 0,
              asc_alt5 = 0,
              delta_a = 0.03,
              gamma_fullydetached_a = 0,
              gamma_majorurban_a = 0,
              gamma_female_a = 0,
              gamma_yeshavelived_a = 0,
              delta_b = 0,
              gamma_fullydetached_b = 0,
              gamma_majorurban_b = 0,
              gamma_female_b = 0,
              gamma_yeshavelived_b = 0,
              mu_worst       = 1)

### 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("asc_alt5","bnoconsider_a","bnoconsider_b","delta_b","gamma_fullydetached_b","gamma_majorurban_b","gamma_female_b","gamma_yeshavelived_b")


# ################################################################# #
#### LATENT CLASS COMPONENTS                                  ####
# ################################################################# #

apollo_lcPars=function(apollo_beta,apollo_inputs){
  lcpars = list()
  lcpars[["bafford"]] = list(bafford_a, bafford_b)
  lcpars[["bpayback"]] = list(bpayback_a, bpayback_b)
  lcpars[["bmargins"]] = list(bmargins_a, bmargins_b)
  lcpars[["brentalmarket"]] = list(brentalmarket_a, brentalmarket_b)
  lcpars[["bsavings"]] = list(bsavings_a, bsavings_b)
  lcpars[["brent"]] = list(brent_a, brent_b)
  lcpars[["bresale"]] = list(bresale_a, bresale_b)
  lcpars[["bnoinstall"]] = list(bnoinstall_a, bnoinstall_b)
  lcpars[["btime"]] = list(btime_a, btime_b)
  lcpars[["bnoconsider"]] = list(bnoconsider_a, bnoconsider_b)
  lcpars[["bgov"]] = list(bgov_a, bgov_b)
  
  ### Utilities of class allocation model
  V=list()
  V[["class_a"]] = delta_a + gamma_fullydetached_a*fullydetached + gamma_majorurban_a*majorurban + gamma_female_a*female + gamma_yeshavelived_a*yeshavelived
  V[["class_b"]] = delta_b + gamma_fullydetached_b*fullydetached + gamma_majorurban_b*majorurban + gamma_female_b*female + gamma_yeshavelived_b*yeshavelived

  
  mnl_settings = list(
    alternatives = c(class_a=1, class_b=2),
    avail = 1,
    choiceVar = NA,
    V = V
  )
  
lcpars[["pi_values"]]=apollo_mnl(mnl_settings, functionality="raw")

lcpars[["pi_values"]]=apollo_firstRow(lcpars[["pi_values"]],apollo_inputs)
  
  return(lcpars)
}

# ################################################################# #
#### GROUP AND VALIDATE INPUTS                                   ####
# ################################################################# #

apollo_inputs = apollo_validateInputs()

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

apollo_probabilities=function(apollo_beta, apollo_inputs, functionality="estimate"){
  
  ### 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()
  P_bw = list()
  
  ## Define settings for MNL model component that are generic across classes
  mnl_settings = list(
    alternatives = c(alt1=1,alt2=2,alt3=3,alt4=4,alt5=5)
  )
  
  ## Loop over classes
  
  s=1
  while (s<=2){
  
  ### List of utilities: these must use the same names as in mnl_settings, order is irrelevant
  V = list()
  V[["alt1"]] = asc_alt1 + bafford[[s]]*afford1 + bpayback[[s]]*payback1 + bmargins[[s]]*margins1 + brentalmarket[[s]]*rentalmarket1 + bsavings[[s]]*savings1 + brent[[s]]*rent1 + bresale[[s]]*resale1 + bnoinstall[[s]]*noinstall1 + btime[[s]]*time1 + bnoconsider[[s]]*noconsider1 + bgov[[s]]*gov1
  
  V[["alt2"]] = asc_alt2 + bafford[[s]]*afford2 + bpayback[[s]]*payback2 + bmargins[[s]]*margins2 + brentalmarket[[s]]*rentalmarket2 + bsavings[[s]]*savings2 + brent[[s]]*rent2 + bresale[[s]]*resale2 + bnoinstall[[s]]*noinstall2 + btime[[s]]*time2 + bnoconsider[[s]]*noconsider2 + bgov[[s]]*gov2
  
  V[["alt3"]] = asc_alt3 + bafford[[s]]*afford3 + bpayback[[s]]*payback3 + bmargins[[s]]*margins3 + brentalmarket[[s]]*rentalmarket3 + bsavings[[s]]*savings3 + brent[[s]]*rent3 + bresale[[s]]*resale3 + bnoinstall[[s]]*noinstall3 + btime[[s]]*time3 + bnoconsider[[s]]*noconsider3 + bgov[[s]]*gov3
  
  V[["alt4"]] = asc_alt4 + bafford[[s]]*afford4 + bpayback[[s]]*payback4 + bmargins[[s]]*margins4 + brentalmarket[[s]]*rentalmarket4 + bsavings[[s]]*savings4 + brent[[s]]*rent4 + bresale[[s]]*resale4 + bnoinstall[[s]]*noinstall4 + btime[[s]]*time4 + bnoconsider[[s]]*noconsider4 + bgov[[s]]*gov4
  
  V[["alt5"]] = asc_alt5 + bafford[[s]]*afford5 + bpayback[[s]]*payback5 + bmargins[[s]]*margins5 + brentalmarket[[s]]*rentalmarket5 + bsavings[[s]]*savings5 + brent[[s]]*rent5 + bresale[[s]]*resale5 + bnoinstall[[s]]*noinstall5 + btime[[s]]*time5 + bnoconsider[[s]]*noconsider5 + bgov[[s]]*gov5
  
  ### Compute probabilities for "best" choice using MNL model
  mnl_settings$avail = list(alt1=1, alt2=1, alt3=1, alt4=1, alt5=1)
  mnl_settings$choiceVar = best1
  mnl_settings$V = V
  mnl_settings$componentName = paste0("Class_",s)
  
  ## Compute within-class choice probabilities using MNL model
  P_bw[["best"]]=apollo_mnl(mnl_settings,functionality)
  
  ### Take product across observations for same individual
  P_bw[["best"]]=apollo_panelProd(P_bw[["best"]],apollo_inputs,functionality)
  
  ### Compute probabilities for "worst" choice using MNL model
  mnl_settings$avail = list(alt1=(best1!=1), alt2=(best1!=2), alt3=(best1!=3), alt4=(best1!=4), alt5=(best1!=5))
  mnl_settings$choiceVar = worst1
  mnl_settings$V = lapply(V,"*",-mu_worst)
  mnl_settings$componentName = paste0("Class_",s)
  
  ###Compute within-class choice probabilities using MNL model
  P_bw[["worst"]]=apollo_mnl(mnl_settings,functionality)
  
  ##Take product across observations for same individual
  P_bw[["worst"]] = apollo_panelProd(P_bw[["worst"]], apollo_inputs ,functionality)
  
  P[[paste0("Class_",s)]] = apollo_combineModels(P_bw, apollo_inputs, functionality)$model
  
  s=s+1
  }
  
  ###Compute latent class model probabilities
  lc_settings = list(inClassProb = P, classProb = pi_values)
  P[["model"]]=apollo_lc(lc_settings, 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)               ----
  # ----------------------------------------------------------------- #
  
  apollo_saveOutput(model)