tests/testthat/test-indicators-group3.R

In rbiogeme: Interface to 'Biogeme'

test_that("the indicators group 3 examples are self-contained", {
  directory <- rbiogeme_example_path( "indicators")
  files <- file.path(
    directory,
    c(
      "indicator_utils.R",
      "optima.R",
      "plot_b03simulation.R",
      "plot_b04market_shares.R",
      "plot_b05revenues.R"
    )
  )
  expect_true(all(file.exists(files)))
  for (file in files) expect_silent(parse(file))
  expect_true(file.exists(file.path(directory, "optima.dat")))
})

indicator_group3_r_specification <- function(factor = 1.0) {
  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_scenario <- variable("MarginalCostPT") * factor

  v_pt <- asc_pt + beta_time_fulltime * time_pt_scaled * fulltime +
    beta_time_other * time_pt_scaled * not_fulltime +
    beta_cost * (marginal_cost_scenario / 10)
  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 <- list(`0` = v_pt, `1` = v_car, `2` = v_sm)
  nests <- nested_nests(
    choice_set = c(0, 1, 2),
    nests = list(
      nested_nest(mu_no_car, alternatives = c(0, 2), name = "no_car"),
      nested_nest(1, alternatives = 1, name = "car")
    )
  )
  list(
    log_probability = nested_log_probability(
      utilities,
      availability = NULL,
      nests = nests,
      alternative = variable("Choice")
    ),
    v_pt = v_pt,
    v_car = v_car,
    v_sm = v_sm,
    prob_pt = nested_probability(utilities, NULL, nests, 0),
    prob_car = nested_probability(utilities, NULL, nests, 1),
    prob_sm = nested_probability(utilities, NULL, nests, 2),
    marginal_cost_scenario = marginal_cost_scenario
  )
}

native_indicator_group3_specification <- function(factor = 1.0) {
  optima <- reticulate::import("biogeme.data.optima", convert = FALSE)
  expressions <- reticulate::import("biogeme.expressions", convert = FALSE)
  nests_module <- reticulate::import("biogeme.nests", convert = FALSE)
  models <- reticulate::import("biogeme.models", convert = FALSE)

  beta <- expressions$Beta
  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)

  male <- optima$Gender == 1
  female <- optima$Gender == 2
  unreported_gender <- optima$Gender == -1
  fulltime <- optima$OccupStat == 1
  not_fulltime <- optima$OccupStat != 1
  marginal_cost_scenario <- optima$MarginalCostPT * factor
  v_pt <- asc_pt + beta_time_fulltime * (optima$TimePT / 200) * fulltime +
    beta_time_other * (optima$TimePT / 200) * not_fulltime +
    beta_cost * (marginal_cost_scenario / 10)
  v_car <- asc_car + beta_time_fulltime * (optima$TimeCar / 200) * fulltime +
    beta_time_other * (optima$TimeCar / 200) * not_fulltime +
    beta_cost * (optima$CostCarCHF / 10)
  v_sm <- asc_sm + beta_dist_male * (optima$distance_km / 5) * male +
    beta_dist_female * (optima$distance_km / 5) * female +
    beta_dist_unreported * (optima$distance_km / 5) * 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)
  )
  list(
    optima = optima,
    models = models,
    utilities = utilities,
    nests = nests,
    v_pt = v_pt,
    v_car = v_car,
    v_sm = v_sm,
    log_probability = models$lognested(
      utilities,
      NULL,
      nests,
      optima$Choice
    ),
    prob_pt = models$nested(utilities, NULL, nests, 0L),
    prob_car = models$nested(utilities, NULL, nests, 1L),
    prob_sm = models$nested(utilities, NULL, nests, 2L),
    marginal_cost_scenario = marginal_cost_scenario
  )
}

test_that("indicators b03-b05 match native simulation and confidence intervals", {
  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"
  )

  optima_file <- rbiogeme_example_path( "indicators", "optima.dat")
  source(rbiogeme_example_path( "indicators", "indicator_utils.R"))
  source(rbiogeme_example_path( "indicators", "optima.R"))
  data <- read.delim(optima_file, check.names = FALSE, stringsAsFactors = FALSE)
  r_database <- optima_database(data, name = "indicators_group3_r")
  expect_equal(biogeme_database_nrow(r_database), 1899L)

  r_spec <- indicator_group3_r_specification(1.0)
  r_model <- biogeme_model(r_database, formula = r_spec$log_probability)
  r_fit_directory <- tempfile("rbiogeme-indicators-group3-r-")
  native_directory <- tempfile("rbiogeme-indicators-group3-native-")
  dir.create(r_fit_directory, recursive = TRUE)
  dir.create(native_directory, recursive = TRUE)

  original_directory <- getwd()
  setwd(r_fit_directory)
  on.exit(setwd(original_directory), add = TRUE)
  r_fit <- estimate(
    r_model,
    model_name = "b02estimation_group3_r",
    control = biogeme_control(
      model_name = "b02estimation_group3_r",
      bootstrap_samples = 3L,
      generate_html = FALSE,
      generate_yaml = FALSE,
      save_iterations = FALSE
    ),
    run_bootstrap = TRUE
  )
  r_simulation_model <- biogeme_model(
    r_database,
    simulations = list(
      weight = variable("normalized_weight"),
      `Utility PT` = r_spec$v_pt,
      `Utility car` = r_spec$v_car,
      `Utility SM` = r_spec$v_sm,
      `Prob. PT` = r_spec$prob_pt,
      `Prob. car` = r_spec$prob_car,
      `Prob. SM` = r_spec$prob_sm
    )
  )
  r_values <- as.data.frame(simulate(r_simulation_model, beta = r_fit), check.names = FALSE)
  r_bootstrap <- indicator_bootstrap_parameter_draws(r_fit)
  r_intervals <- biogeme_confidence_intervals(
    r_simulation_model,
    beta_values = r_bootstrap,
    interval_size = 0.9
  )

  setwd(native_directory)
  native_spec <- native_indicator_group3_specification(1.0)
  native_database <- native_spec$optima$read_data()
  native_module <- reticulate::import("biogeme.biogeme", convert = FALSE)
  native_fit_object <- native_module$BIOGEME(
    native_database,
    native_spec$log_probability,
    bootstrap_samples = 3L,
    number_of_threads = 1L,
    generate_html = FALSE,
    generate_yaml = FALSE,
    save_iterations = FALSE
  )
  native_fit_object$model_name <- "b02estimation_group3_native"
  native_fit <- native_fit_object$estimate(run_bootstrap = TRUE)
  native_simulator <- native_module$BIOGEME(
    native_database,
    reticulate::dict(
      weight = native_spec$optima$normalized_weight,
      `Utility PT` = native_spec$v_pt,
      `Utility car` = native_spec$v_car,
      `Utility SM` = native_spec$v_sm,
      `Prob. PT` = native_spec$prob_pt,
      `Prob. car` = native_spec$prob_car,
      `Prob. SM` = native_spec$prob_sm
    ),
    generate_html = FALSE,
    generate_yaml = FALSE,
    save_iterations = FALSE
  )
  native_values <- reticulate::py_to_r(native_simulator$simulate(
    the_beta_values = native_fit$get_beta_values()
  ))
  native_beta_values <- reticulate::py_to_r(native_fit$get_beta_values())
  expect_equal(names(r_values), names(native_values))
  expect_equal(
    unname(as.matrix(r_values)),
    unname(as.matrix(native_values)),
    tolerance = 1e-12
  )
  expect_equal(unname(coef(r_fit)), unname(unlist(native_beta_values)), tolerance = 1e-8)
  expect_equal(
    r_fit$final_log_likelihood,
    as.numeric(reticulate::py_to_r(native_fit$final_log_likelihood)),
    tolerance = 1e-8
  )

  native_intervals <- native_simulator$confidence_intervals(
    reticulate::r_to_py(lapply(r_bootstrap, as.list)),
    0.9
  )
  native_left_right <- reticulate::py_to_r(native_intervals)
  expect_equal(
    unname(as.matrix(r_intervals$left)),
    unname(as.matrix(native_left_right[[1L]])),
    tolerance = 1e-12
  )
  expect_equal(
    unname(as.matrix(r_intervals$right)),
    unname(as.matrix(native_left_right[[2L]])),
    tolerance = 1e-12
  )

  r_market_shares <- c(
    PT = mean(r_values$weight * r_values$`Prob. PT`),
    car = mean(r_values$weight * r_values$`Prob. car`),
    SM = mean(r_values$weight * r_values$`Prob. SM`)
  )
  native_market_shares <- c(
    PT = mean(native_values$weight * native_values$`Prob. PT`),
    car = mean(native_values$weight * native_values$`Prob. car`),
    SM = mean(native_values$weight * native_values$`Prob. SM`)
  )
  expect_equal(r_market_shares, native_market_shares, tolerance = 1e-12)

  # Rebuild the public-transportation scenario at factor 1.2 and compare the
  # revenue point estimate and native confidence intervals using identical
  # parameter draws on both sides of the bridge.
  r_revenue_spec <- indicator_group3_r_specification(1.2)
  r_revenue_model <- biogeme_model(
    r_database,
    simulations = list(
      weight = variable("normalized_weight"),
      `Revenue public transportation` =
        r_revenue_spec$prob_pt * r_revenue_spec$marginal_cost_scenario
    )
  )
  r_revenue_values <- as.data.frame(simulate(r_revenue_model, beta = r_fit), check.names = FALSE)
  r_revenue_intervals <- biogeme_confidence_intervals(
    r_revenue_model,
    beta_values = r_bootstrap,
    interval_size = 0.9
  )
  native_revenue_spec <- native_indicator_group3_specification(1.2)
  native_revenue_simulator <- native_module$BIOGEME(
    native_database,
    reticulate::dict(
      weight = native_revenue_spec$optima$normalized_weight,
      `Revenue public transportation` =
        native_revenue_spec$prob_pt * native_revenue_spec$marginal_cost_scenario
    ),
    generate_html = FALSE,
    generate_yaml = FALSE,
    save_iterations = FALSE
  )
  native_revenue_values <- reticulate::py_to_r(native_revenue_simulator$simulate(
    the_beta_values = native_fit$get_beta_values()
  ))
  expect_equal(
    unname(as.matrix(r_revenue_values)),
    unname(as.matrix(native_revenue_values)),
    tolerance = 1e-12
  )
  expect_equal(
    sum(r_revenue_values$weight * r_revenue_values$`Revenue public transportation`),
    sum(native_revenue_values$weight * native_revenue_values$`Revenue public transportation`),
    tolerance = 1e-12
  )
  native_revenue_intervals <- reticulate::py_to_r(native_revenue_simulator$confidence_intervals(
    reticulate::r_to_py(lapply(r_bootstrap, as.list)),
    0.9
  ))
  expect_equal(
    unname(as.matrix(r_revenue_intervals$left)),
    unname(as.matrix(native_revenue_intervals[[1L]])),
    tolerance = 1e-12
  )
  expect_equal(
    unname(as.matrix(r_revenue_intervals$right)),
    unname(as.matrix(native_revenue_intervals[[2L]])),
    tolerance = 1e-12
  )
})