# ============================================================================
# RWE DATA PREPARATION FROM SYNTHEA-DERIVED RDS FILES
# ============================================================================
# Developer: Planimeter Kft.
# Version:   1.0.0
# Date:      2026-09-23
#
# Copyright (C) 2026 Planimeter Kft.
# SPDX-License-Identifier: GPL-3.0-or-later
#
# 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.
#
# 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:
#   Document the transformation of four Synthea-derived source tables into the
#   two compact analysis-ready CSV files distributed with the public analysis.
#
# Expected source files (not included in the public package):
#   raw_data/synthea_patients.rds
#   raw_data/synthea_conditions.rds
#   raw_data/synthea_medications.rds
#   raw_data/synthea_encounters.rds
#
# Outputs:
#   data/rwe_patient_features.csv
#   data/rwe_htn_outcome.csv
# ============================================================================

required_packages <- c("dplyr", "tidyr", "lubridate", "readr")
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(lubridate)
library(readr)

RAW_DIR  <- "raw_data"
DATA_DIR <- "data"

dir.create(DATA_DIR, recursive = TRUE, showWarnings = FALSE)

source_files <- c(
  patients = file.path(RAW_DIR, "synthea_patients.rds"),
  conditions = file.path(RAW_DIR, "synthea_conditions.rds"),
  medications = file.path(RAW_DIR, "synthea_medications.rds"),
  encounters = file.path(RAW_DIR, "synthea_encounters.rds")
)

if (!all(file.exists(source_files))) {
  stop(
    paste0(
      "One or more source files are missing from raw_data/. ",
      "The public analysis can still be reproduced directly from the two ",
      "analysis-ready CSV files supplied in data/."
    )
  )
}

patients    <- readRDS(source_files[["patients"]])
conditions  <- readRDS(source_files[["conditions"]])
medications <- readRDS(source_files[["medications"]])
encounters  <- readRDS(source_files[["encounters"]])

# A fixed, data-derived analysis date avoids dependence on the date on which
# the script happens to be executed.
analysis_date <- as.Date(max(encounters$STOP, na.rm = TRUE))

patients2 <- patients %>%
  mutate(
    BIRTHDATE = as.Date(BIRTHDATE),
    DEATHDATE = as.Date(DEATHDATE),
    age = floor(interval(BIRTHDATE, analysis_date) / years(1)),
    alive = is.na(DEATHDATE)
  )

comorbidity <- conditions %>%
  group_by(PATIENT) %>%
  summarise(
    n_conditions = n_distinct(DESCRIPTION),
    .groups = "drop"
  )

med_count <- medications %>%
  group_by(PATIENT) %>%
  summarise(
    n_medications = n_distinct(DESCRIPTION),
    .groups = "drop"
  )

enc_summary <- encounters %>%
  mutate(
    START = as.Date(START),
    STOP = as.Date(STOP)
  ) %>%
  group_by(PATIENT) %>%
  summarise(
    n_encounters = n(),
    n_emergency = sum(ENCOUNTERCLASS == "emergency", na.rm = TRUE),
    n_inpatient = sum(ENCOUNTERCLASS == "inpatient", na.rm = TRUE),
    first_visit = min(START, na.rm = TRUE),
    last_visit = max(STOP, na.rm = TRUE),
    followup_years = as.numeric(
      difftime(last_visit, first_visit, units = "days")
    ) / 365.25,
    .groups = "drop"
  )

patient_features <- patients2 %>%
  left_join(comorbidity, by = c("Id" = "PATIENT")) %>%
  left_join(med_count, by = c("Id" = "PATIENT")) %>%
  left_join(enc_summary, by = c("Id" = "PATIENT")) %>%
  mutate(
    across(
      c(
        n_conditions,
        n_medications,
        n_encounters,
        n_emergency,
        n_inpatient,
        followup_years
      ),
      ~ replace_na(.x, 0)
    )
  )

cv_flags <- conditions %>%
  mutate(desc_low = tolower(DESCRIPTION)) %>%
  group_by(PATIENT) %>%
  summarise(
    hypertension = any(grepl("hypertension|high blood pressure", desc_low)),
    hyperlipidemia = any(
      grepl("hyperlipidemia|hypercholesterolemia|cholesterol", desc_low)
    ),
    coronary_heart_disease = any(
      grepl("coronary|ischemic heart|ischaemic heart|angina", desc_low)
    ),
    stroke_history = any(grepl("stroke|cerebrovascular", desc_low)),
    myocardial_infarction_history = any(
      grepl("myocardial infarction|heart attack", desc_low)
    ),
    heart_failure = any(grepl("heart failure|congestive", desc_low)),
    chronic_kidney_disease = any(
      grepl("chronic kidney|renal disease|kidney disease", desc_low)
    ),
    cancer_history = any(grepl("cancer|malignant|neoplasm", desc_low)),
    .groups = "drop"
  )

patient_features <- patient_features %>%
  left_join(cv_flags, by = c("Id" = "PATIENT"))

flag_vars <- c(
  "hypertension",
  "hyperlipidemia",
  "coronary_heart_disease",
  "stroke_history",
  "myocardial_infarction_history",
  "heart_failure",
  "chronic_kidney_disease",
  "cancer_history"
)

patient_features <- patient_features %>%
  mutate(across(all_of(flag_vars), ~ replace_na(.x, FALSE)))

condition_cut <- 25
medication_cut <- 8

patient_features <- patient_features %>%
  mutate(
    crit_age = age >= 45 & age <= 75,
    crit_alive = alive,
    crit_followup = followup_years >= 2,
    crit_cv_risk =
      hypertension |
      hyperlipidemia |
      coronary_heart_disease,
    crit_comorbidity = n_conditions <= condition_cut,
    crit_polypharmacy = n_medications <= medication_cut,
    crit_no_stroke = !stroke_history,
    crit_no_mi = !myocardial_infarction_history,
    crit_no_hf = !heart_failure,
    crit_no_ckd = !chronic_kidney_disease,
    crit_no_cancer = !cancer_history,
    crit_utilization = n_emergency <= 2 & n_inpatient <= 1,
    trial_eligible =
      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
  )

cond2 <- conditions %>%
  mutate(
    START = as.Date(START),
    desc_low = tolower(DESCRIPTION)
  )

htn_index <- cond2 %>%
  filter(grepl("hypertension", desc_low)) %>%
  group_by(PATIENT) %>%
  summarise(
    index_date = min(START, na.rm = TRUE),
    .groups = "drop"
  )

pred_dates <- cond2 %>%
  filter(grepl("prediabetes", desc_low)) %>%
  group_by(PATIENT) %>%
  summarise(
    pred_date = min(START, na.rm = TRUE),
    .groups = "drop"
  )

diabetes_events <- cond2 %>%
  filter(
    grepl("diabetes", desc_low),
    !grepl("prediabetes", desc_low),
    !grepl("gestational", desc_low)
  ) %>%
  select(PATIENT, diabetes_date = START)

htn_cohort <- patient_features %>%
  inner_join(htn_index, by = c("Id" = "PATIENT")) %>%
  left_join(pred_dates, by = c("Id" = "PATIENT")) %>%
  mutate(
    age_index = floor(interval(BIRTHDATE, index_date) / years(1)),
    pred_at_index = !is.na(pred_date) & pred_date <= index_date
  ) %>%
  filter(age_index >= 40, age_index <= 85)

htn_outcome <- htn_cohort %>%
  left_join(diabetes_events, by = c("Id" = "PATIENT")) %>%
  mutate(
    diabetes_5y =
      !is.na(diabetes_date) &
      diabetes_date > index_date &
      diabetes_date <= index_date + years(5)
  ) %>%
  group_by(Id) %>%
  summarise(
    index_date = first(index_date),
    trial_eligible = first(trial_eligible),
    pred_at_index = first(pred_at_index),
    age_index = first(age_index),
    gender = first(GENDER),
    n_conditions = first(n_conditions),
    n_medications = first(n_medications),
    n_emergency = first(n_emergency),
    n_inpatient = first(n_inpatient),
    diabetes_5y = any(diabetes_5y, na.rm = TRUE),
    .groups = "drop"
  )

# Public identifiers are sequential synthetic IDs. Original Synthea UUIDs are
# not required for the analysis distributed on the website.
patient_map <- patients %>%
  transmute(
    Id,
    patient_id = sprintf("P%05d", row_number())
  )

patient_features_public <- patient_features %>%
  left_join(patient_map, by = "Id") %>%
  transmute(
    patient_id,
    birth_date = as.Date(BIRTHDATE),
    gender = GENDER,
    age_at_analysis_date = age,
    n_conditions,
    n_medications,
    n_encounters,
    n_emergency,
    n_inpatient,
    followup_years,
    hypertension,
    hyperlipidemia,
    coronary_heart_disease,
    stroke_history,
    myocardial_infarction_history,
    heart_failure,
    chronic_kidney_disease,
    cancer_history,
    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
  )

htn_outcome_public <- htn_outcome %>%
  left_join(patient_map, by = "Id") %>%
  transmute(
    patient_id,
    index_date,
    trial_eligible,
    pred_at_index,
    age_index,
    gender,
    n_conditions,
    n_medications,
    n_emergency,
    n_inpatient,
    diabetes_5y
  )

write_csv(
  patient_features_public,
  file.path(DATA_DIR, "rwe_patient_features.csv")
)

write_csv(
  htn_outcome_public,
  file.path(DATA_DIR, "rwe_htn_outcome.csv")
)

cat("\nData preparation complete.\n")
cat("Analysis date: ", format(analysis_date), "\n", sep = "")
cat("Patients: ", nrow(patient_features_public), "\n", sep = "")
cat(
  "Trial eligible: ",
  sum(patient_features_public$trial_eligible),
  "\n",
  sep = ""
)
cat("HTN outcome cohort: ", nrow(htn_outcome_public), "\n", sep = "")
