inst/doc/apollofirstexample.R

In apollo: Tools for Choice Model Estimation and Application

## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## -----------------------------------------------------------------------------
# ####################################################### #
#### 1. Definition of core settings                        
# ####################################################### #

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

### Load libraries
library(apollo)

### Initialise code
apollo_initialise()

### Set core controls
apollo_control = list(
  modelName  ="MNL",
  modelDescr ="Simple MNL model on mode choice SP data",
  indivID    ="ID"
)

# ####################################################### #
#### 2. Data loading                                   ####
# ####################################################### #

data("apollo_modeChoiceData", package="apollo")
database = apollo_modeChoiceData
rm(apollo_modeChoiceData)

### Use only SP data
database = subset(database,database$SP==1)

# ####################################################### #
#### 3. Parameter definition                           ####
# ####################################################### #

### Vector of parameters, including any that are kept fixed 
### during estimation
apollo_beta=c(asc_car  = 0,
              asc_bus  = 0,
              asc_air  = 0,
              asc_rail = 0,
              b_tt_car = 0,
              b_tt_bus = 0,
              b_tt_air = 0,
              b_tt_rail= 0,
              b_c      = 0)

### Vector with names (in quotes) of parameters to be
###  kept fixed at their starting value in apollo_beta.
### Use apollo_beta_fixed = c() for no fixed parameters.
apollo_fixed = c("asc_car")

# ####################################################### #
#### 4. Input validation                               ####
# ####################################################### #

apollo_inputs = apollo_validateInputs()

# ####################################################### #
#### 5. Likelihood definition                          ####
# ####################################################### #

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[['car']] = asc_car + b_tt_car *time_car + b_c*cost_car
  V[['bus']] = asc_bus + b_tt_bus *time_bus + b_c*cost_bus
  V[['air']] = asc_air + b_tt_air *time_air + b_c*cost_air
  V[['rail']]= asc_rail+ b_tt_rail*time_rail+ b_c*cost_rail
  
  ### Define settings for MNL model component
  mnl_settings = list(
    alternatives  = c(car=1, bus=2, air=3, rail=4), 
    avail         = list(car=av_car, bus=av_bus, 
                         air=av_air, rail=av_rail), 
    choiceVar     = choice,
    V             = V
  )
  
  ### Compute probabilities using MNL model
  P[['model']] = apollo_mnl(mnl_settings, functionality)

  ### Take product across observation for same individual
  P = apollo_panelProd(P, apollo_inputs, functionality)

  ### Prepare and return outputs of function
  P = apollo_prepareProb(P, apollo_inputs, functionality)
  return(P)
}

# ####################################################### #
#### 6. Model estimation and reporting                 ####
# ####################################################### #

model = apollo_estimate(apollo_beta, apollo_fixed, 
                        apollo_probabilities, 
                        apollo_inputs,
                        list(writeIter=FALSE))

apollo_modelOutput(model)

#apollo_saveOutput(model)

# ####################################################### #
#### 7. Postprocessing of results                      ####
# ####################################################### #

### Use the estimated model to make predictions
predictions_base = apollo_prediction(model, 
                                     apollo_probabilities, 
                                     apollo_inputs)

### Now imagine the cost for rail increases by 10% 
### and predict again
database$cost_rail = 1.1*database$cost_rail
apollo_inputs   = apollo_validateInputs()
predictions_new = apollo_prediction(model, 
                                    apollo_probabilities, 
                                    apollo_inputs)

### Compare predictions
change=(predictions_new-predictions_base)/predictions_base
### Not interested in chosen alternative now, 
### so drop last column
change=change[,-ncol(change)]
### Summary of changes (possible presence of NAs due to
### unavailable alternatives)
summary(change)


## ----eval=FALSE---------------------------------------------------------------
#  # ####################################################### #
#  #### 1. Definition of core settings
#  # ####################################################### #
#  
#  ### Clear memory
#  rm(list = ls())
#  
#  ### Load libraries
#  library(apollo)
#  
#  ### Initialise code
#  apollo_initialise()
#  
#  ### Set core controls
#  apollo_control = list(
#    modelName  ="MMNL",
#    modelDescr ="Simple MMNL model on mode choice SP data",
#    indivID    ="ID",
#    mixing     = TRUE,
#    nCores     = 2
#  )
#  
#  # ####################################################### #
#  #### 2. Data loading                                   ####
#  # ####################################################### #
#  
#  data("apollo_modeChoiceData", package="apollo")
#  database = apollo_modeChoiceData
#  rm(apollo_modeChoiceData)
#  
#  ### Use only SP data
#  database = subset(database,database$SP==1)
#  
#  ### Create new variable with average income
#  database$mean_income = mean(database$income)
#  
#  # ####################################################### #
#  #### 3. Parameter definition                           ####
#  # ####################################################### #
#  
#  ### Vector of parameters, including any that are kept fixed
#  ### during estimation
#  apollo_beta=c(asc_car  = 0,
#                asc_bus  =-2,
#                asc_air  =-1,
#                asc_rail =-1,
#                mu_tt    =-4,
#                sigma_tt = 0,
#                b_c      = 0)
#  
#  ### Vector with names (in quotes) of parameters to be
#  ###  kept fixed at their starting value in apollo_beta.
#  ### Use apollo_beta_fixed = c() for no fixed parameters.
#  apollo_fixed = c("asc_car")
#  
#  ### Set parameters for generating draws
#  apollo_draws = list(
#    interDrawsType = "halton",
#    interNDraws    = 500,
#    interUnifDraws = c(),
#    interNormDraws = c("draws_tt")
#  )
#  
#  ### Create random parameters
#  apollo_randCoeff = function(apollo_beta, apollo_inputs){
#    randcoeff = list()
#  
#    randcoeff[["b_tt"]] = -exp(mu_tt + sigma_tt*draws_tt)
#  
#    return(randcoeff)
#  }
#  
#  # ####################################################### #
#  #### 4. Input validation                               ####
#  # ####################################################### #
#  
#  apollo_inputs = apollo_validateInputs()
#  
#  # ####################################################### #
#  #### 5. Likelihood definition                          ####
#  # ####################################################### #
#  
#  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[['car']]  = asc_car  + b_tt*time_car  + b_c*cost_car
#    V[['bus']]  = asc_bus  + b_tt*time_bus  + b_c*cost_bus
#    V[['air']]  = asc_air  + b_tt*time_air  + b_c*cost_air
#    V[['rail']] = asc_rail + b_tt*time_rail + b_c*cost_rail
#  
#    ### Define settings for MNL model component
#    mnl_settings = list(
#      alternatives  = c(car=1, bus=2, air=3, rail=4),
#      avail         = list(car=av_car, bus=av_bus,
#                           air=av_air, rail=av_rail),
#      choiceVar     = choice,
#      V             = V
#    )
#  
#    ### Compute probabilities using MNL model
#    P[['model']] = apollo_mnl(mnl_settings, functionality)
#  
#    ### Take product across observation for same individual
#    P = apollo_panelProd(P, apollo_inputs, functionality)
#  
#    ### Average draws
#    P = apollo_avgInterDraws(P, apollo_inputs, functionality)
#  
#    ### Prepare and return outputs of function
#    P = apollo_prepareProb(P, apollo_inputs, functionality)
#    return(P)
#  }
#  
#  # ####################################################### #
#  #### 6. Model estimation and reporting                 ####
#  # ####################################################### #
#  
#  model = apollo_estimate(apollo_beta, apollo_fixed,
#                          apollo_probabilities,
#                          apollo_inputs,
#                          list(writeIter=FALSE))
#  
#  apollo_modelOutput(model)
#  
#  #apollo_saveOutput(model)
#  
#  # ####################################################### #
#  #### 7. Postprocessing of results                      ####
#  # ####################################################### #
#  
#  ### Use the estimated model to make predictions
#  predictions_base = apollo_prediction(model,
#                                       apollo_probabilities,
#                                       apollo_inputs)
#  
#  ### Now imagine the cost for rail increases by 10%
#  ### and predict again
#  database$cost_rail = 1.1*database$cost_rail
#  apollo_inputs   = apollo_validateInputs()
#  predictions_new = apollo_prediction(model,
#                                      apollo_probabilities,
#                                      apollo_inputs)
#  
#  ### Compare predictions
#  change=(predictions_new-predictions_base)/predictions_base
#  ### Not interested in chosen alternative now,
#  ### so drop last column
#  change=change[,-ncol(change)]
#  ### Summary of changes (possible presence of NAs due to
#  ### unavailable alternatives)
#  summary(change)