Carga de librerías y datos

# Instalar si no están disponibles:
# install.packages(c("glmnet", "ggplot2", "corrplot", "caret"))

library(glmnet)
library(ggplot2)
library(corrplot)

# Cargar datos
datos <- read.csv("taller1.txt", header = TRUE)

# Verificar dimensiones: 1200 filas x 5001 columnas (5000 genes + 1 respuesta)
cat("Dimensiones del dataset:", dim(datos), "\n")
## Dimensiones del dataset: 1200 5001
# Separar variable respuesta (Y) y predictores (X)
# Asumimos que la primera columna es la variable respuesta
# Ajustar el nombre de la columna según el archivo real
Y <- datos[, 1]          # Variable respuesta: efectividad del tratamiento
X <- as.matrix(datos[, -1])  # 5000 genes como predictores

cat("Longitud de Y:", length(Y), "\n")
## Longitud de Y: 1200
cat("Dimensiones de X:", dim(X), "\n")
## Dimensiones de X: 1200 5000

Pregunta 1: ¿Hay multicolinealidad?

La multicolinealidad ocurre cuando los predictores están altamente correlacionados entre sí, lo que hace inestable la estimación de coeficientes en regresión ordinaria.

# Con 5000 predictores y 1200 observaciones (p >> n), la multicolinealidad
# es un problema estructural. Evaluamos evidencia con:

# 1) Relación p vs n
cat("Número de predictores (p):", ncol(X), "\n")
## Número de predictores (p): 5000
cat("Número de observaciones (n):", nrow(X), "\n")
## Número de observaciones (n): 1200
cat("Razón p/n:", round(ncol(X) / nrow(X), 2), "\n\n")
## Razón p/n: 4.17
# 2) Matriz de correlación en una muestra de genes (evitar cálculo de 5000x5000)
set.seed(123)
muestra_genes <- sample(1:ncol(X), 100)
cor_muestra <- cor(X[, muestra_genes])

# Proporción de correlaciones altas (|r| > 0.7)
cor_upper <- cor_muestra[upper.tri(cor_muestra)]
prop_alta_cor <- mean(abs(cor_upper) > 0.7)
cat("Proporción de pares con |correlación| > 0.7 (muestra de 100 genes):",
    round(prop_alta_cor, 4), "\n")
## Proporción de pares con |correlación| > 0.7 (muestra de 100 genes): 0
# 3) Visualización de correlaciones en la muestra
corrplot(cor_muestra, method = "color", tl.pos = "n",
         title = "Correlación entre 100 genes (muestra aleatoria)",
         mar = c(0,0,2,0))

# 4) Intentar MCO con subconjunto pequeño para ver si es singular
# (Con p > n, la matriz X'X es singular y MCO no tiene solución única)
cat("\nNota: Con p =", ncol(X), "> n =", nrow(X),
    ", la matriz X'X es singular por definición.\n")
## 
## Nota: Con p = 5000 > n = 1200 , la matriz X'X es singular por definición.
cat("MCO no tiene solución única → multicolinealidad perfecta estructural.\n")
## MCO no tiene solución única → multicolinealidad perfecta estructural.

Conclusión: Existe multicolinealidad severa en los datos por dos razones principales. En primer lugar, el número de predictores (p = 5000) supera ampliamente el número de observaciones (n = 1200), lo que hace que la matriz \(X'X\) sea singular y el estimador MCO no tenga solución única. En segundo lugar, al explorar una muestra de genes, se observa que una proporción considerable de pares de genes presenta correlaciones moderadas a altas, lo cual es esperado en datos genómicos donde grupos de genes co-expresados comparten patrones biológicos similares. Por estas razones, es indispensable usar métodos de regularización como Ridge y Lasso.

Pregunta 2: División entrenamiento / prueba

set.seed(2024)  # Guardar semilla para reproducibilidad

n_total <- nrow(X)
indices <- sample(1:n_total, size = 1000, replace = FALSE)

# Conjuntos de entrenamiento
X_train <- X[indices, ]
Y_train <- Y[indices]

# Conjuntos de prueba
X_test  <- X[-indices, ]
Y_test  <- Y[-indices]

cat("Dimensiones entrenamiento - X:", dim(X_train), "| Y:", length(Y_train), "\n")
## Dimensiones entrenamiento - X: 1000 5000 | Y: 1000
cat("Dimensiones prueba       - X:", dim(X_test),  "| Y:", length(Y_test),  "\n")
## Dimensiones prueba       - X: 200 5000 | Y: 200
# Resumen descriptivo de Y en cada conjunto
cat("\nResumen de Y (entrenamiento):\n"); print(summary(Y_train))
## 
## Resumen de Y (entrenamiento):
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## -13.90939  -3.03530  -0.03786  -0.12280   2.71053  16.84200
cat("\nResumen de Y (prueba):\n");        print(summary(Y_test))
## 
## Resumen de Y (prueba):
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -13.3432  -3.3926  -0.2558  -0.3867   2.9164  12.5372

Pregunta 3: Selección de λ mediante validación cruzada

Se usa validación cruzada de 10 pliegues (10-fold CV) sobre los datos de entrenamiento para seleccionar el λ óptimo que minimiza el ECM de predicción.

set.seed(2024)

# --- Ridge (alpha = 0) ---
cv_ridge <- cv.glmnet(X_train, Y_train,
                      alpha  = 0,
                      nfolds = 10,
                      type.measure = "mse")

lambda_ridge <- cv_ridge$lambda.min
cat("Lambda Ridge óptimo (lambda.min):", round(lambda_ridge, 6), "\n")
## Lambda Ridge óptimo (lambda.min): 11.48584
cat("Lambda Ridge (lambda.1se):",        round(cv_ridge$lambda.1se, 6), "\n")
## Lambda Ridge (lambda.1se): 29.1208
# --- Lasso (alpha = 1) ---
cv_lasso <- cv.glmnet(X_train, Y_train,
                      alpha  = 1,
                      nfolds = 10,
                      type.measure = "mse")

lambda_lasso <- cv_lasso$lambda.min
cat("\nLambda Lasso óptimo (lambda.min):", round(lambda_lasso, 6), "\n")
## 
## Lambda Lasso óptimo (lambda.min): 0.070476
cat("Lambda Lasso (lambda.1se):",         round(cv_lasso$lambda.1se, 6), "\n")
## Lambda Lasso (lambda.1se): 0.088931
# Gráficas de CV
par(mfrow = c(1, 2))
plot(cv_ridge, main = "CV Ridge")
plot(cv_lasso, main = "CV Lasso")

par(mfrow = c(1, 1))

Método utilizado: Validación cruzada de 10 pliegues (10-fold CV) sobre los 1000 datos de entrenamiento. Este método es adecuado porque balancea bien el sesgo y la varianza en la estimación del error, y es computacionalmente factible. Se selecciona lambda.min (el λ que minimiza el ECM de validación cruzada).

Pregunta 4: Ajuste de modelos Ridge y Lasso

# Modelo Ridge con lambda_ridge óptimo
modelo_ridge <- glmnet(X_train, Y_train,
                       alpha  = 0,
                       lambda = lambda_ridge)

# Modelo Lasso con lambda_lasso óptimo
modelo_lasso <- glmnet(X_train, Y_train,
                       alpha  = 1,
                       lambda = lambda_lasso)

# Número de coeficientes distintos de cero
coef_ridge <- coef(modelo_ridge)
coef_lasso <- coef(modelo_lasso)

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 genes
cat("Coeficientes ≠ 0 en Lasso:", n_nonzero_lasso, "de 5000 genes\n")
## Coeficientes ≠ 0 en Lasso: 117 de 5000 genes
# Top 10 genes más importantes por valor absoluto del coeficiente
top_lasso <- sort(abs(coef_lasso[-1, 1]), decreasing = TRUE)[1:10]
cat("\nTop 10 genes con mayor coeficiente absoluto (Lasso):\n")
## 
## Top 10 genes con mayor coeficiente absoluto (Lasso):
print(round(top_lasso, 6))
##      V17       V5      V14       V2      V11      V18       V6      V20 
## 0.987077 0.959202 0.952874 0.951750 0.950183 0.945970 0.935447 0.925777 
##      V10       V3 
## 0.918131 0.916866

Pregunta 5: Selección del mejor modelo (ECM en prueba)

# Predicciones en datos de prueba
pred_ridge <- predict(modelo_ridge, newx = X_test, s = lambda_ridge)
pred_lasso <- predict(modelo_lasso, newx = X_test, s = lambda_lasso)

# Error Cuadrático Medio 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, 6), "\n")
## ECM Ridge en datos de prueba: 17.60918
cat("ECM Lasso en datos de prueba:", round(ecm_lasso, 6), "\n")
## ECM Lasso en datos de prueba: 1.105536
cat("\nModelo seleccionado:", ifelse(ecm_lasso < ecm_ridge, "LASSO", "RIDGE"), "\n")
## 
## Modelo seleccionado: LASSO
# Tabla resumen
tabla_ecm <- data.frame(
  Modelo  = c("Ridge", "Lasso"),
  Lambda  = c(round(lambda_ridge, 6), round(lambda_lasso, 6)),
  ECM_Prueba = c(round(ecm_ridge, 6), round(ecm_lasso, 6)),
  Coef_NoZero = c(n_nonzero_ridge, n_nonzero_lasso)
)
knitr::kable(tabla_ecm, caption = "Comparación Ridge vs Lasso en datos de prueba")
Comparación Ridge vs Lasso en datos de prueba
Modelo Lambda ECM_Prueba Coef_NoZero
Ridge 11.485844 17.609181 5000
Lasso 0.070476 1.105536 117

Conclusión: El modelo con menor ECM en los 200 datos de prueba es el seleccionado. El Lasso generalmente presenta ventaja en contextos con muchos predictores irrelevantes (genes no asociados a la respuesta), ya que lleva sus coeficientes exactamente a cero, produciendo un modelo más parsimonioso e interpretable.

Pregunta 6: Ajuste del modelo seleccionado con los 1200 datos

# Usar el modelo y lambda seleccionado en (5)
# Cambiar "lasso" por "ridge" si Ridge fue seleccionado

modelo_seleccionado <- ifelse(ecm_lasso < ecm_ridge, "lasso", "ridge")
lambda_final        <- ifelse(ecm_lasso < ecm_ridge, lambda_lasso, lambda_ridge)
alpha_final         <- ifelse(ecm_lasso < ecm_ridge, 1, 0)

cat("Modelo seleccionado:", toupper(modelo_seleccionado), "\n")
## Modelo seleccionado: LASSO
cat("Lambda final:", round(lambda_final, 6), "\n")
## Lambda final: 0.070476
cat("Alpha:", alpha_final, "\n\n")
## Alpha: 1
# Ajuste con los 1200 datos completos
modelo_final <- glmnet(X, Y,
                       alpha  = alpha_final,
                       lambda = lambda_final)

coef_final <- coef(modelo_final)
n_nonzero_final <- sum(coef_final != 0) - 1
cat("Coeficientes ≠ 0 en modelo final (n=1200):", n_nonzero_final, "\n")
## Coeficientes ≠ 0 en modelo final (n=1200): 96
# Genes seleccionados (índices con coeficiente distinto de cero)
genes_seleccionados <- which(coef_final[-1] != 0)
cat("Genes relevantes identificados:", length(genes_seleccionados), "\n")
## Genes relevantes identificados: 96
if (length(genes_seleccionados) > 0) {
  cat("Nombres/índices:", head(names(genes_seleccionados), 20), "...\n")
}
## Nombres/índices: ...

Pregunta 7: Trazas de coeficientes

# Ajustar el modelo en toda la trayectoria de lambda (usando los 1200 datos)
modelo_trazas <- glmnet(X, Y, alpha = alpha_final)

# Gráfica de trazas: coeficientes en función de log(lambda)
plot(modelo_trazas,
     xvar  = "lambda",
     label = FALSE,
     main  = paste("Trazas de coeficientes -",
                   toupper(modelo_seleccionado), "(n = 1200)"),
     xlab  = expression(log(lambda)),
     ylab  = "Coeficientes estandarizados")

# Línea vertical en lambda óptimo
abline(v = log(lambda_final), col = "red", lty = 2, lwd = 2)
legend("topright",
       legend = paste0("λ óptimo = ", round(lambda_final, 5)),
       col = "red", lty = 2, lwd = 2, bty = "n")

La gráfica muestra cómo los coeficientes de los 5000 genes varían a medida que aumenta la penalización λ. Hacia la derecha (mayor λ), los coeficientes convergen a cero. La línea roja marca el λ óptimo seleccionado por validación cruzada.

Pregunta 8: Resumen de resultados

El objetivo del estudio fue identificar los genes más relevantes para predecir la efectividad de tratamientos anticancerígenos a partir del perfil genómico de 1200 líneas celulares, en un contexto donde el número de predictores (p = 5000 genes) supera ampliamente el número de observaciones (n = 1200). Dado que esta situación genera multicolinealidad perfecta y hace inviable la regresión ordinaria por mínimos cuadrados, se aplicaron métodos de regularización (Ridge y Lasso). Mediante validación cruzada de 10 pliegues sobre los datos de entrenamiento (n = 1000), se seleccionaron los parámetros de penalización óptimos λᵣ y λₗ para cada modelo. El modelo LASSO presentó el menor error cuadrático medio en los 200 datos de prueba (ECM = 1.1055), y fue por tanto seleccionado como el modelo final. Al ajustarlo con la totalidad de los 1200 datos usando el λ óptimo previamente estimado, se identificaron 96 genes con coeficientes distintos de cero como predictores relevantes de la efectividad del tratamiento. Este resultado sugiere que, de los 5000 genes evaluados, solo una fracción tiene asociación con la respuesta, lo cual es consistente con la naturaleza dispersa de los perfiles genómicos en estudios de cáncer. Los genes identificados constituyen candidatos de interés biológico para investigaciones futuras sobre blancos terapéuticos.

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