tests/testthat/test-tabts.R

In R4VN: Health Data Analysis and Publication-Ready Reporting

test_that("tabts basic structure is created", {
  skip_if_not_installed("forecast")

  d <- data.frame(
    month = seq.Date(as.Date("2018-01-01"), by = "month", length.out = 60),
    y = 100 + sin(2*pi*(1:60)/12) * 20 + seq_len(60)/5
  )

  x <- tabts(
    y, time = month, data = d,
    show = FALSE, plot = FALSE, interpret = TRUE
  )

  expect_s3_class(x, "r4vn_tabts")
  expect_true(is.data.frame(x$summary))
  expect_true(is.list(x$tables))
  expect_true(!is.null(x$model))
  expect_true(is.data.frame(x$diagnostics))
})

test_that("forecast table has separate interval columns", {
  skip_if_not_installed("forecast")

  set.seed(1)
  d <- data.frame(
    month = seq.Date(as.Date("2018-01-01"), by = "month", length.out = 72),
    y = 50 + sin(2*pi*(1:72)/12) * 10 + rnorm(72)
  )

  x <- tabts(
    y, time = month, data = d,
    forecast = 6, level = c(.80, .95),
    show = FALSE, plot = FALSE
  )

  expect_equal(nrow(x$forecast), 6)
  expect_true(all(c(
    "PI80_lower", "PI80_upper",
    "PI95_lower", "PI95_upper"
  ) %in% names(x$forecast)))
})

test_that("manual SARIMA works", {
  skip_if_not_installed("forecast")

  set.seed(2)
  d <- data.frame(
    month = seq.Date(as.Date("2017-01-01"), by = "month", length.out = 84),
    y = 100 + sin(2*pi*(1:84)/12) * 15 + rnorm(84, 0, 3)
  )

  x <- tabts(
    y, time = month, data = d,
    model = "arima",
    order = c(1, 0, 1),
    seasonal = c(0, 1, 1),
    period = 12,
    show = FALSE, plot = FALSE
  )

  expect_s3_class(x, "r4vn_tabts")
  expect_equal(x$model_name, "arima")
  expect_true(is.data.frame(x$coefficients))
})

test_that("Poisson ITS returns separated CI columns and interpretation evidence", {
  set.seed(3)
  n <- 72
  month <- seq.Date(as.Date("2019-01-01"), by = "month", length.out = n)
  tt <- 0:(n-1)
  post <- as.integer(tt >= 36)
  after <- pmax(0, tt - 36)
  mu <- exp(log(60) + .004*tt - .25*post - .01*after)
  d <- data.frame(
    month = month,
    y = rpois(n, mu),
    population = rep(100000, n)
  )

  x <- tabts(
    y, time = month, data = d,
    model = "its",
    intervention = as.Date("2022-01-01"),
    family = "poisson",
    population = population,
    correlation = "none",
    interpret = TRUE,
    show = FALSE, plot = FALSE
  )

  expect_s3_class(x, "r4vn_tabts")
  expect_true(all(c("CI_lower", "CI_upper") %in% names(x$coefficients)))
  expect_true(is.data.frame(x$interpretation))
  expect_true(all(c("Finding", "Evidence", "Source") %in% names(x$interpretation)))
})

test_that("missing time interpolation regularizes series", {
  skip_if_not_installed("forecast")

  set.seed(4)
  d <- data.frame(
    month = seq.Date(as.Date("2020-01-01"), by = "month", length.out = 48),
    y = rnorm(48)
  )
  d <- d[-10, ]

  x <- tabts(
    y, time = month, data = d,
    missing_time = "interpolate",
    show = FALSE, plot = FALSE
  )

  expect_equal(nrow(x$data), 48)
  expect_false(anyNA(x$data$.y))
})

test_that("grouped ordinary analysis returns one result per group", {
  skip_if_not_installed("forecast")

  set.seed(5)
  n <- 48
  d <- rbind(
    data.frame(
      month = seq.Date(as.Date("2020-01-01"), by = "month", length.out = n),
      group = "A", y = rnorm(n, 10, 1)
    ),
    data.frame(
      month = seq.Date(as.Date("2020-01-01"), by = "month", length.out = n),
      group = "B", y = rnorm(n, 12, 1)
    )
  )

  x <- tabts(
    y, time = month, group = group, data = d,
    show = FALSE, plot = FALSE
  )

  expect_s3_class(x, "r4vn_tabts_grouped")
  expect_equal(sort(names(x$results)), c("A", "B"))
})


test_that("fixed SARIMA with d + D >= 2 does not request an invalid drift", {
  skip_if_not_installed("forecast")

  set.seed(6)
  d <- data.frame(
    month = seq.Date(as.Date("2017-01-01"), by = "month", length.out = 84),
    y = 100 + sin(2*pi*(1:84)/12) * 15 + cumsum(rnorm(84, 0, 2))
  )

  expect_no_warning(
    x <- tabts(
      y, time = month, data = d,
      model = "arima",
      order = c(1, 1, 1),
      seasonal = c(0, 1, 1),
      period = 12,
      show = FALSE, plot = FALSE
    )
  )

  expect_s3_class(x, "r4vn_tabts")
  expect_equal(x$model_name, "arima")
})

test_that("automatic ITS family can fit negative binomial counts", {
  skip_if_not_installed("MASS")

  set.seed(7)
  n <- 84
  month <- seq.Date(as.Date("2018-01-01"), by = "month", length.out = n)
  tt <- 0:(n - 1)
  post <- as.integer(tt >= 48)
  after <- pmax(0, tt - 48)
  mu <- exp(log(65) + .002 * tt - .20 * post - .006 * after)
  y <- rnbinom(n, size = .7, mu = mu)

  d <- data.frame(
    month = month,
    y = y,
    population = rep(100000, n)
  )

  x <- tabts(
    y, time = month, data = d,
    model = "its",
    intervention = as.Date("2022-01-01"),
    family = "auto",
    population = population,
    correlation = "none",
    show = FALSE, plot = FALSE
  )

  expect_s3_class(x, "r4vn_tabts")
  expect_identical(x$its$family, "negativebinomial")
  expect_s3_class(x$model, "negbin")
  expect_true(all(c("CI_lower", "CI_upper") %in% names(x$coefficients)))
})


test_that("tabts user-facing interpretation is English by default", {
  skip_if_not_installed("forecast")

  set.seed(8)
  d <- data.frame(
    month = seq.Date(as.Date("2019-01-01"), by = "month", length.out = 60),
    y = 80 + 8 * sin(2 * pi * (1:60) / 12) + rnorm(60, 0, 2)
  )

  x <- tabts(
    y, time = month, data = d,
    show = FALSE, plot = FALSE, interpret = TRUE
  )

  expect_s3_class(x, "r4vn_tabts")
  expect_true(is.data.frame(x$interpretation))
  expect_true(all(nzchar(x$interpretation$Finding)))
  expect_error(
    tabts(
      y, time = month, data = d,
      language = "vi", show = FALSE, plot = FALSE
    ),
    "currently supports only"
  )
})


test_that("negative-binomial ITS uses the MASS default log link safely", {
  skip_if_not_installed("MASS")

  set.seed(9)
  n <- 72
  month <- seq.Date(as.Date("2019-01-01"), by = "month", length.out = n)
  tt <- 0:(n - 1)
  post <- as.integer(tt >= 36)
  after <- pmax(0, tt - 36)
  mu <- exp(log(55) + .003 * tt - .18 * post - .005 * after)

  d <- data.frame(
    month = month,
    y = rnbinom(n, size = .6, mu = mu),
    population = rep(100000, n)
  )

  expect_no_error(
    x <- tabts(
      y, time = month, data = d,
      model = "its",
      intervention = as.Date("2022-01-01"),
      family = "negativebinomial",
      population = population,
      correlation = "none",
      show = FALSE, plot = FALSE
    )
  )

  expect_s3_class(x, "r4vn_tabts")
  expect_identical(x$its$family, "negativebinomial")
  expect_s3_class(x$model, "negbin")
})

test_that("tabts uses the canonical R4VN active data", {
  skip_if_not_installed("forecast")
  skip_if_not(exists("usedf", mode = "function"))

  set.seed(10)
  d <- data.frame(
    month = seq.Date(as.Date("2020-01-01"), by = "month", length.out = 60),
    y = 60 + 8 * sin(2 * pi * (1:60) / 12) + rnorm(60, 0, 2)
  )

  usedf(d, quiet = TRUE)
  on.exit(
    try(usedf(clear = TRUE, quiet = TRUE), silent = TRUE),
    add = TRUE
  )

  x <- tabts(
    y,
    time = month,
    forecast = 3,
    show = FALSE,
    plot = FALSE
  )

  expect_s3_class(x, "r4vn_tabts")
  expect_equal(nrow(x$data), nrow(d))
  expect_equal(nrow(x$forecast), 3L)
})

test_that("tabts v6 revision is loaded", {
  expect_true(exists(".r4vn_tabts_revision", inherits = TRUE))
  expect_identical(get(".r4vn_tabts_revision", inherits = TRUE), "2026-09-02-v6")
})



test_that("interpretation is off by default", {
  d <- data.frame(
    month = seq.Date(as.Date("2020-01-01"), by = "month", length.out = 48),
    y = 30 + sin(2 * pi * (1:48) / 12)
  )
  x <- tabts(
    y, time = month, data = d,
    model = "arima", order = c(1, 0, 0),
    show = FALSE, plot = FALSE
  )
  expect_null(x$interpretation)
  expect_false("interpretation" %in% names(x$tables))
})


test_that("dependency-light automatic ARIMA helper fits a finite model", {
  set.seed(20260902)
  y <- stats::ts(50 + 0.15 * (1:48) + sin(2 * pi * (1:48) / 12) + rnorm(48), frequency = 12)
  fit <- R4VN:::.r4vn_tabts_auto_arima_base(
    yts = y,
    xmat = NULL,
    frequency = 12,
    seasonal = NULL,
    options = list(include_mean = TRUE)
  )
  expect_s3_class(fit, "Arima")
  expect_true(isTRUE(attr(fit, "r4vn_auto_base")))
  expect_true(is.finite(attr(fit, "r4vn_auto_aicc")))
  expect_equal(length(attr(fit, "r4vn_auto_spec")), 7L)
})


test_that("model information reports the final specification", {
  set.seed(123)
  d <- data.frame(
    month = seq.Date(as.Date("2020-01-01"), by = "month", length.out = 48),
    y = 20 + sin(2 * pi * (1:48) / 12) + rnorm(48)
  )
  x <- tabts(
    y, time = month, data = d,
    model = "arima", order = c(1, 0, 1),
    show = FALSE, plot = FALSE
  )
  expect_true(is.data.frame(x$model_info))
  expect_true(all(c("Model", "N", "Seasonal_period", "AIC", "AICc", "BIC") %in% names(x$model_info)))
  expect_match(x$model_info$Model[1], "ARIMA")
  expect_true("model" %in% names(x$tables))
})


test_that("ITS pre and post counts include the correct intervention boundary", {
  n <- 72L
  tt <- 0:(n - 1L)
  d <- data.frame(
    month = seq.Date(as.Date("2019-01-01"), by = "month", length.out = n),
    # Add a small deterministic oscillation so the test exercises the
    # intervention boundary without creating an essentially perfect fit.
    y = 40 + 0.1 * tt - 5 * as.integer(tt >= 36) + 0.2 * sin(tt / 3)
  )
  x <- tabts(
    y, time = month, data = d,
    model = "its", intervention = as.Date("2022-01-01"),
    family = "gaussian", correlation = "none",
    show = FALSE, plot = FALSE
  )
  expect_equal(x$its$n_pre, 36L)
  expect_equal(x$its$n_post, 36L)
})


test_that("base plot objects can render without ggplot2", {
  d <- data.frame(
    .time_original = seq.Date(as.Date("2021-01-01"), by = "month", length.out = 24),
    .y = 1:24,
    .fitted = 1:24 + 0.1,
    .residual = rep(c(-0.1, 0.1), 12)
  )
  p <- R4VN:::.r4vn_tabts_plot_spec(
    "series", d = d, forecast_tab = NULL, counterfactual = NULL,
    intervention_value = NULL,
    options = R4VN:::.r4vn_ts_defaults()$plot,
    outcome_label = "Outcome", time_label = "Month"
  )
  expect_s3_class(p, "r4vn_tabts_plot")
  f <- tempfile(fileext = ".png")
  expect_no_error(R4VN:::.r4vn_tabts_save_plot_object(p, f, width = 5, height = 3, dpi = 96))
  expect_true(file.exists(f))
  expect_gt(file.info(f)$size, 0)
})


test_that("Viewer HTML embeds all generated figures", {
  set.seed(456)
  d <- data.frame(
    month = seq.Date(as.Date("2020-01-01"), by = "month", length.out = 48),
    y = 40 + 5 * sin(2 * pi * (1:48) / 12) + rnorm(48)
  )
  x <- tabts(
    y, time = month, data = d,
    model = "arima", order = c(1, 0, 1), forecast = 3,
    plots = c("series", "forecast", "residual"),
    show = FALSE, plot = TRUE
  )
  html <- R4VN:::.r4vn_tabts_build_html(x)
  expect_match(html, "Figures", fixed = TRUE)
  expect_match(html, "data:image/png;base64,", fixed = TRUE)
  expect_match(html, "Final model specification", fixed = TRUE)
})


test_that("grouped Viewer HTML contains every group", {
  set.seed(789)
  n <- 36L
  d <- rbind(
    data.frame(month = seq.Date(as.Date("2020-01-01"), by = "month", length.out = n), group = "A", y = rnorm(n, 10)),
    data.frame(month = seq.Date(as.Date("2020-01-01"), by = "month", length.out = n), group = "B", y = rnorm(n, 12))
  )
  x <- tabts(
    y, time = month, group = group, data = d,
    model = "arima", order = c(1, 0, 0),
    plots = "series", show = FALSE, plot = TRUE
  )
  html <- R4VN:::.r4vn_tabts_build_html(x)
  expect_match(html, "Group: A", fixed = TRUE)
  expect_match(html, "Group: B", fixed = TRUE)
  expect_gte(length(gregexpr("data:image/png;base64,", html, fixed = TRUE)[[1L]]), 2L)
})


test_that("controlled ITS reports readable interaction labels and requested effects", {
  set.seed(20260903)
  n <- 60L
  month <- seq.Date(as.Date("2019-01-01"), by = "month", length.out = n)
  tt <- 0:(n - 1L)
  post <- as.integer(tt >= 36L)
  after <- pmax(0, tt - 36L)
  d <- rbind(
    data.frame(
      month = month, group = factor("Control", levels = c("Control", "Intervention")),
      y = 50 + 0.15 * tt + rnorm(n, 0, 1)
    ),
    data.frame(
      month = month, group = factor("Intervention", levels = c("Control", "Intervention")),
      y = 52 + 0.15 * tt - 4 * post - 0.08 * after + rnorm(n, 0, 1)
    )
  )

  x <- tabts(
    y, time = month, group = group, data = d,
    model = "its", intervention = as.Date("2022-01-01"),
    family = "gaussian", correlation = "none",
    its_options = list(controlled = TRUE, effect_at = c(0, 3, 6, 12)),
    show = FALSE, plot = FALSE
  )

  expect_true(isTRUE(x$its$controlled))
  expect_identical(x$its$reference_group, "Control")
  expect_true(is.data.frame(x$effect_at))
  expect_true(all(c("Comparison", "Estimate", "CI_lower", "CI_upper", "p") %in% names(x$effect_at)))
  expect_true(all(x$effect_at$Comparison == "Intervention vs Control"))
  expect_true(any(grepl("Immediate level change × Group: Intervention", x$coefficients$Term_label, fixed = TRUE)))
  expect_false(any(grepl(".int1", x$coefficients$Term_label, fixed = TRUE)))
})


test_that("dependency-light grouped base plots do not connect groups into one line", {
  d <- rbind(
    data.frame(
      .time_original = seq.Date(as.Date("2021-01-01"), by = "month", length.out = 12),
      .y = 1:12, .fitted = 1:12 + .1, .residual = rep(.1, 12), .group = factor("A", levels = c("A", "B"))
    ),
    data.frame(
      .time_original = seq.Date(as.Date("2021-01-01"), by = "month", length.out = 12),
      .y = 11:22, .fitted = 11:22 + .1, .residual = rep(-.1, 12), .group = factor("B", levels = c("A", "B"))
    )
  )
  p <- R4VN:::.r4vn_tabts_plot_spec(
    "series", d = d, forecast_tab = NULL, counterfactual = NULL,
    intervention_value = as.Date("2021-07-01"),
    options = R4VN:::.r4vn_ts_defaults()$plot,
    outcome_label = "Outcome", time_label = "Month"
  )
  f <- tempfile(fileext = ".png")
  expect_no_error(R4VN:::.r4vn_tabts_save_plot_object(p, f, width = 6, height = 5, dpi = 96))
  expect_true(file.exists(f))
  expect_gt(file.info(f)$size, 0)
})