En este laboratorio (Basado en el laboratorio 5.3 del libro ISLR) exploraremos tecnicas de remuestreo.
En el capitulo 5 del libro ISLR se exploran los usos del enfoque de las validaciones para poder estimar la tasa de error, de como entrenar diferentes modelos lineares utilizando data sets Auto.
Para esta laboratorio empesaremos utilizando la funcion sample para generar un dataset con los indices que vamos a utilizar para entrenar el modelo.
library(ISLR)
set.seed(1)
train <- sample(392, 196)
train [1] 105 146 224 354 79 348 365 255 242 24 388 68 262 391 292 188 270 372 143 290 387 382 384 47 99 142
[27] 5 140 317 124 175 217 178 67 297 239 283 39 257 379 289 228 275 194 185 274 9 165 252 238 164 294
[53] 149 83 383 34 107 174 222 136 304 98 152 110 214 85 157 250 28 356 329 376 111 336 330 323 347 123
[79] 245 301 333 334 363 101 234 63 218 38 75 44 73 18 193 263 233 237 135 121 357 360 192 103 371 287
[105] 183 62 305 137 299 170 276 206 100 295 42 4 198 29 315 362 321 296 131 369 203 122 312 55 61 326
[131] 151 21 10 167 240 154 144 271 251 129 173 380 60 65 181 112 303 288 26 211 340 385 373 109 120 43
[157] 125 313 249 50 359 207 291 179 201 94 15 76 163 225 386 186 189 86 339 195 311 160 130 300 307 41
[183] 187 106 314 40 284 370 213 247 256 258 261 375 57 117Ahora calculamos los test index
test
[1] 1 2 3 6 7 8 11 12 13 14 16 17 19 20 22 23 25 27 30 31 32 33 35 36 37 45
[27] 46 48 49 51 52 53 54 56 58 59 64 66 69 70 71 72 74 77 78 80 81 82 84 87 88 89
[53] 90 91 92 93 95 96 97 102 104 108 113 114 115 116 118 119 126 127 128 132 133 134 138 139 141 145
[79] 147 148 150 153 155 156 158 159 161 162 166 168 169 171 172 176 177 180 182 184 190 191 196 197 199 200
[105] 202 204 205 208 209 210 212 215 216 219 220 221 223 226 227 229 230 231 232 235 236 241 243 244 246 248
[131] 253 254 259 260 264 265 266 267 268 269 272 273 277 278 279 280 281 282 285 286 293 298 302 306 308 309
[157] 310 316 318 319 320 322 324 325 327 328 331 332 335 337 338 341 342 343 344 345 346 349 350 351 352 353
[183] 355 358 361 364 366 367 368 374 377 378 381 389 390 392Ahora basado en esto obtenemos un dataset con 196 valores entre 1 y 392, ahora con este dataset vamos a generar una regresion lineal y utilizando el dataset Auto:
str(Auto)
'data.frame': 392 obs. of 9 variables:
$ mpg : num 18 15 18 16 17 15 14 14 14 15 ...
$ cylinders : num 8 8 8 8 8 8 8 8 8 8 ...
$ displacement: num 307 350 318 304 302 429 454 440 455 390 ...
$ horsepower : num 130 165 150 150 140 198 220 215 225 190 ...
$ weight : num 3504 3693 3436 3433 3449 ...
$ acceleration: num 12 11.5 11 12 10.5 10 9 8.5 10 8.5 ...
$ year : num 70 70 70 70 70 70 70 70 70 70 ...
$ origin : num 1 1 1 1 1 1 1 1 1 1 ...
$ name : Factor w/ 304 levels "amc ambassador brougham",..: 49 36 231 14 161 141 54 223 241 2 ...lm_fit <- lm(mpg~horsepower, data=Auto, subset=train)
lm_fit
Call:
lm(formula = mpg ~ horsepower, data = Auto, subset = train)
Coefficients:
(Intercept) horsepower
40.3404 -0.1617 Vemos que para esta regresion obtenemos que \(\beta_{0}=40.3404\) y \(\beta_{1}=-016.17\), ahora procedemo a calcular el MSE utilizando el test:
mean(((Auto$mpg - predict(lm_fit, Auto))[test])^2)
[1] 26.14142Con esto vemos que el fit para la regresion lineal es de 26.14142, ahora utilizando la funcion poly vamos a probrar el fit para una funcion cuadratica:
lm_fit_2 <- lm(Auto$mpg~poly(Auto$horsepower, 2), data=Auto, subset=train)
mean(((Auto$mpg - predict(lm_fit_2, Auto))[test])^2)[1] 19.82259y la cubica:
lm_fit_3 <- lm(Auto$mpg~poly(Auto$horsepower, 3), data=Auto, subset=train)
mean(((Auto$mpg - predict(lm_fit_3, Auto))[test])^2)[1] 19.78252Como vemos los resultados finales son \(mse=26.14142\) \(mse^2=19.82259\) \(mse^3=19.78252\), aparentemente la regresion cubica es la que mejor fit tiene, pero ahora intenemos recalcular el dataset de entrenamiento:
set.seed(2)
train <- sample(392,196)
test <- (1:392)[-train]y calculamos las 3 regresiones y su mse
lm_fit <- lm(mpg~horsepower, data=Auto, subset=train)
mean(((Auto$mpg-predict(lm_fit, Auto))[test])^2)[1] 23.29559lm_fit_2 <- lm(mpg~poly(horsepower,2), data=Auto, subset=train)
mean(((Auto$mpg-predict(lm_fit_2, Auto))[test])^2)[1] 18.90124lm_fit_3 <- lm(mpg~poly(horsepower,3), data=Auto, subset=train)
mean(((Auto$mpg-predict(lm_fit_3, Auto))[test])^2)[1] 19.2574Ahora al recalcular la muestra los resultados son \(mse=23.29559\) \(mse^2=18.90124\) \(mse^3=19.2574\), con lo que la regresion cuadratica nos da un mejor fit que las anteriores.
Para utilizar el LOOCV utlizaremos la funcion glm, cuando utilizamos esta funcion sin utilizar el parametro familiy="binomial", la funcion glm crea una regresion lineal como la funcion lm, podemos hacer la prueba:
glm_fit <- glm(mpg~horsepower, data=Auto)
coef(glm_fit)(Intercept) horsepower
39.9358610 -0.1578447 ahora utilizando lm
lm_fit <- lm(mpg~horsepower, data=Auto)
coef(lm_fit)(Intercept) horsepower
39.9358610 -0.1578447 Vemos que para ambas funciones \(\beta_{0}=39.9358610\) y \(\beta_{1}=-0.1578447\). Por esta razon utilizaremos glm y utilizaremos cv.glm para realizar la validacion cruzada:
library(boot)
glm_fit <- glm(mpg~horsepower, data=Auto)
cv_error <- cv.glm(Auto, glm_fit)
cv_error$delta[1] 24.23151 24.23114Como vemos cv_error nos da dos valores, los cuales son identicos, ahora procederemos a calcular estos indicadores delta para los polinimos desde 1 hasta 5.
glm_error_function <- function(grado, dataset){
glm_fit <- glm(mpg~poly(horsepower, grado), data=dataset)
return(cv.glm(dataset, glm_fit)$delta[1])
}
cv_error <- lapply(1:5, glm_error_function, dataset=Auto)
cv_error[[1]]
[1] 24.23151
[[2]]
[1] 19.24821
[[3]]
[1] 19.33498
[[4]]
[1] 19.42443
[[5]]
[1] 19.03321Esta funcion nos permite ver que hay una gran mejora entre el polinomio lineal y el cuadratico, pero no hay una mejor considerable pasando de ahi.
Utilizando la funcion cv.glm podemos tambien implementar la funcion Kfold, para ello haremos la misma implementacion de la funcion anterior, utilizando \(k=10\) el cual es un valor comun a utilizar en este tipo de ejercicios:
set.seed(17)
glm_error_k <- function(grado, dataset){
glm_fit <- glm(mpg~poly(horsepower, grado), data=dataset)
return (cv.glm(dataset, glm_fit, K=10)$delta[1])
}
cv_error_10 <- lapply(1:5, glm_error_k, dataset=Auto)
cv_error_10[[1]]
[1] 24.2052
[[2]]
[1] 19.18924
[[3]]
[1] 19.30662
[[4]]
[1] 19.33799
[[5]]
[1] 18.87911Al utilizar esta nueva funcion vemos que los resultados son muy similares con la ventaja que tarda mucho menos en computar la validacion. Pero seguimos viendo que el mejor fit es el cuadratico.
Una de las ventajas de hacer Bootstrap es que no requiere complicadas funciones matematicas, para ello haremos uso de la funcion boot() presente en la libreria boot. Comenzaremos definiendo la funcion alpha_fn:
Esta funcion calcula las \(\alpha\) de los indices que se le provea segun un dataset:
alpha_fn(Portfolio, 1:100)[1] 0.5758321Ahora lo haremos utilizando la funcion sample para utilizar los indices sin orden:
set.seed(1)
alpha_fn(Portfolio, sample(100, 100, replace = T))[1] 0.5963833Ahora podemos calcular el Bootstrap utilizando la funcion boot y la funcion alpha_fn y procedemos a ejecutar esta prueba 1000 veces:
boot(Portfolio, alpha_fn, R=1000)
ORDINARY NONPARAMETRIC BOOTSTRAP
Call:
boot(data = Portfolio, statistic = alpha_fn, R = 1000)
Bootstrap Statistics :
original bias std. error
t1* 0.5758321 -7.315422e-05 0.08861826Con esto vemos que \(\hat{\alpha}=0.5758\) y que el \(SE(\hat{\alpha})=0.08861\).
La tecnica de bootstrap puede ser muy efectiva para medir la variabilidad de los coeficientes estimados y de prediccion desde un punto de vista de aprendisaje estadistico. Para ello estimaremos la variabilidad de \(\beta_{0}\) y \(\beta_{1}\) de la funcion \(Y=\beta_{0} + \beta_{1}X\) donde \(Y=mpg\) y \(X=horsepower\) del dataset Auto y compararemos los estimados obtenidos usando el bootstrap de las formulas de \(SE(\hat{\beta}_{0})\) y \(SE(\hat{\beta}_{1})\).
boot_fn <- function(data, index) {
return (coef(lm(mpg~horsepower, data=data, subset = index)))
}
boot_fn(Auto, 1:392)(Intercept) horsepower
39.9358610 -0.1578447 Ahora usaremos la funcion con los indices desordenados:
set.seed(1)
boot_fn(Auto, sample(392, 392, replace=T))(Intercept) horsepower
38.7387134 -0.1481952 boot_fn(Auto, sample(392, 392, replace=T))
(Intercept) horsepower
40.0383086 -0.1596104 Ahora usaremos la funcion boot para computar los errores standard de 1000 estimaciones de bootstrap:
boot(Auto, boot_fn, 1000)
ORDINARY NONPARAMETRIC BOOTSTRAP
Call:
boot(data = Auto, statistic = boot_fn, R = 1000)
Bootstrap Statistics :
original bias std. error
t1* 39.9358610 0.02972191 0.860007896
t2* -0.1578447 -0.00030823 0.007404467Como vemos \(SE(\hat{\beta_{0}})=0.8600\) y \(SE(\hat{\beta_{1}})=0.0074\), ahora revisemos el resumen:
summary(lm(mpg~horsepower, data=Auto))$coef Estimate Std. Error t value Pr(>|t|)
(Intercept) 39.9358610 0.717498656 55.65984 1.220362e-187
horsepower -0.1578447 0.006445501 -24.48914 7.031989e-81El resumen nos indica que \(SE(\hat{\beta_{0}})=0.7175\) y \(SE(\hat{\beta_{1}})=0.0064\) distinto de los valores que obtuvimos al realizar el bootstrapping, esto es un problema con el algoritmod e bootstrap? No realmente ya que la ecuacion de \(\alpha^2\) que se plantean en el libro, lo que depende de varias asumciones y ruido que conlleva dicha varaible. Ahora realicemos la misma prueba, utilizando la ecuacion cuadratica para lo que calcularemos \(SE(\hat{\beta_{0}})\), \(SE(\hat{\beta_{1}})\) y \(SE(\hat{\beta_{2}})\)
boot_fn <- function(data, index){
return(coefficients(lm(mpg~horsepower+I(horsepower^2), data=data, subset=index)))
}
set.seed(1)
boot(Auto, boot_fn, 1000)
ORDINARY NONPARAMETRIC BOOTSTRAP
Call:
boot(data = Auto, statistic = boot_fn, R = 1000)
Bootstrap Statistics :
original bias std. error
t1* 56.900099702 6.098115e-03 2.0944855842
t2* -0.466189630 -1.777108e-04 0.0334123802
t3* 0.001230536 1.324315e-06 0.0001208339summary(lm(mpg~horsepower+I(horsepower^2), data=Auto))$coef Estimate Std. Error t value Pr(>|t|)
(Intercept) 56.900099702 1.8004268063 31.60367 1.740911e-109
horsepower -0.466189630 0.0311246171 -14.97816 2.289429e-40
I(horsepower^2) 0.001230536 0.0001220759 10.08009 2.196340e-21Como ya lo hemos visto, al definir un modelo cuadratico los coeficientes si dan el mismo resultando, con lo cual terminamos de definir que el modelo cuadratico es el que mejor se ajusta.