tests/testthat/test-equivalence.R

In rbiogeme: Interface to 'Biogeme'

test_that("R estimation is numerically equivalent to native Biogeme", {
  skip_if_not(
    identical(Sys.getenv("RBIOGEME_RUN_INTEGRATION"), "1"),
    "Set RBIOGEME_RUN_INTEGRATION=1 to run estimation equivalence tests"
  )
  skip_if_not(
    rbiogeme_test_configure_python(),
    "A compatible Python Biogeme environment is not available"
  )
  original_directory <- getwd()
  temporary_directory <- tempfile("rbiogeme-equivalence-")
  dir.create(temporary_directory, recursive = TRUE)
  setwd(temporary_directory)
  on.exit(setwd(original_directory), add = TRUE)

  data <- data.frame(
    choice = c(1, 2, 1, 2, 1, 2, 1, 2),
    x = c(1, 2, 0, 1, 2, 3, 0, 4)
  )
  controls <- list(
    calculating_second_derivatives = "never",
    max_iterations = 20L
  )

  database <- biogeme_database("toy_r", data)
  beta <- biogeme_beta("b", start = 0)
  r_model <- logit_model(
    database,
    choice = "choice",
    utilities = list(`1` = beta * variable("x"), `2` = 0)
  )
  r_fit <- estimate(r_model, model_name = "r_reference", controls = controls)

  expression_module <- reticulate::import("biogeme.expressions", convert = FALSE)
  models_module <- reticulate::import("biogeme.models", convert = FALSE)
  database_module <- reticulate::import("biogeme.database", convert = FALSE)
  biogeme_module <- reticulate::import("biogeme.biogeme", convert = FALSE)

  native_database <- database_module$Database(
    "toy_native",
    reticulate::r_to_py(data)
  )
  native_beta <- expression_module$Beta("b", 0.0, NULL, NULL, 0L)
  native_utilities <- reticulate::dict(
    `1` = native_beta * expression_module$Variable("x"),
    `2` = 0
  )
  native_log_likelihood <- models_module$loglogit(
    native_utilities,
    NULL,
    expression_module$Variable("choice")
  )
  native_biogeme <- biogeme_module$BIOGEME(
    native_database,
    native_log_likelihood,
    calculating_second_derivatives = "never",
    max_iterations = 20L,
    generate_html = FALSE,
    generate_yaml = FALSE,
    save_iterations = FALSE
  )
  native_biogeme$model_name <- "native_reference"
  native_results <- native_biogeme$estimate()
  bridge <- rbiogeme:::biogeme_bridge()
  native_result <- reticulate::py_to_r(
    bridge$extract_estimation_results(native_results)
  )

  expect_equal(unname(coef(r_fit)), native_result$beta_values, tolerance = 1e-8)
  expect_equal(
    as.numeric(logLik(r_fit)),
    native_result$final_log_likelihood,
    tolerance = 1e-8
  )
  expect_identical(isTRUE(r_fit$convergence), isTRUE(native_result$convergence))
  native_vcov <- matrix(
    as.numeric(unlist(native_result$variance_covariance, use.names = FALSE)),
    nrow = length(native_result$beta_names)
  )
  expect_equal(unname(vcov(r_fit)), native_vcov, tolerance = 1e-8)
})