inst/examples/bayesian_swissmetro/plot_b03_scale.R

In rbiogeme: Interface to 'Biogeme'

#!/usr/bin/env Rscript

# b03. Bayesian moneymetric and heteroscedastic specification.
#
# The cost coefficient is fixed at -1, so the utility scale is estimated. The
# group-specific scale expression is compiled into native Biogeme and the
# positive lower bound is enforced by native parameter handling.

library(rbiogeme)

# prepare_swissmetro_example() is defined in ../swissmetro/example_utils.R.
# It makes the script independent of the user's current working directory.
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_model <- function(database) {
  positive_lower_bound <- 1e-5
  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, upper = 0)
  b_cost <- biogeme_beta("b_cost", start = -1, fixed = TRUE)
  scale_not_group3 <- biogeme_beta(
    "scale_not_group3", start = 1, lower = positive_lower_bound
  )
  scale_group3 <- biogeme_beta(
    "scale_group3", start = 1, lower = positive_lower_bound
  )

  scale <- (variable("GROUP") != 3) * scale_not_group3 +
    (variable("GROUP") == 3) * scale_group3
  train <- asc_train + b_time * variable("TRAIN_TT_SCALED") +
    b_cost * variable("TRAIN_COST_SCALED")
  swissmetro <- asc_sm + b_time * variable("SM_TT_SCALED") +
    b_cost * variable("SM_COST_SCALED")
  car <- asc_car + b_time * variable("CAR_TT_SCALED") +
    b_cost * variable("CAR_CO_SCALED")

  logit_model(
    database = database,
    choice = "CHOICE",
    utilities = list(`1` = scale * train, `2` = scale * swissmetro, `3` = scale * car),
    availability = list(
      `1` = variable("TRAIN_AV_SP"),
      `2` = variable("SM_AV"),
      `3` = variable("CAR_AV_SP")
    )
  )
}

prepared <- prepare_swissmetro_example(
  commandArgs(trailingOnly = TRUE),
  default_model = "b03_scale"
)
database <- swissmetro_data(prepared$data)
unlink(file.path(prepared$output, c("b03_scale.yaml", "b03_scale.nc", "b03_scale.html")))

model <- build_model(database)
fit <- bayesian_estimate(
  model,
  model_name = "b03_scale",
  control = biogeme_control(
    output_directory = prepared$output,
    user_notes = paste0(
      "Illustrates a moneymetric heteroscedastic specification. A different scale is",
      " associated with different segments of the sample."
    ),
    generate_yaml = TRUE,
    generate_html = TRUE,
    generate_netcdf = TRUE
  )
)
print(summary(fit))
print(coef(fit))
invisible(fit)