Los básicos

En R tenemos los comandos set.seed(...) y sample(...) * set.seed(...) nos permite definir una semilla para el generador de números aleatorios

  • sample(...) nos permite tomar muestras aleatorias (muestrear con probabilidades predefinidas) ?sample

  • También, podemos usar el paquete ISLR que implementa algunos flujos que facilitan el uso de modelos de aprendizaje máquina en R

install.packages('ISLR')

Train, test split y lm

Usaremos los datos de autos en la base mpg que son parte de la paquetería tidyverse

library(ISLR2)
library(tidyverse)

mpg = read.table('data/Auto.data', header = TRUE)
mpg <- mpg |> 
          mutate(horsepower = as.numeric(horsepower)) |> 
          select(c(mpg, horsepower)) |> 
          na.omit()
glimpse(mpg)
set.seed(1)
n_obs = nrow(mpg)
1train = sample(n_obs, size = 196)
1
Aquí hacemos el muestreo de índices, aunque pudimos hacer sample(mpg, size = 196). Esto nos daría una base de datos completa con 196 observaciones seleccionadas al azar. Pero hacer el muestreo por índices facilita el trabajo con la función lm 
Rows: 392
Columns: 2
$ mpg        <dbl> 18, 15, 18, 16, 17, 15, 14, 14, 14, 15, 15, 14, 15, 14, 24,…
$ horsepower <dbl> 130, 165, 150, 150, 140, 198, 220, 215, 225, 190, 170, 160,…

lm_model <- lm(mpg ~ horsepower, data = mpg)
1all_error <- mpg$mpg - predict(lm_model, mpg)
2test_error <- all_error[-train]
3test_mse <- mean((test_error)^2)
print(paste('Validation MSE: ', test_mse)) # con una partición aleatoria, pero fija, de la muestra
1
Aplicamos predict con datos mpg para obtener \(\hat{y}\) para cada observación presente en la muestra
2
Calculamos el error total y seleccionamos omitiendo los índices de entrenamiento: -train
3
El argumento na.rm=TRUE evita que un NA presente en los datos genere NA por respuesta 
[1] "Validation MSE:  22.6119681573627"

Podemos escribir una función para calcular el MSE de validación con objetos lm:

mse_validation <- function(model, all_data, train_indexes){
  formula = model$call$formula 
1  y_variable = all.vars(formula)[1]
2  all_error = all_data[[y_variable]] - predict(model, all_data)
  test_error = all_error[-train]
  test_mse <- mean(test_error^2, na.rm = TRUE)
  return(test_mse)
}
1
Estas primeras líneas obtienen el nombre de la variable dependient
2
El uso de all_data[[y_variable]] y cómo R base entiende el subsetting…

Ya armados de esta función podemos analizar distintos modelos polinomiales más fácilmente

for (degree in 1:3){
2  lm_degree_fit <- lm(mpg ~ poly(horsepower, degree), data = mpg, subset = train)
  mse_degree <- mse_validation(lm_degree_fit, mpg)
  print(paste('Grado', degree, 'MSE: ', mse_degree))
}
2
poly(var, degree) especifica un polinomio de grado degree en la variable var 
[1] "Grado 1 MSE:  23.2660086465003"
[1] "Grado 2 MSE:  18.7164594933828"
[1] "Grado 3 MSE:  18.7940067973945"

Cambiando los datos de referencia cambiamos el error

set.seed(2)
train <- sample(n_obs, size = 196)
for (degree in 1:3){
  lm_degree_fit <- lm(mpg ~ poly(horsepower, degree), data = mpg, subset = train)
  mse_degree <- mse_validation(lm_degree_fit, mpg)
  print(paste('Grado', degree, 'MSE: ', mse_degree))
}
[1] "Grado 1 MSE:  25.7265106448139"
[1] "Grado 2 MSE:  20.4303642741464"
[1] "Grado 3 MSE:  20.3853268638776"

Obsérvese que se usó set.seed(2) cambiando los datos de entrenamiento

LOOCV

La función glm también puede ajustar una regresión lineal común si se deja el argumento family en su valor por defecto. Se puede usar la función cv.glm() del paquete boot para realizar cálculos de validación cruzada

library(boot)
lm_fit <- glm(mpg ~  horsepower, data = mpg)
cv_error <- cv.glm(mpg, lm_fit)
str(cv_error)
List of 4
 $ call : language cv.glm(data = mpg, glmfit = lm_fit)
 $ K    : num 392
 $ delta: num [1:2] 24.2 24.2
 $ seed : int [1:626] 10403 292 241038711 -724551010 -982165160 -1828904773 1219649113 639708648 -249104362 -1532738671 ...

?cv.glm

compute_poly_cv <- function(degree){
1  fit_degree <- glm(mpg ~ poly(horsepower, degree), data = mpg)
  cv_degree <- cv.glm(mpg, fit_degree)
  return(cv_degree)
}

2purrr::map(1:5, compute_poly_cv)  |>
  purrr::map("delta")
1
Aquí no estamos usando subset = train porque la función cv.glm() hace lo necesario
2
La función purrr::map(...) es un substituto sencillo del ciclo for en lenguaje funcional. También tenemos lapply 
[[1]]
[1] 24.23151 24.23114

[[2]]
[1] 19.24821 19.24787

[[3]]
[1] 19.33498 19.33448

[[4]]
[1] 19.42443 19.42371

[[5]]
[1] 19.03321 19.03242

compute_poly_cv_numeric <- function(degree){
  fit_degree <- glm(mpg ~ poly(horsepower, degree), data = mpg)
  cv_degree <- cv.glm(mpg, fit_degree)
  return(cv_degree$delta[1])
}

purrr::map_dbl(1:5, compute_poly_cv_numeric) 
[1] 24.23151 19.24821 19.33498 19.42443 19.03321

K-Fold Cross Validation

compute_poly_cv_numeric_k <- function(degree, k){
1  fit_degree <- glm(mpg ~ poly(horsepower, degree), data = mpg)
  cv_degree <- cv.glm(mpg, fit_degree, K = k)
  return(cv_degree$delta[1])
}

set.seed(17)
purrr::map_dbl(1:10, compute_poly_cv_numeric_k, k = 10)
1
Lo único que varía es que en la función cv.glm(...) usamos el argumento K para decir cuántos folds se debn usar. En este caso, usamos k = 10 en el llamado a compute_poly_cv_numeric_k 
 [1] 24.27207 19.26909 19.34805 19.29496 19.03198 18.89781 19.12061 19.14666
 [9] 18.87013 20.95520