# File: S6_Statistical_Analysis.R
# Study: Diagnostic Accuracy of Haematological Indices in Yemen
# Author:  Naif Taleb Ali 
# Date: 2026
# Software: R version 4.3.1

# Load required packages
library(tidyverse)
library(pROC)
library(caret)
library(ggplot2)
library(forestplot)
library(gtsummary)

# 1. DATA IMPORT AND CLEANING
# Note: Ensure 'yemen_hematological_data.csv' is in your working directory
if (file.exists("yemen_hematological_data.csv")) {
  data <- read_csv("yemen_hematological_data.csv") %>%
    mutate(
      infection_status = as.factor(infection_status),
      sex = as.factor(sex),
      governorate = as.factor(governorate),
      fever = as.factor(fever),
      prior_abx = as.factor(prior_abx),
      # Calculate derived indices
      nlr = neutrophils / lymphocytes,
      plr = platelets / lymphocytes,
      # Create anemia flag
      anemia = ifelse(hemoglobin < 12, 1, 0),
      # Age groups
      age_group = cut(age_years, 
                      breaks = c(5, 18, 50, 70),
                      labels = c("<18", "18-50", ">50"))
    )

  # 2. DESCRIPTIVE STATISTICS
  # Table 1: Baseline characteristics
  table1 <- tbl_summary(
    data,
    by = infection_status,
    include = c(age_years, sex, governorate, fever, prior_abx),
    statistic = list(all_continuous() ~ "{mean} ± {sd}",
                     all_categorical() ~ "{n} ({p}%)"),
    digits = all_continuous() ~ 1
  )

  # Table 2: Hematological indices
  table2 <- data %>%
    select(infection_status, neutrophils, lymphocytes, platelets,
           hemoglobin, rdw, nlr, plr) %>%
    tbl_summary(
      by = infection_status,
      statistic = list(all_continuous() ~ "{median} ({p25}, {p75})"),
      digits = all_continuous() ~ 2
    ) %>%
    add_p(test = list(all_continuous() ~ "wilcox.test"))

  # 3. ROC CURVE ANALYSIS
  # Calculate AUC for each index
  roc_nlr <- roc(data$infection_status, data$nlr)
  roc_plr <- roc(data$infection_status, data$plr)
  roc_rdw <- roc(data$infection_status, data$rdw)
  roc_hb <- roc(data$infection_status, data$hemoglobin)

  # Get optimal cut-off using Youden Index
  optimal_nlr <- coords(roc_nlr, "best", best.method = "youden")
  optimal_plr <- coords(roc_plr, "best", best.method = "youden")
  optimal_rdw <- coords(roc_rdw, "best", best.method = "youden")
  optimal_hb <- coords(roc_hb, "best", best.method = "youden")

  # Create ROC curve plot
  roc_data <- bind_rows(
    tibble(
      sensitivity = roc_nlr$sensitivities,
      specificity = roc_nlr$specificities,
      index = "NLR",
      auc = auc(roc_nlr)
    ),
    tibble(
      sensitivity = roc_plr$sensitivities,
      specificity = roc_plr$specificities,
      index = "PLR",
      auc = auc(roc_plr)
    ),
    tibble(
      sensitivity = roc_rdw$sensitivities,
      specificity = roc_rdw$specificities,
      index = "RDW",
      auc = auc(roc_rdw)
    ),
    tibble(
      sensitivity = roc_hb$sensitivities,
      specificity = roc_hb$specificities,
      index = "Hemoglobin",
      auc = auc(roc_hb)
    )
  )

  # 4. MULTIVARIABLE LOGISTIC REGRESSION
  model <- glm(infection_status ~ nlr + plr + rdw + hemoglobin + 
                 age_years + sex + prior_abx,
               data = data, family = "binomial")

  # Calculate odds ratios and confidence intervals
  or_ci <- exp(cbind(OR = coef(model), confint(model)))

  # 5. SUBGROUP ANALYSIS
  # By governorate
  roc_nlr_aldhalea <- roc(data$infection_status[data$governorate == "Al-Dhalea"],
                          data$nlr[data$governorate == "Al-Dhalea"])
  roc_nlr_lahj <- roc(data$infection_status[data$governorate == "Lahj"],
                      data$nlr[data$governorate == "Lahj"])

  # Compare AUCs using DeLong test
  delong_test <- roc.test(roc_nlr_aldhalea, roc_nlr_lahj)

  # 6. SENSITIVITY ANALYSIS (Microscopy only)
  data_microscopy <- data %>% filter(!is.na(microscopy))
  if (nrow(data_microscopy) > 0) {
    roc_nlr_micro <- roc(data_microscopy$microscopy, data_microscopy$nlr)
  }

  # 7. CREATE SUPPLEMENTARY FIGURES
  # Figure S1: Distribution plots
  p1 <- ggplot(data, aes(x = infection_status, y = nlr, fill = infection_status)) +
    geom_boxplot() +
    scale_y_log10() +
    labs(title = "NLR Distribution by Infection Status")

  # Figure S2: Forest plot
  forest_data <- data.frame(
    variable = c("NLR", "PLR", "RDW", "Hemoglobin", "Age", "Sex (Female)", "Prior Antibiotic"),
    or = c(5.60, 3.40, 1.02, 0.94, 1.00, 0.92, 1.00),
    lower = c(3.87, 2.48, 0.95, 0.84, 0.99, 0.66, 0.66),
    upper = c(8.11, 4.66, 1.10, 1.05, 1.01, 1.28, 1.52)
  )

  # 8. EXPORT RESULTS
  write_csv(roc_data, "roc_curve_data.csv")
  write_csv(forest_data, "forest_plot_data.csv")
  saveRDS(model, "logistic_model.rds")

  # Generate report (if .Rmd exists)
  if (file.exists("analysis_report.Rmd")) {
    rmarkdown::render("analysis_report.Rmd")
  }

  print("Analysis complete. Check output files.")
} else {
  print("Error: 'yemen_hematological_data.csv' not found. Please provide the data file.")
}