Carga de datos

En esta sección cargamos el archivo de texto que contiene las observaciones

datos <- read.table("taller1.txt",
                    header = TRUE,
                    sep = ",",
                    dec = ".")

Realizamos descriptivos basicos para ver como se comportan nuestras datos y darle una primera ojeada a ese

## descriptivos basicos 

str(datos)
## 'data.frame':    1200 obs. of  5001 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 ...
##  $ V6   : num  -0.699 0.996 -0.693 -0.103 0.604 ...
##  $ V7   : num  0.928 0.359 0.187 -0.42 -1.166 ...
##  $ V8   : num  0.524 -0.742 1.071 -1.29 2.098 ...
##  $ V9   : num  -1.183 -0.376 -1.467 -0.295 -0.563 ...
##  $ V10  : num  -0.635 -0.511 0.944 0.97 0.771 ...
##  $ V11  : num  -0.151 0.801 -1.187 0.431 0.217 ...
##  $ V12  : num  -1.5652 -0.2584 -0.0937 -1.4607 0.2471 ...
##  $ V13  : num  -0.406 -0.883 -0.379 0.722 0.207 ...
##  $ V14  : num  -0.0649 0.194 1.8154 -1.0807 2.035 ...
##  $ V15  : num  0.3233 0.3191 0.0686 0.838 -0.208 ...
##  $ V16  : num  0.278 -1.603 -0.321 -1.226 -0.777 ...
##  $ V17  : num  0.7861 -1.3184 -0.0459 0.0779 -0.7658 ...
##  $ V18  : num  -0.0388 -0.5845 1.2893 0.3011 0.5539 ...
##  $ V19  : num  1.78 1.352 1.523 0.148 -0.15 ...
##  $ V20  : num  -1.56327 -1.16939 0.00337 1.46111 0.81096 ...
##  $ V21  : num  -0.8115 1.3911 -0.0119 -0.8961 -1.0879 ...
##  $ V22  : num  -0.982 -0.382 -0.426 1.691 -0.334 ...
##  $ V23  : num  0.963 -0.713 0.53 0.826 -1.407 ...
##  $ V24  : num  1.45 1.756 -0.622 0.503 -0.613 ...
##  $ V25  : num  0.505 -1.006 0.755 1.352 0.641 ...
##  $ V26  : num  -0.194 0.258 -0.538 -1.179 0.901 ...
##  $ V27  : num  -0.797 2.186 -1.521 0.777 -0.789 ...
##  $ V28  : num  -0.0294 -0.3975 -0.9008 -1.4824 0.1352 ...
##  $ V29  : num  -0.406 0.38 -0.244 -2.189 1.216 ...
##  $ V30  : num  -0.2072 -1.5692 0.8002 -0.0384 0.768 ...
##  $ V31  : num  0.672 0.552 -0.297 0.517 0.275 ...
##  $ V32  : num  -0.054 0.0371 -0.9252 1.0066 1.0011 ...
##  $ V33  : num  0.7735 0.2906 0.8592 -0.7595 -0.0765 ...
##  $ V34  : num  0.378 0.598 0.577 -0.579 1.587 ...
##  $ V35  : num  0.621 0.569 -0.506 -0.719 -1.031 ...
##  $ V36  : num  -0.869 -0.557 0.983 -0.309 0.133 ...
##  $ V37  : num  2.509 -1.649 1.851 1.659 -0.684 ...
##  $ V38  : num  -0.8007 0.432 -0.0697 -1.1298 0.9909 ...
##  $ V39  : num  -0.238 -0.169 -1.33 -1.949 1.451 ...
##  $ V40  : num  -0.0161 0.6645 -3.0278 -0.4405 0.825 ...
##  $ V41  : num  1.688491 -0.02114 -0.0205 -1.332804 -0.000505 ...
##  $ V42  : num  -0.8897 0.1875 0.5075 0.2903 0.0319 ...
##  $ V43  : num  -0.803 -2.459 0.826 2.377 -0.47 ...
##  $ V44  : num  0.581 0.188 0.371 0.378 -0.531 ...
##  $ V45  : num  -0.929 1.27 -0.166 -0.439 0.336 ...
##  $ V46  : num  -2.041 0.733 2.298 -1.119 1.391 ...
##  $ V47  : num  0.65 0.574 1.025 -0.18 0.566 ...
##  $ V48  : num  1.411 -0.448 1.395 -0.404 -1.494 ...
##  $ V49  : num  0.677 0.378 0.261 0.766 1.736 ...
##  $ V50  : num  2.0016 0.5263 0.0307 0.313 1.1568 ...
##  $ V51  : num  -0.3883 0.0274 -0.2761 -0.0867 2.1477 ...
##  $ V52  : num  -1.042 -0.317 0.61 -0.881 1.822 ...
##  $ V53  : num  0.675 -0.697 -0.303 0.822 0.359 ...
##  $ V54  : num  0.3602 0.0088 0.4656 -0.2431 0.5458 ...
##  $ V55  : num  -0.0953 0.5758 0.444 -0.7138 -0.7698 ...
##  $ V56  : num  0.4344 0.2179 -0.4274 0.0636 1.408 ...
##  $ V57  : num  -1.3637 -1.1178 0.5863 0.0356 1.625 ...
##  $ V58  : num  -2.201 1.362 0.518 -1.204 -0.518 ...
##  $ V59  : num  0.514 -0.648 1.727 -0.778 0.987 ...
##  $ V60  : num  0.6205 -0.0438 -0.2588 -1.6446 -0.8323 ...
##  $ V61  : num  -1.211 -1.001 0.235 1.034 -0.392 ...
##  $ V62  : num  -0.6274 0.9167 0.6837 -0.0141 -0.344 ...
##  $ V63  : num  0.50569 2.16501 -0.00697 0.51938 1.71001 ...
##  $ V64  : num  0.0931 -0.1291 0.9979 0.7016 -1.3228 ...
##  $ V65  : num  0.693 -0.735 2.059 -0.505 -0.394 ...
##  $ V66  : num  1.5969 0.844 -0.0487 -1.1216 -0.2766 ...
##  $ V67  : num  0.515 0.11 1.273 -0.481 1.302 ...
##  $ V68  : num  -0.3482 0.8303 0.0248 -0.0231 0.196 ...
##  $ V69  : num  0.781 -0.133 1.927 1.159 1.385 ...
##  $ V70  : num  0.257 -2.536 0.321 -0.635 1.611 ...
##  $ V71  : num  1.312 0.47 1.545 0.827 -0.427 ...
##  $ V72  : num  1.667 1.368 -0.951 0.737 -0.121 ...
##  $ V73  : num  1.253 1.155 -0.377 0.237 0.445 ...
##  $ V74  : num  -0.873 -0.593 -0.661 -1.946 0.43 ...
##  $ V75  : num  -0.584 1.279 -1.487 -0.241 -1.008 ...
##  $ V76  : num  0.5732 0.0183 -0.022 -0.4278 -0.4776 ...
##  $ V77  : num  0.838 -0.974 -1.884 -0.605 -1.142 ...
##  $ V78  : num  -1.024 0.822 0.335 -1.691 -0.82 ...
##  $ V79  : num  -2.1596 1.2989 0.0894 -0.6549 0.4251 ...
##  $ V80  : num  1.481 -1.232 1.108 -0.374 -0.105 ...
##  $ V81  : num  0.352 -0.654 -0.522 -0.418 -0.286 ...
##  $ V82  : num  -0.388 1.148 0.348 0.648 -2.347 ...
##  $ V83  : num  1.186 0.847 -1.201 1.447 -0.222 ...
##  $ V84  : num  0.069 0.5856 -0.5185 -0.0754 -0.2746 ...
##  $ V85  : num  1.998 0.543 -0.771 0.265 1.095 ...
##  $ V86  : num  -0.628 -0.419 -0.174 -0.906 -1.522 ...
##  $ V87  : num  0.509 -0.541 1.634 -0.283 0.249 ...
##  $ V88  : num  0.4983 0.0892 -0.6019 -3.6607 0.0666 ...
##  $ V89  : num  0.3684 0.0613 -2.1403 -0.5057 -0.1787 ...
##  $ V90  : num  -1.04054 -0.63697 0.00132 -0.39905 0.55501 ...
##  $ V91  : num  -0.684 1.059 0.405 -0.467 -0.299 ...
##  $ V92  : num  0.282 0.456 0.935 -1.418 -0.931 ...
##  $ V93  : num  0.131 -0.173 0.385 0.217 0.833 ...
##  $ V94  : num  0.112 1.522 -0.787 0.773 -1.872 ...
##  $ V95  : num  -0.00274 -0.9983 -0.80048 0.73148 -0.5522 ...
##  $ V96  : num  -0.317 0.715 -0.348 0.452 0.724 ...
##  $ V97  : num  -0.3684 -0.5157 0.6402 0.0344 0.0505 ...
##  $ V98  : num  -0.63296 -0.00101 -0.93924 -2.16108 -0.59711 ...
##   [list output truncated]
summary(datos$y)
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -13.9094  -3.0931  -0.0759  -0.1668   2.7338  16.8420
hist(datos$y, breaks=30)

1 ¿Hay multicolinealidad en los datos?

Sí, es muy probable que exista multicolinealidad en este conjunto de datos. En primer lugar, tenemos 5000 genes como variables predictoras y solo 1200 observaciones. Esto implica que el número de variables es mucho mayor que el número de datos disponibles. En este escenario, es casi inevitable que existan relaciones lineales fuertes entre varias de las variables.

Además, en el contexto biológico, los genes no actúan de forma independiente. Muchos participan en las mismas rutas metabólicas o procesos celulares, lo que genera correlaciones altas entre ellos. Esto significa que varias columnas de la matriz de diseño pueden estar fuertemente relacionadas.

Desde el punto de vista matemático, cuando el número de variables supera el número de observaciones, la matriz 𝑋 𝑇 𝑋 X T X no es invertible. Esto hace que la estimación por mínimos cuadrados ordinarios no sea única y que los coeficientes sean inestables. Esa es una manifestación clara de multicolinealidad severa.

Por estas razones, se puede concluir que sí existe multicolinealidad en los datos, tanto por la estructura dimensional del problema como por la naturaleza biológica de las variables. Esto justifica el uso de métodos de regularización como ridge y lasso en lugar de regresión lineal clásica.

#### Definimos nuestras X y Y 
y = datos$y

X = as.matrix(datos[,-1])

2

Ahora vamos a realizar la separacion de nuestro conjunto de datos

# Usamos la semilla 49 por elección personal

set.seed(49)  # Fijamos semilla

n <- nrow(datos)

# Seleccionamos los índices para el entrenamiento (1000 filas)
indices_train <- sample(1:n, 1000)

# Creamos los subsets
X_train <- X[indices_train, ]
y_train <- y[indices_train]

X_test <- X[-indices_train, ]
y_test <- y[-indices_train]

Ahora verificamos que las dimensiones sean correctas

dim(X_train)
## [1] 1000 5000
dim(X_test)
## [1]  200 5000

3

set.seed(49)

cv_ridge <- cv.glmnet(
  X_train,
  y_train,
  alpha = 0,           # ridge
  nfolds = 10,
  standardize = TRUE
)

lambda_r <- cv_ridge$lambda.min
set.seed(49)

cv_lasso <- cv.glmnet(
  X_train,
  y_train,
  alpha = 1,           # lasso
  nfolds = 10,
  standardize = TRUE
)

lambda_l <- cv_lasso$lambda.min

Se utilizó validación cruzada de 10 particiones sobre el conjunto de entrenamiento para seleccionar el valor de λ que minimiza el error cuadrático medio (ECM). Este procedimiento permite estimar el desempeño predictivo del modelo en datos no utilizados durante el ajuste y evita utilizar el conjunto de prueba para la selección del hiperparámetro. Se estimaron los valores óptimos de λ para los modelos ridge y lasso utilizando la función cv.glmnet

# Visualización del error en función de lambda
par(mfrow = c(1, 2))
plot(cv_ridge, main = "Ridge (ECM)")
plot(cv_lasso, main = "Lasso (ECM)")

Vamos a validar cuanto es el resultado del lambda y del ECM promedio en validacion cruzada

lambda_r
## [1] 12.44889
lambda_l
## [1] 0.0800224
min(cv_ridge$cvm)
## [1] 17.03078
min(cv_lasso$cvm)
## [1] 1.238432

Para el modelo ridge se obtuvo un valor óptimo de λᵣ = 12.4489, con un ECM mínimo de 17.03078. En el caso del modelo lasso, el valor óptimo fue λₗ = 0.0800224, con un ECM mínimo de 1.238432.

Estos resultados indican que, dentro del conjunto de entrenamiento, el modelo lasso alcanza un error de predicción considerablemente menor que el modelo ridge bajo el esquema de validación cruzada utilizado. Esto sugiere que la penalización L1 podría estar capturando mejor la estructura del problema, posiblemente al realizar selección automática de variables y eliminar genes con poca contribución predictiva.

4

Con los valores óptimos de λ obtenidos mediante validación cruzada en el punto anterior, se ajustaron los modelos de regresión ridge y lasso utilizando las 1000 observaciones del conjunto de entrenamiento. Para ridge se utilizó λᵣ = 12.4489 y para lasso λₗ = 0.0800224. Estos modelos serán utilizados posteriormente para evaluar su capacidad predictiva sobre el conjunto de prueba.

Primero Ridge :

ridge_model <- glmnet(
  X_train,
  y_train,
  alpha = 0,
  lambda = lambda_r,
  standardize = TRUE
)

Segundo Lasso

lasso_model <- glmnet(
  X_train,
  y_train,
  alpha = 1,
  lambda = lambda_l,
  standardize = TRUE
)

5

Para comparar el desempeño predictivo de los modelos ridge y lasso, se calcularon las predicciones sobre el conjunto de prueba de 200 observaciones y se evaluó el error cuadrático medio (ECM). El modelo con menor ECM en el conjunto de prueba se considera el más apropiado para propósitos de predicción.

primero se genera las predicciones

pred_ridge <- predict(ridge_model, newx = X_test)
pred_lasso <- predict(lasso_model, newx = X_test)

Ahora pasamos a calcular el ECM

ecm_ridge <- mean((y_test - pred_ridge)^2)
ecm_lasso <- mean((y_test - pred_lasso)^2)

ecm_ridge
## [1] 17.83095
ecm_lasso
## [1] 1.250956

6 Dado que el modelo lasso fue seleccionado como el más apropiado para predicción en el punto anterior, se ajustó nuevamente este modelo utilizando la totalidad de las 1200 observaciones disponibles. En este paso se mantuvo el valor de λ previamente estimado mediante validación cruzada. El objetivo es obtener una estimación final de los coeficientes utilizando toda la información disponible en los datos.

Primero reconstruimos la matriz completa

X_all <- as.matrix(datos[, -1])
y_all <- datos$y

Y ahora ajustamos el modelo final usando el λ que encontramos previamente

final_model <- glmnet(
  X_all,
  y_all,
  alpha = 1,
  lambda = lambda_l,
  standardize = TRUE
)
sum(coef(final_model) != 0)
## [1] 60

7

Se graficaron las trazas de los coeficientes del modelo lasso en función del parámetro de penalización λ utilizando las 1200 observaciones. Cada curva representa la evolución del coeficiente de un gen a medida que cambia el nivel de penalización. Para valores grandes de λ la mayoría de coeficientes permanecen en cero, mientras que al disminuir λ algunos genes comienzan a ingresar al modelo. La línea azul punteada indica el valor óptimo de λ seleccionado previamente mediante validación cruzada

library(glmnet)

X_all <- as.matrix(datos[, -1])
y_all <- datos$y

lasso_path <- glmnet(
  X_all,
  y_all,
  alpha = 1,
  standardize = TRUE
)
plot(lasso_path, xvar = "lambda", label = FALSE)

8

El objetivo del estudio fue identificar qué genes son relevantes para predecir la efectividad de un tratamiento anticancerígeno a partir de un conjunto de datos con 1200 líneas celulares y 5000 genes. Dado que el número de variables es mucho mayor que el número de observaciones, es esperable la presencia de multicolinealidad, por lo que se utilizaron métodos de regresión penalizada. Primero se dividió el conjunto de datos en entrenamiento con 1000 observaciones y prueba con 200 observaciones. Sobre el conjunto de entrenamiento se aplicó validación cruzada de diez particiones para estimar el valor óptimo del parámetro de regularización lambda en los modelos ridge y lasso, utilizando el error cuadrático medio como criterio de selección. Posteriormente se ajustaron ambos modelos y se evaluó su capacidad predictiva en el conjunto de prueba. Los resultados mostraron que el modelo lasso obtuvo un error cuadrático medio de alrededor de 1.25, mientras que el modelo ridge presentó un error de alrededor de 17.83. Esto indica que el modelo lasso ofrece un desempeño predictivo considerablemente mejor en este conjunto de datos. Este resultado sugiere que la penalización L1 es más adecuada en este problema, ya que permite realizar selección automática de variables y conservar únicamente los genes con mayor capacidad predictiva. Finalmente, el modelo lasso seleccionado se ajustó utilizando las 1200 observaciones disponibles y se analizaron las trazas de los coeficientes para observar cómo los genes entran al modelo a medida que disminuye el nivel de penalización. En conjunto, los resultados indican que solo un subconjunto de los genes disponibles contribuye de forma relevante a la predicción de la efectividad del tratamiento.