Laboratorio 5

5.3 Lab: Cross-Validation and the Bootstrap (RESAMPLING)

5.3.1 The Validation Set Aproach

Utilizamos el dataset de Auto. Tambien la libreria ISLR.

setwd("~/CICLO 6/FIABILIDAD")
getwd()

[1] "C:/Users/Ericka/Documents/CICLO 6/FIABILIDAD"

library(dplyr)
library(ggplot2)
library(ISLR)

Auto <- read.csv(file = "auto.csv")
names(Auto)

[1] "mpg"          "cylinders"    "displacement" "horsepower"   "weight"       "acceleration" "year"        
[8] "origin"       "name"        

dim(Auto)

[1] 397   9

Auto$horsepower <- as.numeric(Auto$horsepower)
print(Auto)

Realizamos un split seleccionando 196 observaciones en forma aleatoria de las 397.

set.seed (1) ## Generando un numero aleatorio
train=sample(397, 196)
length(train)

[1] 196

Realizamos una regresion linear lm utilizando la informacion de train

lm.fit =lm(mpg~horsepower, data=Auto, subset= train )
summary(lm.fit)

Call:
lm(formula = mpg ~ horsepower, data = Auto, subset = train)

Residuals:
     Min       1Q   Median       3Q      Max 
-13.4692  -5.3321  -0.4667   4.8524  21.7856 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 17.89059    1.00936  17.725  < 2e-16 ***
horsepower   0.11167    0.01709   6.535 5.48e-10 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 7.172 on 194 degrees of freedom
Multiple R-squared:  0.1804,    Adjusted R-squared:  0.1762 
F-statistic:  42.7 on 1 and 194 DF,  p-value: 5.479e-10

Utilizamos la funcion predict para estimar la respuesta de los 397 observaciones, y utilizamos mean para calcular la media o MSE de los 196

#attach(Auto)
mean((mpg-predict(lm.fit, Auto))[-train ]^2) # -train seleciona los que nos estan en el training set

[1] 49.47809

Por lo tanto el estimado de la media en test es MSE = 49.47809

Utilizando la funcion poly para estimar el error en el test para las regresiones polinomiales y cubicas

lm.fit2=lm(mpg~poly(horsepower,2) ,data=Auto ,subset =train ) ## Regresion polinomial
mean((mpg-predict(lm.fit2,Auto))[-train ]^2)

[1] 49.28182

lm.fit3=lm(mpg~poly(horsepower,3) ,data=Auto ,subset =train ) ## Regresion cubica
mean((mpg-predict(lm.fit3 ,Auto))[-train ]^2)

[1] 31.69728

Los errores son de 49.28 y 31.69, asi que es mejor utilizar otro train

Calculamos de nuevo las regresiones polinomiales y cubicas para ver si mejora el error

set.seed (2)
train=sample(397,196)
lm.fit=lm(mpg~horsepower,subset=train) ## Lineal
mean((mpg-predict (lm.fit,Auto))[-train]^2)

[1] 44.79806

lm.fit2=lm(mpg~poly(horsepower,2),data=Auto ,subset =train ) ## Polinomial
mean((mpg-predict (lm.fit2 ,Auto))[-train]^2)

[1] 44.51045

lm.fit3=lm(mpg~poly(horsepower,3),data=Auto ,subset =train ) ## Cubica
mean((mpg-predict (lm.fit3 ,Auto))[-train]^2)

[1] 30.70154

Y si vemos que mejora un poco pero no significativamente, y si son consistentes a los calculos anteriores

5.3.2 Leave-One-Out Cross-Validation

Abreviado como estimacion de LOOCV, se puede calcular utilizando la funcion glm, este puede realizar la regresion lineal sin incluir la familia family=“binomial”, y al relizarlo nos da como resultado lo que nos da el lm

glm.fit=glm(mpg~horsepower, data=Auto)
coef(glm.fit)

(Intercept)  horsepower 
 17.8076120   0.1108047 

lm.fit=lm(mpg~horsepower,data=Auto)
coef(lm.fit)

(Intercept)  horsepower 
 17.8076120   0.1108047 

Podemos ver que son los mismos resultados

En este laboratorio utilizamos glm, y con ello la funcion cv.glm(), con la libreria boot, cv es cross-validation

library(boot)
glm.fit=glm(mpg~horsepower,data=Auto)
cv.err=cv.glm(Auto,glm.fit)
cv.err$delta

[1] 50.59100 50.59047

La funcion cv.glm() produce varios componentes, y los dos calores en el delta son los resultados de cross-validation

En este caso son iguales y corresponden a LOOCV, nuestro estimado de cv es 50.59100

Podemos realizar los calculos de nuevo corriendo un for de 1 a 5 con el fin de buscar una mejora en la regresion polinomial

cv.error=rep (0,5)
for (i in 1:5){
  glm.fit=glm(mpg~poly(horsepower,i),data=Auto)
  cv.error[i]=cv.glm (Auto,glm.fit)$delta[1] }
cv.error

[1] 50.59100 50.65400 32.51864 32.72261 27.32256

5.3.3 k-Fold Cross-Validation

La funcion cv.glm tambien nos ayuda a implementar k-Fold Cross-Validation

En este ejemplo realizamos la busqueda de 1 a 10 iteracciones para regresiones polinomiales

set.seed (17)
cv.error.10= rep(0 ,10)
for (i in 1:10) {
glm.fit=glm(mpg~poly(horsepower,i),data=Auto)
cv.error.10[i]=cv.glm (Auto,glm.fit,K=10)$delta [1] }
cv.error.10

 [1] 50.37182 50.65488 32.46844 32.49507 27.20315 27.43767 24.16583 24.35135 23.54375 23.80250

Como principio el tiempo de computo de LOOCV para modelos cuadraticos debe ser mas rapido que k-fold CV

5.3.4 The Bootstrap

Aplicable a casi todas las situaciones, no es complicado y tiene solo dos pasos

1: creamos la funcion que deseamos utilizadar en la estadisticas, 2: utilizamos la funcion boot que esta en la libreria boot

Utilizando el dataset de Portfolio en ISLR

library(ISLR)
dim(Portfolio)

[1] 100   2

Paso 1: creamos la funcion alpha.fn - el cual ingresa la data x, y, asi como el vector de las observaciones para estimar el ??

alpha.fn=function(data ,index){
X=data$X [index]
Y=data$Y [index]
return((var(Y)-cov(X,Y))/(var(X)+var(Y)-2*cov(X,Y)))}

Entonces alfa, colocamos la data, y el indice de 1 a 100 obsevaciones

alpha.fn(Portfolio,1:100)

[1] 0.5758321

Alfa de Portfolio de las cien observaciones es 0.5758

Utilizamos aqui la funcion sample para seleccionar un grupo de observaciones aleatoreamente de 1 a 100, y se reemplazan. Con la finalidad de reconstruir un nuevo bootstrap dataset y recalcular el alfa

set.seed (1)
alpha.fn(Portfolio,sample(100,100,replace=T))

[1] 0.5963833

El nuevo alfa el 0.5963

Utilizando la funcion boot - Realizamos un bootstrap analisis con calculo de alfa por 1000 veces, y computando la desviacion standard.

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.08861826

El resultado es que utilizando la data original es decir Portfolio ?? = 0.5758, y la desviacion standard SE(^??) es 0.0886.

Ahora estimamos la exactitud del modelos de la regresion lineal para el dataset de Auto

Utilizando la funcion creada para boot

boot.fn=function(data,index)
return (coef(lm(mpg~horsepower,data=data,subset =index)))
boot.fn(Auto,1:397)

(Intercept)  horsepower 
 17.8076120   0.1108047 

Creando bootstraps para auto

set.seed(1)
boot.fn(Auto,sample(397,397, replace =T))

(Intercept)  horsepower 
 17.4674883   0.1125574 

boot.fn(Auto,sample(397,397, replace =T))

(Intercept)  horsepower 
 16.7417737   0.1315611 

E igual como hicimos para portfolio aqui realizamos un boot para Auto con 1000 corridas de desviacion standards

boot(Auto,boot.fn,1000)

ORDINARY NONPARAMETRIC BOOTSTRAP


Call:
boot(data = Auto, statistic = boot.fn, R = 1000)


Bootstrap Statistics :
      original       bias    std. error
t1* 17.8076120 0.0011771179 0.581449335
t2*  0.1108047 0.0001112272 0.009547989

Estimade de SE Bo para el boostrap es de 0.58, y el SE B1 es de 0.00954

Igualmente se puede utilizar summary para ver las SE

summary(lm(mpg~horsepower,data=Auto))$coef

              Estimate Std. Error  t value     Pr(>|t|)
(Intercept) 17.8076120  0.7112886 25.03571 1.542017e-83
horsepower   0.1108047  0.0119490  9.27314 1.197223e-18

Pero los resultados dan un poco diferentes al calculo con strapboot

boot.fn=function(data,index)
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* 16.8389455597 -0.0261845064 0.8982609661
t2*  0.1801991803  0.0013580196 0.0560723042
t3* -0.0007355174 -0.0000106461 0.0005759085

summary(lm(mpg~horsepower+I(horsepower^2),data=Auto))$coef

                     Estimate   Std. Error   t value     Pr(>|t|)
(Intercept)     16.8389455597 0.9890343218 17.025643 4.087631e-49
horsepower       0.1801991803 0.0507204803  3.552789 4.273549e-04
I(horsepower^2) -0.0007355174 0.0005224972 -1.407696 1.600095e-01

Y aqui vemos los resultados iguales del uno al otro