# =============================================================================
# CONFIDENCE WITHOUT CERTAINTY
# NHANES 2017-2018 Prediabetes Uncertainty Demonstration
#
# Developed by: Planimeter Kft.
# Version:      1.0.0
# Date:         2026-09-24
# Copyright:    (c) 2026 Planimeter Kft.
#
# License:
#   This program is free software: you may redistribute it and/or modify it
#   under the terms of the GNU General Public License as published by the
#   Free Software Foundation, either version 3 of the License, or
#   (at your option) any later version.
#
# SPDX-License-Identifier: GPL-3.0-or-later
#
# AI use statement:
#   Generative AI tools were used to support code development, review, and
#   validation. The analytical design, final code, interpretation, and
#   responsibility for the released version remain with Planimeter Kft.
# =============================================================================


required_packages <- c("dplyr", "ggplot2", "readr", "boot", "scales")
missing_packages <- required_packages[
  !vapply(required_packages, requireNamespace, logical(1), quietly = TRUE)
]
if (length(missing_packages) > 0) {
  stop(paste("Missing required packages:", paste(missing_packages, collapse = ", ")))
}

library(dplyr)
library(ggplot2)
library(readr)
library(boot)
library(scales)

DATA_FILE <- file.path("data", "nhanes_prediabetes_analysis.csv")
REFERENCE_FILE <- file.path("reference", "confidence_summary.csv")
OUTPUT_DIR <- "output"
FIGURE_DIR <- "figures"
dir.create(OUTPUT_DIR, recursive = TRUE, showWarnings = FALSE)
dir.create(FIGURE_DIR, recursive = TRUE, showWarnings = FALSE)

# -----------------------------------------------------------------------------
# 1. Analysis data
# -----------------------------------------------------------------------------

dat <- read_csv(DATA_FILE, show_col_types = FALSE)
required_columns <- c("SEQN", "age", "hba1c", "bmi", "prediabetes")
stopifnot(all(required_columns %in% names(dat)))
stopifnot(!anyNA(dat[, required_columns]))

n <- nrow(dat)
y <- sum(dat$prediabetes)
prev <- mean(dat$prediabetes)

# -----------------------------------------------------------------------------
# 2. Participant-level nonparametric bootstrap
# -----------------------------------------------------------------------------

boot_prev <- function(data, indices) {
  mean(data$prediabetes[indices])
}

set.seed(123)
boot_res <- boot(dat, statistic = boot_prev, R = 3000)
boot_df <- data.frame(
  replicate = seq_len(nrow(boot_res$t)),
  prevalence = boot_res$t[, 1]
)
boot_ci <- quantile(boot_df$prevalence, c(0.025, 0.975), na.rm = TRUE)
write_csv(boot_df, file.path(OUTPUT_DIR, "bootstrap_prevalence_recalculated.csv"))

# -----------------------------------------------------------------------------
# 3. Beta-binomial posterior illustration
#    Beta(1,1) prior -> Beta(1+y, 1+n-y) posterior
# -----------------------------------------------------------------------------

set.seed(123)
posterior_prev <- rbeta(
  n = 10000,
  shape1 = 1 + y,
  shape2 = 1 + n - y
)
posterior_df <- data.frame(
  draw = seq_along(posterior_prev),
  prevalence = posterior_prev
)
post_ci <- quantile(posterior_prev, c(0.025, 0.5, 0.975), na.rm = TRUE)
write_csv(posterior_df, file.path(OUTPUT_DIR, "posterior_prevalence_recalculated.csv"))

# -----------------------------------------------------------------------------
# 4. Numerical summary and validation
# -----------------------------------------------------------------------------

summary_table <- data.frame(
  metric = c(
    "Analysis sample size",
    "Prediabetes cases",
    "Unweighted prevalence",
    "Bootstrap 95% lower",
    "Bootstrap 95% upper",
    "Posterior 95% lower",
    "Posterior median",
    "Posterior 95% upper",
    "Posterior P(prevalence > 25%)",
    "Posterior P(prevalence > 30%)",
    "Posterior P(prevalence > 35%)",
    "Median age",
    "Median BMI"
  ),
  value = c(
    n,
    y,
    prev,
    unname(boot_ci[1]),
    unname(boot_ci[2]),
    unname(post_ci[1]),
    unname(post_ci[2]),
    unname(post_ci[3]),
    mean(posterior_prev > 0.25),
    mean(posterior_prev > 0.30),
    mean(posterior_prev > 0.35),
    median(dat$age),
    median(dat$bmi)
  )
)
write_csv(summary_table, file.path(OUTPUT_DIR, "confidence_summary_recalculated.csv"))

if (file.exists(REFERENCE_FILE)) {
  ref <- read_csv(REFERENCE_FILE, show_col_types = FALSE)
  validation <- ref %>%
    rename(reference_value = value) %>%
    left_join(summary_table %>% rename(recalculated_value = value), by = "metric") %>%
    mutate(
      absolute_difference = abs(reference_value - recalculated_value),
      status = if_else(absolute_difference <= 1e-8, "PASS", "CHECK")
    )
  write_csv(validation, file.path(OUTPUT_DIR, "validation_checks_runtime.csv"))
}

# -----------------------------------------------------------------------------
# 5. Publication figures
# -----------------------------------------------------------------------------

bg_dark <- "#07111F"
panel_bg <- "#0B1728"
blue <- "#36A3FF"
cyan <- "#5DEBFF"
gold <- "#F6C85F"
white <- "#F7FAFC"
muted <- "#9FB3C8"
grid <- "#243447"

theme_web_dark <- function(base_size = 15) {
  theme_minimal(base_size = base_size) +
    theme(
      plot.background = element_rect(fill = bg_dark, color = NA),
      panel.background = element_rect(fill = bg_dark, color = NA),
      panel.grid.major = element_line(color = grid, linewidth = 0.25),
      panel.grid.minor = element_blank(),
      axis.text = element_text(color = muted),
      axis.title = element_text(color = muted),
      plot.title = element_text(color = white, face = "bold", size = base_size + 7),
      plot.subtitle = element_text(color = muted, size = base_size),
      plot.caption = element_text(color = muted, size = base_size - 2),
      legend.position = "none"
    )
}

save_png <- function(plot, filename, width = 9, height = 5) {
  ggsave(
    filename = file.path(FIGURE_DIR, filename),
    plot = plot,
    width = width,
    height = height,
    dpi = 320,
    bg = bg_dark
  )
}

# Figure 1: conventional point estimate + bootstrap interval
fig1 <- data.frame(estimate = prev, low = boot_ci[1], high = boot_ci[2]) %>%
  ggplot(aes(x = "Prediabetes prevalence", y = estimate)) +
  geom_col(width = 0.45, fill = "#BFC7D5") +
  geom_errorbar(aes(ymin = low, ymax = high), width = 0.12, color = "#4B5563") +
  geom_text(aes(label = percent(estimate, accuracy = 0.1)), vjust = -1.1,
            size = 7, fontface = "bold", color = "#111827") +
  scale_y_continuous(labels = percent_format(), limits = c(0, max(boot_ci) + 0.08)) +
  labs(title = "The conventional answer",
       subtitle = "A single point estimate with a bootstrap interval",
       x = NULL, y = "Estimated prevalence") +
  theme_minimal(base_size = 15) +
  theme(plot.background = element_rect(fill = "white", color = NA),
        panel.grid.minor = element_blank(), axis.text.x = element_blank(),
        axis.ticks.x = element_blank())
save_png(fig1, "01_conventional_prevalence.png")

# Figure 2: bootstrap distribution
fig2 <- ggplot(boot_df, aes(x = prevalence)) +
  geom_density(fill = blue, alpha = 0.25, color = cyan, linewidth = 1.2) +
  geom_vline(xintercept = prev, color = white, linewidth = 1.2) +
  geom_vline(xintercept = boot_ci, color = gold, linetype = "dashed", linewidth = 0.9) +
  labs(title = "One estimate became a distribution",
       subtitle = "Bootstrap resampling reveals sampling uncertainty",
       x = "Prediabetes prevalence", y = NULL,
       caption = "NHANES 2017-2018 analysed sample; unweighted demonstration") +
  scale_x_continuous(labels = percent_format()) +
  theme_web_dark()
save_png(fig2, "02_bootstrap_distribution.png")

# Figure 3: posterior cloud
set.seed(123)
cloud_df <- posterior_df[sample(seq_len(nrow(posterior_df)), 2500), ]
cloud_df$y <- rnorm(nrow(cloud_df), 0, 0.16)
cloud_df$distance <- abs(cloud_df$prevalence - median(posterior_prev))
fig3 <- ggplot(cloud_df, aes(x = prevalence, y = y)) +
  geom_point(aes(alpha = 1 - rescale(distance)), color = cyan, size = 1.8) +
  geom_vline(xintercept = post_ci[2], color = white, linewidth = 1.2) +
  geom_vline(xintercept = post_ci[c(1,3)], color = gold, linetype = "dashed") +
  scale_alpha(range = c(0.08, 0.9)) +
  scale_x_continuous(labels = percent_format()) +
  labs(title = "10,000 plausible realities",
       subtitle = "The posterior distribution makes model-based uncertainty visible",
       x = "Prediabetes prevalence", y = NULL) +
  theme_web_dark() +
  theme(axis.text.y = element_blank(), axis.ticks.y = element_blank(),
        panel.grid.major.y = element_blank())
save_png(fig3, "03_posterior_cloud.png")

# Figure 4: decision probability
threshold <- 0.30
p_decision <- mean(posterior_prev > threshold)
gauge_df <- data.frame(x = seq(0, 1, length.out = 300)) %>%
  mutate(included = x <= p_decision)
fig4 <- ggplot(gauge_df, aes(x = x, y = 1)) +
  geom_tile(aes(fill = included), height = 0.22, width = 1/300) +
  scale_fill_manual(values = c("TRUE" = cyan, "FALSE" = "#1B2A3A")) +
  annotate("text", x = 0, y = 1.42, hjust = 0,
           label = percent(p_decision, accuracy = 0.1),
           color = white, size = 18, fontface = "bold") +
  annotate("text", x = 0, y = 1.18, hjust = 0,
           label = paste0("Posterior probability that prevalence exceeds ", percent(threshold)),
           color = muted, size = 5) +
  labs(title = "From estimate to decision probability",
       subtitle = "A model-based threshold probability reframes the reporting question",
       x = NULL, y = NULL) +
  coord_cartesian(xlim = c(0,1), ylim = c(0.55,1.55), clip = "off") +
  theme_web_dark() +
  theme(axis.text = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank())
save_png(fig4, "04_probability_gauge.png")

# Figure 5: credible interval strip
thresholds <- data.frame(threshold = c(0.25, 0.30, 0.35), label = c("25%", "30%", "35%"))
fig5 <- ggplot() +
  geom_segment(aes(x = post_ci[1], xend = post_ci[3], y = 1, yend = 1),
               color = cyan, linewidth = 8, lineend = "round", alpha = 0.85) +
  geom_point(aes(x = post_ci[2], y = 1), color = white, size = 6) +
  geom_vline(data = thresholds, aes(xintercept = threshold),
             color = gold, linetype = "dashed", linewidth = 0.8) +
  geom_text(data = thresholds, aes(x = threshold, y = 1.32, label = label),
            color = gold, size = 4.5, fontface = "bold") +
  annotate("text", x = post_ci[2], y = 0.72,
           label = paste0("Median ", percent(post_ci[2], 0.1),
                          "   |   95% credible interval ",
                          percent(post_ci[1], 0.1), "-", percent(post_ci[3], 0.1)),
           color = white, size = 4.5, fontface = "bold") +
  scale_x_continuous(labels = percent_format()) +
  coord_cartesian(ylim = c(0.6,1.42)) +
  labs(title = "The answer is a range of plausible realities",
       subtitle = "Decision thresholds become visible against the posterior interval",
       x = "Prediabetes prevalence", y = NULL) +
  theme_web_dark() +
  theme(axis.text.y = element_blank(), axis.ticks.y = element_blank(),
        panel.grid.major.y = element_blank())
save_png(fig5, "05_decision_interval_strip.png")

# Figure 6: bootstrap logistic risk curve (illustrative, not a clinical model)
set.seed(123)
grid_bmi <- seq(quantile(dat$bmi, 0.05), quantile(dat$bmi, 0.95), length.out = 120)
curves <- vector("list", 400)
for (i in seq_len(400)) {
  idx <- sample(seq_len(nrow(dat)), replace = TRUE)
  d <- dat[idx, ]
  model <- glm(prediabetes ~ bmi + age, data = d, family = binomial())
  pred_data <- data.frame(bmi = grid_bmi, age = median(dat$age))
  curves[[i]] <- data.frame(
    iteration = i,
    bmi = grid_bmi,
    risk = predict(model, newdata = pred_data, type = "response")
  )
}
curves <- bind_rows(curves)
curve_summary <- curves %>%
  group_by(bmi) %>%
  summarise(low = quantile(risk, 0.05), mid = median(risk),
            high = quantile(risk, 0.95), .groups = "drop")
fig6 <- ggplot() +
  geom_line(data = curves %>% filter(iteration <= 120),
            aes(x = bmi, y = risk, group = iteration),
            color = cyan, alpha = 0.055, linewidth = 0.55) +
  geom_ribbon(data = curve_summary, aes(x = bmi, ymin = low, ymax = high),
              fill = blue, alpha = 0.18) +
  geom_line(data = curve_summary, aes(x = bmi, y = mid), color = white, linewidth = 1.4) +
  labs(title = "Risk is not a single curve",
       subtitle = "Bootstrap logistic models show uncertainty around the BMI association",
       x = "BMI", y = "Estimated probability of prediabetes") +
  scale_y_continuous(labels = percent_format()) +
  theme_web_dark()
save_png(fig6, "06_bmi_uncertainty_ribbon.png")

cat("Confidence Without Certainty analysis complete.\n")
print(summary_table)