📍 Maestría en Bioestadística. Departamento de Epidemiología Clínica y Bioestadística. Facultad de Medicina. Pontificia Universidad Javeriana, Bogotá.

📚 Créditos. Este material fue desarrollado para el taller “Clasificación binaria en salud: modelos para pronóstico y diagnóstico clínico utilizando R” en el marco.La bases de datos utilizadas durante el taller han sido modificadas exclusivamente con fines académicos.

1 Introducción

Este taller tiene como propósito fortalecer las habilidades para la construcción de modelos de predicción en Salud, específicamente modelos de clasificación binaria. A través de un ejercicio práctico basado en un problema de la vida real, los asistentes aprenderán a seleccionar el paso a paso de los aspectos más relevantes a considerar durante la construcción de un modelo de predicción, así como evaluar su desempeño mediante métricas de discriminación y calibración.

El caso de estudio se centrará en niños con dengue, considerando dos posibles preguntas de investigación:

  • Caso 1 — Diagnóstico: En un niño febril que consulta a urgencias, ¿cuál es la probabilidad de que tenga dengue confirmado en el momento de la evaluación (sí/no)?
  • Caso 2 — Pronóstico: Entre los pacientes ya hospitalizados con dengue confirmado, ¿cuál es la probabilidad de que desarrollen dengue grave durante la hospitalización?

2 Paquetes

Para cargar los paquetes de R corra el siguiente código:

library(Hmisc)        # Funciones avanzadas para análisis estadístico.
library(rms)          # Modelos de regresión avanzada.
library(data.table)   # Manipulación de datos.
library(nlme)         # Modelos lineales y no lineales mixtos.
library(tidyverse)    # Para manipulación, visualización y análisis de datos.
library(car)          # Herramientas para análisis de regresión y diagnóstico.
library(caret)        # Marco unificado para entrenamiento y validación de modelos predictivos.
library(skimr)        # Resúmenes estadísticos.
library(kableExtra)   # Mejora la visualización de tablas en reportes.
library(lmtest)       # Pruebas estadísticas y diagnóstico para modelos lineales.
library(knitr)        # Generación de reportes dinámicos.
library(ROCR)         # Curva ROC
library(ranger)       # RF

3 Caso 1 - Diagnóstico: Dengue confirmado

Objetivo del análisis: Desarrollar un modelo de predicción diagnóstica que estime la probabilidad de que un paciente con síndrome febril agudo tenga dengue confirmado en el momento de la consulta.

Diseño del estudio: Estudio transversal.

Lugar: Hospital público, ubicado en una zona con transmisión activa de dengue en un país de ingreso mediano bajo, durante un período de vigilancia epidemiológica de varios años.

Pacientes: Niños de 15 años o menos, que consultaron de manera consecutiva por fiebre de inicio agudo (1 a 10 días de evolución) sin foco infeccioso evidente, y a quienes se les realizó la prueba de referencia (PCR y/o antígeno NS1) para confirmar o descartar dengue.

Variable de resultado: La variable de resultado es la confirmación de dengue mediante la prueba de referencia, codificada como variable binaria (Sí=1, No=0).

3.1 Datos

El conjunto de datos contiene 700 niños que consultaron de manera consecutiva por fiebre de inicio agudo (1 a 10 días de evolución) sin foco infeccioso evidente, con las siguientes variables:

Diccionario de variables
Variable Descripción Valores
id Identificador del paciente Único
dengue_confirmado Dengue confirmado por PCR/NS1 (desenlace) Sí=1 / No=0
edad Edad Numérico (años)
sexo Sexo Masculino / Femenino
dias_fiebre Días de fiebre antes de la consulta Numérico (días)
zona_endemica Reside o viajó a zona con transmisión activa (14 días) No / Sí
cefalea Cefalea No / Sí
mialgia_artralgia Mialgias o artralgias No / Sí
dolor_retroocular Dolor retroocular No / Sí
exantema Exantema/sarpullido No / Sí
nauseas_vomito Náuseas o vómito No / Sí
petequias Petequias/manchas puntiformes, punto hemorrágico o microhemorragia cutánea No / Sí
torniquete_pos Prueba del torniquete positiva No / Sí
leucocitos Leucocitos al ingreso (x10^3/mm3) Numérico
plaquetas Plaquetas al ingreso (x10^3/mm3) Numérico
hematocrito Hematocrito al ingreso (%) Numérico

4 Pasos del modelamiento

4.1 Paso 1: Lectura de datos

Para cargar la base de datos ejecute el siguiente código:

datos <- readr::read_csv("bd_dengue.csv") # cargar datos

rmarkdown::paged_table(head(datos))

4.1.1 Tamaño de muestra / EPV

La muestra simulada tiene 1600 pacientes con una prevalencia de dengue confirmado de 27.7%, es decir 443 eventos.

El rango mínimo recomendado de número de eventos por parámetro (EPV) es 20-50 para un modelo de regresión logística múltiple. No obstante, se recomienda realizar cálculos formales (Riley et al., BMJ 2020).

Con 14 predictores candidatos, el EPV es de aproximadamente 31.6, dentro del rango mínimo recomendado (20-50 EPV) para un modelo de regresión logística múltiple.

Consultar: Riley R D, Ensor J, Snell K I E, Harrell F E, Martin G P, Reitsma J B et al. Calculating the sample size required for developing a clinical prediction model. BMJ 2020; 368:m441.

4.1.2 Variable de respuesta

En este estudio, se encontró una incidencia de dengue del 27.7% (31.6/1600) en niños que consultaron a urgencias.

datos %>%
  count(dengue_confirmado) %>%
  mutate(porcentaje = round(n / sum(n) * 100, 1))

4.1.3 Potenciales candidatos

Basados en el conocimiento clínico y en una revisión de la literatura, los investigadores del estudio primario decidieron considerar 14 predictores candidatos. Para crear un subconjunto de la base de datos que contenga estas variables junto con la variable de resultado, ejecute el siguiente código:

datos.mod<-datos%>%
  dplyr::select(
    -id # Eliminar la variable identificador
  )
datos.mod <- as.data.frame(datos.mod)

🧠 ¿Cuál es el número máximo de variables predictoras que podrían considerarse en el modelo inicial?

4.2 Paso 2: Evaluación de la calidad de los datos

La inspección general de cada variable puede realizarse mediante la función describe() del paquete rms. Esta función proporciona información sobre el porcentaje de valores faltantes y los principales estadísticos descriptivos de cada variable. A continuación, exploraremos individualmente cada una de ellas.

4.2.1 Variable de respuesta (y)

Para inspeccionar la variable de resultado, ejecute:

# Resumen de la variable respuesta
describe(datos.mod$dengue_confirmado)
## datos.mod$dengue_confirmado 
##        n  missing distinct     Info      Sum     Mean 
##     1600        0        2    0.601      443   0.2769

4.2.2 Predictores categóricos

Para inspeccionar los predictores categóricos, ejecute:

#  Identificar variables categóricas
cat_vars <- datos.mod %>%
   dplyr::select(where(~is.factor(.) | is.character(.))) %>%
  names()

# Resumen de variables categóricas
describe(datos.mod[,cat_vars ])
## datos.mod[, cat_vars] 
## 
##  9  Variables      1600  Observations
## --------------------------------------------------------------------------------
## sexo 
##        n  missing distinct 
##     1600        0        2 
##                               
## Value       Femenino Masculino
## Frequency        783       817
## Proportion     0.489     0.511
## --------------------------------------------------------------------------------
## zona_endemica 
##        n  missing distinct 
##     1600        0        2 
##                       
## Value         No    Si
## Frequency    757   843
## Proportion 0.473 0.527
## --------------------------------------------------------------------------------
## cefalea 
##        n  missing distinct 
##     1560       40        2 
##                     
## Value        No   Si
## Frequency   453 1107
## Proportion 0.29 0.71
## --------------------------------------------------------------------------------
## mialgia_artralgia 
##        n  missing distinct 
##     1600        0        2 
##                       
## Value         No    Si
## Frequency    712   888
## Proportion 0.445 0.555
## --------------------------------------------------------------------------------
## dolor_retroocular 
##        n  missing distinct 
##     1600        0        2 
##                       
## Value         No    Si
## Frequency   1105   495
## Proportion 0.691 0.309
## --------------------------------------------------------------------------------
## exantema 
##        n  missing distinct 
##     1600        0        2 
##                       
## Value         No    Si
## Frequency   1194   406
## Proportion 0.746 0.254
## --------------------------------------------------------------------------------
## nauseas_vomito 
##        n  missing distinct 
##     1600        0        2 
##                     
## Value        No   Si
## Frequency   944  656
## Proportion 0.59 0.41
## --------------------------------------------------------------------------------
## petequias 
##        n  missing distinct 
##     1600        0        2 
##                       
## Value         No    Si
## Frequency   1421   179
## Proportion 0.888 0.112
## --------------------------------------------------------------------------------
## torniquete_pos 
##        n  missing distinct 
##     1600        0        2 
##                       
## Value         No    Si
## Frequency   1271   329
## Proportion 0.794 0.206
## --------------------------------------------------------------------------------

4.2.3 Predictores continuos

Para inspeccionar los predictores continuos, ejecute:

# Identificar variables continuas
cont_vars <- datos.mod %>%
  dplyr::select(-dengue_confirmado) %>%         # Excluye la variable respuesta
  dplyr::select(where(is.numeric)) %>%  # Mantiene solo las numéricas
  names()    # Extrae los nombres de las columnas

# Resumen de variables continuas
describe(datos.mod[,cont_vars ])
## datos.mod[, cont_vars] 
## 
##  5  Variables      1600  Observations
## --------------------------------------------------------------------------------
## edad 
##        n  missing distinct     Info     Mean  pMedian      Gmd      .05 
##     1600        0       15     0.99    7.566      7.5    3.286        3 
##      .10      .25      .50      .75      .90      .95 
##        4        5        8       10       12       12 
##                                                                             
## Value          1     2     3     4     5     6     7     8     9    10    11
## Frequency      5    41    90   117   172   181   182   206   175   157   110
## Proportion 0.003 0.026 0.056 0.073 0.108 0.113 0.114 0.129 0.109 0.098 0.069
##                                   
## Value         12    13    14    15
## Frequency     94    52    17     1
## Proportion 0.059 0.032 0.011 0.001
## 
## For the frequency table, variable is rounded to the nearest 0
## --------------------------------------------------------------------------------
## dias_fiebre 
##        n  missing distinct     Info     Mean  pMedian      Gmd      .05 
##     1600        0       10    0.966    3.959        4    1.903        1 
##      .10      .25      .50      .75      .90      .95 
##        2        3        4        5        6        7 
##                                                                       
## Value          1     2     3     4     5     6     7     8     9    10
## Frequency     82   240   372   381   237   149    86    34    11     8
## Proportion 0.051 0.150 0.232 0.238 0.148 0.093 0.054 0.021 0.007 0.005
## 
## For the frequency table, variable is rounded to the nearest 0
## --------------------------------------------------------------------------------
## leucocitos 
##        n  missing distinct     Info     Mean  pMedian      Gmd      .05 
##     1600        0      118        1    7.266     7.25    2.353      3.8 
##      .10      .25      .50      .75      .90      .95 
##      4.6      5.9      7.3      8.7     10.0     10.8 
## 
## lowest : 1.2  1.4  1.5  1.6  1.7 , highest: 13.2 13.4 13.9 14.5 16  
## --------------------------------------------------------------------------------
## plaquetas 
##        n  missing distinct     Info     Mean  pMedian      Gmd      .05 
##     1600        0      300        1    247.1    247.5    71.89      141 
##      .10      .25      .50      .75      .90      .95 
##      164      205      248      292      329      348 
## 
## lowest :  20  56  62  74  75, highest: 438 446 448 453 455
## --------------------------------------------------------------------------------
## hematocrito 
##        n  missing distinct     Info     Mean  pMedian      Gmd      .05 
##     1574       26      200        1    34.96    34.95    4.593    28.27 
##      .10      .25      .50      .75      .90      .95 
##    29.70    32.20    34.90    37.68    40.20    41.60 
## 
## lowest : 22   22.4 23.6 23.7 24.3, highest: 46.2 46.6 47.5 47.8 48.4
## --------------------------------------------------------------------------------

4.2.4 Valores faltantes

Para finalizar la exploración, es importante analizar la presencia de valores faltantes en cada variable. Para obtener esta información, ejecute el siguiente código:

tabla_faltantes <- data.frame(
  Variable = names(datos.mod),
  N_faltantes = sapply(datos.mod, function(x) sum(is.na(x))),
  Porcentaje = round(sapply(datos.mod, function(x) mean(is.na(x)) * 100), 2))%>%
  dplyr::filter(N_faltantes > 0)

knitr::kable(tabla_faltantes, caption = "Valores faltantes por variable", align = "l")
Valores faltantes por variable
Variable N_faltantes Porcentaje
cefalea cefalea 40 2.50
hematocrito hematocrito 26 1.62

Con base en la salida anterior, las variables con mayor proporción de valores faltantes son cefalea (2.5%) y hematocrito (1.6%). En este punto, resulta importante considerar la estrategia más adecuada para manejar los valores faltantes en los análisis posteriores. De manera general, las recomendaciones incluyen:

  • Para variables con un bajo porcentaje de valores faltantes (<5%), se podría considerar la eliminación de las observaciones con datos faltantes si estamos seguros que el mecanismo de faltantes se puede atribuir completamente al azar MCAR y no se pone en riesgo el tamaño de muestra ni la validez del análisis. En caso contrario, se recomienda considerar métodos de imputación simples, como la imputación por regresión o el algoritmo de \(k\) vecinos más cercanos (KNN).

  • Para variables con un porcentaje moderado a alto de valores faltantes (>5% y \(\leq\) 20%), se recomienda considerar métodos de imputación más sofisticados, como la imputación múltiple (MI).

Por simplicidad vamos a utilizar la estrategia de eliminación de observaciones con datos faltantes para este taller, aunque en la práctica se recomienda evaluar cuidadosamente el mecanismo de faltantes y considerar métodos de imputación más robustos.

datos.mod<-na.omit(datos.mod) # Observaciones con información completa

Total de eventos en la base de datos:

datos.mod %>%
  count(dengue_confirmado) %>%
  mutate(porcentaje = round(n / sum(n) * 100, 1))

🧠 ¿Podemos concluir que eliminar estos casos es una estrategia adecuada?

4.3 Paso 3: Manejo de predictores

  • ¿Es necesario recodificar algún predictor categórico?

  • ¿Conviene aplicar alguna transformación a las variables continuas para facilitar su interpretación?

  • ¿Los predictores continuos presentan una relación lineal con el logit de la variable respuesta?

  • ¿Hay presencia de valores atípicos en los predictores? ¿Cuál sería la mejor forma de tratarlos?

  • ¿Se observa colinealidad entre los predictores? ¿Cómo podría abordarse?

  • ¿Cuál es el número máximo de predictores que sería adecuado considerar al construir el modelo?

4.3.1 Predictores categóricos

En el caso de los predictores categóricos, es importante evaluar la frecuencia y proporción de observaciones en cada categoría antes de incorporarlos al modelo.

arsenal::tableby( ~ cefalea + mialgia_artralgia+dolor_retroocular+exantema+nauseas_vomito+petequias+torniquete_pos, data = datos.mod) %>%
  summary(text = TRUE)
Overall (N=1534)
cefalea
- No 444 (28.9%)
- Si 1090 (71.1%)
mialgia_artralgia
- No 683 (44.5%)
- Si 851 (55.5%)
dolor_retroocular
- No 1057 (68.9%)
- Si 477 (31.1%)
exantema
- No 1145 (74.6%)
- Si 389 (25.4%)
nauseas_vomito
- No 910 (59.3%)
- Si 624 (40.7%)
petequias
- No 1367 (89.1%)
- Si 167 (10.9%)
torniquete_pos
- No 1215 (79.2%)
- Si 319 (20.8%)

🧠 ¿Sería necesario agrupar categorías o crear nuevas variables? ¿Qué criterios utilizarían para tomar esta decisión?

4.3.2 Predictores continuos

En los predictores continuos siempre es recomendable visualizar su distribución. Esto lo podemos realizar por medio del siguiente código:

cont_vars<-cont_vars[-10] # eliminar AST de la lista

# Convertir los datos a formato largo
datos.mod %>% 
  dplyr::select(all_of(cont_vars)) %>%
  pivot_longer(cols = everything(), names_to = "variable", values_to = "valor") %>%
ggplot(., aes(x = variable, y = valor)) +
  geom_boxplot(fill = "skyblue", color = "gray30") +
  facet_wrap(~ variable, scales = "free", ncol = 2) +
  theme_minimal(base_size = 13) +
  theme(
    axis.text.x = element_blank(),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    strip.text = element_text(face = "bold"),
    panel.spacing = unit(1, "lines")
  ) +
  labs(title = "Distribución de predictores continuos", y = "Valor")

4.3.3 Predictores continuos - Linealidad con el logit

Un aspecto importante durante la construcción del modelo de predicción es modelar la verdadera relación entre los predictores continuos y la variable respuesta. En el caso de la regresión logística, se asume que existe una relación lineal entre cada predictor continuo y el logit de la variable respuesta.

Lo anterior, es posible realizarlo por medio del siguiente código, encontrando los porcentajes observados del evento de interés dengue confirmado de acuerdo a los deciles de cada predictor continuo y graficar la relación entre el promedio del predictor y el logit de la proporción del evento. Para esto, corra el siguiente código:

#  Crear una lista de tablas resumen por cada variable continua
tablas_deciles <- map(cont_vars, function(var) {
  
  # Crear cortes en deciles (puede fallar si hay pocos valores únicos)
  decile_breaks <- unique(quantile(datos.mod[[var]], 
                                   probs = seq(0, 1, by = 0.1), 
                                   na.rm = TRUE))
  
  # Generar tabla resumen
  datos.mod %>%
    mutate(decile = cut(.data[[var]], 
                        breaks = decile_breaks, 
                        include.lowest = TRUE)) %>%
    group_by(decile) %>%
    summarise(
      variable = var,
      mean_x = mean(.data[[var]], na.rm = TRUE),
      prop = mean(dengue_confirmado, na.rm = TRUE),
      logit = log(prop / (1 - prop)),
      n = n(),
      .groups = "drop"
    )
})

#  Unir todas las tablas en un solo data frame
tabla_deciles_total <- bind_rows(tablas_deciles)

#  Visualizar la relación entre el promedio del predictor y el logit por deciles

ggplot(tabla_deciles_total, aes(x = mean_x, y = logit)) +
  geom_point( alpha = 0.6) +
  geom_smooth(method = "loess", color = "blue", se = FALSE, span = 0.9) + 
  facet_wrap(~ variable, scales = "free_x") +
  theme_bw() +
  labs(
    title = "Relación entre predictores continuos y logit por deciles",
    x = "Promedio del predictor por decil",
    y = "Logit de la proporción del evento"
  )

4.3.4 Ejemplo para modelar una relación no lineal de un predictor continuo

Inicialmente, se considera el siguiente modelo logístico simple con la variable dias_fiebre como predictor continuo en su forma lineal:

\[ \ln\!\left(\frac{P(y=1|dias.fiebre)}{1-P(y=1|dias.fiebre)}\right) = \beta_0 + \beta_1 \, \text{hct} \]

En R es posible estimar el modelo logístico por medio de la función glm() de la siguiente manera:

model.lineal<- datos.mod%>%
  glm(dengue_confirmado ~ dias_fiebre, data=., family="binomial")

Una posible solución para abordar la no linealidad es categorizar la variable continua. Por ejemplo, se podría categorizar dias_fiebre de acuerdo a los cuartiles:

# 1. Calcular los puntos de corte según los cuartiles
cuartiles <- quantile(datos.mod$dias_fiebre, probs = c(0, 0.25, 0.5, 0.75, 1))

# 2. Crear la variable categórica según los cuartiles
datos.mod$dfiebrecat <- cut(datos.mod$dias_fiebre,
                        breaks = cuartiles,
                        include.lowest = TRUE,
                        labels = c("Q1", "Q2", "Q3", "Q4"))

# 3. Ajustar el modelo de regresión logística
model.cat <- datos.mod %>%
  glm(dengue_confirmado ~ dfiebrecat, data = ., family = "binomial")

No obstante, la categorización de predictores continuos, como dividirlos en cuartiles, terciles o puntos de corte arbitrarios, no es recomendable desde el punto de vista estadístico. Este procedimiento implica una pérdida de información y puede generar puntos de corte artificiales que no tienen una justificación clínica ni biológica.

En la actualidad existen métodos más flexibles como los splines cubicos restringidos (RCS, por su nombre en inglés, restricted cubic spline). Este enfoque será utilizado para modelar la relación entre dias_fiebre y el logit de dengue severo.

4.3.4.1 Splines cubicos restringidos

  • Son funciones flexibles que permiten modelar relaciones no lineales entre una variable continua y la variable de respuesta.
  • Forzan que la función sea lineal más allá de los puntos extremos (nodos).
  • Requieren la selección de un número adecuado de nodos (knots).

Matemáticamente una función spline \(k\) nodos en \(t_1, t_2,..., t_k\), se define como:

\[f(x)= \beta_0 + \beta_1 x_1 + \beta_2 x_2+ ...+\beta_{k-1}x_{k-1},\]

donde, \(x_1=x\) y para \(j=2, ...,k-1\), se tiene que

\[x_{j}=(x-t_{j-1})_+^3 -\frac{(x-t_{k-1})_+^3(t_k-t_{j-1})}{(t_k-t_{k-1})} + \frac{(x-t_{k})_+^3(t_{k-1}-t_{j-1})}{(t_k-t_{k-1})},\] \(x_j\) será lineal cuando \(x \geq t_k\),

donde,

\[(u)_+=\begin{cases} u, & \text{si } u > 0, \\ 0, & \text{si } u \leq 0 \end{cases}\]

Aplicando este enfoque en la variable dias_fiebre, se pueden considerar diferentes números de nodos. A continuación, se ajustan dos modelos utilizando RCS con 3 y 4 nodos, respectivamente:

model.rcs.k3 <-datos.mod%>%
  dplyr::select(dengue_confirmado, dias_fiebre)%>%
  glm(dengue_confirmado ~ rcs(dias_fiebre,3), data = .,family="binomial")

model.rcs.k4 <-datos.mod%>%
  dplyr::select(dengue_confirmado, dias_fiebre)%>%
  glm(dengue_confirmado ~ rcs(dias_fiebre,4), data = ., family="binomial")

Para facilitar la comparación es posible graficar los tres enfoques de modelamiento por medio del siguiente código:

# Muestra: División en deciles y % del evento por decil

tabladfiebre<- tabla_deciles_total %>%
  filter(variable == "dias_fiebre")

datos.dfiebre <-datos.mod%>%
  dplyr::select(dengue_confirmado, dias_fiebre)

#Predicción por cada modelo
datos.dfiebre$lineal<- predict(model.lineal, data=datos.hct, type="link")
datos.dfiebre$cat <-predict(model.cat,data=datos.hct, type="link")
datos.dfiebre$rcs3 <-predict(model.rcs.k3, data=datos.hct, type="link")
datos.dfiebre$rcs4 <-predict(model.rcs.k4, data=datos.hct, type="link")

# Unificación de datos en formato largo

datos.dfiebre<-datos.dfiebre%>%
   pivot_longer(
    cols = c(lineal, cat, rcs3, rcs4),  # nombres reales de tus columnas"
    names_to = "modelo",               # nuevo nombre de columna para el nombre del modelo
    values_to = "pred"                 # nuevo nombre de columna para el valor predicho
)

ggplot() +
  # LOESS para lineal y rcs
  geom_smooth(
    data = datos.dfiebre %>% filter(modelo %in% c("lineal", "rcs3", "rcs4")),
    aes(x = dias_fiebre, y = pred, color = modelo),
    method = "loess", se = FALSE
  ) +
  # Línea directa para cat
  geom_line(
    data = datos.dfiebre %>% filter(modelo == "cat"),
    aes(x = dias_fiebre, y = pred, color = modelo)
  ) +
  # Puntos adicionales de tablehct
  geom_point(
    data = tabladfiebre, aes(x = mean_x, y = logit),
    inherit.aes = FALSE, shape = 21, size = 2, fill = "black"
  ) +
  labs(x = "Dias de fiebre", y = "logit", color = " ") + theme_gray()

Tambien se recomienda calcular el criterio de información de Akaike (AIC) de cada enfoque. Recordemos que este criterio se define como:

\[AIC=−2×LL+2k,\]

donde:

La log-verosimilitud (LL) mide el ajuste de los datos al modelo y \(k\) es el número de paramétros estimados. Un AIC más bajo indica un modelo mejor (balance entre ajuste y simplicidad).

Comparación de modelos (LogLik y AIC)
Modelo LogLik AIC
Lineal -909.904 1823.808
Categórico -909.898 1827.795
RCS (3 nudos) -908.628 1823.256
RCS (4 nudos) -908.251 1824.502
Note:
El AIC penaliza la complejidad del modelo; un valor menor indica mejor ajuste relativo.
1 LogLik: Log-likelihood; AIC: Criterio de Información de Akaike.
* Modelos ajustados con familia binomial (logit).

4.3.5 Correlación entre predictores

Otro aspecto relevante es evaluar la correlación entre los predictores continuos. La presencia de alta correlación entre predictores puede indicar multicolinealidad, lo que puede afectar el modelo.

corrplot::corrplot(cor(datos.mod[, cont_vars], 
          use = "pairwise.complete.obs"), # Usar valores excluyendo NA's
          method="number" # Mostrar valores numéricos
          )

🧠 Antes de ajustar un Random Forest, ¿es necesario preocuparnos por la no linealidad y la correlación entre los predictores?

4.4 Paso 4: Especificación del modelo completo

4.4.1 Herramienta 1: Regresión logística

El modelo inicial tiene un total de 9 predictores, como se muestra a continuación:

\[ \begin{array}{lcl} \text{logit}\big(P(Y = 1)\big) &=& \beta_0 + \beta_1 \text{edad} + \beta_2 \text{sexo} + \beta_3 \text{dias_fiebre} \\ &+& \beta_4 \text{zona_endemica} + \beta_5 \text{cefalea} + \beta_6 \text{mialgia_artralgia} \\ &+& \beta_7 \text{dolor_retroocular} + \beta_8 \text{exantema} \\ &+& \beta_{9} \text{nauseas_vomito} + \beta_{10} \text{petequias} + \beta_{11} \text{torniquete_pos} \\ &+& \beta_{12} \text{leucocitos} + \beta_{13} \text{plaquetas} + \beta_{14} \text{hematocrito} \end{array} \] Vamos definir el conjunto de datos con los predictores seleccionados para el modelo inicial. Para esto, ejecute el siguiente código:

datos.mod<-datos.mod%>% 
    dplyr:: select(-dfiebrecat) # Eliminar variable categórica creada para el ejemplo

El modelo completo inicial (M0) es el siguiente:

dd <- datadist(datos.mod) # Configurar datadist para rms
options(datadist = "dd")

model.MO.lg <- lrm(dengue_confirmado~edad +
              sexo+
              dias_fiebre+ 
              zona_endemica+
              cefalea +
              mialgia_artralgia +
              dolor_retroocular +
              exantema +
              nauseas_vomito +
              petequias +
              torniquete_pos +
              leucocitos +
              plaquetas +
              hematocrito,
              data=datos.mod, y=T, x=T)

model.MO.lg
## Logistic Regression Model
## 
## lrm(formula = dengue_confirmado ~ edad + sexo + dias_fiebre + 
##     zona_endemica + cefalea + mialgia_artralgia + dolor_retroocular + 
##     exantema + nauseas_vomito + petequias + torniquete_pos + 
##     leucocitos + plaquetas + hematocrito, data = datos.mod, x = T, 
##     y = T)
## 
##                        Model Likelihood       Discrimination    Rank Discrim.    
##                              Ratio Test              Indexes          Indexes    
## Obs          1534    LR chi2     474.44       R2       0.383    C       0.834    
##  0           1103    d.f.            14     R2(14,1534)0.259    Dxy     0.668    
##  1            431    Pr(> chi2) <0.0001    R2(14,929.7)0.391    gamma   0.668    
## max |deriv| 3e-06                             Brier    0.143    tau-a   0.270    
## 
##                      Coef    S.E.   Wald Z Pr(>|Z|)
## Intercept             3.1787 0.7450   4.27 <0.0001 
## edad                 -0.0173 0.0246  -0.70 0.4814  
## sexo=Masculino       -0.0189 0.1366  -0.14 0.8901  
## dias_fiebre          -0.0483 0.0397  -1.22 0.2242  
## zona_endemica=Si      0.6850 0.1390   4.93 <0.0001 
## cefalea=Si           -0.0044 0.1499  -0.03 0.9768  
## mialgia_artralgia=Si  1.0941 0.1443   7.58 <0.0001 
## dolor_retroocular=Si  0.3607 0.1469   2.46 0.0140  
## exantema=Si           0.9634 0.1539   6.26 <0.0001 
## nauseas_vomito=Si     0.0938 0.1385   0.68 0.4983  
## petequias=Si          0.4623 0.2103   2.20 0.0280  
## torniquete_pos=Si     1.4029 0.1598   8.78 <0.0001 
## leucocitos           -0.3447 0.0360  -9.57 <0.0001 
## plaquetas            -0.0158 0.0012 -12.82 <0.0001 
## hematocrito           0.0150 0.0169   0.89 0.3736

🧠 ¿Es este el momento adecuado para realizar selección de variables? ¿Por qué?

4.4.1.1 Validación interna del modelo inicial

La validación del modelo se realizará por medio de la técnica de remuestreo Boostrap con el principal objetivo de corregir las medidas de desempeño del modelo por el optimismo. Esto puede realizarse por con la función validate:

# 200 remuestreos 
vinterna.M0.lg<- rms::validate(model.MO.lg,  B = 200)

Aunque la anterior función reporta múltiples indicadores, vamos a enfocarnos en las siguientes métricas de ajuste global, discriminación y calibración:

  • B: Brier Score: Menor puntuación, mejor el modelo. Para facilitar su interpretación se recomienda trabajar con el Brier escalado, el cual toma valores entre 0 y 100% (mejor modelo) y está dado por:

\[Brier_{escalado}=1-\frac{Brier}{Brier_{max}},\] donde, \(Brier_{max}=\hat{p}\times(1-\hat p)^2 + (1-\hat p) \times \hat p^2\).

  • R2: Nagelkerke \(R^2\).

  • Dxy: Somer’s \(D_{xy}\) index. Datos Binarios:\(D_{xy}=2 \times(AUC-1/2)\). Valores cercanos a 1 discriminan perfectamente. El AUC puede ser obtenido por \(AUC: Dxy/2 + 0.5\).

  • Intercept: Intercepto de la curva de calibración. Valores cercanos a 0 son mejores. No corresponde a calibration-in-the-large

  • Slope: Pendiente de la curva de calibración. Valores cercanos a 1 son mejores. Slope < 1 \(\rightarrow\) sobreajuste (overfitting)

vinterna.M0.lg[c(1,2,3,4,9),]
##           index.orig  training        test     optimism index.corrected   n
## Dxy        0.6675340 0.6784316  0.65955532  0.018876288      0.64865773 200
## R2         0.3827286 0.3961035  0.37383683  0.022266700      0.36046185 200
## Intercept  0.0000000 0.0000000 -0.03064858  0.030648577     -0.03064858 200
## Slope      1.0000000 1.0000000  0.94926238  0.050737616      0.94926238 200
## B          0.1425695 0.1399845  0.14443466 -0.004450157      0.14701969 200

4.4.2 Herramienta 2: Random Forest

Para usar la implementación del paquete caret, se requiere que el desenlace sea un factor con etiquetas válidas (no 0/1), por lo tanto, vamos a declarar el significado de estos números en el programa:

datos.rf <- datos.mod %>%
  mutate(
    dengue_confirmado = factor(dengue_confirmado,levels = c(1, 0),
      labels = c("Si", "No")
    )
  )
levels(datos.rf$dengue_confirmado)
## [1] "Si" "No"

Dado que el tamaño muestral es moderado (1534 pacientes y 431 eventos) este algoritmo podría no ser la mejor alternativa inicialmente. El interés principal de este ejercicio es ilustrar el proceso de optimización, utilizando toda la muestra disponible combinada con remuestreo (validación cruzada para la optimización de hiperparámetros, y bootstrap para la estimación del optimismo del modelo final.

Aspectos importantes a considerar en el RF

  • El desbalanceo de clases o la baja proporción de eventos puede afectar la capacidad del modelo para aprender patrones de la clase minoritaria. Esto puede llevar a un sesgo hacia la clase mayoritaria, resultando en un modelo que predice principalmente la clase mayoritaria y tiene un rendimiento deficiente en la clase minoritaria.

  • En estos casos no se debe elegir los hiperparámetros del modelo basados en la exactitud/accuracy del modelo (clasificación correcta), podría ser una mejor opción el área bajo la curva. Se recomienda el indicador AUC (Area Under the Curve) que considera la sensibilidad y la especificidad para identificar casos con y sin dengue.

  • Los hiperparametros más importantes para optimizar son:

    • mtry: número de predictores considerados aleatoriamente en cada partición. Es, con diferencia, el hiperparámetro más influyente; el valor por defecto de Breiman (\(\sqrt(p)\) para clasificación) es un buen punto de partida.
    • min.node.size (tamaño mínimo del nodo terminal): controla la profundidad/complejidad de cada árbol; valores más altos generan árboles más simples y regularizados (mayor sesgo, menor varianza).

Para el proceso de optimización de hiperparámetros, se recomienda utilizar validación cruzada repetida (repeated cross-validation). Esto permite evaluar el rendimiento del modelo en diferentes particiones de los datos y obtener una estimación de su desempeño. En R puede implementarse con:

set.seed(2026)

ctrl_rf <- trainControl(
  method          = "repeatedcv",
  number          = 5,
  repeats         = 2,
  classProbs      = TRUE,
  summaryFunction = twoClassSummary,
  savePredictions = "final",
  verboseIter     = FALSE
)

grid_rf <- expand.grid(
  mtry           = 2:8,
  min.node.size  = c(1, 5, 10),
  splitrule      = "gini"
)

Ahora, vamos a proceder a correr el modelo con el fin de identificar los mejores hiperparámetros por medio de la función train. Este proceso puede tardar varios minutos dependiendo de la capacidad de procesamiento de su computadora.

set.seed(2026) # fijar semilla
modelo_rf <- train(
  dengue_confirmado ~ .,
  data          = datos.rf,
  method        = "ranger",
  metric        = "ROC",       #  ROC AUC como métrica de optimización
  tuneGrid      = grid_rf,
  trControl     = ctrl_rf,
  num.trees     = 500,
  importance    = "permutation"
)

Los mejores hiperparámetros encontrados por el algoritmo de optimización son:

modelo_rf$bestTune
kable(modelo_rf$results %>% arrange(desc(ROC)) %>% head(10),
      digits = 3, caption = "Diez mejores combinaciones de hiperparámetros (por AUC-ROC en validación cruzada)")
Diez mejores combinaciones de hiperparámetros (por AUC-ROC en validación cruzada)
mtry min.node.size splitrule ROC Sens Spec ROCSD SensSD SpecSD
2 10 gini 0.801 0.261 0.954 0.023 0.034 0.020
3 10 gini 0.800 0.370 0.927 0.024 0.031 0.015
2 5 gini 0.799 0.284 0.954 0.023 0.030 0.017
3 5 gini 0.799 0.381 0.922 0.023 0.049 0.014
4 10 gini 0.798 0.402 0.913 0.025 0.042 0.014
5 10 gini 0.797 0.428 0.908 0.024 0.037 0.012
3 1 gini 0.796 0.369 0.927 0.023 0.044 0.013
4 5 gini 0.796 0.407 0.912 0.024 0.037 0.011
2 1 gini 0.796 0.276 0.949 0.022 0.037 0.016
4 1 gini 0.794 0.408 0.908 0.024 0.038 0.009 
plot(modelo_rf)

Finalmente, podemos ver la importancia de variables calculada por el modelo de Random Forest. La importancia de las variables se puede evaluar mediante la métrica de importancia de permutación, que mide la disminución en el rendimiento del modelo cuando se permuta aleatoriamente una variable específica. Una mayor disminución indica que la variable es más importante para el modelo.

importancia_rf <- varImp(modelo_rf, scale = TRUE)$importance %>%
  rownames_to_column("Variable") %>%
  arrange(desc(Overall))

ggplot(importancia_rf, aes(x = reorder(Variable, Overall), y = Overall)) +
  geom_col(fill = "steelblue") +
  coord_flip() +
  theme_minimal() +
  labs(title = "Importancia de variables (permutación) — Random Forest optimizado",
       x = NULL, y = "Importancia relativa (0-100)")

4.4.2.1 Validación interna del modelo

La función rms::validate() no tiene un método implementado para objetos ranger/train (random forest).

En su lugar, debemos reimplementar el mismo algoritmo de remuestreo bootstrap con corrección por optimismo. Esta es exactamente la lógica que usa validate() por dentro: (1) ajustar el modelo en cada muestra bootstrap, (2) evaluar su desempeño dentro de esa muestra bootstrap (“entrenamiento”), (3) evaluar ese mismo modelo bootstrap en la muestra original (“prueba”), (4) el optimismo es la diferencia entre ambos, y (5) el desempeño corregido es el aparente menos el optimismo promedio en los \(B\) remuestreos.

discriminacion_calibracion <- function(y, prob) {
  if (is.factor(y)) y <- ifelse(y == "Si", 1, 0) # admite y numérico (0/1) o factor Si/No

  prob <- pmin(pmax(prob, 1e-4), 1 - 1e-4)
  logit_p <- qlogis(prob)

  auc_val    <- as.numeric(pROC::auc(pROC::roc(y, prob, quiet = TRUE)))
    # glm(y ~ logit_p), exactamente como lo hace rms::validate() 
  recalibracion <- glm(y ~ logit_p, family = binomial)
  intercept_val <- unname(coef(recalibracion)[1])
  slope_val     <- unname(coef(recalibracion)[2])
  b_val <- mean((prob - y)^2) # Brier score, igual que la columna "B" de rms::validate()

  c(AUC = auc_val, Intercept = intercept_val, Slope = slope_val, B = b_val)
}

validar_bootstrap <- function(datos, y_var, ajustar_fun, predecir_fun, B = 200, semilla = 2026) {
  set.seed(semilla)
  n <- nrow(datos)
  y_orig <- datos[[y_var]]

  modelo_orig <- ajustar_fun(datos)
  prob_orig   <- predecir_fun(modelo_orig, datos)
  aparente    <- discriminacion_calibracion(y_orig, prob_orig)

  optimismos <- matrix(NA_real_, nrow = B, ncol = length(aparente),
                        dimnames = list(NULL, names(aparente)))

  for (b in 1:B) {
    idx   <- sample.int(n, n, replace = TRUE)
    dboot <- datos[idx, ]
    mboot <- tryCatch(ajustar_fun(dboot), error = function(e) NULL)
    if (is.null(mboot)) next
    y_boot <- dboot[[y_var]]

    perf_boot_boot <- discriminacion_calibracion(y_boot, predecir_fun(mboot, dboot))
    perf_boot_orig <- discriminacion_calibracion(y_orig, predecir_fun(mboot, datos))
    optimismos[b, ] <- perf_boot_boot - perf_boot_orig
  }
  optimismo_medio <- colMeans(optimismos, na.rm = TRUE)
  corregido <- aparente - optimismo_medio
  list(aparente = aparente, optimismo = optimismo_medio, corregido = corregido,
       B_efectivo = sum(!is.na(optimismos[, 1])))
}
ajustar_lr <- function(datos) {
  dd_tmp <- datadist(datos)
  assign("dd_tmp", dd_tmp, envir = .GlobalEnv)
  old <- options(datadist = "dd_tmp")
  on.exit(options(old), add = TRUE)
  lrm(dengue_confirmado ~ edad + sexo + dias_fiebre + zona_endemica + cefalea +
        mialgia_artralgia + dolor_retroocular + exantema + nauseas_vomito +
        petequias + torniquete_pos + leucocitos + plaquetas + hematocrito,
      data = datos, x = TRUE, y = TRUE)
}
predecir_lr <- function(modelo, datos) as.numeric(predict(modelo, newdata = datos, type = "fitted"))
val_lr <- validar_bootstrap(datos.mod, "dengue_confirmado", ajustar_lr, predecir_lr, B = 200)
auc_rms   <- vinterna.M0.lg["Dxy", "index.corrected"] / 2 + 0.5
slope_rms <- vinterna.M0.lg["Slope", "index.corrected"]
int_rms   <- vinterna.M0.lg["Intercept", "index.corrected"]
b_rms     <- vinterna.M0.lg["B", "index.corrected"]

tabla_comparacion <- data.frame(
  Metrica = c("AUC", "Slope", "Intercept", "B"),
  `rms::validate()` = c(auc_rms, slope_rms, int_rms, b_rms),
  `validar_bootstrap()` = c(val_lr$corregido["AUC"], val_lr$corregido["Slope"],
                             val_lr$corregido["Intercept"], val_lr$corregido["B"]),
  check.names = FALSE
)
kable(tabla_comparacion, digits = 4,
      caption = "Regresión logística: comparación entre rms::validate() y validar_bootstrap() (B=200, misma base n=1534)") %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE) %>%
  row_spec(0, bold = TRUE, color = "white", background = "#0D47A1")
Regresión logística: comparación entre rms::validate() y validar_bootstrap() (B=200, misma base n=1534)
Metrica rms::validate() validar_bootstrap()
AUC 0.8243 0.8262
Slope 0.9493 NA
Intercept -0.0306 NA
B 0.1470 NA

Usamos los hiperparámetros óptimos ya encontrados (modelo_rf$bestTune) y reajustamos un random forest en cada muestra bootstrap considerando las predicciones realizadas fueran de la bolsa OOB:

# Random Forest evaluado con sus propias predicciones OOB (out-of-bag)
set.seed(2026)
mejor_rf <- modelo_rf$bestTune
rf_oob <- ranger(
  dengue_confirmado ~ ., data = datos.rf,
  mtry = mejor_rf$mtry, min.node.size = mejor_rf$min.node.size,
  splitrule = as.character(mejor_rf$splitrule),
  num.trees = 500, probability = TRUE
)
prob_rf_oob <- rf_oob$predictions[, "Si"]  # predicciones OOB, no las del ajuste completo

oof_rf <- discriminacion_calibracion(datos.mod$dengue_confirmado, prob_rf_oob)

prob_lr <- as.numeric(predict(model.MO.lg, type = "fitted"))
Comparación de las formas de estimar el desempeño: aparente vs. fuera de muestra
Random Forest
Regresión logística
Metrica RF - OOB (fuera de muestra) RL - Bootstrap corregido
AUC 0.7984 0.8262
Intercepto 0.4242 NA
Pendiente 1.5872 NA
Brier 0.1584 NA
Note:
RL: regresión logística. RF: random forest

🧠 ¿Cuál de los dos modelos LR o RF tiene mejor desempeño predictivo para el diagnóstico de dengue en niños menores de 15 años?

4.4.3 Curvas de calibración

Para explorar la calibración del modelo de una manera más exhaustiva se recomienda la construcción de curvas de calibración. Estas se construyen estimando un modelo de regresión logística condicionado a las probabilidades predichas (o estimadas) por el modelo propuesto.

# --- Curva de calibración (loess flexible) ---
graficar_calibracion <- function(y, prob, color, add = FALSE, ...) {
  lo <- loess(y ~ prob, degree = 2)
  xs <- seq(min(prob), max(prob), length.out = 200)
  pred_lo <- pmin(pmax(predict(lo, newdata = data.frame(prob = xs)), 0), 1)
  if (!add) {
    plot(xs, pred_lo, type = "l", col = color, lwd = 2, xlim = c(0, 1), ylim = c(0, 1),
         xlab = "Probabilidad predicha", ylab = "Proporción observada",
         main = "Curva de calibración", ...)
    abline(0, 1, lty = 2, col = "gray40")
  } else {
    lines(xs, pred_lo, col = color, lwd = 2)
  }
}
graficar_calibracion(datos.mod$dengue_confirmado, prob_lr, "#0D47A1")
graficar_calibracion(datos.mod$dengue_confirmado, prob_rf_oob, "#C62828", add = TRUE)
rug(prob_lr, col = adjustcolor("#0D47A1", alpha.f = 0.4), side = 1)
rug(prob_rf_oob, col = adjustcolor("#C62828", alpha.f = 0.4), side = 3)
legend("topleft", legend = c("Regresión logística", "Random Forest (OOB)", "Calibración ideal"),
       col = c("#0D47A1", "#C62828", "gray40"), lty = c(1, 1, 2), lwd = 2, bty = "n", cex = 0.85)

par(mfrow = c(1, 1))

4.4.4 Curva ROC

Otro gráfico de gran relevancia es la curva de característica operativa (ROC, por sus siglas en inglés), la cual permite visualizar el área bajo la curva (AUC) y, con ello, evaluar la capacidad discriminativa del modelo. La curva ROC muestra la relación entre la tasa de verdaderos positivos y la tasa de falsos positivos. Una curva ubicada por encima de la línea roja refleja una mayor capacidad de discriminación del modelo. En R esto puede realizarse con:

roc_lr     <- pROC::roc(datos.mod$dengue_confirmado, prob_lr, quiet = TRUE)
roc_rf_oob <- pROC::roc(datos.mod$dengue_confirmado, prob_rf_oob, quiet = TRUE)

par(mfrow = c(1, 2))

# --- Curva ROC ---
plot(roc_lr, col = "#0D47A1", lwd = 2, main = "Curva ROC")
plot(roc_rf_oob, col = "#C62828", lwd = 2, add = TRUE)
legend("bottomright",
       legend = c(paste0("Regresión logística (AUC=", round(as.numeric(pROC::auc(roc_lr)), 3), ")"),
                  paste0("Random Forest, OOB (AUC=", round(as.numeric(pROC::auc(roc_rf_oob)), 3), ")")),
       col = c("#0D47A1", "#C62828"), lwd = 2, bty = "n", cex = 0.85)

4.5 Paso 5: Reducción del modelo

**Objetivo:* identificar si existe un subconjunto de los predictores que ofrezca un buen desempeño predictivo, mientras se conservan aquellos con relevancia clínica.

La selección del modelo final debe realizarse conservando las decisiones realizadas para modelar o categorizar cada variable predictora. En este punto podría considerarse la posibilidad de eliminar predictores que no aporten información adicional al modelo, basándose en criterios estadísticos y clínicos. Aquí, se realizará este proceso con el modelo de regresión logística, utilizando la función fastbw() del paquete rms, que implementa un procedimiento de selección hacia atrás (backward selection) basado en el criterio de información de Akaike (AIC). Este procedimiento permite identificar un subconjunto de predictores que optimiza el balance entre ajuste y complejidad del modelo. Adiciomalmente, se realizará un análisis de remuestreo bootstrap para evaluar la estabilidad de la selección de variables, es decir, qué variables tienden a ser seleccionadas con mayor frecuencia en diferentes muestras bootstrap.

# número de remuestreos bootstrap
B <- 200

# guardar las variables retenidas en cada remuestreo
vars <- list()
set.seed(123)
for (i in 1:B) {
  # 1. muestra bootstrap
  datos_boot <- datos.mod[sample(1:nrow(datos.mod), replace = TRUE), ]

  # 2. ajustar el modelo completo (mismas variables que model.MO.lg)
  model_boot <- lrm(
    dengue_confirmado ~ edad + sexo + dias_fiebre + zona_endemica +
      cefalea + mialgia_artralgia + dolor_retroocular + exantema +
      nauseas_vomito + petequias + torniquete_pos +
      leucocitos + plaquetas + hematocrito,
    data = datos_boot, x = TRUE, y = TRUE
  )

  # 3. selección backward (AIC)
  bw_boot <- fastbw(model_boot, rule = "aic")

  # 4. guardar variables seleccionadas en esta réplica
  vars[[i]] <- bw_boot$names.kept
}

# juntar resultados de las B réplicas
vars_all <- unlist(vars)

# frecuencia de selección
freq <- table(vars_all)

# porcentaje (as.numeric() evita que quede como objeto "table", que es lo que
# generaba la columna espuria "Porcentaje.vars_all")
freq_percent <- round(100 * as.numeric(freq) / B, 1)

# tabla final
tabla_final <- data.frame(
  Variable   = names(freq),
  Frecuencia = as.numeric(freq),
  Porcentaje = freq_percent
)

tabla_final %>%
  arrange(desc(Porcentaje)) %>%
  knitr::kable(
    caption = "Frecuencia de selección de variables en el bootstrap (B = 200)",
    align = "lcc", row.names = FALSE
  ) %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE) %>%
  row_spec(0, bold = TRUE, color = "white", background = "#0D47A1")
Frecuencia de selección de variables en el bootstrap (B = 200)
Variable Frecuencia Porcentaje
exantema 200 100.0
leucocitos 200 100.0
mialgia_artralgia 200 100.0
plaquetas 200 100.0
torniquete_pos 200 100.0
zona_endemica 199 99.5
dolor_retroocular 101 50.5
petequias 65 32.5
hematocrito 19 9.5
dias_fiebre 17 8.5
nauseas_vomito 13 6.5
edad 9 4.5
cefalea 6 3.0
sexo 6 3.0

🧠 Con los resultados obtenidos, ¿qué variables considera que deberían incluirse en el modelo final? ¿Por qué?

Cualquiera que sea el camino elegido, el modelo final debe reajustarse sobre la totalidad de los datos (no sobre una submuestra) y luego validarse internamente con la misma lógica de corrección por optimismo que ya aplicamos a model.MO.lg.

4.6 Paso 6: Selección de un punto de corte mediante curvas de decisión

Hasta ahora hemos evaluado los modelos con medidas estadísticas de discriminación (AUC) y calibración (Intercept, Slope). Sin embargo, en la práctica clínica es necesario decidir: ¿a partir de qué probabilidad predicha se considera a un paciente “positivo” para dengue (por ejemplo, para iniciar manejo hospitalario u órdenes adicionales)?

4.6.1 Análisis de curvas de decisión (Decision Curve Analysis, DCA)

El beneficio neto (net benefit) propuesto por Vickers y Elkin (2006) permite comparar estrategias de decisión (usar el modelo vs. tratar/estudiar a todos vs. no tratar/estudiar a nadie) incorporando explícitamente el costo relativo entre falsos positivos y falsos negativos, a través de una probabilidad umbral (\(p_t\)) que refleja ese costo relativo:

\[\text{Beneficio neto}(p_t) = \frac{VP}{n} - \frac{FP}{n}\times\left(\frac{p_t}{1-p_t}\right)\]

donde \(VP\) y \(FP\) son los verdaderos y falsos positivos al clasificar como “positivo” a todo paciente con probabilidad predicha \(\geq p_t\), y \(n\) es el tamaño de la muestra. El término \(p_t/(1-p_t)\) es el odds del umbral, y representa cuántos falsos positivos un clínico está dispuesto a aceptar por cada verdadero positivo adicional identificado. Un \(p_t=0.20\), por ejemplo, implica que se aceptarían hasta 4 falsos positivos (pacientes sin dengue estudiados innecesariamente) por cada verdadero positivo (paciente con dengue correctamente identificado).

Se compara el beneficio neto del modelo contra dos estrategias de referencia:

  • Tratar/estudiar a todos: \(\text{BN}_{todos}(p_t) = \text{prevalencia} - (1-\text{prevalencia})\times\frac{p_t}{1-p_t}\)
  • Tratar/estudiar a ninguno: \(\text{BN}_{ninguno}=0\) (por definición, no se generan ni verdaderos ni falsos positivos)

Un modelo es clínicamente útil en un umbral \(p_t\) solo si su beneficio neto supera ambas estrategias de referencia en ese punto.

net_benefit_modelo <- function(y, prob, pt) {
  n <- length(y)
  vp <- sum(prob >= pt & y == 1)
  fp <- sum(prob >= pt & y == 0)
  vp / n - (fp / n) * (pt / (1 - pt))
}
net_benefit_tratar_todos <- function(y, pt) {
  prev <- mean(y)
  prev - (1 - prev) * (pt / (1 - pt))
}
pts <- seq(0.01, 0.70, by = 0.01)
dca <- data.frame(
  pt = pts,
  `Regresión logística` = sapply(pts, function(p) net_benefit_modelo(datos.mod$dengue_confirmado, prob_lr, p)),
  `Random Forest`       = sapply(pts, function(p) net_benefit_modelo(datos.mod$dengue_confirmado, prob_rf_oob, p)),
  `Estudiar a todos`    = sapply(pts, function(p) net_benefit_tratar_todos(datos.mod$dengue_confirmado, p)),
  `Estudiar a ninguno`  = 0,
  check.names = FALSE
)

Un umbral específico solo puede justificarse clínicamente, no estadísticamente (Vickers & Elkin, 2006): debe reflejar cuánto más grave se considera no identificar un caso de dengue frente al costo de estudiar innecesariamente a un niño sano (una prueba confirmatoria adicional, una consulta de seguimiento).

4.7 Paso 7: Validación externa

No se dispone de datos para realizar este paso en el presente ejercicio práctico. Sin embargo, el propósito de esta etapa es evaluar el rendimiento del modelo propuesto en datos distintos a los utilizados para su construcción o en nuevas poblaciones. Para más información, se recomienda consultar:Steyerberg EW. Validation of prediction models. In: Clinical prediction models: a practical approach to development, validation, and updating. 2nd ed. Cham: Springer; 2019. Chapter 17, p. 329-343.

4.8 Paso 8: Presentación del modelo

El desarrollo de aplicaciones interactivas en Shiny ofrece la posibilidad de implementar el modelo en una plataforma web, permitiendo que los profesionales de la salud ingresen valores individuales de los predictores y obtengan de inmediato la probabilidad estimada del evento. Un ejemplo sencillo se muestra a continuación:

5 Referencias

  • Breiman L. Random Forests. Machine Learning. 2001;45(1):5-32.
  • Harrell FE Jr. Regression Modeling Strategies: With Applications to Linear Models, Logistic and Ordinal Regression, and Survival Analysis. 2nd ed. Springer; 2015.
  • Steyerberg EW. Clinical Prediction Models: A Practical Approach to Development, Validation, and Updating. 2nd ed. Springer; 2019.