1 Introducción

El problema es que el modelo se ajusta muy bien a los datos que hemos utilizado para estimarlo, luego es razonable pensar que va a predecir bien dichos datos. Pero en la práctica, queremos el modelo para predecir datos que no conocemos. El objetivo es analizar la predicción del modelo frente a datos no conocidos.

2 El método del subconjunto de validación (validation set approach)

  • Dividimos en dos partes los datos (50%-50%, 40%-60%,…).
  • En el training set estimamos el modelo.
  • En el test set (validation set) predecimos la respuesta y calculamos el Mean Squared Error (MSE) entre la variable respuesta observada y la predicha.

\[ MSE = \frac{1}{n}\sum _{i=1}^{n}{(y_i - \hat y_i)^2} \]

datos = read.csv("Advertising.csv")
str(datos)
## 'data.frame':    200 obs. of  5 variables:
##  $ X        : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ TV       : num  230.1 44.5 17.2 151.5 180.8 ...
##  $ radio    : num  37.8 39.3 45.9 41.3 10.8 48.9 32.8 19.6 2.1 2.6 ...
##  $ newspaper: num  69.2 45.1 69.3 58.5 58.4 75 23.5 11.6 1 21.2 ...
##  $ sales    : num  22.1 10.4 9.3 18.5 12.9 7.2 11.8 13.2 4.8 10.6 ...
  • Dividimos los datos en training set y test set:
n = nrow(datos)
n_train = 100
n_test = n - n_train
set.seed(115)

v = 1:n
pos_train = sample(v, n_train, replace = F) # muestreo sin reemplazamiento
pos_test = v[-pos_train]

# dividimos los datos en training set y validation set
datos_train = datos[pos_train,]
datos_test = datos[pos_test,]

plot(datos_train$TV, datos_train$sales, pch = 19, col = "blue")
points(datos_test$TV, datos_test$sales, pch = 19, col = "red")
legend(x=0,y=25, legend=c("train","test"), fill=c("blue","red"))

  • En el training set estimamos el modelo:
m1 = lm(sales ~ TV, data = datos_train)

# error en el training set
y_train = datos_train$sales
y_train_p = predict(m1,datos_train)
  • Error en el validation set calculamos el MSE:
source("MSE.R")
MSE # se muestra la definición de la funcion
## function (y, yp) 
## {
##     y = as.numeric(y)
##     yp = as.numeric(yp)
##     n = length(y)
##     d = (y - yp)^2
##     suma = sum(d)
##     MSE = 1/n * suma
##     return(MSE)
## }
y_test = datos_test$sales
y_test_p = predict(m1,datos_test)

(mse_test = MSE(y_test,y_test_p))
## [1] 11.6203

Problemas: - El MSE cambia en función de como se elija el training set y el validation set. - Sólo se incluye una parte de los datos para estimar el modelo, por lo que la estimación es peor que incluyéndolos a todos.

3 Leave-One-Out Cross-Validation

  • Se dividen los datos en dos partes
    • training set: \(\{(y_2,x_2), (y_3,x_3), \ldots, (y_n,x_n) \}\). Aquí se estima el modelo
    • test set: \((y_1,x_1)\). Calculamos el error de predicción

\[ MSE_1 = (y_1 - \hat y_1)^2 \] - Repetimos el proceso con: - training set: \(\{(y_1,x_1), (y_3,x_3), \ldots, (y_n,x_n) \}\). Aquí se estima el modelo - test set: \((y_2,x_2)\). Calculamos el error de predicción

\[ MSE_2 = (y_2 - \hat y_2)^2 \] - Repetimos el proceso para calcular \(MSE_3, MSE_4, \ldots, MSE_n\).

  • el error final es la media

\[ CV = \frac{1}{n}\sum _{i=1}^{n}{MSE_i} \]

error_i = rep(0,n)
for (i in 1:n){
  datos_train = datos[-i,]
  datos_test = datos[i,]
  mi = lm(sales ~ TV, data = datos_train)
  
  yi = datos_test$sales
  yi_p = predict(mi,datos_test)
  
  error_i[i] = MSE(yi,yi_p)
}
(error_i_media = mean(error_i))
## [1] 10.74109
plot(error_i)
lines(error_i_media*rep(1,n),col='red')

4 k-fold Cross-Validation

Se divide aleatoriamente los datos en k partes (o folds). Cada parte tiene n1 = n/k datos.

  • Para i = 1:k
    • la parte i constituye el test set.
    • las otras partes constituyen el train set.
    • se calcula el MSE(i)
  • El error total es:

\[ MSE_{TOTAL} = \frac{1}{k} \sum _{i=1}^k MSE_i \]

# numero de folds
num_folds = 5
n = nrow(datos)
n1 = trunc(n/num_folds)

# se definen las posiciones de test de cada fold
set.seed(342)
v = sample(1:n,n,replace = F) # vector auxiliar
pos = list()
for (k in 1:(num_folds-1)){
  pos[[k]] = v[((k-1)*n1+1):(k*n1)]
}
# el ultimo fold puede tener un numero diferente de datos
# (por redondeo)
pos[[num_folds]] = v[(n-n1+1):n]

v = 1:n
error_k = rep(0,num_folds)
for (k in 1:num_folds){
  # datos de training y de validation de cada fold
  pos_test = pos[[k]]
  pos_train = v[-pos_test]
  datos_train = datos[pos_train,]
  datos_test = datos[pos_test,]
  
  # estimamos el modelo
  mk = lm(sales ~ TV, data = datos_train)
  
  # error de prediccion
  y_va = datos_test$sales
  y_va_p = predict(mk,datos_test)
  
  error_k[k] = MSE(y_va,y_va_p)
}
(errork_media = mean(error_k))
## [1] 10.60946

4.1 Funciones para el cross-validation

Utilizando las definiciones anteriores, se van a definir funciones para poder aplicar cross-validation en temas posteriores:

  • funcion que calcula las posiciones de train y test dado el numero de folds (descargar):
source("cross_val_pos.R")
cross_val_pos # se muestra el contenido de la función
## function (num_datos, num_folds) 
## {
##     n1 = trunc(num_datos/num_folds)
##     v = sample(1:num_datos, num_datos, replace = F)
##     train = list()
##     test = list()
##     for (k in 1:(num_folds - 1)) {
##         pos_test = ((k - 1) * n1 + 1):(k * n1)
##         test[[k]] = v[pos_test]
##         train[[k]] = v[-pos_test]
##     }
##     pos_test = ((num_folds - 1) * n1 + 1):num_datos
##     test[[num_folds]] = v[pos_test]
##     train[[num_folds]] = v[-pos_test]
##     return(list(train = train, test = test))
## }

Aplicamos estas funciones a los datos para comprobar que funcionan correctamente:

set.seed(342)
num_datos = 10
num_folds = 3
cross_val_pos(num_datos,num_folds)
## $train
## $train[[1]]
## [1]  5  6  7  8 10  4  1
## 
## $train[[2]]
## [1]  3  9  2  8 10  4  1
## 
## $train[[3]]
## [1] 3 9 2 5 6 7
## 
## 
## $test
## $test[[1]]
## [1] 3 9 2
## 
## $test[[2]]
## [1] 5 6 7
## 
## $test[[3]]
## [1]  8 10  4  1