inst/examples/assisted/plot_b02nonlinear.R

In rbiogeme: Interface to 'Biogeme'

#!/usr/bin/env Rscript

# b02nonlinear. Catalog of nonlinear travel-time specifications.
#
# The linear, Box--Cox, and power-series alternatives are symbolic Biogeme
# expressions. Native Python Biogeme evaluates and estimates the three
# synchronized catalog configurations; R does not implement a second
# likelihood or transformation engine.

library(rbiogeme)

# prepare_swissmetro_example() is defined in ../swissmetro/example_utils.R.
# It supplies data and run configuration only. The full model specification is
# kept below, as in the native documentation example.
script_path <- commandArgs(trailingOnly = FALSE)
script_path <- sub("^--file=", "", script_path[startsWith(script_path, "--file=")][[1L]])
source(file.path(dirname(normalizePath(script_path)), "..", "swissmetro", "example_utils.R"))

build_b02nonlinear_model <- function(database) {
  # These parameter names, bounds, and starts match native b02nonlinear.
  asc_car <- biogeme_beta("asc_car", start = 0)
  asc_train <- biogeme_beta("asc_train", start = 0)
  b_time <- biogeme_beta("b_time", start = 0, upper = 0)
  b_cost <- biogeme_beta("b_cost", start = 0, upper = 0)
  lambda_travel_time <- biogeme_beta(
    "lambda_travel_time", start = 1, lower = -10, upper = 10
  )
  square_tt_coef <- biogeme_beta("square_tt_coef", start = 0)
  cube_tt_coef <- biogeme_beta("cube_tt_coef", start = 0)

  # This is the same degree-three power_series() expression as native
  # b02nonlinear, including its exact symbolic multiplication order.
  power_series <- function(the_variable) {
    the_variable + square_tt_coef * the_variable^2 +
      cube_tt_coef * the_variable * the_variable^3
  }

  # A shared controller synchronizes the transformation selected for Train,
  # Swissmetro, and Car.
  time_controller <- catalog_controller(
    "train_tt_catalog",
    c("linear", "boxcox", "power")
  )
  train_tt_catalog <- catalog(
    "train_tt_catalog",
    list(
      linear = variable("TRAIN_TT_SCALED"),
      boxcox = boxcox(variable("TRAIN_TT_SCALED"), lambda_travel_time),
      power = power_series(variable("TRAIN_TT_SCALED"))
    ),
    controller = time_controller
  )
  sm_tt_catalog <- catalog(
    "sm_tt_catalog",
    list(
      linear = variable("SM_TT_SCALED"),
      boxcox = boxcox(variable("SM_TT_SCALED"), lambda_travel_time),
      power = power_series(variable("SM_TT_SCALED"))
    ),
    controller = time_controller
  )
  car_tt_catalog <- catalog(
    "car_tt_catalog",
    list(
      linear = variable("CAR_TT_SCALED"),
      boxcox = boxcox(variable("CAR_TT_SCALED"), lambda_travel_time),
      power = power_series(variable("CAR_TT_SCALED"))
    ),
    controller = time_controller
  )

  utilities <- list(
    `1` = asc_train + b_time * train_tt_catalog +
      b_cost * variable("TRAIN_COST_SCALED"),
    `2` = b_time * sm_tt_catalog + b_cost * variable("SM_COST_SCALED"),
    `3` = asc_car + b_time * car_tt_catalog +
      b_cost * variable("CAR_CO_SCALED")
  )
  log_probability <- logit_log_probability(
    utilities = utilities,
    availability = list(
      `1` = variable("TRAIN_AV_SP"),
      `2` = variable("SM_AV"),
      `3` = variable("CAR_AV_SP")
    ),
    alternative = variable("CHOICE")
  )
  biogeme_model(
    database = database,
    formula = log_probability,
    control = biogeme_control(
    output_directory = prepared$output,
      model_name = "b02nonlinear",
      generate_html = FALSE,
      generate_yaml = FALSE,
      save_iterations = FALSE
    )
  )
}

prepared <- prepare_swissmetro_example(
  commandArgs(trailingOnly = TRUE),
  default_model = "b02nonlinear"
)
database <- swissmetro_data(prepared$data, filter_purpose = FALSE)
model <- build_b02nonlinear_model(database)

fit <- estimate_catalog(
  model,
  model_name = "b02nonlinear",
  control = model$control,
  force = TRUE
)

cat("A total of ", length(fit$results), " models have been estimated.\n", sep = "")
for (configuration in names(fit$results)) {
  result <- fit$results[[configuration]]
  cat(
    configuration,
    ": LL=",
    formatC(result$final_log_likelihood, digits = 2, format = "f"),
    " K=",
    length(result$beta_names),
    "\n",
    sep = ""
  )
}
print(fit$summary)
for (name in names(fit$description)) {
  if (!identical(name, unname(fit$description[[name]]))) {
    cat(name, "\t", fit$description[[name]], "\n", sep = "")
  }
}
cat("Non dominated models:\n")
for (configuration in fit$non_dominated) cat(configuration, "\n", sep = "")
print(fit$non_dominated_summary)
cat(fit$latex, "\n")

invisible(fit)