inst/examples/bayesian_swissmetro/plot_b02_weight.R

In rbiogeme: Interface to 'Biogeme'

#!/usr/bin/env Rscript

# b02. Bayesian logit estimation with a weighted sample.
#
# The weight expression is symbolic and is passed to native Biogeme as the
# observation weight.  It is not evaluated by R before estimation.

library(rbiogeme)

# prepare_swissmetro_example() is defined in ../swissmetro/example_utils.R;
# it parses --data/--python/--output and prepares a reproducible input session.
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) {
  # Fixed and free parameters retain the native Python names and order.
  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)
  weight <- 0.8890991 * (variable("GROUP") == 2) +
    1.2 * (variable("GROUP") == 3)

  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")
    ),
    weight = weight
  )
}

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

model <- build_model(database)
fit <- bayesian_estimate(
  model,
  model_name = "b02_weight",
  control = biogeme_control(
    output_directory = prepared$output,
    user_notes = paste0(
      "Example of a logit model with three alternatives: Train, Car and ",
      "Swissmetro. Weighted Exogenous Sample Maximum Likelihood estimator (WESML)"
    ),
    mcmc_sampling_strategy = "pymc",
    generate_yaml = TRUE,
    generate_html = TRUE,
    generate_netcdf = TRUE
  )
)
print(summary(fit))
print(coef(fit))
invisible(fit)