1 Introducción

Para llevar a cabo un afianciamiento del contenido de la materia Estádistica Aplicada de la Pontificia Universidad Javeriana Cali en el periodo 2025-1 se lleva a cabo una implementación de dos métodos diferentes: clustering jerárquico y regresión multilineal, cada uno para abordar un problema abierto a escogencia del equipo de trabajo.

El análisis se enmarca en la metodología CRISP-DM(Cross-Industry Standard Process for Data Mining), que guía el desarrollo de proyectos de ciencia de datos mediante seis fases: comprensión del negocio, comprensión de los datos, preparación de los datos, modelado, evaluación y despliegue. Seguido por su flexibilidad, fácil personalización, por el enfásis que hace en los aspectos fundamentales para el planteamiento y desarrollo de un proyecto y por recomendación del profesor.

2 Clustering para la Caracterización de Ciudades

2.1 Entendimiento del Problema

En años recientes (2020-2023) Colombia ha tenido una reducción en su competitividad relativa de acuerdo con IMD World Competitiveness Ranking, además, de acuerdo con el informe Trayectorias: Prosperidad y reducción de la pobreza en el territorio colombiano del World Bank Group las ciudades más pobres suelen estar más aisladas, considerando la accidentada topografía de Colombia, por lo tanto se presentan grandes desafios en infrastructura de conectividad terrestre relativa a otros paises, también indica que se establece una concentración de la actividad económica en las ciudades que representan mayores aglomeraciones poblacionales.

Por ello, entender cómo agrupar las ciudades en función de similitudes estructurales es clave para mejorar la efectividad de las intervenciones gubernamentales.

2.1.1 Objetivo del estudio

  • Reducir la complejidad del análisis nacional para el diseño de politicas públicas.

2.2 Datos

Se toma la base de datos del informe que evalúa el desempeño competitivo de las 32 ciudades capitales de departamento en Colombia realizado por la Cámara de Comercio de la ciudad de Bucaramanga, el Consejo Privado de Competitividad (CPC) y la Universidad del Rosario, llevado a cabo por ofrecer una herramienta clave para el análisis y la toma de decisiones en materia de desarrollo regional. También, de dicho informe se logra obtener el Índice de Competitividad de Ciudades (ICC) 2024 y la población de cada capital de departamento en el 2024. Con Google Maps se halla la distancia y tiempo de viaje a media noche del día 13/5/2025 desde las ciudades más pobladas del país Bogotá, Medellín AM, Cali AM, Barranquilla AM.

# ——————————————
# 1. Preparar datos
# ——————————————
library(tidyverse)
library(cluster)
library(factoextra)
library(readxl)

# Leer la base desde archivo Excel
Base <- read_excel("C:/Users/mario/Downloads/Base.xlsx")

# Seleccionar las variables de interés
vars <- c("ICC",
          "Índice_Tiempo_Distancia",
          "Población")

# Extraer el subconjunto con esas columnas
x <- Base[, vars]

# Estandarizar las variables numéricas
x <- scale(x)

# Convertir a data frame para compatibilidad con clustering
x <- as.data.frame(x)

# Usar el nombre de las ciudades como nombre de fila
row.names(x) <- Base[[1]]
head(x)
##                        ICC Índice_Tiempo_Distancia  Población
## Arauca          -0.9080571              0.07873301 -0.4828798
## Armenia          0.6179615             -0.50384697 -0.3480476
## Barranquilla AM  0.8267980             -0.50579987  0.7542448
## Bogotá D.C.      2.0226467             -0.50579987  4.5729713
## Bucaramanga AM   0.9608666             -0.28697126  0.3091366
## Cali AM          1.0899095             -0.50579987  0.9973042

2.3 Entendimiento de los datos

Se usa directamente el índice de Competitividad entre Ciudades (ICC) que toma en cuenta los 4 factores Condiciones habilitantes, Capital Humano, Eficiencia de los Mercados y Ecosistema Innovador, la población de las ciudades se toma, como ya se menciono para considerar el efecto de las aglomeraciones sobre la caracterización de las ciudades y las distancias y tiempos de las ciudades a las cuatro más grandes se toma para tener en cuenta el factor de aglomeración y distancia juntos y así considerar el aspecto geográfico del pais.

2.4 Preparación de los datos

Inicialmente no hay que hacer mucho en este sentido, ya que, se trata de una base de datos que ya ha procesado los datos el ICC, la población se toma como tal y con las distancias y tiempos se crea un indicador de tiempo distancia como “peso logístico” para penalizar trayectos largos y lentos.

2.5 Modelado

Se va a implementar le método de clustering jerárquico para obtener la representación en dendograma de las observaciones sin importar el número de variables usadas (en este caso 3). Para la disimilitud se usa la distancia euclidiana sobre la correlación entre variables (que sigue patrones y halla perfiles entre ciudades de acuerdo a esto) debido al factor de aglomeración que se esta teniendo en cuenta en este estudio, lo cual hace que las magnitudes sean importantes.

Se realiza inicialmente el clustering con los tres tipos de enlaces más populares entre los estadísticos:

Completo: Máxima disimilitud entre conglomerados. Calcule todas las disimilitudes por pares entre las observaciones del conglomerado A y las del conglomerado B, y registre la mayor de estas disimilitudes.

Simple: Disimilitud mínima entre conglomerados. Calcule todas las disimilitudes por pares entre las observaciones del conglomerado A y las del conglomerado B, y registre la menor de estas disimilitudes. El enlace simple puede resultar en conglomerados extensos y rezagados, en los que las observaciones individuales se fusionan una a una.

Promedio: Disimilitud media entre conglomerados. Calcule todas las disimilitudes por pares entre las observaciones del conglomerado A y las del conglomerado B, y registre el promedio de estas disimilitudes.

# ——————————————
# 2. Distancias
# ——————————————
dist_euc <- dist(x, method = "euclidean")

# ——————————————
# 3. Dendrogramas
# ——————————————
par(mfrow = c(1, 3), mar = c(2, 4, 2, 2))

hc.euc.complete <- hclust(dist_euc, method = "complete")
hc.euc.average  <- hclust(dist_euc, method = "average")
hc.euc.single   <- hclust(dist_euc, method = "single")

plot(hc.euc.complete, main = "Euclidiana / Complete", xlab = "", sub = "", cex = .8)
plot(hc.euc.average , main = "Euclidiana / Average",  xlab = "", sub = "", cex = .8)
plot(hc.euc.single  , main = "Euclidiana / Single",   xlab = "", sub = "", cex = .8)

Para determinar cuanto clusters se deberían tomar dadas las observaciones y sus factores primero se usa la métrica del ancho medio de silueta:

# ——————————————
# 4. Análisis de silueta
# ——————————————
sil_width <- function(dist_mat, hc, k) {
  clus <- cutree(hc, k = k)
  ss <- silhouette(clus, dist_mat)
  mean(ss[, "sil_width"])
}

ks <- 2:31
res <- tibble(
  k = rep(ks, times = 3),
  method = rep(c("Euc-Complete", "Euc-Average", "Euc-Single"), each = length(ks)),
  sil = NA_real_
)

for(i in seq_along(res$method)) {
  m <- res$method[i]; k <- res$k[i]
  hc_obj <- switch(m,
                   "Euc-Complete" = hc.euc.complete,
                   "Euc-Average"  = hc.euc.average,
                   "Euc-Single"   = hc.euc.single)
  res$sil[i] <- sil_width(dist_euc, hc_obj, k)
}

# Visualizar silueta
ggplot(res, aes(x = k, y = sil, color = method)) +
  geom_line() + geom_point() +
  labs(title = "Silueta media vs k", x = "Número de clusters", y = "Ancho medio de silueta") +
  theme_minimal()

Esta gráfica compara el número de clusters con respecto al ancho medio de silueta como medida de cuán bien se ajusta cada elemento a su propio clúster, comparado con otros clústeres. La idea central al “maximizar la silueta media”es elegir el número de grupos k que optimiza, en promedio, esa medida de ajuste.

El ancho de silueta se calcula así:

Considerando que se quieren más de cuatro clusteres para que estudio pueda llegar a tener en cuenta las disimilitudes de las ciudades al momento de plantear políticas se aprecia que el enlace simple tiene valores muy bajos en el rango hasta entorno a los 20 clusteres, lo cual ya se considera un número muy elevado de grupos para establecer políticas que puedan orientarse a un fin general sin mucha experiencia previa, por lo tanto se descarta.

En los tipos de enlace completo y promedio se nota una región en el gráfico que es parecido a una meseta por lo que se toma un valor \(k=9\) que se encuentra aproximadamente a la mitad de la misma.

Para complementar la métrica del ancho medio de silueta se usa el método del codo usando el pseudo-WCSS (Within-Cluster Sum Squares), ya que la primera penaliza fuertemente clústeres grandes con subestructura (puede preferir muchos clústeres pequeños) y no evalúa la geometría global de los clústeres, ni su estabilidad. Pseudo, ya que normalmente se hace uso de el en el método de clustering por k-means.

El pseudo-WCSS es la suma de las distancias al centroide dentro de cada clúster, mide qué tan compactos son los cluústeres.

# ——————————————
# 5. Método del codo (WCSS)
# ——————————————
wcss <- function(x, hc, k) {
  clus <- cutree(hc, k = k)
  sum(sapply(split(as.data.frame(x), clus), function(sub) {
    sum(rowSums((scale(sub, scale = FALSE))^2))
  }))
}

ks <- 2:20
wcss_df <- expand.grid(method = c("Euc-Complete", "Euc-Average", "Euc-Single"), k = ks) %>%
  rowwise() %>%
  mutate(WCSS = wcss(x,
                     switch(method,
                            "Euc-Complete" = hc.euc.complete,
                            "Euc-Average"  = hc.euc.average,
                            "Euc-Single"   = hc.euc.single),
                     k)) %>%
  ungroup()

# Visualizar codo
ggplot(wcss_df, aes(x = k, y = WCSS, color = method)) +
  geom_line() + geom_point() +
  labs(title = "Método del codo (WCSS)", x = "Número de clusters", y = "Suma intra-cluster (WCSS)") +
  theme_minimal()

De la gráfica anterior se logra apreciar que para las formas de enlace promedio y completa a partir de \(k=9\) la pendiente se atenua significativamente más visualmente que en los cambios de \(k\) anteriores, con esto se define que se hará uso de este número de clústeres.

Se debe destacar que con el enlace completo y el promedio se han agrupado las mismas ciudades en los nueve clústeres elegidos, a través de la visualización de los siguientes dendogramas se considera más sensato el enlace completo, por la estructuración interna de sus clústeres.

# ——————————————
# 6. Dendrogramas con cortes óptimos
# (elige los "k" más razonables luego de ver gráficas)
# ——————————————
k1 <- 9
k2 <- 9
k3 <- 15

par(mfrow = c(1, 3), mar = c(2, 4, 2, 2))

plot(hc.euc.complete, main = "Euc / Complete", xlab = "", sub = "")
rect.hclust(hc.euc.complete, k = k1, border = 2:4)

plot(hc.euc.average, main = "Euc / Average", xlab = "", sub = "")
rect.hclust(hc.euc.average, k = k2, border = 2:4)

plot(hc.euc.single, main = "Euc / Single", xlab = "", sub = "")
rect.hclust(hc.euc.single, k = k3, border = 2:4)

# ——————————————
# 7. Clusters finales
# ——————————————
clusters <- list(
  EucComplete = cutree(hc.euc.complete, k1),
  EucAverage  = cutree(hc.euc.average,  k2),
  EucSingle   = cutree(hc.euc.single,   k3)
)

# Ver resultados de un método específico
print(clusters$EucComplete)
##                Arauca               Armenia       Barranquilla AM 
##                     1                     2                     3 
##           Bogotá D.C.        Bucaramanga AM               Cali AM 
##                     4                     2                     3 
##             Cartagena             Cúcuta AM             Florencia 
##                     2                     5                     1 
##                Ibagué               Inírida               Leticia 
##                     2                     6                     7 
##          Manizales AM           Medellín AM                  Mitú 
##                     2                     8                     6 
##                 Mocoa              Montería                 Neiva 
##                     1                     5                     2 
##                 Pasto            Pereira AM               Popayán 
##                     2                     2                     2 
##        Puerto Carreño                Quibdó              Riohacha 
##                     6                     1                     1 
##            San Andrés San José del Guaviare           Santa Marta 
##                     9                     1                     2 
##             Sincelejo                 Tunja            Valledupar 
##                     1                     2                     5 
##         Villavicencio                 Yopal 
##                     5                     5
# ——————————————
# 8. WCSS final para los k seleccionados
# ——————————————
wcss_values <- tibble(
  method = c("Euc-Complete", "Euc-Average", "Euc-Single"),
  k = c(k1, k2, k3),
  WCSS = c(wcss(x, hc.euc.complete, k1),
           wcss(x, hc.euc.average,  k2),
           wcss(x, hc.euc.single,   k3))
)
print(wcss_values)
## # A tibble: 3 × 3
##   method           k  WCSS
##   <chr>        <dbl> <dbl>
## 1 Euc-Complete     9  3.19
## 2 Euc-Average      9  3.13
## 3 Euc-Single      15  3.11

También se ha calculado el WCSS (a menor valor para los tres tipos de enlace con el número de clústeres elegidos, de lo que se evidencia si bien el WCSS del enlace simple es menor al de los otros dos tipos de enlace, cuenta con muchos más clusteres, lo cual no es deseado en la clasificación deseada.

2.6 Evaluación

Para finalizar, se muestra el dendograma del tipo de enlace promedio elegido en el que se puede apreciar la distinción de:

Grupo 1: Florencia, Sincelejo, Riohacha, San José del Guaviare, Mocoa, Arauca, Quibdó (Población baja, baja alta competitividad, acceso difícil a ciudad de aglomeración)

Grupo 2 : Armenia, Popayán, Tunja, Manizales AM, Pereira AM, Bucaramanga AM, Cartagena, Ibagué, Neiva, Pasto, Santa Marta (Población intermedia, alta competitividad, acceso fácil a ciudad de aglomeración)

Grupo 3: Barranquilla AM, Cali AM (Alta población, alta competitividad, ciudad de aglomeración)

Grupo 4: Bogotá D.C. (Muy alta población, muy alta competitividad, ciudad de aglomeración)

Grupo 5: Cúcuta AM, Sincelejo, Yopal, Montería, Villavicencio, Valledupar (Población intermedia, medio-baja competitividad, acceso intermedio a ciudad de aglomeración)

Grupo 6: Mitú, Inírida, Puerto Carreño (Baja población, baja competitividad, acceso difícil a ciudad de aglomeración)

Grupo 7: Leticia (Población intermedia, media competitividad, acceso fácil a ciudad de aglomeración)

Grupo 8: Medellín AM (Alta población, muy alta competitividad, ciudad de aglomeración)

Grupo 9: San Andrés (Baja población, medio-baja competitividad, acceso muy difícil a ciudad de aglomeración).

# Asignar clústeres
k1 <- 9
grupos <- cutree(hc.euc.complete, k = k1)

# Plot del dendrograma
plot(hc.euc.complete, main = "Euc / Complete", xlab = "", sub = "", cex = 0.7)
rect.hclust(hc.euc.complete, k = k1, border = rainbow(k1))

Con la breve descripción de los grupos se puede apreciar que el modelo empleado tcumple los objetivos de reducir la complejidad del análisis nacional para el diseño de politicas públicas, ya que las ciudades que se han juntado tienen características similares.

2.7 Despliegue

Se calculan los centroides para tener forma de determinar una distancia de la ciudad ingresada a los grupos establecidos:

# Obtener las asignaciones de cluster
cluster_assignments <- cutree(hc.euc.complete, k = 9)

# Calcular los centroides para cada cluster
centroids <- x %>%
  mutate(Cluster = cluster_assignments) %>%
  group_by(Cluster) %>%
  summarise(across(everything(), mean))

# Mostrar los centroides
print(centroids)
## # A tibble: 9 × 4
##   Cluster     ICC Índice_Tiempo_Distancia Población
##     <int>   <dbl>                   <dbl>     <dbl>
## 1       1 -0.777                   -0.257    -0.447
## 2       2  0.680                   -0.422    -0.173
## 3       3  0.958                   -0.506     0.876
## 4       4  2.02                    -0.506     4.57 
## 5       5 -0.0891                  -0.271    -0.168
## 6       6 -1.75                     1.12     -0.527
## 7       7 -1.52                     4.24     -0.512
## 8       8  1.47                    -0.506     2.15 
## 9       9 -0.238                    2.23     -0.512

Clasificación de Ciudad según ICC, Población, Distancia y Tiempo

2.8 Conclusión

Este estudio ha demostrado que la clasificación de ciudades colombianas basada en indicadores como el ICC, población y factores de distancia-tiempo permite identificar grupos con características socioeconómicas y geográficas claramente diferenciadas. La segmentación en nueve clusters refleja la realidad heterogénea del país, donde ciudades como Bogotá y Medellín destacan por su alta competitividad y conectividad, mientras que otras localidades más aisladas o con menor población presentan desafíos significativos para el desarrollo.

La integración de métodos estadísticos con herramientas interactivas facilita no solo el análisis sino también la interpretación dinámica de los resultados, permitiendo a los interesados —desde planificadores hasta académicos— tomar decisiones informadas que impulsen el desarrollo regional. En suma, esta clasificación no es solo un ejercicio académico, sino una base sólida para orientar políticas públicas focalizadas y estrategias de crecimiento territorial más eficientes.

Todo esto con el uso del clustering

3 Regresión Logística para Predicción en Salud

3.1 Entendimiento del Problema

El modelo de riesgos proporcionales de Cox es una herramienta estadística ampliamente utilizada en el análisis de supervivencia. Este modelo permite analizar el tiempo hasta que ocurre un evento de interés (como la muerte, el fallo de un dispositivo o la recaída de un paciente), considerando al mismo tiempo múltiples variables explicativas.

En este proyecto se aplica el modelo de Cox para predecir la probabilidad de muerte en pacientes con insuficiencia cardíaca, utilizando un conjunto de datos clínicos. Este conjunto incluye información como la edad del paciente, presión arterial, fracción de eyección, niveles de creatinina y otros factores relevantes.

La insuficiencia cardíaca es una condición médica crónica en la que el corazón no puede bombear sangre de manera eficiente para satisfacer las necesidades del cuerpo. Esta patología es una de las principales causas de hospitalización y mortalidad en todo el mundo, especialmente en pacientes de edad avanzada o con comorbilidades.

Desde el punto de vista clínico y operativo, es crucial poder predecir el riesgo de muerte de un paciente con insuficiencia cardíaca a partir de sus características médicas y demográficas. Esta predicción puede apoyar la toma de decisiones médicas, como priorizar tratamientos, asignar recursos hospitalarios o definir estrategias de seguimiento personalizado.

El objetivo de este análisis es utilizar un modelo estadístico adecuado para analizar el tiempo hasta la muerte de los pacientes, considerando múltiples variables clínicas. Para ello, se empleará el modelo de regresión de Cox, que permite estudiar cómo diferentes factores (como la edad, fracción de eyección, niveles de creatinina, etc.) influyen sobre el riesgo de muerte a lo largo del tiempo.

Este modelo será aplicado sobre un conjunto de datos clínicos reales, que contiene registros de 299 pacientes con insuficiencia cardíaca, recolectados durante su estancia hospitalaria. A través del análisis, se busca:

3.1.1 Objetivos del estudio

  • Identificar las variables más asociadas al riesgo de muerte.

  • Estimar el efecto de cada variable sobre el riesgo.

  • Predecir probabilidades relativas de mortalidad entre diferentes perfiles de pacientes.

Con este enfoque, se contribuye a una medicina más proactiva y basada en datos, que puede mejorar los resultados clínicos y optimizar los procesos hospitalarios.

3.2 Datos

Para llevar a cabo el análisis de supervivencia con el modelo de Cox, se utiliza un conjunto de datos clínicos que contiene información sobre 299 pacientes con insuficiencia cardíaca. Este dataset incluye variables médicas, demográficas y de seguimiento, todas ellas recolectadas durante el ingreso hospitalario de cada paciente.

## 📄 Descripción general del dataset:

- **Observaciones (filas):** 299 pacientes

- **Variables (columnas):** 13 en total

- **Variable objetivo (evento):** `DEATH_EVENT` (1 si el paciente falleció, 0 si sobrevivió durante el seguimiento)

- **Variable de tiempo:** `time` (días de seguimiento hasta el evento o censura)

## 🧬 Variables disponibles:

library(knitr)
library(kableExtra)
library(readr)

hfcrd <- read_excel("C:/Users/mario/OneDrive/Desktop/Estadística_aplicada/Proyecto Final/heart_failure_clinical_records_dataset.xlsx")

# Crear la tabla
variables <- data.frame(
  Variable = c(
    "age", "anaemia", "creatinine_phosphokinase", "diabetes",
    "ejection_fraction", "high_blood_pressure", "platelets",
    "serum_creatinine", "serum_sodium", "sex", "smoking",
    "time", "DEATH_EVENT"
  ),
  Tipo = c(
    "Numérica", "Binaria", "Numérica", "Binaria",
    "Numérica", "Binaria", "Numérica",
    "Numérica", "Numérica", "Binaria", "Binaria",
    "Numérica", "Binaria"
  ),
  Descripción = c(
    "Edad del paciente (años)",
    "1 si el paciente tiene anemia, 0 si no",
    "Nivel de CPK en sangre (enzima cardíaca)",
    "1 si el paciente es diabético, 0 si no",
    "Porcentaje de sangre bombeada por el corazón en cada latido",
    "1 si el paciente tiene hipertensión, 0 si no",
    "Número de plaquetas en sangre (mil/mL)",
    "Nivel de creatinina en sangre (mg/dL), indicador de función renal",
    "Nivel de sodio en sangre (mEq/L)",
    "1 para hombre, 0 para mujer",
    "1 si el paciente fuma, 0 si no",
    "Días de seguimiento",
    "1 si el paciente murió durante el seguimiento, 0 si no"
  )
)

# Mostrar la tabla en formato bonito
kable(variables, format = "html", caption = "Descripción de las variables del dataset") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = F)
Descripción de las variables del dataset
Variable Tipo Descripción
age Numérica Edad del paciente (años)
anaemia Binaria 1 si el paciente tiene anemia, 0 si no
creatinine_phosphokinase Numérica Nivel de CPK en sangre (enzima cardíaca)
diabetes Binaria 1 si el paciente es diabético, 0 si no
ejection_fraction Numérica Porcentaje de sangre bombeada por el corazón en cada latido
high_blood_pressure Binaria 1 si el paciente tiene hipertensión, 0 si no
platelets Numérica Número de plaquetas en sangre (mil/mL)
serum_creatinine Numérica Nivel de creatinina en sangre (mg/dL), indicador de función renal
serum_sodium Numérica Nivel de sodio en sangre (mEq/L)
sex Binaria 1 para hombre, 0 para mujer
smoking Binaria 1 si el paciente fuma, 0 si no
time Numérica Días de seguimiento
DEATH_EVENT Binaria 1 si el paciente murió durante el seguimiento, 0 si no

3.3 📊 Resumen estadístico

# Librería para estadísticas más detalladas (opcional)
# install.packages("psych")
library(psych)

# Estadísticas detalladas: media, sd, min, max, etc.
describe(hfcrd)
##                          vars   n      mean       sd   median   trimmed
## age                         1 299     60.83    11.89     60.0     60.22
## anaemia                     2 299      0.43     0.50      0.0      0.41
## creatinine_phosphokinase    3 299    581.84   970.29    250.0    365.49
## diabetes                    4 299      0.42     0.49      0.0      0.40
## ejection_fraction           5 299     38.08    11.83     38.0     37.43
## high_blood_pressure         6 299      0.35     0.48      0.0      0.32
## platelets                   7 299 263358.03 97804.24 262000.0 256730.09
## serum_creatinine            8 299      1.39     1.03      1.1      1.19
## serum_sodium                9 299    136.63     4.41    137.0    136.82
## sex                        10 299      0.65     0.48      1.0      0.68
## smoking                    11 299      0.32     0.47      0.0      0.28
## time                       12 299    130.26    77.61    115.0    129.28
## DEATH_EVENT                13 299      0.32     0.47      0.0      0.28
##                               mad     min      max    range  skew kurtosis
## age                         14.83    40.0     95.0     55.0  0.42    -0.22
## anaemia                      0.00     0.0      1.0      1.0  0.28    -1.93
## creatinine_phosphokinase   269.83    23.0   7861.0   7838.0  4.42    24.53
## diabetes                     0.00     0.0      1.0      1.0  0.33    -1.90
## ejection_fraction           11.86    14.0     80.0     66.0  0.55     0.00
## high_blood_pressure          0.00     0.0      1.0      1.0  0.62    -1.62
## platelets                65234.40 25100.0 850000.0 824900.0  1.45     6.03
## serum_creatinine             0.30     0.5      9.4      8.9  4.41    25.19
## serum_sodium                 4.45   113.0    148.0     35.0 -1.04     3.98
## sex                          0.00     0.0      1.0      1.0 -0.62    -1.62
## smoking                      0.00     0.0      1.0      1.0  0.76    -1.42
## time                       105.26     4.0    285.0    281.0  0.13    -1.22
## DEATH_EVENT                  0.00     0.0      1.0      1.0  0.76    -1.42
##                               se
## age                         0.69
## anaemia                     0.03
## creatinine_phosphokinase   56.11
## diabetes                    0.03
## ejection_fraction           0.68
## high_blood_pressure         0.03
## platelets                5656.17
## serum_creatinine            0.06
## serum_sodium                0.26
## sex                         0.03
## smoking                     0.03
## time                        4.49
## DEATH_EVENT                 0.03
library(dplyr)
library(corrplot)

# Seleccionar solo las variables clave numéricas
vars_relevantes <- c(
  "age", "creatinine_phosphokinase", "ejection_fraction", 
  "platelets", "serum_creatinine", "serum_sodium", 
  "time", "DEATH_EVENT"
)

hfcrd_cl <- hfcrd %>%
  select(all_of(vars_relevantes)) %>%
  na.omit()

# Renombrar para claridad en español
colnames(hfcrd_cl) <- c(
  "Edad", "CPK", "Fracción", "Plaquetas",
  "Creatinina", "Sodio", "Tiempo", "Muerte"
)

# Calcular matriz de correlación
cor_matrix <- cor(hfcrd_cl)

# Graficar matriz sin números, solo colores y etiquetas claras
corrplot(cor_matrix, method = "color",
         col = colorRampPalette(c("blue", "white", "red"))(200),
         tl.col = "black", tl.srt = 45,
         number.cex = 0,  # sin números en cuadritos
         cl.pos = "r",
         title = "Matriz de Correlación de Variables Clave",
         mar = c(0, 0, 2, 0))

3.4 Entendimiento de los datos/Preparación de los datos

3.4.1 Interpretación de la matriz de correlación

  • Los colores azules indican correlaciones negativas, rojos positivas, y blancos ausencia de correlación.
  • Por ejemplo, la fracción de eyección muestra correlación negativa con la muerte, es decir, a menor fracción mayor riesgo.
  • Variables como la creatinina sérica y edad muestran correlación positiva con la muerte, señalando que valores altos aumentan el riesgo.
  • Estas asociaciones sugieren que el modelo de Cox debe incluir variables clínicamente relevantes y estadísticamente asociadas con el evento para mejorar la predicción.

3.5 Preparación de los Datos

Selección, limpieza, construcción de variables y formateo.

## Cargar librerías necesarias
library(readr)
library(dplyr)

## Leer los datos
hfcrd <- read_excel("C:/Users/mario/OneDrive/Desktop/Estadística_aplicada/Proyecto Final/heart_failure_clinical_records_dataset.xlsx")

## Recodificar variables binarias a factores con etiquetas en español
hfcrd <- hfcrd %>%
  mutate(
    anaemia = factor(ifelse(anaemia == 1, "Sí", "No")),
    diabetes = factor(ifelse(diabetes == 1, "Sí", "No")),
    high_blood_pressure = factor(ifelse(high_blood_pressure == 1, "Sí", "No")),
    sex = factor(ifelse(sex == 1, "Hombre", "Mujer")),
    smoking = factor(ifelse(smoking == 1, "Sí", "No")),
    DEATH_EVENT = factor(ifelse(DEATH_EVENT == 1, "Muerte", "Censurado"))
  )

## Crear variable transformada (logaritmo de CPK)
hfcrd <- hfcrd %>%
  mutate(log_cpk = log(creatinine_phosphokinase + 1))  # Se suma 1 para evitar log(0)

## Convertir a data.frame para evitar imprimir los atributos de readr
hfcrd <- as.data.frame(hfcrd)

## Visualizar la estructura final del conjunto de datos
str(hfcrd)
## 'data.frame':    299 obs. of  14 variables:
##  $ age                     : num  75 55 65 50 65 90 75 60 65 80 ...
##  $ anaemia                 : Factor w/ 2 levels "No","Sí": 1 1 1 2 2 2 2 2 1 2 ...
##  $ creatinine_phosphokinase: num  582 7861 146 111 160 ...
##  $ diabetes                : Factor w/ 2 levels "No","Sí": 1 1 1 1 2 1 1 2 1 1 ...
##  $ ejection_fraction       : num  20 38 20 20 20 40 15 60 65 35 ...
##  $ high_blood_pressure     : Factor w/ 2 levels "No","Sí": 2 1 1 1 1 2 1 1 1 2 ...
##  $ platelets               : num  265000 263358 162000 210000 327000 ...
##  $ serum_creatinine        : num  1.9 1.1 1.3 1.9 2.7 2.1 1.2 1.1 1.5 9.4 ...
##  $ serum_sodium            : num  130 136 129 137 116 132 137 131 138 133 ...
##  $ sex                     : Factor w/ 2 levels "Hombre","Mujer": 1 1 1 1 2 1 1 1 2 1 ...
##  $ smoking                 : Factor w/ 2 levels "No","Sí": 1 1 2 1 1 2 1 2 1 2 ...
##  $ time                    : num  4 6 7 7 8 8 10 10 10 10 ...
##  $ DEATH_EVENT             : Factor w/ 2 levels "Censurado","Muerte": 2 2 2 2 2 2 2 2 2 2 ...
##  $ log_cpk                 : num  6.37 8.97 4.99 4.72 5.08 ...

3.6 Modelado

# Ajustamos un modelo de Cox proporcional de riesgos
library(survival)

cox_model <- coxph(Surv(time, DEATH_EVENT == "Muerte") ~ age + serum_creatinine + ejection_fraction,data = hfcrd)

summary(cox_model)
## Call:
## coxph(formula = Surv(time, DEATH_EVENT == "Muerte") ~ age + serum_creatinine + 
##     ejection_fraction, data = hfcrd)
## 
##   n= 299, number of events= 96 
## 
##                        coef exp(coef)  se(coef)      z Pr(>|z|)    
## age                0.044164  1.045154  0.008951  4.934 8.05e-07 ***
## serum_creatinine   0.357834  1.430228  0.068240  5.244 1.57e-07 ***
## ejection_fraction -0.048876  0.952300  0.010177 -4.803 1.57e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                   exp(coef) exp(-coef) lower .95 upper .95
## age                  1.0452     0.9568    1.0270    1.0637
## serum_creatinine     1.4302     0.6992    1.2512    1.6349
## ejection_fraction    0.9523     1.0501    0.9335    0.9715
## 
## Concordance= 0.717  (se = 0.03 )
## Likelihood ratio test= 66.51  on 3 df,   p=2e-14
## Wald test            = 71.68  on 3 df,   p=2e-15
## Score (logrank) test = 72.54  on 3 df,   p=1e-15

Este modelo usa como predictores:

  • age (edad):

    A mayor edad, mayor riesgo de mortalidad. Es una variable común en todos los modelos de supervivencia médica.

  • serum_creatinine (creatinina sérica):

    Altos niveles indican disfunción renal, que es un fuerte predictor de mal pronóstico en pacientes con insuficiencia cardíaca.

  • ejection_fraction (fracción de eyección):

    Mide qué porcentaje de sangre bombea el corazón con cada contracción. Valores bajos indican insuficiencia cardíaca grave y se asocian con mayor riesgo de muerte.

3.6.1 🔹 Coeficientes estimados

  • Indican el efecto logarítmico de cada predictor sobre el riesgo de muerte.

  • Si exponencias estos valores (exp(coef)), obtienes el hazard ratio (HR):

    • HR > 1 → mayor riesgo

    • HR < 1 → menor riesgo

    • HR = 1 → sin efecto

3.6.2 🔹 Error estándar, z y valor p

  • Sirven para evaluar si un predictor tiene un efecto estadísticamente significativo.

  • Un valor p < 0.05 suele indicar que el predictor es importante.

3.6.3 🔹 Estadísticas globales del modelo:

  • Concordance: mide cuán bien el modelo discrimina entre pacientes que mueren antes que otros. Un valor cercano a 1 es mejor (ideal: > 0.7).

  • Likelihood ratio test, Wald test, logrank test: pruebas estadísticas que evalúan si el modelo globalmente es significativo.

  • n / events: número de observaciones totales y cuántas muertes (eventos) ocurrieron.

3.6.4 📊 Gráfica 1

# Cargar librerías necesarias
library(survival)
library(survminer)
library(readr)

# Leer el dataset original
hfcrd <- read_excel("C:/Users/mario/OneDrive/Desktop/Estadística_aplicada/Proyecto Final/heart_failure_clinical_records_dataset.xlsx")

# Asegurarse de que las columnas estén correctas
hfcrd$DEATH_EVENT <- as.numeric(hfcrd$DEATH_EVENT)
hfcrd$sex <- factor(hfcrd$sex, levels = c(0, 1), labels = c("Mujer", "Hombre"))

# Crear objeto de supervivencia de Kaplan-Meier por sexo
fit_km <- survfit(Surv(time, DEATH_EVENT) ~ sex, data = hfcrd)

# Graficar curvas de supervivencia
ggsurvplot(
  fit_km,
  data = hfcrd,
  pval = TRUE,
  conf.int = TRUE,
  risk.table = TRUE,
  risk.table.height = 0.5,
  legend.title = "Sexo",
  legend.labs = c("Mujer", "Hombre"),
  palette = c("#D55E00", "#0072B2"),
  xlab = "Tiempo (días)",
  ylab = "P. de Supervivencia",
  break.time.by = 100,
  surv.median.line = "hv"
)

3.6.5 🔹 Análisis

La curva de supervivencia Kaplan-Meier muestra la probabilidad estimada de que los pacientes sobrevivan a lo largo del tiempo, estratificada por sexo.

-La línea naranja corresponde a las mujeres y la azul a los hombres.

-El valor p = 0.95 indica que no hay una diferencia estadísticamente significativa en las tasas de supervivencia entre hombres y mujeres en esta muestra.

-Las bandas sombreadas representan los intervalos de confianza al 95%, mostrando la incertidumbre en la estimación de la supervivencia.

-La tabla de riesgo debajo del gráfico muestra el número de pacientes que permanecen en riesgo (sin evento) en cada punto de tiempo clave.

En conclusión, el sexo no parece ser un factor discriminante importante en la supervivencia de pacientes con insuficiencia cardíaca dentro de este conjunto de datos. Este hallazgo puede orientar el enfoque clínico para considerar otros factores más relevantes en el pronóstico.

3.6.6 📊 Gráfica 2

library(survival)
library(survminer)
library(readr)
library(dplyr)

# Leer datos
hfcrd  <- read_excel("C:/Users/mario/OneDrive/Desktop/Estadística_aplicada/Proyecto Final/heart_failure_clinical_records_dataset.xlsx")

# Convertir variables necesarias
hfcrd$DEATH_EVENT <- as.numeric(hfcrd$DEATH_EVENT)

# Crear variable categórica para creatinina
hfcrd <- hfcrd %>%
  mutate(creatinina_cat = ifelse(serum_creatinine > 1.0, "Alta", "Normal"),
         creatinina_cat = factor(creatinina_cat, levels = c("Normal", "Alta")))

# Modelo de Kaplan-Meier por niveles de creatinina
fit_km_crea <- survfit(Surv(time, DEATH_EVENT) ~ creatinina_cat, data = hfcrd)

# Graficar curvas
ggsurvplot(
  fit_km_crea,
  data = hfcrd,
  pval = TRUE,
  conf.int = TRUE,
  risk.table = TRUE,
  risk.table.height = 0.3,
  legend.title = "Creatinina",
  legend.labs = c("Normal (≤1.0)", "Alta (>1.0)"),
  palette = c("#009E73", "#D55E00"),
  xlab = "Tiempo (días)",
  ylab = "Probabilidad de Supervivencia",
  break.time.by = 100,
  surv.median.line = "hv"
)

3.6.7 🔹 Análisis

En esta curva de Kaplan-Meier se evalúa el impacto de los niveles de creatinina sérica sobre la supervivencia:

  • Se observa que los pacientes con creatinina alta (>1.0 mg/dL) tienen una probabilidad de supervivencia notablemente menor en comparación con quienes presentan niveles normales.

  • El valor p reportado en el gráfico permite concluir que esta diferencia es estadísticamente significativa (si p < 0.05), lo cual refuerza la importancia clínica de esta variable como predictor de mortalidad.

  • Esto es consistente con el conocimiento médico: valores elevados de creatinina indican disfunción renal, un factor que suele agravar el pronóstico en insuficiencia cardíaca.

  • La tabla inferior muestra cómo el número de pacientes en riesgo disminuye con el tiempo en cada grupo, reforzando visualmente la divergencia en las curvas.

3.6.8 📊 Gráfica 3

library(survival)
library(survminer)
library(readr)
library(dplyr)

# Leer datos
hfcrd <- read_excel("C:/Users/mario/OneDrive/Desktop/Estadística_aplicada/Proyecto Final/heart_failure_clinical_records_dataset.xlsx")

# Convertir variables necesarias
hfcrd$DEATH_EVENT <- as.numeric(hfcrd$DEATH_EVENT)

# Clasificar fracción de eyección
hfcrd <- hfcrd %>%
  mutate(feyeccion_cat = ifelse(ejection_fraction <= 40, "Reducida", "Preservada"),
         feyeccion_cat = factor(feyeccion_cat, levels = c("Preservada", "Reducida")))

# Modelo de Kaplan-Meier por fracción de eyección
fit_km_ef <- survfit(Surv(time, DEATH_EVENT) ~ feyeccion_cat, data = hfcrd)

# Gráfico
ggsurvplot(
  fit_km_ef,
  data = hfcrd,
  pval = TRUE,
  conf.int = TRUE,
  risk.table = TRUE,
  risk.table.height = 0.4,
  legend.title = "Fracción eyección",
  legend.labs = c("Preservada (>40%)", "Reducida (≤40%)"),
  palette = c("#0072B2", "#D55E00"),
  xlab = "Tiempo (días)",
  ylab = "P. Supervivencia",
  break.time.by = 100,
  surv.median.line = "hv"
)

3.6.9 🔹 Análisis

Esta curva de supervivencia muestra cómo la fracción de eyección afecta la mortalidad:

  • Los pacientes con fracción de eyección reducida (≤ 40%) tienen menor probabilidad de supervivencia a lo largo del tiempo, como se refleja en la separación clara de las curvas.

  • El valor p generado por la prueba log-rank indica si esta diferencia es estadísticamente significativa. Un valor p < 0.05 confirmaría que la fracción de eyección es un predictor significativo de mortalidad.

  • Esto es coherente con la fisiopatología: una fracción de eyección baja implica peor función ventricular, lo que suele asociarse a mayor riesgo de eventos adversos y mortalidad.

  • La visualización con bandas de confianza y tabla de pacientes en riesgo aporta claridad y precisión al análisis.

3.7 Conclusiones

El análisis de supervivencia utilizando el modelo de riesgos proporcionales de Cox ha permitido identificar y cuantificar los factores clínicos más relevantes asociados con el riesgo de muerte en pacientes con insuficiencia cardíaca. A partir de los resultados obtenidos, se destacan las siguientes conclusiones clave:

  • Edad, creatinina sérica y fracción de eyección resultaron ser predictores estadísticamente significativos del riesgo de muerte en un momento dado. Esto no solo valida conocimientos médicos previos, sino que ofrece un marco cuantitativo para estimar el pronóstico de los pacientes.

  • La edad avanzada se asocia con un mayor riesgo, reflejando el deterioro fisiológico acumulado. A cada año adicional se incrementa el riesgo relativo de muerte, lo cual debe considerarse en decisiones clínicas.

  • Niveles elevados de creatinina sérica fueron fuertemente predictivos de mortalidad, indicando que la función renal deteriorada agrava significativamente el pronóstico. Esta variable mostró diferencias claras y significativas en las curvas de Kaplan-Meier.

  • Una fracción de eyección reducida (≤40%) se relaciona con menor supervivencia, confirmando su papel como indicador crítico del desempeño cardíaco. Esta categoría mostró una separación clara en las curvas de supervivencia y un valor p estadísticamente significativo.

  • Por el contrario, variables como el sexo no mostraron un efecto significativo sobre la supervivencia en este conjunto de datos, lo que sugiere que en este contexto clínico específico, el género no debería ser un criterio principal para estratificación del riesgo.

  • El modelo de Cox demostró buena capacidad de discriminación (índice de concordancia adecuado) y coherencia clínica, lo que lo convierte en una herramienta robusta para apoyar la toma de decisiones médicas y la planificación del seguimiento de pacientes.

Este estudio evidencia el poder del análisis de supervivencia para transformar datos clínicos en información accionable. Con un enfoque basado en datos, se pueden diseñar estrategias de tratamiento más personalizadas y eficientes, priorizando recursos hacia los pacientes con mayor riesgo y promoviendo una medicina verdaderamente preventiva y proactiva.

4 Referencias

[1] Consejo Privado de Competitividad & Universidad del Rosario. (2024). Índice de Competitividad de Ciudades 2024. https://compite.com.co/indice-de-competitividad-de-ciudades/

[2] James, G., Witten, D., Hastie, T., & Tibshirani, R. (2023). An introduction to statistical learning: With applications in R (2nd ed., corrected version). Springer. https://www.statlearning.com