# ============================================================================
# TRIAL-LIKE VS REAL-WORLD EVIDENCE ANALYSIS
# ============================================================================
# Developer: Planimeter Kft.
# Version:   1.0.0
# Date:      2026-09-23
#
# Copyright (C) 2026 Planimeter Kft.
#
# License:
#   This program is free software: you can 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.
#
#   This program is distributed in the hope that it will be useful, but
#   WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
#   or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
#   for more details.
#
#   SPDX-License-Identifier: GPL-3.0-or-later
#
# AI use statement:
#   Generative AI was used as a support tool during code development, review
#   and validation. The analytical design, final code, interpretation and
#   responsibility for the released materials remain with Planimeter Kft.
#
# Purpose:
#   Reproduce the public Planimeter analysis comparing a trial-like subset
#   with a broader real-world hypertensive population using synthetic data.
#
# Required input files (relative paths):
#   data/rwe_patient_features.csv
#   data/rwe_htn_outcome.csv
#
# Main outputs:
#   figures/01_synthea_ecosystem.png
#   figures/02_trial_eligibility_funnel.png
#   figures/03_trial_vs_rwe_baseline.png
#   figures/04_prediabetes_diabetes_risk.png
#   figures/05_forest_trial_vs_rwe.png
#   figures/06_transportability_concept.png
#   output/eligibility_funnel.csv
#   output/prediabetes_diabetes_risk.csv
#   output/model_results.csv
#   output/interaction_result.csv
#   output/adjusted_model_results.csv
# ============================================================================

# ---- 0. Packages -----------------------------------------------------------

required_packages <- c(
  "dplyr",
  "tidyr",
  "ggplot2",
  "readr",
  "scales"
)

missing_packages <- required_packages[
  !vapply(required_packages, requireNamespace, logical(1), quietly = TRUE)
]

if (length(missing_packages) > 0) {
  stop(
    paste0(
      "Missing required package(s): ",
      paste(missing_packages, collapse = ", "),
      ". Install them before running this script."
    )
  )
}

library(dplyr)
library(tidyr)
library(ggplot2)
library(readr)
library(scales)

# ---- 1. Relative paths -----------------------------------------------------

DATA_DIR   <- "data"
FIGURE_DIR <- "figures"
OUTPUT_DIR <- "output"

dir.create(FIGURE_DIR, recursive = TRUE, showWarnings = FALSE)
dir.create(OUTPUT_DIR, recursive = TRUE, showWarnings = FALSE)

patient_features_path <- file.path(DATA_DIR, "rwe_patient_features.csv")
htn_outcome_path      <- file.path(DATA_DIR, "rwe_htn_outcome.csv")

if (!file.exists(patient_features_path) || !file.exists(htn_outcome_path)) {
  stop(
    paste0(
      "Required input files were not found. Run the script from the project root ",
      "and keep the two CSV files in the data/ directory."
    )
  )
}

# ---- 2. Read analysis-ready data ------------------------------------------

patient_features <- read_csv(
  patient_features_path,
  show_col_types = FALSE
)

htn_outcome <- read_csv(
  htn_outcome_path,
  show_col_types = FALSE
)

required_patient_columns <- c(
  "patient_id",
  "age_at_analysis_date",
  "n_conditions",
  "n_medications",
  "n_encounters",
  "n_emergency",
  "n_inpatient",
  "followup_years",
  "crit_age",
  "crit_alive",
  "crit_followup",
  "crit_cv_risk",
  "crit_comorbidity",
  "crit_polypharmacy",
  "crit_no_stroke",
  "crit_no_mi",
  "crit_no_hf",
  "crit_no_ckd",
  "crit_no_cancer",
  "crit_utilization",
  "trial_eligible"
)

required_outcome_columns <- c(
  "patient_id",
  "trial_eligible",
  "pred_at_index",
  "age_index",
  "gender",
  "n_conditions",
  "n_medications",
  "n_emergency",
  "n_inpatient",
  "diabetes_5y"
)

if (!all(required_patient_columns %in% names(patient_features))) {
  stop("rwe_patient_features.csv does not contain all required columns.")
}

if (!all(required_outcome_columns %in% names(htn_outcome))) {
  stop("rwe_htn_outcome.csv does not contain all required columns.")
}

# ---- 3. Visual theme -------------------------------------------------------

bg     <- "#07111f"
grid   <- "#24384F"
txt    <- "white"
subtxt <- "#A9BED6"
cyan   <- "#55DDE0"
purple <- "#7B61FF"
gold   <- "#FFD166"
muted  <- "#334A63"
red    <- "#FF6B6B"

theme_rwe <- function(base_size = 18) {
  theme_minimal(base_size = base_size) +
    theme(
      plot.background = element_rect(fill = bg, colour = NA),
      panel.background = element_rect(fill = bg, colour = NA),
      panel.grid.major = element_line(colour = grid, linewidth = 0.35),
      panel.grid.minor = element_blank(),
      plot.title = element_text(
        colour = txt,
        face = "bold",
        size = base_size + 10
      ),
      plot.subtitle = element_text(
        colour = subtxt,
        size = base_size + 1
      ),
      axis.title = element_text(colour = subtxt),
      axis.text = element_text(colour = subtxt),
      legend.background = element_rect(fill = bg, colour = NA),
      legend.key = element_rect(fill = bg, colour = NA),
      legend.text = element_text(colour = subtxt),
      legend.title = element_text(colour = txt),
      plot.margin = margin(30, 40, 30, 40)
    )
}

theme_rwe_void <- function(base_size = 18) {
  theme_void(base_size = base_size) +
    theme(
      plot.background = element_rect(fill = bg, colour = NA),
      panel.background = element_rect(fill = bg, colour = NA),
      plot.title = element_text(
        colour = txt,
        face = "bold",
        size = base_size + 10
      ),
      plot.subtitle = element_text(
        colour = subtxt,
        size = base_size + 1
      ),
      plot.margin = margin(35, 45, 35, 45)
    )
}

save_rwe <- function(p, filename, width = 14, height = 8) {
  ggsave(
    filename = file.path(FIGURE_DIR, filename),
    plot = p,
    width = width,
    height = height,
    dpi = 300,
    bg = bg
  )
}

# ---- 4. Eligibility funnel -------------------------------------------------

eligibility_funnel <- tibble(
  step = factor(
    c(
      "All patients",
      "Age window",
      "CV risk",
      "Comorbidity",
      "Medication burden",
      "No major exclusions",
      "Trial eligible"
    ),
    levels = c(
      "All patients",
      "Age window",
      "CV risk",
      "Comorbidity",
      "Medication burden",
      "No major exclusions",
      "Trial eligible"
    )
  ),
  count = c(
    nrow(patient_features),
    sum(patient_features$crit_age, na.rm = TRUE),
    sum(
      patient_features$crit_age &
        patient_features$crit_alive &
        patient_features$crit_followup &
        patient_features$crit_cv_risk,
      na.rm = TRUE
    ),
    sum(
      patient_features$crit_age &
        patient_features$crit_alive &
        patient_features$crit_followup &
        patient_features$crit_cv_risk &
        patient_features$crit_comorbidity,
      na.rm = TRUE
    ),
    sum(
      patient_features$crit_age &
        patient_features$crit_alive &
        patient_features$crit_followup &
        patient_features$crit_cv_risk &
        patient_features$crit_comorbidity &
        patient_features$crit_polypharmacy,
      na.rm = TRUE
    ),
    sum(
      patient_features$crit_age &
        patient_features$crit_alive &
        patient_features$crit_followup &
        patient_features$crit_cv_risk &
        patient_features$crit_comorbidity &
        patient_features$crit_polypharmacy &
        patient_features$crit_no_stroke &
        patient_features$crit_no_mi &
        patient_features$crit_no_hf &
        patient_features$crit_no_ckd &
        patient_features$crit_no_cancer,
      na.rm = TRUE
    ),
    sum(patient_features$trial_eligible, na.rm = TRUE)
  )
) %>%
  mutate(
    percent_of_all = count / first(count),
    width = sqrt(count / max(count)) * 10,
    y = rev(seq_len(n())),
    xmin = -width / 2,
    xmax = width / 2,
    ymin = y - 0.33,
    ymax = y + 0.33,
    label = paste0(comma(count), "\n", step)
  )

write_csv(
  eligibility_funnel %>% select(step, count, percent_of_all),
  file.path(OUTPUT_DIR, "eligibility_funnel.csv")
)

# ---- 5. Epidemiological analyses ------------------------------------------

risk_df <- htn_outcome %>%
  mutate(
    pred_group = if_else(
      pred_at_index,
      "Hypertension + Prediabetes",
      "Hypertension only"
    )
  ) %>%
  group_by(pred_group) %>%
  summarise(
    n = n(),
    events = sum(diabetes_5y, na.rm = TRUE),
    risk = mean(diabetes_5y, na.rm = TRUE),
    .groups = "drop"
  )

write_csv(
  risk_df,
  file.path(OUTPUT_DIR, "prediabetes_diabetes_risk.csv")
)

fit_all <- glm(
  diabetes_5y ~ pred_at_index,
  data = htn_outcome,
  family = binomial()
)

fit_trial <- glm(
  diabetes_5y ~ pred_at_index,
  data = htn_outcome %>% filter(trial_eligible),
  family = binomial()
)

fit_rwe <- glm(
  diabetes_5y ~ pred_at_index,
  data = htn_outcome %>% filter(!trial_eligible),
  family = binomial()
)

fit_interaction <- glm(
  diabetes_5y ~ pred_at_index * trial_eligible,
  data = htn_outcome,
  family = binomial()
)

fit_adjusted <- glm(
  diabetes_5y ~
    pred_at_index +
    age_index +
    gender +
    n_conditions +
    n_medications,
  data = htn_outcome,
  family = binomial()
)

extract_or <- function(fit, term, population, n_obs) {
  beta <- coef(fit)[term]
  se <- sqrt(diag(vcov(fit)))[term]
  p <- summary(fit)$coefficients[term, "Pr(>|z|)"]

  tibble(
    population = population,
    n = n_obs,
    odds_ratio = exp(beta),
    ci_lower = exp(beta - 1.96 * se),
    ci_upper = exp(beta + 1.96 * se),
    p_value = p
  )
}

model_results <- bind_rows(
  extract_or(
    fit_all,
    "pred_at_indexTRUE",
    "Full HTN cohort",
    nrow(htn_outcome)
  ),
  extract_or(
    fit_trial,
    "pred_at_indexTRUE",
    "Trial-like",
    sum(htn_outcome$trial_eligible)
  ),
  extract_or(
    fit_rwe,
    "pred_at_indexTRUE",
    "RWE-only",
    sum(!htn_outcome$trial_eligible)
  )
)

write_csv(
  model_results,
  file.path(OUTPUT_DIR, "model_results.csv")
)

interaction_result <- extract_or(
  fit_interaction,
  "pred_at_indexTRUE:trial_eligibleTRUE",
  "Prediabetes x trial eligibility",
  nrow(htn_outcome)
)

write_csv(
  interaction_result,
  file.path(OUTPUT_DIR, "interaction_result.csv")
)

extract_all_terms <- function(fit) {
  coefs <- coef(fit)
  ses <- sqrt(diag(vcov(fit)))
  pvals <- summary(fit)$coefficients[, "Pr(>|z|)"]

  tibble(
    term = names(coefs),
    odds_ratio = exp(coefs),
    ci_lower = exp(coefs - 1.96 * ses),
    ci_upper = exp(coefs + 1.96 * ses),
    p_value = pvals
  )
}

write_csv(
  extract_all_terms(fit_adjusted),
  file.path(OUTPUT_DIR, "adjusted_model_results.csv")
)

# ---- 6. Figure 1: source-data ecosystem -----------------------------------

ecosystem_nodes <- data.frame(
  label = c(
    "Patients",
    "Conditions",
    "Medications",
    "Encounters",
    "Patient-level\nanalysis dataset"
  ),
  x = c(1, 2.7, 4.4, 6.1, 3.55),
  y = c(3, 3, 3, 3, 1.4),
  w = c(1.25, 1.35, 1.45, 1.35, 2.6),
  h = c(0.55, 0.55, 0.55, 0.55, 0.7),
  type = c("source", "source", "source", "source", "derived")
)

p1 <- ggplot(ecosystem_nodes) +
  geom_rect(
    aes(
      xmin = x - w / 2,
      xmax = x + w / 2,
      ymin = y - h / 2,
      ymax = y + h / 2,
      fill = type
    ),
    alpha = 0.95,
    colour = subtxt,
    linewidth = 0.7
  ) +
  geom_text(
    aes(x = x, y = y, label = label),
    colour = txt,
    fontface = "bold",
    size = 5.5,
    lineheight = 0.9
  ) +
  geom_segment(
    data = data.frame(
      x = c(1, 2.7, 4.4, 6.1),
      xend = rep(3.55, 4),
      y = rep(2.68, 4),
      yend = rep(1.8, 4)
    ),
    aes(x = x, xend = xend, y = y, yend = yend),
    arrow = grid::arrow(length = grid::unit(0.16, "inches")),
    colour = cyan,
    alpha = 0.7,
    linewidth = 0.8
  ) +
  scale_fill_manual(values = c("source" = muted, "derived" = gold)) +
  coord_cartesian(xlim = c(0.2, 6.9), ylim = c(0.7, 3.6), clip = "off") +
  labs(
    title = "Real-world evidence starts as linked patient histories",
    subtitle = paste(
      "Patients, diagnoses, treatments, and healthcare encounters are",
      "transformed into an analysis-ready cohort"
    )
  ) +
  theme_rwe_void() +
  theme(legend.position = "none")

save_rwe(p1, "01_synthea_ecosystem.png")

# ---- 7. Figure 2: eligibility funnel --------------------------------------

p2 <- ggplot(eligibility_funnel) +
  geom_rect(
    aes(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax, fill = step),
    alpha = 0.96
  ) +
  geom_text(
    aes(x = 0, y = y, label = label),
    colour = txt,
    fontface = "bold",
    size = 4.8,
    lineheight = 0.9
  ) +
  scale_fill_manual(values = c(muted, cyan, cyan, purple, purple, gold, red)) +
  coord_cartesian(xlim = c(-5.8, 5.8), ylim = c(0.3, 7.7), clip = "off") +
  labs(
    title = "Eligibility criteria reshape the real-world population",
    subtitle = "A broad synthetic population is reduced to a trial-like subset"
  ) +
  theme_rwe_void() +
  theme(legend.position = "none")

save_rwe(p2, "02_trial_eligibility_funnel.png")

# ---- 8. Figure 3: trial-like vs RWE baseline ------------------------------

baseline_df <- htn_outcome %>%
  mutate(population = if_else(trial_eligible, "Trial-like", "RWE-only")) %>%
  group_by(population) %>%
  summarise(
    Age = mean(age_index, na.rm = TRUE),
    Conditions = mean(n_conditions, na.rm = TRUE),
    Medications = mean(n_medications, na.rm = TRUE),
    Emergency = mean(n_emergency, na.rm = TRUE),
    Inpatient = mean(n_inpatient, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  pivot_longer(
    cols = -population,
    names_to = "metric",
    values_to = "value"
  ) %>%
  group_by(metric) %>%
  mutate(value_scaled = value / max(value, na.rm = TRUE)) %>%
  ungroup()

p3 <- ggplot(
  baseline_df,
  aes(metric, value_scaled, fill = population)
) +
  geom_col(
    position = position_dodge(width = 0.65),
    width = 0.55,
    alpha = 0.95
  ) +
  geom_text(
    aes(label = round(value, 1)),
    position = position_dodge(width = 0.65),
    vjust = -0.35,
    colour = txt,
    fontface = "bold",
    size = 4.5
  ) +
  scale_y_continuous(
    labels = percent_format(accuracy = 1),
    expand = expansion(mult = c(0, 0.18))
  ) +
  scale_fill_manual(values = c("Trial-like" = gold, "RWE-only" = cyan)) +
  labs(
    title = "The trial-like population is not simply a smaller copy",
    subtitle = "Baseline characteristics differ after applying eligibility criteria",
    x = NULL,
    y = "Scaled within metric",
    fill = NULL
  ) +
  theme_rwe() +
  theme(legend.position = "bottom")

save_rwe(p3, "03_trial_vs_rwe_baseline.png")

# ---- 9. Figure 4: five-year diabetes risk ---------------------------------

p4 <- ggplot(
  risk_df,
  aes(pred_group, risk * 100, fill = pred_group)
) +
  geom_col(width = 0.58, alpha = 0.96) +
  geom_text(
    aes(
      label = paste0(
        round(risk * 100, 1),
        "%\n",
        events,
        "/",
        n
      )
    ),
    vjust = -0.35,
    colour = txt,
    fontface = "bold",
    size = 6,
    lineheight = 0.9
  ) +
  scale_fill_manual(
    values = c(
      "Hypertension only" = cyan,
      "Hypertension + Prediabetes" = gold
    )
  ) +
  scale_y_continuous(
    labels = function(x) paste0(x, "%"),
    expand = expansion(mult = c(0, 0.18))
  ) +
  labs(
    title = "Baseline risk matters",
    subtitle = paste(
      "Among hypertensive patients, prediabetes was associated with higher",
      "5-year diabetes incidence"
    ),
    x = NULL,
    y = "5-year diabetes incidence"
  ) +
  theme_rwe() +
  theme(legend.position = "none")

save_rwe(p4, "04_prediabetes_diabetes_risk.png")

# ---- 10. Figure 5: effect estimates ---------------------------------------

forest_df <- model_results %>%
  mutate(
    population = factor(
      population,
      levels = c("RWE-only", "Trial-like", "Full HTN cohort")
    ),
    label = paste0(
      "OR ",
      round(odds_ratio, 2),
      "  (",
      round(ci_lower, 2),
      "-",
      round(ci_upper, 2),
      ")"
    )
  )

p5 <- ggplot(forest_df, aes(x = odds_ratio, y = population)) +
  geom_vline(
    xintercept = 1,
    colour = subtxt,
    linetype = "dashed",
    linewidth = 0.9
  ) +
  geom_errorbarh(
    aes(xmin = ci_lower, xmax = ci_upper),
    height = 0.18,
    colour = cyan,
    linewidth = 1.2
  ) +
  geom_point(colour = gold, size = 5) +
  geom_text(
    aes(label = label),
    x = max(forest_df$ci_upper, na.rm = TRUE) * 1.05,
    colour = txt,
    hjust = 0,
    fontface = "bold",
    size = 4.6
  ) +
  scale_x_log10(
    breaks = c(0.25, 0.5, 1, 2, 4, 8),
    labels = c("0.25", "0.5", "1", "2", "4", "8")
  ) +
  coord_cartesian(
    xlim = c(0.25, max(forest_df$ci_upper, na.rm = TRUE) * 1.8),
    clip = "off"
  ) +
  labs(
    title = "The direction of association is consistent",
    subtitle = paste(
      "Trial-like estimates are less precise, but not directionally",
      "contradictory to RWE"
    ),
    x = "Odds ratio for 5-year diabetes, log scale",
    y = NULL
  ) +
  theme_rwe() +
  theme(panel.grid.major.y = element_blank())

save_rwe(p5, "05_forest_trial_vs_rwe.png")

# ---- 11. Figure 6: transportability concept -------------------------------

transport_nodes <- data.frame(
  label = c(
    "Real-world\nhypertensive population",
    "Eligibility\ncriteria",
    "Trial-like\ncohort",
    "Estimate\nassociation",
    "Generalize\nback"
  ),
  x = c(1, 2.6, 4.1, 5.6, 7.1),
  y = c(2, 2, 2, 2, 2),
  w = c(1.65, 1.25, 1.35, 1.25, 1.25),
  h = c(0.72, 0.72, 0.72, 0.72, 0.72),
  fill = c(muted, purple, gold, cyan, red)
)

p6 <- ggplot(transport_nodes) +
  geom_segment(
    data = data.frame(
      x = c(1.85, 3.25, 4.8, 6.25),
      xend = c(2.0, 3.45, 5.0, 6.45),
      y = 2,
      yend = 2
    ),
    aes(x = x, xend = xend, y = y, yend = yend),
    arrow = grid::arrow(length = grid::unit(0.18, "inches")),
    colour = subtxt,
    linewidth = 1
  ) +
  geom_rect(
    aes(
      xmin = x - w / 2,
      xmax = x + w / 2,
      ymin = y - h / 2,
      ymax = y + h / 2,
      fill = label
    ),
    alpha = 0.96,
    colour = subtxt,
    linewidth = 0.7
  ) +
  geom_text(
    aes(x = x, y = y, label = label),
    colour = txt,
    fontface = "bold",
    size = 4.8,
    lineheight = 0.9
  ) +
  scale_fill_manual(values = setNames(transport_nodes$fill, transport_nodes$label)) +
  coord_cartesian(xlim = c(0.1, 8), ylim = c(1.2, 2.8), clip = "off") +
  labs(
    title = "RWE asks whether evidence travels",
    subtitle = paste(
      "The key question is not only whether an association exists, but whether",
      "it remains interpretable across populations"
    )
  ) +
  theme_rwe_void() +
  theme(legend.position = "none")

save_rwe(p6, "06_transportability_concept.png")

# ---- 12. Console summary ---------------------------------------------------

cat("\n")
cat("============================================================\n")
cat("TRIAL-LIKE VS REAL-WORLD EVIDENCE ANALYSIS COMPLETE\n")
cat("============================================================\n\n")

cat("Synthetic population: ", nrow(patient_features), "\n", sep = "")
cat(
  "Trial eligible: ",
  sum(patient_features$trial_eligible),
  " (",
  percent(mean(patient_features$trial_eligible), accuracy = 0.1),
  ")\n",
  sep = ""
)
cat("Hypertensive outcome cohort: ", nrow(htn_outcome), "\n", sep = "")
cat("Trial-like HTN cohort: ", sum(htn_outcome$trial_eligible), "\n", sep = "")
cat("RWE-only HTN cohort: ", sum(!htn_outcome$trial_eligible), "\n\n", sep = "")

print(risk_df)
cat("\n")
print(model_results)
cat("\nInteraction model:\n")
print(interaction_result)

cat("\nFigures saved to: ", FIGURE_DIR, "/\n", sep = "")
cat("Tables saved to: ", OUTPUT_DIR, "/\n", sep = "")
cat("============================================================\n")
