EpiSimple

# ============================================================================
# ANÁLISIS DE PREVALENCIA DE LESIÓN RENAL AGUDA EN SEPSIS - UCI
# Hospital Federico Lleras Acosta, Ibagué, Colombia
# 
# Autor: [Tu nombre]
# Fecha: [Fecha]
# Versión: 1.0
# 
# Este código realiza:
# 1. Importación interactiva de datos (CSV, Excel, SPSS, Stata)
# 2. Limpieza y validación de datos
# 3. Análisis descriptivo completo
# 4. Estimación de prevalencias con IC95%
# 5. Análisis bivariado exploratorio
# 6. Gráficos listos para publicación (estilo The Lancet)
# 7. Tablas formato APA/Vancouver
# ============================================================================

# ============================================================================
# SECCIÓN 1: INSTALACIÓN Y CARGA DE PAQUETES
# ============================================================================

# Lista de paquetes necesarios
paquetes_necesarios <- c(
  "tidyverse",      # Manipulación de datos (dplyr, ggplot2, tidyr, etc.)
  "readxl",         # Leer archivos Excel (.xlsx, .xls)
  "haven",          # Leer archivos SPSS (.sav) y Stata (.dta)
  "janitor",        # Limpieza de nombres de variables
  "tableone",       # Tabla 1 de características (clásica en artículos médicos)
  "gtsummary",      # Tablas elegantes listas para publicación
  "flextable",      # Tablas exportables a Word/PowerPoint
  "epiR",           # Funciones epidemiológicas (prevalencias, IC95%)
  "DescTools",      # Herramientas descriptivas adicionales
  "ggplot2",        # Gráficos de alta calidad
  "ggpubr",         # Publicación: p-values, comparaciones en gráficos
  "patchwork",      # Combinar múltiples gráficos
  "scales",         # Formateo de ejes (porcentajes, números)
  "RColorBrewer",   # Paletas de colores profesionales
  "viridis",        # Paletas accesibles para daltónicos
  "naniar",         # Análisis de datos faltantes
  "skimr",          # Resumen rápido y completo de datos
  "here"            # Manejo de rutas de archivos
)

# Función para instalar paquetes faltantes
instalar_si_falta <- function(paquete) {
  if (!require(paquete, character.only = TRUE)) {
    install.packages(paquete, dependencies = TRUE)
    library(paquete, character.only = TRUE)
  }
}

# Instalar y cargar todos los paquetes
cat("Instalando/cargando paquetes necesarios...\n")
invisible(lapply(paquetes_necesarios, instalar_si_falta))
cat("✓ Todos los paquetes cargados exitosamente.\n\n")


# ============================================================================
# SECCIÓN 2: IMPORTACIÓN INTERACTIVA DE DATOS
# ============================================================================

# Función para importar datos de múltiples formatos
importar_datos <- function() {
  cat("\n========================================\n")
  cat("IMPORTACIÓN DE DATOS\n")
  cat("========================================\n\n")
  
  # Abrir cuadro de diálogo para seleccionar archivo
  cat("Por favor, seleccione su archivo de datos...\n")
  ruta_archivo <- file.choose()
  
  # Detectar extensión del archivo
  extension <- tolower(tools::file_ext(ruta_archivo))
  
  # Importar según el tipo de archivo
  datos <- switch(extension,
    "csv"  = read.csv(ruta_archivo, stringsAsFactors = FALSE),
    "xlsx" = read_excel(ruta_archivo),
    "xls"  = read_excel(ruta_archivo),
    "sav"  = read_sav(ruta_archivo),
    "dta"  = read_dta(ruta_archivo),
    stop("Formato no soportado. Use: .csv, .xlsx, .xls, .sav, o .dta")
  )
  
  cat(sprintf("✓ Archivo importado: %s\n", basename(ruta_archivo)))
  cat(sprintf("✓ Dimensiones: %d filas × %d columnas\n\n", 
              nrow(datos), ncol(datos)))
  
  return(datos)
}

# Ejecutar importación
datos_raw <- importar_datos()


# ============================================================================
# SECCIÓN 3: LIMPIEZA Y VALIDACIÓN DE DATOS
# ============================================================================

cat("\n========================================\n")
cat("LIMPIEZA DE DATOS\n")
cat("========================================\n\n")

# 3.1 Limpiar nombres de variables (sin espacios, minúsculas, sin caracteres especiales)
cat("Paso 1: Estandarizando nombres de variables...\n")
datos <- datos_raw %>%
  clean_names()  # Convierte a snake_case: edad_anos, apache_ii, etc.

cat("✓ Nombres estandarizados a formato snake_case\n")
cat("Nombres de variables:\n")
print(names(datos))

# 3.2 Inspección inicial de datos
cat("\n\nPaso 2: Inspección inicial de datos...\n")
cat("Resumen de la estructura de datos:\n")
glimpse(datos)

# 3.3 Análisis de datos faltantes
cat("\n\nPaso 3: Análisis de datos faltantes...\n")
datos_faltantes <- datos %>%
  summarise(across(everything(), 
                      ~sum(is.na(.)) / length(.) * 100)) %>%
  pivot_longer(everything(), 
               names_to = "variable", 
               values_to = "porcentaje_faltante") %>%
  filter(porcentaje_faltante > 0) %>%
  arrange(desc(porcentaje_faltante))

if (nrow(datos_faltantes) > 0) {
  cat("\n⚠ Variables con datos faltantes:\n")
  print(datos_faltantes, n = Inf)
} else {
  cat("✓ No hay datos faltantes\n")
}

# Gráfico de patrón de datos faltantes (si hay alguno)
if (any(is.na(datos))) {
  cat("\nGenerando gráfico de datos faltantes...\n")
  gg_miss_var(datos) +
    labs(title = "Missing Data Pattern",
         subtitle = sprintf("Dataset: %d observations", nrow(datos)),
         x = "Variable",
         y = "Number of Missing Values") +
    theme_minimal(base_size = 12) +
    theme(axis.text.x = element_text(angle = 45, hjust = 1))
}

# 3.4 Conversión de tipos de variables y etiquetado
cat("\n\nPaso 4: Conversión de variables categóricas a factores...\n")

# IMPORTANTE: Adaptar estos nombres a los de TU base de datos
# Esta es una plantilla - debes modificar según tus variables
datos <- datos %>%
  mutate(
    # Variables demográficas
    sexo = factor(sexo, 
                   levels = c(1, 2), 
                   labels = c("Male", "Female")),
    
    # Variable de desenlace principal
    lra = factor(lra, 
                 levels = c(0, 1), 
                 labels = c("No AKI", "AKI")),
    
    estadio_lra = factor(estadio_lra,
                          levels = c(0, 1, 2, 3),
                          labels = c("No AKI", "Stage 1", "Stage 2", "Stage 3"),
                          ordered = TRUE),
    
    # Variables clínicas
    choque_septico = factor(choque_septico,
                               levels = c(0, 1),
                               labels = c("No Septic Shock", "Septic Shock")),
    
    tipo_sepsis = factor(tipo_sepsis,
                          levels = c(1, 2, 3, 4, 5),
                          labels = c("Pulmonary", "Abdominal", "Urinary", 
                                    "Skin/Soft Tissue", "Primary Bacteremia")),
    
    # Comorbilidades
    diabetes = factor(diabetes,
                       levels = c(0, 1),
                       labels = c("No Diabetes", "Diabetes")),
    
    hipertension = factor(hipertension,
                           levels = c(0, 1),
                           labels = c("No Hypertension", "Hypertension")),
    
    # Desenlaces
    mortalidad_uci = factor(mortalidad_uci,
                               levels = c(0, 1),
                               labels = c("Survivor", "Non-survivor")),
    
    # Asegurar que variables continuas sean numéricas
    edad = as.numeric(edad),
    apache_ii = as.numeric(apache_ii),
    sofa = as.numeric(sofa),
    creatinina_basal = as.numeric(creatinina_basal),
    creatinina_max = as.numeric(creatinina_max),
    lactato = as.numeric(lactato)
  )

cat("✓ Variables categóricas convertidas a factores con etiquetas\n")
cat("✓ Variables continuas verificadas como numéricas\n")

# 3.5 Detección de valores atípicos (outliers) en variables continuas
cat("\n\nPaso 5: Detección de valores atípicos...\n")

detectar_outliers <- function(datos, variable) {
  q1 <- quantile(datos[[variable]], 0.25, na.rm = TRUE)
  q3 <- quantile(datos[[variable]], 0.75, na.rm = TRUE)
  iqr <- q3 - q1
  limite_inferior <- q1 - 3 * iqr
  limite_superior <- q3 + 3 * iqr
  
  outliers <- datos[[variable]] < limite_inferior | datos[[variable]] > limite_superior
  n_outliers <- sum(outliers, na.rm = TRUE)
  
  if (n_outliers > 0) {
    cat(sprintf("  %s: %d outliers extremos (%.1f%%)\n", 
                variable, n_outliers, n_outliers/length(datos[[variable]])*100))
  }
}

variables_continuas <- c("edad", "apache_ii", "sofa", "creatinina_basal", 
                           "creatinina_max", "lactato")
invisible(lapply(variables_continuas, function(x) detectar_outliers(datos, x)))

cat("\n✓ Limpieza de datos completada\n")


# ============================================================================
# SECCIÓN 4: ANÁLISIS DESCRIPTIVO
# ============================================================================

cat("\n\n========================================\n")
cat("ANÁLISIS DESCRIPTIVO\n")
cat("========================================\n\n")

# 4.1 Resumen general de datos
cat("Resumen estadístico de todas las variables:\n\n")
skim(datos)

# 4.2 Evaluación de normalidad para variables continuas
cat("\n\nEvaluación de normalidad (Test de Shapiro-Wilk):\n")
cat("(Valores p < 0.05 indican distribución NO normal)\n\n")

test_normalidad <- sapply(variables_continuas, function(var) {
  if (sum(!is.na(datos[[var]])) < 5000) {
    test <- shapiro.test(datos[[var]][!is.na(datos[[var]])])
    c(Statistic = test$statistic, P_value = test$p.value)
  } else {
    c(Statistic = NA, P_value = NA)  # Shapiro no funciona bien con n>5000
  }
})

print(t(test_normalidad))

# 4.3 Crear Tabla 1 de características basales (estilo publicación)
cat("\n\nGenerando Tabla 1 de características...\n")

variables_tabla1 <- c("edad", "sexo", "apache_ii", "sofa", 
                        "choque_septico", "tipo_sepsis", "diabetes", 
                        "hipertension", "creatinina_basal", "lactato")

tabla1 <- datos %>%
  select(all_of(variables_tabla1)) %>%
  tbl_summary(
    statistic = list(
      all_continuous() ~ "{median} ({p25}, {p75})",  # Mediana (RIQ) para continuas
      all_categorical() ~ "{n} ({p}%)"              # n (%) para categóricas
    ),
    digits = list(
      all_continuous() ~ 1,
      all_categorical() ~ c(0, 1)
    ),
    label = list(
      edad ~ "Age, years",
      sexo ~ "Sex",
      apache_ii ~ "APACHE II score",
      sofa ~ "SOFA score",
      choque_septico ~ "Septic shock",
      tipo_sepsis ~ "Source of sepsis",
      diabetes ~ "Diabetes mellitus",
      hipertension ~ "Hypertension",
      creatinina_basal ~ "Baseline creatinine, mg/dL",
      lactato ~ "Lactate, mmol/L"
    ),
    missing_text = "Missing data"
  ) %>%
  modify_header(label = "**Variable**") %>%
  modify_caption("Table 1. Baseline Characteristics") %>%
  bold_labels()

print(tabla1)

# Exportar Tabla 1 a Word (opcional)
# tabla1 %>% as_flex_table() %>% save_as_docx(path = "Tabla1.docx")


# ============================================================================
# SECCIÓN 5: ESTIMACIÓN DE PREVALENCIA
# ============================================================================

cat("\n\n========================================\n")
cat("ESTIMACIÓN DE PREVALENCIA DE LRA\n")
cat("========================================\n\n")

# 5.1 Prevalencia global con IC95% (método de Wilson)
n_total <- sum(!is.na(datos$lra))
n_lra <- sum(datos$lra == "AKI", na.rm = TRUE)

prev_global <- epi.conf(
  dat = matrix(c(n_lra, n_total - n_lra), nrow = 1),
  ctype = "prevalence",
  method = "wilson",
  conf.level = 0.95
)

cat("PREVALENCIA GLOBAL DE LESIÓN RENAL AGUDA:\n")
cat(sprintf("  Casos: %d / %d\n", n_lra, n_total))
cat(sprintf("  Prevalencia: %.1f%% (95%% CI: %.1f%% - %.1f%%)\n\n",
            prev_global$est * 100,
            prev_global$lower * 100,
            prev_global$upper * 100))

cat("📝 CÓMO REDACTAR EN EL ARTÍCULO:\n")
cat(sprintf("'The prevalence of AKI in patients with sepsis was %.1f%% (95%% CI %.1f%%-%.1f%%, n=%d/%d).'\n\n",
            prev_global$est * 100, prev_global$lower * 100, 
            prev_global$upper * 100, n_lra, n_total))

# 5.2 Prevalencia por estadio KDIGO
cat("\nPREVALENCIA POR ESTADIO KDIGO:\n")

prev_kdigo <- datos %>%
  count(estadio_lra) %>%
  mutate(
    prevalencia = n / sum(n) * 100,
    # Calcular IC95% para cada estadio
    ic_lower = mapply(function(casos, total) {
      epi.conf(matrix(c(casos, total - casos), nrow = 1),
               ctype = "prevalence", method = "wilson")$lower * 100
    }, n, sum(n)),
    ic_upper = mapply(function(casos, total) {
      epi.conf(matrix(c(casos, total - casos), nrow = 1),
               ctype = "prevalence", method = "wilson")$upper * 100
    }, n, sum(n))
  )

print(prev_kdigo)


# ============================================================================
# SECCIÓN 6: GRÁFICOS LISTOS PARA PUBLICACIÓN
# ============================================================================

cat("\n\n========================================\n")
cat("GENERANDO GRÁFICOS (ESTILO THE LANCET)\n")
cat("========================================\n\n")

# Tema personalizado para publicación
tema_publicacion <- theme_minimal(base_size = 14, base_family = "sans") +
  theme(
    plot.title = element_text(size = 16, face = "bold", hjust = 0),
    plot.subtitle = element_text(size = 12, color = "grey30", hjust = 0),
    axis.title = element_text(size = 12, face = "bold"),
    axis.text = element_text(size = 11, color = "black"),
    legend.title = element_text(size = 12, face = "bold"),
    legend.text = element_text(size = 11),
    legend.position = "bottom",
    panel.grid.minor = element_blank(),
    panel.grid.major = element_line(color = "grey90", size = 0.3),
    plot.caption = element_text(size = 9, color = "grey50", hjust = 0)
  )

# Paleta de colores profesional (The Lancet style)
colores_lancet <- c("#00468B", "#ED0000", "#42B540", "#0099B4", 
                     "#925E9F", "#FDAF91", "#AD002A")

# 6.1 Gráfico de barras: Prevalencia por estadio KDIGO
cat("Generando Figure 1: Prevalencia por estadio KDIGO...\n")

fig1_prevalencia_kdigo <- ggplot(prev_kdigo, aes(x = estadio_lra, y = prevalencia, 
                                                fill = estadio_lra)) +
  geom_col(alpha = 0.9, width = 0.7) +
  geom_errorbar(aes(ymin = ic_lower, ymax = ic_upper), 
                width = 0.2, size = 0.8) +
  geom_text(aes(label = sprintf("%.1f%%\n(n=%d)", prevalencia, n)),
            vjust = -0.5, size = 4, fontface = "bold") +
  scale_fill_manual(values = colores_lancet) +
  scale_y_continuous(limits = c(0, max(prev_kdigo$ic_upper) * 1.2),
                     expand = expansion(mult = c(0, 0.05))) +
  labs(
    title = "Figure 1. Prevalence of Acute Kidney Injury by KDIGO Stage",
    subtitle = sprintf("Patients with sepsis admitted to ICU (N=%d)", n_total),
    x = "KDIGO Stage",
    y = "Prevalence (%)",
    caption = "Error bars represent 95% confidence intervals (Wilson method).\nAKI = Acute Kidney Injury; KDIGO = Kidney Disease: Improving Global Outcomes"
  ) +
  tema_publicacion +
  theme(legend.position = "none")

print(fig1_prevalencia_kdigo)

# Guardar en alta resolución para publicación
ggsave("Figure1_AKI_Prevalence_by_KDIGO.png", 
       fig1_prevalencia_kdigo,
       width = 10, height = 7, dpi = 300, bg = "white")
ggsave("Figure1_AKI_Prevalence_by_KDIGO.pdf", 
       fig1_prevalencia_kdigo,
       width = 10, height = 7, device = "pdf")

cat("✓ Figura 1 guardada: Figure1_AKI_Prevalence_by_KDIGO.png (300 DPI)\n")
cat("✓ Figura 1 guardada: Figure1_AKI_Prevalence_by_KDIGO.pdf (vectorial)\n\n")

# ============================================================================
# SECCIÓN 7: GRÁFICOS ADICIONALES PARA PUBLICACIÓN
# ============================================================================

# 7.1 Distribución de edad por presencia de LRA (boxplot + violin)
cat("Generando Figure 2: Edad por estado de LRA...\n")

fig2_edad_lra <- ggplot(datos, aes(x = lra, y = edad, fill = lra)) +
  geom_violin(alpha = 0.5, trim = FALSE) +
  geom_boxplot(width = 0.2, alpha = 0.8, outlier.shape = 21, 
               outlier.size = 2, outlier.alpha = 0.5) +
  stat_summary(fun = mean, geom = "point", shape = 23, 
               size = 3, fill = "white", color = "black") +
  scale_fill_manual(values = c("No AKI" = colores_lancet[3], 
                                 "AKI" = colores_lancet[2])) +
  stat_compare_means(method = "wilcox.test", 
                      label = "p.format",
                      label.x = 1.5, label.y = max(datos$edad, na.rm = TRUE) * 0.95) +
  labs(
    title = "Figure 2. Age Distribution by AKI Status",
    subtitle = sprintf("N=%d patients with sepsis", n_total),
    x = "",
    y = "Age (years)",
    caption = "Box represents IQR with median line. Diamond represents mean.\nWilcoxon rank-sum test for group comparison."
  ) +
  tema_publicacion +
  theme(legend.position = "none")

print(fig2_edad_lra)
ggsave("Figure2_Age_by_AKI.png", fig2_edad_lra, 
       width = 8, height = 7, dpi = 300, bg = "white")


# 7.2 APACHE II y SOFA por presencia de LRA (comparación múltiple)
cat("Generando Figure 3: Scores de severidad por LRA...\n")

# Preparar datos en formato largo
datos_severidad <- datos %>%
  select(lra, apache_ii, sofa) %>%
  pivot_longer(cols = c(apache_ii, sofa),
               names_to = "score",
               values_to = "value") %>%
  mutate(score = factor(score, 
                          levels = c("apache_ii", "sofa"),
                          labels = c("APACHE II", "SOFA")))

fig3_scores <- ggplot(datos_severidad, aes(x = lra, y = value, fill = lra)) +
  geom_boxplot(alpha = 0.7, outlier.alpha = 0.3) +
  stat_summary(fun = mean, geom = "point", shape = 23, 
               size = 2.5, fill = "white") +
  facet_wrap(~score, scales = "free_y", ncol = 2) +
  scale_fill_manual(values = c("No AKI" = colores_lancet[3], 
                                 "AKI" = colores_lancet[2])) +
  stat_compare_means(method = "wilcox.test", 
                      label = "p.format",
                      label.y.npc = 0.95) +
  labs(
    title = "Figure 3. Illness Severity Scores by AKI Status",
    subtitle = "Comparison of APACHE II and SOFA scores",
    x = "",
    y = "Score",
    caption = "Diamond represents mean. P-values from Wilcoxon rank-sum test."
  ) +
  tema_publicacion +
  theme(legend.position = "none")

print(fig3_scores)
ggsave("Figure3_Severity_Scores.png", fig3_scores,
       width = 10, height = 6, dpi = 300, bg = "white")


# ============================================================================
# SECCIÓN 8: ANÁLISIS BIVARIADO COMPLETO
# ============================================================================

cat("\n\n========================================\n")
cat("ANÁLISIS BIVARIADO AVANZADO\n")
cat("========================================\n\n")

# 8.1 Tabla comparativa estratificada (Tabla 2 del artículo)
cat("Generando Tabla 2: Características según presencia de LRA...\n")

tabla2_bivariado <- datos %>%
  select(all_of(variables_tabla1), lra) %>%
  tbl_summary(
    by = lra,
    statistic = list(
      all_continuous() ~ "{median} ({p25}, {p75})",
      all_categorical() ~ "{n} ({p}%)"
    ),
    digits = list(
      all_continuous() ~ 1,
      all_categorical() ~ c(0, 1)
    ),
    label = list(
      edad ~ "Age, years",
      sexo ~ "Sex",
      apache_ii ~ "APACHE II score",
      sofa ~ "SOFA score",
      choque_septico ~ "Septic shock",
      tipo_sepsis ~ "Source of sepsis",
      diabetes ~ "Diabetes mellitus",
      hipertension ~ "Hypertension",
      creatinina_basal ~ "Baseline creatinine, mg/dL",
      lactato ~ "Lactate, mmol/L"
    ),
    missing_text = "Missing"
  ) %>%
  add_p(
    test = list(
      all_continuous() ~ "wilcox.test",  # Mann-Whitney U para continuas
      all_categorical() ~ "chisq.test"   # Chi-cuadrado para categóricas
    ),
    pvalue_fun = function(x) style_pvalue(x, digits = 3)
  ) %>%
  add_overall() %>%
  modify_header(
    label = "**Variable**",
    stat_0 = "**Overall**\nN = {N}",
    stat_1 = "**No AKI**\nN = {n}",
    stat_2 = "**AKI**\nN = {n}"
  ) %>%
  modify_caption("Table 2. Baseline Characteristics Stratified by AKI Status") %>%
  bold_labels() %>%
  italicize_levels()

print(tabla2_bivariado)

# Exportar a Word
tabla2_bivariado %>% 
  as_flex_table() %>% 
  save_as_docx(path = "Table2_Bivariate_Analysis.docx")

cat("✓ Tabla 2 exportada: Table2_Bivariate_Analysis.docx\n\n")


# 8.2 Calcular Razones de Prevalencia (RP) para factores de riesgo
cat("Calculando Razones de Prevalencia con IC95%...\n\n")

# Función mejorada para calcular RP
calcular_rp_completo <- function(datos, variable_exposicion, variable_desenlace = "lra") {
  
  # Crear tabla 2x2
  tabla <- table(datos[[variable_exposicion]], datos[[variable_desenlace]])
  
  # Verificar que sea 2x2
  if (nrow(tabla) != 2 | ncol(tabla) != 2) {
    return(list(error = "Variable no es dicotómica"))
  }
  
  # Calcular con epiR
  resultado <- epi.2by2(tabla, 
                         method = "cross.sectional",
                         conf.level = 0.95)
  
  # Extraer RP e IC95%
  rp <- resultado$massoc.detail$PR.strata.wald
  
  return(list(
    variable = variable_exposicion,
    rp = rp["est"],
    ic_lower = rp["lower"],
    ic_upper = rp["upper"],
    p_valor = resultado$massoc.detail$chi2.strata.uncor["p.value"]
  ))
}

# Calcular RP para variables categóricas dicotómicas
variables_categoricas_dicho <- c("sexo", "choque_septico", "diabetes", "hipertension")

resultados_rp <- lapply(variables_categoricas_dicho, function(var) {
  calcular_rp_completo(datos, var)
})

# Crear tabla de resultados RP
tabla_rp <- do.call(rbind, lapply(resultados_rp, function(x) {
  if (is.null(x$error)) {
    data.frame(
      Variable = x$variable,
      PR = sprintf("%.2f", x$rp),
      CI95 = sprintf("%.2f - %.2f", x$ic_lower, x$ic_upper),
      P_value = sprintf("%.3f", x$p_valor)
    )
  }
}))

cat("\nRAZONES DE PREVALENCIA (RP) CRUDAS:\n")
print(tabla_rp)

cat("\n📝 INTERPRETACIÓN:\n")
cat("- RP > 1: Mayor prevalencia de LRA en expuestos\n")
cat("- RP < 1: Menor prevalencia de LRA en expuestos\n")
cat("- IC95% que no incluye 1.0: Asociación estadísticamente significativa\n\n")

cat("📝 EJEMPLO DE REDACCIÓN:\n")
cat("'Patients with septic shock had a XX% higher prevalence of AKI\n")
cat("compared to those without septic shock (PR X.XX, 95% CI X.XX-X.XX, P=X.XXX).'\n\n")


# ============================================================================
# SECCIÓN 9: ANÁLISIS MULTIVARIADO - REGRESIÓN DE POISSON
# ============================================================================

cat("\n\n========================================\n")
cat("ANÁLISIS MULTIVARIADO\n")
cat("========================================\n\n")

cat("IMPORTANTE: Usaremos REGRESIÓN DE POISSON con varianza robusta\n")
cat("para estimar RP ajustados (NO regresión logística que daría OR).\n\n")

# 9.1 Modelo de Regresión de Poisson con varianza robusta
# Esto es MEJOR que regresión logística para desenlaces frecuentes (>10%)
# porque OR sobreestima el verdadero riesgo relativo

library(sandwich)  # Para errores estándar robustos
library(lmtest)    # Para coeftest

# Convertir LRA a numérico para Poisson
datos <- datos %>%
  mutate(lra_num = as.numeric(lra == "AKI"))

# Modelo completo (ajustar variables según tu estudio)
modelo_poisson <- glm(
  lra_num ~ edad + sexo + apache_ii + choque_septico + diabetes + 
            hipertension + creatinina_basal,
  data = datos,
  family = poisson(link = "log")
)

# Obtener RR ajustados con errores robustos
coef_robustos <- coeftest(modelo_poisson, vcov = vcovHC(modelo_poisson, type = "HC0"))

# Crear tabla de resultados
resultados_multivariado <- data.frame(
  Variable = rownames(coef_robustos),
  Coef = coef_robustos[, "Estimate"],
  SE = coef_robustos[, "Std. Error"],
  Z = coef_robustos[, "z value"],
  P = coef_robustos[, "Pr(>|z|)"]
) %>%
  mutate(
    PR = exp(Coef),
    IC_lower = exp(Coef - 1.96 * SE),
    IC_upper = exp(Coef + 1.96 * SE)
  ) %>%
  filter(Variable != "(Intercept)") %>%
  mutate(
    PR_formatted = sprintf("%.2f (%.2f-%.2f)", PR, IC_lower, IC_upper),
    P_formatted = ifelse(P < 0.001, "<0.001", sprintf("%.3f", P))
  )

cat("\nRESULTADOS DEL MODELO MULTIVARIADO:\n")
cat("Razones de Prevalencia (PR) Ajustadas con IC95%:\n\n")
print(resultados_multivariado[, c("Variable", "PR_formatted", "P_formatted")])

# Crear tabla publicable con gtsummary
modelo_tabla <- tbl_regression(
  modelo_poisson,
  exponentiate = TRUE,  # Convertir a PR (exp de coeficientes)
  label = list(
    edad ~ "Age (per year)",
    sexo ~ "Sex (Female vs Male)",
    apache_ii ~ "APACHE II (per point)",
    choque_septico ~ "Septic Shock",
    diabetes ~ "Diabetes Mellitus",
    hipertension ~ "Hypertension",
    creatinina_basal ~ "Baseline Creatinine (per mg/dL)"
  )
) %>%
  modify_header(
    estimate = "**PR**",
    ci = "**95% CI**",
    p.value = "**P-value**"
  ) %>%
  modify_caption("Table 3. Multivariable Analysis - Adjusted Prevalence Ratios for AKI") %>%
  bold_labels() %>%
  bold_p(t = 0.05)

print(modelo_tabla)

# Exportar
modelo_tabla %>%
  as_flex_table() %>%
  save_as_docx(path = "Table3_Multivariable_Model.docx")

cat("\n✓ Tabla 3 exportada: Table3_Multivariable_Model.docx\n\n")

cat("📝 CÓMO REDACTAR:\n")
cat("'In multivariable analysis adjusting for age, sex, illness severity,\n")
cat("and comorbidities, septic shock remained independently associated with\n")
cat("higher prevalence of AKI (adjusted PR X.XX, 95% CI X.XX-X.XX, P=X.XXX).'\n\n")


# 9.2 Evaluación del modelo
cat("Evaluación del ajuste del modelo:\n")

# Pseudo R-squared (McFadden)
modelo_nulo <- glm(lra_num ~ 1, data = datos, family = poisson())
pseudo_r2 <- 1 - (logLik(modelo_poisson) / logLik(modelo_nulo))
cat(sprintf("Pseudo R² (McFadden): %.3f\n", pseudo_r2))

# AIC
cat(sprintf("AIC: %.1f\n\n", AIC(modelo_poisson)))


# ============================================================================
# SECCIÓN 10: FOREST PLOT DE RAZONES DE PREVALENCIA
# ============================================================================

cat("\n\n========================================\n")
cat("GENERANDO FOREST PLOT\n")
cat("========================================\n\n")

# Preparar datos para forest plot
forest_data <- resultados_multivariado %>%
  mutate(
    Variable_label = case_when(
      Variable == "edad" ~ "Age (per year)",
      Variable == "sexoFemale" ~ "Female sex",
      Variable == "apache_ii" ~ "APACHE II (per point)",
      Variable == "choque_septicoSeptic Shock" ~ "Septic shock",
      Variable == "diabetesDiabetes" ~ "Diabetes mellitus",
      Variable == "hipertensionHypertension" ~ "Hypertension",
      Variable == "creatinina_basal" ~ "Baseline creatinine (per mg/dL)",
      TRUE ~ Variable
    )
  ) %>%
  arrange(desc(PR))

fig4_forest <- ggplot(forest_data, aes(x = PR, y = reorder(Variable_label, PR))) +
  geom_vline(xintercept = 1, linetype = "dashed", color = "gray50", size = 0.8) +
  geom_errorbarh(aes(xmin = IC_lower, xmax = IC_upper), 
                  height = 0.3, size = 0.8, color = colores_lancet[1]) +
  geom_point(size = 4, color = colores_lancet[1], shape = 18) +
  geom_text(aes(label = sprintf("%.2f (%.2f-%.2f)", PR, IC_lower, IC_upper)),
            hjust = -0.1, size = 3.5, fontface = "bold") +
  scale_x_continuous(
    trans = "log10",
    breaks = c(0.5, 0.75, 1, 1.5, 2, 3, 4),
    limits = c(0.4, max(forest_data$IC_upper) * 1.5)
  ) +
  labs(
    title = "Figure 4. Adjusted Prevalence Ratios for Acute Kidney Injury",
    subtitle = "Multivariable Poisson regression with robust variance",
    x = "Adjusted Prevalence Ratio (95% CI)",
    y = "",
    caption = "Vertical dashed line represents PR = 1.0 (no association).\nError bars represent 95% confidence intervals."
  ) +
  tema_publicacion +
  theme(
    axis.text.y = element_text(size = 11, hjust = 1),
    panel.grid.major.y = element_blank(),
    panel.grid.major.x = element_line(color = "grey90")
  )

print(fig4_forest)

ggsave("Figure4_Forest_Plot.png", fig4_forest,
       width = 12, height = 7, dpi = 300, bg = "white")
ggsave("Figure4_Forest_Plot.pdf", fig4_forest,
       width = 12, height = 7, device = "pdf")

cat("✓ Figure 4 guardada en PNG y PDF\n\n")


# ============================================================================
# SECCIÓN 11: RESUMEN FINAL Y EXPORTACIÓN COMPLETA
# ============================================================================

cat("\n\n========================================\n")
cat("RESUMEN DEL ANÁLISIS COMPLETO\n")
cat("========================================\n\n")

cat("ARCHIVOS GENERADOS:\n")
cat("  Figuras (PNG 300 DPI + PDF vectorial):\n")
cat("    ✓ Figure1_AKI_Prevalence_by_KDIGO.png/.pdf\n")
cat("    ✓ Figure2_Age_by_AKI.png\n")
cat("    ✓ Figure3_Severity_Scores.png\n")
cat("    ✓ Figure4_Forest_Plot.png/.pdf\n\n")
cat("  Tablas (formato Word):\n")
cat("    ✓ Table2_Bivariate_Analysis.docx\n")
cat("    ✓ Table3_Multivariable_Model.docx\n\n")

cat("RESULTADOS PRINCIPALES:\n")
cat(sprintf("  • Prevalencia global de LRA: %.1f%% (IC95%% %.1f%%-%.1f%%)\n",
            prev_global$est * 100, prev_global$lower * 100, prev_global$upper * 100))
cat(sprintf("  • Casos analizados: %d\n", n_total))
cat("  • Factores independientes en modelo multivariado: Ver Tabla 3\n\n")

cat("📊 TODOS LOS GRÁFICOS Y TABLAS ESTÁN LISTOS PARA ENVIAR A:\n")
cat("  The Lancet, JAMA, NEJM, BMJ, Critical Care Medicine, Intensive Care Medicine\n\n")

cat("✅ ANÁLISIS COMPLETADO EXITOSAMENTE\n")
cat("========================================\n")

# ============================================
# CÁLCULO DE TAMAÑO DE MUESTRA - COHORTE
# Comparación de dos proporciones
# ============================================

# Cargar librería
library(pwr)  # Para cálculos de poder

# Método 1: Usando pwr.2p.test
# --------------------------------------------

# Parámetros
p0 <- 0.30  # Incidencia en NO expuestos (sedación superficial)
p1 <- 0.50  # Incidencia en EXPUESTOS (sedación profunda)
alpha <- 0.05  # Error tipo I
power <- 0.80  # Poder estadístico (1-β)

# Calcular tamaño de efecto (h de Cohen)
h <- ES.h(p1, p0)
cat(sprintf("Tamaño de efecto h de Cohen: %.3f\n", h))

# Calcular tamaño de muestra
resultado <- pwr.2p.test(
  h = h,
  sig.level = alpha,
  power = power,
  alternative = "two.sided"
)

print(resultado)

n_por_grupo <- ceiling(resultado$n)
cat(sprintf("\n✓ Tamaño de muestra (sin ajustar): %d por grupo\n", n_por_grupo))

# Ajustar por pérdidas
tasa_perdida <- 0.10  # 10% de pérdidas esperadas
n_ajustado <- ceiling(n_por_grupo / (1 - tasa_perdida))
n_total <- n_ajustado * 2

cat(sprintf("\n✓ Tamaño ajustado por %.0f%% pérdidas: %d por grupo\n", 
            tasa_perdida * 100, n_ajustado))
cat(sprintf("✓ Muestra total: %d pacientes\n", n_total))


# Método 2: Usando fórmula manual (más control)
# --------------------------------------------

calcular_muestra_cohorte <- function(p0, p1, alpha = 0.05, power = 0.80, 
                                       bilateral = TRUE, perdidas = 0) {
  # Valores Z
  z_alpha <- if (bilateral) qnorm(1 - alpha/2) else qnorm(1 - alpha)
  z_beta <- qnorm(power)
  
  # Fórmula de Fleiss
  numerador <- (z_alpha + z_beta)^2 * (p1*(1-p1) + p0*(1-p0))
  denominador <- (p1 - p0)^2
  
  n <- numerador / denominador
  n <- ceiling(n)
  
  # Ajustar por pérdidas
  n_ajustado <- if (perdidas > 0) ceiling(n / (1 - perdidas)) else n
  
  # Resultados
  cat("\n========================================\n")
  cat("CÁLCULO DE TAMAÑO DE MUESTRA - COHORTE\n")
  cat("========================================\n\n")
  
  cat("PARÁMETROS:\n")
  cat(sprintf("  p0 (No expuestos): %.1f%%\n", p0*100))
  cat(sprintf("  p1 (Expuestos): %.1f%%\n", p1*100))
  cat(sprintf("  Diferencia absoluta: %.1f%%\n", (p1-p0)*100))
  cat(sprintf("  Riesgo Relativo: %.2f\n", p1/p0))
  cat(sprintf("  Alpha: %.2f (%s)\n", alpha, 
              if(bilateral) "bilateral" else "unilateral"))
  cat(sprintf("  Poder: %.0f%%\n", power*100))
  cat(sprintf("  Z alpha: %.3f\n", z_alpha))
  cat(sprintf("  Z beta: %.3f\n\n", z_beta))
  
  cat("RESULTADOS:\n")
  cat(sprintf("  n por grupo (sin ajustar): %d\n", n))
  cat(sprintf("  Pérdidas esperadas: %.0f%%\n", perdidas*100))
  cat(sprintf("  n por grupo (ajustado): %d\n", n_ajustado))
  cat(sprintf("  MUESTRA TOTAL: %d\n\n", n_ajustado * 2))
  
  # Devolver lista con resultados
  return(list(
    n_sin_ajustar = n,
    n_ajustado = n_ajustado,
    n_total = n_ajustado * 2,
    parametros = list(p0 = p0, p1 = p1, alpha = alpha, power = power,
                       diferencia_absoluta = p1 - p0,
                       riesgo_relativo = p1/p0)
  ))
}

# Ejecutar cálculo
resultado <- calcular_muestra_cohorte(
  p0 = 0.30,
  p1 = 0.50,
  alpha = 0.05,
  power = 0.80,
  bilateral = TRUE,
  perdidas = 0.10
)


# ============================================
# ANÁLISIS DE SENSIBILIDAD
# ¿Qué pasa si cambio los parámetros?
# ============================================

# Crear tabla de diferentes escenarios
escenarios <- expand.grid(
  p0 = c(0.25, 0.30, 0.35),
  diferencia = c(0.15, 0.20, 0.25),
  power = c(0.80, 0.90)
)

escenarios$p1 <- escenarios$p0 + escenarios$diferencia

# Calcular n para cada escenario
escenarios$n_por_grupo <- sapply(1:nrow(escenarios), function(i) {
  res <- calcular_muestra_cohorte(
    p0 = escenarios$p0[i],
    p1 = escenarios$p1[i],
    power = escenarios$power[i],
    perdidas = 0.10
  )
  res$n_ajustado
})

escenarios$n_total <- escenarios$n_por_grupo * 2

print(escenarios[order(escenarios$n_total), ])

# ============================================
# CÁLCULO DE TAMAÑO DE MUESTRA
# CASOS Y CONTROLES
# ============================================

library(epiR)  # Excelente para cálculos epidemiológicos

# ============================================
# MÉTODO 1: Casos-Controles NO PAREADOS
# ============================================

calcular_muestra_cc_no_pareado <- function(p0, OR, m = 1, 
                                                 alpha = 0.05, power = 0.80) {
  
  # Calcular p1 desde OR y p0
  # OR = [p1/(1-p1)] / [p0/(1-p0)]
  # Despejando: p1 = (OR × p0) / (1 - p0 + OR × p0)
  p1 <- (OR * p0) / (1 - p0 + OR * p0)
  
  # Prevalencia promedio ponderada
  p_bar <- (p1 + m * p0) / (m + 1)
  
  # Valores Z
  z_alpha <- qnorm(1 - alpha/2)
  z_beta <- qnorm(power)
  
  # Fórmula
  numerador <- (z_alpha + z_beta)^2 * (m + 1) * p_bar * (1 - p_bar)
  denominador <- m * (p1 - p0)^2
  
  n_casos <- ceiling(numerador / denominador)
  n_controles <- n_casos * m
  n_total <- n_casos + n_controles
  
  # Imprimir resultados
  cat("\n========================================\n")
  cat("CASOS Y CONTROLES - NO PAREADO\n")
  cat("========================================\n\n")
  
  cat("PARÁMETROS:\n")
  cat(sprintf("  p0 (exposición en controles): %.1f%%\n", p0*100))
  cat(sprintf("  p1 (exposición en casos): %.1f%%\n", p1*100))
  cat(sprintf("  Odds Ratio: %.2f\n", OR))
  cat(sprintf("  Ratio controles:casos: %d:1\n", m))
  cat(sprintf("  Alpha: %.2f (bilateral)\n", alpha))
  cat(sprintf("  Poder: %.0f%%\n\n", power*100))
  
  cat("RESULTADOS:\n")
  cat(sprintf("  Casos necesarios: %d\n", n_casos))
  cat(sprintf("  Controles necesarios: %d\n", n_controles))
  cat(sprintf("  MUESTRA TOTAL: %d\n\n", n_total))
  
  return(list(n_casos = n_casos, n_controles = n_controles, 
              n_total = n_total, p1 = p1))
}

# Ejemplo: IBP y NAV
resultado_nopareado <- calcular_muestra_cc_no_pareado(
  p0 = 0.40,    # 40% de controles usan IBP
  OR = 2.25,    # OR esperado
  m = 2,        # 2 controles por caso
  alpha = 0.05,
  power = 0.80
)


# ============================================
# MÉTODO 2: Usando paquete epiR (MÁS PRECISO)
# ============================================

# Para diseño NO pareado
resultado_epir <- epi.ccsize(
  OR = 2.25,
  p0 = 0.40,
  n = NA,          # Lo que queremos calcular
  power = 0.80,
  r = 2,           # 2 controles por caso
  rho = 0,          # 0 = no pareado
  design = 1,      # 1 = no pareado
  sided.test = 2,  # bilateral
  conf.level = 0.95
)

print(resultado_epir)


# ============================================
# MÉTODO 3: Casos-Controles PAREADOS
# ============================================

# Para pareado, usar epiR con rho y design diferentes
resultado_pareado <- epi.ccsize(
  OR = 2.25,
  p0 = 0.40,
  n = NA,
  power = 0.80,
  r = 2,
  rho = 0.2,       # Correlación intra-par (típicamente 0.1-0.3)
  design = 2,      # 2 = pareado
  sided.test = 2,
  conf.level = 0.95
)

cat("\n========================================\n")
cat("CASOS Y CONTROLES - PAREADO 1:2\n")
cat("========================================\n")
print(resultado_pareado)


# ============================================
# ANÁLISIS DE SENSIBILIDAD
# ============================================

# Evaluar diferentes escenarios de OR
escenarios_or <- data.frame(
  OR = c(1.5, 2.0, 2.5, 3.0),
  p0 = 0.40,
  m = 2,
  power = 0.80
)

escenarios_or$n_casos <- sapply(1:nrow(escenarios_or), function(i) {
  res <- calcular_muestra_cc_no_pareado(
    p0 = escenarios_or$p0[i],
    OR = escenarios_or$OR[i],
    m = escenarios_or$m[i],
    power = escenarios_or$power[i]
  )
  res$n_casos
})

escenarios_or$n_total <- escenarios_or$n_casos * (1 + escenarios_or$m)

cat("\n========================================\n")
cat("ANÁLISIS DE SENSIBILIDAD - Diferentes OR\n")
cat("========================================\n")
print(escenarios_or)

# Gráfico
library(ggplot2)
ggplot(escenarios_or, aes(x = OR, y = n_total)) +
  geom_line(color = "steelblue", size = 1.5) +
  geom_point(color = "red", size = 3) +
  geom_text(aes(label = n_total), vjust = -1) +
  theme_minimal() +
  labs(title = "Tamaño de muestra según OR esperado",
       subtitle = "Diseño casos-controles 1:2, p0=40%, poder=80%",
       x = "Odds Ratio",
       y = "Muestra total (casos + controles)")