tests/testthat/test-indicators-group4.R

In rbiogeme: Interface to 'Biogeme'

test_that("the indicators group 4 examples are self-contained", {
  directory <- rbiogeme_example_path( "indicators")
  files <- file.path(
    directory,
    c(
      "indicator_utils.R",
      "optima.R",
      "plot_b06point_elasticities.R",
      "plot_b07cross_elasticities.R",
      "plot_b08arc_elasticities.R"
    )
  )
  expect_true(all(file.exists(files)))
  for (file in files) expect_silent(parse(file))
})

indicator_group4_r_specification <- function(factor = 1.0, nests = NULL) {
  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 <- variable("TimePT")
  time_car <- variable("TimeCar")
  marginal_cost_pt <- variable("MarginalCostPT")
  cost_car <- variable("CostCarCHF")
  distance_km <- variable("distance_km")
  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 <- marginal_cost_pt * factor
  v_pt <- asc_pt + beta_time_fulltime * (time_pt / 200) * fulltime +
    beta_time_other * (time_pt / 200) * not_fulltime +
    beta_cost * (marginal_cost_scenario / 10)
  v_car <- asc_car + beta_time_fulltime * (time_car / 200) * fulltime +
    beta_time_other * (time_car / 200) * not_fulltime +
    beta_cost * (cost_car / 10)
  v_sm <- asc_sm + beta_dist_male * (distance_km / 5) * male +
    beta_dist_female * (distance_km / 5) * female +
    beta_dist_unreported * (distance_km / 5) * unreported_gender
  utilities <- list(`0` = v_pt, `1` = v_car, `2` = v_sm)
  if (is.null(nests)) {
    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")
      )
    )
  }
  prob_pt <- nested_probability(utilities, NULL, nests, 0)
  prob_car <- nested_probability(utilities, NULL, nests, 1)
  prob_sm <- nested_probability(utilities, NULL, nests, 2)
  list(
    utilities = utilities,
    nests = nests,
    v_pt = v_pt,
    v_car = v_car,
    v_sm = v_sm,
    log_probability = nested_log_probability(
      utilities,
      NULL,
      nests,
      variable("Choice")
    ),
    prob_pt = prob_pt,
    prob_car = prob_car,
    prob_sm = prob_sm,
    direct = list(
      pt_time = Derive(prob_pt, "TimePT") * time_pt / prob_pt,
      pt_cost = Derive(prob_pt, "MarginalCostPT") * marginal_cost_pt / prob_pt,
      car_time = Derive(prob_car, "TimeCar") * time_car / prob_car,
      car_cost = Derive(prob_car, "CostCarCHF") * cost_car / prob_car,
      sm_distance = Derive(prob_sm, "distance_km") * distance_km / prob_sm
    ),
    cross = list(
      pt_time = Derive(prob_pt, "TimeCar") * time_car / prob_pt,
      pt_cost = Derive(prob_pt, "CostCarCHF") * cost_car / prob_pt,
      car_time = Derive(prob_car, "TimePT") * time_pt / prob_car,
      car_cost = Derive(prob_car, "MarginalCostPT") * marginal_cost_pt / prob_car
    ),
    marginal_cost_scenario = marginal_cost_scenario
  )
}

native_indicator_group4_specification <- function(factor = 1.0, nests = NULL) {
  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)
  if (is.null(nests)) {
    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)
    )
  }
  prob_pt <- models$nested(utilities, NULL, nests, 0L)
  prob_car <- models$nested(utilities, NULL, nests, 1L)
  prob_sm <- models$nested(utilities, NULL, nests, 2L)
  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 = prob_pt,
    prob_car = prob_car,
    prob_sm = prob_sm,
    direct = list(
      pt_time = expressions$Derive(prob_pt, "TimePT") * optima$TimePT / prob_pt,
      pt_cost = expressions$Derive(prob_pt, "MarginalCostPT") * optima$MarginalCostPT / prob_pt,
      car_time = expressions$Derive(prob_car, "TimeCar") * optima$TimeCar / prob_car,
      car_cost = expressions$Derive(prob_car, "CostCarCHF") * optima$CostCarCHF / prob_car,
      sm_distance = expressions$Derive(prob_sm, "distance_km") * optima$distance_km / prob_sm
    ),
    cross = list(
      pt_time = expressions$Derive(prob_pt, "TimeCar") * optima$TimeCar / prob_pt,
      pt_cost = expressions$Derive(prob_pt, "CostCarCHF") * optima$CostCarCHF / prob_pt,
      car_time = expressions$Derive(prob_car, "TimePT") * optima$TimePT / prob_car,
      car_cost = expressions$Derive(prob_car, "MarginalCostPT") * optima$MarginalCostPT / prob_car
    ),
    marginal_cost_scenario = marginal_cost_scenario
  )
}

test_that("indicators b06-b08 match native elasticity calculations", {
  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"
  )

  directory <- rbiogeme_example_path( "indicators")
  source(file.path(directory, "optima.R"))
  data <- read.delim(
    file.path(directory, "optima.dat"),
    check.names = FALSE,
    stringsAsFactors = FALSE
  )
  r_database <- optima_database(data, name = "indicators_group4_r")
  r_spec <- indicator_group4_r_specification(1.0)
  r_model <- biogeme_model(r_database, formula = r_spec$log_probability)
  r_directory <- tempfile("rbiogeme-indicators-group4-r-")
  native_directory <- tempfile("rbiogeme-indicators-group4-native-")
  dir.create(r_directory, recursive = TRUE)
  dir.create(native_directory, recursive = TRUE)

  original_directory <- getwd()
  setwd(r_directory)
  on.exit(setwd(original_directory), add = TRUE)
  r_fit <- estimate(
    r_model,
    model_name = "b02estimation_group4_r",
    control = biogeme_control(
      model_name = "b02estimation_group4_r",
      generate_html = FALSE,
      generate_yaml = FALSE,
      save_iterations = FALSE
    )
  )
  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. car` = r_spec$prob_car,
      `Prob. public transportation` = r_spec$prob_pt,
      `Prob. slow modes` = r_spec$prob_sm,
      direct_elas_pt_time = r_spec$direct$pt_time,
      direct_elas_pt_cost = r_spec$direct$pt_cost,
      direct_elas_car_time = r_spec$direct$car_time,
      direct_elas_car_cost = r_spec$direct$car_cost,
      direct_elas_sm_dist = r_spec$direct$sm_distance,
      cross_elas_pt_time = r_spec$cross$pt_time,
      cross_elas_pt_cost = r_spec$cross$pt_cost,
      cross_elas_car_time = r_spec$cross$car_time,
      cross_elas_car_cost = r_spec$cross$car_cost
    )
  )
  r_values <- as.data.frame(simulate(r_simulation_model, beta = r_fit), check.names = FALSE)

  setwd(native_directory)
  native_spec <- native_indicator_group4_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,
    generate_html = FALSE,
    generate_yaml = FALSE,
    save_iterations = FALSE
  )
  native_fit_object$model_name <- "b02estimation_group4_native"
  native_fit <- native_fit_object$estimate()
  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. car` = native_spec$prob_car,
      `Prob. public transportation` = native_spec$prob_pt,
      `Prob. slow modes` = native_spec$prob_sm,
      direct_elas_pt_time = native_spec$direct$pt_time,
      direct_elas_pt_cost = native_spec$direct$pt_cost,
      direct_elas_car_time = native_spec$direct$car_time,
      direct_elas_car_cost = native_spec$direct$car_cost,
      direct_elas_sm_dist = native_spec$direct$sm_distance,
      cross_elas_pt_time = native_spec$cross$pt_time,
      cross_elas_pt_cost = native_spec$cross$pt_cost,
      cross_elas_car_time = native_spec$cross$car_time,
      cross_elas_car_cost = native_spec$cross$car_cost
    ),
    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()
  ))
  expect_identical(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(reticulate::py_to_r(native_fit$get_beta_values()), use.names = FALSE)),
    tolerance = 1e-8
  )

  r_direct <- c(
    car_time = sum(r_values$`Prob. car` * r_values$direct_elas_car_time * r_values$weight) /
      sum(r_values$`Prob. car` * r_values$weight),
    car_cost = sum(r_values$`Prob. car` * r_values$direct_elas_car_cost * r_values$weight) /
      sum(r_values$`Prob. car` * r_values$weight),
    pt_time = sum(r_values$`Prob. public transportation` * r_values$direct_elas_pt_time * r_values$weight) /
      sum(r_values$`Prob. public transportation` * r_values$weight),
    pt_cost = sum(r_values$`Prob. public transportation` * r_values$direct_elas_pt_cost * r_values$weight) /
      sum(r_values$`Prob. public transportation` * r_values$weight),
    sm_distance = sum(r_values$`Prob. slow modes` * r_values$direct_elas_sm_dist * r_values$weight) /
      sum(r_values$`Prob. slow modes` * r_values$weight)
  )
  native_direct <- c(
    car_time = sum(native_values$`Prob. car` * native_values$direct_elas_car_time * native_values$weight) /
      sum(native_values$`Prob. car` * native_values$weight),
    car_cost = sum(native_values$`Prob. car` * native_values$direct_elas_car_cost * native_values$weight) /
      sum(native_values$`Prob. car` * native_values$weight),
    pt_time = sum(native_values$`Prob. public transportation` * native_values$direct_elas_pt_time * native_values$weight) /
      sum(native_values$`Prob. public transportation` * native_values$weight),
    pt_cost = sum(native_values$`Prob. public transportation` * native_values$direct_elas_pt_cost * native_values$weight) /
      sum(native_values$`Prob. public transportation` * native_values$weight),
    sm_distance = sum(native_values$`Prob. slow modes` * native_values$direct_elas_sm_dist * native_values$weight) /
      sum(native_values$`Prob. slow modes` * native_values$weight)
  )
  expect_equal(r_direct, native_direct, tolerance = 1e-12)

  # The arc expression is a separate post-estimation scenario using the same
  # base nest structure as the native b08 example.
  r_after <- indicator_group4_r_specification(1.2, nests = r_spec$nests)
  r_arc <- (
    r_after$prob_pt - r_spec$prob_pt
  ) * variable("MarginalCostPT") /
    (
      r_spec$prob_pt *
        (r_after$marginal_cost_scenario - r_spec$marginal_cost_scenario)
    )
  r_arc_model <- biogeme_model(
    r_database,
    simulations = list(
      weight = variable("normalized_weight"),
      `Prob. PT` = r_spec$prob_pt,
      direct_elas_pt = r_arc
    )
  )
  r_arc_values <- as.data.frame(simulate(r_arc_model, beta = r_fit), check.names = FALSE)
  native_after <- native_indicator_group4_specification(1.2, nests = native_spec$nests)
  native_arc_simulator <- native_module$BIOGEME(
    native_database,
    reticulate::dict(
      weight = native_spec$optima$normalized_weight,
      `Prob. PT` = native_spec$prob_pt,
      direct_elas_pt = (
        native_after$prob_pt - native_spec$prob_pt
      ) * native_spec$optima$MarginalCostPT /
        (
          native_spec$prob_pt *
            (native_after$marginal_cost_scenario - native_spec$marginal_cost_scenario)
        )
    ),
    generate_html = FALSE,
    generate_yaml = FALSE,
    save_iterations = FALSE
  )
  native_arc_values <- reticulate::py_to_r(native_arc_simulator$simulate(
    the_beta_values = native_fit$get_beta_values()
  ))
  expect_equal(
    unname(as.matrix(r_arc_values)),
    unname(as.matrix(native_arc_values)),
    tolerance = 1e-12
  )
  expect_equal(
    sum(
      r_arc_values$weight * r_arc_values$`Prob. PT` * r_arc_values$direct_elas_pt /
        sum(r_arc_values$weight * r_arc_values$`Prob. PT`),
      na.rm = TRUE
    ),
    sum(
      native_arc_values$weight * native_arc_values$`Prob. PT` * native_arc_values$direct_elas_pt /
        sum(native_arc_values$weight * native_arc_values$`Prob. PT`),
      na.rm = TRUE
    ),
    tolerance = 1e-12
  )
})