tests/testthat/manual-checks.R

In bage: Bayesian Estimation and Forecasting of Age-Specific Rates

## Checks that involve inspection by a human, and so are not run by testthat

library(bage)
library(dplyr)
library(poputils)
library(ggplot2)

## Datasets -------------------------------------------------------------------

p_isl_deaths <- ggplot(isl_deaths, aes(x = age_mid(age), y = deaths / popn, col = sex)) +
  facet_wrap(vars(time)) +
  geom_line() +
  scale_y_log10() +
  ggtitle("isl_deaths - rates")

p_kor_births <- ggplot(kor_births, aes(x = time, y = births / popn)) +
  facet_grid(vars(age), vars(region)) +
  geom_line() +
  ggtitle("kor_births - rates")

p_nld_expenditure <- ggplot(nld_expenditure,
                            aes(x = age_mid(age), y = value)) +
  facet_grid(vars(diag), vars(year), scale = "free_y") +
  geom_line() +
  ggtitle("nld_expenditure - expenditures")

p_nzl_divorces <- ggplot(nzl_divorces,
                         aes(x = age_mid(age), y = divorces / population, color = sex)) +
  facet_wrap(vars(time)) +
  geom_line() +
  ggtitle("nzl_divorces - rates")

p_nzl_households <- ggplot(nzl_households,
                           aes(x = age_mid(age),
                               y = oneperson / total,
                               color = factor(year))) +
  facet_wrap(vars(region)) +
  geom_line() +
  ggtitle("nzl_households - rates one-person")

p_nzl_injuries <- ggplot(nzl_injuries,
                         aes(x = age_mid(age),
                             y = injuries / popn,
                             linetype = sex,
                             color = ethnicity)) +
  facet_wrap(vars(year)) +
  geom_line() +
  ggtitle("nzl_injuries - rates")

p_prt_deaths <- ggplot(prt_deaths,
                         aes(x = age_mid(age),
                             y = deaths / exposure,
                             color = sex)) +
  facet_wrap(vars(time)) +
  geom_line() +
  scale_y_log10() +
  ggtitle("prt_deaths - rates")

p_swe_infant <- ggplot(swe_infant,
                       aes(x = time,
                           y = deaths / births)) +
  facet_wrap(vars(county)) +
  geom_line() +
  ggtitle("swe_infant - rates")


p_usa_deaths <- ggplot(usa_deaths,
         aes(x = as.Date(paste0(month, "-01"), format = "%Y-%b-%d"),
             y = deaths)) +
  geom_line() +
  ggtitle("usa_deaths - counts")


graphics.off()
pdf(file = "manual-checks-datasets.pdf",
    w = 10,
    h = 10,
    onefile = TRUE)
plot(p_isl_deaths)
plot(p_kor_births)
plot(p_nld_expenditure)
plot(p_nzl_divorces)
plot(p_nzl_households)
plot(p_nzl_injuries)
plot(p_prt_deaths)
plot(p_swe_infant)
plot(p_usa_deaths)
dev.off()

file.remove("manual-checks-datasets.pdf")


## Scaled SVDs ----------------------------------------------------------------

nms_svd <- c("CSA",
             "HFD",
             "HIMD_R", "HIMD_P1", "HIMD_P5",
             "HMD",
             "LFP",
             "WMD_C", "WMD_E")

transforms <- list(invlogit,
                   exp,
                   exp, invlogit, invlogit,
                   identity,
                   invlogit,
                   invlogit, invlogit)

has_sex_var <- function(x) "indep" %in% x$data$type

n_comp <- 5

## components

graphics.off()
pdf(file = "manual-checks-svd-components.pdf",
    w = 8,
    h = 8,
    onefile = TRUE)

for (nm in nms_svd) {
  obj <- get(nm)
  has_sex <- has_sex_var(obj)
  indep <- if (has_sex) FALSE else NULL
  comp <- components(obj, n_comp = n_comp, indep = indep)
  if (has_sex) {
    p <- ggplot(comp, aes(x = age_mid(age), y = value, color = sex))
  } else {
    p <- ggplot(comp, aes(x = age_mid(age), y = value))
  }
  p <- p +
    facet_wrap(vars(component)) +
    geom_line() +
    ggtitle(nm)
  plot(p)
}

dev.off()


## generate

set.seed(0)
graphics.off()
pdf(file = "manual-checks-svd-generate.pdf",
    w = 12,
    h = 12,
    onefile = TRUE)

for (i in seq_along(nms_svd)) {
  nm <- nms_svd[[i]]
  transform <- transforms[[i]]
  obj <- get(nm)
  has_sex <- has_sex_var(obj)
  indep <- if (has_sex) FALSE else NULL
  gen <- generate(obj, indep = indep, n_comp = n_comp)
  gen$value <- transform(gen$value)
  if (has_sex) {
    p <- ggplot(gen, aes(x = age_mid(age), y = value, color = sex))
  } else {
    p <- ggplot(gen, aes(x = age_mid(age), y = value))
  }
  p <- p +
    facet_wrap(vars(draw)) +
    geom_line() +
    ggtitle(nm)
  plot(p)
}

dev.off()

file.remove("manual-checks-svd-components.pdf")
file.remove("manual-checks-svd-generate.pdf")