ApolloChoiceModelling forum

Code: Select all

#working directory
setwd("C:\\Users\\calde\\Google Drive\\Mestrado 2019_Gabriel Caldeira\\Modelos\\OD\\Banco de dados")

library(apollo)
library(tidyverse)

rm(list = ls())

parallel::detectCores() #cores available

#read database
database <- read.csv(file = paste0(getwd(),'\\cleaning\\database_calib.csv'))

database <- database %>% 
  mutate(picomanha = case_when(h_saida >= 6 & h_saida <= 9 ~1,TRUE~0))

###Preparing environment
apollo_initialise()

###Options controlling the running of the code
apollo_control = list(
  modelName  ="mnl_lc18",
  modelDescr = "latent class model",
  indivID    ="id_pess",
  mixing     = T,
  nCores     = 10
)

#parameters
apollo_beta <- c(asc1 = 0,
                 asc2 = 0,
                 asc3 = 0,
                 #class 1
                 b_ttauto_1 = 0,
                 b_ttbus_1 = 0,
                 b_ttmetro_1 = 0,
                 b_coauto_1 = 0,
                 b_cometro_1 = 0,
                 #class 2
                 b_ttauto_2 = 0,
                 b_ttbus_2 = 0,
                 b_ttmetro_2 = 0,
                 b_coauto_2 = 0,
                 b_cometro_2 = 0,
                 #class membership
                 s_1 = 0,
                 s_2 = 0,
                 b_sexo_al = 0,
                 b_idade_al = 0,
                 b_nomorad_al = 0,
                 b_totviag_al = 0,
                 b_grau1_al = 0,
                 b_grau2_al = 0,
                 b_pico_al = 0,
                 param = 0)

#fixed parameters
apollo_fixed <- c('asc1','s_1')

#starting values
setwd('C:\\Users\\calde\\Google Drive\\Mestrado 2019_Gabriel Caldeira\\JOCM\\modelo')

apollo_beta <- apollo_readBeta(apollo_beta,apollo_fixed,'mnl_lc10',overwriteFixed = FALSE)

#criando draws
apollo_draws = list(
  interDrawsType = 'halton',
  interNDraws = 300,
  interNormDraws = c('eta')
)

apollo_randCoeff <- function(apollo_beta,apollo_inputs){
  randcoeff = list()
  
  randcoeff[['LV']] = b_sexo_al*sexod + b_idade_al*idade + b_nomorad_al*no_morad + b_totviag_al*tot_viag +
    b_grau1_al*grau1 + b_grau2_al*grau2 + b_pico_al*picomanha + eta
  return(randcoeff)
}

###Grouping latent class parameters
apollo_lcPars = function(apollo_beta, apollo_inputs){
  lcpars = list()
  
  lcpars[['b_ttauto']] = list(b_ttauto_1,b_ttauto_2)
  
  lcpars[['b_ttbus']] = list(b_ttbus_1,b_ttbus_2)
  
  lcpars[['b_ttmetro']] = list(b_ttmetro_1,b_ttmetro_2)
  
  lcpars[['b_coauto']] = list(b_coauto_1,b_coauto_2)
  
  lcpars[['b_cometro']] = list(b_cometro_1,b_cometro_2)
  
  ###Class membership probabilities based on s_1, s_2: use of MNL fomrula
  V=list()
  
  V[["class_1"]] = s_1
  V[["class_2"]] = s_2 + param*LV
  
  ###Settings for class membership probabilities 
  mnl_settings = list(
    alternatives = c(class_1=1, class_2=2), 
    avail        = 1, 
    choiceVar    = NA, 
    V            = V
  )
  ###Class membership probabilities
  lcpars[["pi_values"]] = apollo_mnl(mnl_settings, functionality="raw")
  return(lcpars)
}

###busca pelos dados de entrada necessarios 
apollo_inputs = apollo_validateInputs()

#cosntruindo a funcao de verossimilhanca
apollo_probabilities <-  function(apollo_beta,apollo_inputs,functionality = 'estimate'){
  
  ###Attaches parameters and data so that variables can be referred by name
  apollo_attach(apollo_beta, apollo_inputs)
  on.exit(apollo_detach(apollo_beta, apollo_inputs))
  
  ###Create list for probabilities
  P = list()

  ### Define settings for MNL model component
  mnl_settings = list(
    alternatives = c(auto=3, bus=2, metro=1),
    avail        = list(auto=1, bus=1, metro=1),
    choiceVar    = choice
   )
  
  ### Compute class-specific utilities
  V=list()
  
  V[['auto']]  = asc1 + b_ttauto[[1]]*tt_auto + b_coauto[[1]]*co_auto 
  
  V[['bus']]  = asc2 + b_ttbus[[1]]*tt_bus_otp 
  
  V[['metro']]  = asc3 + b_ttmetro[[1]]*tt_metro + b_cometro[[1]]*co_metro 
  
  ###Calculating probabilities based on MNL function for class 1
  mnl_settings$V = V
  P[[1]] = apollo_mnl(mnl_settings, functionality)
  
  ### Compute class-specific regrets
  R=list()
  
  R[['auto']]  = asc1 - log(1 + exp(b_ttbus[[2]]*tt_bus_otp-b_ttauto[[2]]*tt_auto)) -
    log(1 + exp(b_ttmetro[[2]]*tt_metro-b_ttauto[[2]]*tt_auto)) - log(1 + exp(b_cometro[[2]]*co_metro-b_coauto[[2]]*co_auto)) 
  
  R[['bus']]  = asc2 - log(1 + exp(b_ttauto[[2]]*tt_auto-b_ttbus[[2]]*tt_bus_otp)) - 
    log(1 + exp(b_ttmetro[[2]]*tt_metro-b_ttbus[[2]]*tt_bus_otp))
  
  R[['metro']]  = asc3 - log(1 + exp(b_ttauto[[2]]*tt_auto-b_ttmetro[[2]]*tt_metro)) - 
    log(1 + exp(b_ttbus[[2]]*tt_bus_otp-b_ttmetro[[2]]*tt_metro)) - log(1 + exp(b_coauto[[2]]*co_auto-b_cometro[[2]]*co_metro)) 
  
  ###Calculating probabilities based on MNL function for class 2
  mnl_settings$V = R 
  P[[2]] = apollo_mnl(mnl_settings, functionality)
  
  ###Calculating choice probabilities using class membership and conditional probabilities 
  lc_settings   = list(inClassProb = P, classProb=pi_values)
  P[["model"]] = apollo_lc(lc_settings, apollo_inputs, functionality)
  
  P = apollo_avgInterDraws(P,apollo_inputs,functionality)
  
  ### Prepare and return outputs of function
  P = apollo_prepareProb(P, apollo_inputs, functionality)
  
  return(P)
}