La validación cruzada como selector de modelos. Aplicación a los datos de prostate.data
Pablo Gregori
05/12/2016
Introducción
En la prevención del cáncer de próstata, existe un antígeno que se puede analizar y servir de alerta temprana. Se llama PSA (Prostate Specific Antigen), y se puede medir su concentración en la sangre de los pacientes. Se mide en logaritmos pues son valores muy pequeños, y el valor se almacena en la variable lpsa.
Existen otros indicadores en pacientes que pueden influir sobre el valor de lpsa. Estos indicadores los facilitan los expertos, gracias a su experiencia.
Para investigar la influencia de esos otros factores en el lpsa, se toma una muestra de individuos a los que se les ha medido todos esos rasgos, para que esa muestra nos enseñe la relación que hay.
Vamos a aplicar el modelo más sencillo que hay, el modelo lineal, para tratar de predecir variable lpsa a partir de las otras variables recogidas en los datos (descargar datos y ver explicación de los mismos).
El modelo lineal es demasiado sencillo para ser cierto en situaciones normales, pero si su aproximación es suficiente, podemos intuir la importancia de las distintas variables independientes en el resultado final.
La nube de puntos que dibuja la muestra de datos, en un espacio de más de 3 dimensiones, debe mostrar una tendencia general, pero los datos orbitan en torno a esa tendencia, porque dependen de otros factores que no han recogido los expertos, y que deben influir, aunque sea en menor medida, sobre el lpsa. Esa dispersión la veremos como ruido en torno a la verdadera tendencia de los datos.
El modelo lineal intenta encontrar esa tendencia, con el inconveniente de la rigidez que tienen las rectas, planos o hiperplanos.
Cuantas más variables incluya en el modelo lineal, más conseguiremos adaptar el hiperplano de regresión a la nube de puntos, lo cual parece bueno, pero corremos el peligro de adaptarnos al “ruido” de los datos, más que la tendencia.
La bondad del modelo no consiste en aproximar los datos de la muestra, sino aproximar datos futuros, de fuera de la muestra, de los que no conocemos el valor de lpsa.
El error fuera de la muestra es imposible de conocer, pero se puede estimar: la mejor manera de hacerlo es usar una parte de la muestra para “entrenar el modelo”, y otra parte para “comprobar el modelo” (medir el error).
Leemos los datos
x = read.table(file="prostate.data", header=TRUE)
str(x)## 'data.frame': 97 obs. of 10 variables:
## $ lcavol : num -0.58 -0.994 -0.511 -1.204 0.751 ...
## $ lweight: num 2.77 3.32 2.69 3.28 3.43 ...
## $ age : int 50 58 74 58 62 50 64 58 47 63 ...
## $ lbph : num -1.39 -1.39 -1.39 -1.39 -1.39 ...
## $ svi : int 0 0 0 0 0 0 0 0 0 0 ...
## $ lcp : num -1.39 -1.39 -1.39 -1.39 -1.39 ...
## $ gleason: int 6 6 7 6 6 6 6 6 6 6 ...
## $ pgg45 : int 0 0 20 0 0 0 0 0 0 0 ...
## $ lpsa : num -0.431 -0.163 -0.163 -0.163 0.372 ...
## $ train : logi TRUE TRUE TRUE TRUE TRUE TRUE ...summary(x)## lcavol lweight age lbph
## Min. :-1.3471 Min. :2.375 Min. :41.00 Min. :-1.3863
## 1st Qu.: 0.5128 1st Qu.:3.376 1st Qu.:60.00 1st Qu.:-1.3863
## Median : 1.4469 Median :3.623 Median :65.00 Median : 0.3001
## Mean : 1.3500 Mean :3.629 Mean :63.87 Mean : 0.1004
## 3rd Qu.: 2.1270 3rd Qu.:3.876 3rd Qu.:68.00 3rd Qu.: 1.5581
## Max. : 3.8210 Max. :4.780 Max. :79.00 Max. : 2.3263
## svi lcp gleason pgg45
## Min. :0.0000 Min. :-1.3863 Min. :6.000 Min. : 0.00
## 1st Qu.:0.0000 1st Qu.:-1.3863 1st Qu.:6.000 1st Qu.: 0.00
## Median :0.0000 Median :-0.7985 Median :7.000 Median : 15.00
## Mean :0.2165 Mean :-0.1794 Mean :6.753 Mean : 24.38
## 3rd Qu.:0.0000 3rd Qu.: 1.1787 3rd Qu.:7.000 3rd Qu.: 40.00
## Max. :1.0000 Max. : 2.9042 Max. :9.000 Max. :100.00
## lpsa train
## Min. :-0.4308 Mode :logical
## 1st Qu.: 1.7317 FALSE:30
## Median : 2.5915 TRUE :67
## Mean : 2.4784 NA's :0
## 3rd Qu.: 3.0564
## Max. : 5.5829N = dim(x)[1]Tenemos \(97\) individuos en la muestra.
Observar que la columna train es una columna artificial, creada en el libro de Hastie para separar los datos de la muestra en las dos partes comentadas.
Ajuste del modelo lineal con una variable independiente (por ejemplo, lweight)
r1 = lm(formula=lpsa ~ lweight, data=x) # la función que hace el trabajo
summary(r1) # información de la regresión lineal realizada##
## Call:
## lm(formula = lpsa ~ lweight, data = x)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.27976 -0.67507 -0.03503 0.53984 2.93649
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.7586 0.9103 -1.932 0.0564 .
## lweight 1.1676 0.2491 4.686 9.28e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.046 on 95 degrees of freedom
## Multiple R-squared: 0.1878, Adjusted R-squared: 0.1792
## F-statistic: 21.96 on 1 and 95 DF, p-value: 9.276e-06Los Residuals son los errores cometidos por la recta de regresión al no pasar por los puntos de la muestra.
La ecuación de la recta de regresión es
\[\text{lpsa} = -1.7585989 + 1.1675538 \times \text{lweight}.\]
En este sencillo caso, podemos representar gráficamente los datos junto a su recta de regresión (y de paso señalar los errores).
plot(x[,c('lweight', 'lpsa')]) # dibuja los puntos
abline(r1) # dibuja la recta
for(i in 1:N) {
points(x=rep(x[i, 'lweight'],2),
y=c(r1$fitted.values[i], x[i, 'lpsa']),
type='l', col="red") # para pintar los errores
}
El error cuadrático total lo podesmo tomar como 1.0822824, y podemos usarlo como indicador de lo que falla el modelo lineal ajustado, al intentar aproximar los datos.
La predicción de lpsa para un nuevo individuo cuya lweight es conocida, se haría a partir de la ecuación de la recta encontrada.
Modelos lineales anidados: añadir cada vez una variable independiente (la de mejor ajuste)
Con 1 sola variable:
x2 = x[,1:8] # solo variables independientes
ErrorCuadratico = numeric(0)
ErrorCuadraticoAcum = ErrorCuadratico
for(i in names(x2)) { # recorremos var. indep.
reg = lm(formula=paste("lpsa ~ ", i), data=x) # ajustamos
ErrorCuadratico[i] = var(reg$residuals) # calculamos error
}
VarElegida = names(which.min(ErrorCuadratico)) # var. elegida
VarQuedan = names(x2)[ names(x2) != VarElegida ] # la quitamos
Formulas = character(0) # aquí guardaremos las fórmulas
ErrorCuadraticoAcum = c(ErrorCuadraticoAcum, min(ErrorCuadratico))La variable que ofrece una recta de regresión con menor error cuadrático es… lcavol, y el error cuadrático es 0.6136957.
Ahora añadimos a ésta, una segunda variable, para obtener el mejor modelo con 2 variables independientes, y así sucesivamente hasta incorporar todas:
VarQuedan = VarQuedan[ VarQuedan != VarElegida ]
FormulaAnt = paste("lpsa ~ ", VarElegida)
Formulas = c(Formulas, FormulaAnt)
while( length(VarQuedan)>0 ) {
ErrorCuadratico = numeric(0)
for( i in VarQuedan ) {
reg = lm(formula=paste(FormulaAnt, i, sep="+"), data=x)
# EXTRAER ERROR CUADRATICO DEL MODELO
ErrorCuadratico[i] = var(reg$residuals)
}
VarElegida = names(which.min(ErrorCuadratico))
VarQuedan = VarQuedan[ VarQuedan != VarElegida ]
FormulaAnt = paste(FormulaAnt, VarElegida, sep="+")
Formulas = c(Formulas, FormulaAnt)
ErrorCuadraticoAcum = c(ErrorCuadraticoAcum, min(ErrorCuadratico))
}
Formulas## [1] "lpsa ~ lcavol"
## [2] "lpsa ~ lcavol+lweight"
## [3] "lpsa ~ lcavol+lweight+svi"
## [4] "lpsa ~ lcavol+lweight+svi+lbph"
## [5] "lpsa ~ lcavol+lweight+svi+lbph+age"
## [6] "lpsa ~ lcavol+lweight+svi+lbph+age+pgg45"
## [7] "lpsa ~ lcavol+lweight+svi+lbph+age+pgg45+lcp"
## [8] "lpsa ~ lcavol+lweight+svi+lbph+age+pgg45+lcp+gleason"Los errores cuadráticos que se cometen en cada modelo son: 0.6136957, 0.538981, 0.4850879, 0.4749528, 0.4628821, 0.4559997, 0.4490371, 0.4485252. Como se ve, cada vez se disminuye.
Validación cruzada: estimar el error fuera de la muestra
El objetivo es predecir el valor de lpsa para individuos de fuera de la muestra, con el menor error posible, lo cual es inaccesible ahora, ya que no tenemos individuos fuera de la muestra, de los que conozcamos el valor de lpsa.
Para estimar ese error fuera de la muestra, se debe tomar la muestra actual y:
- Elegir una parte como conjunto de entrenamiento
- Elegir otra parte como conjunto de validación
- Ajustar el modelo al conjunto de entrenamiento
- Valorar el error en el conjunto de validación
El mejor modelo será aquel que dé menor error en la parte de validación.
Hay muchas formas de realizar esa división de la muestra:
- Dejar uno fuera (leave-one-out): se entrena el modelo con todos los datos menos uno, y se valida con el dato extraído.
- Validación en k pliegues (k-fold validation): se divide la muestra en \(k\) partes de igual tamaño, y se prodece a entrenar el modelo con \(k-1\) partes y ajustar con la restante. Elecciones habituales de \(k\) son \(5\), \(10\) y \(1/N\) (= leave-one-out).
En todos los casos, se repite el proceso variando el conjunto de validación, y se promedian los errores cuadráticos calculados para eliminar el posible sesgo de la elección.
En este caso vamos a tomar \(k = 5\) grupos, por lo que los tamaños de los grupos serán 19, 19, 19, 19, 21, y la asignación será completamente aleatoria.
Para poder “trazar” en caso de errores, usamos una semilla (que suprimiremos en la última versión del documento).
# asignamos los individuos a los grupos aleatoriamente
set.seed(20161205) # fecha de hoy
grupo = sample(x=rep(1:k, times=Nj), size=N, rep=FALSE)
# el vector grupo servirá para seleccionar con corchetes
# los individuos que entran en el ajuste, y los que
# entran en la validaciónAhora comenzamos con cada uno de los modelos lineales ajustados en la sección anterior.
- For (cada modelo lineal)
- For (cada conjunto de validación)
- Ajustar el modelo al conjunto de entrenamiento
- Calcular el error cuadrático sobre el conjunto de validación
- Promediar y guardar el error cuadrático
- For (cada conjunto de validación)
- Comparar los errores cuadráticos de todos los modelos
- Elegir el modelo con menor error cuadrático