ApolloChoiceModelling forum

Code: Select all

# ################################################################# #
#### LOAD LIBRARY AND DEFINE CORE SETTINGS                       ####
# ################################################################# #

### Clear memory
rm(list = ls())

#to check directory (folder where we are working)
getwd()

#Then we load them up 
library(apollo)

### Initialise code
apollo_initialise()

### Set core controls
apollo_control = list(
  modelName       = "MMNL_uncorrelated",
  modelDescr      = "Mixed logit model in WTP space",
  indivID         = "RID",  
  nCores          = 11,
  outputDirectory = "output"
)


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

setwd("//nash/mtec-home/gleduc/My Documents/Postdoc Bordeaux/Choice Experiment/Apollo")
database = read.csv("widedata.csv",header=TRUE)


###Data transformation###

##Rescaling the cost variable to divide it by 10,000##

database[c("COST1", "COST2")] <- database[c("COST1", "COST2")] / 10000


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

### Vector of parameters, including any that are kept fixed in estimation
### starting values = starting values taken from estimated with mixed logit in utility space. See end of code for starting values of WTP with LC that were divided by 1000
##Laure:OECD MXL WTPspace_fixedprice, Date: 2022-12-20 16:42:47 (which used MNL coefs as starting values)

apollo_beta = c(b_asc_alt3            = -12.011,
                sd_asc_alt3           = 3.0157,
                mean_b_suppnone       = 0,
                sd_b_suppnone         = 0,
                mean_b_suppindiv      = 0.5208,
                sd_b_suppindiv        = 0.5557,
                mean_b_suppcoll       = 0.6493,
                sd_b_suppcoll         = 0.8402,
                mean_b_tvxindiv       = 0,
                sd_b_tvxindiv         = 0,
                mean_b_tvxindgp       = -0.0054,
                sd_b_tvxindgp         = 0.6998,
                mean_b_tvxcoop        = -0.1110,
                sd_b_tvxcoop          = 0.4753,
                mean_b_pulvindiv      = 0,
                sd_b_pulvindiv        = 0,
                mean_b_pulvindgp      = -0.6841,
                sd_b_pulvindgp        = 0.9320,
                mean_b_pulvcoop       = -0.9684,
                sd_b_pulvcoop         = 0.7094,
                mean_b_techindiv      = 0, 
                sd_b_techindiv        = 0,
                mean_b_techcoll       = -0.3077,
                sd_b_techcoll         = 0.7801,
                b_cost                = -1.431)
                #sd_log_b_cost         = 0.01)

### 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("mean_b_suppnone", "sd_b_suppnone", "mean_b_tvxindiv", "sd_b_tvxindiv", "mean_b_pulvindiv", "sd_b_pulvindiv", "mean_b_techindiv","sd_b_techindiv")


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

### Set parameters for generating draws

apollo_draws = list(
  interDrawsType = "sobol",
  interNDraws    = 3000,
  interNormDraws = c("draws_suppindiv","draws_suppcoll","draws_tvxindgp","draws_tvxcoop", "draws_pulvindgp", "draws_pulvcoop", "draws_techcoll", "draws_asc_alt3"),
  intraDrawsType = "sobol",
  intraNDraws = 0,
  intraNormDraws = c()
)

### Create random parameters

apollo_randCoeff = function(apollo_beta, apollo_inputs){
  randcoeff = list()
  randcoeff[["asc_alt3"]] = b_asc_alt3 + sd_asc_alt3 * draws_asc_alt3
  randcoeff[["b_suppindiv"]] = mean_b_suppindiv + sd_b_suppindiv * draws_suppindiv
  randcoeff[["b_suppcoll"]] = mean_b_suppcoll + sd_b_suppcoll * draws_suppcoll
  randcoeff[["b_tvxindgp"]] = mean_b_tvxindgp + sd_b_tvxindgp * draws_tvxindgp
  randcoeff[["b_tvxcoop"]] = mean_b_tvxcoop + sd_b_tvxcoop * draws_tvxcoop
  #randcoeff[["b_cost"]] = -exp(mean_log_b_cost + sd_log_b_cost * draws_cost )
  randcoeff[["b_pulvindgp"]] = mean_b_pulvindgp + sd_b_pulvindgp * draws_pulvindgp
  randcoeff[["b_pulvcoop"]] = mean_b_pulvcoop + sd_b_pulvcoop * draws_pulvcoop
  randcoeff[["b_techcoll"]] = mean_b_techcoll + sd_b_techcoll * draws_techcoll
  
  
  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()
  
  ### List of utilities: these must use the same names as in mnl_settings, order is irrelevant
  V = list()
  
  V[["alt1"]] = b_cost*(mean_b_suppnone*(SUPP1==1) + b_suppindiv*(SUPP1==2) + b_suppcoll*(SUPP1==3) + mean_b_tvxindiv*(TVX1==1) + b_tvxindgp*(TVX1==2) + b_tvxcoop*(TVX1==3) + COST1 + mean_b_pulvindiv*(PULV1==1) + b_pulvindgp*(PULV1==2) + b_pulvcoop*(PULV1==3) + mean_b_techindiv*(TECH1==1) + b_techcoll*(TECH1==2))                  
  V[["alt2"]] = b_cost*(mean_b_suppnone*(SUPP2==1) + b_suppindiv*(SUPP2==2) + b_suppcoll*(SUPP2==3) + mean_b_tvxindiv*(TVX2==1) + b_tvxindgp*(TVX2==2) + b_tvxcoop*(TVX2==3) + COST2 + mean_b_pulvindiv*(PULV2==1) + b_pulvindgp*(PULV2==2) + b_pulvcoop*(PULV2==3) + mean_b_techindiv*(TECH2==1) + b_techcoll*(TECH2==2))                     
  V[["alt3"]] = asc_alt3 + mean_b_suppnone*(SUPP3==1)
  
  
  ### Define settings for MNL model component
  mnl_settings = list(
    alternatives  = c("alt1"=1, "alt2"=2, "alt3"=3), 
    avail = 1,
    choiceVar     = prefapollo,
    utilities     = V
  )
  
  ### Compute probabilities using MNL model
  P[['model']] = apollo_mnl(mnl_settings, functionality)
  
  ### Average across intra???individual draws
  #P = apollo_avgIntraDraws(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, estimate_settings=list(maxIterations=300))

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

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

apollo_modelOutput(model, list(printPVal = 2))