tests/testthat/test-stan_model.R

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

test_that("stan_model() returns the expected string", {

  m1      <- "y ~ normal(net_flow(C), tau)"
  meas_mdl <- list(m1)

  estimated_params <- list(
    sd_prior("par_beta", "lognormal", c(0, 1)),
    sd_prior("par_rho", "beta", c(2, 2)),
    sd_prior("I0", "lognormal", c(0, 1), "init"),
    list(par_name = "tau", dist = "exponential", beta = 0.1, min = 0,
         type = "meas_par"))

  lvl_names <- c("S", "E", "I", "R", "C")

  actual <- stan_model(estimated_params, meas_mdl, lvl_names)

  expected <- paste(
    "model {",
    "  par_beta ~ lognormal(0, 1);",
    "  par_rho ~ beta(2, 2);",
    "  I0 ~ lognormal(0, 1);",
    "  tau ~ exponential(0.1);",
    "  y ~ normal(delta_x_1, tau);",
    "}", sep = "\n")

  expect_equal(actual, expected)
})

test_that("stan_model() returns the expected string for a vectorised model", {

  ag <- c("A", "B", "C", "D") # age_groups

  measurements <- stringr::str_glue("y_{ag} ~ poisson(net_flow(C_{ag}))")
  meas_mdl     <- as.list(measurements)

  estimated_params <- list(sd_prior("par_rho", "beta", c(2, 2)))

  filepath  <- system.file("models/", "SEIR_age.stmx", package = "readsdr")
  mdl       <- read_xmile(filepath)
  lvl_names <- sd_stocks(mdl)$name

  actual <- stan_model(estimated_params, meas_mdl, lvl_names)

  expected <- paste(
    "model {",
    "  par_rho ~ beta(2, 2);",
    "  y_A ~ poisson(delta_x_1);",
    "  y_B ~ poisson(delta_x_2);",
    "  y_C ~ poisson(delta_x_3);",
    "  y_D ~ poisson(delta_x_4);",
    "}", sep = "\n")

  expect_equal(actual, expected)
})

test_that("stan_model() handles the lognormal distribution", {

  estimated_params  <- list(sd_prior("par_alpha", "normal",
                                     c(1, 0.5), min_0 = TRUE))

  lvl_names <- c("Hares", "Lynx")

  meas_mdl <- list("y1 ~ lognormal(Lynx, sigma_1)")
  actual   <- stan_model(estimated_params, meas_mdl, lvl_names)

  expected <- paste(
    "model {",
    "  par_alpha ~ normal(1, 0.5);",
    "  y1 ~ lognormal(x[:, 2], sigma_1);",
    "}", sep = "\n")

  expect_equal(actual, expected)
})

test_that("stan_model() handles measurements of init values", {

  estimated_params  <- list(
    list(
      par_name  = "sigma_1",
      dist      = "lognormal",
      mu        = -1,
      sigma     = 1,
      min       = 0,
      type      = "meas_par"),
    list(
      par_name = "H0",
      dist     = "lognormal",
      mu       = log(10),
      sigma    = 1,
      min      = 0,
      type     = "init"),
    list(
      par_name = "par_gamma",
      dist     = "normal",
      mu       = 1,
      sigma    = 0.5,
      type     = "constant",
      min      = 0))

  meas_mdl <- list("y ~ lognormal(Hares, sigma_1)",
                   "y0 ~ lognormal(log(Hares[0]), sigma_1)")

  lvl_names <- c("Hares", "Lynx")
  actual    <- stan_model(estimated_params, meas_mdl, lvl_names)

  expected <- paste(
    "model {",
    "  sigma_1 ~ lognormal(-1, 1);",
    "  H0 ~ lognormal(2.30258509299405, 1);",
    "  par_gamma ~ normal(1, 0.5);",
    "  y ~ lognormal(x[:, 1], sigma_1);",
    "  y0 ~ lognormal(log(x0[1]), sigma_1);",
    "}", sep = "\n")

  expect_equal(actual, expected)
})

# Construct prior line----------------------------------------------------------
test_that("construct_prior_line() returns the expected string", {

  prior_obj <- list(par_name = "par_beta",
                    dist     = "lognormal",
                    mu       =  0,
                    sigma    =  1,
                    min      =  0,
                    type     =  "constant")

  actual   <- construct_prior_line(prior_obj)

  expected <- "  par_beta ~ lognormal(0, 1);"

  expect_equal(actual, expected)
})

test_that("construct_likelihood_line() returns the expected string", {

  meas_obj <- "y ~ neg_binomial_2(net_flow(C), phi)"

  actual   <- construct_likelihood_line(meas_obj, 1)

  expected <- list(line          = "  y ~ neg_binomial_2(delta_x_1, phi);",
                   delta_counter = 2)

  expect_equal(actual, expected)

  meas_obj <- "y ~ poisson(C)"

  lvl_names <- c("S", "E", "I", "R", "C")

  actual   <- construct_likelihood_line(meas_obj, 1, lvl_names)

  expected <- list(line          = "  y ~ poisson(x[:, 5]);",
                   delta_counter = 1)

  expect_equal(actual, expected)
})

test_that("construct_likelihood_line() handles log transformations", {

  meas_obj <- "y1 ~ lognormal(log(Lynx), sigma_1)"
  actual   <- construct_likelihood_line(meas_obj, 1, c("Hares", "Lynx"))

  expected <- list(line          = "  y1 ~ lognormal(log(x[:, 2]), sigma_1);",
                   delta_counter = 1)

  expect_equal(actual, expected)
})