tests/testthat/test-swissmetro-b18a.R

In rbiogeme: Interface to 'Biogeme'

native_swissmetro_b18a <- 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)
  bridge <- rbiogeme:::biogeme_bridge()

  database <- database_module$Database(
    "swissmetro_native_b18a",
    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")
  train_cost <- database$define_variable(
    "TRAIN_COST",
    variable("TRAIN_CO") * (ga == 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
  )

  beta <- expressions$Beta
  b_time <- beta("b_time", 0, NULL, NULL, 0)
  b_cost <- beta("b_cost", 0, NULL, NULL, 0)
  tau1 <- beta("tau1", -1, NULL, 0, 0)
  delta2 <- beta("delta2", 2, 0, NULL, 0)
  tau2 <- tau1 + delta2
  utility <- b_time * train_tt_scaled + b_cost * train_cost_scaled
  log_probability <- expressions$OrderedLogLogit(
    eta = utility,
    cutpoints = reticulate::r_to_py(list(tau1, tau2)),
    y = choice,
    categories = reticulate::r_to_py(c(1, 2, 3)),
    neutral_labels = reticulate::r_to_py(numeric())
  )
  biogeme <- biogeme_module$BIOGEME(
    database,
    log_probability,
    generate_html = FALSE,
    generate_yaml = FALSE,
    save_iterations = FALSE
  )
  biogeme$model_name <- "b18a_ordinal_logit"
  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("b18a Swissmetro ordered logit matches 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)
  b_time <- biogeme_beta("b_time", start = 0)
  b_cost <- biogeme_beta("b_cost", start = 0)
  tau1 <- biogeme_beta("tau1", start = -1, upper = 0)
  delta2 <- biogeme_beta("delta2", start = 2, lower = 0)
  log_probability <- ordered_logit_log_probability(
    eta = b_time * variable("TRAIN_TT_SCALED") +
      b_cost * variable("TRAIN_COST_SCALED"),
    cutpoints = list(tau1, tau1 + delta2),
    alternative = variable("CHOICE"),
    categories = c(1, 2, 3),
    neutral_labels = numeric()
  )
  model <- biogeme_model(
    database = database,
    formula = log_probability
  )

  temporary_directory <- tempfile("rbiogeme-b18a-")
  dir.create(temporary_directory, recursive = TRUE)
  original_directory <- getwd()
  setwd(temporary_directory)
  on.exit(setwd(original_directory), add = TRUE)
  controls <- biogeme_control(
    model_name = "b18a_ordinal_logit",
    generate_html = FALSE,
    generate_yaml = FALSE,
    save_iterations = FALSE
  )
  r_fit <- estimate(
    model,
    model_name = "b18a_ordinal_logit",
    control = controls
  )
  native <- native_swissmetro_b18a(data)
  native_results <- native$results

  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_equal(r_fit$number_of_excluded_data, native_results$number_of_excluded_data)
  expect_identical(isTRUE(r_fit$convergence), isTRUE(native_results$convergence))
})