1 DATA SET

Contexto: El dataset injury_data.csv contiene información sobre jugadores deportivos y factores asociados a la probabilidad de sufrir una lesión. Cada fila corresponde a un jugador

# cargamos los datos, nos aseguramos de cargarlo con dos opciones: una es cargando el archivo del data set en la carpeta en la que está guardado el script, o cargandolo directamente en environment > import dataset > from text(readr)
datos <- read.csv("injury_data.csv", header = TRUE, stringsAsFactors = FALSE)

# Convertimos la variable objetivo a factor (categórica binaria)
datos$Likelihood_of_Injury <- factor(datos$Likelihood_of_Injury,
                                      levels = c(0, 1),
                                      labels = c("No Lesión", "Lesión"))

# Vemos las primeras 6 filas del dataset
kable(head(datos, 6),
      caption = "Primeras 6 observaciones del dataset",
      digits = 2) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                full_width = FALSE, font_size = 12)
Primeras 6 observaciones del dataset
Player_Age Player_Weight Player_Height Previous_Injuries Training_Intensity Recovery_Time Likelihood_of_Injury
24 66.25 175.73 1 0.46 5 No Lesión
37 71.00 174.58 0 0.23 6 Lesión
32 80.09 186.33 0 0.61 2 Lesión
28 87.47 175.50 1 0.25 4 Lesión
25 84.66 190.18 0 0.58 1 Lesión
38 75.82 206.63 1 0.36 4 No Lesión

Variables que tenemos:

  • Player_Age: Edad del jugador (años).
  • Player_Weight: Peso del jugador (kg).
  • Player_Height: Altura del jugador (cm).
  • Previous_Injuries: Lesiones previas (0 = No, 1 = Sí).
  • Training_Intensity: Intensidad de entrenamiento [0–1].
  • Recovery_Time: Tiempo de recuperación (días).
  • Likelihood_of_Injury: Variable objetivo (0 = No lesión, 1 = Lesión).

2 Estadísticas Descriptivas

2.1 Resumen General

# estadistica descriptiva (univariada)
# media, mediana, desviación estándar,
# mínimo, máximo y cuartiles por variable


# Variables numéricas para el análisis
numericas <- c("Player_Age", "Player_Weight", "Player_Height",
               "Previous_Injuries", "Training_Intensity", "Recovery_Time")

resumen <- data.frame(
  Variable   = numericas,
  N          = sapply(datos[, numericas], function(x) sum(!is.na(x))),
  Media      = sapply(datos[, numericas], mean, na.rm = TRUE),
  Mediana    = sapply(datos[, numericas], median, na.rm = TRUE),
  DE         = sapply(datos[, numericas], sd, na.rm = TRUE),
  Min        = sapply(datos[, numericas], min, na.rm = TRUE),
  Max        = sapply(datos[, numericas], max, na.rm = TRUE),
  Q1         = sapply(datos[, numericas], quantile, 0.25, na.rm = TRUE),
  Q3         = sapply(datos[, numericas], quantile, 0.75, na.rm = TRUE)
)
rownames(resumen) <- NULL

kable(resumen,
      caption = "Estadistica descriptiva por Variable",
      digits = 3) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                full_width = FALSE, font_size = 12)
Estadistica descriptiva por Variable
Variable N Media Mediana DE Min Max Q1 Q3
Player_Age 1000 28.231 28.000 6.538 18.000 39.000 22.000 34.000
Player_Weight 1000 74.794 75.021 9.893 40.192 104.650 67.944 81.303
Player_Height 1000 179.751 180.034 9.889 145.286 207.309 173.037 186.558
Previous_Injuries 1000 0.515 1.000 0.500 0.000 1.000 0.000 1.000
Training_Intensity 1000 0.491 0.484 0.286 0.000 0.998 0.241 0.730
Recovery_Time 1000 3.466 4.000 1.701 1.000 6.000 2.000 5.000

Interpretación: Se observan escalas muy distintas entre variables ( Player_Height en cm vs Training_Intensity en [0,1]). Esto justifica la estandarización antes de hacer el analisis multi. ## Distribución por Variable: Histogramas

# histogramas
# usamos la función multi.hist del paquete psych para hacer los distintos histogramas que queremos por variable
# superpone la curva de densidad y la normal teórica


multi.hist(datos[, numericas], global = FALSE,
           main = "Histogramas con densidad por variable")

2.2 Boxplots: Detección Visual de Atípicos

#boxplot
# Los puntos fuera de los bigotes (±1.5*IQR)
# son candidatos a valores atípicos univariados
#vemos que las escalas son heterogéneas
# Player_Height domina visualmente por su escala (cm).
# Se requiere estandarización para análisis comparativos

boxplot(datos[, numericas],
        main  = "Boxplots - Escala Original",
        xlab  = "Variables",
        ylab  = "Valor",
        col   = "steelblue",
        frame = FALSE,
        las   = 2,
        cex.axis = 0.7)
grid()

2.3 Distribución de la Variable Objetivo

# Variable objetivo, Conteo y proporción de cada clase.
# Un desbalance severo (>80/20) puede afectar modelos
# supervisados y la interpretación de clusters
#
# Si las clases están balanceadas, los
# clusters obtenidos deberían reflejar ambas categorías
# de forma proporcional

tabla_target <- table(datos$Likelihood_of_Injury)
prop_target  <- round(prop.table(tabla_target) * 100, 1)

df_target <- data.frame(
  Clase      = names(tabla_target),
  Frecuencia = as.integer(tabla_target),
  Porcentaje = as.numeric(prop_target)
)

kable(df_target, caption = "Distribución de la Variable Objetivo") %>%
  kable_styling(bootstrap_options = c("striped", "hover"),
                full_width = FALSE, font_size = 12)
Distribución de la Variable Objetivo
Clase Frecuencia Porcentaje
No Lesión 500 50
Lesión 500 50 
ggplot(df_target, aes(x = Clase, y = Frecuencia, fill = Clase)) +
  geom_bar(stat = "identity", width = 0.5) +
  geom_text(aes(label = paste0(Porcentaje, "%")), vjust = -0.5, size = 4) +
  scale_fill_manual(values = c("steelblue", "tomato")) +
  labs(title = "Distribución de Likelihood_of_Injury",
       x = "Clase", y = "Frecuencia") +
  theme_minimal() +
  theme(legend.position = "none")

2.4 Matriz de Correlaciones

# mapa de calor de correlaciones (si recuerdan es como el de la vez pasada, basicamente entre mas color mas correlación hay entre las variables comparadas)
# Valores cercanos a ±1 indican
# relaciones lineales fuertes
#Las correlaciones altas entre variables
# sugieren redundancia de información, lo cual es relevante
# para PCA (que busca maximizar varianza explicada reduciendo
# dimensionalidad)

cor_matrix <- cor(datos[, numericas], use = "complete.obs")

corrplot(cor_matrix,
         method  = "color",
         type    = "upper",
         addCoef.col = "black",
         tl.col  = "black",
         tl.srt  = 45,
         number.cex = 0.7,
         title   = "Correlaciones entre Variables",
         mar     = c(0, 0, 1, 0))

3 Preprocesamiento

# Estandarizacion (z-score)
# Técnicamente: Transforma cada variable a media=0 y DE=1.
# Fórmula: z = (x - media) / DE
#
#Necesario porque las variables tienen escalas
# muy distintas. Sin estandarizar, variables como Player_Height
# (~170-210 cm) dominarían sobre Training_Intensity (~0-1)

datos.std <- scale(datos[, numericas])
rownames(datos.std) <- rownames(datos)


boxplot(datos.std,
        main    = "Boxplots - Variables Estandarizadas",
        xlab    = "Variables",
        ylab    = "Escala Z",
        col     = "lightgreen",
        frame   = FALSE,
        las     = 2,
        cex.axis = 0.7)
grid()

4 Detección de Observaciones Atípicas

4.1 Atípicos Univariados (Regla IQR)

# detección univariada de atípicos
# Técnicamente: Una observación es atípica si su valor
# supera Q3 + 1.5*IQR o es menor que Q1 - 1.5*IQR.
# Se itera sobre cada variable numérica estandarizada
#Los índices listados son observaciones que se desvían significativamente en al menos una variable

atipicos      <- list()
atipicos.unicos <- c()

for (variable in colnames(datos.std)) {
  vals <- boxplot(datos.std[, variable], plot = FALSE)$out
  idx  <- which(datos.std[, variable] %in% vals)
  atipicos[[variable]] <- idx
  atipicos.unicos <- union(atipicos.unicos, idx)
}

# Resumen por variable
resumen_atip <- data.frame(
  Variable      = names(atipicos),
  N_atipicos    = sapply(atipicos, length)
)

kable(resumen_atip,
      caption = "Atípicos Univariados por Variable",
      row.names = FALSE) %>%
  kable_styling(bootstrap_options = c("striped", "hover"),
                full_width = FALSE, font_size = 12)
Atípicos Univariados por Variable
Variable N_atipicos
Player_Age 0
Player_Weight 7
Player_Height 5
Previous_Injuries 0
Training_Intensity 0
Recovery_Time 0 
cat("Total de observaciones atípicas únicas (univariado):", length(atipicos.unicos), "\n")
## Total de observaciones atípicas únicas (univariado): 12

4.2 Atípicos Multivariados: Distancia de Mahalanobis

# Distancia de Mahalanobis
# La distancia de Mahalanobis considera la
# covarianza entre variables. Un punto es atípico multivariado
# si su distancia supera el cuantil 97.5% de una distribución
# Chi-cuadrado con p grados de libertad (donde p = núm. vars).
#
# Un punto con distancia de Mahalanobis alta
# es inusual considerando simultáneamente todas las variables

maha_dist <- mahalanobis(datos.std, colMeans(datos.std), cov(datos.std))

# Umbral Chi-cuadrado al 97.5%
umbral_maha <- qchisq(0.975, df = ncol(datos.std))

# Hist. de distancias
hist(maha_dist,
     freq    = FALSE,
     col     = "steelblue",
     border  = 0,
     main    = "Distancias de Mahalanobis",
     xlab    = "Distancia",
     ylab    = "Densidad")
curve(dchisq(x, df = ncol(datos.std)), add = TRUE, col = "red", lwd = 2)
abline(v = umbral_maha, col = "orange", lty = 2, lwd = 2)
legend("topright", legend = c("Chi-cuadrado teórico", "Umbral 97.5%"),
       col = c("red", "orange"), lty = c(1, 2), bty = "n")
grid()

# Identificamos atípicos multivariados
idx_maha <- which(maha_dist > umbral_maha)
atipicos[["Mahalanobis"]] <- idx_maha

cat("Atípicos multivariados detectados:", length(idx_maha), "\n")
## Atípicos multivariados detectados: 6
cat("Índices:", paste(idx_maha, collapse = ", "), "\n")
## Índices: 226, 281, 365, 435, 658, 961
# Univariados n Multivariados (intersección)
# Un atípico que aparece tanto en el análisis
# univariado como en el multivariado tiene mayor evidencia
# de ser "verdaderamente" atípico

# Estos son los casos que con mayor seguridad
# requieren revisión o tratamiento especial en el análisis

implicados <- intersect(atipicos.unicos, idx_maha)

cat("Observaciones atípicas en ambos métodos:", length(implicados), "\n")
## Observaciones atípicas en ambos métodos: 5
cat("Índices:", paste(implicados, collapse = ", "), "\n")
## Índices: 281, 435, 961, 365, 658
# Etiquetas para gráficos 
etiquetas <- rep("Típico", nrow(datos.std))
etiquetas[idx_maha] <- "Atípico"
etiquetas <- factor(etiquetas)

5 Clusterización

5.1 Selección del Método: Jerárquico Divisivo (DIANA)

# CLUSTERIZACIÓN JERÁRQUICA DIVISIVA (DIANA)
 
# DIANA (DIvisive ANAlysis) parte de un único
# cluster que contiene todos los datos y los va dividiendo
# recursivamente. El coeficiente de divisibilidad mide la
# cohesión: valores cercanos a 1 = clusters homogéneos
#
#  El dendrograma muestra cómo se forman los
# grupos. Las ramas largas indican separaciones claras entre
# clusters (alta heterogeneidad entre grupos)

# Matriz de disimilaridades Euclidiana sobre datos estandarizados
D <- daisy(datos.std, metric = "euclidean", stand = FALSE)

# Clusterizacion
datos.diana <- diana(x = D, diss = TRUE, stand = FALSE)

cat("Coeficiente de Divisibilidad DIANA:", round(datos.diana$dc, 4), "\n")
## Coeficiente de Divisibilidad DIANA: 0.8817
cat("(Cercano a 1 = clusters bien definidos)\n")
## (Cercano a 1 = clusters bien definidos)
# Dendrograma
fviz_dend(datos.diana,
          k          = 3,
          k_colors   = c("steelblue", "tomato", "forestgreen"),
          rect       = TRUE,
          rect_fill  = TRUE,
          main       = "Dendrograma DIANA - Injury Data",
          ylab       = "Altura",
          cex        = 0.3)

5.2 Número Óptimo de Clusters: Índices Dunn y Davies-Bouldin

# num. óptimo de clusters
# se evalúan soluciones de K=2 a K=7 usando
# dos índices internos de validación:
#
#   - indice DUNN: Rango [0, ∞]. MAYOR = mejor.
#   - indice Davies-Bouldin (DB): Rango [0, ∞]. MENOR = mejor
#
#El K óptimo es aquel donde DUNN es máximo
# y DB es mínimo

dunn_vals <- numeric()
db_vals   <- numeric()

for (k in 2:7) {
  clusters_k <- cutree(datos.diana, k = k)
  dunn_vals  <- c(dunn_vals,
                  as.numeric(intCriteria(datos.std, clusters_k, "Dunn")))
  db_vals    <- c(db_vals,
                  as.numeric(intCriteria(datos.std, clusters_k, "Davies_Bouldin")))
}

names(dunn_vals) <- names(db_vals) <- 2:7

# resumen
df_indices <- data.frame(
  K              = 2:7,
  Dunn           = round(dunn_vals, 4),
  Davies_Bouldin = round(db_vals, 4)
)

kable(df_indices,
      caption = "Índices de Validación por Número de Clusters",
      row.names = FALSE) %>%
  kable_styling(bootstrap_options = c("striped", "hover"),
                full_width = FALSE, font_size = 12)
Índices de Validación por Número de Clusters
K Dunn Davies_Bouldin
2 0.0560 2.1906
3 0.0538 2.0867
4 0.0551 1.9684
5 0.0559 2.0878
6 0.0565 1.9946
7 0.0607 1.8834 
# Graficas 
par(mfrow = c(1, 2))

plot(2:7, dunn_vals,
     type  = "b",
     pch   = 16,
     col   = "steelblue",
     bty   = "n",
     main  = "Índice DUNN",
     xlab  = "Número de Clusters",
     ylab  = "Dunn")
abline(v = which.max(dunn_vals) + 1, col = "red", lty = 2)
grid()

plot(2:7, db_vals,
     type  = "b",
     pch   = 16,
     col   = "tomato",
     bty   = "n",
     main  = "Davies-Bouldin",
     xlab  = "Número de Clusters",
     ylab  = "DB")
abline(v = which.min(db_vals) + 1, col = "red", lty = 2)
grid()

par(mfrow = c(1, 1))

k_optimo_dunn <- which.max(dunn_vals) + 1
k_optimo_db   <- which.min(db_vals) + 1
cat("K óptimo según DUNN:", k_optimo_dunn, "\n")
## K óptimo según DUNN: 7
cat("K óptimo según Davies-Bouldin:", k_optimo_db, "\n")
## K óptimo según Davies-Bouldin: 7
# Metodo del silhouette (validación adicional)

# El coeficiente de silueta mide qué tan similar
# es cada observación a su propio cluster vs el vecino más
# cercano. Rango [-1, 1]: valores altos = buena asignación
#
#Un pico en el promedio de silueta indica el
# número ideal de grupos para este dataset

fviz_nbclust(datos.std,
             FUN      = hcut,
             method   = "silhouette",
             k.max    = 7,
             hc_method = "ward.D2") +
  labs(title    = "Método Silhouette - Número Óptimo de Clusters",
       subtitle = "Máximo silhouette = K óptimo") +
  theme_minimal()

5.3 Clusterización Final con K Óptimo

# DENDROGRAMA FINAL CON K ÓPTIMO

# Se corta el árbol jerárquico de DIANA en el
# nivel que produce el número óptimo de clusters determinado
# por los índices de validación

# el color de cada rama identifica el cluster
# al que pertenece cada observación

k_final <- k_optimo_dunn  # Usamos el mejor según DUNN

clusters_finales <- cutree(datos.diana, k = k_final)

fviz_dend(datos.diana,
          k         = k_final,
          rect      = TRUE,
          rect_fill = TRUE,
          main      = paste("Dendrograma DIANA - K =", k_final),
          cex       = 0.3,
          ylab      = "Altura")

cat("Distribución de observaciones por cluster:\n")
## Distribución de observaciones por cluster:
print(table(clusters_finales))
## clusters_finales
##   1   2   3   4   5   6   7 
## 306 237 127  57 110 114  49

5.4 Identificación de Atípicos mediante Clusterización

#aitpicos por clusterizaicon
# En métodos jerárquicos, los clusters con
# muy pocas observaciones (< 2% del total) suelen contener
# outliers
#
# Los jugadores en clusters pequeños tienen
# perfiles únicos que los separan del resto

tabla_clusters <- table(clusters_finales)
n_total        <- nrow(datos.std)
umbral_cluster <- max(3, floor(0.02 * n_total))  # 2% del total o mínimo 3

cat("Tamaño de cada cluster:\n")
## Tamaño de cada cluster:
print(tabla_clusters)
## clusters_finales
##   1   2   3   4   5   6   7 
## 306 237 127  57 110 114  49
cat("\nUmbral para considerar cluster 'pequeño':", umbral_cluster, "observaciones\n")
## 
## Umbral para considerar cluster 'pequeño': 20 observaciones
# Identificar posibles atípicos
clusters_pequenos <- names(tabla_clusters[tabla_clusters <= umbral_cluster])
obs_atipicas_cluster <- which(clusters_finales %in% as.integer(clusters_pequenos))

cat("\nClusters pequeños identificados:", paste(clusters_pequenos, collapse = ", "), "\n")
## 
## Clusters pequeños identificados:
cat("Observaciones en clusters pequeños:", length(obs_atipicas_cluster), "\n")
## Observaciones en clusters pequeños: 0
# Gráfica
barplot(tabla_clusters,
        col  = rainbow(k_final),
        main = "Tamaño de Clusters",
        xlab = "Cluster",
        ylab = "Número de Observaciones",
        border = NA)
abline(h = umbral_cluster, col = "red", lty = 2, lwd = 2)
legend("topright", legend = "Umbral atípico", col = "red", lty = 2, bty = "n")
grid()

6 Escalamiento Multidimensional (MDS)

6.1 MDS Métrico (cmdscale)

# escalamiento multidimensional metrico (cmsdscale)

# MDS busca una representación en baja
# dimensión (2D o 3D) que preserve las distancias originales
# entre observaciones
# Observaciones cercanas en el plano MDS
# tienen perfiles similares. Los clusters DIANA se colorean
# para ver si MDS los separa visualmente

datos.centrados <- scale(datos[, numericas], center = TRUE, scale = FALSE)
distancias      <- dist(datos.centrados)

# MDS a 2 dimensiones
mds_coords <- cmdscale(distancias, k = 2, eig = TRUE)

# Proporción de varianza explicada por cada dimensión
varianza_mds <- mds_coords$eig / sum(abs(mds_coords$eig))
cat("Varianza explicada por Dim1:", round(varianza_mds[1] * 100, 1), "%\n")
## Varianza explicada por Dim1: 41.7 %
cat("Varianza explicada por Dim2:", round(varianza_mds[2] * 100, 1), "%\n")
## Varianza explicada por Dim2: 39.4 %
cat("Total explicado (2D):", round(sum(varianza_mds[1:2]) * 100, 1), "%\n")
## Total explicado (2D): 81.1 %
# Gráfico MDS coloreado por cluster
df_mds <- data.frame(
  Dim1    = mds_coords$points[, 1],
  Dim2    = mds_coords$points[, 2],
  Cluster = factor(clusters_finales),
  Lesion  = datos$Likelihood_of_Injury
)

ggplot(df_mds, aes(x = Dim1, y = Dim2, color = Cluster, shape = Lesion)) +
  geom_point(alpha = 0.7, size = 2) +
  scale_color_brewer(palette = "Set1") +
  labs(title    = "MDS Métrico - Injury Data",
       subtitle = paste0("Dim1: ", round(varianza_mds[1]*100,1),
                         "% | Dim2: ", round(varianza_mds[2]*100,1), "%"),
       x = "Dimensión 1", y = "Dimensión 2",
       color = "Cluster", shape = "Lesión") +
  theme_minimal() +
  theme(legend.position = "right")

6.2 MDS No Métrico (isoMDS)

# ============================================================
# Escalamiento multidimensioanl no metrico (isoMDS)

# A diferencia del MDS métrico, el no métrico
# solo preserva el ORDEN de las distancias (ranking).
# El STRESS mide el ajuste: stress < 0.10 = bueno,
# stress < 0.05 = excelente
#
# Útil cuando las distancias no son
# perfectamente euclidianas

# Se trabaja con una muestra para eficiencia computacional
set.seed(42)
n_muestra  <- min(300, nrow(datos.std))
idx_muestra <- sample(1:nrow(datos.std), n_muestra)
datos_muestra <- datos.std[idx_muestra, ]

dist_muestra  <- dist(datos_muestra)
fit_iso <- isoMDS(dist_muestra, k = 2)
## initial  value 35.464899 
## iter   5 value 31.064631
## final  value 30.174144 
## converged
cat("Stress del MDS No Métrico:", round(fit_iso$stress, 2), "%\n")
## Stress del MDS No Métrico: 30.17 %
cat("(< 10% = bueno, < 5% = excelente)\n")
## (< 10% = bueno, < 5% = excelente)
# Gráfico
df_iso <- data.frame(
  x       = fit_iso$points[, 1],
  y       = fit_iso$points[, 2],
  Cluster = factor(clusters_finales[idx_muestra]),
  Lesion  = datos$Likelihood_of_Injury[idx_muestra]
)

ggplot(df_iso, aes(x = x, y = y, color = Cluster, shape = Lesion)) +
  geom_point(alpha = 0.7, size = 2) +
  scale_color_brewer(palette = "Set2") +
  labs(title    = paste0("MDS No Métrico (isoMDS) - n = ", n_muestra),
       subtitle = paste0("Stress = ", round(fit_iso$stress, 2), "%"),
       x = "Dimensión 1", y = "Dimensión 2",
       color = "Cluster", shape = "Lesión") +
  theme_minimal()

7 Componentes Principales (PCA)

7.1 Cálculo y Varianza Explicada

#Ánalisis de componentes principales PCA

# PCA descompone la matriz de covarianzas en
# vectores propios (direcciones de máxima varianza) y valores
# propios (cantidad de varianza en esa dirección)
# Se usa princomp() con datos estandarizados y cor=FALSE
# (ya estan en misma escala)
#
# Se busca retener suficientes componentes
# para explicar al menos el 80% de la varianza total

datos.pca <- princomp(datos.std, cor = FALSE)
resumen_pca <- summary(datos.pca)

# Varianza acumulada
varianza_acum <- cumsum(datos.pca$sdev^2) / sum(datos.pca$sdev^2)

df_pca_var <- data.frame(
  Componente        = paste0("PC", 1:length(datos.pca$sdev)),
  Varianza_Prop     = round(datos.pca$sdev^2 / sum(datos.pca$sdev^2) * 100, 2),
  Varianza_Acum     = round(varianza_acum * 100, 2)
)

kable(df_pca_var,
      caption = "Varianza Explicada por Componente Principal",
      row.names = FALSE) %>%
  kable_styling(bootstrap_options = c("striped", "hover"),
                full_width = FALSE, font_size = 12)
Varianza Explicada por Componente Principal
Componente Varianza_Prop Varianza_Acum
PC1 18.62 18.62
PC2 18.09 36.71
PC3 17.14 53.85
PC4 16.10 69.95
PC5 15.44 85.39
PC6 14.61 100.00 
# Número de componentes para 80% de varianza
n_comp_80 <- which(varianza_acum >= 0.80)[1]
cat("\nComponentes necesarios para explicar >= 80% de varianza:", n_comp_80, "\n")
## 
## Componentes necesarios para explicar >= 80% de varianza: 5

7.2 Gráfica de Sedimentación (Scree Plot)

# Scree plot
#Grafica el valor propio (eigenvalue) de cada
# componente en orden descendente. El codo de la
# curva señala el número óptimo de componentes a retener
#
# La línea roja discontinua marca el umbral del 80%

eigenvalues <- datos.pca$sdev^2

par(mfrow = c(1, 2))

# Scree plot
plot(eigenvalues,
     type  = "b",
     pch   = 16,
     col   = "steelblue",
     bty   = "n",
     main  = "Scree Plot",
     xlab  = "Componente",
     ylab  = "Valor Propio")
abline(h = 1, col = "red", lty = 2)
grid()

# Varianza acumulada
plot(varianza_acum * 100,
     type  = "b",
     pch   = 16,
     col   = "steelblue",
     bty   = "n",
     main  = "Varianza Acumulada",
     xlab  = "Número de Componentes",
     ylab  = "% Varianza Acumulada",
     ylim  = c(0, 105))
abline(h = 80, col = "red", lty = 2, lwd = 2)
abline(v = n_comp_80, col = "orange", lty = 2, lwd = 2)
text(n_comp_80 + 0.2, 75,
     labels = paste("PC", n_comp_80, "\n>=80%"),
     col = "orange", cex = 0.8, adj = 0)
grid()

par(mfrow = c(1, 1))

7.3 Biplot y Cargas (Loadings)

# BIPLOT PCA
#El biplot representa simultáneamente las
# observaciones (puntos) y las variables originales (flechas)
# en el espacio de las primeras dos componentes

#  Flechas en la misma dirección = variables correlacionadas positivamente.
#  Flechas opuestas = correlación negativa
#  Flechas perpendiculares = variables independientes

biplot(datos.pca,
       main   = "Biplot - PCA Injury Data",
       col    = c("gray50", "steelblue"),
       cex    = 0.5,
       scale  = 0,
       xlab   = paste0("PC1 (", round(df_pca_var$Varianza_Prop[1], 1), "%)"),
       ylab   = paste0("PC2 (", round(df_pca_var$Varianza_Prop[2], 1), "%)"))
grid()

# cargas de las componentes principales
# Las cargas muestran la correlación entre
# cada variable original y cada componente principal.
# Valores absolutos altos (> 0.3) indican contribución
# significativa a la componente
#
# Permite "nombrar" cada componente principal
# en términos de las variables que la definen

loadings_df <- as.data.frame(unclass(datos.pca$loadings)[, 1:n_comp_80])
colnames(loadings_df) <- paste0("PC", 1:n_comp_80)
rownames(loadings_df) <- numericas

kable(round(loadings_df, 3),
      caption = paste0("Cargas (Loadings) de las Primeras ", n_comp_80, " Componentes")) %>%
  kable_styling(bootstrap_options = c("striped", "hover"),
                full_width = FALSE, font_size = 12)
Cargas (Loadings) de las Primeras 5 Componentes
PC1 PC2 PC3 PC4 PC5
Player_Age 0.405 0.528 0.117 0.510 0.110
Player_Weight 0.374 -0.595 -0.053 -0.274 0.416
Player_Height 0.387 0.115 -0.715 0.118 0.339
Previous_Injuries -0.108 -0.535 -0.283 0.679 -0.387
Training_Intensity 0.600 0.025 -0.070 -0.328 -0.726
Recovery_Time -0.419 0.258 -0.622 -0.288 -0.148

7.4 Visualización PCA con Clusters

# Scatter plot PC1 vs PC2 coloreado por cluster

# Se proyectan las observaciones al espacio
# de los dos primeros componentes y se colorean según el
# cluster asignado por DIANA
#
# Si los clusters DIANA se separan bien en
# el plano PCA, esto sugiere que PCA y la clusterización
# capturan la misma estructura en los datos

scores_pca <- as.data.frame(datos.pca$scores[, 1:2])
colnames(scores_pca) <- c("PC1", "PC2")
scores_pca$Cluster <- factor(clusters_finales)
scores_pca$Lesion  <- datos$Likelihood_of_Injury
scores_pca$Atipico <- ifelse(1:nrow(datos.std) %in% idx_maha, "Atípico", "Típico")

ggplot(scores_pca, aes(x = PC1, y = PC2, color = Cluster,
                        shape = Atipico, size = Atipico)) +
  geom_point(alpha = 0.7) +
  scale_color_brewer(palette = "Set1") +
  scale_shape_manual(values = c("Atípico" = 17, "Típico" = 16)) +
  scale_size_manual(values  = c("Atípico" = 3, "Típico" = 1.5)) +
  labs(title    = "PCA: PC1 vs PC2 (Clusters + Atípicos)",
       subtitle = paste0("PC1: ", round(df_pca_var$Varianza_Prop[1], 1),
                         "% | PC2: ", round(df_pca_var$Varianza_Prop[2], 1), "%"),
       x = paste0("PC1 (", round(df_pca_var$Varianza_Prop[1], 1), "%)"),
       y = paste0("PC2 (", round(df_pca_var$Varianza_Prop[2], 1), "%)"),
       color = "Cluster", shape = "Tipo obs.", size = "Tipo obs.") +
  theme_minimal()

8 Comparación MDS vs PCA

# Comparación: MDS MÉTRICO vs PCA (Primeras 2 Dimensiones)

# Tanto MDS como PCA producen representaciones
# en baja dimensión. La equivalencia teórica es:
#   - PCA en datos centrados = MDS con distancias Euclidianas
#
#  Si ambas representaciones son similares,
# la estructura encontrada es robusta

# equivalencia con MDS
pca_centrado <- princomp(datos.centrados, cor = FALSE)

# Comparación visual
par(mfrow = c(1, 2))

plot(mds_coords$points[, 1], mds_coords$points[, 2],
     col  = rainbow(k_final)[clusters_finales],
     pch  = 16,
     cex  = 0.5,
     bty  = "n",
     main = "MDS Métrico",
     xlab = paste0("Dim1 (", round(varianza_mds[1]*100,1), "%)"),
     ylab = paste0("Dim2 (", round(varianza_mds[2]*100,1), "%)"))
grid()

plot(pca_centrado$scores[, 1], pca_centrado$scores[, 2],
     col  = rainbow(k_final)[clusters_finales],
     pch  = 16,
     cex  = 0.5,
     bty  = "n",
     main = "PCA (datos centrados)",
     xlab = paste0("PC1 (", round(summary(pca_centrado)$sdev[1]^2/
                     sum(summary(pca_centrado)$sdev^2)*100, 1), "%)"),
     ylab = paste0("PC2 (", round(summary(pca_centrado)$sdev[2]^2/
                     sum(summary(pca_centrado)$sdev^2)*100, 1), "%)"))
grid()

par(mfrow = c(1, 1))
# Tabla comparativa mds vs pca
#La tabla resume las características
# fundamentales de cada enfoque

df_comp <- data.frame(
  Criterio         = c("Objetivo", "Input", "Output",
                       "Equivalencia", "Escala requerida",
                       "Varianza explicada (2D)",
                       "Ventaja principal"),
  MDS_Metrico      = c("Preservar distancias",
                       "Matriz de distancias",
                       "Coordenadas en baja dim.",
                       "Con PCA (dist. Euclidiana)",
                       "No requerida",
                       paste0(round(sum(varianza_mds[1:2]) * 100, 1), "%"),
                       "Funciona con cualquier distancia"),
  PCA              = c("Maximizar varianza",
                       "Matriz de datos",
                       "Componentes principales",
                       "Con MDS (datos centrados)",
                       "Estandarización recomendada",
                       paste0(round(sum(df_pca_var$Varianza_Prop[1:2]), 1), "%"),
                       "Interpretación de loadings")
)

kable(df_comp,
      caption = "Comparación Técnica: MDS vs PCA",
      col.names = c("Criterio", "MDS Métrico", "PCA")) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                full_width = TRUE, font_size = 12)
Comparación Técnica: MDS vs PCA
Criterio MDS Métrico PCA
Objetivo Preservar distancias Maximizar varianza
Input Matriz de distancias Matriz de datos
Output Coordenadas en baja dim. Componentes principales
Equivalencia Con PCA (dist. Euclidiana) Con MDS (datos centrados)
Escala requerida No requerida Estandarización recomendada
Varianza explicada (2D) 81.1% 36.7%
Ventaja principal Funciona con cualquier distancia Interpretación de loadings

9 Feature Importance con Random Forest

9.1 Entrenamiento del Modelo

# FEATURE IMPORTANCE CON RANDOM FOREST

#Random Forest es un ensamble de árboles de
# decisión. La importancia de cada variable se mide por:
#   . Mean Decrease in Gini (MDG): reducción promedio del
#      índice de impureza Gini al dividir por esa variable
#   . Mean Decrease in Accuracy (MDA): reducción del accuracy
#      cuando se permuta aleatoriamente esa variable
# Mayor MDG/MDA = variable más importante para predecir
#
# Variables con alta importancia en RF son
# las que más contribuyen a distinguir jugadores con/sin lesión

set.seed(42)

# Variable objetivo 
y_rf <- datos$Likelihood_of_Injury

# Entrenamiento del Random Forest
rf_model <- randomForest(
  x          = datos[, numericas],
  y          = y_rf,
  ntree      = 500,
  importance = TRUE,
  random.seed = 42
)

cat("Random Forest - OOB Error Rate:", 
    round(rf_model$err.rate[500, "OOB"] * 100, 2), "%\n")
## Random Forest - OOB Error Rate: 50.3 %
cat("(Out-of-Bag: estimación sin necesidad de conjunto de prueba)\n\n")
## (Out-of-Bag: estimación sin necesidad de conjunto de prueba)
importance_df <- as.data.frame(importance(rf_model))
importance_df$Variable <- rownames(importance_df)

# Renombrar columnas según lo que devuelva importance()
colnames(importance_df)[1:2] <- c("MDA", "MDG")
importance_df <- importance_df[order(-importance_df$MDG), ]

tabla_imp <- importance_df[, c("Variable", "MDA", "MDG")]
tabla_imp$MDA <- round(tabla_imp$MDA, 4)
tabla_imp$MDG <- round(tabla_imp$MDG, 4)

kable(tabla_imp,
      caption = "Feature Importance - Random Forest",
      row.names = FALSE) %>%
  kable_styling(bootstrap_options = c("striped", "hover"),
                full_width = FALSE, font_size = 12)
Feature Importance - Random Forest
Variable MDA MDG
Previous_Injuries -1.7241 4.2598
Recovery_Time 3.1825 -0.2438
Player_Age -0.4875 -0.3096
Player_Weight 1.3844 -0.5173
Training_Intensity -0.1112 -0.7322
Player_Height -1.5623 -4.1006

9.2 Gráfica de Importancia

#feature importance visualiacion
# se grafican ambas métricas de importancia
# (MDA y MDG) en graficas de barras ordenadas
#
# 
# - MDA (Accuracy): más relevante para predicción
# - MDG (Gini): más relevante para separación de clases

importance_df_plot <- importance_df[order(-importance_df$MDA), ]

p1 <- ggplot(importance_df_plot,
             aes(x = reorder(Variable, MDA), y = MDA)) +
  geom_bar(stat = "identity", fill = "steelblue") +
  coord_flip() +
  labs(title = "Feature Importance\n(Mean Decrease Accuracy)",
       x = NULL, y = "MDA") +
  theme_minimal()

p2 <- ggplot(importance_df_plot,
             aes(x = reorder(Variable, MDG), y = MDG)) +
  geom_bar(stat = "identity", fill = "tomato") +
  coord_flip() +
  labs(title = "Feature Importance\n(Mean Decrease Gini)",
       x = NULL, y = "MDG") +
  theme_minimal()

grid.arrange(p1, p2, ncol = 2)

9.3 Comparación: PCA Loadings vs Feature Importance

# 
# COMPARACIÓN: PCA (varianza) vs FEATURE IMPORTANCE (predicción)
# 
# - PCA Loadings: contribución de cada variable a la primera
#   componente (estructura de varianza)
# - RF Importance: contribución de cada variable a predecir
#   la probabilidad de lesion (utilidad predictiva)
#
# 
# - Alta carga PCA Y alta importancia RF = variable informativa y predictiva
# - Alta carga PCA pero baja importancia RF = varianza sin predictibilidad
# - Baja carga PCA pero alta importancia RF = variable discriminante pese a poca varianza

# Carga absoluta en PC1 (primera componente de mayor varianza)
loadings_pc1 <- abs(datos.pca$loadings[, 1])

# Normalizar ambas métricas al rango [0,1]
norm01 <- function(x) (x - min(x)) / (max(x) - min(x))

comp_df <- data.frame(
  Variable        = numericas,
  PCA_Loading_PC1 = norm01(loadings_pc1[numericas]),
  RF_Importance   = norm01(importance_df[match(numericas, importance_df$Variable), "MDG"])
)
comp_df <- comp_df[order(-comp_df$RF_Importance), ]

# comparativa
comp_long <- reshape(comp_df,
                     varying   = c("PCA_Loading_PC1", "RF_Importance"),
                     v.names   = "Valor",
                     timevar   = "Metodo",
                     times     = c("PCA Loading |PC1|", "RF Importance (MDG)"),
                     direction = "long")

ggplot(comp_long, aes(x = reorder(Variable, Valor),
                      y = Valor, fill = Metodo)) +
  geom_bar(stat = "identity", position = "dodge") +
  coord_flip() +
  scale_fill_manual(values = c("steelblue", "tomato")) +
  labs(title    = "Comparación: PCA Loadings vs Feature Importance (RF)",
       subtitle = "Valores normalizados al rango [0,1] para comparabilidad",
       x = NULL, y = "Importancia Normalizada",
       fill = "Método") +
  theme_minimal() +
  theme(legend.position = "bottom")

# pca vs feautre importance
# La tabla consolida el ranking de cada
# variable según ambos enfoques

comp_df$Rank_PCA <- rank(-comp_df$PCA_Loading_PC1)
comp_df$Rank_RF  <- rank(-comp_df$RF_Importance)
comp_df$Diferencia_Ranking <- abs(comp_df$Rank_PCA - comp_df$Rank_RF)

kable(comp_df[order(comp_df$Rank_RF), ],
      caption  = "Comparación de Rankings: PCA vs Random Forest",
      row.names = FALSE, digits = 3) %>%
  kable_styling(bootstrap_options = c("striped", "hover"),
                full_width = FALSE, font_size = 12)
Comparación de Rankings: PCA vs Random Forest
Variable PCA_Loading_PC1 RF_Importance Rank_PCA Rank_RF Diferencia_Ranking
Previous_Injuries 0.000 1.000 6 1 5
Recovery_Time 0.632 0.461 2 2 0
Player_Age 0.603 0.453 3 3 0
Player_Weight 0.541 0.429 5 4 1
Training_Intensity 1.000 0.403 1 5 4
Player_Height 0.567 0.000 4 6 2

10 Resumen globak

# resuen
cat("=== RESUMEN DEL ANÁLISIS ===\n\n")
## === RESUMEN DEL ANÁLISIS ===
cat("1. DATOS:\n")
## 1. DATOS:
cat("   -", nrow(datos), "jugadores,", length(numericas), "variables predictoras\n")
##    - 1000 jugadores, 6 variables predictoras
cat("   - Variable objetivo: Likelihood_of_Injury\n")
##    - Variable objetivo: Likelihood_of_Injury
cat("   - Distribución:", table(datos$Likelihood_of_Injury)[1], "sin lesión /",
    table(datos$Likelihood_of_Injury)[2], "con lesión\n\n")
##    - Distribución: 500 sin lesión / 500 con lesión
cat("2. ATÍPICOS:\n")
## 2. ATÍPICOS:
cat("   - Univariados (IQR):", length(atipicos.unicos), "observaciones\n")
##    - Univariados (IQR): 12 observaciones
cat("   - Multivariados (Mahalanobis):", length(idx_maha), "observaciones\n")
##    - Multivariados (Mahalanobis): 6 observaciones
cat("   - Ambos métodos (intersección):", length(implicados), "observaciones\n\n")
##    - Ambos métodos (intersección): 5 observaciones
cat("3. CLUSTERIZACIÓN (DIANA):\n")
## 3. CLUSTERIZACIÓN (DIANA):
cat("   - K óptimo (DUNN):", k_optimo_dunn, "clusters\n")
##    - K óptimo (DUNN): 7 clusters
cat("   - K óptimo (Davies-Bouldin):", k_optimo_db, "clusters\n")
##    - K óptimo (Davies-Bouldin): 7 clusters
cat("   - Coef. divisibilidad:", round(datos.diana$dc, 3), "\n\n")
##    - Coef. divisibilidad: 0.882
cat("4. MDS:\n")
## 4. MDS:
cat("   - Varianza explicada (MDS métrico 2D):",
    round(sum(varianza_mds[1:2]) * 100, 1), "%\n")
##    - Varianza explicada (MDS métrico 2D): 81.1 %
cat("   - Stress MDS no métrico:", round(fit_iso$stress, 2), "%\n\n")
##    - Stress MDS no métrico: 30.17 %
cat("5. PCA:\n")
## 5. PCA:
cat("   - Componentes para >=80% varianza:", n_comp_80, "\n")
##    - Componentes para >=80% varianza: 5
cat("   - Varianza acumulada con", n_comp_80, "CPs:",
    round(varianza_acum[n_comp_80] * 100, 1), "%\n\n")
##    - Varianza acumulada con 5 CPs: 85.4 %
cat("6. FEATURE IMPORTANCE (RF):\n")
## 6. FEATURE IMPORTANCE (RF):
cat("   - Variable más importante (MDG):",
    importance_df$Variable[1], "\n")
##    - Variable más importante (MDG): Previous_Injuries
cat("   - OOB Error Rate:",
    round(rf_model$err.rate[500, "OOB"] * 100, 2), "%\n")
##    - OOB Error Rate: 50.3 %

10.1 Conclusiones

Estadísticas Descriptivas: Las variables presentan escalas heterogéneas, lo que requiere estandarización obligatoria. No se detectaron porcentajes elevados de atípicos extremos, pero la distancia de Mahalanobis reveló casos multivariados que conviene monitorear.

Clusterización (DIANA): El método jerárquico divisivo permitió identificar grupos de jugadores con perfiles de riesgo diferenciados. Los índices Dunn y Davies-Bouldin confirmaron el número óptimo de clusters, y los clusters pequeños concentran los perfiles más atípicos.

MDS: La representación bidimensional captura una porción razonable de la varianza original. La concordancia entre MDS métrico y PCA (datos centrados) confirma la equivalencia teórica entre ambos métodos.

PCA: Con los primeros componentes se alcanza el umbral del 80% de varianza. Los loadings permiten identificar qué variables dominan cada dimensión latente.

Feature Importance vs PCA: La comparación revela que algunas variables con alta varianza (alta carga en PCA) no necesariamente son las más predictivas de lesión (RF importance), y viceversa. Esto ilustra la diferencia fundamental entre análisis exploratorio (no supervisado) y predictivo (supervisado).