Estimation of predictions for a joint RP-SP MNL model

[Giansoldati]

Dear All,

I estimated a joint RP-SP multinomial logit model and a joint RP-SP mixed logit model to describe students' preferences for alternative train feeder modes.

I succeeded in the the estimations of both models and in the computation of the predictions_base and the predictions_new following a hypothetical change in the length of cycling lanes, relying on the examples 22 and 3 provided on Apollo's website. Yet, I could calculate the change in the predictions for the RP and the SP components, only separately and not for the whole joint RP-SP MNL (and MXL) models.

I read Apollo's manual, latest version, in different sections andI believe that the issue lies at page 87, where it is stated the following:

The probability for the combined model is obtained by multiplying the RP and the SP components together in P[["model"]], which is the component used for estimation. Rather than doing this manually, we use P = apollo_combineModels(P,apollo_inputs, functionality) as this function also prepares different output depending on the settings of functionality, allowing the use of joint models also for prediction, for example. The multiplication of probabilities across all observations for the same respondent happens after combining the two model components using apollo_panelProd
apollo_combineModels

Moreover, I searched on the R documentation for the

apollo_prediction {apollo}

and noticed that the argument

prediction_settings

can have the following elements:

• *modelComponent Character. Name of component of apollo_probabilities output to calculate predictions for. Default is

  model

  , i.e. the whole model.
  *silent Boolean. If TRUE, this function won't print any output to screen.

Now, I tried to replicate what suggested in the R documentation in my sintax which derives from example 22, with a number of failed attempts, where the following one seemed the most obvious:

Code: Select all

predictions_base = apollo_prediction(model, 
                                     apollo_probabilities, 
                                     apollo_inputs,
                                     prediction_settings=list(modelComponent = "model"))

It does not work and I assume this is the case because if I followed exactly what suggested in the example 22 there is no "model", bust just P, which has no name. Thus I can perform the abovementioned command only for RP and SP model components separately. In that case it worked.

As I said at the beginning of this post I would like to compute the difference between

predictions_new

and

predictions_base

for the whole model and not only for its separate components.

I suppose that I have to abandon the

P = apollo_combineModels(P,apollo_inputs, functionality)

and work manually to get a P[["model"]], but I do not know how to do it.

I would be very thankful if anyone could help me.

Please find attached this post the R file (in txt extension" I used for my estimations which works perfectly fine.

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  ="MXL_student_mobility_RP_e_SP_07_06_prova",
  modelDescr ="Simple MNL model on mode choice RP data",
  indivID    ="ID",
  mixing= TRUE,
  ncores=8
)

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

database = read.csv("Database_06_06_2020.csv", sep=";", header=TRUE)

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


### Create new variable with average income
#database$mean_income = mean(database$income)


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

### Vector of parameters, including any that are kept fixed in estimation
apollo_beta=c(asc_auto_p                      = 0,
              asc_auto_c                      = 0,
              asc_bus                         = 0,
              asc_bse                         = 0,
              asc_piedi                       = 0,
              asc_bici                        = 0,
              asc_scooter                     = 0, 
              #b_fr_bus                       = 0,
              #b_dist                          = 0,
              b_pr_ped                        = 0,
              b_pr_cicl                       = 0,
              b_temp_park                     = 0,
              #b_dist_bus                     = 0,
              b_pend                          = 0,
              mu_b_tempo    =0,
              sigma_b_tempo = 0,
              mu_b_costo = 0,
              sigma_b_costo = 0,
              mu_b_distanza = 0,
              sigma_b_distanza = 0, 
              mu_RP                           = 1,
              mu_SP                           = 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_auto_p","mu_RP")

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

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

### Create random parameters
apollo_randCoeff = function(apollo_beta, apollo_inputs){
  randcoeff = list()
  
  randcoeff[["b_tempo"]] = mu_b_tempo + sigma_b_tempo * draws_tempo
  randcoeff[["b_costo"]] = mu_b_costo + sigma_b_costo * draws_costo
  randcoeff[["b_dist"]] = mu_b_distanza + sigma_b_distanza * draws_distanza
  
  
  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[['auto_p']]   = asc_auto_p   + b_tempo  * Tempo_auto + b_costo * Costo_auto 
  V[['auto_c']]   = asc_auto_c   + b_tempo  * Tempo_auto + b_costo * Costo_auto + b_temp_park * Temp_park
  V[['bus']]      = asc_bus      + b_tempo  * Tempo_autobus + b_costo * Costo_autobus 
  V[['bse']]      = asc_bse      + b_tempo * Tempo_bike_sharing_elettrico  + b_costo * Cents_bike_sharing_elect_minute + b_dist*Distanza +  b_pr_cicl * Perc_ciclabile + b_pend*Pendolare
  V[['piedi']]    = asc_piedi    + b_tempo * Tempo_piedi + b_dist*Distanza    + b_costo * Costo_piedi + b_pr_ped * Perc_pedonale + b_pend*Pendolare
  V[['bici']]     = asc_bici     + b_tempo * Tempo_bici + b_dist*Distanza +  b_pr_cicl * Perc_ciclabile + b_costo * Costo_bici + b_pend*Pendolare 
  V[['scooter']]  = asc_scooter  + b_tempo * Tempo_scooter  + b_costo * Costo_scooter  
  
  
  ### Compute probabilities for the RP part of the data using MNL model
  mnl_settings = list(
    alternatives  = c(auto_p=1, auto_c=7, bus=2, bse=3, piedi=4, bici=5, scooter=6), 
    avail         = list(auto_p=AV_auto_p, auto_c=AV_auto_c, bus=AV_autobus, bse=AV_bse, piedi=AV_piedi, bici=AV_bici, scooter=AV_scooter), 
    choiceVar     = Scelta_2,
    V             = lapply(V, "*", mu_RP),
    rows          = (RP==1)
  )
  
  P[['RP']] = apollo_mnl(mnl_settings, functionality)
  
  ### Compute probabilities for the SP part of the data using MNL model
  mnl_settings$V = lapply(V, "*", mu_SP)
  mnl_settings$rows = (SP==1)
  
  P[['SP']] = apollo_mnl(mnl_settings, functionality)
  
  ### Combined 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                                            ####
# ################################################################# #

model = apollo_estimate(apollo_beta, apollo_fixed, apollo_probabilities, apollo_inputs)

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

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

apollo_modelOutput(model, list(printPVal = TRUE))

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

apollo_saveOutput(model, list(printPVal = TRUE))


# ----------------------------------------------------------------- #
#---- MODEL PREDICTIONS AND ELASTICITY CALCULATIONS              ----
# ----------------------------------------------------------------- #


### RP elasticities

### Use the estimated model to make predictions
predictions_base = apollo_prediction(model, 
                                     apollo_probabilities, 
                                     apollo_inputs,
                                     prediction_settings=list(modelComponent = "RP"))

### Look at a summary of the predicted choice probabilities
summary(predictions_base)

### Now imagine for Perc_ciclabile increases by 50%
database$Perc_ciclabile = 1.50*database$Perc_ciclabile

### Rerun predictions with the new data, and save into a separate matrix
predictions_new = apollo_prediction(model, 
                                    apollo_probabilities, 
                                    apollo_inputs,
                                    prediction_settings=list(modelComponent = "RP"))

### Look at a summary of the predicted choice probabilities
summary(predictions_new)

### Return to original data
database$Perc_ciclabile = 1/1.50*database$Perc_ciclabile


### Compute own elasticity for bike use:
log(sum(predictions_new[,8],na.rm=TRUE)/sum(predictions_base[,8],na.rm=TRUE))/log(1.50)


### SP elasticities

### Use the estimated model to make predictions
predictions_base = apollo_prediction(model, 
                                     apollo_probabilities, 
                                     apollo_inputs,
                                     prediction_settings=list(modelComponent = "SP"))

### Look at a summary of the predicted choice probabilities
summary(predictions_base)

### Now imagine cycling lanes increases by 50%
database$Perc_ciclabile = 1.50*database$Perc_ciclabile

### Rerun predictions with the new data, and save into a separate matrix
predictions_new = apollo_prediction(model, 
                                    apollo_probabilities, 
                                    apollo_inputs,
                                    prediction_settings=list(modelComponent = "SP"))

### Look at a summary of the predicted choice probabilities
summary(predictions_new)

### Return to original data
database$Perc_ciclabile = 1/1.50*database$Perc_ciclabile

### Compute own elasticity for bike:
log(sum(predictions_new[,8],na.rm=TRUE)/sum(predictions_base[,8],na.rm=TRUE))/log(1.50)



########## AUMENTO PERC_PEDONALE


### Use the estimated model to make predictions
predictions_base = apollo_prediction(model, 
                                     apollo_probabilities, 
                                     apollo_inputs,
                                     prediction_settings=list(modelComponent = "RP"))

### Look at a summary of the predicted choice probabilities
summary(predictions_base)

### Now imagine Perc_pedonale increases by 50%
database$Perc_pedonale = 1.50*database$Perc_pedonale

### Rerun predictions with the new data, and save into a separate matrix
predictions_new = apollo_prediction(model, 
                                    apollo_probabilities, 
                                    apollo_inputs,
                                    prediction_settings=list(modelComponent = "RP"))

### Look at a summary of the predicted choice probabilities
summary(predictions_new)

### Return to original data
database$Perc_pedonale = 1/1.50*database$Perc_pedonale


### Compute own elasticity for walk use:
log(sum(predictions_new[,7],na.rm=TRUE)/sum(predictions_base[,7],na.rm=TRUE))/log(1.50)


### SP elasticities

### Use the estimated model to make predictions
predictions_base = apollo_prediction(model, 
                                     apollo_probabilities, 
                                     apollo_inputs,
                                     prediction_settings=list(modelComponent = "SP"))

### Look at a summary of the predicted choice probabilities
summary(predictions_base)

### Now imagine the walking lanes increase by 50%
database$Perc_pedonale = 1.50*database$Perc_pedonale

### Rerun predictions with the new data, and save into a separate matrix
predictions_new = apollo_prediction(model, 
                                    apollo_probabilities, 
                                    apollo_inputs,
                                    prediction_settings=list(modelComponent = "SP"))

### Look at a summary of the predicted choice probabilities
summary(predictions_new)

### Return to original data
database$Perc_pedonale = 1/1.50*database$Perc_pedonale

### Compute own elasticity for walk:
log(sum(predictions_new[,7],na.rm=TRUE)/sum(predictions_base[,7],na.rm=TRUE))/log(1.50)

[stephanehess]

Hi

in estimation, you maximise the log-likelihood of the joint model, i.e. the object P[["model"]] which is given as the product of the RP and SP choices. But in prediction, you obviously predict individual choices, not products of RP and SP choices. That's why in prediction, if your model uses multiple components, you cannot make a prediction from the component "model" but only from the components "RP" and "SP". Predictions look at the probabilities of individual alternatives in individual scenarios, and that's it would not make sense to make a prediction of a joint RP-SP choice. You will of course get probabilities for the RP alternatives and for the SP alternatives, and you could thus calculate the probability of a set of RP&SP choices by multiplying those probabilities together after prediction

Best wishes

Stephane

[Giansoldati]

Dear Prof. Hess,

I do thank you very much for your kind and prompt reply.

Many many thanks.

Best wishes

Marco