tests/testthat/test-indicators-group5.R

In rbiogeme: Interface to 'Biogeme'

test_that("the indicators group 5 example is self-contained", {
  directory <- rbiogeme_example_path( "indicators")
  files <- file.path(
    directory,
    c(
      "indicator_utils.R",
      "optima.R",
      "plot_b09wtp.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_group5_r_specification <- function() {
  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

  v_pt <- asc_pt + beta_time_fulltime * (time_pt / 200) * fulltime +
    beta_time_other * (time_pt / 200) * not_fulltime +
    beta_cost * (marginal_cost_pt / 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)
  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,
    wtp_pt_time = Derive(v_pt, "TimePT") / Derive(v_pt, "MarginalCostPT"),
    wtp_car_time = Derive(v_car, "TimeCar") / Derive(v_car, "CostCarCHF")
  )
}

native_indicator_group5_specification <- function() {
  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
  v_pt <- asc_pt + beta_time_fulltime * (optima$TimePT / 200) * fulltime +
    beta_time_other * (optima$TimePT / 200) * not_fulltime +
    beta_cost * (optima$MarginalCostPT / 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,
    log_probability = models$lognested(
      util = utilities,
      availability = NULL,
      nests = nests,
      choice = optima$Choice
    ),
    v_pt = v_pt,
    v_car = v_car,
    wtp_pt_time = expressions$Derive(v_pt, "TimePT") /
      expressions$Derive(v_pt, "MarginalCostPT"),
    wtp_car_time = expressions$Derive(v_car, "TimeCar") /
      expressions$Derive(v_car, "CostCarCHF")
  )
}

test_that("indicators b09 WTP and native confidence intervals match", {
  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, "indicator_utils.R"))
  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_group5_r")
  expect_equal(biogeme_database_nrow(r_database), 1899L)

  r_spec <- indicator_group5_r_specification()
  r_model <- biogeme_model(r_database, formula = r_spec$log_probability)
  temporary_directory <- tempfile("rbiogeme-indicators-group5-")
  r_directory <- file.path(temporary_directory, "r")
  native_directory <- file.path(temporary_directory, "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 = "b09wtp_group5_r",
    control = biogeme_control(
      model_name = "b09wtp_group5_r",
      bootstrap_samples = 3L,
      number_of_threads = 1L,
      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"),
      `WTP PT time` = r_spec$wtp_pt_time,
      `WTP CAR time` = r_spec$wtp_car_time
    )
  )
  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_group5_specification()
  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 <- "b09wtp_group5_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,
      `WTP PT time` = native_spec$wtp_pt_time,
      `WTP CAR time` = native_spec$wtp_car_time
    ),
    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_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(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 <- reticulate::py_to_r(native_simulator$confidence_intervals(
    reticulate::r_to_py(lapply(r_bootstrap, as.list)),
    0.9
  ))
  expect_equal(
    unname(as.matrix(r_intervals$left)),
    unname(as.matrix(native_intervals[[1L]])),
    tolerance = 1e-12
  )
  expect_equal(
    unname(as.matrix(r_intervals$right)),
    unname(as.matrix(native_intervals[[2L]])),
    tolerance = 1e-12
  )

  r_average <- mean(60 * r_values$`WTP CAR time` * r_values$weight)
  native_average <- mean(60 * native_values$`WTP CAR time` * native_values$weight)
  expect_equal(r_average, native_average, tolerance = 1e-12)
  expect_equal(
    unique(60 * r_values$`WTP CAR time`),
    unique(60 * native_values$`WTP CAR time`),
    tolerance = 1e-12
  )

  subgroup_wtp <- function(values, left, right, filter) {
    filter <- as.logical(filter)
    size <- sum(filter)
    subgroup_weight <- values$weight[filter]
    subgroup_weight <- subgroup_weight * size / sum(subgroup_weight)
    c(
      value = mean(60 * values$`WTP CAR time`[filter] * subgroup_weight),
      lower = mean(60 * left$`WTP CAR time`[filter] * subgroup_weight),
      upper = mean(60 * right$`WTP CAR time`[filter] * subgroup_weight)
    )
  }
  groups <- list(
    workers = r_database$data$OccupStat == 1,
    females = r_database$data$Gender == 2,
    males = r_database$data$Gender == 1
  )
  for (group_name in names(groups)) {
    r_group <- subgroup_wtp(r_values, r_intervals$left, r_intervals$right, groups[[group_name]])
    native_group <- subgroup_wtp(
      native_values,
      native_intervals[[1L]],
      native_intervals[[2L]],
      groups[[group_name]]
    )
    expect_equal(r_group, native_group, tolerance = 1e-12, info = group_name)
  }
})