tests/testthat/test-swissmetro-b01b.R

In rbiogeme: Interface to 'Biogeme'

native_swissmetro_b01b <- function(data, bootstrap_samples = 3L, user_notes) {
  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)
  segmentation_module <- reticulate::import("biogeme.segmentation", convert = FALSE)

  database <- database_module$Database(
    "swissmetro_native_b01b",
    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)

  gender <- database$generate_segmentation(
    variable = variable("MALE"),
    mapping = reticulate::dict(`0` = "female", `1` = "male")
  )
  ga_segmentation <- database$generate_segmentation(
    variable = ga,
    mapping = reticulate::dict(`0` = "without_ga", `1` = "with_ga")
  )
  segmentations <- list(gender, ga_segmentation)

  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", -1.28, NULL, NULL, 0)
  b_cost <- beta("b_cost", -1.08, NULL, NULL, 0)
  segmented_asc_car <- segmentation_module$Segmentation(asc_car, segmentations)$segmented_beta()
  segmented_asc_train <- segmentation_module$Segmentation(asc_train, segmentations)$segmented_beta()
  term <- expressions$LinearTermTuple
  utility <- expressions$LinearUtility
  v_train <- segmented_asc_train + utility(list(
    term(b_time, train_tt_scaled),
    term(b_cost, train_cost_scaled)
  ))
  v_sm <- utility(list(
    term(b_time, sm_tt_scaled),
    term(b_cost, sm_cost_scaled)
  ))
  v_car <- segmented_asc_car + utility(list(
    term(b_time, car_tt_scaled),
    term(b_cost, car_co_scaled)
  ))
  utilities <- reticulate::dict(`1` = v_train, `2` = v_sm, `3` = v_car)
  availability <- reticulate::dict(
    `1` = train_av_sp,
    `2` = variable("SM_AV"),
    `3` = car_av_sp
  )
  biogeme <- biogeme_module$BIOGEME(
    database,
    models$loglogit(utilities, availability, choice),
    user_notes = user_notes,
    save_iterations = FALSE,
    bootstrap_samples = as.integer(bootstrap_samples),
    number_of_threads = 1L,
    calculating_second_derivatives = "never",
    generate_html = FALSE,
    generate_yaml = FALSE
  )
  biogeme$model_name <- "b01b_native"
  biogeme$calculate_null_loglikelihood(availability)
  results <- biogeme$estimate(run_bootstrap = TRUE)
  bridge <- rbiogeme:::biogeme_bridge()
  list(
    results = reticulate::py_to_r(
      bridge$extract_estimation_results(results, variance_covariance_type = "BHHH")
    ),
    number_of_rows = nrow(reticulate::py_to_r(database$dataframe))
  )
}

test_that("b01b segmented Swissmetro 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)
  model <- swissmetro_b01b_model(database)
  expected_names <- c(
    "asc_train_ref", "asc_train_diff_male", "asc_train_diff_with_ga",
    "b_time", "b_cost",
    "asc_car_ref", "asc_car_diff_male", "asc_car_diff_with_ga"
  )
  expect_equal(biogeme_model_parameters(model)$name, expected_names)

  user_notes <- "b01b equivalence test"
  temporary_directory <- tempfile("rbiogeme-b01b-")
  dir.create(temporary_directory, recursive = TRUE)
  original_directory <- getwd()
  setwd(temporary_directory)
  on.exit(setwd(original_directory), add = TRUE)
  control <- biogeme_control(
    second_derivatives = "never",
    variance_covariance_type = "BHHH",
    bootstrap_samples = 3L,
    number_of_threads = 1L,
    seed = 123L,
    user_notes = user_notes,
    generate_html = FALSE,
    generate_yaml = FALSE,
    save_iterations = FALSE
  )
  r_fit <- estimate(
    model,
    model_name = "b01b_r",
    control = control,
    run_bootstrap = TRUE
  )
  native <- native_swissmetro_b01b(data, bootstrap_samples = 3L, user_notes = user_notes)
  native_results <- native$results

  expect_equal(nobs(r_fit), native$number_of_rows)
  expect_identical(r_fit$beta_names, expected_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(r_fit$variance_covariance_type, "BHHH")
  expect_identical(native_results$variance_covariance_type, "BHHH")
  native_vcov <- matrix(
    as.numeric(unlist(native_results$variance_covariance, use.names = FALSE)),
    nrow = length(native_results$beta_names)
  )
  expect_equal(unname(vcov(r_fit)), native_vcov, tolerance = 1e-8)
  expect_identical(r_fit$user_notes, user_notes)
  expect_true(isTRUE(r_fit$bootstrap_complete))
  expect_true(isTRUE(native_results$bootstrap_complete))
  expect_length(r_fit$bootstrap, 3L)
  expect_length(native_results$bootstrap, 3L)
  expect_true(all(vapply(r_fit$bootstrap, length, integer(1)) == length(expected_names)))
})