tests/testthat/test-tabscore.R

In R4VN: Health Data Analysis and Publication-Ready Reporting

testthat::test_that("tabscore builds a complete logistic clinical score", {
  set.seed(1001)
  n <- 450
  d <- data.frame(
    age = rnorm(n, 50, 11),
    htn = factor(rbinom(n, 1, .30), 0:1, c("No", "Yes")),
    smoke = factor(rbinom(n, 1, .25), 0:1, c("No", "Yes"))
  )
  lp <- -4.7 + .04*d$age + .8*(d$htn == "Yes") + .6*(d$smoke == "Yes")
  d$event <- rbinom(n, 1, plogis(lp))

  s <- tabscore("event", c("age", "htn", "smoke"), data = d,
                validate = "none", show = FALSE, plot = FALSE)

  testthat::expect_s3_class(s, "r4vn_tabscore")
  testthat::expect_true(all(c("Predictor", "Category", "Point") %in% names(s$tables$scorecard)))
  testthat::expect_equal(tail(s$tables$scorecard$Predictor, 1), "Total score")
  testthat::expect_true(nrow(s$tables$risk) >= 2)
  testthat::expect_true(all(s$attainable_scores %in% s$tables$risk_raw$Score))
  testthat::expect_true(is.finite(s$selected_cutoff))
  testthat::expect_true(all(c("Original model", "Model score", "Clinical score") %in% s$tables$comparison$Model))
})

testthat::test_that("manual cut points accept ASCII aliases for publication labels", {
  set.seed(1002)
  n <- 400
  d <- data.frame(
    age = round(rnorm(n, 50, 10)),
    htn = factor(rbinom(n, 1, .3), 0:1, c("No", "Yes"))
  )
  d$y <- rbinom(n, 1, plogis(-4 + .05*d$age + .9*(d$htn == "Yes")))

  s <- tabscore(
    "y", c("age", "htn"), data = d,
    cuts = list(age = c(40, 50, 60)),
    points = list(
      age = c("<40" = 0, "40-49" = 1, "50-59" = 2, ">=60" = 3),
      htn = c("No" = 0, "Yes" = 2)
    ),
    validate = "none", show = FALSE, plot = FALSE
  )

  a <- s$.dictionary[s$.dictionary$predictor == "age", ]
  testthat::expect_equal(a$clinical_point, 0:3)
  testthat::expect_equal(unname(s$score_range), c(0, 5))
})

testthat::test_that("probability threshold maps to an integer score cutoff", {
  set.seed(1003)
  n <- 500
  d <- data.frame(x = rnorm(n), z = factor(rbinom(n, 1, .4)))
  d$y <- rbinom(n, 1, plogis(-1.7 + .9*d$x + .7*(d$z == "1")))

  s <- tabscore("y", c("x", "z"), data = d,
                cutoff = "risk", riskcut = .20,
                validate = "none", show = FALSE, plot = FALSE)

  testthat::expect_match(s$selected_cutoff_method, "Clinical risk probability")
  testthat::expect_true(s$selected_cutoff %in% s$attainable_scores)
  testthat::expect_true(all(c("Sensitivity", "Specificity") %in% s$tables$cutoff_performance$Measure))
})

testthat::test_that("reference probability comparison is returned", {
  set.seed(1004)
  n <- 350
  d <- data.frame(x = rnorm(n), z = factor(rbinom(n, 1, .4)))
  d$y <- rbinom(n, 1, plogis(-1.2 + .7*d$x + .5*(d$z == "1")))
  d$ref <- plogis(-1.1 + .65*d$x + .45*(d$z == "1"))

  s <- tabscore("y", c("x", "z"), data = d,
                refprob = "ref", refcut = .20, cutoff = "refprob",
                validate = "none", show = FALSE, plot = FALSE)

  testthat::expect_true(nrow(s$tables$reference_probability) == 4)
  testthat::expect_match(s$selected_cutoff_method, "Reference probability")
})

testthat::test_that("predict applies the frozen scorecard to new patients", {
  set.seed(1005)
  n <- 300
  d <- data.frame(
    age = rnorm(n, 50, 10),
    htn = factor(rbinom(n, 1, .3), 0:1, c("No", "Yes"))
  )
  d$y <- rbinom(n, 1, plogis(-4 + .05*d$age + .8*(d$htn == "Yes")))
  s <- tabscore("y", c("age", "htn"), data = d,
                cuts = list(age = c(40, 50, 60)),
                validate = "none", show = FALSE, plot = FALSE)

  nd <- data.frame(age = c(45, 65), htn = factor(c("No", "Yes"), levels = c("No", "Yes")))
  p <- predict(s, nd, type = "all")
  testthat::expect_equal(nrow(p), 2)
  testthat::expect_true(all(c("Score", "Predicted_risk", "Risk_group") %in% names(p)))
  testthat::expect_true(all(p$Predicted_risk >= 0 & p$Predicted_risk <= 1))
})

testthat::test_that("an already-fitted logistic model can be converted", {
  set.seed(1006)
  n <- 300
  d <- data.frame(age = rnorm(n, 50, 10), smoke = factor(rbinom(n, 1, .3)))
  d$y <- rbinom(n, 1, plogis(-3 + .04*d$age + .7*(d$smoke == "1")))
  m <- stats::glm(y ~ age + smoke, data = d, family = stats::binomial())
  s <- tabscore(m, validate = "none", show = FALSE, plot = FALSE)
  testthat::expect_s3_class(s, "r4vn_tabscore")
  testthat::expect_true(nrow(s$tables$scorecard) > 1)
})

testthat::test_that("Poisson count score returns expected values", {
  set.seed(1007)
  n <- 350
  d <- data.frame(age = rnorm(n, 50, 10), smoke = factor(rbinom(n, 1, .3)))
  d$count <- rpois(n, exp(-.8 + .015*d$age + .4*(d$smoke == "1")))
  s <- tabscore("count", c("age", "smoke"), data = d, family = "poisson",
                validate = "none", show = FALSE, plot = FALSE)
  testthat::expect_true("Expected value (95% CI)" %in% names(s$tables$risk))
  testthat::expect_true(all(c("RMSE", "MAE") %in% names(s$tables$comparison)))
})

testthat::test_that("Cox scorecard produces time-specific risk", {
  testthat::skip_if_not_installed("survival")
  set.seed(1008)
  n <- 450
  d <- data.frame(age = rnorm(n, 50, 10), htn = factor(rbinom(n, 1, .3)))
  true_t <- rexp(n, rate = exp(-3 + .02*d$age + .5*(d$htn == "1")))
  cen <- rexp(n, rate = .08)
  d$status <- as.integer(true_t <= cen)
  d$ftime <- pmin(true_t, cen)

  s <- tabscore("status", c("age", "htn"), data = d,
                family = "cox", time = "ftime", times = c(1, 3, 5), cutoff_time = 5,
                validate = "none", show = FALSE, plot = FALSE)
  testthat::expect_true(all(c(1, 3, 5) %in% s$tables$risk_raw$Time))
  testthat::expect_true("C_index" %in% names(s$tables$comparison))
  testthat::expect_true(is.finite(s$selected_cutoff))
})

testthat::test_that("R4VN vars() specifications and forced predictors are accepted", {
  set.seed(1009)
  n <- 320
  d <- data.frame(
    age = round(rnorm(n, 50, 10)),
    smoke = factor(rbinom(n, 1, .3), 0:1, c("No", "Yes")),
    bmi = rnorm(n, 24, 4)
  )
  d$y <- rbinom(n, 1, plogis(-4 + .05*d$age + .7*(d$smoke == "Yes")))
  s <- tabscore("y", vars(c.age, b2.smoke, c.bmi), data = d,
                select = "purposeful", force = vars(c.age, b2.smoke),
                validate = "none", show = FALSE, plot = FALSE)
  testthat::expect_true("age" %in% s$selected_predictors)
  testthat::expect_equal(levels(stats::model.frame(s$models$original)$smoke)[1], "Yes")
  testthat::expect_true(all(c("Final model", "Clinical scorecard", "Score to risk") %in%
                              names(s$publication_tables)))
})

testthat::test_that("R4VN logistic result objects can be converted through raw$model", {
  testthat::skip_if_not(exists("logistic", mode = "function", inherits = TRUE))
  set.seed(1010)
  n <- 300
  d <- data.frame(
    age = round(rnorm(n, 50, 10)),
    smoke = factor(rbinom(n, 1, .3), 0:1, c("No", "Yes"))
  )
  d$y <- factor(rbinom(n, 1, plogis(-3.5 + .045*d$age + .7*(d$smoke == "Yes"))),
                0:1, c("No", "Yes"))
  mr4 <- logistic(y, age, smoke, data = d, event = "Yes", show = FALSE)
  s <- tabscore(mr4, validate = "none", show = FALSE, plot = FALSE)
  testthat::expect_s3_class(s, "r4vn_tabscore")
  testthat::expect_equal(s$event, "Yes")
  testthat::expect_true(all(c("age", "smoke") %in% s$selected_predictors))
})


testthat::test_that("external validation keeps the frozen binary score and cutoff", {
  set.seed(1011)
  n <- 500
  d <- data.frame(age = round(rnorm(n, 50, 10)), htn = factor(rbinom(n, 1, .3), 0:1, c("No", "Yes")))
  d$y <- rbinom(n, 1, plogis(-4 + .05*d$age + .8*(d$htn == "Yes")))
  dev <- d[1:300, ]
  val <- d[301:500, ]
  s <- tabscore("y", c("age", "htn"), data = dev,
                cuts = list(age = c(40, 50, 60)),
                validation = val, validate = "none", show = FALSE, plot = FALSE)
  testthat::expect_true(is.data.frame(s$tables$external_validation))
  testthat::expect_true(all(c("AUC", "Brier", "Calibration_intercept",
                              "Sensitivity", "Specificity") %in% names(s$tables$external_validation)))
  testthat::expect_equal(s$tables$external_validation$N[1], nrow(val))
})

testthat::test_that("riskonly rebases a protective logistic contrast to an add-only score", {
  set.seed(1101)
  n <- 1400
  d <- data.frame(
    age = rnorm(n, 50, 10),
    sex = factor(rbinom(n, 1, .50), 0:1, c("Female", "Male"))
  )
  d$y <- rbinom(n, 1, plogis(-4 + .035*d$age + .90*(d$sex == "Male")))

  s <- tabscore(y, vars(c.age, b2.sex), data = d,
                cuts = list(age = c(40, 50, 60)),
                riskonly = TRUE, scoreref = "lowest",
                validate = "none", show = FALSE, plot = FALSE)

  sx <- s$.dictionary[s$.dictionary$predictor == "sex", , drop = FALSE]
  female <- which(sx$category == "Female")
  male <- which(sx$category == "Male")

  testthat::expect_true(sx$model_reference[male])
  testthat::expect_lt(sx$model_ratio[female], 1)
  testthat::expect_true(sx$scoring_reference[female])
  testthat::expect_equal(sx$scoring_ratio[female], 1, tolerance = 1e-8)
  testthat::expect_gt(sx$scoring_ratio[male], 1)
  testthat::expect_true(all(s$.dictionary$scoring_ratio >= 1 - 1e-10))
  testthat::expect_true(all(s$.dictionary$clinical_point >= 0))
  testthat::expect_true("Risk-oriented OR" %in% names(s$tables$risk_orientation))
})

testthat::test_that("riskonly handles a multi-level predictor without using abs(beta)", {
  set.seed(1102)
  n <- 1800
  d <- data.frame(
    grp = factor(sample(c("A", "B", "C"), n, TRUE), levels = c("A", "B", "C"))
  )
  # B has the highest risk, C the lowest. b2.grp makes B the model reference.
  lp <- -2.8 + .40*(d$grp == "A") + 1.00*(d$grp == "B") + 0*(d$grp == "C")
  d$y <- rbinom(n, 1, plogis(lp))

  s <- tabscore(y, vars(b2.grp), data = d,
                riskonly = TRUE, validate = "none",
                show = FALSE, plot = FALSE)
  g <- s$.dictionary[s$.dictionary$predictor == "grp", , drop = FALSE]

  testthat::expect_equal(g$category[g$scoring_reference], "C")
  testthat::expect_true(all(g$scoring_ratio >= 1 - 1e-10))
  testthat::expect_true(all(g$clinical_point >= 0))
  testthat::expect_gt(g$scoring_ratio[g$category == "B"],
                      g$scoring_ratio[g$category == "A"])
})

testthat::test_that("negative manual protective points are shifted to add-only form by default", {
  set.seed(1103)
  n <- 1000
  d <- data.frame(
    exercise = factor(rbinom(n, 1, .50), 0:1, c("No", "Yes"))
  )
  d$y <- rbinom(n, 1, plogis(-1.6 - .85*(d$exercise == "Yes")))

  s <- tabscore(y, c(exercise), data = d,
                points = list(exercise = c("No" = 0, "Yes" = -2)),
                riskonly = TRUE,
                validate = "none", show = FALSE, plot = FALSE)
  e <- s$.dictionary[s$.dictionary$predictor == "exercise", , drop = FALSE]
  testthat::expect_equal(e$clinical_point[e$category == "Yes"], 0)
  testthat::expect_equal(e$clinical_point[e$category == "No"], 2)
  testthat::expect_true(all(e$clinical_point >= 0))

  s_signed <- tabscore(y, c(exercise), data = d,
                       points = list(exercise = c("No" = 0, "Yes" = -2)),
                       riskonly = FALSE, scoreref = "model",
                       validate = "none", show = FALSE, plot = FALSE)
  es <- s_signed$.dictionary[s_signed$.dictionary$predictor == "exercise", , drop = FALSE]
  testthat::expect_equal(es$clinical_point[es$category == "Yes"], -2)
})

testthat::test_that("riskonly rejects a non-lowest custom scoring reference", {
  set.seed(1104)
  n <- 900
  d <- data.frame(sex = factor(rbinom(n, 1, .50), 0:1, c("Female", "Male")))
  d$y <- rbinom(n, 1, plogis(-2 + .9*(d$sex == "Male")))

  testthat::expect_error(
    tabscore(y, c(sex), data = d,
             riskonly = TRUE, scoreref = list(sex = "Male"),
             validate = "none", show = FALSE, plot = FALSE),
    "lowest-risk category"
  )
})

testthat::test_that("Cox riskonly reorients protective HRs to HR >= 1", {
  testthat::skip_if_not_installed("survival")
  set.seed(1105)
  n <- 1200
  d <- data.frame(sex = factor(rbinom(n, 1, .50), 0:1, c("Female", "Male")))
  true_t <- rexp(n, rate = exp(-2.7 + .8*(d$sex == "Male")))
  cen <- rexp(n, rate = .08)
  d$status <- as.integer(true_t <= cen)
  d$ftime <- pmin(true_t, cen)

  s <- tabscore(status, vars(b2.sex), data = d, family = "cox", time = ftime,
                times = c(1, 3), cutoff_time = 3,
                riskonly = TRUE, validate = "none",
                show = FALSE, plot = FALSE)
  sx <- s$.dictionary[s$.dictionary$predictor == "sex", , drop = FALSE]
  testthat::expect_true(all(sx$scoring_ratio >= 1 - 1e-10))
  testthat::expect_true(all(sx$clinical_point >= 0))
  testthat::expect_equal(s$effect_measure, "HR")
})

testthat::test_that("plot metadata contains all available publication figures", {
  set.seed(1201)
  n <- 500
  d <- data.frame(
    age = rnorm(n, 50, 10),
    bmi = rnorm(n, 24, 4),
    smoke = factor(rbinom(n, 1, .30), 0:1, c("No", "Yes"))
  )
  d$y <- rbinom(n, 1, plogis(-4 + .045*d$age + .05*d$bmi + .70*(d$smoke == "Yes")))

  s <- tabscore(y, c(age, bmi, smoke), data = d,
                cuts = list(age = c(40, 50, 60), bmi = c(20, 25, 30)),
                validate = "none", plot = TRUE, show = FALSE)

  testthat::expect_true(isTRUE(s$plot_enabled))
  testthat::expect_true(all(c("risk", "roc", "calibration", "decision", "distribution") %in% s$plots))
  testthat::expect_true(all(s$plots %in% names(s$plot_titles)))
  testthat::expect_true(all(vapply(s$plot_data[s$plots], is.data.frame, logical(1))))
})

testthat::test_that("Viewer embeds figures when plot TRUE and omits them when plot FALSE", {
  if (!(capabilities("png") || capabilities("cairo"))) {
    testthat::skip("No raster/vector graphics capability for Viewer smoke test")
  }
  set.seed(1202)
  n <- 420
  d <- data.frame(
    age = rnorm(n, 48, 11),
    htn = factor(rbinom(n, 1, .35), 0:1, c("No", "Yes"))
  )
  d$y <- rbinom(n, 1, plogis(-3.5 + .045*d$age + .8*(d$htn == "Yes")))

  s1 <- tabscore(y, c(age, htn), data = d,
                 cuts = list(age = c(40, 50, 60)),
                 validate = "none", plot = TRUE, show = FALSE)
  p1 <- R4VN:::.r4vn_score_show_html(s1, show = FALSE)
  on.exit(unlink(p1), add = TRUE)
  h1 <- paste(readLines(p1, warn = FALSE, encoding = "UTF-8"), collapse = "\n")
  testthat::expect_true(grepl("Figure 1", h1, fixed = TRUE))
  testthat::expect_true(grepl("Score-to-risk curve", h1, fixed = TRUE))
  testthat::expect_true(grepl("data:image/png;base64|<svg", h1))

  s0 <- tabscore(y, c(age, htn), data = d,
                 cuts = list(age = c(40, 50, 60)),
                 validate = "none", plot = FALSE, show = FALSE)
  p0 <- R4VN:::.r4vn_score_show_html(s0, show = FALSE)
  on.exit(unlink(p0), add = TRUE)
  h0 <- paste(readLines(p0, warn = FALSE, encoding = "UTF-8"), collapse = "\n")
  testthat::expect_false(grepl("<section class='figure'>", h0, fixed = TRUE))
})

testthat::test_that("ROC and calibration plotting are constrained to probability axes", {
  set.seed(1203)
  n <- 450
  d <- data.frame(age = rnorm(n, 52, 12))
  d$y <- rbinom(n, 1, plogis(-3 + .05*d$age))
  s <- tabscore(y, c(age), data = d,
                cuts = list(age = c(40, 50, 60)),
                validate = "none", plot = FALSE, show = FALSE)

  f <- tempfile(fileext = ".pdf")
  grDevices::pdf(f)
  on.exit({
    try(grDevices::dev.off(), silent = TRUE)
    unlink(f)
  }, add = TRUE)
  testthat::expect_silent(plot(s, which = "roc"))
  testthat::expect_silent(plot(s, which = "calibration"))
  grDevices::dev.off()
})