R/swissmetro.R

In rbiogeme: Interface to 'Biogeme'

#' Prepare the Swissmetro database used by Biogeme's examples
#'
#' The transformation is deliberately expressed as native database
#' operations. It removes rows with `PURPOSE` other than 1 or 3 and rows with
#' `CHOICE == 0`, then defines the cost, availability, and `/100` scaled
#' variables used by `plot_b01a_logit.py`. Row order and original row IDs are
#' retained. Set `panel = TRUE` to declare `ID` as the panel identifier.
#'
#' @param data Numeric Swissmetro data frame.
#' @param name Database name.
#' @param panel Whether to declare the `ID` panel identifier.
#' @param filter_purpose Whether to retain only observations with `PURPOSE`
#'   equal to 1 or 3, in addition to removing `CHOICE == 0`. Set this to
#'   `FALSE` for native examples whose `read_data()` keeps all purposes.
#' @return A `biogeme_database` specification.
#' @export
swissmetro_data <- function(
    data,
    name = "swissmetro",
    panel = FALSE,
    filter_purpose = TRUE
) {
  required <- c(
    "PURPOSE", "CHOICE", "GA", "SM_CO", "TRAIN_CO", "CAR_AV", "SP",
    "TRAIN_AV", "TRAIN_TT", "SM_TT", "CAR_TT", "CAR_CO", "SM_AV"
  )
  if (!is.data.frame(data) || any(!required %in% names(data))) {
    stop("data must contain the Swissmetro columns: ", paste(required, collapse = ", "), call. = FALSE)
  }
  if (!is.logical(filter_purpose) || length(filter_purpose) != 1L ||
      is.na(filter_purpose)) {
    stop("filter_purpose must be one non-missing logical value.", call. = FALSE)
  }
  database <- biogeme_database(name, data)
  exclude <- if (isTRUE(filter_purpose)) {
    ((variable("PURPOSE") != 1) & (variable("PURPOSE") != 3)) |
      (variable("CHOICE") == 0)
  } else {
    variable("CHOICE") == 0
  }
  database <- biogeme_database_remove(database, exclude)
  derived <- list(
    SM_COST = variable("SM_CO") * (variable("GA") == 0),
    TRAIN_COST = variable("TRAIN_CO") * (variable("GA") == 0),
    CAR_AV_SP = variable("CAR_AV") * (variable("SP") != 0),
    TRAIN_AV_SP = variable("TRAIN_AV") * (variable("SP") != 0),
    TRAIN_TT_SCALED = variable("TRAIN_TT") / 100,
    TRAIN_COST_SCALED = variable("TRAIN_COST") / 100,
    SM_TT_SCALED = variable("SM_TT") / 100,
    SM_COST_SCALED = variable("SM_COST") / 100,
    CAR_TT_SCALED = variable("CAR_TT") / 100,
    CAR_CO_SCALED = variable("CAR_CO") / 100
  )
  for (name_i in names(derived)) {
    database <- biogeme_database_define_variable(database, name_i, derived[[name_i]])
  }
  if (isTRUE(panel)) database <- biogeme_database_panel(database, "ID")
  database
}

#' @rdname swissmetro_data
#' @export
prepare_swissmetro <- swissmetro_data

#' Build the baseline Swissmetro MNL specification
#' @param database A database prepared by [swissmetro_data()].
#' @return A `biogeme_logit_model`.
#' @export
swissmetro_mnl_model <- function(database) {
  validate_biogeme_database(database)
  asc_car <- biogeme_beta("asc_car", start = 0)
  asc_train <- biogeme_beta("asc_train", start = 0)
  asc_sm <- biogeme_beta("asc_sm", start = 0, fixed = TRUE)
  b_time <- biogeme_beta("b_time", start = 0)
  b_cost <- biogeme_beta("b_cost", start = 0)
  model <- logit_model(
    database = database,
    choice = "CHOICE",
    utilities = list(
      `1` = asc_train + b_time * variable("TRAIN_TT_SCALED") + b_cost * variable("TRAIN_COST_SCALED"),
      `2` = asc_sm + b_time * variable("SM_TT_SCALED") + b_cost * variable("SM_COST_SCALED"),
      `3` = asc_car + b_time * variable("CAR_TT_SCALED") + b_cost * variable("CAR_CO_SCALED")
    ),
    availability = list(
      `1` = variable("TRAIN_AV_SP"),
      `2` = variable("SM_AV"),
      `3` = variable("CAR_AV_SP")
    )
  )
}

#' Build the segmented linear-utility Swissmetro specification from b01b
#'
#' This uses native `LinearUtility` and `Segmentation` during bridge
#' compilation.  The generated parameter names are `asc_*_ref` and
#' `asc_*_diff_*`, matching the Python example.
#'
#' @param database A database prepared by [swissmetro_data()].
#' @param user_notes Optional notes stored in the native estimation result.
#' @return A `biogeme_logit_model`.
#' @export
swissmetro_b01b_model <- function(
    database,
    user_notes = paste0(
      "Example of a logit model with three alternatives: Train, Car and ",
      "Swissmetro. Same as 01logit and introducing LinearUtility and ",
      "automatic segmentation of parameters."
    )
) {
  validate_biogeme_database(database)
  male <- biogeme_database_segmentation(
    database,
    "MALE",
    c(`0` = "female", `1` = "male")
  )
  ga <- biogeme_database_segmentation(
    database,
    "GA",
    c(`0` = "without_ga", `1` = "with_ga")
  )
  segmentations <- list(male, ga)

  asc_car <- biogeme_beta("asc_car", start = 0)
  asc_train <- biogeme_beta("asc_train", start = 0)
  b_time <- biogeme_beta("b_time", start = -1.28)
  b_cost <- biogeme_beta("b_cost", start = -1.08)
  segmented_asc_car <- segment_beta(asc_car, segmentations)
  segmented_asc_train <- segment_beta(asc_train, segmentations)

  term <- function(beta, name) linear_term(beta, variable(name))
  model <- logit_model(
    database = database,
    choice = "CHOICE",
    utilities = list(
      `1` = segmented_asc_train + linear_utility(list(
        term(b_time, "TRAIN_TT_SCALED"),
        term(b_cost, "TRAIN_COST_SCALED")
      )),
      `2` = linear_utility(list(
        term(b_time, "SM_TT_SCALED"),
        term(b_cost, "SM_COST_SCALED")
      )),
      `3` = segmented_asc_car + linear_utility(list(
        term(b_time, "CAR_TT_SCALED"),
        term(b_cost, "CAR_CO_SCALED")
      ))
    ),
    availability = list(
      `1` = variable("TRAIN_AV_SP"),
      `2` = variable("SM_AV"),
      `3` = variable("CAR_AV_SP")
    )
  )
  model$control <- biogeme_control(
    user_notes = user_notes,
    save_iterations = FALSE,
    bootstrap_samples = 100,
    second_derivatives = "never"
  )
  model
}