test-swissmetro-b15a.R
In rbiogeme: Interface to 'Biogeme'

native_swissmetro_b15a <- function(data, number_of_draws = 128L, seed = 1223L) {
  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)
  bridge <- rbiogeme:::biogeme_bridge()

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

  beta <- expressions$Beta
  classes <- 0:1
  b_cost <- lapply(classes, function(class) beta(paste0("b_cost_class", class), 0, NULL, NULL, 0))
  b_time <- lapply(classes, function(class) beta(paste0("b_time_class", class), 0, NULL, NULL, 0))
  b_time_s <- lapply(classes, function(class) beta(paste0("b_time_s_class", class), 1, NULL, NULL, 0))
  b_time_rnd <- Map(
    function(class, location, scale) location + scale * expressions$Draws(
      paste0("b_time_rnd_class", class), "NORMAL_ANTI"
    ),
    classes,
    b_time,
    b_time_s
  )

  asc_car <- lapply(classes, function(class) beta(paste0("asc_car_class", class), 0, NULL, NULL, 0))
  asc_car_s <- lapply(classes, function(class) beta(paste0("asc_car_s_class", class), 1, NULL, NULL, 0))
  asc_car_rnd <- Map(
    function(class, location, scale) location + scale * expressions$Draws(
      paste0("asc_car_rnd_class", class), "NORMAL_ANTI"
    ),
    classes,
    asc_car,
    asc_car_s
  )

  asc_train <- lapply(classes, function(class) beta(paste0("asc_train_class", class), 0, NULL, NULL, 0))
  asc_train_s <- lapply(classes, function(class) beta(paste0("asc_train_s_class", class), 1, NULL, NULL, 0))
  asc_train_rnd <- Map(
    function(class, location, scale) location + scale * expressions$Draws(
      paste0("asc_train_rnd_class", class), "NORMAL_ANTI"
    ),
    classes,
    asc_train,
    asc_train_s
  )

  asc_sm <- lapply(classes, function(class) beta(paste0("asc_sm_class", class), 0, NULL, NULL, 1))
  asc_sm_s <- lapply(classes, function(class) beta(paste0("asc_sm_s_class", class), 1, NULL, NULL, 0))
  asc_sm_rnd <- Map(
    function(class, location, scale) location + scale * expressions$Draws(
      paste0("asc_sm_rnd_class", class), "NORMAL_ANTI"
    ),
    classes,
    asc_sm,
    asc_sm_s
  )
  b_time_rnd[[1L]] <- 0

  score_class_0 <- beta("score_class_0", -1.7, NULL, NULL, 0)
  probability_class_0 <- models$logit(
    reticulate::dict(`0` = score_class_0, `1` = 0),
    NULL,
    0L
  )
  probability_class_1 <- models$logit(
    reticulate::dict(`0` = score_class_0, `1` = 0),
    NULL,
    1L
  )

  utilities <- lapply(seq_along(classes), function(index) {
    reticulate::dict(
      `1` = asc_train_rnd[[index]] + b_time_rnd[[index]] * train_tt_scaled +
        b_cost[[index]] * train_cost_scaled,
      `2` = asc_sm_rnd[[index]] + b_time_rnd[[index]] * sm_tt_scaled +
        b_cost[[index]] * sm_cost_scaled,
      `3` = asc_car_rnd[[index]] + b_time_rnd[[index]] * car_tt_scaled +
        b_cost[[index]] * car_co_scaled
    )
  })
  availability <- reticulate::dict(
    `1` = train_av_sp,
    `2` = variable("SM_AV"),
    `3` = car_av_sp
  )
  trajectories <- lapply(utilities, function(class_utilities) {
    expressions$PanelLikelihoodTrajectory(models$logit(
      class_utilities,
      availability,
      choice
    ))
  })
  choice_probability <- probability_class_0 * trajectories[[1L]] +
    probability_class_1 * trajectories[[2L]]
  log_probability <- expressions$log(expressions$MonteCarlo(choice_probability))

  parameter_groups <- reticulate::r_to_py(list(
    `Class 0` = c(
      "b_cost_class0", "asc_car_class0", "asc_car_s_class0",
      "asc_train_class0", "asc_train_s_class0", "asc_sm_s_class0"
    ),
    `Class 1` = c(
      "b_cost_class1", "b_time_class1", "b_time_s_class1",
      "asc_car_class1", "asc_car_s_class1", "asc_train_class1",
      "asc_train_s_class1", "asc_sm_s_class1"
    )
  ))
  biogeme <- biogeme_module$BIOGEME(
    database,
    log_probability,
    number_of_draws = as.integer(number_of_draws),
    seed = as.integer(seed),
    calculating_second_derivatives = "never",
    group_of_parameters = parameter_groups,
    generate_html = FALSE,
    generate_yaml = FALSE,
    save_iterations = FALSE
  )
  biogeme$model_name <- "b15a_panel_discrete"
  results <- biogeme$estimate()
  native_data <- reticulate::py_to_r(database$dataframe)
  list(
    results = reticulate::py_to_r(bridge$extract_estimation_results(results)),
    number_of_rows = nrow(native_data),
    number_of_individuals = length(unique(native_data$ID))
  )
}

test_that("b15a panel discrete mixture 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, panel = TRUE)
  classes <- 0:1
  b_cost <- lapply(classes, function(class) biogeme_beta(paste0("b_cost_class", class), start = 0))
  b_time <- lapply(classes, function(class) biogeme_beta(paste0("b_time_class", class), start = 0))
  b_time_s <- lapply(classes, function(class) biogeme_beta(paste0("b_time_s_class", class), start = 1))
  b_time_rnd <- Map(
    function(class, location, scale) location + scale * draw(paste0("b_time_rnd_class", class), "NORMAL_ANTI"),
    classes, b_time, b_time_s
  )
  asc_car <- lapply(classes, function(class) biogeme_beta(paste0("asc_car_class", class), start = 0))
  asc_car_s <- lapply(classes, function(class) biogeme_beta(paste0("asc_car_s_class", class), start = 1))
  asc_car_rnd <- Map(
    function(class, location, scale) location + scale * draw(paste0("asc_car_rnd_class", class), "NORMAL_ANTI"),
    classes, asc_car, asc_car_s
  )
  asc_train <- lapply(classes, function(class) biogeme_beta(paste0("asc_train_class", class), start = 0))
  asc_train_s <- lapply(classes, function(class) biogeme_beta(paste0("asc_train_s_class", class), start = 1))
  asc_train_rnd <- Map(
    function(class, location, scale) location + scale * draw(paste0("asc_train_rnd_class", class), "NORMAL_ANTI"),
    classes, asc_train, asc_train_s
  )
  asc_sm <- lapply(classes, function(class) biogeme_beta(paste0("asc_sm_class", class), start = 0, fixed = TRUE))
  asc_sm_s <- lapply(classes, function(class) biogeme_beta(paste0("asc_sm_s_class", class), start = 1))
  asc_sm_rnd <- Map(
    function(class, location, scale) location + scale * draw(paste0("asc_sm_rnd_class", class), "NORMAL_ANTI"),
    classes, asc_sm, asc_sm_s
  )
  b_time_rnd[[1L]] <- 0
  score_class_0 <- biogeme_beta("score_class_0", start = -1.7)
  probability_class_0 <- logit_probability(list(`0` = score_class_0, `1` = 0), alternative = 0)
  probability_class_1 <- logit_probability(list(`0` = score_class_0, `1` = 0), alternative = 1)
  utilities <- lapply(seq_along(classes), function(index) list(
    `1` = asc_train_rnd[[index]] + b_time_rnd[[index]] * variable("TRAIN_TT_SCALED") + b_cost[[index]] * variable("TRAIN_COST_SCALED"),
    `2` = asc_sm_rnd[[index]] + b_time_rnd[[index]] * variable("SM_TT_SCALED") + b_cost[[index]] * variable("SM_COST_SCALED"),
    `3` = asc_car_rnd[[index]] + b_time_rnd[[index]] * variable("CAR_TT_SCALED") + b_cost[[index]] * variable("CAR_CO_SCALED")
  ))
  availability <- list(`1` = variable("TRAIN_AV_SP"), `2` = variable("SM_AV"), `3` = variable("CAR_AV_SP"))
  trajectories <- lapply(utilities, function(class_utilities) panel_likelihood_trajectory(logit_probability(
    class_utilities, availability, variable("CHOICE")
  )))
  formula <- log(monte_carlo(probability_class_0 * trajectories[[1L]] + probability_class_1 * trajectories[[2L]]))
  draws <- list(
    biogeme_draws("b_time_rnd_class1", "NORMAL_ANTI", 128L, 1223L),
    biogeme_draws("asc_car_rnd_class0", "NORMAL_ANTI", 128L, 1223L),
    biogeme_draws("asc_car_rnd_class1", "NORMAL_ANTI", 128L, 1223L),
    biogeme_draws("asc_train_rnd_class0", "NORMAL_ANTI", 128L, 1223L),
    biogeme_draws("asc_train_rnd_class1", "NORMAL_ANTI", 128L, 1223L),
    biogeme_draws("asc_sm_rnd_class0", "NORMAL_ANTI", 128L, 1223L),
    biogeme_draws("asc_sm_rnd_class1", "NORMAL_ANTI", 128L, 1223L)
  )
  model <- biogeme_model(
    database = database,
    formula = formula,
    draws = draws,
    control = biogeme_control(
      model_name = "b15a_panel_discrete",
      number_of_draws = 128L,
      seed = 1223L,
      second_derivatives = "never",
      group_of_parameters = list(
        `Class 0` = c("b_cost_class0", "asc_car_class0", "asc_car_s_class0", "asc_train_class0", "asc_train_s_class0", "asc_sm_s_class0"),
        `Class 1` = c("b_cost_class1", "b_time_class1", "b_time_s_class1", "asc_car_class1", "asc_car_s_class1", "asc_train_class1", "asc_train_s_class1", "asc_sm_s_class1")
      ),
      generate_html = FALSE,
      generate_yaml = FALSE,
      save_iterations = FALSE
    )
  )

  temporary_directory <- tempfile("rbiogeme-b15a-")
  dir.create(temporary_directory, recursive = TRUE)
  original_directory <- getwd()
  setwd(temporary_directory)
  on.exit(setwd(original_directory), add = TRUE)

  r_fit <- estimate(model, model_name = "b15a_panel_discrete", control = model$control)
  native <- native_swissmetro_b15a(data, 128L, 1223L)
  expect_equal(nobs(r_fit), native$number_of_individuals)
  expect_identical(r_fit$beta_names, native$results$beta_names)
  expect_equal(unname(coef(r_fit)), native$results$beta_values, tolerance = 1e-7)
  expect_equal(as.numeric(logLik(r_fit)), native$results$final_log_likelihood, tolerance = 1e-7)
  expect_equal(r_fit$number_of_excluded_data, native$results$number_of_excluded_data)
  expect_identical(isTRUE(r_fit$convergence), isTRUE(native$results$convergence))
})