Las técnicas de re-sampling son bastante utiles y poderosas cuando la cantidad de información no es limitada o si se quiren hacer modelos con resultados más robustos.

library("dplyr")
library("boot")
library("purrr")
library("ggplot2")
library("MASS")
auto <- as_tibble(ISLR::Auto)

Existen desventaja con el enfoque de train-test para modelos predictivos, y es que cada una de estas particiones pueder ser muy variable arrojando resultados distintos dependiente el muestreo que se realice.

ggplot(auto, aes(horsepower, mpg)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE) +
  geom_smooth(method = "lm", formula = y ~ poly(x, 2), se = FALSE, linetype = 2) +
  geom_smooth(method = "lm", formula = y ~ poly(x, 3), se = FALSE, linetype = 3) +
  geom_smooth(method = "lm", formula = y ~ poly(x, 4), se = FALSE, linetype = 4) +
  theme_minimal()
## `geom_smooth()` using formula 'y ~ x'

Hagamos una revisión de la metodología tradicional train-test sobre el dataset de autos:

# Plantamos semilla 
set.seed(1)

# Creamos los datasets de entrenamiento y testeo
sample <- sample(c(TRUE, FALSE), nrow(auto), replace = T, prob = c(0.6,0.4))
train <- auto[sample, ]
test <- auto[!sample, ]


# Creamos un DF para evaluar el error con polinomios de distinto orden
mse.df <- tibble(degree = 1:10, mse = NA)

# Estimamos una regresión, una por cada orden del polinomio creado
for(i in 1:10) {
  lm.fit <- lm(mpg ~ poly(horsepower, i), data = train)
  # MSE
  mse.df[i, 2] <- mean((test$mpg - predict(lm.fit, test))^2)
}

ggplot(mse.df, aes(degree, mse)) +
  geom_line() +
  geom_point() +
  ylim(c(10, 30))+
  theme_classic()

Existe una disminución considerable del error cuando se aumenta en 2 el grado del polinomio en la regresión. Sin embargo, este resultado sigue dependiendo de cómo se hagan los splits de la base de datos, por lo que si intentamos recrear nuevos train-test sets tendremos un error distinto:

mse.df.2 <- tibble(sample = vector("integer", 100), 
                   degree = vector("integer", 100), 
                   mse = vector("double", 100))
counter <- 1

for(i in 1:10) {
  # Cambios la semilla para que cada train-test split sea distinto
  set.seed(i)
  
  # Hacemos el split 
  sample <- sample(c(TRUE, FALSE), nrow(auto), replace = T, prob = c(0.6,0.4))
  train <- auto[sample, ]
  test <- auto[!sample, ]
  
  # Por cada split debemos correr 10 modelos distintos
  for(j in 1:10) {
    lm.fit <- lm(mpg ~ poly(horsepower, j), data = train)
    
   # Agregamos al dataframe los datos correspondientes al error y el grado polinomial
    mse.df.2[counter, 2] <- j
    mse.df.2[counter, 3] <- mean((test$mpg - predict(lm.fit, test))^2)
    
    # Identificado de muestra
    mse.df.2[counter, 1] <- i
    counter <- counter + 1
  }
  next
}

ggplot(mse.df.2, aes(degree, mse, color = factor(sample))) +
  geom_line(show.legend = FALSE) +
  geom_point(show.legend = FALSE) +
  ylim(c(10, 30)) + theme_classic()

Efectivamente, el error varia de muestra a muestra con un rango de aproximadamente 10, que es bastante amplio.

Cross validation

Leave-one-out Cross-validation (LOOCV)

El método consiste en partir los datos en dos grupos, entrenamiento y validación. La principal diferencia con el enfoque train-test es que el tamaño del set de validación es de 1, es decir, que entrenamos el modelo con n-1 observaciones. Por lo que el MSE de validación será:

\[ MSE_{1} = (y_1-\hat{y_1})^2 \]

Sin embargo la estimación es altamente variable, dado que estamos basando todo el error del modelo en una única observación.Aquí es donde entra la validación cruzada. Este proceso no se hace una unica vez, por el contrario la hacemos \(n\) veces, cada vez cambiando los datos de entrenamiento y los datos de evaluación. En total tendríamos n valores del MSE, por lo que sacamos su promedio para evaluar su desempeño total:

\[ CV_{(n)} = \frac{1}{k}\sum_{i=1}^2{MSE_i} \] Podemos realizar esto de la siguiente manera:

# Regresión linea habitual
glm_fit <- glm(mpg ~ horsepower, data = auto)
coef(glm_fit)
## (Intercept)  horsepower 
##  39.9358610  -0.1578447

Usamos el la función cv.glm para realizar LOOCV (debemos tener previamente estimado el modelo al que deseamos hacer cross-validation)

loocv_err <- cv.glm(auto, glm_fit)

str(loocv_err)
## List of 4
##  $ call : language cv.glm(data = auto, glmfit = glm_fit)
##  $ K    : num 392
##  $ delta: num [1:2] 24.2 24.2
##  $ seed : int [1:626] 10403 392 -1703707781 1994959178 434562476 -1277611857 -1105401243 1020654108 526650482 -1538305299 ...
# Estamos interesados en el error medio del procediento

paste0("El MSE de LOOCV es ", loocv_err$delta[1])
## [1] "El MSE de LOOCV es 24.2315135179293"

Ahora podemos revisar también por medio de validación cruzada el modelo con polinomios de distinto orden

# create function that computes LOOCV MSE based on specified polynomial degree
loocv_error <- function(x) {
  glm.fit <- glm(mpg ~ poly(horsepower, x), data = auto)
  cv.glm(auto, glm.fit)$delta[1]
}

# compute LOOCV MSE for polynomial degrees 1-5
1:5 %>% map_dbl(loocv_error)
## [1] 24.23151 19.24821 19.33498 19.42443 19.03321

Hay que notar es questa metodología podría no ser computacionalmente eficiente si hay que estimar n modelos cuando n es demasiado grande. Por ello podría ser más viable trabajar con menos iteraciones y con grupos más grandes. A esta metodología se le llama K-fold cross validation.

K-fold cross validation

LOOCV en realidad es un caso especial de k-fold (o k-pliegues) cross validation, en donde k (el número de pliegues deseados para partir los grupos) es igual a n. En este caso entonces se estimaran k modelos y se partiran los datos en k muestras de igual tamaño donde cada una hacer parte del conjunto de entrenamiento

\[ CV_{(k)} = \frac{1}{n}\sum_{i=1}^k{MSE_i} \]

Los valores de k más frecuentes son 5 y 10 pliegues pero su elección puede depender de más factores. La implementación en R sería:

# Crea,ps una función que estime un modelo de regresión, 10-fold cross-validation y
kfcv_error <- function(x) {
  glm.fit <- glm(mpg ~ poly(horsepower, x), data = auto)
  cv.glm(auto, glm.fit, K = 10)$delta[1]
}

# Calculamos los errores para la regresión con polinomios entre 1 y 5
1:5 %>% map_dbl(kfcv_error)
## [1] 24.24098 19.16217 19.17898 19.43967 19.11506
# Comparamos los resultados con la metodología LOOCV
1:5 %>% map_dbl(loocv_error)
## [1] 24.23151 19.24821 19.33498 19.42443 19.03321

El mismo enfoque puede ser aplicado a problemas de clasificación

\[ CV_{(n)} = \frac{1}{n}\sum_{i=1}^k{Err_i} \] donde \(Err_i = I(y_i - \hat{y_i})\), es decir, la veces que el modelo estimado se equivocó con la predicción

stock <- ISLR::Smarket

# Estimamos una regresión logística
glm.fit <- glm(Direction ~ Lag1 + Lag2, family = binomial, data = stock)

# Creamos una función costo que se ajuste a un modelo de classificatión 
cost <- function(r, pi = 0) mean(abs(r - pi) > 0.5)

# Estimamos el error con cross validation
cv.glm(stock, glm.fit, cost, K = 10)$delta[1]
## [1] 0.4856

Ejemplos con análisis linear discriminante y análisis cuadratico discrimiante, basados en el mismo ejemplo que la regresión logística

library(MASS)

# fit discriminant analysis models with CV = TRUE for LOOCV
lda.fit <- lda(Direction ~ Lag1 + Lag2, CV = TRUE, data = stock)
qda.fit <- qda(Direction ~ Lag1 + Lag2, CV = TRUE, data = stock)

# compute estimated test error based on cross validation
mean(lda.fit$class != stock$Direction)
## [1] 0.4816
mean(qda.fit$class != stock$Direction)
## [1] 0.4872

Bootstrapping

Herramienta estadística que puede ser usada para cuantificar la incertidumbre asociada con un estimador dado o un metodo estadistico.

Bootstrapping repetidamente hace “sorteos” para sacar muestras independientes de los datos originales para crear datos “bootstrap”. Este muestreo es realizado con reemplazo (cada observación puede aparecer más de una vez en el muestreo).

Pasos para realizar boostrap

  1. Crear una función que calcule el estadístico de interés (media, varianza, o lo que sea)
  2. Usar la función boot() para realizar boostrapping
  3. Revisamos el error estandar

Hagamos un ejemplo rápido de cómo funciona la técnica:

x <- c(60, 75, 80, 85, 90)
mean(x)
## [1] 78
sd(x)
## [1] 11.51086
(str_err <- sd(x)/sqrt(length(x)))
## [1] 5.147815
B <- 100 # Veces que deseo realizar el remuestreo

st_resampled <- vector(mode = "numeric", 100) # Distribución boostrap
for (i in 1:B) {
  resampled_x <- sample(x,5, replace=TRUE)
  st_resampled[i] <- mean(resampled_x)
  
}

sd(st_resampled) #Error estandar boostrap de la media
## [1] 4.233536
hist(st_resampled,breaks = 10)

Incrementos de B, nos darán una mejor estimación de la distribución muestral y más confiabilidad en el error estandar

Ahora un ejemplo más complejo: La idea del ejemplo es minimizar el riesgo de inversión de un portafolio con montos de los activos X,Y

# Creamos el estadístico que estmoa interesados en estimar
statistic <- function(data, index) {
  x <- data$X[index]
  y <- data$Y[index]
  (var(y) - cov(x, y)) / (var(x) + var(y) - 2* cov(x, y))
}

# cargamos la data
portfolio <- ISLR::Portfolio

# El estadístico inicialmente nos da este valor
statistic(portfolio, 1:100)
## [1] 0.5758321
# Probemos si hacemos un muestreo con reemplazo del estadístico sobre los mismos datos
# En este caso es bootstrapping con B=1
statistic(portfolio, sample(100, 100, replace = TRUE))
## [1] 0.4385693
# La función boot nos permite recrear el paso de arriba la veces que consideremos necesarias

#También nos da el error estandar
set.seed(123)
result <- boot(portfolio, statistic, R = 1000)
result
## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = portfolio, statistic = statistic, R = 1000)
## 
## 
## Bootstrap Statistics :
##      original      bias    std. error
## t1* 0.5758321 0.007544609  0.08952496
# Podemos ver también el intervalo de confianza
boot.ci(result, type = "basic")
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = result, type = "basic")
## 
## Intervals : 
## Level      Basic         
## 95%   ( 0.3918,  0.7449 )  
## Calculations and Intervals on Original Scale
plot(result)

Este documento es una copia del artículo publicado por la universidad de Cincinnati http://uc-r.github.io/resampling_methods