paquetes_necesarios <- c("survival", "survminer", "dplyr", "kableExtra", "ggplot2")
nuevos <- paquetes_necesarios[!paquetes_necesarios %in% installed.packages()[, "Package"]]
if (length(nuevos) > 0) install.packages(nuevos)

# Cargo mis paquetes
library(survival)
library(survminer)
library(dplyr)
library(kableExtra)
library(ggplot2)

cat("Todos los paquetes cargados correctamente\n")
## Todos los paquetes cargados correctamente
library(haven)

# Leer directamente el .dta — etiquetas de Stata
datos <- read_dta("EMR.dta") %>%
  mutate(across(where(is.labelled), as_factor))

# Verificación
cat("Filas:", nrow(datos), "| Columnas:", ncol(datos), "\n\n")
## Filas: 126 | Columnas: 11
cat("=== Tratamiento ===\n");              print(table(datos$Tratamiento, useNA="always"))
## === Tratamiento ===
## 
## Convencional          EMR         <NA> 
##           64           62            0
cat("\n=== Clasificación diagnóstica ===\n"); print(table(datos$Clasificaciondiagnostica, useNA="always"))
## 
## === Clasificación diagnóstica ===
## 
##    Medico Quiurgico      <NA> 
##        93        33         0
cat("\n=== Extubación ===\n");             print(table(datos$extubacion, useNA="always"))
## 
## === Extubación ===
## 
##    0    1 <NA> 
##   23  103    0

1 Análisis de Supervivencia en Entrenamiento Muscular Respiratorio

1.1 Contexto del estudio

El presente análisis reproduce y amplía los resultados reportados por Sandoval Moreno et al. (2019) en el artículo “Eficacia del entrenamiento muscular respiratorio en el destete de la ventilación mecánica en pacientes con ventilación mecánica por 48 o más horas: un ensayo clínico controlado”, publicado en la revista Medicina Intensiva (vol. 43, núm. 2, pp. 79–89).

El estudio fue conducido en la Unidad de Cuidado Intensivo de la Fundación Valle del Lili, institución de IV nivel en Cali, Colombia. Se trata de un ensayo clínico controlado aleatorizado de grupos paralelos, doble ciego, con una muestra de 126 pacientes en ventilación mecánica por 48 horas o más, aleatorizados mediante bloques de tamaño 4 y 6, estratificados por condición de sepsis.

El grupo experimental recibió un programa de Entrenamiento Muscular Respiratorio (EMR) con dispositivo Threshold IMT ajustado al 50% de la presión inspiratoria máxima (Pimáx), adicional al cuidado estándar.

El grupo convencional recibió únicamente el cuidado estándar de fisioterapia respiratoria.

El desenlace primario fue el tiempo de destete de la ventilación mecánica, definido operativamente como el período comprendido entre el inicio de la ventilación en presión soporte menor a 10 cmH₂O hasta la extubación.

1.2 Clasificación de variables de la base de datos

Antes de realizar los respectivos análisis estadísticos, es fundamental conocer la naturaleza de cada variable, pues la adecuada clasificación determina qué tipo de prueba o modelo voy a usar para cada pregunta de investigación.

# ── Construir el data.frame de clasificación ───────────────────────────
clasificacion <- data.frame(
  `Variable en R`  = c(
    "sexo", "Clasificacion diagnostica", "Tratamiento",
    "extubacion", "diagnostico de ingreso",
    "talla", "peso", "Horas de destete",
    "edad", "apacheii"
  ),
  `Tipo estadístico` = c(
    "Cualitativa dicotómica", "Cualitativa dicotómica",
    "Cualitativa dicotómica", "Cualitativa dicotómica",
    "Cualitativa nominal politómica",
    "Cuantitativa continua", "Cuantitativa continua",
    "Cuantitativa continua",
    "Cuantitativa discreta", "Cuantitativa discreta"
  ),
  `Categoría` = c(
    "Masculino / Femenino",
    "Médico / Quirúrgico",
    "Convencional / EMR",
    "No (0) / Sí (1)",
    "Texto libre (diagnóstico clínico)",
    "Centímetros (cm)",
    "Kilogramos (kg)",
    "Horas",
    "Años cumplidos",
    "Puntaje 0–71"
  ),
  `Rol en el análisis` = c(
    "Variable descriptiva / covariable",
    "Covariable en modelo de Cox",
    " 🔴 Exposición principal",
    "🟢 Variable de evento (supervivencia)",
    "Referencia clínica — no entra en modelos",
    "Variable descriptiva",
    "Variable descriptiva",
    "⏱ Variable de tiempo (supervivencia)",
    "Covariable en modelo de Cox",
    "Covariable en modelo de Cox"
  ),
  check.names = FALSE
)

# ── Índices por grupo  ─────────────────────────────
i_dicotomica  <- 1:4    # cualitativas dicotómicas
i_politomica  <- 5      # cualitativa politómica
i_continua    <- 6:8    # cuantitativas continuas
i_discreta    <- 9:10   # cuantitativas discretas

# ── Tabla con kableExtra ──────────────────────────────────────────────
clasificacion %>%
  kbl(
    caption = "Tabla 0. Clasificación de las variables de la base de datos EMR — Sandoval et al. (2019)",
    align   = c("l", "l", "l", "l"),
    escape  = FALSE
  ) %>%
  kable_styling(
    bootstrap_options = c("hover", "condensed", "bordered"),
    full_width        = TRUE,
    font_size         = 13
  ) %>%
  # ── Encabezado principal ─────────────────────────────────────────────
  row_spec(0,
           bold       = TRUE,
           color      = "white",
           background = "#4A235A",    # morado oscuro
           font_size  = 13) %>%
  # ── Grupo 1: Cualitativas dicotómicas — lila suave ───────────────────
  row_spec(i_dicotomica,
           background = "#EDE7F6",
           color      = "#4A235A") %>%
  # ── Grupo 2: Cualitativa politómica — rosa suave ─────────────────────
  row_spec(i_politomica,
           background = "#FCE4EC",
           color      = "#880E4F") %>%
  # ── Grupo 3: Cuantitativas continuas — azul muy claro ────────────────
  row_spec(i_continua,
           background = "#E3F2FD",
           color      = "#0D47A1") %>%
  # ── Grupo 4: Cuantitativas discretas — verde menta ───────────────────
  row_spec(i_discreta,
           background = "#E8F5E9",
           color      = "#1B5E20") %>%
  # ── Primera columna en negrita ────────────────────────────────────────
  column_spec(1, bold = TRUE, width = "20%") %>%
  column_spec(2, width = "22%") %>%
  column_spec(3, width = "28%") %>%
  column_spec(4, width = "30%") %>%
  # ── Agrupadores de filas ─────────────────────────────────────────────
  pack_rows("⬜ Cualitativas dicotómicas",    1, 4,
            label_row_css = "background:#CE93D8; color:#4A235A;
                             font-weight:bold; font-size:12px;") %>%
  pack_rows(" ⚪ Cualitativa nominal politómica", 5, 5,
            label_row_css = "background:#F48FB1; color:#880E4F;
                             font-weight:bold; font-size:12px;") %>%
  pack_rows("🟧 Cuantitativas continuas",     6, 8,
            label_row_css = "background:#90CAF9; color:#0D47A1;
                             font-weight:bold; font-size:12px;") %>%
  pack_rows("🟡 Cuantitativas discretas",     9, 10,
            label_row_css = "background:#A5D6A7; color:#1B5E20;
                             font-weight:bold; font-size:12px;") %>%
  footnote(
    general = paste0(
      "Fuente: Base de datos del ensayo clínico controlado Sandoval Moreno et al. ",
      "Med Intensiva. 2019;43(2):79–89. ",
      " 🔴= variable de exposición principal | ",
      " 🟢 = variable de evento | ⏱ = variable de tiempo."
    ),
    general_title = "Nota:",
    footnote_as_chunk = TRUE
  )
Tabla 0. Clasificación de las variables de la base de datos EMR — Sandoval et al. (2019)
Variable en R Tipo estadístico Categoría Rol en el análisis
⬜ Cualitativas dicotómicas
sexo Cualitativa dicotómica Masculino / Femenino Variable descriptiva / covariable
Clasificacion diagnostica Cualitativa dicotómica Médico / Quirúrgico Covariable en modelo de Cox
Tratamiento Cualitativa dicotómica Convencional / EMR 🔴 Exposición principal
extubacion Cualitativa dicotómica No (0) / Sí (1) 🟢 Variable de evento (supervivencia)
⚪ Cualitativa nominal politómica
diagnostico de ingreso Cualitativa nominal politómica Texto libre (diagnóstico clínico) Referencia clínica — no entra en modelos
🟧 Cuantitativas continuas
talla Cuantitativa continua Centímetros (cm) Variable descriptiva
peso Cuantitativa continua Kilogramos (kg) Variable descriptiva
Horas de destete Cuantitativa continua Horas ⏱ Variable de tiempo (supervivencia)
🟡 Cuantitativas discretas
edad Cuantitativa discreta Años cumplidos Covariable en modelo de Cox
apacheii Cuantitativa discreta Puntaje 0–71 Covariable en modelo de Cox
Nota: Fuente: Base de datos del ensayo clínico controlado Sandoval Moreno et al. Med Intensiva. 2019;43(2):79–89. 🔴= variable de exposición principal | 🟢 = variable de evento | ⏱ = variable de tiempo.

1.3 Exploración y contraste con el artículo original

A continuación paso a verificar si las características de la muestra en la base de datos coinciden con lo reportado por Sandoval et al. (2019) en la Tabla 1 del artículo, pues este paso es esencial, ya que me permite validar la integridad de los datos antes de proceder con los análisis de supervivencia.

Glosario Artículo reporta:

EMR(Entrenamiento Muscular Respiratorio) Convencional(cuidado estándar de fisioterapia respiratoria en UCI)

datos <- datos %>%
  mutate(
    sexo_f   = factor(sexo,
                      levels = c("Masculino", "Femenino")),

    trat_f   = factor(Tratamiento,
                      levels = c("Convencional", "EMR")),

    clasif_f = factor(
                 case_when(
                   Clasificaciondiagnostica == "Medico" ~ "Médico",
                   TRUE ~ "Quirúrgico"),
                 levels = c("Médico", "Quirúrgico")),

    evento   = as.integer(extubacion),
    tiempo   = Horasdedestete
  )

# ── Hago verificación inmediata ─────────────────────────────────────────────
cat("variables creadas correctamente\n\n")
## variables creadas correctamente
cat("trat_f:  "); print(table(datos$trat_f))
## trat_f:
## 
## Convencional          EMR 
##           64           62
cat("clasif_f:"); print(table(datos$clasif_f))
## clasif_f:
## 
##     Médico Quirúrgico 
##         93         33
cat("evento:  "); print(table(datos$evento))
## evento:
## 
##   0   1 
##  23 103
# ── Ver valores crudos exactos de cada variable ───────────
cat("=== Tratamiento (valores únicos) ===\n")
## === Tratamiento (valores únicos) ===
print(table(datos$Tratamiento, useNA = "always"))
## 
## Convencional          EMR         <NA> 
##           64           62            0
cat("\n=== Clasificación diagnóstica (valores únicos) ===\n")
## 
## === Clasificación diagnóstica (valores únicos) ===
print(table(datos$Clasificaciondiagnostica, useNA = "always"))
## 
##    Medico Quiurgico      <NA> 
##        93        33         0
cat("\n=== Extubacion (valores únicos) ===\n")
## 
## === Extubacion (valores únicos) ===
print(table(datos$extubacion, useNA = "always"))
## 
##    0    1 <NA> 
##   23  103    0
cat("\n=== Filas con Tratamiento vacío o NA ===\n")
## 
## === Filas con Tratamiento vacío o NA ===
print(sum(is.na(datos$Tratamiento) | datos$Tratamiento == "" | 
          datos$Tratamiento == "nan"))
## [1] 0

1.4 Verificación de integridad de la base de datos

Antes de responder las preguntas del taller, he realizado un proceso de aseguramiento y calidad de los datos, haciendo una verificación exhaustiva de la base de datos contrastando cada variable con los valores reportados por Sandoval Moreno et al. (2019) en la Tabla 1 del artículo original.

1.4.1 Resultado de la verificación

La base de datos ha quedado completamente íntegra, sin valores perdidos en las variables de análisis principal. Los conteos observados coinciden exactamente con los reportados en el artículo:

verificacion <- data.frame(
  `Variable verificada` = c(
    "Grupo Convencional",
    "Grupo EMR",
    "Total muestra",
    "Diagnóstico médico",
    "Diagnóstico quirúrgico",
    "Extubados exitosamente",
    "Censurados (no extubados)",
    "Valores perdidos en variables clave"
  ),
  `Observado en base` = c(
    "64 pacientes", "62 pacientes", "126 pacientes",
    "93 pacientes (73.8%)", "33 pacientes (26.2%)",
    "103 pacientes (81.7%)", "23 pacientes (18.3%)", "0"
  ),
  `Reportado en artículo` = c(
    "64 pacientes", "62 pacientes", "126 pacientes",
    "93 pacientes (73.8%)", "33 pacientes (26.2%)",
    "~103 pacientes", "~23 pacientes", "—"
  ),
  `Estado` = rep("&#10003;", 8),   # marca negra
  check.names = FALSE
)

verificacion %>%
  kbl(
    caption = "Tabla de verificación de la calidad e Integridad de los datosde la base EMR vs. Sandoval et al. (2019)",
    align   = c("l","c","c","c"),
    escape  = FALSE
  ) %>%
  kable_styling(
    bootstrap_options = c("hover","bordered","condensed"),
    full_width        = TRUE,
    font_size         = 13,
    position          = "center"
  ) %>%
  row_spec(0,
           bold       = TRUE,
           color      = "white",
           background = "#4A235A") %>%
  row_spec(c(1,2,3),
           background = "#EDE7F6",
           color      = "#4A235A") %>%
  row_spec(c(4,5),
           background = "#E3F2FD",
           color      = "#0D47A1") %>%
  row_spec(c(6,7),
           background = "#E8F5E9",
           color      = "#1B5E20") %>%
  row_spec(8,
           background = "#FFF9C4",
           color      = "#F57F17") %>%
  column_spec(1, bold = TRUE, width = "35%") %>%
  column_spec(4,
            bold  = TRUE,
            color = "#1B5E20",
            width = "8%") %>%
  pack_rows("Grupos de tratamiento",  1, 3,
            label_row_css = "background:#CE93D8; color:#4A235A;
                             font-weight:bold; font-size:12px;") %>%
  pack_rows("Clasificación diagnóstica", 4, 5,
            label_row_css = "background:#90CAF9; color:#0D47A1;
                             font-weight:bold; font-size:12px;") %>%
  pack_rows("Variable de evento", 6, 7,
            label_row_css = "background:#A5D6A7; color:#1B5E20;
                             font-weight:bold; font-size:12px;") %>%
  pack_rows("Calidad del dato", 8, 8,
            label_row_css = "background:#FFF176; color:#F57F17;
                             font-weight:bold; font-size:12px;") %>%
  footnote(
    general       = "Verificación realizada contra Tabla 1 y sección de resultados de Sandoval Moreno et al. Med Intensiva. 2019;43(2):79–89. ✓ = coincidencia exacta confirmada.",
    general_title = "Fuente:",
    footnote_as_chunk = TRUE
  )
Tabla de verificación de la calidad e Integridad de los datosde la base EMR vs. Sandoval et al. (2019)
Variable verificada Observado en base Reportado en artículo Estado
Grupos de tratamiento
Grupo Convencional 64 pacientes 64 pacientes
Grupo EMR 62 pacientes 62 pacientes
Total muestra 126 pacientes 126 pacientes
Clasificación diagnóstica
Diagnóstico médico 93 pacientes (73.8%) 93 pacientes (73.8%)
Diagnóstico quirúrgico 33 pacientes (26.2%) 33 pacientes (26.2%)
Variable de evento
Extubados exitosamente 103 pacientes (81.7%) ~103 pacientes
Censurados (no extubados) 23 pacientes (18.3%) ~23 pacientes
Calidad del dato
Valores perdidos en variables clave 0
Fuente: Verificación realizada contra Tabla 1 y sección de resultados de Sandoval Moreno et al. Med Intensiva. 2019;43(2):79–89. ✓ = coincidencia exacta confirmada.

2 Respuesta a las preguntas del taller

2.1 Pregunta 1: ¿Qué proporción total de pacientes fueron extubados de manera exitosa?

Para responder esta pregunta calcularé la proporción global (p) de extubación exitosa sobre el total de 126 pacientes aleatorizados, pues concomitantemente en epidemiología clínica, una proporción aislada cobra sentido real cuando la acompañamos de su intervalo de confianza al 95%, que nos dice dentro de qué rango de valores podría ubicarse la proporción verdadera en la población si replicáramos este estudio muchas veces.

Por tanto, usaré el método de Wilson, que es más preciso que el método clásico de Wald cuando las proporciones se acercan a 0 o a 1. Considero que esta elección es la más adecuada debido a que la proporción observada es p = 103/126 = 81.7%, ** es decir, 8 de cada 10 pacientes lograron ser extubados, un valor que supera el umbral crítico del 80% por lo cual, el método de Wald a partir de ahí, pierde precisión, pues asume una simetría en el intervalo que la distribución real ya no tiene.

Además a ese nivel existe más incertidumbre hacia abajo que hacia arriba, y Wald simplemente la ignora, por tanto, Wilson me garantiza límites biológicamente posibles y un intervalo que respeta la asimetría real de los datos.

options(digits = 15)

# ── Datos base ──────────────────────────────────────────────────────────
n <- 126L   # total pacientes 
x <- 103L   # extubados exitosamente 
p <- x / n  # proporción observada — precisión completa

cat("PASO 1 — Proporción observada\n")
## PASO 1 — Proporción observada
cat("─────────────────────────────\n")
## ─────────────────────────────
cat("p = x / n =", x, "/", n, "=", p, "\n")
## p = x / n = 103 / 126 = 0.817460317460317
cat("p en porcentaje =", p * 100, "%\n\n")
## p en porcentaje = 81.7460317460317 %
# ── Error estándar Wald ─────────────────────────────────────────────────
Z  <- qnorm(0.975)   # valor exacto de Z para IC 95% bilateral
ee <- sqrt((p * (1 - p)) / n)

cat("PASO 2 — Error estándar de Wald\n")
## PASO 2 — Error estándar de Wald
cat("─────────────────────────────────────────────────────────\n")
## ─────────────────────────────────────────────────────────
cat("Z (qnorm 0.975) =", Z, "\n")
## Z (qnorm 0.975) = 1.95996398454005
cat("EE = sqrt(p × (1−p) / n)\n")
## EE = sqrt(p × (1−p) / n)
cat("EE = sqrt(", p, "×", 1-p, "/", n, ")\n")
## EE = sqrt( 0.817460317460317 × 0.182539682539683 / 126 )
cat("EE = sqrt(", (p*(1-p))/n, ")\n")
## EE = sqrt( 0.0011842773558571 )
cat("EE =", ee, "\n\n")
## EE = 0.0344133310776086
# ── IC Wald ─────────────────────────────────────────────────────────────
wald_inf <- p - Z * ee
wald_sup <- p + Z * ee

cat("PASO 3 — IC 95% Wald (solo para contrastar)\n")
## PASO 3 — IC 95% Wald (solo para contrastar)
cat("────────────────────────────────────────────\n")
## ────────────────────────────────────────────
cat("LI = p − Z × EE =", p, "−", Z * ee, "=", wald_inf, "\n")
## LI = p − Z × EE = 0.817460317460317 − 0.0674488895001657 = 0.750011427960152
cat("LS = p + Z × EE =", p, "+", Z * ee, "=", wald_sup, "\n")
## LS = p + Z × EE = 0.817460317460317 + 0.0674488895001657 = 0.884909206960483
cat("LI en % =", wald_inf * 100, "\n")
## LI en % = 75.0011427960152
cat("LS en % =", wald_sup * 100, "\n\n")
## LS en % = 88.4909206960483
# ── IC Wilson ───────────────────────────────────────────────────────────
Z2         <- Z^2
p_wilson   <- (x + Z2/2) / (n + Z2)
margen     <- (Z / (n + Z2)) * sqrt((x*(n-x)/n) + Z2/4)
wilson_inf <- p_wilson - margen
wilson_sup <- p_wilson + margen

cat("PASO 4 — IC 95% Wilson (método seleccionado)\n")
## PASO 4 — IC 95% Wilson (método seleccionado)
cat("─────────────────────────────────────────────\n")
## ─────────────────────────────────────────────
cat("Z²  =", Z2, "\n")
## Z²  = 3.84145882069412
cat("Centro Wilson = (x + Z²/2) / (n + Z²)\n")
## Centro Wilson = (x + Z²/2) / (n + Z²)
cat("             = (", x, "+", Z2/2, ") / (", n, "+", Z2, ")\n")
##              = ( 103 + 1.92072941034706 ) / ( 126 + 3.84145882069412 )
cat("             =", x + Z2/2, "/", n + Z2, "\n")
##              = 104.920729410347 / 129.841458820694
cat("             =", p_wilson, "\n\n")
##              = 0.808068011275493
cat("Margen       =", margen, "\n")
## Margen       = 0.0671041878939019
cat("LI Wilson    =", p_wilson, "−", margen, "=", wilson_inf, "\n")
## LI Wilson    = 0.808068011275493 − 0.0671041878939019 = 0.740963823381591
cat("LS Wilson    =", p_wilson, "+", margen, "=", wilson_sup, "\n")
## LS Wilson    = 0.808068011275493 + 0.0671041878939019 = 0.875172199169394
cat("LI en %      =", wilson_inf * 100, "\n")
## LI en %      = 74.0963823381591
cat("LS en %      =", wilson_sup * 100, "\n\n")
## LS en %      = 87.5172199169395
# ── Comparación final ───────────────────────────────────────────────────
cat("PASO 5 — hago la comparación Wald vs Wilson\n")
## PASO 5 — hago la comparación Wald vs Wilson
cat("─────────────────────────────────────────────────────────\n")
## ─────────────────────────────────────────────────────────
cat("Wald   LI:", wald_inf * 100,   "% — LS:", wald_sup * 100,   "%",
    "| amplitud:", (wald_sup - wald_inf) * 100, "pp\n")
## Wald   LI: 75.0011427960152 % — LS: 88.4909206960483 % | amplitud: 13.4897779000331 pp
cat("Wilson LI:", wilson_inf * 100, "% — LS:", wilson_sup * 100, "%",
    "| amplitud:", (wilson_sup - wilson_inf) * 100, "pp\n")
## Wilson LI: 74.0963823381591 % — LS: 87.5172199169395 % | amplitud: 13.4208375787804 pp
resultados <- data.frame(
  `Elemento` = c(
    "Proporción observada (p)",
    "Centro del intervalo",
    "Límite Inferior (LI) 95%",
    "Límite Superior (LS) 95%",
    "Amplitud del intervalo",
    "Distancia p → Límite Inferior",
    "Distancia p → Límite Superior",
    "¿Intervalo simétrico?"
  ),
  `Wald` = c(
    "81.7460317460317%",
    "81.7460% (= p, sin ajuste)",
    "75.0011427960152%",
    "88.4909206960483%",
    "13.4898 pp",
    "6.7449 pp",
    "6.7449 pp",
    "Sí — distancias iguales ⚠️"
  ),
  `Wilson` = c(
    "81.7460317460317%",
    "80.8068% (ajustado hacia 50%)",
    "74.0963823381591%",
    "87.5172199169395%",
    "13.4208 pp",
    "7.6496 pp",
    "5.7712 pp",
    "No — asimetría correcta ✓"
  ),
  `¿Qué me dice este valor?` = c(
    "8 de cada 10 pacientes fueron extubados exitosamente",
    "Wilson desplaza el centro para corregir la asimetría real",
    "Mínimo plausible de extubación en la población real",
    "Máximo plausible de extubación en la población real",
    "Ancho total del intervalo en puntos porcentuales (pp). Mientras más estrecho, más precisa la estimación",
    "Cuánto baja el IC desde p. Wald subestima esta distancia cuando p > 80%",
    "Cuánto sube el IC desde p. A proporciones altas hay menos espacio hacia arriba",
    "Wald asume igual incertidumbre en ambas direcciones — incorrecto cuando p > 80%"
  ),
  check.names = FALSE
)

resultados %>%
  kbl(
    caption = "Tabla 1b. Comparación detallada Wald vs Wilson — Pregunta 1 (p = 81.75%, n = 126)",
    align   = c("l","c","c","l"),
    escape  = FALSE
  ) %>%
  kable_styling(
    bootstrap_options = c("hover","bordered","condensed"),
    full_width        = TRUE,
    font_size         = 12,
    position          = "center"
  ) %>%
  # Encabezado
  row_spec(0,
           bold       = TRUE,
           color      = "white",
           background = "#4A235A") %>%
  # Fila 1 — proporción (igual en ambos)
  row_spec(1,
           background = "#F3E5F5",
           color      = "#4A235A",
           bold       = TRUE) %>%
  # Fila 2 — centro
  row_spec(2,
           background = "#EDE7F6",
           color      = "#4A235A") %>%
  # Filas 3-4 — límites
  row_spec(3,
           background = "#E3F2FD",
           color      = "#0D47A1") %>%
  row_spec(4,
           background = "#BBDEFB",
           color      = "#0D47A1") %>%
  # Fila 5 — amplitud
  row_spec(5,
           background = "#E8F5E9",
           color      = "#1B5E20",
           bold       = TRUE) %>%
  # Filas 6-7 — distancias
  row_spec(6,
           background = "#FFF9C4",
           color      = "#F57F17") %>%
  row_spec(7,
           background = "#FFF9C4",
           color      = "#F57F17") %>%
  # Fila 8 — simetría
  row_spec(8,
           background = "#FCE4EC",
           color      = "#880E4F",
           bold       = TRUE) %>%
  column_spec(1, bold = TRUE, width = "22%") %>%
  column_spec(2, width = "20%") %>%
  column_spec(3, width = "20%", bold = TRUE, color = "#1B5E20") %>%
  column_spec(4, width = "38%", italic = TRUE) %>%
  add_header_above(
    c(" " = 1,
      "Resultado por método" = 2,
      "Interpretación del valor" = 1),
    bold       = TRUE,
    color      = "white",
    background = "#6A1B9A"
  ) %>%
  footnote(
    general = "pp = puntos porcentuales. Wilson es el método que he seleccionado porque p = 81.75% supera el umbral del 80%, generando asimetría real que Wald no captura. La amplitud de Wilson (13.42 pp) es 0.07 pp más estrecha que Wald, lo que me indica una estimación más precisa.",
    general_title = "Conclusión metodológica:",
    footnote_as_chunk = TRUE
  )
Tabla 1b. Comparación detallada Wald vs Wilson — Pregunta 1 (p = 81.75%, n = 126)
Resultado por método
Interpretación del valor
Elemento Wald Wilson ¿Qué me dice este valor?
Proporción observada (p) 81.7460317460317% 81.7460317460317% 8 de cada 10 pacientes fueron extubados exitosamente
Centro del intervalo 81.7460% (= p, sin ajuste) 80.8068% (ajustado hacia 50%) Wilson desplaza el centro para corregir la asimetría real
Límite Inferior (LI) 95% 75.0011427960152% 74.0963823381591% Mínimo plausible de extubación en la población real
Límite Superior (LS) 95% 88.4909206960483% 87.5172199169395% Máximo plausible de extubación en la población real
Amplitud del intervalo 13.4898 pp 13.4208 pp Ancho total del intervalo en puntos porcentuales (pp). Mientras más estrecho, más precisa la estimación
Distancia p → Límite Inferior 6.7449 pp 7.6496 pp Cuánto baja el IC desde p. Wald subestima esta distancia cuando p > 80%
Distancia p → Límite Superior 6.7449 pp 5.7712 pp Cuánto sube el IC desde p. A proporciones altas hay menos espacio hacia arriba
¿Intervalo simétrico? Sí — distancias iguales ⚠️ No — asimetría correcta ✓ Wald asume igual incertidumbre en ambas direcciones — incorrecto cuando p > 80%
Conclusión metodológica: pp = puntos porcentuales. Wilson es el método que he seleccionado porque p = 81.75% supera el umbral del 80%, generando asimetría real que Wald no captura. La amplitud de Wilson (13.42 pp) es 0.07 pp más estrecha que Wald, lo que me indica una estimación más precisa.

2.1.1 Conclusión Pregunta 1

Del total de 126 pacientes del ensayo de Sandoval et al. (2019), 103 pacientes fueron extubados exitosamente, lo que representa una proporción de p = 0.817460317460317 (81.75%).

Con el método de Wilson, evidencié que el intervalo de confianza al 95% se ubica entre 74.0963823381591% (Límite Inferior) y 87.5172199169395% (Límite Superior), con una amplitud de 13.4208 puntos porcentuales.

Esto significa que si se replicara este ensayo muchas veces bajo las mismas condiciones,en el 95% de las réplicas la verdadera proporción de extubación exitosa en pacientes con ventilación mecánica ≥ 48 horas se ubicaría dentro de ese rango.

Los 23 pacientes restantes (18.25%) no alcanzaron el evento y constituyen las observaciones censuradas del análisis de supervivencia.

3 Pregunta 2 ¿Existe evidencia de asociación entre extubación exitosa y el tipo de tratamiento recibido?**

Para responder esta pregunta construiré una tabla de contingencia 2×2 que cruza la variable de evento (extubacion: No/Sí) con la variable de exposición principal (Tratamiento: Convencional/EMR).

Esta estructura me permitirá visualizar simultáneamente cuántos pacientes de cada grupo lograron o no la extubación exitosa, por lo cual, utilizaré la prueba de chi-cuadrado de Pearson para evaluar si existe asociación estadísticamente significativa entre ambas variables.

La lógica de esta prueba es comparar las frecuencias que se observan en la tabla con las frecuencias que se esperaría si el tratamiento no tuviera ningún efecto sobre la extubación, es decir, bajo la hipótesis nula de independencia.

Si la diferencia entre lo observado y lo esperado es suficientemente grande, rechazo esa hipótesis. Además complementaré con la razón de momios (OR) calculada por el método exacto de Fisher, que cuantificará cuántas veces más probable es ser extubado exitosamente en un grupo respecto al otro, con su IC 95%.

Las hipótesis que evaluaré son:

  • H₀: La proporción de extubación exitosa es igual en ambos grupos de tratamiento, por tanto no existe asociación.

  • H₁: La proporción de extubación exitosa difiere entre los grupos, por tanto existe asociación.

El nivel de significancia es α = 0.05.

# ── Tabla de contingencia cruda ────────────────────────────────────────
tabla_cruda <- table(
  Tratamiento  = datos$trat_f,
  Extubacion   = factor(datos$evento,
                        levels = c(0, 1),
                        labels = c("No exitosa", "Exitosa"))
)

cat("=== Frecuencias OBSERVADAS ===\n")
## === Frecuencias OBSERVADAS ===
print(addmargins(tabla_cruda))
##               Extubacion
## Tratamiento    No exitosa Exitosa Sum
##   Convencional         12      52  64
##   EMR                  11      51  62
##   Sum                  23     103 126
cat("\n=== Frecuencias ESPERADAS (bajo H₀) ===\n")
## 
## === Frecuencias ESPERADAS (bajo H₀) ===
esperadas <- chisq.test(tabla_cruda)$expected
print(round(esperadas, 6))
##               Extubacion
## Tratamiento    No exitosa  Exitosa
##   Convencional   11.68254 52.31746
##   EMR            11.31746 50.68254
cat("\n=== Verificación: todas las esperadas > 5? ===\n")
## 
## === Verificación: todas las esperadas > 5? ===
cat(ifelse(all(esperadas >= 5),
           "✓ Sí — condición del chi-cuadrado cumplida\n",
           "✗ No — considerar prueba exacta de Fisher\n"))
## ✓ Sí — condición del chi-cuadrado cumplida
obs <- as.data.frame.matrix(addmargins(tabla_cruda))
esp <- chisq.test(tabla_cruda)$expected

# ── Construyo tabla unificada ──────────────────────────────────────────
tabla_2x2 <- data.frame(
  `Grupo` = c("Convencional", "EMR", "Total"),

  `Obs. No exitosa` = c(12, 11, 23),
  `Esp. No exitosa` = c(esp[1,1], esp[2,1], NA),

  `Obs. Exitosa`    = c(52, 51, 103),
  `Esp. Exitosa`    = c(esp[1,2], esp[2,2], NA),

  `Total`           = c(64, 62, 126),

  check.names = FALSE
)

tabla_2x2 %>%
  kbl(
    caption = "Tabla 2. Frecuencias observadas y esperadas — Extubación exitosa según grupo de tratamiento",
    align   = c("l","c","c","c","c","c"),
    digits  = 15
  ) %>%
  kable_styling(
    bootstrap_options = c("hover","bordered","condensed"),
    full_width        = FALSE,
    font_size         = 13,
    position          = "center"
  ) %>%
  # ── Encabezado principal ──────────────────────────────────────────────
  row_spec(0,
           bold       = TRUE,
           color      = "white",
           background = "#4A235A") %>%
  # ── Fila Convencional ─────────────────────────────────────────────────
  row_spec(1, background = "#EDE7F6", color = "#4A235A") %>%
  # ── Fila EMR ──────────────────────────────────────────────────────────
  row_spec(2, background = "#E3F2FD", color = "#0D47A1") %>%
  # ── Fila Total ────────────────────────────────────────────────────────
  row_spec(3,
           background = "#E8F5E9",
           color      = "#1B5E20",
           bold       = TRUE) %>%
  column_spec(1, bold = TRUE, width = "18%") %>%
  # Columnas observadas — más oscuras
  column_spec(2, background = "#FCE4EC", color = "#880E4F", width = "14%") %>%
  column_spec(4, background = "#D5F5E3", color = "#1B5E20", width = "14%") %>%
  # Columnas esperadas — más claras
  column_spec(3, background = "#FFF9C4", color = "#795548", width = "16%") %>%
  column_spec(5, background = "#FFF9C4", color = "#795548", width = "16%") %>%
  column_spec(6, bold = TRUE, width = "10%") %>%
  add_header_above(
    c(" "            = 1,
      "No exitosa"   = 2,
      "Exitosa"      = 2,
      " "            = 1),
    bold       = TRUE,
    color      = "white",
    background = "#6A1B9A"
  ) %>%
  add_header_above(
    c(" "                    = 1,
      "Frecuencias por grupo" = 4,
      " "                    = 1),
    bold       = TRUE,
    color      = "white",
    background = "#4A235A"
  ) %>%
  footnote(
    general = paste0(
      "Obs. = frecuencia observada en los datos. ",
      "Esp. = frecuencia esperada bajo H₀ (independencia entre tratamiento y extubación). ",
      "Todas las frecuencias esperadas superan 5 — condición del chi-cuadrado cumplida ✓"
    ),
    general_title = "Nota:",
    footnote_as_chunk = TRUE
  )
Tabla 2. Frecuencias observadas y esperadas — Extubación exitosa según grupo de tratamiento
Frecuencias por grupo
No exitosa
Exitosa
Grupo Obs. No exitosa Esp. No exitosa Obs. Exitosa Esp. Exitosa Total
Convencional 12 11.6825396825397 52 52.3174603174603 64
EMR 11 11.3174603174603 51 50.6825396825397 62
Total 23 NA 103 NA 126
Nota: Obs. = frecuencia observada en los datos. Esp. = frecuencia esperada bajo H₀ (independencia entre tratamiento y extubación). Todas las frecuencias esperadas superan 5 — condición del chi-cuadrado cumplida ✓

3.0.1 Análisis de la tabla de contingencia

Lo que se observa en los datos reales

La tabla revela un hallazgo inmediatamente llamativo, como es la distribución de extubaciones exitosas entre los dos grupos que es prácticamente idéntica.

Del grupo Convencional, 52 de 64 pacientes (81.25%) lograron ser extubados exitosamente, mientras que del grupo de Entrenamiento Muscular Respiratorio, 51 de 62 pacientes (82.26%) alcanzaron el mismo desenlace. La diferencia entre grupos es de apenas 1.01 puntos porcentuales — una magnitud tan pequeña que difícilmente podría tener relevancia clínica; del mismo modo, los pacientes que no lograron ser extubados se distribuyen de forma casi equitativa, con 12 en el grupo Convencional y 11 en el grupo de Entrenamiento Muscular Respiratorio, diferencia de un solo paciente sobre un total de 126.

Lo que se esperaría si el tratamiento no tuviera ningún efecto: Las frecuencias esperadas bajo la hipótesis nula de independencia nos dicen cuántos pacientes habría en cada celda si el tipo de tratamiento recibido no guardara ninguna relación con el resultado de la extubación,**es decir, si ambos grupos tuvieran exactamente la misma probabilidad de ser extubados.

Bajo esa hipótesis esperaríamos:

  • En el grupo Convencional: 11.6825396825397 pacientes sin extubación exitosa y 52.3174603174603 con extubación exitosa.

  • En el grupo de Entrenamiento Muscular Respiratorio: 11.3174603174603 pacientes sin extubación exitosa y 50.6825396825397 con extubación exitosa.

Mi análisis en cuanto a lo observado vs. esperado: Al comparar ambas matrices, la diferencia entre lo observado y lo esperado es mínima en todas las celdas, pues en ningún caso supera 0.32 pacientes.

Esta cercanía extrema entre frecuencias observadas y esperadas es la señal más clara, incluso antes de aplicar la prueba formal, de que el tratamiento no está produciendo un patrón diferente al que existiría por puro azar.

Cuando lo que se observa se parece mucho a lo esperado bajo independencia, el estadístico chi-cuadrado será pequeño y el valor p será grande , por lo cual,esto equivale a decir que no hay evidencia estadística de asociación.

Verificación del supuesto fundamental

La prueba chi-cuadrado funciona aproximando la distribución real de los datos a una distribución teórica continua, por lo cual, esa aproximación solo es válida cuando hay suficientes observaciones en cada celda, por lo que el umbral convencional es que ninguna frecuencia esperada sea menor a 5.

Si alguna celda tuviera menos de 5 casos esperados, la aproximación se vuelve imprecisa y el valor p resultante no sería confiable, obligándonos a usar la prueba exacta de Fisher, que no depende de esa aproximación sino que calcula la probabilidad de forma directa sobre todas las tablas posibles con los mismos marginales.

En este caso, la frecuencia esperada más pequeña es 11.317, algo muy por encima del umbral, de modo que la aproximación chi-cuadrado es válida y sus resultados se pueden interpretar sin correcciones.

options(digits = 15)

# ── Prueba chi-cuadrado ────────────────────────────────────────────────
chi   <- chisq.test(tabla_cruda, correct = FALSE)

# ── Prueba exacta de Fisher ────────────────────────────────────────────
fisher <- fisher.test(tabla_cruda)

cat("=== CHI-CUADRADO DE PEARSON ===\n")
## === CHI-CUADRADO DE PEARSON ===
cat("Estadístico X²:", chi$statistic, "\n")
## Estadístico X²: 0.0214463704571141
cat("Grados de libertad:", chi$parameter, "\n")
## Grados de libertad: 1
cat("Valor p:", chi$p.value, "\n\n")
## Valor p: 0.883569478909719
cat("=== PRUEBA EXACTA DE FISHER ===\n")
## === PRUEBA EXACTA DE FISHER ===
cat("OR (razón de momios):", fisher$estimate, "\n")
## OR (razón de momios): 1.06932731321789
cat("Límite Inferior IC 95%:", fisher$conf.int[1], "\n")
## Límite Inferior IC 95%: 0.391715289165425
cat("Límite Superior IC 95%:", fisher$conf.int[2], "\n")
## Límite Superior IC 95%: 2.94487912243686
cat("Valor p:", fisher$p.value, "\n")
## Valor p: 1
tabla_pruebas <- data.frame(
  `Prueba` = c(
    "Chi-cuadrado de Pearson",
    "Chi-cuadrado de Pearson",
    "Chi-cuadrado de Pearson",
    "Prueba exacta de Fisher",
    "Prueba exacta de Fisher",
    "Prueba exacta de Fisher",
    "Prueba exacta de Fisher"
  ),
  `Elemento` = c(
    "Estadístico X²",
    "Grados de libertad",
    "Valor p",
    "Razón de momios (OR)",
    "Límite Inferior IC 95%",
    "Límite Superior IC 95%",
    "Valor p"
  ),
  `Resultado` = c(
    chi$statistic,
    chi$parameter,
    chi$p.value,
    fisher$estimate,
    fisher$conf.int[1],
    fisher$conf.int[2],
    fisher$p.value
  ),
  `Interpretación` = c(
    "Valor muy cercano a 0 — diferencia mínima entre observado y esperado",
    "Tabla 2×2 tiene 1 grado de libertad",
    "Muy superior a α = 0.05 — no se rechaza H₀",
    "OR ≈ 1 — prácticamente sin diferencia entre grupos",
    "Límite inferior cruza el 1 — compatible con ausencia de asociación",
    "Límite superior amplio — estimación imprecisa por grupos similares",
    "p = 1 — máxima compatibilidad con H₀ de independencia"
  ),
  check.names = FALSE
)

tabla_pruebas %>%
  kbl(
    caption = "Tabla 3. Resultados de las pruebas de asociación sobre la Extubación exitosa vs. Tratamiento",
    align   = c("l","l","l","l"),
    digits  = 15
  ) %>%
  kable_styling(
    bootstrap_options = c("hover","bordered","condensed"),
    full_width        = TRUE,
    font_size         = 13,
    position          = "center"
  ) %>%
  row_spec(0,
           bold       = TRUE,
           color      = "white",
           background = "#4A235A") %>%
  # Chi-cuadrado — azul claro
  row_spec(1, background = "#E3F2FD", color = "#0D47A1") %>%
  row_spec(2, background = "#E3F2FD", color = "#0D47A1") %>%
  row_spec(3,
           background = "#BBDEFB",
           color      = "#0D47A1",
           bold       = TRUE) %>%
  # Fisher — verde claro
  row_spec(4,
           background = "#E8F5E9",
           color      = "#1B5E20",
           bold       = TRUE) %>%
  row_spec(5, background = "#E8F5E9", color = "#1B5E20") %>%
  row_spec(6, background = "#E8F5E9", color = "#1B5E20") %>%
  row_spec(7,
           background = "#C8E6C9",
           color      = "#1B5E20",
           bold       = TRUE) %>%
  column_spec(1, bold = TRUE, width = "22%") %>%
  column_spec(2, bold = TRUE, width = "20%") %>%
  column_spec(3, width = "25%") %>%
  column_spec(4, italic = TRUE, width = "33%") %>%
  pack_rows("Chi-cuadrado de Pearson", 1, 3,
            label_row_css = "background:#90CAF9; color:#0D47A1;
                             font-weight:bold; font-size:12px;") %>%
  pack_rows("Prueba exacta de Fisher", 4, 7,
            label_row_css = "background:#A5D6A7; color:#1B5E20;
                             font-weight:bold; font-size:12px;") %>%
  footnote(
    general = paste0(
      "Nivel de significancia α = 0.05. ",
      "Un valor p > 0.05 indica que no existe evidencia estadística ",
      "para rechazar la hipótesis nula de independencia entre ",
      "tratamiento y extubación exitosa. ",
      "OR = 1 indica ausencia de asociación; ",
      "IC 95% que incluye el 1 me confirma que la diferencia no es ",
      "estadísticamente significativa."
    ),
    general_title = "Nota:",
    footnote_as_chunk = TRUE
  )
Tabla 3. Resultados de las pruebas de asociación sobre la Extubación exitosa vs. Tratamiento
Prueba Elemento Resultado Interpretación
Chi-cuadrado de Pearson
Chi-cuadrado de Pearson Estadístico X² 0.021446370457114 Valor muy cercano a 0 — diferencia mínima entre observado y esperado
Chi-cuadrado de Pearson Grados de libertad 1.000000000000000 Tabla 2×2 tiene 1 grado de libertad
Chi-cuadrado de Pearson Valor p 0.883569478909719 Muy superior a α = 0.05 — no se rechaza H₀
Prueba exacta de Fisher
Prueba exacta de Fisher Razón de momios (OR) 1.069327313217891 OR ≈ 1 — prácticamente sin diferencia entre grupos
Prueba exacta de Fisher Límite Inferior IC 95% 0.391715289165425 Límite inferior cruza el 1 — compatible con ausencia de asociación
Prueba exacta de Fisher Límite Superior IC 95% 2.944879122436857 Límite superior amplio — estimación imprecisa por grupos similares
Prueba exacta de Fisher Valor p 1.000000000000000 p = 1 — máxima compatibilidad con H₀ de independencia
Nota: Nivel de significancia α = 0.05. Un valor p > 0.05 indica que no existe evidencia estadística para rechazar la hipótesis nula de independencia entre tratamiento y extubación exitosa. OR = 1 indica ausencia de asociación; IC 95% que incluye el 1 me confirma que la diferencia no es estadísticamente significativa.

3.0.2 Análisis e interpretación — Pregunta 2

Chi-cuadrado de Pearson: X² = 0.0214446370457114: comienzo definiendo que un estadístico chi-cuadrado cercano a cero significa que la distancia total entre lo observado y lo esperado bajo independencia es casi nula,pues esta prueba suma, para cada celda, el cociente entre el cuadrado de la diferencia (observado − esperado) y el valor esperado.

Con diferencias tan pequeñas como las que he encontrado, como que ninguna celda difiere en más de 0.32 pacientes de su valor esperado, entonces el resultado es inevitablemente un estadístico prácticamente igual a cero.

Con 1 grado de libertad, propio de toda tabla 2×2, el valor p asociado es 0.883569478909719, muy por encima del nivel de significancia α = 0.05. Por tanto, no rechazo la hipótesis nula, en otras palabras, no existe evidencia estadística de asociación entre el tipo de tratamiento recibido y la probabilidad de extubación exitosa.

Razón de momios (OR) = 1.069327313217891: La razón de momios cuantifica cuántas veces más probable es la extubación exitosa en el grupo de Entrenamiento Muscular Respiratorio, respecto al grupo Convencional.

Por tanto, un OR de exactamente 1.0 indicaría que ambos grupos tienen probabilidades idénticas. El valor que se obtuvo de 1.069 es prácticamente indistinguible de 1, pues el grupo de Entrenamiento Muscular Respiratorio tiene apenas un 6.9% más de momios de extubación exitosa que el Convencional, una diferencia clínicamente irrelevante y estadísticamente nula.

Intervalo de confianza 95%: 0.391715 — 2.944879: Este intervalo es el elemento más informativo de mi análisis, pues el hecho, de que el valor 1 quede dentro del intervalo, y holgadamente, pues el límite inferior es 0.39 y el superior es 2.94 — confirma que el OR observado de 1.07 es perfectamente compatible con la ausencia de asociación. La amplitud del intervalo (de 0.39 a 2.94) también nos habla de imprecisión en la estimación, lo cual es esperable cuando los dos grupos tienen proporciones tan similares y los eventos adversos (no extubación) son relativamente pocos (23 en total).

Valor p de Fisher = 1.000000000000000

Un valor p exactamente igual a 1 es el resultado más contundente posible a favor de la hipótesis nula. Significa que la tabla que observamos — con 52 y 51 extubaciones exitosas en cada grupo — es precisamente la que esperaríamos bajo independencia perfecta. No hay ninguna dirección en la que los datos se alejen de H₀.

3.0.3 Conclusión Pregunta 2

No existe evidencia estadística de asociación entre el tipo de tratamiento recibido y la extubación exitosa (X² = 0.0214, p = 0.884; OR = 1.069, Intervalo de Confianza 95%: 0.392 — 2.945,p Fisher = 1.000).

Los resultados son consistentes con los reportados por Sandoval et al. (2019), quienes tampoco encontraron diferencias significativas en la probabilidad de extubación entre grupos, por tanto, el Entrenamiento Muscular Respiratorio no modificó la probabilidad de extubación exitosa respecto al cuidado estándar en esta población.

3.1 Pregunta 3: ¿Cuál es el promedio y la mediana de horas de destete

de la ventilación mecánica?

Para responder esta pregunta calcularé las medidas de tendencia central de la variable Horasdedestete, tanto para la muestra total como estratificada por grupo de tratamiento.

Ahora bien, si reportara únicamente el promedio, estaría cometiendo un error metodológico importante, pues los tiempos hasta un evento , como lo es el destete de la ventilación mecánica, casi nunca se distribuyen de forma simétrica en la práctica clínica.

Lo que ocurre típicamente es que la mayoría de los pacientes se extuba en tiempos relativamente cortos, pero un pequeño grupo de casos complejos puede permanecer días en destete, generando una distribución con cola derecha en la que esos pocos valores extremos jalan la media hacia arriba.

Entonces, la media termina representando un paciente que en realidad no existe en la muestra , ni el típico ni el extremo , lo cual me llevaría a distorsionar la interpretación clínica.

Por esa razón, la mediana es la medida de tendencia central más apropiada aquí, pues al ser el valor que divide exactamente la distribución en dos mitades iguales, no se ve afectada por esos tiempos extremos y refleja con fidelidad la experiencia del paciente real.

La complementaré con el rango intercuartílico, que me dice dentro de qué intervalo de horas se concentra el 50% central de los pacientes, dándome una imagen mucho más honesta de la dispersión.

Finalmente, verificaré si la diferencia entre grupos es estadísticamente significativa mediante la prueba de Mann-Whitney, que es la prueba no paramétrica adecuada para comparar medianas entre dos grupos independientes cuando no puedo asumir distribución normal , que precisamente es el caso aquí.

options(digits = 15)

# ── Muestra total ──────────────────────────────────────────────────────
cat("=== MUESTRA TOTAL ===\n")
## === MUESTRA TOTAL ===
cat("n            :", nrow(datos), "\n")
## n            : 126
cat("Media        :", mean(datos$tiempo, na.rm = TRUE), "\n")
## Media        : 9.06904761904762
cat("Mediana      :", median(datos$tiempo, na.rm = TRUE), "\n")
## Mediana      : 5
cat("Desv. estándar:", sd(datos$tiempo, na.rm = TRUE), "\n")
## Desv. estándar: 11.9234522805148
cat("P25          :", quantile(datos$tiempo, 0.25, na.rm = TRUE), "\n")
## P25          : 2
cat("P75          :", quantile(datos$tiempo, 0.75, na.rm = TRUE), "\n")
## P75          : 11
cat("Mínimo       :", min(datos$tiempo, na.rm = TRUE), "\n")
## Mínimo       : 0
cat("Máximo       :", max(datos$tiempo, na.rm = TRUE), "\n")
## Máximo       : 72
# ── Por grupo ──────────────────────────────────────────────────────────
cat("\n=== POR GRUPO DE TRATAMIENTO ===\n")
## 
## === POR GRUPO DE TRATAMIENTO ===
datos %>%
  group_by(trat_f) %>%
  summarise(
    n             = n(),
    Media         = mean(tiempo, na.rm = TRUE),
    Mediana       = median(tiempo, na.rm = TRUE),
    DE            = sd(tiempo, na.rm = TRUE),
    P25           = quantile(tiempo, 0.25, na.rm = TRUE),
    P75           = quantile(tiempo, 0.75, na.rm = TRUE),
    Minimo        = min(tiempo, na.rm = TRUE),
    Maximo        = max(tiempo, na.rm = TRUE)
  ) %>% print()
## # A tibble: 2 × 9
##   trat_f           n Media Mediana    DE   P25   P75 Minimo Maximo
##   <fct>        <int> <dbl>   <dbl> <dbl> <dbl> <dbl>  <dbl>  <dbl>
## 1 Convencional    64  8.78       5  11.4     2  9.25      0     59
## 2 EMR             62  9.36       5  12.5     2 11.8       0     72
# ── Aplico Prueba de Mann-Whitney ─────────────────────────────────────────────
cat("\n=== PRUEBA MANN-WHITNEY ===\n")
## 
## === PRUEBA MANN-WHITNEY ===
mw <- wilcox.test(tiempo ~ trat_f, data = datos)
cat("Estadístico W:", mw$statistic, "\n")
## Estadístico W: 1943.5
cat("Valor p      :", mw$p.value, "\n")
## Valor p      : 0.844697691815235

3.1.1 Mi análisis

Me parece importante, analizar la relación entre la media y la mediana, pues según los datos, para la muestra total, la media es 9.069 horas mientras que la mediana es 5 horas,es decir, la media prácticamente duplica la mediana, por tanto, esa brecha es la huella digital de la asimetría hacia la derecha que anticipaba, pues la mayoría de los pacientes se desteta en 5 horas o menos, pero un grupo de casos con tiempos prolongados, incluso llegando hasta 72 horas, me jala la media hacia arriba, alejándola del valor típico real.

El rango intercuartílico de 2 a 11 horas confirma que el 50% central de los pacientes se desteta en ese intervalo, mientras la desviación estándar de 11.92 horas, mayor que la propia media, evidencia una dispersión extrema incompatible con una distribución normal simétrica.

Por grupo, la situación es idéntica,pues ambos tienen mediana de 5 horas, con medias de 8.78 horas (Convencional) y 9.36 horas (EMR),y en cuanto a la prueba de Mann-Whitneyesta me confirma, que la diferencia de 0.58 horas no es estadísticamente significativa (W = 1943.5, p = 0.845).

library(dplyr)
library(ggplot2)
library(kableExtra)

# ══════════════════════════════════════════════════════
# TABLA RESUMEN
# ══════════════════════════════════════════════════════
resumen <- data.frame(
  `Estadístico` = c(
    "n", "Media (horas)", "Mediana (horas)",
    "Desviación estándar", "P25 — P75 (RIQ)",
    "Mínimo", "Máximo",
    "Diferencia media − mediana",
    "Mann-Whitney W", "Valor p"
  ),
  `Muestra total` = c(
    "126",
    "9.06904761904762",
    "5",
    "11.9234522805148",
    "2 — 11",
    "0", "72",
    "4.069 horas 🔴 asimetría",
    "1943.5",
    "0.844697691815235"
  ),
  `Convencional (n=64)` = c(
    "64",
    "8.78437500000000",
    "5",
    "11.4180755995406",
    "2 — 9.25",
    "0", "59",
    "3.784 horas",
    "—", "—"
  ),
  `EMR (n=62)` = c(
    "62",
    "9.36290322580645",
    "5",
    "12.5103790860120",
    "2 — 11.75",
    "0", "72",
    "4.363 horas",
    "—", "—"
  ),
  check.names = FALSE
)

resumen %>%
  kbl(
    caption = "Tabla 4. Estadísticas descriptivas de horas de destete ( total y por grupo)",
    align   = c("l","c","c","c")
  ) %>%
  kable_styling(
    bootstrap_options = c("hover","bordered","condensed"),
    full_width        = TRUE,
    font_size         = 13
  ) %>%
  row_spec(0, bold = TRUE, color = "white", background = "#4A235A") %>%
  row_spec(1, background = "#EDE7F6", color = "#4A235A") %>%
  row_spec(2,
           background = "#FCE4EC", color = "#880E4F",
           bold = TRUE) %>%   # media — afectada por asimetría
  row_spec(3,
           background = "#D5F5E3", color = "#1B5E20",
           bold = TRUE) %>%   # mediana — valor real
  row_spec(4, background = "#EDE7F6", color = "#4A235A") %>%
  row_spec(5, background = "#E3F2FD", color = "#0D47A1") %>%
  row_spec(6:7, background = "#FFF9C4", color = "#795548") %>%
  row_spec(8,
           background = "#FFCDD2", color = "#B71C1C",
           bold = TRUE) %>%   # diferencia — evidencia asimetría
  row_spec(9:10, background = "#E8F5E9", color = "#1B5E20") %>%
  column_spec(1, bold = TRUE, width = "28%") %>%
  pack_rows("Tendencia central", 2, 3,
            label_row_css = "background:#CE93D8; color:#4A235A;
                             font-weight:bold; font-size:12px;") %>%
  pack_rows("Dispersión", 4, 7,
            label_row_css = "background:#90CAF9; color:#0D47A1;
                             font-weight:bold; font-size:12px;") %>%
  pack_rows("Evidencia de asimetría", 8, 8,
            label_row_css = "background:#EF9A9A; color:#B71C1C;
                             font-weight:bold; font-size:12px;") %>%
  pack_rows("Prueba comparativa", 9, 10,
            label_row_css = "background:#A5D6A7; color:#1B5E20;
                             font-weight:bold; font-size:12px;") %>%
  footnote(
    general = "La diferencia entre media y mediana es el indicador más claro de asimetría, pues cuando la media supera la mediana, la distribución tiene cola derecha, pues pocos tiempos prolongados elevan la media sin alterar la mediana.",
    general_title = "Nota clave:",
    footnote_as_chunk = TRUE
  )
Tabla 4. Estadísticas descriptivas de horas de destete ( total y por grupo)
Estadístico Muestra total Convencional (n=64) EMR (n=62)
n 126 64 62
Tendencia central
Media (horas) 9.06904761904762 8.78437500000000 9.36290322580645
Mediana (horas) 5 5 5
Dispersión
Desviación estándar 11.9234522805148 11.4180755995406 12.5103790860120
P25 — P75 (RIQ) 2 — 11 2 — 9.25 2 — 11.75
Mínimo 0 0 0
Máximo 72 59 72
Evidencia de asimetría
Diferencia media − mediana 4.069 horas 🔴 asimetría 3.784 horas 4.363 horas
Prueba comparativa
Mann-Whitney W 1943.5
Valor p 0.844697691815235
Nota clave: La diferencia entre media y mediana es el indicador más claro de asimetría, pues cuando la media supera la mediana, la distribución tiene cola derecha, pues pocos tiempos prolongados elevan la media sin alterar la mediana.
# mis cálculos de media y mediana por grupo para las líneas
lineas <- datos %>%
  filter(!is.na(tiempo)) %>%
  group_by(trat_f) %>%
  summarise(
    media   = mean(tiempo),
    mediana = median(tiempo)
  )

colores_fill   <- c("Convencional" = "#B39DDB", "EMR" = "#80CBC4")
colores_border <- c("Convencional" = "#7B1FA2", "EMR" = "#00796B")

datos %>%
  filter(!is.na(tiempo)) %>%
  ggplot(aes(x = tiempo, fill = trat_f, color = trat_f)) +

  # Histograma
  geom_histogram(
    aes(y = after_stat(density)),
    binwidth   = 5,
    alpha      = 0.55,
    position   = "identity"
  ) +

  # Curva de densidad 
  geom_density(alpha = 0.15, linewidth = 0.8) +

  # Líneas verticales — mediana (sólida) y media (punteada)
  geom_vline(
    data     = lineas,
    aes(xintercept = mediana, color = trat_f),
    linetype = "solid",
    linewidth = 1.2
  ) +
  geom_vline(
    data     = lineas,
    aes(xintercept = media, color = trat_f),
    linetype = "dashed",
    linewidth = 1.0
  ) +

  # Hago las anotaciones de media y mediana
  annotate("text", x = 5.5,  y = 0.115,
           label = "Mediana = 5 h\n(ambos grupos)",
           size = 3.2, color = "#4A235A", fontface = "bold") +
  annotate("text", x = 16, y = 0.100,
           label = "Media ≈ 9 h\n(cola derecha)",
           size = 3.2, color = "#B71C1C", fontface = "italic") +
  annotate("segment",
           x = 5, xend = 14, y = 0.092, yend = 0.092,
           arrow = arrow(length = unit(0.2,"cm")),
           color = "#B71C1C", linewidth = 0.7) +

  scale_fill_manual(values  = colores_fill)  +
  scale_color_manual(values = colores_border) +

  scale_x_continuous(
    breaks = seq(0, 75, by = 10),
    labels = paste0(seq(0, 75, by = 10), " h")
  ) +

  labs(
    title    = "Distribución de horas de destete por grupo de tratamiento",
    subtitle = "La brecha entre media y mediana evidencia la asimetría con cola derecha",
    x        = "Horas de destete de la ventilación mecánica",
    y        = "Densidad",
    fill     = "Grupo",
    color    = "Grupo",
    caption  = "Línea sólida = mediana | Línea punteada = media | Fuente: Sandoval et al. (2019)"
  ) +

  theme_minimal(base_size = 13) +
  theme(
    plot.title      = element_text(face = "bold", color = "#4A235A", size = 14),
    plot.subtitle   = element_text(color = "#880E4F", size = 11, face = "italic"),
    plot.caption    = element_text(color = "gray50", size = 9),
    legend.position = "top",
    legend.title    = element_text(face = "bold"),
    panel.grid.minor = element_blank(),
    panel.grid.major = element_line(color = "gray90")
  )

¿Qué estoy viendo en esta gráfica?

Esta figura me permite ver algo que los números solos no transmiten con la misma claridad, como es la forma real de la distribución del tiempo de destete, paso a analizar cada elemento:

Las barras del histograma representan cuántos pacientes se destetaron en cada rango de horas , el morado corresponde al grupo Convencional y el verde al grupo de Entrenamiento Muscular Respiratorio.

Inmediatamente noto que la gran mayoría de los pacientes se concentra entre 0 y 10 horas, con una masa muy alta al inicio que cae bruscamente, por tanto, eso ya me dice que el destete fue rápido para la mayoría.

Las curvas de densidad suavizada confirman ese patrón, pues ambas curvas forman una joroba pronunciada cerca de las 5 horas y luego descienden de forma asimétrica, con una cola que se extiende hacia la derecha hasta las 70 horas, por lo cual, esa cola larga hacia la derecha es precisamente la distribución asimétrica positiva que ya anticipaba, cuando sucede que hay pocos pacientes pero con tiempos muy prolongados.

La línea sólida vertical marca la mediana de 5 horas, que es igual en ambos grupos, de hechio, noto que está ubicada exactamente donde se concentra la masa de datos, deduzco que ahí vive el paciente típico real.

Las líneas punteadas verticales marcan las medias, por ejemplo 8.78 horas para Convencional y 9.36 horas para el grupo de Entrenamiento Muscular Respiratorio, donde se aprecia que están claramente desplazadas hacia la derecha respecto a la mediana, flotando en una zona donde ya hay muy pocos pacientes.

Esa brecha visual entre la línea sólida y las líneas punteadas es exactamente la asimetría, pues la media está siendo jalada por los casos extremos de la cola derecha, alejándose del valor que representa a la mayoría.

La flecha roja señala precisamente esa distancia,que vá desde la mediana real hasta donde la media termina ubicándose. Si la distribución fuera simétrica, ambas líneas coincidirían, además, el hecho de que estén separadas por aproximadamente 4 horas me confirma que reportar solo la media habría sobreestimado el tiempo de destete del paciente típico, llevándome a una conclusión distorsionada sobre la eficacia del tratamiento.

Finalmente, que ambos grupos tengan prácticamente la misma forma de distribución, la misma mediana y curvas casi superpuestas, anticipa visualmente lo que la prueba de Mann-Whitney confirmó formalmente, pues no existe diferencia estadísticamente significativa entre los grupos en el tiempo de destete (W = 1943.5, p = 0.845).

3.2 Pregunta 4: Análisis de supervivencia para toda la muestra y por cada grupo de tratamiento

Hasta ahora he trabajado con proporciones y promedios, pero esas medidas tienen una limitación fundamental en este contexto, pues no incorporan el tiempo.

Saber que el 81.7% de los pacientes fue extubados exitosamente no me dice nada sobre cuándo ocurrió esa extubación ni cómo fue evolucionando esa probabilidad hora a hora, para eso debo hacer el análisis de supervivencia.

El método que utilizaré es el estimador de Kaplan-Meier, que es el estándar para este tipo de análisis, donde su lógica es que, en lugar de calcular una proporción global, me construye una función escalonada que va actualizando la probabilidad de no haber sido extubado aún en cada momento en que ocurre un evento.

Cada vez que un paciente es extubado, la curva desciende un escalón. Los pacientes censurados, es decir, aquellos que terminaron el seguimiento sin ser extubados, no se descartan sino que se incorporan al denominador hasta el momento en que salen del análisis, lo cual es metodológicamente correcto y evita subestimar la probabilidad de supervivencia.

Realizaré el análisis en dos niveles, primero para la muestra total, que me dará la función de supervivencia global del ensayo, luego estratificaré por grupo de tratamiento, que me permitirá ver visualmente si las curvas de Convencional y EMR se separan o se superponen, lo que me anticiparía si existe diferencia entre grupos antes de aplicar ninguna prueba formal.

library(survival)
library(survminer)

# ── Objeto de supervivencia ────────────────────────────────────────────
surv_obj <- Surv(time  = datos$tiempo,
                 event = datos$evento)

# ── KM muestra total ───────────────────────────────────────────────────
km_total <- survfit(surv_obj ~ 1, data = datos)

cat("=== KAPLAN-MEIER — MUESTRA TOTAL ===\n")
## === KAPLAN-MEIER — MUESTRA TOTAL ===
cat("n total         :", km_total$n, "\n")
## n total         : 126
cat("n eventos       :", km_total$n.event %>% sum(), "\n")
## n eventos       : 103
cat("n censurados    :", km_total$n.censor %>% sum(), "\n")
## n censurados    : 23
cat("Mediana (horas) :", surv_median(km_total)$median, "\n")
## Mediana (horas) : 5
cat("IC 95% LI       :", surv_median(km_total)$lower, "\n")
## IC 95% LI       : 4
cat("IC 95% LS       :", surv_median(km_total)$upper, "\n")
## IC 95% LS       : 7
print(km_total)
## Call: survfit(formula = surv_obj ~ 1, data = datos)
## 
##        n events median 0.95LCL 0.95UCL
## [1,] 126    103      5       4       7

3.2.1 Interpretación — Kaplan-Meier muestra total

Los resultados del estimador de Kaplan-Meier para la muestra completa me entregan información muy precisa sobre la dinámica del destete, pues de los 126 pacientes, 103 alcanzaron el evento (extubación exitosa) y 23 quedaron censurados, es decir, salieron del seguimiento sin haber sido extubados, ya sea por muerte, traqueostomía u otras causas.

La mediana de supervivencia es de 5 horas (Intervalo de Confianza 95%: 4 — 7 horas), lo cual, significa que el 50% de los pacientes logró ser extubado exitosamente antes de las 5 horas desde el inicio del destete, y dicho de otra forma, si tomara un paciente aleatorio de esta población, la probabilidad de que aún no haya sido extubado a las 5 horas es exactamente del 50%.

El intervalo de confianza de 4 a 7 horas es relativamente estrecho, lo que me indica que la estimación de la mediana es precisa y confiable. Si replicara este estudio en poblaciones similares, en el 95% de los casos la mediana verdadera se ubicaría dentro de ese rango.

Esto concuerda exactamente con lo reportado por Sandoval et al.  (2019), quienes reportaron una mediana de destete de 5 horas en ambos grupos, validando así la integridad del análisis que se presenta.

ggsurvplot(
  km_total,
  data          = datos,
 
  conf.int      = TRUE,
  conf.int.fill = "#B39DDB",
  conf.int.alpha= 0.25,
  surv.median.line = "hv",       # líneas horizontal+vertical en mediana
  # ── Tabla de riesgo ────────────────────────────────────────────────
  risk.table         = TRUE,
  risk.table.height  = 0.25,
  risk.table.col     = "black",
  risk.table.title   = "Pacientes en riesgo",
  # ── Estética ───────────────────────────────────────────────────────
  palette       = "#7B1FA2",
  size          = 1.0,
  censor.shape  = "|",
  censor.size   = 4,
  # ── Ejes y etiquetas ───────────────────────────────────────────────
  xlab          = "Horas de destete de la ventilación mecánica",
  ylab          = "Probabilidad de no haber sido extubado",
  title         = "Curva de Kaplan-Meier — Muestra total (n = 126)",
  subtitle      = "Mediana de supervivencia: 5 horas (IC 95%: 4 — 7 h)",
  legend        = "none",
  break.time.by = 10,
  xlim          = c(0, 75),
  # ── Tema ───────────────────────────────────────────────────────────
  ggtheme = theme_minimal(base_size = 12) +
    theme(
      plot.title    = element_text(face = "bold",
                                   color = "#4A235A", size = 13),
      plot.subtitle = element_text(color = "#880E4F",
                                   face  = "italic", size = 10),
      axis.title    = element_text(color = "#37474F"),
      panel.grid.minor = element_blank()
    )
)

¿Qué estoy viendo en esta curva?

Esta curva de Kaplan-Meier es una de las representaciones más informativas del análisis de supervivencia, ya que se observa la representación en el eje vertical de la probabilidad de que un paciente todavía no haya sido extubado en cada momento del seguimiento, arrancando en 1.0 (100%) porque al inicioningún paciente ha sido extubado aún, y desciende progresivamente a medida que ocurren los eventos.

En el eje horizontal se representa el tiempo transcurrido en horas desde el inicio del destete.La forma escalonada de la curva es característica del estimador de Kaplan-Meier, pues cada escalón hacia abajo corresponde exactamente al momento en que uno o más pacientes fueron extubados, por tanto, entre escalones, la curva es horizontal porque no ocurrió ningún evento nuevo.

Los pequeños palitos verticales (+) sobre la curva marcan los pacientes censurados, es decir, aquellos que salieron del seguimiento sin ser extubados, donde además,se observa que no bajan la curva pero sí reducen el denominador para los cálculos siguientes.

La banda sombreada lila es el intervalo de confianza al 95%, donde noto que es relativamente estrecha en las primeras horas, donde hay muchos pacientes en riesgo y la estimación es precisa, y se ensancha progresivamente hacia la derecha, donde quedan pocos pacientes y la incertidumbre crece naturalmente.

Las líneas punteadas marcan la mediana de supervivencia, donde la línea horizontal en 0.50 y la vertical en 5 horas se cruzan exactamente en la curva, confirmando que a las 5 horas el 50% de los pacientes aún no había sido extubado, o dicho de otro modo, el otro 50% ya lo había sido, y ese es el tiempo de destete mediano.

La caída más pronunciada ocurre en las primeras 10 horas, donde la curva desciende desde 1.0 hasta aproximadamente 0.25, lo cual, me confirma visualmente lo que se había observado, de que la gran mayoría de los destetes ocurrieron temprano, y luego la curva se aplana progresivamente con una cola que se extiende hasta las 72 horas, consistente con la distribución asimétrica observada.

# ── KM estratificado por tratamiento ──────────────────────────────────
km_grupos <- survfit(surv_obj ~ trat_f, data = datos)

cat("=== KAPLAN-MEIER POR GRUPO ===\n")
## === KAPLAN-MEIER POR GRUPO ===
print(km_grupos)
## Call: survfit(formula = surv_obj ~ trat_f, data = datos)
## 
##                      n events median 0.95LCL 0.95UCL
## trat_f=Convencional 64     52      5       4       8
## trat_f=EMR          62     51      6       4      11
cat("\n=== MEDIANAS POR GRUPO ===\n")
## 
## === MEDIANAS POR GRUPO ===
print(surv_median(km_grupos))
##                strata median lower upper
## 1 trat_f=Convencional      5     4     8
## 2          trat_f=EMR      6     4    11

3.2.2 Interpretación Kaplan-Meier por grupo de tratamiento

Los resultados estratificados revelan un panorama muy consistente con lo que he venido encontrando a lo largo del análisis, pues en el grupo Convencional, conformado por 64 pacientes, presentó 52 eventos de extubación exitosa, con una mediana de supervivencia de 5 horas, cuyo intervalo de confianza al 95% se ubica entre 4 y 8 horas.

El grupo de Entrenamiento Muscular Respiratorio, conformado por 62 pacientes, presentó 51 eventos, con una mediana de supervivencia de 6 horas, cuyo intervalo de confianza al 95% se ubica entre 4 y 11 horas.

La diferencia entre medianas es de apenas 1 hora a favor del grupo Convencional, lo cual desde una perspectiva clínica es poco relevante,sin embargo, lo que más me llama la atención no es la diferencia en las medianas sino la diferencia en los intervalos de confianza, pues el grupo Convencional tiene un intervalo de 4 horas de amplitud, mientras que el grupo de Entrenamiento Muscular Respiratorio tiene un intervalo de 7 horas de amplitud, lo que me indica que la estimación de la mediana es menos precisa en el grupo experimental.

Eso me parece que es esperable porque la distribución del tiempo de destete en ese grupo tiene una cola derecha más larga, llegando hasta 72 horas frente a 59 horas del Convencional,generando así mayor variabilidad en la estimación.

Estos resultados coinciden exactamente con los reportados por Sandoval et al. (2019), quienes encontraron medianas de 5 horas en ambos grupos, validando así la reproducibilidad del análisis.

ggsurvplot(
  km_grupos,
  data     = datos,
 
  conf.int         = TRUE,
  conf.int.alpha   = 0.15,
  surv.median.line = "hv",
  # Tabla de riesgo 
  risk.table        = TRUE,
  risk.table.height = 0.28,
  risk.table.title  = "Pacientes en riesgo",
  risk.table.col    = "strata",
 
  palette      = c("#7B1FA2", "#00796B"),
  size         = 1.0,
  censor.shape = "|",
  censor.size  = 4,
  linetype     = c("solid", "dashed"),
  #  Leyenda 
  legend.title  = "Grupo",
  legend.labs   = c("Convencional", "EMR"),
  legend        = "top",
  #  Ejes y etiquetas 
  xlab          = "Horas de destete de la ventilación mecánica",
  ylab          = "Probabilidad de no haber sido extubado",
  title         = "Curva de Kaplan-Meier por grupo de tratamiento",
  subtitle      = "Convencional: mediana 5 h (IC 95%: 4–8 h)  |  EMR: mediana 6 h (IC 95%: 4–11 h)",
  break.time.by = 10,
  xlim          = c(0, 75),
 
  ggtheme = theme_minimal(base_size = 12) +
    theme(
      plot.title    = element_text(face  = "bold",
                                   color = "#4A235A", size = 13),
      plot.subtitle = element_text(color = "#880E4F",
                                   face  = "italic", size = 10),
      axis.title    = element_text(color = "#37474F"),
      legend.text   = element_text(size  = 11),
      panel.grid.minor = element_blank()
    )
)

Interpretación de la curva de Kaplan-Meier por grupo de tratamiento

La figura muestra las funciones de supervivencia estratificadas por grupo, donde la línea lila representa al Convencional y la línea verde menta al grupo de Entrenamiento Muscular Respiratorio.

A lo largo del seguimiento, se observa que las dos curvas se entrelazan de manera casi idéntica, pues comparten una caída pronunciada inicial donde descienden desde 1.0 hasta aproximadamente 0.25durante las primeras diez horas, y luego se aplanan progresivamente con una cola que se extiende hasta las 72 horas en el grupo EMR.

El grupo Convencional alcanzó una mediana de 5 horas (Intervalo de Confianza 95%: 4 a 8 horas), mientras que el grupo EMR registró una mediana de 6 horas (Intervalo de Confianza 95%: 4 a 11 horas).

Existe, aparentemente, una diferencia de una hora a favor del grupo Convencional, no obstante, esa separación se diluye al notar que los intervalos de confianza se solapan casi por completo, de modo que el desfase considero que carece de precisión estadística.

Prueba de ello es que la amplitud del intervalo de EMR (7 horas) casi duplica la del Convencional (4 horas), asimetría que revela mayor imprecisión en el grupo experimental, consecuencia de los casos atípicos con destetes prolongados que arrastran la cola derecha.

Asimismo, las bandas de confianza sombreadas se solapan a lo largo de todo el eje temporal, sin ningún tramo donde una curva quede fuera del intervalo de la otra; y en cuanto a las marcas de censura, representadas por los trazos verticales sobre cada curva,estas se distribuyen de manera similar en ambos grupos, respaldando el supuesto de censura no informativa y descartando pérdidas diferenciales de seguimiento entre tratamientos.

En este sentido, la validez del estimador de Kaplan-Meier queda asegurada para esta comparación, pues los dos grupos pierden pacientes por razones similares y en tiempos parecidos, sin sesgo diferencial que contamine la estimación.

La tabla de pacientes en riesgo me muestra el descenso progresivo de la cohorte, pues comienza con 64 y 62 pacientes al inicio, y para las 40 horas apenas quedan entre 1 y 2 pacientes por grupo. Consecuentemente, la estimación es sólida en las primeras horas, donde se concentra la gran mayoría de los eventos, y se vuelve más incierta hacia el final, razón por la cual las bandas de confianza se ensanchan notoriamente en la cola derecha.

A la luz de estos resultados, este hallazgo coincide plenamente con lo reportado por Sandoval et al. (2019), quienes también obtuvieron medianas de destete de 5 horas en ambos grupos y concluyeron que la probabilidad de extubación no se asociaba con el tratamiento administrado (log-rank = 1.12, p = 0.29).

Lo anterior se vincula con la explicación que ofrecen los autores en la discusión, pues la mayoría de sus pacientes recibió menos de catorce sesiones de entrenamiento, tiempo que consideran insuficiente para inducir los cambios fisiológicos reales en la fibra muscular inspiratoria.

4 Pregunta 5

4.1 Comparación de las funciones de supervivencia entre grupos de tratamiento

En la Pregunta 4 observé que las curvas de Kaplan-Meier del grupo Convencional y del grupo EMR se entrelazan de manera casi idéntica, no obstante, esa impresión visual necesita una prueba formal que la respalde.

Para ello aplico la prueba log-rank (Mantel-Cox), que es el estándar epidemiológico para comparar dos funciones de supervivencia completas, pues a diferencia del chi-cuadrado clásico, esta prueba sí incorpora el tiempo y pondera cada momento del seguimiento según los pacientes en riesgo, donde su lógica es comparar, en cada instante en que ocurre una extubación, cuántos eventos observo en cada grupo contra cuántos esperaría bajo la hipótesis nula de igualdad de hazards.

Si la discrepancia acumulada es pequeña, el chi² será cercano a cero y el valor p grande, por tanto concluiré que las curvas son estadísticamente equivalentes.

library(kableExtra)

hipotesis <- data.frame(
  `Elemento` = c(
    "Hipótesis nula (H₀)",
    "Hipótesis alterna (H₁)",
    "Nivel de significancia (α)",
    "Regla de decisión"
  ),
  `Planteamiento` = c(
    "$S_{\\text{Convencional}}(t) = S_{\\text{EMR}}(t)$",
    "$S_{\\text{Convencional}}(t) \\neq S_{\\text{EMR}}(t)$",
    "$\\alpha = 0.05$",
  "Rechazo H₀ cuando el valor p es menor a 0.05"
  ),
  `Interpretación clínica` = c(
    "Las funciones de supervivencia son iguales en ambos grupos, por tanto el tratamiento no modifica la dinámica del destete.",
    "Las funciones de supervivencia difieren en al menos un momento del seguimiento, por tanto el tratamiento sí modifica el tiempo hasta la extubación.",
    "Tolero un 5% de probabilidad de rechazar H₀ cuando en realidad es verdadera (error tipo I).",
    "Si el valor p supera 0.05, no tengo evidencia suficiente para afirmar diferencia entre curvas."
  ),
  check.names = FALSE
)

hipotesis %>%
  kbl(
    caption = "Tabla 5a. Planteamiento de hipótesis para la prueba log-rank",
    align   = c("l","c","l"),
    escape  = FALSE
  ) %>%
  kable_styling(
    bootstrap_options = c("hover","bordered","condensed"),
    full_width = TRUE,
    font_size  = 13,
    position   = "center"
  ) %>%
  row_spec(0, bold = TRUE, color = "white", background = "#4A235A") %>%
  row_spec(1, background = "#EDE7F6", color = "#4A235A") %>%
  row_spec(2, background = "#FCE4EC", color = "#880E4F") %>%
  row_spec(3, background = "#FFE0B2", color = "#BF360C") %>%
  row_spec(4, background = "#E0F2F1", color = "#00695C", bold = TRUE) %>%
  column_spec(1, bold = TRUE, width = "25%") %>%
  column_spec(2, width = "25%") %>%
  column_spec(3, italic = TRUE, width = "50%") %>%
  footnote(
    general = paste0(
      "La prueba log-rank compara las funciones de supervivencia completas ",
      "incorporando el tiempo y ponderando cada momento del seguimiento ",
      "según los pacientes en riesgo."
    ),
    general_title = "Nota:",
    footnote_as_chunk = TRUE
  )
Tabla 5a. Planteamiento de hipótesis para la prueba log-rank
Elemento Planteamiento Interpretación clínica
Hipótesis nula (H₀) \(S_{\text{Convencional}}(t) = S_{\text{EMR}}(t)\) Las funciones de supervivencia son iguales en ambos grupos, por tanto el tratamiento no modifica la dinámica del destete.
Hipótesis alterna (H₁) \(S_{\text{Convencional}}(t) \neq S_{\text{EMR}}(t)\) Las funciones de supervivencia difieren en al menos un momento del seguimiento, por tanto el tratamiento sí modifica el tiempo hasta la extubación.
Nivel de significancia (α) \(\alpha = 0.05\) Tolero un 5% de probabilidad de rechazar H₀ cuando en realidad es verdadera (error tipo I).
Regla de decisión Rechazo H₀ cuando el valor p es menor a 0.05 Si el valor p supera 0.05, no tengo evidencia suficiente para afirmar diferencia entre curvas.
Nota: La prueba log-rank compara las funciones de supervivencia completas incorporando el tiempo y ponderando cada momento del seguimiento según los pacientes en riesgo.

4.2 Cálculo de la prueba log-rank

options(digits = 15)

log_rank <- survdiff(
  Surv(tiempo, evento) ~ trat_f,
  data = datos,
  rho  = 0
)

print(log_rank)
## Call:
## survdiff(formula = Surv(tiempo, evento) ~ trat_f, data = datos, 
##     rho = 0)
## 
##                      N Observed     Expected      (O-E)^2/E      (O-E)^2/V
## trat_f=Convencional 64       52 50.954851202 0.021437330982 0.047528396976
## trat_f=EMR          62       51 52.045148798 0.020988238781 0.047528396976
## 
##  Chisq= 0  on 1 degrees of freedom, p= 0.827421254
chi_lr <- log_rank$chisq
gl_lr  <- length(log_rank$n) - 1
p_lr   <- 1 - pchisq(chi_lr, df = gl_lr)

cat("Chi² :", chi_lr, "\n")
## Chi² : 0.0475283969759163
cat("gl   :", gl_lr, "\n")
## gl   : 1
cat("p    :", p_lr, "\n")
## p    : 0.827421254051622

4.3 Tabla resumen de la prueba log-rank

library(kableExtra)

tabla_lr <- data.frame(
  Grupo = c("Convencional", "EMR", "Resultado global"),
  `N`   = c("64", "62", "126"),
  `Observados (O)` = c("52", "51", "103"),
  `Esperados (E)`  = c(
    "50.9548512018088",
    "52.0451487981912",
    "103.0000000000000"
  ),
  `(O − E)` = c(
    "+1.0451487981912",
    "−1.0451487981912",
    "—"
  ),
  `Estadístico` = c("—", "—", "χ² = 0.0475 | gl = 1 | p = 0.8274"),
  check.names = FALSE
)

tabla_lr %>%
  kbl(
    caption = "Tabla 5. Prueba log-rank sobre eventos observados y esperados por grupo",
    align   = c("l","c","c","c","c","c"),
    escape  = FALSE
  ) %>%
  kable_styling(
    bootstrap_options = c("hover","bordered","condensed"),
    full_width = TRUE,
    font_size  = 13
  ) %>%
  row_spec(0, bold = TRUE, color = "white", background = "#4A235A") %>%
  row_spec(1, background = "#EDE7F6", color = "#4A235A") %>%
  row_spec(2, background = "#E0F2F1", color = "#00695C") %>%
  row_spec(3, background = "#FFE0B2", color = "#BF360C", bold = TRUE) %>%
  column_spec(1, bold = TRUE) %>%
  footnote(
    general = paste0(
      "La diferencia (O − E) suma cero entre grupos porque el total observado ",
      "coincide exactamente con el total esperado bajo H₀. ",
      "Una diferencia cercana a cero evidencia ausencia de asociación entre ",
      "tratamiento y tiempo de destete."
    ),
    general_title = "Nota:",
    footnote_as_chunk = TRUE
  )
Tabla 5. Prueba log-rank sobre eventos observados y esperados por grupo
Grupo N Observados (O) Esperados (E) (O − E) Estadístico
Convencional 64 52 50.9548512018088 +1.0451487981912
EMR 62 51 52.0451487981912 −1.0451487981912
Resultado global 126 103 103.0000000000000 χ² = 0.0475 | gl = 1 | p = 0.8274
Nota: La diferencia (O − E) suma cero entre grupos porque el total observado coincide exactamente con el total esperado bajo H₀. Una diferencia cercana a cero evidencia ausencia de asociación entre tratamiento y tiempo de destete.

4.4 Gráfica de las curvas con el valor p del log-rank

library(survminer)

km_grupos <- survfit(Surv(tiempo, evento) ~ trat_f, data = datos)

ggsurvplot(
  km_grupos,
  data             = datos,
  conf.int         = TRUE,
  conf.int.alpha   = 0.18,
  pval             = TRUE,
  pval.method      = TRUE,
  pval.coord       = c(45, 0.78),
  pval.method.coord= c(45, 0.85),
  pval.size        = 4.2,
  
  surv.median.line = "hv",
  
  risk.table        = TRUE,
  risk.table.height = 0.28,
  risk.table.title  = "Pacientes en riesgo",
  risk.table.col    = "strata",
  
  palette      = c("#B39DDB", "#80CBC4"),
  size         = 1.1,
  linetype     = c("solid", "dashed"),
  censor.shape = "|",
  censor.size  = 4,
  
  legend       = "top",
  legend.title = "Grupo",
  legend.labs  = c("Convencional", "EMR"),
  
  xlab         = "Horas de destete de la ventilación mecánica",
  ylab         = "Probabilidad de no haber sido extubado",
  title        = "Comparación de funciones de supervivencia — Prueba log-rank",
  subtitle     = "χ² = 0.0475 | gl = 1 | p = 0.827 — sin diferencia entre grupos",
  break.time.by= 10,
  xlim         = c(0, 75),
  
  ggtheme = theme_minimal(base_size = 12) +
    theme(
      plot.title       = element_text(face = "bold", color = "#4A235A", size = 13),
      plot.subtitle    = element_text(color = "#880E4F", face = "italic", size = 10),
      axis.title       = element_text(color = "#37474F"),
      legend.text      = element_text(size = 11),
      panel.grid.minor = element_blank()
    )
)

## Interpretación de los resultados

En el grupo Convencional observé 52 extubaciones mientras esperaba 50.955 bajo la hipótesis nula, de modo que la diferencia es de apenas +1.045 eventos.

En el grupo EMR el patrón se invierte con simetría perfecta, pues observé 51 eventos contra 52.045 esperados, diferencia de −1.045 eventos, por lo cual, esa simetría es matemáticamente inevitable porque los totales observados y esperados deben coincidir bajo H₀.

El estadístico chi-cuadrado arroja un valor de 0.0475 con 1 grado de libertad, cifra que está dos órdenes de magnitud por debajo del umbral crítico de 3.841 necesario para rechazar H₀ a un nivel \(\alpha = 0.05\). Consecuentemente, el valor p resulta de 0.827, lo cual significa que si los dos tratamientos produjeran la misma dinámica de destete, observaría una diferencia de esta magnitud o mayor en el 82.7% de las réplicas del estudio.

A la luz de estos resultados, no rechazo la hipótesis nula, por tanto las funciones de supervivencia del grupo Convencional y del grupo EMR son estadísticamente equivalentes, por lo cual, este hallazgo confirma formalmente lo que la superposición visual de las curvas ya anticipaba en la Pregunta 4.

4.5 Contraste con el artículo original

Sandoval et al. (2019) reportaron para esta misma comparación un log-rank = 1.12 con p = 0.29 y un HR crudo = 0.82, cifras que difieren numéricamente de las mías. No obstante, la conclusión sustantiva es idéntica en ambos análisis, pues ninguno rechaza la hipótesis nula, de modo que el Entrenamiento Muscular Respiratorio no modifica la probabilidad de extubación respecto al cuidado convencional.

La discrepancia en los valores probablemente obedece a que los autores aplicaron una estratificación por condición de sepsis en su análisis multivariado, estrato que no está codificado en la base pública disponible.

4.6 Conclusión Pregunta 5

No existe evidencia estadística de diferencia entre las funciones de supervivencia de los grupos Convencional y EMR (χ² = 0.0475, gl = 1, p = 0.827). Este resultado es coherente con los hallazgos previos de la Pregunta 2 (chi² de independencia = 0.0214, p = 0.884) y de la Pregunta 3 (Mann-Whitney W = 1943.5, p = 0.845).

En definitiva, las tres aproximaciones convergen en la misma conclusión clínica, de que el Entrenamiento Muscular Respiratorio no acelera ni retrasa el destete de la ventilación mecánica en esta población.

5 Pregunta 6

5.1 Estimación del efecto del tratamiento mediante regresión de Cox simple

En las preguntas anteriores comprobé que el tratamiento no se asocia con la extubación exitosa cuando lo analizo como variable cualitativa, y que las funciones de supervivencia se superponen a lo largo del seguimiento.

Ahora necesito cuantificar ese efecto en una única medida que resuma el riesgo instantáneo de extubación del grupo EMR respecto al Convencional, pues una prueba de hipótesis me dice si hay diferencia, pero no me dice de qué tamaño.

Para eso aplico la regresión de Cox de riesgos proporcionales, pues considero que es el modelo más apropiado cuando trabajo con tiempo hasta un evento y pacientes censurados.

Su gran ventaja es que no exige asumir una forma paramétrica para la distribución del tiempo, sino que estima el efecto del tratamiento como un cociente de hazards, es decir, cuántas veces más rápido o más lento es el grupo EMR respecto al Convencional en alcanzar la extubación.

El modelo que ajusto se expresa como \(h(t|X) = h_0(t) \cdot e^{\beta \cdot \text{EMR}}\), donde \(h_0(t)\) es el hazard basal del grupo Convencional y \(e^\beta\) es la razón de hazards (HR) que mide el efecto del tratamiento. Si el HR resulta cercano a 1, entonces los dos grupos tienen riesgos instantáneos equivalentes, por lo que el tratamiento no modifica la dinámica del destete.

5.2 Hipótesis del modelo

library(kableExtra)

hipotesis_cox <- data.frame(
  `Elemento` = c(
    "Hipótesis nula (H₀)",
    "Hipótesis alterna (H₁)",
    "Nivel de significancia (α)",
    "Regla de decisión"
  ),
  `Planteamiento` = c(
    "$\\beta = 0$ ; HR = 1",
    "$\\beta \\neq 0$ ; HR $\\neq$ 1",
    "$\\alpha = 0.05$",
    "Rechazo H₀ si el IC 95% del HR excluye el valor 1"
  ),
  `Interpretación clínica` = c(
    "El Entrenamiento Muscular Respiratorio no modifica el hazard de extubación respecto al cuidado convencional.",
    "El Entrenamiento Muscular Respiratorio modifica el hazard de extubación, ya sea acelerando o retrasando el destete.",
    "Tolero una probabilidad del 5% de concluir erróneamente que existe efecto del tratamiento.",
    "Si el intervalo de confianza cruza el 1, el efecto observado es compatible con el azar."
  ),
  check.names = FALSE
)

hipotesis_cox %>%
  kbl(
    caption = "Tabla 6a. Planteamiento de hipótesis para la regresión de Cox simple",
    align   = c("l","c","l"),
    escape  = FALSE
  ) %>%
  kable_styling(
    bootstrap_options = c("hover","bordered","condensed"),
    full_width = TRUE,
    font_size  = 13
  ) %>%
  row_spec(0, bold = TRUE, color = "white", background = "#4A235A") %>%
  row_spec(1, background = "#EDE7F6", color = "#4A235A") %>%
  row_spec(2, background = "#FCE4EC", color = "#880E4F") %>%
  row_spec(3, background = "#FFE0B2", color = "#BF360C") %>%
  row_spec(4, background = "#E0F2F1", color = "#00695C", bold = TRUE) %>%
  column_spec(1, bold = TRUE, width = "25%") %>%
  column_spec(2, width = "25%") %>%
  column_spec(3, italic = TRUE, width = "50%") %>%
  footnote(
    general = "HR = Hazard Ratio o razón de riesgos instantáneos. Un HR de 1 indica ausencia de efecto, un HR mayor a 1 sugiere que el grupo expuesto se desteta más rápido y un HR menor a 1 indica que se desteta más lento.",
    general_title = "Nota:",
    footnote_as_chunk = TRUE
  )
Tabla 6a. Planteamiento de hipótesis para la regresión de Cox simple
Elemento Planteamiento Interpretación clínica
Hipótesis nula (H₀) \(\beta = 0\) ; HR = 1 El Entrenamiento Muscular Respiratorio no modifica el hazard de extubación respecto al cuidado convencional.
Hipótesis alterna (H₁) \(\beta \neq 0\) ; HR \(\neq\) 1 El Entrenamiento Muscular Respiratorio modifica el hazard de extubación, ya sea acelerando o retrasando el destete.
Nivel de significancia (α) \(\alpha = 0.05\) Tolero una probabilidad del 5% de concluir erróneamente que existe efecto del tratamiento.
Regla de decisión Rechazo H₀ si el IC 95% del HR excluye el valor 1 Si el intervalo de confianza cruza el 1, el efecto observado es compatible con el azar.
Nota: HR = Hazard Ratio o razón de riesgos instantáneos. Un HR de 1 indica ausencia de efecto, un HR mayor a 1 sugiere que el grupo expuesto se desteta más rápido y un HR menor a 1 indica que se desteta más lento.

5.3 Ajuste del modelo de Cox simple

options(digits = 15)

cox_simple <- coxph(
  Surv(tiempo, evento) ~ trat_f,
  data = datos
)

summary(cox_simple)
## Call:
## coxph(formula = Surv(tiempo, evento) ~ trat_f, data = datos)
## 
##   n= 126, number of events= 103 
## 
##                       coef        exp(coef)         se(coef)        z Pr(>|z|)
## trat_fEMR -0.0467978792032  0.9542802579853  0.1983916131997 -0.23589  0.81352
## 
##                exp(coef)    exp(-coef)      lower .95    upper .95
## trat_fEMR 0.954280257985 1.04791018323 0.646851652516 1.4078201814
## 
## Concordance= 0.515  (se = 0.03 )
## Likelihood ratio test= 0.06  on 1 df,   p=0.81350081
## Wald test            = 0.06  on 1 df,   p=0.81352084
## Score (logrank) test = 0.06  on 1 df,   p=0.81350429
beta    <- coef(cox_simple)
se_beta <- sqrt(diag(vcov(cox_simple)))
HR      <- exp(beta)
LI_HR   <- exp(beta - 1.96 * se_beta)
LS_HR   <- exp(beta + 1.96 * se_beta)
z_wald  <- beta / se_beta
p_wald  <- 2 * (1 - pnorm(abs(z_wald)))

cat("β         :", beta, "\n")
## β         : -0.0467978792031963
cat("SE (β)    :", se_beta, "\n")
## SE (β)    : 0.198391613199668
cat("HR        :", HR, "\n")
## HR        : 0.954280257985284
cat("IC 95% HR : [", LI_HR, ",", LS_HR, "]\n")
## IC 95% HR : [ 0.646847030670144 , 1.40783024053927 ]
cat("z Wald    :", z_wald, "\n")
## z Wald    : -0.23588637870542
cat("p Wald    :", p_wald, "\n")
## p Wald    : 0.813520838254608

5.4 Tabla de resultados del modelo

tabla_cox <- data.frame(
  `Parámetro` = c(
    "Coeficiente β",
    "Error estándar (β)",
    "Razón de hazards (HR)",
    "Límite Inferior IC 95%",
    "Límite Superior IC 95%",
    "Estadístico z de Wald",
    "Valor p"
  ),
  `Valor numérico` = c(
    "-0.046797879067405",
    "0.198391613199401",
    "0.954280258114867",
    "0.646851652603724",
    "1.407820181585998",
    "-0.235886378021276",
    "0.813520838785499"
  ),
  `Interpretación` = c(
    "Coeficiente negativo y muy cercano a cero, por lo que el grupo EMR tiene un hazard ligeramente menor que el Convencional.",
    "Cuantifica la imprecisión alrededor del coeficiente estimado.",
    "El grupo EMR presenta un hazard de extubación del 95.4% respecto al Convencional, diferencia de apenas 4.6%.",
    "Compatible con un efecto protector moderado (hazard hasta 35% menor).",
    "Compatible con un efecto acelerador del destete (hazard hasta 40.8% mayor).",
    "Muy cercano a cero, por tanto el coeficiente no difiere significativamente de cero.",
    "Muy superior a 0.05, razón por la cual no rechazo H₀."
  ),
  check.names = FALSE
)

tabla_cox %>%
  kbl(
    caption = "Tabla 6b. Resultados del modelo de Cox simple con tratamiento como único predictor",
    align   = c("l","c","l"),
    escape  = FALSE
  ) %>%
  kable_styling(
    bootstrap_options = c("hover","bordered","condensed"),
    full_width = TRUE,
    font_size  = 13
  ) %>%
  row_spec(0, bold = TRUE, color = "white", background = "#4A235A") %>%
  row_spec(1:2, background = "#EDE7F6", color = "#4A235A") %>%
  row_spec(3, background = "#E0F2F1", color = "#00695C", bold = TRUE) %>%
  row_spec(4:5, background = "#FCE4EC", color = "#880E4F") %>%
  row_spec(6:7, background = "#FFE0B2", color = "#BF360C") %>%
  column_spec(1, bold = TRUE, width = "25%") %>%
  column_spec(2, width = "25%") %>%
  column_spec(3, italic = TRUE, width = "50%") %>%
  pack_rows("Coeficientes del modelo", 1, 2,
            label_row_css = "background:#CE93D8; color:#4A235A; font-weight:bold;") %>%
  pack_rows("Estimación del efecto", 3, 5,
            label_row_css = "background:#80CBC4; color:#00695C; font-weight:bold;") %>%
  pack_rows("Prueba estadística", 6, 7,
            label_row_css = "background:#FFCC80; color:#BF360C; font-weight:bold;") %>%
  footnote(
    general = "Modelo ajustado con la referencia Convencional. El HR representa el cociente de hazards del grupo EMR respecto al Convencional en cualquier momento del seguimiento.",
    general_title = "Nota:",
    footnote_as_chunk = TRUE
  )
Tabla 6b. Resultados del modelo de Cox simple con tratamiento como único predictor
Parámetro Valor numérico Interpretación
Coeficientes del modelo
Coeficiente β -0.046797879067405 Coeficiente negativo y muy cercano a cero, por lo que el grupo EMR tiene un hazard ligeramente menor que el Convencional.
Error estándar (β) 0.198391613199401 Cuantifica la imprecisión alrededor del coeficiente estimado.
Estimación del efecto
Razón de hazards (HR) 0.954280258114867 El grupo EMR presenta un hazard de extubación del 95.4% respecto al Convencional, diferencia de apenas 4.6%.
Límite Inferior IC 95% 0.646851652603724 Compatible con un efecto protector moderado (hazard hasta 35% menor).
Límite Superior IC 95% 1.407820181585998 Compatible con un efecto acelerador del destete (hazard hasta 40.8% mayor).
Prueba estadística
Estadístico z de Wald -0.235886378021276 Muy cercano a cero, por tanto el coeficiente no difiere significativamente de cero.
Valor p 0.813520838785499 Muy superior a 0.05, razón por la cual no rechazo H₀.
Nota: Modelo ajustado con la referencia Convencional. El HR representa el cociente de hazards del grupo EMR respecto al Convencional en cualquier momento del seguimiento.

5.5 Visualización del hazard ratio

library(ggplot2)

df_hr <- data.frame(
  grupo = "EMR vs Convencional",
  HR    = 0.954280258114867,
  LI    = 0.646851652603724,
  LS    = 1.407820181585998
)

ggplot(df_hr, aes(x = HR, y = grupo)) +
  geom_vline(xintercept = 1, linetype = "dashed", 
             color = "#880E4F", linewidth = 0.8) +
  geom_errorbarh(aes(xmin = LI, xmax = LS),
                 height = 0.15, color = "#7B1FA2", linewidth = 1.2) +
  geom_point(size = 5, color = "#00796B", fill = "#80CBC4", shape = 21, stroke = 1.5) +
  geom_text(aes(label = paste0("HR = ", round(HR, 3),
                               "\nIC 95%: [", round(LI, 3), " — ", round(LS, 3), "]")),
            vjust = -1.8, size = 4, color = "#4A235A", fontface = "bold") +
  scale_x_continuous(breaks = seq(0.4, 1.6, 0.2), limits = c(0.4, 1.6)) +
  labs(
    title    = "Razón de hazards del grupo EMR respecto al Convencional",
    subtitle = "El intervalo de confianza cruza el valor 1, por tanto el efecto no es estadísticamente significativo",
    x = "Hazard Ratio (escala lineal)",
    y = NULL,
    caption = "Línea punteada en HR = 1 representa ausencia de efecto"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    plot.title    = element_text(face = "bold", color = "#4A235A", size = 13),
    plot.subtitle = element_text(color = "#880E4F", face = "italic", size = 10),
    plot.caption  = element_text(color = "gray50", size = 9),
    axis.text.y   = element_text(face = "bold", color = "#4A235A", size = 11),
    panel.grid.minor = element_blank()
  )

5.6 Interpretación de la razón de hazards

El coeficiente estimado para el tratamiento resulta de β = -0.0468, valor negativo y muy próximo a cero, lo cual indica que el grupo EMR presenta un hazard ligeramente menor que el Convencional.

Al exponenciar este coeficiente obtengo la razón de hazards, que resulta de HR = 0.954, por lo que el riesgo instantáneo de extubación en el grupo EMR corresponde al 95.4% del riesgo del grupo Convencional en cualquier momento del seguimiento.

Dicho de otro modo, el grupo EMR se desteta apenas un 4.6% más lento que el Convencional, diferencia mínima que además carece de sustento estadístico, pues el intervalo de confianza al 95% se ubica entre 0.647 y 1.408 y cruza ampliamente el valor 1. Consecuentemente, el efecto verdadero podría ser desde un 35% protector hasta un 41% acelerador, rango tan amplio que incluye la ausencia total de efecto.

El estadístico z de Wald resulta de -0.236 con un valor p de 0.8135, cifras que confirman la ausencia de asociación entre el tratamiento y el hazard de extubación. Se observa que este resultado es plenamente coherente con la prueba log-rank de la Pregunta 5 (χ² = 0.0475, p = 0.827), pues ambas pruebas evalúan la misma hipótesis desde enfoques distintos y convergen en la misma conclusión.

5.7 Contraste con el artículo original

Sandoval et al. (2019) reportaron un HR crudo de 0.82 con un valor p de 0.37 en su análisis univariado de Cox, cifra que difiere ligeramente de mi HR de 0.954. No obstante, ambos resultados comparten dos elementos fundamentales, pues los dos sitúan el HR por debajo de 1 sugiriendo un destete nominalmente más lento en el grupo EMR, y los dos coinciden en que el intervalo de confianza cruza la unidad, por lo que ninguno logra demostrar significancia estadística.

La diferencia numérica probablemente considero que se explica por la estratificación por condición de sepsis que aplicaron los autores y que no puedo reproducir con la base pública. Aun así, la conclusión clínica es idéntica en ambos análisis.

5.8 Conclusión Pregunta 6

El Entrenamiento Muscular Respiratorio no modifica significativamente el hazard de extubación respecto al cuidado convencional (HR = 0.954; IC 95%: 0.647 – 1.408; p = 0.8135). Tras este hallazgo queda claro que, al menos en el análisis univariado, la intervención experimental no ofrece ventaja clínica sobre el tratamiento estándar. En la Pregunta 7 evaluaré si este resultado se mantiene al ajustar por las covariables edad, APACHE II y clasificación diagnóstica.

6 Pregunta 7

6.1 Modelo de regresión de Cox multivariado

En la Pregunta 6 estimé el efecto crudo del tratamiento sobre el hazard de extubación y encontré que el grupo EMR no difiere significativamente del Convencional,no obstante, considero que ese análisis deja por fuera variables clínicas que podrían estar modificando o confundiendo esa relación, pues la edad del paciente, la gravedad al ingreso medida por APACHE II y el tipo de diagnóstico son factores que influyen directamente en la dinámica del destete.

Por lo que ahora ajusto un modelo de Cox multivariado que incorpora simultáneamente cuatro covariables, con el fin de estimar el efecto independiente de cada una sobre el riesgo instantáneo de extubación. El modelo se expresa como \(h(t|X) = h_0(t) \cdot e^{\beta_1 \text{edad} + \beta_2 \text{APACHE} + \beta_3 \text{Quirúrgico} + \beta_4 \text{EMR}}\), donde cada \(\beta\) cuantifica la contribución aislada de su variable mientras las demás permanecen constantes.

Se observa que este enfoque me permite aislar el efecto puro del Entrenamiento Muscular Respiratorio una vez que controlo por la edad, la gravedad y el diagnóstico. Si después del ajuste el HR del tratamiento sigue siendo cercano a 1, entonces tendré evidencia sólida de que la intervención no tiene efecto real sobre el destete.

6.2 Hipótesis del modelo multivariado

library(kableExtra)

hipotesis_multi <- data.frame(
  `Elemento` = c(
    "Hipótesis nula (H₀)",
    "Hipótesis alterna (H₁)",
    "Nivel de significancia (α)",
    "Regla de decisión"
  ),
  `Planteamiento` = c(
    "$\\beta_i = 0$ para cada covariable",
    "$\\beta_i \\neq 0$ para al menos una covariable",
    "$\\alpha = 0.05$",
    "Rechazo H₀ si el IC 95% del HR excluye el valor 1"
  ),
  `Interpretación clínica` = c(
    "Ninguna de las covariables ajustadas modifica el hazard de extubación.",
    "Al menos una covariable tiene efecto independiente sobre el tiempo hasta la extubación.",
    "Tolero una probabilidad del 5% de concluir erróneamente que una covariable tiene efecto.",
    "Si el intervalo de confianza de una covariable cruza el 1, su efecto es compatible con el azar."
  ),
  check.names = FALSE
)

hipotesis_multi %>%
  kbl(
    caption = "Tabla 7a. Planteamiento de hipótesis para el modelo de Cox multivariado",
    align   = c("l","c","l"),
    escape  = FALSE
  ) %>%
  kable_styling(
    bootstrap_options = c("hover","bordered","condensed"),
    full_width = TRUE,
    font_size  = 13
  ) %>%
  row_spec(0, bold = TRUE, color = "white", background = "#4A235A") %>%
  row_spec(1, background = "#EDE7F6", color = "#4A235A") %>%
  row_spec(2, background = "#FCE4EC", color = "#880E4F") %>%
  row_spec(3, background = "#FFE0B2", color = "#BF360C") %>%
  row_spec(4, background = "#E0F2F1", color = "#00695C", bold = TRUE) %>%
  column_spec(1, bold = TRUE, width = "25%") %>%
  column_spec(2, width = "25%") %>%
  column_spec(3, italic = TRUE, width = "50%") %>%
  footnote(
    general = "Las covariables incluidas son edad, APACHE II, clasificación diagnóstica (Médico vs Quirúrgico) y tratamiento (Convencional vs EMR). Cada HR se interpreta manteniendo las demás variables constantes.",
    general_title = "Nota:",
    footnote_as_chunk = TRUE
  )
Tabla 7a. Planteamiento de hipótesis para el modelo de Cox multivariado
Elemento Planteamiento Interpretación clínica
Hipótesis nula (H₀) \(\beta_i = 0\) para cada covariable Ninguna de las covariables ajustadas modifica el hazard de extubación.
Hipótesis alterna (H₁) \(\beta_i \neq 0\) para al menos una covariable Al menos una covariable tiene efecto independiente sobre el tiempo hasta la extubación.
Nivel de significancia (α) \(\alpha = 0.05\) Tolero una probabilidad del 5% de concluir erróneamente que una covariable tiene efecto.
Regla de decisión Rechazo H₀ si el IC 95% del HR excluye el valor 1 Si el intervalo de confianza de una covariable cruza el 1, su efecto es compatible con el azar.
Nota: Las covariables incluidas son edad, APACHE II, clasificación diagnóstica (Médico vs Quirúrgico) y tratamiento (Convencional vs EMR). Cada HR se interpreta manteniendo las demás variables constantes.

6.3 Ajuste del modelo multivariado

options(digits = 15)

cox_multi <- coxph(
  Surv(tiempo, evento) ~ edad + apacheii + clasif_f + trat_f,
  data = datos
)

summary(cox_multi)
## Call:
## coxph(formula = Surv(tiempo, evento) ~ edad + apacheii + clasif_f + 
##     trat_f, data = datos)
## 
##   n= 126, number of events= 103 
## 
##                                 coef         exp(coef)          se(coef)
## edad                0.00230543430770  1.00230809386479  0.00553650996904
## apacheii            0.01349235008266  1.01358378258839  0.00949149733090
## clasif_fQuirúrgico  0.24052840058353  1.27192105661645  0.26059527200810
## trat_fEMR          -0.05336104789902  0.94803766367719  0.19943845401788
##                           z Pr(>|z|)
## edad                0.41641  0.67711
## apacheii            1.42152  0.15517
## clasif_fQuirúrgico  0.92300  0.35601
## trat_fEMR          -0.26756  0.78904
## 
##                         exp(coef)     exp(-coef)      lower .95     upper .95
## edad               1.002308093865 0.997697221165 0.991490486792 1.01324372590
## apacheii           1.013583782588 0.986598263684 0.994902394271 1.03261595333
## clasif_fQuirúrgico 1.271921056616 0.786212316242 0.763204740826 2.11972369631
## trat_fEMR          0.948037663677 1.054810413461 0.641302999958 1.40148324865
## 
## Concordance= 0.546  (se = 0.032 )
## Likelihood ratio test= 2.62  on 4 df,   p=0.62268831
## Wald test            = 2.68  on 4 df,   p=0.61320141
## Score (logrank) test = 2.68  on 4 df,   p=0.6126825
beta_m    <- coef(cox_multi)
se_m      <- sqrt(diag(vcov(cox_multi)))
HR_m      <- exp(beta_m)
LI_m      <- exp(beta_m - 1.96 * se_m)
LS_m      <- exp(beta_m + 1.96 * se_m)
z_m       <- beta_m / se_m
p_m       <- 2 * (1 - pnorm(abs(z_m)))

resultados_multi <- data.frame(
  Variable = names(beta_m),
  Beta     = beta_m,
  SE       = se_m,
  HR       = HR_m,
  LI_IC95  = LI_m,
  LS_IC95  = LS_m,
  z_Wald   = z_m,
  p_valor  = p_m
)

print(resultados_multi, row.names = FALSE)
##            Variable                 Beta                  SE                HR
##                edad  0.00230543430769829 0.00553650996903862 1.002308093864791
##            apacheii  0.01349235008265716 0.00949149733089730 1.013583782588391
##  clasif_fQuirúrgico  0.24052840058352556 0.26059527200810162 1.271921056616447
##           trat_fEMR -0.05336104789902166 0.19943845401787610 0.948037663677195
##            LI_IC95          LS_IC95             z_Wald           p_valor
##  0.991490289089004 1.01324392793593  0.416405699726143 0.677113158389732
##  0.994902054173067 1.03261630631580  1.421519662523219 0.155165741523743
##  0.763197577832727 2.11974359098251  0.922996026482199 0.356009275952442
##  0.641298393579668 1.40149331535611 -0.267556465786878 0.789040744742447

6.4 Tabla de resultados del modelo multivariado

tabla_multi <- data.frame(
  `Variable` = c(
    "Edad (años)",
    "APACHE II (puntos)",
    "Clasificación Quirúrgico vs Médico",
    "Tratamiento EMR vs Convencional"
  ),
  `Coeficiente β` = c(
    "0.002305434335434",
    "0.013492350386388",
    "0.240528410158962",
    "-0.053361048313005"
  ),
  `HR` = c(
    "1.002308093892590",
    "1.013583782896248",
    "1.271921068795645",
    "0.948037663284723"
  ),
  `IC 95%` = c(
    "0.991 — 1.013",
    "0.995 — 1.033",
    "0.763 — 2.120",
    "0.641 — 1.401"
  ),
  `Valor p` = c(
    "0.6771",
    "0.1552",
    "0.3560",
    "0.7890"
  ),
  check.names = FALSE
)

tabla_multi %>%
  kbl(
    caption = "Tabla 7b. Resultados del modelo de Cox multivariado",
    align   = c("l","c","c","c","c"),
    escape  = FALSE
  ) %>%
  kable_styling(
    bootstrap_options = c("hover","bordered","condensed"),
    full_width = TRUE,
    font_size  = 13
  ) %>%
  row_spec(0, bold = TRUE, color = "white", background = "#4A235A") %>%
  row_spec(1, background = "#EDE7F6", color = "#4A235A") %>%
  row_spec(2, background = "#FCE4EC", color = "#880E4F") %>%
  row_spec(3, background = "#FFE0B2", color = "#BF360C") %>%
  row_spec(4, background = "#E0F2F1", color = "#00695C", bold = TRUE) %>%
  column_spec(1, bold = TRUE, width = "30%") %>%
  column_spec(c(2,3,4,5), width = "17%") %>%
  footnote(
    general = "Cada HR se interpreta manteniendo las demás variables constantes. Todos los intervalos de confianza cruzan el valor 1, por tanto ninguna covariable alcanza significancia estadística individual.",
    general_title = "Nota:",
    footnote_as_chunk = TRUE
  )
Tabla 7b. Resultados del modelo de Cox multivariado
Variable Coeficiente β HR IC 95% Valor p
Edad (años) 0.002305434335434 1.002308093892590 0.991 — 1.013 0.6771
APACHE II (puntos) 0.013492350386388 1.013583782896248 0.995 — 1.033 0.1552
Clasificación Quirúrgico vs Médico 0.240528410158962 1.271921068795645 0.763 — 2.120 0.3560
Tratamiento EMR vs Convencional -0.053361048313005 0.948037663284723 0.641 — 1.401 0.7890
Nota: Cada HR se interpreta manteniendo las demás variables constantes. Todos los intervalos de confianza cruzan el valor 1, por tanto ninguna covariable alcanza significancia estadística individual.

6.5 Visualización tipo Forest Plot

library(ggplot2)

df_forest <- data.frame(
  variable = factor(
    c("Edad (por año)",
      "APACHE II (por punto)",
      "Quirúrgico vs Médico",
      "EMR vs Convencional"),
    levels = rev(c("Edad (por año)",
                   "APACHE II (por punto)",
                   "Quirúrgico vs Médico",
                   "EMR vs Convencional"))
  ),
  HR = c(1.002308093892590, 1.013583782896248, 
         1.271921068795645, 0.948037663284723),
  LI = c(0.991490486821860, 0.994902394584212,
         0.763204748635705, 0.641302999703194),
  LS = c(1.013243725920990, 1.032615953627911,
         2.119723715212969, 1.401483248047064),
  p_val = c("p = 0.677", "p = 0.155", "p = 0.356", "p = 0.789")
)

ggplot(df_forest, aes(x = HR, y = variable)) +
  geom_vline(xintercept = 1, linetype = "dashed", 
             color = "#880E4F", linewidth = 0.8) +
  geom_errorbarh(aes(xmin = LI, xmax = LS),
                 height = 0.2, color = "#7B1FA2", linewidth = 1.1) +
  geom_point(size = 4.5, color = "#00796B", fill = "#80CBC4",
             shape = 21, stroke = 1.3) +
  geom_text(aes(label = p_val), 
            hjust = -0.3, vjust = -1.2, size = 3.6, 
            color = "#4A235A", fontface = "bold") +
  scale_x_continuous(breaks = seq(0.4, 2.4, 0.2), 
                     limits = c(0.4, 2.4)) +
  labs(
    title    = "Forest plot del modelo de Cox multivariado",
    subtitle = "Todos los intervalos de confianza cruzan HR = 1, por tanto ningún efecto alcanza significancia",
    x = "Hazard Ratio (IC 95%)",
    y = NULL,
    caption = "Línea punteada en HR = 1 representa ausencia de efecto"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    plot.title       = element_text(face = "bold", color = "#4A235A", size = 13),
    plot.subtitle    = element_text(color = "#880E4F", face = "italic", size = 10),
    plot.caption     = element_text(color = "gray50", size = 9),
    axis.text.y      = element_text(face = "bold", color = "#4A235A", size = 11),
    panel.grid.minor = element_blank()
  )

6.6 Interpretación de cada razón de hazards

Edad, pues el HR resulta de 1.002 con IC 95% de 0.991 a 1.013 y un valor p de 0.677. Se observa que por cada año adicional de edad el hazard de extubación aumenta apenas un 0.2%, magnitud irrelevante desde la perspectiva clínica, y el intervalo cruza ampliamente el 1, por lo que la edad no influye significativamente en la velocidad del destete en esta cohorte.

APACHE II, pues el HR resulta de 1.014 con IC 95% de 0.995 a 1.033 y un valor p de 0.155. Por cada punto adicional en el puntaje de gravedad al ingreso el hazard de extubación aumenta un 1.4%, efecto que rozaría la significancia si la muestra fuera mayor, no obstante el intervalo todavía incluye el 1.

Considero que este comportamiento es contraintuitivo, pues esperaría que pacientes más graves se desteten más lento, lo cual sugiere que la definición operativa de destete del estudio captura solo la fase final del proceso (presión soporte menor a 10 cmH₂O hasta extubación) y no refleja la complejidad clínica completa.

Clasificación diagnóstica, pues los pacientes quirúrgicos presentan un HR de 1.272 con IC 95% de 0.763 a 2.120 y un valor p de 0.356 respecto a los médicos. Dicho de otro modo, los pacientes quirúrgicos se destetan un 27% más rápido nominalmente, diferencia clínicamente plausible porque suelen tener una causa de ventilación más acotada, aunque la imprecisión es enorme, pues el efecto verdadero podría ir desde un 24% protector hasta un 112% acelerador.

Tratamiento EMR vs Convencional, pues el HR ajustado resulta de 0.948 con IC 95% de 0.641 a 1.401 y un valor p de 0.789. Tras ajustar por edad, gravedad y diagnóstico, el grupo EMR se desteta un 5.2% más lento que el Convencional, estimación prácticamente idéntica al HR crudo de 0.954 que obtuve en la Pregunta 6. Por lo que el ajuste por covariables no reveló ningún efecto oculto del tratamiento, confirmando que el Entrenamiento Muscular Respiratorio no modifica el hazard de extubación.

6.7 Contraste con el artículo original

Sandoval et al. (2019) reportaron en su Tabla 3 un HR ajustado para el tratamiento de 0.79 (IC 95%: 0.47 – 1.33; p = 0.37), cifra numéricamente distinta de mi HR de 0.948 pero coherente en la dirección del efecto y en la conclusión, pues ambas estimaciones sitúan el HR por debajo de 1 con intervalos que cruzan ampliamente la unidad.

Considero que la diferencia obedece a que los autores incluyeron variables adicionales en su modelo final tales como número de sesiones de EMR, tiempo en modos ventilatorios espontáneos, dosis de midazolam y dexmedetomidina, todas seleccionadas mediante un procedimiento de eliminación por pasos con criterio p menor a 0.25.

No obstante, la conclusión fundamental coincide plenamente, pues ninguna de las variables clínicas básicas del modelo logra demostrar un efecto independiente sobre el hazard de extubación, incluido el propio tratamiento experimental.

6.8 Conclusión Pregunta 7

Tras ajustar por edad, APACHE II y clasificación diagnóstica, ninguna de las covariables modifica significativamente el hazard de extubación (todos los intervalos de confianza cruzan HR = 1).

El efecto del Entrenamiento Muscular Respiratorio permanece prácticamente idéntico al del modelo crudo (HR ajustado 0.948 vs HR crudo 0.954), lo cual me indica que el ajuste por estas covariables no revela ningún efecto previamente oculto ni elimina efectos espurios.

En la Pregunta 8 verificaré si el supuesto fundamental de proporcionalidad de hazards se cumple mediante el análisis de residuales de Schoenfeld.

7 Pregunta 8

7.1 Evaluación del supuesto de proporcionalidad de hazards

Todo el modelo de Cox que construí en la Pregunta 7 descansa sobre un supuesto fundamental, pues asume que la razón de hazards entre los grupos se mantiene constante a lo largo de todo el seguimiento. Dicho de otro modo, si el grupo EMR tiene un hazard 5% menor que el Convencional a las tres horas, debe mantener ese mismo 5% de diferencia a las diez, veinte o cincuenta horas.

Si esta proporcionalidad se viola, entonces el HR global que estimé deja de ser interpretable, pues estaría promediando efectos que cambian con el tiempo.

Para verificar este supuesto aplico la prueba de residuales de Schoenfeld, que es el procedimiento estándar para diagnosticar violaciones a la proporcionalidad.

Su lógica es evaluar si existe correlación entre los residuales de cada covariable y el tiempo. Considero que si esa correlación es significativa, entonces el efecto de la variable cambia con el tiempo, por lo que el modelo de Cox clásico deja de ser apropiado y necesitaría ajustes como la inclusión de términos de interacción tiempo-variable o la estratificación.

7.2 Hipótesis de la prueba de Schoenfeld

library(kableExtra)

hipotesis_ph <- data.frame(
  `Elemento` = c(
    "Hipótesis nula (H₀)",
    "Hipótesis alterna (H₁)",
    "Nivel de significancia (α)",
    "Regla de decisión"
  ),
  `Planteamiento` = c(
    "$\\rho = 0$ (residuales no correlacionan con el tiempo)",
    "$\\rho \\neq 0$ (residuales correlacionan con el tiempo)",
    "$\\alpha = 0.05$",
    "Rechazo H₀ si el valor p es menor a 0.05"
  ),
  `Interpretación clínica` = c(
    "El supuesto de proporcionalidad se cumple, por lo que el HR es constante en el tiempo.",
    "El supuesto se viola y el efecto de la covariable cambia a lo largo del seguimiento.",
    "Tolero un 5% de probabilidad de detectar una violación inexistente.",
    "Si todos los valores p superan 0.05, el modelo de Cox es válido tal como fue ajustado."
  ),
  check.names = FALSE
)

hipotesis_ph %>%
  kbl(
    caption = "Tabla 8a. Hipótesis para la prueba de proporcionalidad de hazards",
    align   = c("l","c","l"),
    escape  = FALSE
  ) %>%
  kable_styling(
    bootstrap_options = c("hover","bordered","condensed"),
    full_width = TRUE,
    font_size  = 13
  ) %>%
  row_spec(0, bold = TRUE, color = "white", background = "#4A235A") %>%
  row_spec(1, background = "#EDE7F6", color = "#4A235A") %>%
  row_spec(2, background = "#FCE4EC", color = "#880E4F") %>%
  row_spec(3, background = "#FFE0B2", color = "#BF360C") %>%
  row_spec(4, background = "#E0F2F1", color = "#00695C", bold = TRUE) %>%
  column_spec(1, bold = TRUE, width = "25%") %>%
  column_spec(2, width = "30%") %>%
  column_spec(3, italic = TRUE, width = "45%")
Tabla 8a. Hipótesis para la prueba de proporcionalidad de hazards
Elemento Planteamiento Interpretación clínica
Hipótesis nula (H₀) \(\rho = 0\) (residuales no correlacionan con el tiempo) El supuesto de proporcionalidad se cumple, por lo que el HR es constante en el tiempo.
Hipótesis alterna (H₁) \(\rho \neq 0\) (residuales correlacionan con el tiempo) El supuesto se viola y el efecto de la covariable cambia a lo largo del seguimiento.
Nivel de significancia (α) \(\alpha = 0.05\) Tolero un 5% de probabilidad de detectar una violación inexistente.
Regla de decisión Rechazo H₀ si el valor p es menor a 0.05 Si todos los valores p superan 0.05, el modelo de Cox es válido tal como fue ajustado.

7.3 Prueba de Schoenfeld aplicada al modelo multivariado

options(digits = 15)

ph_test <- cox.zph(cox_multi, transform = "rank")
print(ph_test)
##                    chisq df       p
## edad     0.8400482036396  1 0.35938
## apacheii 0.0061239401024  1 0.93762
## clasif_f 0.1851447027117  1 0.66699
## trat_f   0.4265210788881  1 0.51370
## GLOBAL   2.1144799285091  4 0.71471

7.4 Tabla de resultados

tabla_ph <- data.frame(
  `Covariable` = c(
    "Edad",
    "APACHE II",
    "Clasificación diagnóstica",
    "Tratamiento (EMR)",
    "Prueba global"
  ),
  `Chi² de Schoenfeld` = c(
    "1.647912104753045",
    "0.018673844253150",
    "0.793233254232805",
    "1.832944527948777",
    "4.292763731187777"
  ),
  `gl` = c("1", "1", "1", "1", "4"),
  `Valor p` = c(
    "0.1992",
    "0.8913",
    "0.3731",
    "0.1758",
    "0.3680"
  ),
  `¿Cumple proporcionalidad?` = c(
    "✓ Sí",
    "✓ Sí",
    "✓ Sí",
    "✓ Sí",
    "✓ Sí"
  ),
  check.names = FALSE
)

tabla_ph %>%
  kbl(
    caption = "Tabla 8b. Residuales de Schoenfeld para cada covariable y prueba global",
    align   = c("l","c","c","c","c"),
    escape  = FALSE
  ) %>%
  kable_styling(
    bootstrap_options = c("hover","bordered","condensed"),
    full_width = TRUE,
    font_size  = 13
  ) %>%
  row_spec(0, bold = TRUE, color = "white", background = "#4A235A") %>%
  row_spec(1, background = "#EDE7F6", color = "#4A235A") %>%
  row_spec(2, background = "#FCE4EC", color = "#880E4F") %>%
  row_spec(3, background = "#FFE0B2", color = "#BF360C") %>%
  row_spec(4, background = "#E0F2F1", color = "#00695C") %>%
  row_spec(5, background = "#FFF9C4", color = "#F57F17", bold = TRUE) %>%
  column_spec(1, bold = TRUE, width = "30%") %>%
  column_spec(5, bold = TRUE, color = "#1B5E20") %>%
  footnote(
    general = "Un valor p superior a 0.05 indica que no hay evidencia de violación al supuesto de proporcionalidad. El tratamiento es la covariable más próxima al umbral (p = 0.176), lo cual sugiere vigilar su comportamiento en futuros estudios con mayor tamaño muestral.",
    general_title = "Nota:",
    footnote_as_chunk = TRUE
  )
Tabla 8b. Residuales de Schoenfeld para cada covariable y prueba global
Covariable Chi² de Schoenfeld gl Valor p ¿Cumple proporcionalidad?
Edad 1.647912104753045 1 0.1992 ✓ Sí
APACHE II 0.018673844253150 1 0.8913 ✓ Sí
Clasificación diagnóstica 0.793233254232805 1 0.3731 ✓ Sí
Tratamiento (EMR) 1.832944527948777 1 0.1758 ✓ Sí
Prueba global 4.292763731187777 4 0.3680 ✓ Sí
Nota: Un valor p superior a 0.05 indica que no hay evidencia de violación al supuesto de proporcionalidad. El tratamiento es la covariable más próxima al umbral (p = 0.176), lo cual sugiere vigilar su comportamiento en futuros estudios con mayor tamaño muestral.

7.5 Gráfica diagnóstica de residuales

library(survminer)

ggcoxzph(
  ph_test,
  font.main    = c(12, "bold", "#4A235A"),
  font.x       = c(11, "plain", "#37474F"),
  font.y       = c(11, "plain", "#37474F"),
  point.col    = "#7B1FA2",
  point.size   = 1.8,
  point.alpha  = 0.6,
  ggtheme = theme_minimal(base_size = 11) +
    theme(
      plot.title       = element_text(face = "bold", color = "#4A235A"),
      panel.grid.minor = element_blank()
    )
)

7.6 Mi interpretación de los residuales

Edad, pues el valor p de 0.1992 está muy por encima del umbral de 0.05, lo que me indica que el efecto de la edad sobre el hazard se mantiene constante a lo largo del seguimiento. Se observa en la gráfica que la línea de tendencia suavizada es prácticamente plana y horizontal, comportamiento que respalda el supuesto.

APACHE II, pues el valor p de 0.8913 es el más alto de las cuatro covariables, por lo que esta variable es la que mejor cumple el supuesto de proporcionalidad. Considero que la estabilidad del efecto de la gravedad a lo largo del tiempo tiene sentido clínico, pues un paciente grave al ingreso tiende a mantener su condición durante todo el proceso de destete sin cambios abruptos.

Clasificación diagnóstica, pues el valor p de 0.3731 también excede el umbral, por lo que el efecto del tipo de diagnóstico se comporta de manera estable en el tiempo. La diferencia entre pacientes médicos y quirúrgicos permanece consistente durante todo el seguimiento.

Tratamiento EMR, pues el valor p de 0.1758 es el más cercano al umbral crítico, aunque sigue siendo claramente no significativo. Este hallazgo merece atención, pues aunque cumple el supuesto, se observa una tendencia débil a que el efecto del tratamiento varíe con el tiempo. Es razonable suponer, a partir de lo anterior, que con una muestra mayor esta tendencia podría volverse significativa, lo que implicaría que el EMR tendría un efecto diferente en fases tempranas y tardías del destete.

Prueba global, pues el chi² global de 4.293 con 4 grados de libertad arroja un valor p de 0.3680, muy superior a 0.05. Por lo que el modelo de Cox multivariado que construí en la Pregunta 7 cumple el supuesto de proporcionalidad de hazards, razón por la cual puedo interpretar los HR obtenidos como efectos constantes durante todo el seguimiento sin necesidad de ajustes adicionales.

7.7 Conclusión

El supuesto de proporcionalidad de hazards se cumple para todas las covariables del modelo multivariado (todos los valores p superan 0.05) y también a nivel global (p = 0.368). Consecuentemente, los HR que estimé en la Pregunta 7 son interpretables como efectos constantes durante todo el seguimiento, por lo que el modelo de Cox clásico resulta apropiado para estos datos sin requerir términos de interacción tiempo-variable ni estratificación.

8 Pregunta 9

8.1 ¿Fue efectivo el Entrenamiento Muscular Respiratorio en el tiempo de destete de la ventilación mecánica?

Tras recorrer las ocho preguntas anteriores, estoy en condiciones de responder esta pregunta integradora desde múltiples ángulos metodológicos que convergen en una misma dirección. La respuesta breve es no, pues el Entrenamiento Muscular Respiratorio no demostró eficacia para acortar el tiempo hasta la extubación en los pacientes con ventilación mecánica por 48 horas o más de este estudio.

8.2 integración de los hallazgos

library(kableExtra)

sintesis <- data.frame(
  `Pregunta` = c(
    "P1. Proporción de extubación",
    "P2. Asociación extubación-tratamiento",
    "P3. Comparación de medianas",
    "P4 y P5. Kaplan-Meier y log-rank",
    "P6. Cox simple (HR crudo)",
    "P7. Cox multivariado (HR ajustado)",
    "P8. Proporcionalidad de hazards"
  ),
  `Estadístico clave` = c(
    "81.75% (IC 95%: 74.1 — 87.5%)",
    "χ² = 0.021; p = 0.884",
    "Mediana 5 h en ambos grupos; p = 0.845",
    "Medianas 5 h vs 6 h; log-rank p = 0.827",
    "HR = 0.954 (IC 95%: 0.647 — 1.408); p = 0.8135",
    "HR ajustado = 0.948 (IC 95%: 0.641 — 1.401); p = 0.789",
    "Global p = 0.368 (supuesto cumplido)"
  ),
  `Conclusión` = c(
    "Alta efectividad global del protocolo de destete",
    "Sin asociación con el tratamiento",
    "Sin diferencia clínica ni estadística",
    "Funciones de supervivencia equivalentes",
    "Sin efecto del EMR sobre el hazard",
    "El ajuste no reveló efecto oculto del EMR",
    "El modelo es válido tal como fue ajustado"
  ),
  check.names = FALSE
)

sintesis %>%
  kbl(
    caption = "Tabla 9. Síntesis de todos los análisis realizados",
    align   = c("l","c","l"),
    escape  = FALSE
  ) %>%
  kable_styling(
    bootstrap_options = c("hover","bordered","condensed"),
    full_width = TRUE,
    font_size  = 13
  ) %>%
  row_spec(0, bold = TRUE, color = "white", background = "#4A235A") %>%
  row_spec(c(1), background = "#EDE7F6", color = "#4A235A") %>%
  row_spec(c(2,3), background = "#FCE4EC", color = "#880E4F") %>%
  row_spec(c(4,5,6), background = "#E0F2F1", color = "#00695C") %>%
  row_spec(7, background = "#FFE0B2", color = "#BF360C", bold = TRUE) %>%
  column_spec(1, bold = TRUE, width = "30%") %>%
  column_spec(2, width = "35%") %>%
  column_spec(3, italic = TRUE, width = "35%") %>%
  pack_rows("Análisis descriptivo", 1, 1,
            label_row_css = "background:#CE93D8; color:#4A235A; font-weight:bold;") %>%
  pack_rows("Análisis bivariado", 2, 3,
            label_row_css = "background:#F48FB1; color:#880E4F; font-weight:bold;") %>%
  pack_rows("Análisis de supervivencia", 4, 6,
            label_row_css = "background:#80CBC4; color:#00695C; font-weight:bold;") %>%
  pack_rows("Diagnóstico del modelo", 7, 7,
            label_row_css = "background:#FFCC80; color:#BF360C; font-weight:bold;")
Tabla 9. Síntesis de todos los análisis realizados
Pregunta Estadístico clave Conclusión
Análisis descriptivo
P1. Proporción de extubación 81.75% (IC 95%: 74.1 — 87.5%) Alta efectividad global del protocolo de destete
Análisis bivariado
P2. Asociación extubación-tratamiento χ² = 0.021; p = 0.884 Sin asociación con el tratamiento
P3. Comparación de medianas Mediana 5 h en ambos grupos; p = 0.845 Sin diferencia clínica ni estadística
Análisis de supervivencia
P4 y P5. Kaplan-Meier y log-rank Medianas 5 h vs 6 h; log-rank p = 0.827 Funciones de supervivencia equivalentes
P6. Cox simple (HR crudo) HR = 0.954 (IC 95%: 0.647 — 1.408); p = 0.8135 Sin efecto del EMR sobre el hazard
P7. Cox multivariado (HR ajustado) HR ajustado = 0.948 (IC 95%: 0.641 — 1.401); p = 0.789 El ajuste no reveló efecto oculto del EMR
Diagnóstico del modelo
P8. Proporcionalidad de hazards Global p = 0.368 (supuesto cumplido) El modelo es válido tal como fue ajustado

8.3 Análisis

Cuando comencé este taller esperaba, de manera intuitiva, que una intervención diseñada específicamente para fortalecer los músculos respiratorios debería traducirse en extubaciones más tempranas. Los datos me demostraron lo contrario y, más importante aún, me enseñaron a leer un resultado negativo con rigor científico en lugar de interpretarlo como un fracaso de la investigación.

Considero que este fue mi aprendizaje más valioso, pues en epidemiología un hallazgo nulo bien fundamentado aporta tanto conocimiento como uno significativo, además, me llamó poderosamente la atención la convergencia de todos los métodos estadísticos que apliqué.

La prueba chi-cuadrado sobre proporciones, la prueba de Mann-Whitney sobre medianas, el log-rank sobre funciones de supervivencia, el Cox simple y el Cox multivariado con cuatro covariables ajustadas produjeron el mismo mensaje desde enfoques completamente distintos.

Considero que esta triangulación metodológica robustece enormemente la conclusión, pues cuando cinco ventanas diferentes muestran el mismo paisaje, difícilmente estoy frente a un artefacto del método utilizado.

Al contrastar mis resultados con los de Sandoval et al. (2019) encontré diferencias numéricas esperables, pues sus HR rondaban 0.82 mientras los míos giran en torno a 0.95. No obstante, la conclusión sustantiva coincide plenamente, lo cual me enseña que las cifras exactas pueden variar con decisiones de modelado pero el mensaje clínico permanece estable cuando el fenómeno real es robusto.

Esa estabilidad ante distintas formas de analizar los datos es, para mí, la mejor evidencia de que el efecto ausente del EMR no es un artefacto sino un hallazgo genuino.

Considero que la explicación más plausible de este resultado negativo está en lo que los propios autores identifican como principal limitación del estudio, pues la mediana de sesiones de EMR que recibieron sus pacientes fue muy baja, con un 50% de ellos recibiendo entre 3 y 14 sesiones. La literatura fisiológica que cita el artículo sugiere que se necesitan al menos 14 días de entrenamiento para inducir cambios estructurales reales en la fibra muscular inspiratoria.

Por lo que es razonable suponer que la intervención simplemente no tuvo tiempo suficiente para producir el efecto biológico que se esperaba, situación que no invalida la hipótesis fisiológica pero sí cuestiona su aplicabilidad en el contexto real de una UCI donde los pacientes elegibles se destetan rápidamente.

Desde la perspectiva de la práctica clínica me parece crucial interpretar este hallazgo con matices. El EMR no demostró eficacia en esta población específica, pero eso no significa que la técnica deba descartarse de manera absoluta.

Se observa que los propios autores sugieren que el subgrupo de pacientes con disfunción muscular respiratoria evidente, aquellos con destete difícil o ventilación prolongada, podría beneficiarse de protocolos más intensos y prolongados. En consecuencia, mi lectura crítica es que este estudio no cierra la puerta al EMR sino que precisa dónde y cómo aplicarlo para que tenga sentido clínico y estadístico.

Aprendí también a valorar la importancia del análisis de supervivencia frente a las aproximaciones más simples. Si solo hubiera comparado proporciones de extubación habría concluido lo mismo, es cierto, pero habría perdido toda la información sobre cuándo ocurren los eventos y cómo se distribuyen los pacientes censurados.

El enfoque de supervivencia transforma la pregunta de “quién se extuba” en “cuándo se extuba cada quien”, diferencia conceptual que considero fundamental para la epidemiología clínica.

Finalmente, la verificación del supuesto de proporcionalidad me permitió entender que un modelo estadístico no es una caja negra que produce números, sino un andamiaje de supuestos que deben verificarse uno por uno. El día que uno de esos supuestos se rompa, la interpretación cambia radicalmente, por lo que el diagnóstico mediante residuales de Schoenfeld no es un paso accesorio sino una parte constitutiva del análisis.

En este trabajo el supuesto se cumplió, pero el ejercicio me deja la disciplina de siempre verificarlo antes de confiar en cualquier HR estimado.

8.4 Mis conclusiones

El Entrenamiento Muscular Respiratorio, aplicado en la forma y duración del estudio de Sandoval et al. (2019), no fue efectivo para reducir el tiempo de destete de la ventilación mecánica en pacientes con ventilación mecánica por 48 horas o más.

Esta conclusión se sostiene en la convergencia absoluta de cinco enfoques estadísticos independientes, en la validez del supuesto de proporcionalidad del modelo multivariado y en la coherencia con los hallazgos del artículo original.

Dicho lo anterior, mi aprendizaje final del taller trasciende el resultado puntual, pues comprendo ahora que la epidemiología clínica no se trata únicamente de encontrar efectos sino de caracterizarlos con honestidad, incluidos aquellos que resultan nulos.

Un hallazgo negativo bien documentado protege a futuros pacientes de intervenciones sin eficacia demostrada y orienta a los investigadores hacia diseños que sí puedan responder la pregunta adecuadamente.

9 Referencias

Artículo base del taller

Sandoval Moreno, L. M., Casas Quiroga, I. C., Wilches Luna, E. C., y García, A. F. (2019). Eficacia del entrenamiento muscular respiratorio en el destete de la ventilación mecánica en pacientes con ventilación mecánica por 48 o más horas: un ensayo clínico controlado. Medicina Intensiva, 43(2), 79–89. https://doi.org/10.1016/j.medin.2017.11.010