Taller 1 - Análisis Avanzado de Datos
Sifred Mendoza, Juan Camilo Sanmiguel y Alejandro Rubiano
2026-03-11
1 Carga de librerías y datos
library(glmnet)
library(ggplot2)
library(corrplot)
# Carga del archivo (separado por comas)
datos <- read.table("taller1.txt", header = TRUE, sep = ",", dec = ".")
cat("Dimensiones del dataset:", dim(datos), "\n")## Dimensiones del dataset: 1200 5001# Esperado: 1200 5001 (1200 líneas celulares, 5000 genes + 1 variable respuesta)1.1 Análisis descriptivo inicial
Antes de responder las preguntas, realizamos una exploración inicial de los datos para entender su estructura.
# Estructura general
str(datos[, 1:6]) # Primeras 6 columnas para no saturar la salida## 'data.frame': 1200 obs. of 6 variables:
## $ y : num -2.24 -3.53 2.99 -1.42 6.61 ...
## $ V1: num -0.5605 -0.2302 1.5587 0.0705 0.1293 ...
## $ V2: num 0.62 -0.758 0.852 -0.748 0.63 ...
## $ V3: num -0.483 -0.531 -0.588 -0.412 0.709 ...
## $ V4: num -0.214 1.198 0.232 -0.503 0.63 ...
## $ V5: num -0.6327 0.1092 -1.5626 -0.0402 -0.0363 ...# Descriptivos de la variable respuesta Y
cat("\nResumen de la variable respuesta (Y = efectividad del tratamiento):\n")##
## Resumen de la variable respuesta (Y = efectividad del tratamiento):summary(datos$y)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -13.9094 -3.0931 -0.0759 -0.1668 2.7338 16.8420# Histograma de Y
hist(datos$y, breaks = 30,
main = "Distribución de la efectividad del tratamiento (Y)",
xlab = "Efectividad", col = "steelblue", border = "white")
# Separar variable respuesta y predictores
y <- datos$y
X <- as.matrix(datos[, -1])
cat("Dimensiones de X (predictores):", dim(X), "\n")## Dimensiones de X (predictores): 1200 5000cat("Longitud de Y (respuesta):", length(y), "\n")## Longitud de Y (respuesta): 12002 Pregunta 1: ¿Hay multicolinealidad?
En el aprendizaje estadístico supervisado, la multicolinealidad ocurre cuando los predictores presentan correlaciones altas entre sí, lo que hace inestable o imposible la estimación de los parámetros por Mínimos Cuadrados Ordinarios (MCO). Recordando de la Clase 1 que \(\hat{\beta}_{MCO}\) es el BLUE, esto no es necesariamente bueno cuando existe multicolinealidad: aunque es insesgado, su varianza puede ser extremadamente alta, deteriorando el ECM total: \(ECM(\hat{\theta}) = V(\hat{\theta}) + B(\hat{\theta})^2\).
# --- Argumento 1: Dimensionalidad (p >> n) ---
p <- ncol(X)
n <- nrow(X)
cat("Número de predictores (p):", p, "\n")## Número de predictores (p): 5000cat("Número de observaciones (n):", n, "\n")## Número de observaciones (n): 1200cat("Razón p/n:", round(p / n, 2), "\n\n")## Razón p/n: 4.17# Cuando p > n, la matriz X'X (5000x5000) es singular -> no invertible
# -> MCO no tiene solución única -> multicolinealidad perfecta estructural
cat("Como p =", p, "> n =", n, ":\n")## Como p = 5000 > n = 1200 :cat("La matriz X'X es singular y MCO no tiene solución única.\n\n")## La matriz X'X es singular y MCO no tiene solución única.# --- Argumento 2: Correlaciones entre genes (muestra de 100) ---
set.seed(49)
muestra_genes <- sample(1:p, 100)
cor_muestra <- cor(X[, muestra_genes])
cor_upper <- cor_muestra[upper.tri(cor_muestra)]
prop_cor_alta <- mean(abs(cor_upper) > 0.7)
prop_cor_media <- mean(abs(cor_upper) > 0.3)
cat("En una muestra aleatoria de 100 genes:\n")## En una muestra aleatoria de 100 genes:cat(" Proporción de pares con |r| > 0.7:", round(prop_cor_alta, 4), "\n")## Proporción de pares con |r| > 0.7: 0cat(" Proporción de pares con |r| > 0.3:", round(prop_cor_media, 4), "\n")## Proporción de pares con |r| > 0.3: 0# --- Argumento 3: Visualización ---
corrplot(cor_muestra, method = "color", tl.pos = "n",
title = "Mapa de correlaciones entre 100 genes (muestra aleatoria)",
mar = c(0, 0, 2, 0))
Conclusión (P1): Sí existe multicolinealidad severa en los datos, por tres razones:
Estructural: El número de predictores (p = 5000 genes) supera ampliamente el número de observaciones (n = 1200), haciendo que la matriz \(X'X\) sea singular y que MCO no tenga solución única. Este es un caso de multicolinealidad perfecta por construcción.
Biológica: En el contexto genómico, los genes no actúan de forma independiente: muchos participan en las mismas rutas metabólicas o procesos celulares, generando correlaciones altas entre columnas de la matriz de diseño.
Consecuencia estadística: Aunque \(\hat{\beta}_{MCO}\) es el BLUE (Clase 1), su varianza es extremadamente alta en este contexto, deteriorando el ECM total. Esto justifica el uso de métodos de regularización (Ridge y Lasso) que introducen sesgo controlado a cambio de una reducción significativa en la varianza.
3 Pregunta 2: División entrenamiento / prueba
Como se discutió en la Clase 2, cuando el interés es la predicción de nuevas observaciones o la selección del modelo, se requiere un conjunto de prueba para medir el rendimiento del modelo \(P_{\hat{\theta}}\), es decir, una validación externa.
set.seed(49) # Semilla fijada para reproducibilidad
n_total <- nrow(datos)
indices_train <- sample(1:n_total, size = 1000, replace = FALSE)
# Conjuntos de entrenamiento (train)
X_train <- X[indices_train, ]
y_train <- y[indices_train]
# Conjuntos de prueba (test)
X_test <- X[-indices_train, ]
y_test <- y[-indices_train]
cat("Dimensiones entrenamiento - X:", dim(X_train), "| Y:", length(y_train), "\n")## Dimensiones entrenamiento - X: 1000 5000 | Y: 1000cat("Dimensiones prueba - X:", dim(X_test), "| Y:", length(y_test), "\n")## Dimensiones prueba - X: 200 5000 | Y: 200# Verificar distribución similar entre train y test
cat("\nResumen Y entrenamiento:\n"); print(summary(y_train))##
## Resumen Y entrenamiento:## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -13.90939 -3.07947 -0.06371 -0.13811 2.84513 16.84200cat("\nResumen Y prueba:\n"); print(summary(y_test))##
## Resumen Y prueba:## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -11.7506 -3.1766 -0.1332 -0.3102 2.3179 13.8928La variable respuesta presenta distribuciones similares en ambos conjuntos, lo cual indica que la partición aleatoria fue representativa.
4 Pregunta 3: Selección de λ mediante validación cruzada
Como se presentó en la Clase 2, la regresión Ridge y Lasso hacen uso de validación externa con validación cruzada en 10 folds (10-fold CV) para determinar el valor de λ. Este método divide los 1000 datos de entrenamiento en 10 grupos iguales, ajusta el modelo con 9 grupos y evalúa en el restante, rotando hasta cubrir todos, y promedia el ECM de las 10 iteraciones. Se prefiere sobre la validación regular por dos razones:
- Menor sesgo: usa casi todos los datos en entrenamiento en cada iteración.
- Menor varianza: frente a LOOCV, los ECM de cada fold están menos correlacionados.
El valor seleccionado es lambda.min: el λ que minimiza
el ECM de validación cruzada.
set.seed(49)
# --- Ridge (alpha = 0) ---
cv_ridge <- cv.glmnet(X_train, y_train,
alpha = 0,
nfolds = 10,
standardize = TRUE,
type.measure = "mse")
lambda_r <- cv_ridge$lambda.min
cat("=== RIDGE ===\n")## === RIDGE ===cat("Lambda óptimo (lambda.min):", round(lambda_r, 6), "\n")## Lambda óptimo (lambda.min): 12.44889cat("Lambda (lambda.1se):", round(cv_ridge$lambda.1se, 6), "\n")## Lambda (lambda.1se): 23.87586cat("ECM mínimo en CV:", round(min(cv_ridge$cvm), 4), "\n\n")## ECM mínimo en CV: 17.0308# --- Lasso (alpha = 1) ---
cv_lasso <- cv.glmnet(X_train, y_train,
alpha = 1,
nfolds = 10,
standardize = TRUE,
type.measure = "mse")
lambda_l <- cv_lasso$lambda.min
cat("=== LASSO ===\n")## === LASSO ===cat("Lambda óptimo (lambda.min):", round(lambda_l, 6), "\n")## Lambda óptimo (lambda.min): 0.080022cat("Lambda (lambda.1se):", round(cv_lasso$lambda.1se, 6), "\n")## Lambda (lambda.1se): 0.100977cat("ECM mínimo en CV:", round(min(cv_lasso$cvm), 4), "\n")## ECM mínimo en CV: 1.2263# Gráficas de CV (ECM vs log(lambda))
par(mfrow = c(1, 2))
plot(cv_ridge, main = "Ridge: ECM vs log(λ)")
plot(cv_lasso, main = "Lasso: ECM vs log(λ)")
par(mfrow = c(1, 1))Interpretación: Las gráficas muestran la curva de
ECM de validación cruzada en función de log(λ). Las líneas verticales
punteadas indican lambda.min (izquierda) y
lambda.1se (derecha). Se selecciona lambda.min
ya que el objetivo del taller es minimizar el ECM de predicción.
5 Pregunta 4: Ajuste de modelos Ridge y Lasso
Con los λ óptimos obtenidos en (3), se ajustan ambos modelos sobre los 1000 datos de entrenamiento. La diferencia clave entre Ridge (penalización L2) y Lasso (penalización L1) es que Lasso lleva coeficientes exactamente a cero, actuando como selector de variables, mientras que Ridge los encoge hacia cero pero nunca los anula.
# Modelo Ridge
ridge_model <- glmnet(X_train, y_train,
alpha = 0,
lambda = lambda_r,
standardize = TRUE)
# Modelo Lasso
lasso_model <- glmnet(X_train, y_train,
alpha = 1,
lambda = lambda_l,
standardize = TRUE)
# Coeficientes
coef_ridge <- coef(ridge_model)
coef_lasso <- coef(lasso_model)
n_nonzero_ridge <- sum(coef_ridge != 0) - 1 # excluir intercepto
n_nonzero_lasso <- sum(coef_lasso != 0) - 1
cat("Coeficientes ≠ 0 en Ridge:", n_nonzero_ridge, "de 5000 genes\n")## Coeficientes ≠ 0 en Ridge: 5000 de 5000 genescat("Coeficientes ≠ 0 en Lasso:", n_nonzero_lasso, "de 5000 genes\n\n")## Coeficientes ≠ 0 en Lasso: 76 de 5000 genes# Top 10 genes más importantes según Lasso
top10_lasso <- sort(abs(coef_lasso[-1, 1]), decreasing = TRUE)[1:10]
cat("Top 10 genes con mayor coeficiente absoluto (Lasso):\n")## Top 10 genes con mayor coeficiente absoluto (Lasso):print(round(top10_lasso, 6))## V17 V11 V5 V20 V7 V18 V19 V2
## 0.993552 0.976170 0.944601 0.944372 0.936138 0.930369 0.930215 0.928828
## V3 V15
## 0.919577 0.917947Como se esperaba, Ridge conserva todos los genes con coeficientes distintos de cero (encoge pero no elimina), mientras que Lasso selecciona únicamente un subconjunto de genes relevantes.
6 Pregunta 5: Selección del mejor modelo
Para propósitos de predicción, se usa el ECM en los 200 datos de prueba — datos que los modelos nunca vieron durante el entrenamiento. Como se enfatizó en la Clase 2, el ECMTest es la medida apropiada de validación externa para seleccionar entre modelos. Un ECMTrain bajo no garantiza buen desempeño en prueba (sobreajuste).
\[ECM_{Test} = \frac{1}{200}\sum_{i=1}^{200}(y_i - \hat{f}(x_i))^2\]
# Predicciones en datos de prueba
pred_ridge <- predict(ridge_model, newx = X_test)
pred_lasso <- predict(lasso_model, newx = X_test)
# ECM en prueba
ecm_ridge <- mean((y_test - pred_ridge)^2)
ecm_lasso <- mean((y_test - pred_lasso)^2)
cat("ECM Ridge en datos de prueba:", round(ecm_ridge, 4), "\n")## ECM Ridge en datos de prueba: 17.831cat("ECM Lasso en datos de prueba:", round(ecm_lasso, 4), "\n")## ECM Lasso en datos de prueba: 1.251cat("\nModelo seleccionado:", ifelse(ecm_lasso < ecm_ridge, "LASSO", "RIDGE"), "\n")##
## Modelo seleccionado: LASSO# Tabla comparativa
tabla <- data.frame(
Modelo = c("Ridge", "Lasso"),
Lambda = c(round(lambda_r, 6), round(lambda_l, 6)),
ECM_CV = c(round(min(cv_ridge$cvm), 4), round(min(cv_lasso$cvm), 4)),
ECM_Prueba = c(round(ecm_ridge, 4), round(ecm_lasso, 4)),
Genes_activos = c(n_nonzero_ridge, n_nonzero_lasso)
)
knitr::kable(tabla, caption = "Comparación Ridge vs Lasso")| Modelo | Lambda | ECM_CV | ECM_Prueba | Genes_activos |
|---|---|---|---|---|
| Ridge | 12.448893 | 17.0308 | 17.831 | 5000 |
| Lasso | 0.080022 | 1.2263 | 1.251 | 76 |
Conclusión (P5): El modelo Lasso presenta el menor ECM en los datos de prueba (1.251), por lo que es seleccionado. La penalización L1 del Lasso es más adecuada en este problema al realizar selección automática de variables, conservando únicamente los genes con mayor capacidad predictiva.
7 Pregunta 6: Modelo final con los 1200 datos
Una vez seleccionado el modelo y estimado el λ óptimo en (3), se ajusta el modelo final usando la totalidad de los 1200 datos. El λ no se re-estima: se usa el mismo valor obtenido por validación cruzada. Este paso aprovecha toda la información disponible para obtener estimaciones más precisas de los coeficientes.
# Parámetros del modelo ganador
modelo_sel <- ifelse(ecm_lasso < ecm_ridge, "lasso", "ridge")
lambda_fin <- ifelse(ecm_lasso < ecm_ridge, lambda_l, lambda_r)
alpha_fin <- ifelse(ecm_lasso < ecm_ridge, 1, 0)
cat("Modelo seleccionado:", toupper(modelo_sel), "\n")## Modelo seleccionado: LASSOcat("Lambda final (del paso 3):", round(lambda_fin, 6), "\n\n")## Lambda final (del paso 3): 0.080022# Ajuste con los 1200 datos completos
final_model <- glmnet(X, y,
alpha = alpha_fin,
lambda = lambda_fin,
standardize = TRUE)
coef_final <- coef(final_model)
n_nonzero_final <- sum(coef_final != 0) - 1
cat("Genes con coeficiente ≠ 0 en modelo final:", n_nonzero_final, "de 5000\n")## Genes con coeficiente ≠ 0 en modelo final: 59 de 5000# Nombres de genes seleccionados
genes_sel <- which(coef_final[-1] != 0)
cat("Primeros 20 genes seleccionados:\n")## Primeros 20 genes seleccionados:print(head(names(genes_sel), 20))## NULL8 Pregunta 7: Trazas de coeficientes
Las trazas muestran cómo evolucionan los coeficientes de los 5000 genes a medida que varía λ. Para valores grandes de λ (derecha), la penalización es alta y la mayoría de coeficientes permanecen en cero. Al disminuir λ (izquierda), los genes van ingresando al modelo. La línea roja indica el λ óptimo seleccionado por validación cruzada.
# Ajustar en toda la trayectoria de lambda
modelo_path <- glmnet(X, y, alpha = alpha_fin, standardize = TRUE)
# Gráfica de trazas
plot(modelo_path,
xvar = "lambda",
label = FALSE,
main = paste("Trazas de coeficientes —", toupper(modelo_sel), "(n = 1200)"),
xlab = expression(log(lambda)),
ylab = "Coeficientes")
# Línea vertical en lambda óptimo
abline(v = log(lambda_fin), col = "red", lty = 2, lwd = 2)
legend("topright",
legend = paste0("λ óptimo = ", round(lambda_fin, 5)),
col = "red", lty = 2, lwd = 2, bty = "n")
Interpretación: Cada línea de color representa un gen. A medida que λ disminuye (de derecha a izquierda), los genes con mayor asociación con la variable respuesta empiezan a adquirir coeficientes distintos de cero. En el λ óptimo (línea roja), solo 59 genes tienen coeficientes activos de los 5000 disponibles.
9 Pregunta 8: Resumen de resultados
El objetivo del estudio fue determinar cuáles de los 5000 genes del perfil genómico son relevantes para predecir la efectividad de tratamientos anticancerígenos en 1200 líneas celulares. Dado que el número de predictores (p = 5000) supera ampliamente el de observaciones (n = 1200), la matriz \(X'X\) es singular y la estimación por MCO es inviable, evidenciando multicolinealidad perfecta estructural. A esto se suma que, biológicamente, los genes co-expresados en rutas metabólicas comunes presentan correlaciones altas entre sí. Por estas razones, se recurrió a métodos de regularización que introducen un sesgo controlado a cambio de reducir la varianza, mejorando el ECM total. Los datos se dividieron en entrenamiento (n = 1000) y prueba (n = 200), y mediante validación cruzada de 10 pliegues sobre el conjunto de entrenamiento se estimaron los parámetros de penalización óptimos: λᵣ = 12.4489 para Ridge y λₗ = 0.08 para Lasso. Al evaluar ambos modelos en los datos de prueba, el modelo LASSO obtuvo el menor ECM (1.251 frente a 17.831 del modelo alternativo), siendo seleccionado como el más apropiado para predicción. Ajustado finalmente sobre las 1200 observaciones, el modelo identificó 59 genes con coeficientes distintos de cero, lo que sugiere que solo una fracción del perfil genómico contribuye de forma relevante a la predicción de la efectividad del tratamiento. Estos genes constituyen candidatos de interés biológico para investigaciones futuras sobre blancos terapéuticos en oncología.
10 Información de sesión
sessionInfo()## R version 4.5.0 (2025-04-11 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26200)
##
## Matrix products: default
## LAPACK version 3.12.1
##
## locale:
## [1] LC_COLLATE=Spanish_Mexico.utf8 LC_CTYPE=Spanish_Mexico.utf8
## [3] LC_MONETARY=Spanish_Mexico.utf8 LC_NUMERIC=C
## [5] LC_TIME=Spanish_Mexico.utf8
##
## time zone: America/Bogota
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] corrplot_0.95 ggplot2_4.0.2 glmnet_4.1-10 Matrix_1.7-3
##
## loaded via a namespace (and not attached):
## [1] gtable_0.3.6 jsonlite_2.0.0 dplyr_1.1.4 compiler_4.5.0
## [5] tidyselect_1.2.1 Rcpp_1.0.14 jquerylib_0.1.4 splines_4.5.0
## [9] scales_1.4.0 yaml_2.3.10 fastmap_1.2.0 lattice_0.22-6
## [13] R6_2.6.1 generics_0.1.4 shape_1.4.6.1 knitr_1.51
## [17] iterators_1.0.14 tibble_3.2.1 bslib_0.9.0 pillar_1.10.2
## [21] RColorBrewer_1.1-3 rlang_1.1.6 cachem_1.1.0 xfun_0.52
## [25] sass_0.4.10 S7_0.2.1 cli_3.6.5 withr_3.0.2
## [29] magrittr_2.0.3 digest_0.6.37 foreach_1.5.2 grid_4.5.0
## [33] rstudioapi_0.18.0 lifecycle_1.0.4 vctrs_0.6.5 evaluate_1.0.3
## [37] glue_1.8.0 farver_2.1.2 codetools_0.2-20 survival_3.8-3
## [41] rmarkdown_2.30 tools_4.5.0 pkgconfig_2.0.3 htmltools_0.5.8.1