# ==============================================================================
# SCRIPT R: Calculadora de Aproximación HR/IRR a Medidas de impacto
# Autor: E. Ortega. Grupo de Trabajo Pediatría Basada en la Evidencia (GTPBE)
# ==============================================================================

# 1. CARGA DE LIBRERÍAS
# ------------------------------------------------------------------------------
if (!require("rmarkdown")) install.packages("rmarkdown", quiet = TRUE)
if (!require("knitr")) install.packages("knitr", quiet = TRUE)
if (!require("ggplot2")) install.packages("ggplot2", quiet = TRUE)
if (!require("dplyr")) install.packages("dplyr", quiet = TRUE)
if (!require("tidyr")) install.packages("tidyr", quiet = TRUE)

library(rmarkdown)
library(knitr)
library(ggplot2)
library(dplyr)
library(tidyr)

# 2. CONFIGURACIÓN DE PARÁMETROS (INPUTS)
# ------------------------------------------------------------------------------
input_Risk_Pt   <- 1.91       # < 1: Prevenible | > 1: Atribuible
input_Risk_Low  <- 1.40       
input_Risk_Up   <- 2.59       
label_Risk      <- "HR"       

input_Pe        <- 0.0159     # Prevalencia Exposición (Población)
input_N_total   <- 22393      

# --- Configuración Método B (Flexible) ---
input_Pce_obs   <- 0.0282     # Proporción de Casos Expuestos (Observada)
input_N_cases   <- NA         # NA = Pce fija (sin N muestral).
# Número = Pce simulada (con N muestral).

input_n_boot    <- 10000      # Iteraciones Bootstrap

# 3. MOTOR DE CÁLCULO (BOOTSTRAP)
# ------------------------------------------------------------------------------
ejecutar_bootstrap <- function() {
  cat("\nIniciando simulación Bootstrap (n =", input_n_boot, " iteraciones)...\n")
  set.seed(123)
  is_prev <- input_Risk_Pt < 1
  
  # A. Riesgo (Log-Normal)
  log_risk_se <- (log(input_Risk_Up) - log(input_Risk_Low)) / (2 * 1.96)
  sim_risk    <- exp(rnorm(input_n_boot, log(input_Risk_Pt), log_risk_se))
  
  # B. Prevalencia (Beta)
  sim_pe <- rbeta(input_n_boot, round(input_Pe * input_N_total) + 1, 
                  input_N_total - round(input_Pe * input_N_total) + 1)
  
  # C. Pce para Método B
  if (!is.na(input_Pce_obs)) {
    if (!is.na(input_N_cases)) {
      # Con Incertidumbre (Beta)
      sim_pce_b <- rbeta(input_n_boot, round(input_Pce_obs * input_N_cases) + 1, 
                         input_N_cases - round(input_Pce_obs * input_N_cases) + 1)
    } else {
      # Sin Incertidumbre (Fija)
      sim_pce_b <- rep(input_Pce_obs, input_n_boot)
    }
  } else {
    sim_pce_b <- rep(NA, input_n_boot)
  }
  
  if (!is_prev) {
    # --- RIESGO (>1) -> ATRIBUIBLE ---
    sim_fae <- (sim_risk - 1) / sim_risk
    
    # Met A (Levin)
    sim_fap_A <- (sim_pe * (sim_risk - 1)) / (1 + sim_pe * (sim_risk - 1))
    
    # Met B (Miettinen: Pce * FAe)
    sim_fap_B <- sim_pce_b * sim_fae
    
    # Met C (Teórico)
    pce_t <- (sim_pe * sim_risk) / (1 + sim_pe * (sim_risk - 1))
    sim_fap_C <- pce_t * sim_fae
    
  } else {
    # --- PROTECTOR (<1) -> PREVENIBLE ---
    sim_fae <- 1 - sim_risk
    
    # Met A
    sim_fap_A <- sim_pe * (1 - sim_risk)
    
    # Met B (Fórmula ajustada para Preventiva)
    sim_fap_B <- (sim_pce_b * (1 - sim_risk)) / (sim_risk + sim_pce_b * (1 - sim_risk))
    
    # Met C
    pce_t <- (sim_pe * sim_risk) / (sim_pe * sim_risk + (1 - sim_pe))
    sim_fap_C <- (pce_t * (1 - sim_risk)) / (sim_risk + pce_t * (1 - sim_risk))
  }
  
  return(list(df = data.frame(fae=sim_fae, fap_A=sim_fap_A, fap_B=sim_fap_B, fap_C=sim_fap_C), 
              is_prev = is_prev))
}

# 4. GENERACIÓN DE RESULTADOS Y REPORTE
# ------------------------------------------------------------------------------
generar_reporte <- function() {
  res_boot <- ejecutar_bootstrap()
  df <- res_boot$df; is_p <- res_boot$is_prev
  
  etiq_i <- if(is_p) "FPe (Individual)" else "FAe (Individual)"
  etiq_p <- if(is_p) "FPp (Poblacional)" else "FAP (Poblacional)"
  
  # --- TABLA CONSOLA ---
  calc_stats <- function(x, lab) {
    if(all(is.na(x))) return(NULL)
    data.frame(Medida = lab, 
               Media = round(mean(x)*100, 2), 
               IC_2.5 = round(quantile(x, 0.025)*100, 2), 
               IC_97.5 = round(quantile(x, 0.975)*100, 2))
  }
  
  res_tab_consola <- rbind(
    calc_stats(df$fae, etiq_i),
    calc_stats(df$fap_A, paste(etiq_p, "Met. A")),
    calc_stats(df$fap_B, paste(etiq_p, "Met. B")),
    calc_stats(df$fap_C, paste(etiq_p, "Met. C"))
  )
  
  cat("\n================================================================\n")
  cat(" INFORME MEDIDAS IMPACTO - Autor: E. Ortega. GTPBE \n")
  cat("================================================================\n")
  print(kable(res_tab_consola, format = "simple"))
  cat("================================================================\n")
  cat("DETALLE DE MÉTODOS:\n")
  cat("1. Individual: Impacto en expuestos.\n")
  cat("2. Método A: Vía Prevalencia Poblacional (Pe).\n")
  cat(paste0("3. Método B: Vía Pce Observada (", if(is.na(input_N_cases)) "Fija" else paste0("Variable N=", input_N_cases), ").\n"))
  cat("4. Método C: Vía Pce Teórica (Derivada).\n")
  cat("================================================================\n\n")
  
  # --- GENERAR PDF ---
  rmd_file <- "temp_reporte.Rmd"
  tipo_doc <- if(is_p) "PREVENIBLE" else "ATRIBUIBLE"
  pdf_name <- paste0("Informe_Impacto_", tipo_doc, ".pdf")
  
  # Función auxiliar para tabla PDF
  calc_row_pdf <- function(x, nombre) {
    if(all(is.na(x))) return(NULL)
    media <- mean(x)*100; low <- quantile(x, 0.025)*100; up <- quantile(x, 0.975)*100
    data.frame(Concepto = nombre, 
               Estimacion = sprintf("%.2f%% (IC: %.2f%% - %.2f%%)", media, low, up))
  }
  
  # Creamos la tabla del PDF aquí para pasarla como objeto limpio
  tab_pdf_clean <- rbind(
    calc_row_pdf(df$fae, paste("1.", etiq_i, "en expuestos")),
    calc_row_pdf(df$fap_A, paste("2.", etiq_p, "- Método A (Pe)")),
    calc_row_pdf(df$fap_B, paste("3.", etiq_p, "- Método B (Pce Observada)")),
    calc_row_pdf(df$fap_C, paste("4.", etiq_p, "- Método C (Pce Teórica)"))
  )
  
  rmd_content <- c(
    "---",
    paste0("title: 'Informe Metodológico: Medidas de Impacto (", tipo_doc, ")'"),
    "author: 'Autor: E. Ortega. GTPBE.'",
    paste0("date: '", Sys.Date(), "'"),
    "output: pdf_document",
    "params:",
    "  tab_pdf: NA",
    "  df: NA",
    "---",
    "",
    "## 1. Parámetros del Análisis",
    paste0("* **Medida de Asociación (", label_Risk, "):** ", input_Risk_Pt, " (IC 95%: ", input_Risk_Low, "-", input_Risk_Up, ")"),
    paste0("* **Interpretación:** ", ifelse(is_p, "Factor Protector", "Factor de Riesgo")),
    paste0("* **Prevalencia de Exposición (Pe):** ", input_Pe * 100, "%"),
    if(!is.na(input_Pce_obs)) paste0("* **Pce (Casos):** ", input_Pce_obs * 100, "%") else "",
    "",
    "## 2. Resultados de la Simulación (IC 95%)",
    "```{r, echo=FALSE}",
    "knitr::kable(params$tab_pdf, caption = 'Table 1: Estimaciones Bootstrap N=10000')",
    "```",
    "",
    "## 3. Visualización de Distribuciones Poblacionales",
    "Se muestra la distribución de probabilidad para los tres métodos.",
    "",
    "```{r, echo=FALSE, fig.height=5, warning=FALSE, message=FALSE}",
    "library(ggplot2); library(tidyr); library(dplyr)",
    "# Forzamos la selección de columnas para asegurar que Metodo B salga",
    "df_long <- params$df %>% ",
    "  select(fap_A, fap_B, fap_C) %>% ",
    "  pivot_longer(cols=everything(), names_to='Metodo', values_to='Valor') %>%",
    "  filter(!is.na(Valor)) %>%",
    "  mutate(Metodo = recode(Metodo, ",
    "    'fap_A' = 'A. Vía Pe',",
    "    'fap_B' = 'B. Vía Pce (Obs)',",
    "    'fap_C' = 'C. Vía Pce (Teo)'))",
    "",
    "ggplot(df_long, aes(x=Valor*100, fill=Metodo)) +",
    "  geom_density(alpha=0.4, color='black') +",
    "  scale_fill_manual(values=c('A. Vía Pe'='#3498db', 'B. Vía Pce (Obs)'='#e67e22', 'C. Vía Pce (Teo)'='#2ecc71')) +",
    "  labs(x='Impacto Poblacional (%)', y='Densidad', title='Distribución Posterior Comparativa') +",
    "  theme_minimal() + theme(legend.position='bottom')",
    "```",
    "",
    "## 4. Fórmulas Aplicadas",
    "",
    paste0("En este reporte se han aplicado las fórmulas para factores de **", ifelse(is_p, "protección", "riesgo"), "**."),
    "",
    if(is_p) {
      paste0(
        "1. **FPe (Individual):** $$ 1 - ", label_Risk, " $$ \n\n",
        "2. **FPp (Método A):** $$ Pe \\times (1 - ", label_Risk, ") $$ \n\n",
        "3. **FPp (Método B):** $$ \\frac{Pce(1 - ", label_Risk, ")}{", label_Risk, " + Pce(1 - ", label_Risk, ")} $$"
      )
    } else {
      paste0(
        "1. **FAe (Individual):** $$ \\frac{", label_Risk, " - 1}{", label_Risk, "} $$ \n\n",
        "2. **FAP (Método A):** $$ \\frac{Pe(", label_Risk, " - 1)}{1 + Pe(", label_Risk, " - 1)} $$ \n\n",
        "3. **FAP (Método B):** $$ Pce \\times FAe $$"
      )
    }
  )
  
  writeLines(rmd_content, rmd_file)
  rmarkdown::render(rmd_file, output_file = pdf_name, params = list(tab_pdf = tab_pdf_clean, df = df), quiet = TRUE)
  if(file.exists(rmd_file)) file.remove(rmd_file)
  
  cat("[ÉXITO] Informe PDF generado como:", pdf_name, "\n")
}

generar_reporte()