tests/testthat/test-indicators-group2.R

In rbiogeme: Interface to 'Biogeme'

test_that("the indicators b02 example and helpers are self-contained", {
  files <- rbiogeme_example_path( "indicators",
    c("indicator_utils.R", "optima.R", "plot_b02estimation.R")
  )
  expect_true(all(file.exists(files)))
  for (file in files) expect_silent(parse(file))
  expect_true(file.exists(rbiogeme_example_path( "indicators", "optima.dat"
  )))
})

test_that("nested probability compiles to native Biogeme", {
  skip_if_not(
    rbiogeme_test_configure_python(),
    "Set RBIOGEME_PYTHON to a compatible native Biogeme interpreter"
  )
  database <- biogeme_database(
    "nested_probability_toy",
    data.frame(Choice = c(0, 1, 2), x = c(1, 2, 3))
  )
  mu <- biogeme_beta("mu", start = 1, lower = 1, upper = 2)
  nests <- nested_nests(
    choice_set = c(0, 1, 2),
    nests = list(
      nested_nest(mu, c(0, 2), name = "no_car"),
      nested_nest(1, 1, name = "car")
    )
  )
  probability <- nested_probability(
    utilities = list(`0` = variable("x"), `1` = 0, `2` = -variable("x")),
    availability = NULL,
    nests = nests,
    alternative = 0
  )
  expect_equal(
    rbiogeme:::collect_biogeme_parameters(probability),
    "mu"
  )
  expect_match(rbiogeme:::format_biogeme_expression(probability), "nested")
  model <- biogeme_model(database = database, simulations = list(probability = probability))
  compiled <- rbiogeme:::biogeme_compile_model(model)
  expect_true(reticulate::py_has_attr(compiled$formulas$probability, "get_value"))
  expect_match(reticulate::py_repr(compiled$formulas$probability), "LogNested")
})

native_indicators_b02 <- 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)
  nests_module <- reticulate::import("biogeme.nests", convert = FALSE)

  filtered <- data[data$Choice != -1, , drop = FALSE]
  filtered <- filtered[!(filtered$Choice == 1 & filtered$CarAvail == 3), , drop = FALSE]
  database <- database_module$Database(
    "native_indicators_b02",
    reticulate::r_to_py(filtered)
  )
  variable <- expressions$Variable
  beta <- expressions$Beta
  time_pt_scaled <- variable("TimePT") / 200
  time_car_scaled <- variable("TimeCar") / 200
  cost_car_scaled <- variable("CostCarCHF") / 10
  distance_scaled <- variable("distance_km") / 5
  male <- variable("Gender") == 1
  female <- variable("Gender") == 2
  unreported_gender <- variable("Gender") == -1
  fulltime <- variable("OccupStat") == 1
  not_fulltime <- variable("OccupStat") != 1
  marginal_cost_pt_scaled <- variable("MarginalCostPT") / 10

  asc_car <- beta("asc_car", 0, NULL, NULL, 0)
  asc_pt <- beta("asc_pt", 0, NULL, NULL, 1)
  asc_sm <- beta("asc_sm", 0, NULL, NULL, 0)
  beta_time_fulltime <- beta("beta_time_fulltime", 0, NULL, NULL, 0)
  beta_time_other <- beta("beta_time_other", 0, NULL, NULL, 0)
  beta_dist_male <- beta("beta_dist_male", 0, NULL, NULL, 0)
  beta_dist_female <- beta("beta_dist_female", 0, NULL, NULL, 0)
  beta_dist_unreported <- beta("beta_dist_unreported", 0, NULL, NULL, 0)
  beta_cost <- beta("beta_cost", 0, NULL, NULL, 0)
  mu_no_car <- beta("mu_no_car", 1, 1, 2, 0)

  v_pt <- asc_pt + beta_time_fulltime * time_pt_scaled * fulltime +
    beta_time_other * time_pt_scaled * not_fulltime +
    beta_cost * marginal_cost_pt_scaled
  v_car <- asc_car + beta_time_fulltime * time_car_scaled * fulltime +
    beta_time_other * time_car_scaled * not_fulltime +
    beta_cost * cost_car_scaled
  v_sm <- asc_sm + beta_dist_male * distance_scaled * male +
    beta_dist_female * distance_scaled * female +
    beta_dist_unreported * distance_scaled * unreported_gender
  utilities <- reticulate::dict(`0` = v_pt, `1` = v_car, `2` = v_sm)
  no_car <- nests_module$OneNestForNestedLogit(
    nest_param = mu_no_car,
    list_of_alternatives = reticulate::r_to_py(list(0L, 2L)),
    name = "no_car"
  )
  car <- nests_module$OneNestForNestedLogit(
    nest_param = 1.0,
    list_of_alternatives = reticulate::r_to_py(list(1L)),
    name = "car"
  )
  nests <- nests_module$NestsForNestedLogit(
    choice_set = reticulate::r_to_py(list(0L, 1L, 2L)),
    tuple_of_nests = reticulate::tuple(no_car, car)
  )
  log_probability <- models$lognested(
    util = utilities,
    availability = NULL,
    nests = nests,
    choice = variable("Choice")
  )
  biogeme <- biogeme_module$BIOGEME(
    database,
    log_probability,
    bootstrap_samples = as.integer(bootstrap_samples),
    number_of_threads = 1L,
    user_notes = user_notes,
    generate_html = FALSE,
    generate_yaml = FALSE,
    save_iterations = FALSE
  )
  biogeme$model_name <- "b02estimation_native"
  results <- biogeme$estimate(run_bootstrap = TRUE)
  bridge <- rbiogeme:::biogeme_bridge()
  list(
    results = reticulate::py_to_r(bridge$extract_estimation_results(results)),
    number_of_rows = nrow(reticulate::py_to_r(database$dataframe))
  )
}

test_that("Optima nested-logit estimation and get_value_c match native Biogeme", {
  skip_if_not(
    identical(Sys.getenv("RBIOGEME_RUN_INTEGRATION"), "1"),
    "Set RBIOGEME_RUN_INTEGRATION=1 to run indicator equivalence tests"
  )
  skip_if_not(
    rbiogeme_test_configure_python(),
    "Set RBIOGEME_PYTHON to a compatible native Biogeme interpreter"
  )
  data_path <- normalizePath(
    rbiogeme_example_path( "indicators", "optima.dat"),
    mustWork = TRUE
  )
  data <- read.delim(data_path, check.names = FALSE, stringsAsFactors = FALSE)
  original_directory <- getwd()
  temporary_directory <- tempfile("rbiogeme-indicators-b02-")
  dir.create(temporary_directory, recursive = TRUE)
  r_directory <- file.path(temporary_directory, "r")
  native_directory <- file.path(temporary_directory, "native")
  dir.create(r_directory)
  dir.create(native_directory)
  on.exit(setwd(original_directory), add = TRUE)

  source(rbiogeme_example_path( "indicators", "optima.R"))
  setwd(r_directory)
  r_database <- read_optima_database(data_path, name = "r_indicators_b02")
  expect_equal(nrow(r_database$data), 1899L)
  expect_equal(biogeme_database_filtered_row_count(r_database), 1899L)
  expect_length(biogeme_database_row_ids(r_database), 1899L)

  asc_car <- biogeme_beta("asc_car", start = 0)
  asc_pt <- biogeme_beta("asc_pt", start = 0, fixed = TRUE)
  asc_sm <- biogeme_beta("asc_sm", start = 0)
  beta_time_fulltime <- biogeme_beta("beta_time_fulltime", start = 0)
  beta_time_other <- biogeme_beta("beta_time_other", start = 0)
  beta_dist_male <- biogeme_beta("beta_dist_male", start = 0)
  beta_dist_female <- biogeme_beta("beta_dist_female", start = 0)
  beta_dist_unreported <- biogeme_beta("beta_dist_unreported", start = 0)
  beta_cost <- biogeme_beta("beta_cost", start = 0)
  mu_no_car <- biogeme_beta("mu_no_car", start = 1, lower = 1, upper = 2)
  time_pt_scaled <- variable("TimePT") / 200
  time_car_scaled <- variable("TimeCar") / 200
  cost_car_scaled <- variable("CostCarCHF") / 10
  distance_scaled <- variable("distance_km") / 5
  male <- variable("Gender") == 1
  female <- variable("Gender") == 2
  unreported_gender <- variable("Gender") == -1
  fulltime <- variable("OccupStat") == 1
  not_fulltime <- variable("OccupStat") != 1
  marginal_cost_pt_scaled <- variable("MarginalCostPT") / 10
  utilities <- list(
    `0` = asc_pt + beta_time_fulltime * time_pt_scaled * fulltime +
      beta_time_other * time_pt_scaled * not_fulltime + beta_cost * marginal_cost_pt_scaled,
    `1` = asc_car + beta_time_fulltime * time_car_scaled * fulltime +
      beta_time_other * time_car_scaled * not_fulltime + beta_cost * cost_car_scaled,
    `2` = asc_sm + beta_dist_male * distance_scaled * male +
      beta_dist_female * distance_scaled * female + beta_dist_unreported * distance_scaled * unreported_gender
  )
  nests <- nested_nests(
    choice_set = c(0, 1, 2),
    nests = list(
      nested_nest(mu_no_car, c(0, 2), name = "no_car"),
      nested_nest(1, 1, name = "car")
    )
  )
  log_probability <- nested_log_probability(
    utilities = utilities,
    availability = NULL,
    nests = nests,
    alternative = variable("Choice")
  )
  model <- biogeme_model(database = r_database, formula = log_probability)
  user_notes <- "b02 indicators equivalence test"
  control <- biogeme_control(
    bootstrap_samples = 3L,
    number_of_threads = 1L,
    user_notes = user_notes,
    generate_html = FALSE,
    generate_yaml = FALSE,
    save_iterations = FALSE
  )
  r_fit <- estimate(
    model,
    model_name = "b02estimation_r",
    control = control,
    run_bootstrap = TRUE
  )
  r_rowwise <- evaluate_biogeme_expression_c(
    model = model,
    expression = log_probability,
    beta = r_fit,
    aggregation = FALSE,
    number_of_draws = 1000L
  )
  r_aggregate <- evaluate_biogeme_expression_c(
    model = model,
    expression = log_probability,
    beta = r_fit,
    aggregation = TRUE,
    number_of_draws = 1000L
  )

  setwd(native_directory)
  native <- native_indicators_b02(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, 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_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)

  native_expressions <- reticulate::import("biogeme.expressions", convert = FALSE)
  native_jax <- reticulate::import("biogeme.jax_calculator", convert = FALSE)
  native_database <- reticulate::import("biogeme.database", convert = FALSE)$Database(
    "native_indicators_b02_values",
    reticulate::r_to_py(data[data$Choice != -1 & !(data$Choice == 1 & data$CarAvail == 3), , drop = FALSE])
  )
  # Rebuild the native expression with the estimated values solely for the
  # post-estimation get_value_c comparison.
  v <- native_expressions$Variable
  nb <- native_expressions$Beta
  n_asc_car <- nb("asc_car", 0, NULL, NULL, 0)
  n_asc_pt <- nb("asc_pt", 0, NULL, NULL, 1)
  n_asc_sm <- nb("asc_sm", 0, NULL, NULL, 0)
  n_btf <- nb("beta_time_fulltime", 0, NULL, NULL, 0)
  n_bto <- nb("beta_time_other", 0, NULL, NULL, 0)
  n_bdm <- nb("beta_dist_male", 0, NULL, NULL, 0)
  n_bdf <- nb("beta_dist_female", 0, NULL, NULL, 0)
  n_bdu <- nb("beta_dist_unreported", 0, NULL, NULL, 0)
  n_bc <- nb("beta_cost", 0, NULL, NULL, 0)
  n_mu <- nb("mu_no_car", 1, 1, 2, 0)
  n_v <- reticulate::dict(
    `0` = n_asc_pt + n_btf * (v("TimePT") / 200) * (v("OccupStat") == 1) +
      n_bto * (v("TimePT") / 200) * (v("OccupStat") != 1) + n_bc * (v("MarginalCostPT") / 10),
    `1` = n_asc_car + n_btf * (v("TimeCar") / 200) * (v("OccupStat") == 1) +
      n_bto * (v("TimeCar") / 200) * (v("OccupStat") != 1) + n_bc * (v("CostCarCHF") / 10),
    `2` = n_asc_sm + n_bdm * (v("distance_km") / 5) * (v("Gender") == 1) +
      n_bdf * (v("distance_km") / 5) * (v("Gender") == 2) + n_bdu * (v("distance_km") / 5) * (v("Gender") == -1)
  )
  n_no_car <- reticulate::import("biogeme.nests", convert = FALSE)$OneNestForNestedLogit(
    nest_param = n_mu,
    list_of_alternatives = reticulate::r_to_py(list(0L, 2L)),
    name = "no_car"
  )
  n_car <- reticulate::import("biogeme.nests", convert = FALSE)$OneNestForNestedLogit(
    nest_param = 1.0,
    list_of_alternatives = reticulate::r_to_py(list(1L)),
    name = "car"
  )
  n_nests <- reticulate::import("biogeme.nests", convert = FALSE)$NestsForNestedLogit(
    choice_set = reticulate::r_to_py(list(0L, 1L, 2L)),
    tuple_of_nests = reticulate::tuple(n_no_car, n_car)
  )
  n_log_probability <- reticulate::import("biogeme.models", convert = FALSE)$lognested(
    n_v, NULL, n_nests, v("Choice")
  )
  native_beta_values <- as.list(coef(r_fit))
  native_rowwise <- reticulate::py_to_r(native_jax$get_value_c(
    expression = n_log_probability,
    betas = reticulate::r_to_py(native_beta_values),
    database = native_database,
    numerically_safe = FALSE,
    use_jit = TRUE,
    aggregation = FALSE
  ))
  native_aggregate <- reticulate::py_to_r(native_jax$get_value_c(
    expression = n_log_probability,
    betas = reticulate::r_to_py(native_beta_values),
    database = native_database,
    numerically_safe = FALSE,
    use_jit = TRUE,
    aggregation = TRUE
  ))
  expect_equal(r_rowwise, as.numeric(native_rowwise), tolerance = 1e-8)
  expect_equal(r_aggregate, as.numeric(native_aggregate), tolerance = 1e-8)
  expect_equal(r_aggregate, as.numeric(logLik(r_fit)), tolerance = 1e-8)
})