tests/testthat/test-swissmetro-b28.R

In rbiogeme: Interface to 'Biogeme'

native_swissmetro_b28 <- function(data) {
  expressions <- reticulate::import("biogeme.expressions", convert = FALSE)
  database_module <- reticulate::import("biogeme.database", convert = FALSE)
  biogeme_module <- reticulate::import("biogeme.biogeme", convert = FALSE)
  models <- reticulate::import("biogeme.models", convert = FALSE)
  bridge <- rbiogeme:::biogeme_bridge()

  database <- database_module$Database("swissmetro_native_b28", reticulate::r_to_py(data))
  variable <- expressions$Variable
  purpose <- variable("PURPOSE")
  choice <- variable("CHOICE")
  database$remove(((purpose != 1) * (purpose != 3) + (choice == 0)) > 0)
  ga <- variable("GA")
  sp <- variable("SP")
  sm_cost <- database$define_variable("SM_COST", variable("SM_CO") * (ga == 0))
  train_cost <- database$define_variable("TRAIN_COST", variable("TRAIN_CO") * (ga == 0))
  car_av_sp <- database$define_variable("CAR_AV_SP", variable("CAR_AV") * (sp != 0))
  train_av_sp <- database$define_variable("TRAIN_AV_SP", variable("TRAIN_AV") * (sp != 0))
  train_tt_scaled <- database$define_variable("TRAIN_TT_SCALED", variable("TRAIN_TT") / 100)
  train_cost_scaled <- database$define_variable("TRAIN_COST_SCALED", train_cost / 100)
  sm_tt_scaled <- database$define_variable("SM_TT_SCALED", variable("SM_TT") / 100)
  sm_cost_scaled <- database$define_variable("SM_COST_SCALED", sm_cost / 100)
  car_tt_scaled <- database$define_variable("CAR_TT_SCALED", variable("CAR_TT") / 100)
  car_co_scaled <- database$define_variable("CAR_CO_SCALED", variable("CAR_CO") / 100)
  beta <- expressions$Beta
  asc_car <- beta("asc_car", 0, NULL, NULL, 0)
  asc_train <- beta("asc_train", 0, NULL, NULL, 0)
  b_time <- beta("b_time", 0, NULL, NULL, 0)
  b_cost <- beta("b_cost", 0, NULL, NULL, 0)
  log_probability <- models$loglogit(
    reticulate::dict(
      `1` = asc_train + b_time * train_tt_scaled + b_cost * train_cost_scaled,
      `2` = b_time * sm_tt_scaled + b_cost * sm_cost_scaled,
      `3` = asc_car + b_time * car_tt_scaled + b_cost * car_co_scaled
    ),
    reticulate::dict(`1` = train_av_sp, `2` = variable("SM_AV"), `3` = car_av_sp),
    choice
  )
  overrides <- expressions$ParameterOverrides()
  overrides$set("b_cost", beta("b_cost", -1, -10, 0, 1))
  overrides$set("asc_train", expressions$Numeric(0))
  log_probability <- expressions$apply_parameter_overrides(log_probability, overrides)
  biogeme <- biogeme_module$BIOGEME(
    database,
    log_probability,
    generate_html = FALSE,
    generate_yaml = FALSE,
    save_iterations = FALSE
  )
  biogeme$model_name <- "b28_parameter_overrides_native"
  results <- biogeme$estimate()
  list(
    results = reticulate::py_to_r(bridge$extract_estimation_results(results)),
    number_of_rows = nrow(reticulate::py_to_r(database$dataframe))
  )
}

test_that("b28 Swissmetro parameter overrides match native Biogeme", {
  skip_if_not(identical(Sys.getenv("RBIOGEME_RUN_INTEGRATION"), "1"), "Set RBIOGEME_RUN_INTEGRATION=1 to run full Swissmetro equivalence tests")
  skip_if_not(rbiogeme_test_configure_python(), "Set RBIOGEME_PYTHON to a compatible native Biogeme interpreter")
  data_path <- rbiogeme_test_swissmetro_path()
  skip_if(!nzchar(data_path), "Set RBIOGEME_SWISSMETRO_DATA to the Swissmetro .dat file")
  data <- read.delim(data_path, check.names = FALSE, stringsAsFactors = FALSE)
  database <- swissmetro_data(data)
  asc_car <- biogeme_beta("asc_car", start = 0)
  asc_train <- biogeme_beta("asc_train", start = 0)
  b_time <- biogeme_beta("b_time", start = 0)
  b_cost <- biogeme_beta("b_cost", start = 0)
  log_probability <- logit_log_probability(
    utilities = list(
      `1` = asc_train + b_time * variable("TRAIN_TT_SCALED") + b_cost * variable("TRAIN_COST_SCALED"),
      `2` = 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")
    ),
    alternative = variable("CHOICE")
  )
  model <- biogeme_model(
    database = database,
    formula = log_probability,
    parameter_overrides = list(
      b_cost = biogeme_beta("b_cost", start = -1, lower = -10, upper = 0, fixed = TRUE),
      asc_train = 0
    )
  )
  temporary_directory <- tempfile("rbiogeme-b28-")
  dir.create(temporary_directory, recursive = TRUE)
  original_directory <- getwd()
  setwd(temporary_directory)
  on.exit(setwd(original_directory), add = TRUE)
  controls <- biogeme_control(
    model_name = "b28_parameter_overrides",
    generate_html = FALSE,
    generate_yaml = FALSE,
    save_iterations = FALSE
  )
  r_fit <- estimate(model, model_name = "b28_parameter_overrides", control = controls)
  native <- native_swissmetro_b28(data)
  expect_equal(nobs(r_fit), native$number_of_rows)
  expect_identical(r_fit$beta_names, native$results$beta_names)
  expect_equal(unname(coef(r_fit)), native$results$beta_values, tolerance = 1e-8)
  expect_equal(as.numeric(logLik(r_fit)), native$results$final_log_likelihood, tolerance = 1e-8)
  expect_false("asc_train" %in% r_fit$beta_names)
  expect_false("b_cost" %in% r_fit$beta_names)
  definitions <- biogeme_model_parameters(model)
  b_cost_definition <- definitions[definitions$name == "b_cost", , drop = FALSE]
  expect_true(isTRUE(b_cost_definition$fixed[[1L]]))
  expect_equal(b_cost_definition$start[[1L]], -1, tolerance = 1e-12)
})