library (class)
library(MASS)
Attaching package: <U+393C><U+3E31>MASS<U+393C><U+3E32>
The following object is masked from <U+393C><U+3E31>package:dplyr<U+393C><U+3E32>:
select
library(ISLR)
library(dplyr)
library(caret)
Inicializamos fijando un seed para el random luego partiremos la data por el train en como un subconjunot para el subset de lm Por ultimo calcularemos el mse basado en el promedio de 392 observaciones
mean((mpg -predict(lm.fit ,Auto))[-train ]^2)
[1] 26.14142
Ahora utilizaremos una funcion polinomica para un caso cuadratico, como veremos el valor se disminuye
mean((mpg -predict (lm.fit2 ,Auto))[-train ]^2)
[1] 19.82259
La version cubica del polinomio pero veremos que no reduce ya mucho el algoritmo
mean((mpg -predict (lm.fit3 ,Auto))[-train ]^2)
[1] 19.78252
Volveremos hacer la prueba con un seed diferente seleccionaremos nuevamente otra muestra, veremos que el error nuevamente vuelve a incrementarse
mean((mpg -predict (lm.fit ,Auto))[-train ]^2)
[1] 23.29559
PAra el caso del polinomio cuadratico ocrre lo mimso
mean((mpg -predict (lm.fit2 ,Auto))[-train ]^2)
[1] 18.90124
el mejor ajuste esta nuevamente en la versión cubica
lm.fit3=lm(mpg~poly(horsepower ,3) ,data=Auto ,subset =train )
mean((mpg -predict (lm.fit3 ,Auto))[-train ]^2)
[1] 19.2574
5.3.2 LOOCV
Eb los labs del cap4 utilizamos glm para regresion logistica en esta ocacion no procederemos generando una familia binomial en su generación
coef(glm.fit)
(Intercept) horsepower
39.9358610 -0.1578447
Como veremos lm y glm son identicos por el hecho de no especificar la familia
coef(lm.fit)
(Intercept) horsepower
39.9358610 -0.1578447
A continuacion utilizaremos cv de la libreria boot, como podremos ver tendremos dos valores identicos en en la validacion cruzada
library (boot)
glm.fit=glm(mpg~horsepower ,data=Auto)
cv.err =cv.glm(Auto ,glm.fit)
cv.err$delta
[1] 24.23151 24.23114
A continuacion generaremos el listado de errores computados de un bucle de 5 iteraciones donde probaremos cambiar el orden del polinomio, veremos que si existe mejora al comparar la lineal con cuadratica, pero ninguna cuadratica mejora significativamente despues de la cuadratica
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] 24.23151 19.24821 19.33498 19.42443 19.03321
5.3.3 Kfold validation
K-fold esta implementado en la libreria de boot en CV y generaremos diez particiones para poder probar la muestra, tomemos encuenta que el ciclo implica realizar el algoritmo para expresar la funcion entermino de un polinomio de grado i
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] 24.20520 19.18924 19.30662 19.33799 18.87911 19.02103 18.89609 19.71201 18.95140 19.50196
5.3.4 The Bootstrap
Boostrap basicamente puede aplicarse a cualquier estrctura que contenga la funcion, a continuación definiremos una funcion personalizada donde recibiremos como parametro el data set y el indice que deseamos utilizar para el caclulo, en nuestro caso sera un portafolio y nuestro modelo es relativamente simple una variable libre y una variable dependiente
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)))
}
alpha.fn(Portfolio ,1:100)
[1] 0.5758321
al realizar el sample evaluaremos nuestro portafolio
alpha.fn(Portfolio ,sample (100 ,100 , replace =T))
[1] 0.5963833
Ahora utilizaremos nuestra funcion de boot para registrar todos los alpha generados de la desviación que se genera del portafolio, como podemos observar el algoritmo nos calcula el Standar error y el BIAS asociado
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
Estimating the Accuracy of a Linear Regression Model
Inicial mente veremos la regresion lineal y sus coeficientes
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
La funcion boot es utilizada para determinar el intercepto y la pendiente asociada a la muestra
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
38.7387134 -0.1481952
La siguiente funcion estimará el error de 1000 pruebas aplicadas al data set de Auto sobre la funcion personalizada
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.0296667441 0.860440524
t2* -0.1578447 -0.0003113047 0.007411218
Con la funcion de summary podremos extraer los standar error de los coeficientes de los interceptos
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-81
Como pudimos obervar hay variacion entre el std error y los coeficientes del la LM, esto implica que estamos viendo el ruido de la varianza por consiguiente necesitamos adaptar la funcion para operar sobre el cuadrado
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.0001208339
LS0tDQp0aXRsZTogIkxhYm9yYXRvcmlvIDUuMzogQ3Jvc3MtVmFsaWRhdGlvbiBhbmQgdGhlIEJvb3RzdHJhcA0KICINCm91dHB1dDogaHRtbF9ub3RlYm9vaw0KLS0tDQoNCiANCg0KYGBge3J9DQpsaWJyYXJ5IChjbGFzcykNCmxpYnJhcnkoTUFTUykNCmxpYnJhcnkoSVNMUikNCmxpYnJhcnkoZHBseXIpDQpsaWJyYXJ5KGNhcmV0KQ0KYGBgDQpJbmljaWFsaXphbW9zIGZpamFuZG8gdW4gc2VlZCBwYXJhIGVsIHJhbmRvbQ0KbHVlZ28gcGFydGlyZW1vcyBsYSBkYXRhIHBvciBlbCB0cmFpbiBlbiBjb21vIHVuIHN1YmNvbmp1bm90IHBhcmEgZWwgc3Vic2V0IGRlIGxtDQpQb3IgdWx0aW1vIGNhbGN1bGFyZW1vcyBlbCBtc2UgYmFzYWRvIGVuIGVsIHByb21lZGlvIGRlIDM5MiBvYnNlcnZhY2lvbmVzIA0KDQpgYGB7cn0NCiBzZXQuc2VlZCgxKSANCnRyYWluPXNhbXBsZSAoMzkyLDE5NikNCmxtLmZpdD1sbShtcGd+aG9yc2Vwb3dlcixkYXRhPUF1dG8sc3Vic2V0PXRyYWluKQ0KYXR0YWNoKEF1dG8pDQptZWFuKChtcGcgLXByZWRpY3QobG0uZml0ICxBdXRvKSlbLXRyYWluIF1eMikNCmBgYA0KQWhvcmEgdXRpbGl6YXJlbW9zIHVuYSBmdW5jaW9uIHBvbGlub21pY2EgcGFyYSB1biBjYXNvIGN1YWRyYXRpY28sIGNvbW8gdmVyZW1vcyBlbCB2YWxvciBzZSBkaXNtaW51eWUNCmBgYHtyfQ0KbG0uZml0Mj1sbShtcGd+cG9seShob3JzZXBvd2VyICwyKSAsZGF0YT1BdXRvICxzdWJzZXQgPXRyYWluICkNCm1lYW4oKG1wZyAtcHJlZGljdCAobG0uZml0MiAsQXV0bykpWy10cmFpbiBdXjIpDQpgYGANCg0KDQpMYSB2ZXJzaW9uIGN1YmljYSBkZWwgcG9saW5vbWlvIHBlcm8gdmVyZW1vcyBxdWUgbm8gcmVkdWNlIHlhIG11Y2hvIGVsIGFsZ29yaXRtbw0KYGBge3J9DQpsbS5maXQzPWxtKG1wZ35wb2x5KGhvcnNlcG93ZXIgLDMpICxkYXRhPUF1dG8gLHN1YnNldCA9dHJhaW4gKQ0KbWVhbigobXBnIC1wcmVkaWN0IChsbS5maXQzICxBdXRvKSlbLXRyYWluIF1eMikNCmBgYA0KDQpWb2x2ZXJlbW9zIGhhY2VyIGxhIHBydWViYSBjb24gdW4gc2VlZCBkaWZlcmVudGUgc2VsZWNjaW9uYXJlbW9zIG51ZXZhbWVudGUgb3RyYSBtdWVzdHJhLCB2ZXJlbW9zIHF1ZSBlbCBlcnJvciBudWV2YW1lbnRlIHZ1ZWx2ZSBhIGluY3JlbWVudGFyc2UNCg0KYGBge3J9DQpzZXQuc2VlZCAoMikNCnRyYWluPXNhbXBsZSAoMzkyICwxOTYpDQpsbS5maXQgPWxtKG1wZ35ob3JzZXBvd2VyICxzdWJzZXQgPXRyYWluKQ0KbWVhbigobXBnIC1wcmVkaWN0IChsbS5maXQgLEF1dG8pKVstdHJhaW4gXV4yKQ0KYGBgDQpQQXJhIGVsIGNhc28gZGVsIHBvbGlub21pbyBjdWFkcmF0aWNvIG9jcnJlIGxvIG1pbXNvDQoNCmBgYHtyfQ0KbG0uZml0Mj1sbShtcGd+cG9seShob3JzZXBvd2VyICwyKSAsZGF0YT1BdXRvICxzdWJzZXQgPXRyYWluICkNCm1lYW4oKG1wZyAtcHJlZGljdCAobG0uZml0MiAsQXV0bykpWy10cmFpbiBdXjIpDQpgYGANCg0KZWwgbWVqb3IgYWp1c3RlIGVzdGEgbnVldmFtZW50ZSBlbiBsYSB2ZXJzafNuIGN1YmljYSANCmBgYHtyfQ0KbG0uZml0Mz1sbShtcGd+cG9seShob3JzZXBvd2VyICwzKSAsZGF0YT1BdXRvICxzdWJzZXQgPXRyYWluICkNCm1lYW4oKG1wZyAtcHJlZGljdCAobG0uZml0MyAsQXV0bykpWy10cmFpbiBdXjIpDQpgYGANCiMgNS4zLjIgTE9PQ1YNCkViIGxvcyBsYWJzIGRlbCBjYXA0IHV0aWxpemFtb3MgZ2xtIHBhcmEgcmVncmVzaW9uIGxvZ2lzdGljYSAgZW4gZXN0YSBvY2FjaW9uIG5vIHByb2NlZGVyZW1vcyBnZW5lcmFuZG8gdW5hIGZhbWlsaWEgYmlub21pYWwgZW4gc3UgZ2VuZXJhY2nzbiANCmBgYHtyfQ0KZ2xtLmZpdD1nbG0obXBnfmhvcnNlcG93ZXIgLGRhdGE9QXV0bykNCmNvZWYoZ2xtLmZpdCkNCmBgYA0KQ29tbyB2ZXJlbW9zIGxtIHkgZ2xtIHNvbiBpZGVudGljb3MgcG9yIGVsIGhlY2hvIGRlIG5vIGVzcGVjaWZpY2FyIGxhIGZhbWlsaWENCmBgYHtyfQ0KbG0uZml0ID1sbShtcGd+aG9yc2Vwb3dlciAsZGF0YT1BdXRvKQ0KY29lZihsbS5maXQpDQpgYGANCg0KQSBjb250aW51YWNpb24gdXRpbGl6YXJlbW9zIGN2IGRlIGxhIGxpYnJlcmlhIGJvb3QsIGNvbW8gcG9kcmVtb3MgdmVyIHRlbmRyZW1vcyBkb3MgdmFsb3JlcyBpZGVudGljb3MgZW4gZW4gbGEgdmFsaWRhY2lvbiBjcnV6YWRhDQpgYGB7cn0NCmxpYnJhcnkgKGJvb3QpDQpnbG0uZml0PWdsbShtcGd+aG9yc2Vwb3dlciAsZGF0YT1BdXRvKQ0KY3YuZXJyID1jdi5nbG0oQXV0byAsZ2xtLmZpdCkNCmN2LmVyciRkZWx0YQ0KYGBgDQoNCkEgY29udGludWFjaW9uIGdlbmVyYXJlbW9zIGVsIGxpc3RhZG8gZGUgZXJyb3JlcyBjb21wdXRhZG9zIGRlIHVuIGJ1Y2xlIGRlIDUgaXRlcmFjaW9uZXMgZG9uZGUgcHJvYmFyZW1vcyBjYW1iaWFyIGVsIG9yZGVuIGRlbCBwb2xpbm9taW8sIHZlcmVtb3MgcXVlIHNpIGV4aXN0ZSBtZWpvcmEgYWwgY29tcGFyYXIgbGEgbGluZWFsIGNvbiBjdWFkcmF0aWNhLCBwZXJvIG5pbmd1bmEgY3VhZHJhdGljYSBtZWpvcmEgc2lnbmlmaWNhdGl2YW1lbnRlIGRlc3B1ZXMgZGUgbGEgY3VhZHJhdGljYQ0KDQpgYGB7cn0NCmN2LmVycm9yPXJlcCAoMCw1KQ0KDQpmb3IgKGkgaW4gMTo1KSB7DQogIGdsbS5maXQgPSBnbG0obXBnfnBvbHkoaG9yc2Vwb3dlciAsIGkpLCBkYXRhID0gQXV0bykNCiAgY3YuZXJyb3JbaV0gPSBjdi5nbG0gKEF1dG8gLCBnbG0uZml0KSRkZWx0YVsxXQ0KfQ0KDQpjdi5lcnJvcg0KYGBgDQoNCg0KIyA1LjMuMyBLZm9sZCB2YWxpZGF0aW9uDQoNCkstZm9sZCBlc3RhIGltcGxlbWVudGFkbyBlbiBsYSBsaWJyZXJpYSBkZSBib290IGVuIENWIHkgZ2VuZXJhcmVtb3MgZGlleiBwYXJ0aWNpb25lcyBwYXJhIHBvZGVyIHByb2JhciBsYSBtdWVzdHJhLCB0b21lbW9zIGVuY3VlbnRhIHF1ZSBlbCBjaWNsbyBpbXBsaWNhIHJlYWxpemFyIGVsIGFsZ29yaXRtbyBwYXJhIGV4cHJlc2FyIGxhIGZ1bmNpb24gZW50ZXJtaW5vIGRlIHVuIHBvbGlub21pbyBkZSBncmFkbyBpDQpgYGB7cn0NCnNldC5zZWVkICgxNykNCg0KY3YuZXJyb3IuMTAgPSByZXAgKDAgLDEwKQ0KDQpmb3IgKGkgaW4gMToxMCkgew0KICBnbG0uZml0PWdsbShtcGc/Pz9wb2x5KGhvcnNlcG93ZXIgLGkpLGRhdGE9QXV0bykNCiAgY3YuZXJyb3IuMTBbaV09Y3YuZ2xtIChBdXRvICxnbG0uZml0ICxLPTEwKSAkZGVsdGEgWzFdDQp9DQoNCmN2LmVycm9yLjEwDQpgYGANCg0KIyA1LjMuNCBUaGUgQm9vdHN0cmFwDQpCb29zdHJhcCBiYXNpY2FtZW50ZSBwdWVkZSBhcGxpY2Fyc2UgYSBjdWFscXVpZXIgZXN0cmN0dXJhIHF1ZSBjb250ZW5nYSBsYSBmdW5jaW9uLCBhIGNvbnRpbnVhY2nzbiBkZWZpbmlyZW1vcyB1bmEgZnVuY2lvbiBwZXJzb25hbGl6YWRhIGRvbmRlIHJlY2liaXJlbW9zIGNvbW8gcGFyYW1ldHJvIGVsIGRhdGEgc2V0IHkgZWwgaW5kaWNlIHF1ZSBkZXNlYW1vcyB1dGlsaXphciBwYXJhIGVsIGNhY2x1bG8sIGVuIG51ZXN0cm8gY2FzbyBzZXJhIHVuIHBvcnRhZm9saW8geSBudWVzdHJvIG1vZGVsbyBlcyByZWxhdGl2YW1lbnRlIHNpbXBsZSB1bmEgdmFyaWFibGUgbGlicmUgeSB1bmEgdmFyaWFibGUgZGVwZW5kaWVudGUNCmBgYHtyfQ0KYWxwaGEuZm49ZnVuY3Rpb24gKGRhdGEgLGluZGV4KXsNCiAgWD1kYXRhJFggW2luZGV4XQ0KICBZPWRhdGEkWSBbaW5kZXhdDQogIHJldHVybiAoKHZhcihZKS1jb3YgKFgsWSkpLyh2YXIoWCkrdmFyKFkpIC0yKiBjb3YoWCxZKSkpDQp9DQphbHBoYS5mbihQb3J0Zm9saW8gLDE6MTAwKQ0KYGBgDQoNCmFsIHJlYWxpemFyIGVsIHNhbXBsZSBldmFsdWFyZW1vcyBudWVzdHJvIHBvcnRhZm9saW8gDQpgYGB7cn0NCnNldC5zZWVkICgxKQ0KYWxwaGEuZm4oUG9ydGZvbGlvICxzYW1wbGUgKDEwMCAsMTAwICwgcmVwbGFjZSA9VCkpDQpgYGANCkFob3JhIHV0aWxpemFyZW1vcyBudWVzdHJhIGZ1bmNpb24gZGUgYm9vdCBwYXJhIHJlZ2lzdHJhciB0b2RvcyBsb3MgYWxwaGEgZ2VuZXJhZG9zIGRlIGxhIGRlc3ZpYWNp824gcXVlIHNlIGdlbmVyYSBkZWwgcG9ydGFmb2xpbywgY29tbyBwb2RlbW9zIG9ic2VydmFyIGVsIGFsZ29yaXRtbyBub3MgY2FsY3VsYSBlbCBTdGFuZGFyIGVycm9yIHkgZWwgQklBUyBhc29jaWFkbyANCmBgYHtyfQ0KYm9vdChQb3J0Zm9saW8gLGFscGhhLmZuLFI9MTAwMCkNCg0KYGBgDQoNCiMgRXN0aW1hdGluZyB0aGUgQWNjdXJhY3kgb2YgYSBMaW5lYXIgUmVncmVzc2lvbiBNb2RlbA0KSW5pY2lhbCBtZW50ZSB2ZXJlbW9zIGxhIHJlZ3Jlc2lvbiBsaW5lYWwgeSBzdXMgY29lZmljaWVudGVzDQpgYGB7cn0NCmJvb3QuZm49ZnVuY3Rpb24gKGRhdGEgLGluZGV4ICkNCnJldHVybiAoY29lZihsbShtcGc/Pz9ob3JzZXBvd2VyICxkYXRhPWRhdGEgLHN1YnNldCA9aW5kZXgpKSkNCg0KYm9vdC5mbihBdXRvICwxOjM5MikNCmBgYA0KTGEgZnVuY2lvbiBib290IGVzIHV0aWxpemFkYSBwYXJhIGRldGVybWluYXIgZWwgaW50ZXJjZXB0byB5IGxhIHBlbmRpZW50ZSBhc29jaWFkYSBhIGxhIG11ZXN0cmENCmBgYHtyfQ0Kc2V0LnNlZWQgKDEpDQpib290LmZuKEF1dG8gLHNhbXBsZSAoMzkyICwzOTIgLCByZXBsYWNlID1UKSkNCmBgYA0KDQpgYGB7cn0NCmJvb3QuZm4oQXV0byAsc2FtcGxlICgzOTIgLDM5MiAsIHJlcGxhY2UgPVQpKQ0KYGBgDQpMYSBzaWd1aWVudGUgZnVuY2lvbiBlc3RpbWFy4SBlbCBlcnJvciBkZSAxMDAwIHBydWViYXMgYXBsaWNhZGFzIGFsIGRhdGEgc2V0IGRlIEF1dG8gc29icmUgbGEgZnVuY2lvbiBwZXJzb25hbGl6YWRhDQpgYGB7cn0NCmJvb3QoQXV0byAsYm9vdC5mbiAsMTAwMCkNCmBgYA0KDQpDb24gbGEgZnVuY2lvbiBkZSBzdW1tYXJ5IHBvZHJlbW9zIGV4dHJhZXIgbG9zIHN0YW5kYXIgZXJyb3IgZGUgbG9zIGNvZWZpY2llbnRlcyBkZSBsb3MgaW50ZXJjZXB0b3MNCmBgYHtyfQ0Kc3VtbWFyeSAobG0obXBnPz8/aG9yc2Vwb3dlciAsZGF0YT1BdXRvKSkkY29lZg0KYGBgDQpDb21vIHB1ZGltb3Mgb2JlcnZhciBoYXkgdmFyaWFjaW9uIGVudHJlIGVsIHN0ZCBlcnJvciB5IGxvcyBjb2VmaWNpZW50ZXMgZGVsIGxhIExNLCBlc3RvIGltcGxpY2EgcXVlIGVzdGFtb3MgdmllbmRvIGVsIHJ1aWRvIGRlIGxhIHZhcmlhbnphIHBvciBjb25zaWd1aWVudGUgbmVjZXNpdGFtb3MgYWRhcHRhciBsYSBmdW5jaW9uIHBhcmEgb3BlcmFyIHNvYnJlIGVsIGN1YWRyYWRvDQpgYGB7cn0NCmJvb3QuZm4gPSBmdW5jdGlvbiAoZGF0YSAsaW5kZXggKQ0KICAgICAgICAgIGNvZWZmaWNpZW50cyhsbShtcGc/Pz9ob3JzZXBvd2VyICtJKCBob3JzZXBvd2VyIF4yKSAsZGF0YT1kYXRhICwNCiAgICAgICAgICBzdWJzZXQgPWluZGV4KSkNCg0Kc2V0LnNlZWQgKDEpDQpib290KEF1dG8gLGJvb3QuZm4gLDEwMDApDQpgYGANCg0K