Caso completo: análisis estadístico publicable

r
estadistica
Un análisis clínico end-to-end con bajo peso al nacer (birthwt). Descriptiva, comparaciones univariadas, regresión múltiple, diagnóstico, tamaños de efecto y tabla final con gtsummary.

El cierre de la ruta

Has completado 11 tutoriales sobre estadística aplicada. Esta entrega no introduce conceptos nuevos: junta todo lo anterior en un análisis clínico real, desde la pregunta inicial hasta la tabla final publicable.

Usamos el dataset MASS::birthwt, bajo peso al nacer en 189 madres del Baystate Medical Center (Springfield, Massachusetts, 1986). Variables clave:

  • bwt: peso al nacer en gramos (la respuesta principal).
  • low: bajo peso al nacer, < 2500 g (indicador binario, derivado de bwt).
  • age: edad de la madre.
  • lwt: peso de la madre antes del embarazo (libras).
  • smoke: fumó durante el embarazo (0/1).
  • race: raza (1 = blanca, 2 = negra, 3 = otra).
  • ht: antecedentes de hipertensión.
  • ui: irritabilidad uterina.
  • ptl: partos prematuros previos.
  • ftv: visitas médicas en el primer trimestre.

La pregunta clínica

“¿Qué factores maternos predicen el peso al nacer en esta cohorte? ¿En qué magnitud afecta cada uno, controlando por los demás?”

Pregunta clásica de epidemiología: tres niveles de análisis (descriptivo, univariado por exposición, multivariado).

Paso 0: Preparar los datos

library(MASS)
library(dplyr)
library(ggplot2)
library(gtsummary)
library(effectsize)

datos <- birthwt |>
  mutate(
    low   = factor(low,   labels = c("Normal", "Bajo")),
    smoke = factor(smoke, labels = c("No fuma", "Fuma")),
    race  = factor(race,  labels = c("Blanca", "Negra", "Otra")),
    ht    = factor(ht,    labels = c("No", "Sí")),
    ui    = factor(ui,    labels = c("No", "Sí")),
    bwt_kg = bwt / 1000
  )

Convertir los códigos numéricos en factores con etiquetas legibles antes del análisis ahorra trabajo en cada paso siguiente.

Paso 1: Descriptiva inicial

Primero, una mirada al dataset entero:

library(skimr)
skim(datos)

skim() te muestra n, n_missing, distribuciones, el primer comando al recibir cualquier dataset.

Después, descriptiva por exposición principal (fumar). Es el primer paso clásico:

tabla1 <- datos |>
  select(bwt, age, lwt, race, ht, ui, smoke) |>
  tbl_summary(
    by = smoke,
    label = list(
      bwt  ~ "Peso al nacer (g)",
      age  ~ "Edad materna (años)",
      lwt  ~ "Peso materno previo (lb)",
      race ~ "Raza",
      ht   ~ "Hipertensión",
      ui   ~ "Irritabilidad uterina"
    )
  ) |>
  add_p() |>
  add_overall() |>
  bold_labels()

Esta es la “Tabla 1” del paper. Tres columnas: total, no-fuma, fuma. Una fila por covariable. p-value por fila. Mediana (Q1, Q3) para continuas, n (%) para categóricas.

Lectura: ¿hay diferencias basales entre los grupos? Si fumadoras y no-fumadoras difieren en edad, raza o peso materno, hay confusión potencial. Lo controlaremos en el modelo multivariado.

Paso 2: Comparación univariada del outcome principal

Pregunta: “¿Difiere el peso al nacer entre fumadoras y no-fumadoras?”

Inspección visual primero (siempre):

ggplot(datos, aes(smoke, bwt_kg)) +
  geom_jitter(width = 0.2, alpha = 0.5) +
  geom_boxplot(width = 0.3, alpha = 0.5, outlier.shape = NA) +
  labs(
    title = "Peso al nacer por exposición materna al tabaco",
    x = NULL,
    y = "Peso al nacer (kg)"
  ) +
  theme_minimal()

Si las distribuciones se ven razonablemente simétricas y sin outliers extremos (verifica con QQ-plot dentro de cada grupo), t-test de Welch:

test_smoke <- t.test(bwt ~ smoke, data = datos)
test_smoke
#>
#>  Welch Two Sample t-test
#>
#> data:  bwt by smoke
#> t = 2.7299, df = 170.1, p-value = 0.007003
#> alternative hypothesis: true difference in means between group No fuma and group Fuma is not equal to 0
#> 95 percent confidence interval:
#>   78.57486 488.97860
#> sample estimates:
#> mean in group No fuma   mean in group Fuma
#>              3055.696             2771.919

Diferencia de medias: 284 g, IC 95 %: [79, 489] g. Los hijos de fumadoras pesan en media 284 g menos al nacer.

Tamaño de efecto:

cohens_d(bwt ~ smoke, data = datos)
#> Cohen's d |       95% CI
#> ------------------------
#> 0.41      | [0.11, 0.70]

Cohen’s d = 0.41 (efecto pequeño-mediano). El IC excluye cero, así que el efecto es consistente. Reportable:

“Las hijas de mujeres que fumaron durante el embarazo pesaron 284 g menos al nacer (IC 95 %: 79-489. Cohen’s d = 0.41 [IC 95 %: 0.11-0.70]. T = 2.73, p = 0.007).”

Tres datos para decir lo mismo: diferencia bruta, tamaño de efecto, significancia.

Paso 3: Modelo multivariado

Univariado no basta, fumadoras y no-fumadoras pueden diferir en otras características. El modelo multivariado controla por las covariables:

modelo <- lm(bwt ~ smoke + age + lwt + race + ht + ui, data = datos)
summary(modelo)

Coeficiente de smoke:

coef(modelo)["smokeFuma"]
#> smokeFuma
#> -355.6298

Tras controlar por edad, peso materno, raza, hipertensión e irritabilidad uterina, fumar se asocia con un peso al nacer −356 g (más extremo que el univariado, lo que sugiere que algunas covariables ocultaban el efecto verdadero).

IC para el coeficiente:

confint(modelo)["smokeFuma", ]
#>      2.5 %     97.5 %
#> -550.62049 -160.63906

Paso 4: Diagnóstico

par(mfrow = c(2, 2))
plot(modelo)

Verifica los cuatro paneles (ver tutorial de diagnóstico):

  • Residuals vs Fitted: ¿alguna curva o forma de embudo?
  • QQ-plot: ¿residuos aproximadamente normales?
  • Scale-Location: ¿varianza constante?
  • Residuals vs Leverage: ¿algún punto con Cook’s distance > 0.5?

Si todo se ve razonable, el modelo es defendible. Si hay problemas, los abordas (transformar Y, eliminar puntos justificadamente, reformular).

Diagnóstico moderno alternativo:

library(performance)
check_model(modelo)

Genera todos los paneles con interpretación textual.

Paso 5: Tamaño de efecto multivariado

Para reportar la contribución relativa de cada predictor:

standardize_parameters(modelo)

Lectura: cada coeficiente expresado en SD. Permite comparar la importancia relativa de variables con escalas distintas (edad en años vs peso en libras).

Para el R² del modelo entero:

summary(modelo)$adj.r.squared
#> [1] 0.2018

El modelo explica ~20 % de la varianza en peso al nacer. No es altísimo, pero es coherente con la complejidad biológica de la variable (muchos factores no medidos: genética fetal, nutrición materna, etc.).

Paso 6: Tabla final con gtsummary

El producto final, listo para incluir en el manuscrito:

tabla_modelo <- modelo |>
  tbl_regression(
    label = list(
      smoke ~ "Fumar durante embarazo",
      age   ~ "Edad materna (años)",
      lwt   ~ "Peso materno previo (lb)",
      race  ~ "Raza",
      ht    ~ "Hipertensión",
      ui    ~ "Irritabilidad uterina"
    )
  ) |>
  bold_labels()

tabla_modelo

Salida: tabla con coeficiente, IC 95 % y p-value por predictor, con labels legibles. Lista para as_gt() |> gtsave("tabla_modelo.docx").

Para una vista combinada univariada + multivariada (patrón epidemiológico clásico):

uni <- tbl_uvregression(
  datos |> select(bwt, smoke, age, lwt, race, ht, ui),
  method = lm,
  y = bwt
)

tbl_merge(
  list(uni, tabla_modelo),
  tab_spanner = c("**Univariado**", "**Multivariado**")
)

Lo que has hecho

En este caso completo has aplicado los once tutoriales anteriores en orden:

  1. Descriptiva (skim, tbl_summary).
  2. Distribuciones, verificación visual de los supuestos antes del t-test.
  3. t-test de Welch para la comparación univariada principal.
  4. Chi-cuadrado implícito a través de add_p() para las covariables categóricas.
  5. Correlación / regresión simple para entender la relación inicial.
  6. Regresión múltiple para ajustar por confusores.
  7. Diagnóstico con los cuatro paneles de plot(lm).
  8. Tamaño de efecto (Cohen’s d univariado + betas estandarizados multivariado).
  9. Reporte profesional con gtsummary.

Y todo enmarcado en una pregunta clínica concreta, no en un ejercicio abstracto.

¿Quieres ir más a fondo?

El libro Estadística para datos reales con R (en preparación) desarrolla esta misma materia con:

  • Un caso clínico mucho más largo, con datos longitudinales reales.
  • Modelos mixtos (lme4) para diseños con clusters.
  • Análisis de supervivencia con survival y survminer.
  • Capítulo dedicado a la comunicación: cómo defender un análisis en peer review.

Has terminado la Ruta 3

Felicidades. Si quieres seguir aprendiendo R, las rutas que naturalmente continúan son:

  • Machine Learning con tidymodels: del describir e inferir al predecir.
  • Reproducibilidad con Quarto: para que tus análisis sean replicables y publicables.

Ambas están en el listado de Rutas.