tests/testthat/test-loglik_fun.R

In readsdr: Translate Models from System Dynamics Software into 'R'

test_that("get_par_list() returns the expected list", {

  unknown_pars <- list(list(par_name = "par_beta", min = 0),
                       list(par_name = "inv_phi", min = 0, type = "meas_par",
                            par_trans = "inv"))

  actual <- get_par_list(unknown_pars)

  expected <- list(list(par_name = "par_beta", par_trans = "exp"),
                   list(par_name = "phi", par_trans = c("exp", "inv")))

  expect_equal(actual, expected)
})

test_that("transform_pars() returns the expected text", {

  unknown_pars <- list(list(par_name = "par_beta", min = 0),
                       list(par_name = "inv_phi", min = 0))

  actual_text <-  transform_pars(unknown_pars)

  expected_text <- paste(
    "  pars[[1]] <- exp(pars[[1]])",
    "  pars[[2]] <- exp(pars[[2]])", sep = "\n")

  expect_equal(actual_text, expected_text)

  unknown_pars <- list(list(par_name = "par_beta"),
                       list(par_name = "inv_phi"))

  actual_text <- transform_pars(unknown_pars)

  expected_text <- paste(
    "  pars[[1]] <- pars[[1]]",
    "  pars[[2]] <- pars[[2]]", sep = "\n")

  expect_equal(actual_text, expected_text)

  unknown_pars <- list(list(par_name = "par_beta"),
                       list(par_name = "inv_phi", min = 0))

  expected_text <- paste(
    "  pars[[1]] <- pars[[1]]",
    "  pars[[2]] <- exp(pars[[2]])", sep = "\n")

  actual_text <- transform_pars(unknown_pars)

  expect_equal(actual_text, expected_text)

  # Expit

  unknown_pars <- list(list(par_name = "par_beta", min = 0),
                       list(par_name = "par_rho", min = 0, max = 1),
                       list(par_name = "inv_phi", min = 0))

  actual_text <- transform_pars(unknown_pars)

  expected_text <- paste(
    "  pars[[1]] <- exp(pars[[1]])",
    "  pars[[2]] <- expit(pars[[2]])",
    "  pars[[3]] <- exp(pars[[3]])",
    sep = "\n")


  expect_equal(actual_text, expected_text)
})

# assign_pars_text ()-------------------------------------------------------------
test_that("assign_pars_text() returns the expected text", {

  unk_constants <- list(list(par_name = "par_beta", min = 0),
                        list(par_name = "par_rho", min = 0, max = 1))

  actual_text <- assign_pars_text(unk_constants)

  expected_text <- paste('  ds_inputs$consts["par_beta"] <- pars[[1]]',
                         '  ds_inputs$consts["par_rho"] <- pars[[2]]', sep = "\n")

  expect_equal(actual_text, expected_text)
})

# get_model_run_text() ---------------------------------------------------------

test_that("get_model_run_text returns the expected text", {

  actual_text <- get_model_run_text("rk4")

  expected_text <- paste(
    "  readsdr_env <- list2env(as.list(ds_inputs$consts))",
    "  ds_inputs$stocks <- purrr::map_dbl(ds_inputs$stocks, function(x) {",
    "    eval(parse(text = x), envir = readsdr_env)",
    "  })",
    '  o    <- sd_simulate(ds_inputs, integ_method = "rk4")',
    "  o_df <- data.frame(o)",
    sep = "\n")

  expect_equal(actual_text, expected_text)
})

# get_meas_model_text() --------------------------------------------------------

test_that("get_meas_model_text() returns the expected text", {

  # As-is measurement

  meas_data_mdl <- list(list(formula      = "y ~ poisson(C)",
                             measurements = 1:10))

  n_consts <- 1

  unknown_pars <- list(list(par_name = "par_beta", min = 0))

  actual_text <- get_meas_model_text(meas_data_mdl, n_consts, unknown_pars,
                                     FALSE)

  expected_text <- paste(
    'sim_data_1 <- dplyr::filter(o_df, time - trunc(time) == 0)',
    "loglik_1   <- sum(dpois(data_1, lambda = sim_data_1[, 'C'] + 1e-05, log = TRUE))",
    "loglik     <- loglik_1",
    "loglik", sep = "\n")

  expect_equal(actual_text, expected_text)

  # Net flow

  meas_data_mdl <- list(list(formula      = "y ~ neg_binomial_2(net_flow(C), phi)",
                             measurements = 1:10))

  n_consts <- 1

  unknown_pars <- list(list(par_name = "par_beta", min = 0),
                       list(par_name = "inv_phi", min = 0, par_trans = "inv"))

  actual_text <- get_meas_model_text(meas_data_mdl, n_consts, unknown_pars,
                                     FALSE)

  expected_text <- paste(
    'sim_data_1 <- sd_net_change(o_df, "C")',
    "loglik_1   <- sum(dnbinom(data_1, mu = sim_data_1[, 'value'] + 1e-05, size = 1/pars[[2]], log = TRUE))",
    "loglik     <- loglik_1",
    "loglik",
    sep = "\n")

  expect_equal(actual_text, expected_text)

  # Reflect log-lik

  actual_text <- get_meas_model_text(meas_data_mdl, n_consts, unknown_pars, TRUE)

  expected_text <- paste(
    'sim_data_1 <- sd_net_change(o_df, "C")',
    "loglik_1   <- sum(dnbinom(data_1, mu = sim_data_1[, 'value'] + 1e-05, size = 1/pars[[2]], log = TRUE))",
    "loglik     <- loglik_1",
    "-loglik",
    sep = "\n")
})

test_that("get_meas_model_text() handles a known par in the measurement model", {

  meas_data_mdl <- list(list(formula      = "y ~ neg_binomial_2(net_flow(C), 10)",
                             measurements = 1:10))

  n_consts <- 1

  unknown_pars <- list(list(par_name = "par_beta", min = 0))

  actual_text <- get_meas_model_text(meas_data_mdl, n_consts, unknown_pars,
                                     FALSE)

  expected_text <- paste(
    'sim_data_1 <- sd_net_change(o_df, "C")',
    "loglik_1   <- sum(dnbinom(data_1, mu = sim_data_1[, 'value'] + 1e-05, size = 10, log = TRUE))",
    "loglik     <- loglik_1",
    "loglik",
    sep = "\n")

  expect_equal(actual_text, expected_text)
})

test_that("get_meas_model_text() handles fixed pars", {

  meas_data_mdl <- list(list(formula      = "y ~ neg_binomial_2(net_flow(C), phi)",
                             measurements = 1:10))

  n_consts      <- 1

  unknown_pars  <- list(list(par_name = "par_beta", min = 0))

  supplied_pars <- "phi"

  actual_text <- get_meas_model_text(meas_data_mdl, n_consts, unknown_pars,
                                     FALSE, supplied_pars)

  expected_text <- paste(
    'sim_data_1 <- sd_net_change(o_df, "C")',
    "loglik_1   <- sum(dnbinom(data_1, mu = sim_data_1[, 'value'] + 1e-05, size = phi, log = TRUE))",
    "loglik     <- loglik_1",
    "loglik",
    sep = "\n")

  expect_equal(actual_text, expected_text)
})

test_that("get_meas_model_text() handles multiple measurements", {

  n_consts <- 1

  meas_data_mdl <- list(list(formula      = "y_A ~ neg_binomial_2(net_flow(C_A), phi)",
                             measurements = 1:10),
                        list(formula      = "y_B ~ neg_binomial_2(net_flow(C_B), phi)",
                             measurements = 10:20))

  unknown_pars <- list(list(par_name = "par_beta", min = 0),
                       list(par_name = "inv_phi", min = 0, par_trans = "inv"))

  actual_text <- get_meas_model_text(meas_data_mdl, n_consts, unknown_pars,
                                     FALSE)


  expected_text <-      paste(
    'sim_data_1 <- sd_net_change(o_df, "C_A")',
    "loglik_1   <- sum(dnbinom(data_1, mu = sim_data_1[, 'value'] + 1e-05, size = 1/pars[[2]], log = TRUE))",
    'sim_data_2 <- sd_net_change(o_df, "C_B")',
    "loglik_2   <- sum(dnbinom(data_2, mu = sim_data_2[, 'value'] + 1e-05, size = 1/pars[[2]], log = TRUE))",
    'loglik     <- loglik_1 + loglik_2',
    'loglik', sep = "\n")

  expect_equal(actual_text, expected_text)
})

# sd_loglik_fun() --------------------------------------------------------------

test_that("sd_loglik_fun() returns the expected object", {

  filepath      <- system.file("models/", "SEIR.stmx", package = "readsdr")

  unknown_pars  <- list(list(par_name = "par_beta", min = 0))

  meas_data_mdl <- list(list(formula      = "y ~ neg_binomial_2(net_flow(C), phi)",
                             measurements = 1:10))

  fun_obj <- sd_loglik_fun(filepath, unknown_pars, meas_data_mdl, neg_log = FALSE,
                           start_time = 0, stop_time = 10, timestep = 1/32)

  actual_val <- fun_obj$fun(c(1, 0.1))

  expected_val <- -32.47283

  expect_equal(actual_val, expected_val, tolerance = 1e-4)

  actual_list <- fun_obj$par_list

  expected_list <- list(list(par_name = "par_beta", par_trans = "exp"),
                        list(par_name = "phi", par_trans = c("exp", "inv")))

  expect_equal(actual_list, expected_list)

  # Test negative loglik
  fun_obj <- sd_loglik_fun(filepath, unknown_pars, meas_data_mdl,
                           neg_log = TRUE, start_time = 0, stop_time = 10,
                           timestep = 1/32)

  actual_val <- fun_obj$fun(c(1, 0.1))

  expected_val <- 32.47283

  expect_equal(actual_val, expected_val, tolerance = 1e-4)
})

# Multiple meas
test_that("sd_loglik_fun() handles multiple measurements", {

  filepath <- system.file("models/", "SEIR_age.stmx", package = "readsdr")

  unknown_pars  <- list(list(par_name = "k_AA", min = 0),
                        list(par_name = "par_rho", min = 0, max = 1))

  meas_data_mdl <- list(list(formula      = "y_A ~ neg_binomial_2(net_flow(C_A), phi)",
                             measurements = 1:10),
                        list(formula      = "y_B ~ neg_binomial_2(net_flow(C_B), phi)",
                             measurements = 11:20),
                        list(formula      = "y_C ~ neg_binomial_2(net_flow(C_C), phi)",
                             measurements = 21:30),
                        list(formula      = "y_D ~ neg_binomial_2(net_flow(C_D), phi)",
                             measurements = 31:40))

  fun_obj <- sd_loglik_fun(filepath, unknown_pars, meas_data_mdl, neg_log = FALSE,
                           start_time = 0, stop_time = 10, timestep = 1/32)

  actual_val <- fun_obj$fun(c(4, 0.5,0.1))

  expected_val <- -1168.551

  expect_equal(actual_val, expected_val, tolerance = 1e-4)
})

# Fixed pars

test_that("sd_loglik_fun() handles fixed pars", {

  filepath      <- system.file("models/", "SEIR.stmx", package = "readsdr")

  unknown_pars  <- list(list(par_name = "par_beta", min = 0))

  supplied_pars    <- c("par_rho", "I0", "phi")

  meas_data_mdl <- list(list(formula      = "y ~ neg_binomial_2(net_flow(C), phi)",
                             measurements = 1:10))

  fun_obj <- sd_loglik_fun(filepath, unknown_pars, meas_data_mdl, neg_log = FALSE,
                           supplied_pars = supplied_pars, start_time = 0,
                           stop_time = 10, timestep = 1/32)

  actual_val <- fun_obj$fun(c(1), list(par_rho = 0.75,
                                       I0 = 1,
                                       inv_phi = exp(0.1)))

  expected_val <- -32.47283

  expect_equal(actual_val, expected_val, tolerance = 1e-4)

  actual_list <- fun_obj$par_list

  expected_list <- list(list(par_name = "par_beta", par_trans = "exp"),
                        list(par_name = "phi", par_trans = c("exp", "inv")))
})

test_that("sd_loglik_fun() overrides sim params", {

  filepath      <- system.file("models/", "SEIR.stmx", package = "readsdr")
  unknown_pars  <- list(list(par_name = "par_beta", min = 0))
  meas_data_mdl <- list(list(formula      = "y ~ neg_binomial_2(net_flow(C), phi)",
                             measurements = 1:10))

  fun_obj <- sd_loglik_fun(filepath, unknown_pars, meas_data_mdl,
                           start_time = 0, stop_time = 10, timestep = 1/128)

  expect_obj <- list(start = 0,
                     stop  = 10,
                     dt    = 1/128)

  actual_obj <- fun_obj$ds_inputs$sim_params

  expect_equal(actual_obj, expect_obj)
})

test_that("sd_loglik_fun() overrides consts", {

  filepath      <- system.file("models/", "SEIR.stmx", package = "readsdr")

  unknown_pars  <- list(list(par_name = "par_beta", min = 0))
  meas_data_mdl <- list(list(formula      = "y ~ neg_binomial_2(net_flow(C), phi)",
                             measurements = 1:10))

  N_val <- 5234

  const_list <- list(N = N_val)

  fun_obj <- sd_loglik_fun(filepath, unknown_pars, meas_data_mdl,
                           start_time = 0, stop_time = 10, timestep = 1/128,
                           const_list = const_list)

  expect_equal(fun_obj$ds_inputs$consts[["N"]], N_val)
})