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       = "Hybrid_with_OL_8",
  modelDescr      = "Hybrid choice model on drug choice data, using ordered measurement model for indicators",
  indivID         = "Person_ID",
  mixing          = TRUE,
  nCores          = 22, 
  outputDirectory = "New"
)

# ################################################################# #
#### LOAD DATA AND APPLY ANY TRANSFORMATIONS                     ####
# ################################################################# #
setwd("")

getwd()

### 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)
database = read.csv("[size=150][size=50][/size][/size]",header=TRUE)
### for data dictionary, use ?apollo_drugChoiceData

# ################################################################# #
#### ANALYSIS OF CHOICES                                         ####
# ################################################################# #

### Illustration of how to use apollo_choiceAnalysis with user-defined alternatives.
### This is useful in cases where the alternatives in the data differ 
### across tasks. The same approach can then also be used with unlabelled data

# choiceAnalysis_settings <- list(
#   alternatives = c(alt1=1, alt2=2),
#   avail        = with(database,
#                       list(Artemis     = brand_1=="Artemis"     | brand_2=="Artemis"    , 
#                            Novum       = brand_1=="Novum"       | brand_2=="Novum"      ,
#                            BestValue   = brand_3=="BestValue"   | brand_4=="BestValue"  ,
#                            Supermarket = brand_3=="Supermarket" | brand_4=="Supermarket",
#                            PainAway    = brand_3=="PainAway"    | brand_4=="PainAway"   )
#   ),
#   choiceVar    = with(database,
#                       11*((best==1 & brand_1=="Artemis"    ) | (best==2 & brand_2=="Artemis"    )) + 
#                         12*((best==1 & brand_1=="Novum"      ) | (best==2 & brand_2=="Novum"      )) + 
#                         21*((best==3 & brand_3=="BestValue"  ) | (best==4 & brand_4=="BestValue"  )) + 
#                         22*((best==3 & brand_3=="Supermarket") | (best==4 & brand_4=="Supermarket")) + 
#                         23*((best==3 & brand_3=="PainAway"   ) | (best==4 & brand_4=="PainAway"   ))),
#   explanators  = database[,c("regular_user","university_educated","over_50")]
# )
# 
# apollo_choiceAnalysis(choiceAnalysis_settings, apollo_control, database)

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

### Vector of parameters, including any that are kept fixed in estimation
apollo_beta = c(asc_CV               = 0,
                asc_EV               = 0,
                b_pp                 = 0,
                b_oc                 = 0,
                b_ct                 = 0,
                b_r                  = 0,
                b_cf                 = 0,
                b_e                  = 0,
                b_Male               = 0,
                b_female             = 0,
                b_ageAvg             = 0,
                b_dailyDistanceAvg   = 0,
                b_knowledgeLow       = 0,
                b_knowledgeMiddle    = 0,
                b_knowledgeHigh      = 0,
                b_priorityLanes      = 0,
                b_freeParking        = 0,
                b_tollExemption      = 0,  
                lambda               = 1, 
                sigma_eta            = 3,
                zeta_socIma_1        = 1, 
                zeta_socIma_2        = 1, 
                zeta_socIma_3        = 1, 
                zeta_socIma_4        = 1, 
                tau_socIma_1_1       = 0, 
                tau_socIma_1_2       = 2, 
                tau_socIma_1_3       = 3, 
                tau_socIma_1_4       = 4, 
                tau_socIma_2_1       = 0, 
                tau_socIma_2_2       = 2, 
                tau_socIma_2_3       = 3, 
                tau_socIma_2_4       = 4, 
                tau_socIma_3_1       = 0, 
                tau_socIma_3_2       = 2, 
                tau_socIma_3_3       = 3, 
                tau_socIma_3_4       = 4, 
                tau_socIma_4_1       = 0, 
                tau_socIma_4_2       = 2, 
                tau_socIma_4_3       = 3, 
                tau_socIma_4_4       = 4)

### 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_CV","b_female","b_knowledgeHigh","b_knowledgeLow","b_priorityLanes","tau_socIma_1_1","tau_socIma_2_1","tau_socIma_3_1","tau_socIma_4_1","zeta_socIma_1")

# ################################################################# #
#### DEFINE RANDOM COMPONENTS                                    ####
# ################################################################# #

### Set parameters for generating draws
apollo_draws = list(
  interDrawsType="halton", 
  interNDraws=500,          
  interUnifDraws=c(),      
  interNormDraws=c("eta"), 
  
  intraDrawsType="",
  intraNDraws=0,          
  intraUnifDraws=c(),     
  intraNormDraws=c()      
)

### Create random parameters
apollo_randCoeff=function(apollo_beta, apollo_inputs){
  randcoeff = list()
  
  randcoeff[["LV"]] = b_Male*genderCoded + b_ageAvg*ageAvg + b_dailyDistanceAvg*dailyDistanceAvg + b_knowledgeLow*knowledgeLow + b_knowledgeMiddle*knowledgeMiddle  + b_knowledgeHigh*knowledgeHigh + b_priorityLanes*priorityLanes + b_freeParking*freeParking + b_tollExemption*tollExemption  + sigma_eta*eta
  
  return(randcoeff)
}

# ################################################################# #
#### 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()
  
  ### Likelihood of indicators
  ol_settings1 = list(outcomeOrdered = socIma_1, 
                      V              = zeta_socIma_1*LV, 
                      tau            = list(tau_socIma_1_1, tau_socIma_1_2, tau_socIma_1_3, tau_socIma_1_4),
                      #rows           = (task==1),
                      componentName  = "indic_socIma_1")
  ol_settings2 = list(outcomeOrdered = socIma_2, 
                      V              = zeta_socIma_2*LV, 
                      tau            = list(tau_socIma_2_1, tau_socIma_2_2, tau_socIma_2_3, tau_socIma_2_4), 
                      #rows           = (task==1),
                      componentName  = "indic_socIma_2")
  ol_settings3 = list(outcomeOrdered = socIma_3, 
                      V              = zeta_socIma_3*LV, 
                      tau            = list(tau_socIma_3_1, tau_socIma_3_2, tau_socIma_3_3, tau_socIma_3_4), 
                      #rows           = (task==1),
                      componentName  = "indic_socIma_3")
  ol_settings4 = list(outcomeOrdered = socIma_4, 
                      V              = zeta_socIma_4*LV, 
                      tau            = list(tau_socIma_4_1, tau_socIma_4_2, tau_socIma_4_3, tau_socIma_4_4), 
                      #rows           = (task==1),
                      componentName  = "indic_socIma_4")
  
  P[["indic_socIma_1"]]     = apollo_ol(ol_settings1, functionality)
  P[["indic_socIma_2"]]     = apollo_ol(ol_settings2, functionality)
  P[["indic_socIma_3"]]     = apollo_ol(ol_settings3, functionality)
  P[["indic_socIma_4"]]     = apollo_ol(ol_settings4, functionality)
  
  ### Likelihood of choices
  ### List of utilities: these must use the same names as in mnl_settings, order is irrelevant
  V = list()
  V[["alt1"]] = asc_EV + ( b_pp*(ppEV*incomeAvg*6/10000) + b_oc*ocEV + b_ct*ctEV/60 + b_r*rEV/100 + b_cf*cfEV + b_e*eEV)
  V[["alt2"]] = asc_CV + ( b_pp*(ppCV*incomeAvg*6/10000) + b_oc*ocCV                + b_r*rCV/100             + b_e*eCV + lambda*LV )
  
  
  ### Define settings for MNL model component
  mnl_settings = list(
    alternatives  = c(alt1=1, alt2=2),
    avail         = list(alt1=1, alt2=1),
    choiceVar     = choiceVehicle,
    utilities     = V,
    componentName = "choice"
  )
  
  ### Compute probabilities for MNL model component
  P[["choice"]] = apollo_mnl(mnl_settings, functionality)
  
  ### Likelihood of the whole 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                                            ####
# ################################################################# #

### Optional: calculate LL before model estimation
# apollo_llCalc(apollo_beta, apollo_probabilities, apollo_inputs)

### Estimate model
model = apollo_estimate(apollo_beta, apollo_fixed, apollo_probabilities, apollo_inputs, estimate_settings = list(scaleAfterConvergence=FALSE))

# ################################################################# #
#### MODEL OUTPUTS                                               ####
# ################################################################# #

# ----------------------------------------------------------------- #
#---- FORMATTED OUTPUT (TO SCREEN)                               ----
# ----------------------------------------------------------------- #

apollo_modelOutput(model)

# ----------------------------------------------------------------- #
#---- FORMATTED OUTPUT (TO FILE, using model name)               ----
# ----------------------------------------------------------------- #

# apollo_saveOutput(model)
# 
# # ################################################################# #
# ##### POST-PROCESSING                                            ####
# # ################################################################# #
# 
# ### Print outputs of additional diagnostics to new output file (remember to close file writing when complete)
# apollo_sink()
# 
# # ----------------------------------------------------------------- #
# #---- MODEL PREDICTIONS                                          ----
# # ----------------------------------------------------------------- #
# 
# forecast <- apollo_prediction(model, apollo_probabilities, apollo_inputs,
#                               prediction_settings=list(modelComponent="indic_socIma_1"))
# 
# # ----------------------------------------------------------------- #
# #---- CONDITIONALS AND UNCONDITIONALS                            ----
# # ----------------------------------------------------------------- #
# 
# conditionals <- apollo_conditionals(model,apollo_probabilities,apollo_inputs)
# 
# summary(conditionals)
# 
# unconditionals <- apollo_unconditionals(model,apollo_probabilities,apollo_inputs)
# 
# mean(unconditionals[[1]])
# sd(unconditionals[[1]])
# 
# # ----------------------------------------------------------------- #
# #---- switch off writing to file                                 ----
# # ----------------------------------------------------------------- #

#apollo_sink()