Actividad Guiada: Selección de Modelos y Regularización

Datos de Cáncer de Próstata (Elements of Statistical Learning)

Actuaría. Regresión

Introducción

Esta actividad retoma los datos de cáncer de próstata que ya trabajamos en clase y los lleva un paso más allá. Resolverás, en orden, las siguientes tareas:

1. Estandarizar las variables numéricas y preparar los factores (tú mismo/a, sin ver antes la solución).
2. Ajustar el modelo de regresión lineal completo.
3. Repetir la selección de mejores subconjuntos, pero ahora usando todas las variables disponibles (nvmax dinámico) y un criterio objetivo (no solo inspección visual del gráfico de RSS).
4. Ajustar una regresión Ridge con \(\lambda\) seleccionado por validación cruzada.
5. Ajustar una regresión Lasso con \(\lambda\) seleccionado por validación cruzada.
6. Comparar qué variables selecciona cada método.
7. Ajustar Regresión de Componentes Principales (PCR), seleccionando el número de componentes por validación cruzada.
8. Ajustar Mínimos Cuadrados Parciales (PLS), seleccionando el número de componentes por validación cruzada.
9. Comparar el error de predicción de todos los métodos en el conjunto de prueba.

Cada sección tiene un bloque de código con una plantilla (# TODO) para que lo completes. Justo después encontrarás un bloque colapsable “Solución” — intenta resolver el ejercicio antes de abrirlo.

Paquetes

library(tidyverse)
library(leaps)
library(glmnet)
library(pls)
library(knitr)

Paso 0: Carga de datos (dado)

Note

El paquete ElemStatLearn, que originalmente distribuía este conjunto de datos, ya no está disponible en CRAN. Por ello, leemos los datos directamente desde la fuente original en el sitio de Stanford, lo cual además garantiza reproducibilidad total del ejercicio.
Los datos se distribuyen como parte del libro ``The Elements of Statistical Learning’’ disponible en https://hastie.su.domains/ElemStatLearn/

Se refieren a un estudio originalmente publicado como:

Stamey, T.A., Kabalin, J.N., McNeal, J.E., Johnstone, I.M., Freiha, F., Redwine, E.A., & Yang, N. (1989). Prostate specific antigen in the diagnosis and treatment of adenocarcinoma of the prostate II: radical prostatectomy treated patients. Journal of Urology, 16(5), 1076–1083.

prostate <- read_tsv(
  "https://hastie.su.domains/ElemStatLearn/datasets/prostate.data",
  col_select = -1
)

glimpse(prostate)

Rows: 97
Columns: 10
$ lcavol  <dbl> -0.5798185, -0.9942523, -0.5108256, -1.2039728, 0.7514161, -1.…
$ lweight <dbl> 2.769459, 3.319626, 2.691243, 3.282789, 3.432373, 3.228826, 3.…
$ age     <dbl> 50, 58, 74, 58, 62, 50, 64, 58, 47, 63, 65, 63, 63, 67, 57, 66…
$ lbph    <dbl> -1.3862944, -1.3862944, -1.3862944, -1.3862944, -1.3862944, -1…
$ svi     <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
$ lcp     <dbl> -1.3862944, -1.3862944, -1.3862944, -1.3862944, -1.3862944, -1…
$ gleason <dbl> 6, 6, 7, 6, 6, 6, 6, 6, 6, 6, 6, 6, 7, 7, 7, 6, 7, 6, 6, 6, 6,…
$ pgg45   <dbl> 0, 0, 20, 0, 0, 0, 0, 0, 0, 0, 0, 0, 30, 5, 5, 0, 30, 0, 0, 0,…
$ lpsa    <dbl> -0.4307829, -0.1625189, -0.1625189, -0.1625189, 0.3715636, 0.7…
$ train   <lgl> TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, FALSE, TRUE, FALSE, FALSE,…

dim(prostate)

[1] 97 10

Las columnas son las mismas que vimos en clase: lcavol, lweight, age, lbph, svi, lcp, gleason, pgg45 como predictores; lpsa como variable dependiente; y train como indicador lógico de partición entrenamiento/prueba. Lee el abstract del paper e identifica con ello las variables de la base de datos.

Abstract

Serum prostate specific antigen was determined (Yang polyclonal radioimmunoassay) in 102 men before hospitalization for radical prostatectomy. Prostate specimens were subjected to detailed histological and morphometric analysis. Levels of prostate specific antigen were significantly different between patients with and without a Gleason score of 7 or greater (p less than 0.001), capsular penetration greater than 1 cm. in linear extent (p less than 0.001), seminal vesicle invasion (p less than 0.001) and pelvic lymph node metastasis (p less than 0.005). Prostate specific antigen was strongly correlated with volume of prostate cancer (r equals 0.70). Bivariate and multivariate analyses indicate that cancer volume is the primary determinant of serum prostate specific antigen levels. Prostate specific antigen was elevated 3.5 ng. per ml. for every cc of cancer, a level at least 10 times that observed for benign prostatic hyperplasia. Prostate specific antigen is useful as a preoperative marker because no patient with lymph node metastasis had serum levels of less than 10 ng. per ml. (4 times the upper limit of normal range). Of the patients with greater than 50 ng. per ml. two-thirds had microscopic lymph node metastasis and 90 per cent had seminal vesicle invasion. Serum prostatic acid phosphatase levels showed a significantly weaker correlation with cancer volume (r equals 0.51) and every other pathological parameter. Of the patients 73 per cent had serum prostatic acid phosphatase levels in the normal range (0 to 2.1 ng. per ml.), including 7 per cent who had pelvic lymph node metastasis. Postoperative prostate specific antigen values were available in 97 of 102 patients, with a mean and maximum followup of 12 and 38 months. No patient with pelvic lymph node metastasis achieved an undetectable prostate specific antigen level without adjunctive therapy (hormonal or radiation). No difference in preoperative or postoperative prostate specific antigen levels, cancer volume, seminal vesicle invasion or incidence of pelvic lymph node metastasis was seen between patients with no capsular penetration and those with minimal capsular penetration (1 cm. or less total linear extent of full thickness penetration), providing the first quantitative evidence that small amounts of capsular penetration may not be of biological or prognostic significance.

Paso 1: Estandarización — Tu turno

Instrucciones: crea un objeto df_scaled a partir de prostate que:

• Convierta svi en factor.
• Convierta gleason en factor ordenado.
• Estandarice (resta la media, divide entre la desviación estándar) todas las variables numéricas, excepto la variable dependiente lpsa (y, por supuesto, sin tocar train).

y_variable_name <- "lpsa"

df_scaled <- prostate |>
  mutate(
    svi = ___,
    gleason = ___
  ) |>
  mutate(
    across(
      ___,          # ¿qué columnas? (numéricas y distintas de lpsa)
      ___           # ¿qué función de estandarización?
    )
  )

df_scaled |>
  select(-c(train, lpsa)) |>
  summary()

TipSolución

y_variable_name <- "lpsa"

df_scaled <- prostate |>
  mutate(
    svi = factor(svi),
    gleason = factor(gleason, ordered = TRUE)
  ) |>
  mutate(
    across(
      where(is.numeric) & -all_of(y_variable_name),
      \(x) (x - mean(x)) / sd(x)
    )
  )

df_scaled |>
  select(-c(train, lpsa)) |>
  summary()

     lcavol            lweight              age               lbph        
 Min.   :-2.28833   Min.   :-2.92718   Min.   :-3.0713   Min.   :-1.0247  
 1st Qu.:-0.71031   1st Qu.:-0.59070   1st Qu.:-0.5193   1st Qu.:-1.0247  
 Median : 0.08222   Median :-0.01386   Median : 0.1523   Median : 0.1377  
 Mean   : 0.00000   Mean   : 0.00000   Mean   : 0.0000   Mean   : 0.0000  
 3rd Qu.: 0.65927   3rd Qu.: 0.57761   3rd Qu.: 0.5553   3rd Qu.: 1.0048  
 Max.   : 2.09651   Max.   : 2.68770   Max.   : 2.0327   Max.   : 1.5343  
 svi         lcp          gleason     pgg45        
 0:76   Min.   :-0.8632   6:35    Min.   :-0.8645  
 1:21   1st Qu.:-0.8632   7:56    1st Qu.:-0.8645  
        Median :-0.4428   8: 1    Median :-0.3326  
        Mean   : 0.0000   9: 5    Mean   : 0.0000  
        3rd Qu.: 0.9712           3rd Qu.: 0.5538  
        Max.   : 2.2053           Max.   : 2.6811  

Paso 2: Modelo completo — Tu turno

Instrucciones: ajusta lpsa contra todas las demás variables (excepto train), usando únicamente las observaciones de entrenamiento (train == TRUE). Muestra el summary().

full_model <- lm(
  ___,
  data = ___,
  subset = ___
)
summary(full_model)

TipSolución

full_model <- lm(
  lpsa ~ . - train,
  data = df_scaled,
  subset = train == TRUE
)
summary(full_model)

Call:
lm(formula = lpsa ~ . - train, data = df_scaled, subset = train == 
    TRUE)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.63386 -0.29322 -0.08569  0.44765  1.54455 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  2.32130    0.22602  10.270 1.72e-14 ***
lcavol       0.66959    0.12967   5.164 3.32e-06 ***
lweight      0.27682    0.09569   2.893  0.00543 ** 
age         -0.16515    0.10105  -1.634  0.10780    
lbph         0.19140    0.10369   1.846  0.07020 .  
svi1         0.70716    0.30509   2.318  0.02413 *  
lcp         -0.34327    0.16812  -2.042  0.04590 *  
gleason.L   -0.21194    0.45168  -0.469  0.64074    
gleason.Q   -0.70358    0.46424  -1.516  0.13526    
gleason.C   -0.47598    0.58087  -0.819  0.41602    
pgg45        0.31176    0.17189   1.814  0.07509 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.7031 on 56 degrees of freedom
Multiple R-squared:  0.7124,    Adjusted R-squared:  0.6611 
F-statistic: 13.87 on 10 and 56 DF,  p-value: 6.847e-12

Paso 3: Selección de mejor subconjunto, con todas las variables — Tu turno

En clase elegimos visualmente \(p = 3\) variables viendo el gráfico de RSS. Ahora harás algo más riguroso:

Instrucciones:

1. Calcula cuántas columnas tendría la matriz de diseño con todas las variables (usa model.matrix() sobre la fórmula completa y réstale 1, por el intercepto). Ese número es tu nvmax.
2. Ajusta regsubsets() con ese nvmax (en vez de un valor fijo pequeño).
3. Extrae summary() del objeto y grafica tres criterios en un panel: \(C_p\), BIC y \(R^2\) ajustado, marcando en cada uno el tamaño óptimo según ese criterio (which.min o which.max, según corresponda).
4. Decide con cuántas y cuáles variables te quedas (puedes usar el criterio BIC, que suele ser el más conservador) y reajusta el modelo final con lm().

n_predictors <- ncol(model.matrix(___, data = df_scaled)) - ___

best_subsets <- regsubsets(
  lpsa ~ .,
  data = df_scaled,
  subset = train == TRUE,
  nvmax = ___,
  method = "exhaustive"
)

subsets_summary <- summary(best_subsets)

par(mfrow = c(1, 3))
# Cp
plot(subsets_summary$cp, xlab = "Número de variables", ylab = "Cp", type = "b")
points(___, subsets_summary$cp[___], col = "red", pch = 19)

# BIC
plot(subsets_summary$bic, xlab = "Número de variables", ylab = "BIC", type = "b")
points(___, subsets_summary$bic[___], col = "red", pch = 19)

# Adjusted R^2
plot(subsets_summary$adjr2, xlab = "Número de variables", ylab = "Adj. R2", type = "b")
points(___, subsets_summary$adjr2[___], col = "red", pch = 19)
par(mfrow = c(1, 1))

# ¿Qué variables selecciona el criterio que elegiste?
coef(best_subsets, id = ___)

TipSolución

n_predictors <- ncol(model.matrix(lpsa ~ . - train, data = df_scaled)) - 1
n_predictors

[1] 10

best_subsets <- regsubsets(
  lpsa ~ .,
  data = df_scaled,
  subset = train == TRUE,
  nvmax = n_predictors,
  method = "exhaustive"
)

subsets_summary <- summary(best_subsets)
subsets_summary

Subset selection object
Call: regsubsets.formula(lpsa ~ ., data = df_scaled, subset = train == 
    TRUE, nvmax = n_predictors, method = "exhaustive")
11 Variables  (and intercept)
          Forced in Forced out
lcavol        FALSE      FALSE
lweight       FALSE      FALSE
age           FALSE      FALSE
lbph          FALSE      FALSE
svi1          FALSE      FALSE
lcp           FALSE      FALSE
gleason.L     FALSE      FALSE
gleason.Q     FALSE      FALSE
gleason.C     FALSE      FALSE
pgg45         FALSE      FALSE
trainTRUE     FALSE      FALSE
1 subsets of each size up to 10
Selection Algorithm: exhaustive
          lcavol lweight age lbph svi1 lcp gleason.L gleason.Q gleason.C pgg45
1  ( 1 )  "*"    " "     " " " "  " "  " " " "       " "       " "       " "  
2  ( 1 )  "*"    "*"     " " " "  " "  " " " "       " "       " "       " "  
3  ( 1 )  "*"    "*"     " " " "  " "  " " " "       " "       "*"       " "  
4  ( 1 )  "*"    "*"     " " " "  "*"  " " " "       " "       "*"       " "  
5  ( 1 )  "*"    "*"     " " "*"  "*"  " " " "       "*"       " "       " "  
6  ( 1 )  "*"    "*"     "*" "*"  "*"  " " " "       "*"       " "       " "  
7  ( 1 )  "*"    "*"     "*" "*"  "*"  "*" " "       " "       " "       "*"  
8  ( 1 )  "*"    "*"     "*" "*"  "*"  "*" " "       "*"       " "       "*"  
9  ( 1 )  "*"    "*"     "*" "*"  "*"  "*" " "       "*"       "*"       "*"  
10  ( 1 ) "*"    "*"     "*" "*"  "*"  "*" "*"       "*"       "*"       "*"  
          trainTRUE
1  ( 1 )  " "      
2  ( 1 )  " "      
3  ( 1 )  " "      
4  ( 1 )  " "      
5  ( 1 )  " "      
6  ( 1 )  " "      
7  ( 1 )  " "      
8  ( 1 )  " "      
9  ( 1 )  " "      
10  ( 1 ) " "      

best_cp    <- which.min(subsets_summary$cp)
best_bic   <- which.min(subsets_summary$bic)
best_adjr2 <- which.max(subsets_summary$adjr2)

c(best_cp = best_cp, best_bic = best_bic, best_adjr2 = best_adjr2)

   best_cp   best_bic best_adjr2 
         8          3          8 

par(mfrow = c(1, 3))
plot(subsets_summary$cp, xlab = "Número de variables", ylab = "Cp", type = "b")
points(best_cp, subsets_summary$cp[best_cp], col = "red", pch = 19)

plot(subsets_summary$bic, xlab = "Número de variables", ylab = "BIC", type = "b")
points(best_bic, subsets_summary$bic[best_bic], col = "red", pch = 19)

plot(subsets_summary$adjr2, xlab = "Número de variables", ylab = "Adj. R2", type = "b")
points(best_adjr2, subsets_summary$adjr2[best_adjr2], col = "red", pch = 19)

par(mfrow = c(1, 1))

Usamos el criterio BIC (el más parsimonioso) para elegir el modelo final y lo reajustamos con lm():

# regsubsets() expands factors into dummy/contrast columns -- our ordered
# factor `gleason` becomes "gleason.L", "gleason.Q", "gleason.C" (linear,
# quadratic, cubic contrasts). Those raw coefficient names are NOT valid
# variable names for reformulate()/lm(): R would look for a column
# literally called "gleason.C", which doesn't exist in df_scaled (only
# "gleason" does). We map each selected coefficient back to its original
# predictor name before rebuilding the formula.

coef_names <- names(coef(best_subsets, id = best_bic))[-1]
coef_names

[1] "lcavol"    "lweight"   "gleason.C"

predictor_names <- setdiff(names(df_scaled), c("lpsa", "train"))

selected_vars <- predictor_names[
  sapply(predictor_names, function(v) any(startsWith(coef_names, v)))
]
selected_vars

[1] "lcavol"  "lweight" "gleason"

subset_formula <- reformulate(selected_vars, response = "lpsa")
subset_formula

lpsa ~ lcavol + lweight + gleason

subset_final_model <- lm(subset_formula, data = df_scaled, subset = train == TRUE)
summary(subset_final_model)

Call:
lm(formula = subset_formula, data = df_scaled, subset = train == 
    TRUE)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.53127 -0.45955 -0.04718  0.51205  1.85745 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  2.30871    0.21984  10.502 2.65e-15 ***
lcavol       0.62368    0.10637   5.863 1.98e-07 ***
lweight      0.33344    0.08687   3.838 0.000297 ***
gleason.L   -0.07376    0.35861  -0.206 0.837717    
gleason.Q   -0.17579    0.44028  -0.399 0.691091    
gleason.C    0.42212    0.51757   0.816 0.417914    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.7426 on 61 degrees of freedom
Multiple R-squared:  0.6506,    Adjusted R-squared:  0.622 
F-statistic: 22.72 on 5 and 61 DF,  p-value: 8.533e-13

Paso 4: Regresión Ridge con CV — Tu turno

Instrucciones: construye las matrices X y y de entrenamiento (recuerda quitar la columna de intercepto con [,-1]), selecciona \(\lambda\) por validación cruzada con alpha = 0, extrae los coeficientes en lambda.min, y grafica la trayectoria completa de los coeficientes en función de \(\log(\lambda)\).

X_train <- model.matrix(
  ___,
  data = filter(df_scaled, train == TRUE)
)[, ___]

y_train <- df_scaled |>
  filter(train == TRUE) |>
  pull(___)

set.seed(123)
cv_ridge <- cv.glmnet(X_train, y_train, alpha = ___)
plot(cv_ridge)

best_lambda_ridge <- cv_ridge$___

ridge_coeffs <- coef(cv_ridge, s = "lambda.min")
print(ridge_coeffs)

ridge_path <- glmnet(X_train, y_train, alpha = ___)
plot(ridge_path, xvar = "lambda", label = TRUE)

TipSolución

X_train <- model.matrix(
  lpsa ~ . - train,
  data = filter(df_scaled, train == TRUE)
)[, -1]

y_train <- df_scaled |>
  filter(train == TRUE) |>
  pull(lpsa)

set.seed(123)
cv_ridge <- cv.glmnet(X_train, y_train, alpha = 0)
plot(cv_ridge)

best_lambda_ridge <- cv_ridge$lambda.min
best_lambda_ridge

[1] 0.08788804

ridge_coeffs <- coef(cv_ridge, s = "lambda.min")
print(ridge_coeffs)

11 x 1 sparse Matrix of class "dgCMatrix"
             lambda.min
(Intercept)  2.27316445
lcavol       0.55042243
lweight      0.27532301
age         -0.12447206
lbph         0.17266600
svi1         0.60771023
lcp         -0.16359728
gleason.L   -0.08937738
gleason.Q   -0.35171726
gleason.C    0.01567191
pgg45        0.17684936

ridge_path <- glmnet(X_train, y_train, alpha = 0)
plot(ridge_path, xvar = "lambda", label = TRUE)
abline(v = log(best_lambda_ridge), lty = 2)

Paso 5: Regresión Lasso con CV — Tu turno

Instrucciones: repite exactamente el mismo procedimiento del paso anterior, cambiando únicamente alpha al valor correspondiente a Lasso. Observa cuáles coeficientes quedan en cero exacto.

set.seed(123)
cv_lasso <- cv.glmnet(X_train, y_train, alpha = ___)
plot(cv_lasso)

best_lambda_lasso <- cv_lasso$___

lasso_coeffs <- coef(cv_lasso, s = "lambda.min")
print(lasso_coeffs)

lasso_path <- glmnet(X_train, y_train, alpha = ___)
plot(lasso_path, xvar = "lambda", label = TRUE)

TipSolución

set.seed(123)
cv_lasso <- cv.glmnet(X_train, y_train, alpha = 1)
plot(cv_lasso)

best_lambda_lasso <- cv_lasso$lambda.min
best_lambda_lasso

[1] 0.01109021

lasso_coeffs <- coef(cv_lasso, s = "lambda.min")
print(lasso_coeffs)

11 x 1 sparse Matrix of class "dgCMatrix"
             lambda.min
(Intercept)  2.29425705
lcavol       0.61252661
lweight      0.27955187
age         -0.13508483
lbph         0.16835536
svi1         0.61249772
lcp         -0.21262216
gleason.L   -0.04636499
gleason.Q   -0.34077176
gleason.C    .         
pgg45        0.17796172

lasso_path <- glmnet(X_train, y_train, alpha = 1)
plot(lasso_path, xvar = "lambda", label = TRUE)
abline(v = log(best_lambda_lasso), lty = 2)

Paso 6: ¿Qué variables selecciona cada método? — Tu turno

Instrucciones: construye una tabla (una fila por variable, una columna por método) que indique si cada variable fue conservada por: (a) el modelo de mejor subconjunto del Paso 3, (b) Ridge (en la práctica, casi siempre todas, ya que Ridge no hace selección exacta), y (c) Lasso (coeficiente distinto de cero en lambda.min).

# Dado: igual que en el Paso 3, los coeficientes de Ridge/Lasso vienen
# nombrados por columna dummy/contraste (p.ej. "gleason.L"), no por
# variable original. Esta función mapea cada uno de vuelta a su variable.
predictor_names <- setdiff(names(df_scaled), c("lpsa", "train"))

base_var_of <- function(coef_name) {
  predictor_names[sapply(predictor_names, \(v) startsWith(coef_name, v))][1]
}

ridge_names   <- rownames(ridge_coeffs)[-1]
ridge_nonzero <- as.vector(ridge_coeffs)[-1] != 0
lasso_names   <- rownames(lasso_coeffs)[-1]
lasso_nonzero <- as.vector(lasso_coeffs)[-1] != 0

ridge_base <- sapply(ridge_names, base_var_of)
lasso_base <- sapply(lasso_names, base_var_of)

comparison_table <- tibble(
  variable = predictor_names,
  subconjuntos = predictor_names %in% ___,
  ridge = predictor_names %in% ___[___],
  lasso = predictor_names %in% ___[___]
)

comparison_table |> kable()

TipSolución

predictor_names <- setdiff(names(df_scaled), c("lpsa", "train"))

base_var_of <- function(coef_name) {
  predictor_names[sapply(predictor_names, \(v) startsWith(coef_name, v))][1]
}

ridge_names   <- rownames(ridge_coeffs)[-1]
ridge_nonzero <- as.vector(ridge_coeffs)[-1] != 0
lasso_names   <- rownames(lasso_coeffs)[-1]
lasso_nonzero <- as.vector(lasso_coeffs)[-1] != 0

ridge_base <- sapply(ridge_names, base_var_of)
lasso_base <- sapply(lasso_names, base_var_of)

comparison_table <- tibble(
  variable = predictor_names,
  subconjuntos = predictor_names %in% selected_vars,
  ridge = predictor_names %in% ridge_base[ridge_nonzero],
  lasso = predictor_names %in% lasso_base[lasso_nonzero]
)

comparison_table |> kable()

variable	subconjuntos	ridge	lasso
lcavol	TRUE	TRUE	TRUE
lweight	TRUE	TRUE	TRUE
age	FALSE	TRUE	TRUE
lbph	FALSE	TRUE	TRUE
svi	FALSE	TRUE	TRUE
lcp	FALSE	TRUE	TRUE
gleason	TRUE	TRUE	TRUE
pgg45	FALSE	TRUE	TRUE

Observa que Ridge casi nunca produce ceros exactos (reduce los coeficientes, pero no los elimina), mientras que Lasso sí puede coincidir con el subconjunto de variables que eligió la selección de subconjuntos — aunque no necesariamente son idénticos, porque ambos métodos penalizan de forma distinta.

Paso 7: Regresión de Componentes Principales (PCR) — Tu turno

Instrucciones: ajusta un modelo PCR sobre los datos de entrenamiento con pls::pcr(), usando validation = "CV" (validación cruzada interna). Grafica el error de validación (validationplot(..., val.type = "MSEP")) y elige el número de componentes que minimiza el error (puedes usar selectNcomp()).

set.seed(123)
pcr_model <- pcr(
  ___,
  data = df_scaled,
  subset = train == TRUE,
  scale = ___,          # ¿ya estandarizamos los datos? ¿hace falta volver a escalar?
  validation = ___
)

validationplot(pcr_model, val.type = "MSEP")

ncomp_pcr <- selectNcomp(pcr_model, method = "onesigma", plot = TRUE)
ncomp_pcr

summary(pcr_model)

TipSolución

set.seed(123)
pcr_model <- pcr(
  lpsa ~ . - train,
  data = df_scaled,
  subset = train == TRUE,
  scale = FALSE,     # las variables numéricas ya están estandarizadas desde el Paso 1
  validation = "CV"
)

validationplot(pcr_model, val.type = "MSEP")

ncomp_pcr <- selectNcomp(pcr_model, method = "onesigma", plot = TRUE)

ncomp_pcr

[1] 2

summary(pcr_model)

Data:   X dimension: 67 10 
    Y dimension: 67 1
Fit method: svdpc
Number of components considered: 10

VALIDATION: RMSEP
Cross-validated using 10 random segments.
       (Intercept)  1 comps  2 comps  3 comps  4 comps  5 comps  6 comps
CV           1.217   0.8588   0.8450   0.8206   0.8122   0.8133   0.7595
adjCV        1.217   0.8557   0.8426   0.8149   0.8094   0.8094   0.7545
       7 comps  8 comps  9 comps   10 comps
CV      0.7585   0.7490   0.7749  1.215e+14
adjCV   0.7543   0.7437   0.7683  1.150e+14

TRAINING: % variance explained
      1 comps  2 comps  3 comps  4 comps  5 comps  6 comps  7 comps  8 comps
X       40.66    64.46    75.89    84.44    91.32    95.37    98.20    99.33
lpsa    52.16    54.00    59.86    60.75    62.26    68.01    68.25    70.32
      9 comps  10 comps
X       99.81    100.00
lpsa    70.36     71.24

Paso 8: Mínimos Cuadrados Parciales (PLS) — Tu turno

Instrucciones: repite el procedimiento anterior usando pls::plsr() en vez de pcr(). Todo lo demás es análogo.

set.seed(123)
pls_model <- plsr(
  ___,
  data = df_scaled,
  subset = train == TRUE,
  scale = ___,
  validation = ___
)

validationplot(pls_model, val.type = "MSEP")

ncomp_pls <- selectNcomp(pls_model, method = "onesigma", plot = TRUE)
ncomp_pls

summary(pls_model)

TipSolución

set.seed(123)
pls_model <- plsr(
  lpsa ~ . - train,
  data = df_scaled,
  subset = train == TRUE,
  scale = FALSE,
  validation = "CV"
)

validationplot(pls_model, val.type = "MSEP")

ncomp_pls <- selectNcomp(pls_model, method = "onesigma", plot = TRUE)

ncomp_pls

[1] 1

summary(pls_model)

Data:   X dimension: 67 10 
    Y dimension: 67 1
Fit method: kernelpls
Number of components considered: 10

VALIDATION: RMSEP
Cross-validated using 10 random segments.
       (Intercept)  1 comps  2 comps  3 comps  4 comps  5 comps  6 comps
CV           1.217   0.8119   0.7746   0.7649   0.7562   0.7572   0.7615
adjCV        1.217   0.8100   0.7705   0.7610   0.7515   0.7522   0.7556
       7 comps  8 comps  9 comps   10 comps
CV      0.7618   0.7774   0.7814  9.166e+10
adjCV   0.7561   0.7702   0.7734  8.674e+10

TRAINING: % variance explained
      1 comps  2 comps  3 comps  4 comps  5 comps  6 comps  7 comps  8 comps
X       40.11    53.76    71.29    79.89    85.61    89.28     96.1    98.13
lpsa    58.69    66.29    68.33    69.69    70.25    70.72     70.8    71.03
      9 comps  10 comps
X       99.53    100.00
lpsa    71.24     71.24

Note

Nota pedagógica: PLS suele necesitar menos componentes que PCR para alcanzar un error de validación comparable, porque construye las componentes usando información de lpsa, no solo de la varianza de los predictores (a diferencia de PCR, que es completamente no supervisado en la construcción de componentes).

Paso 9: Comparación final en el conjunto de prueba — Tu turno

Esta es la pregunta que de verdad importa: ¿qué tan bien predice cada método en datos que no vio durante el entrenamiento?

Instrucciones: usando las observaciones con train == FALSE, calcula el error cuadrático medio (MSE) de predicción para cada uno de los seis modelos que ajustamos: modelo completo, mejor subconjunto, Ridge, Lasso, PCR y PLS. Preséntalos en una tabla ordenada de menor a mayor error, y en una gráfica de barras.

test_data <- df_scaled |> filter(train == FALSE)

X_test <- model.matrix(___, data = test_data)[, -1]
y_test <- test_data |> pull(___)

pred_full     <- predict(___, newdata = ___)
pred_subset   <- predict(___, newdata = ___)
pred_ridge    <- predict(cv_ridge, newx = ___, s = "lambda.min")
pred_lasso    <- predict(cv_lasso, newx = ___, s = "lambda.min")
pred_pcr      <- predict(pcr_model, newdata = ___, ncomp = ___)
pred_pls      <- predict(pls_model, newdata = ___, ncomp = ___)

mse <- function(pred, actual) mean((actual - pred)^2)

results <- tibble(
  metodo = c("OLS completo", "Mejor subconjunto", "Ridge", "Lasso", "PCR", "PLS"),
  test_mse = c(
    mse(pred_full, y_test),
    mse(pred_subset, y_test),
    mse(pred_ridge, y_test),
    mse(pred_lasso, y_test),
    mse(pred_pcr, y_test),
    mse(pred_pls, y_test)
  )
) |> arrange(test_mse)

results |> kable(digits = 4)

results |>
  ggplot(aes(x = reorder(metodo, test_mse), y = test_mse)) +
  geom_col(fill = "#2c3e50") +
  coord_flip() +
  labs(
    title = "Error de predicción (MSE) en el conjunto de prueba",
    x = NULL, y = "MSE de prueba"
  ) +
  theme_minimal()

TipSolución

test_data <- df_scaled |> filter(train == FALSE)

X_test <- model.matrix(lpsa ~ . - train, data = test_data)[, -1]
y_test <- test_data |> pull(lpsa)

pred_full   <- predict(full_model, newdata = test_data)
pred_subset <- predict(subset_final_model, newdata = test_data)
pred_ridge  <- predict(cv_ridge, newx = X_test, s = "lambda.min")
pred_lasso  <- predict(cv_lasso, newx = X_test, s = "lambda.min")
pred_pcr    <- predict(pcr_model, newdata = test_data, ncomp = ncomp_pcr)
pred_pls    <- predict(pls_model, newdata = test_data, ncomp = ncomp_pls)

mse <- function(pred, actual) mean((actual - pred)^2)

results <- tibble(
  metodo = c("OLS completo", "Mejor subconjunto", "Ridge", "Lasso", "PCR", "PLS"),
  test_mse = c(
    mse(pred_full, y_test),
    mse(pred_subset, y_test),
    mse(pred_ridge, y_test),
    mse(pred_lasso, y_test),
    mse(pred_pcr, y_test),
    mse(pred_pls, y_test)
  )
) |> arrange(test_mse)

results |> kable(digits = 4)

metodo	test_mse
Lasso	0.4944
Ridge	0.4988
Mejor subconjunto	0.5369
OLS completo	0.5512
PLS	0.5797
PCR	0.7004

results |>
  ggplot(aes(x = reorder(metodo, test_mse), y = test_mse)) +
  geom_col(fill = "#2c3e50") +
  coord_flip() +
  labs(
    title = "Error de predicción (MSE) en el conjunto de prueba",
    x = NULL, y = "MSE de prueba"
  ) +
  theme_minimal()

Conclusión

Responde brevemente (2–3 líneas cada una):

1. ¿Cuál método tuvo el menor error de prueba? ¿Coincide con el que hubieras esperado antes de hacer el ejercicio?
2. ¿Hubo acuerdo entre las variables seleccionadas por mejor subconjunto y por Lasso? ¿Qué nos dice eso sobre la estabilidad de la selección de variables?
3. PCR y PLS suelen necesitar distinto número de componentes para un desempeño similar. ¿Por qué ocurre esto, dado cómo construye cada método sus componentes?
4. Con un conjunto de datos tan pequeño (menos de 100 observaciones), ¿qué tan confiable te parece una sola partición entrenamiento/prueba para comparar métodos? ¿Qué alternativa propondrías?