Capitulo 7: Mas alla de la linealidad
Regresión polinomial
library(ISLR)
attach(Wage)
head(Wage)
## year age sex maritl race education
## 231655 2006 18 1. Male 1. Never Married 1. White 1. < HS Grad
## 86582 2004 24 1. Male 1. Never Married 1. White 4. College Grad
## 161300 2003 45 1. Male 2. Married 1. White 3. Some College
## 155159 2003 43 1. Male 2. Married 3. Asian 4. College Grad
## 11443 2005 50 1. Male 4. Divorced 1. White 2. HS Grad
## 376662 2008 54 1. Male 2. Married 1. White 4. College Grad
## region jobclass health health_ins logwage
## 231655 2. Middle Atlantic 1. Industrial 1. <=Good 2. No 4.318
## 86582 2. Middle Atlantic 2. Information 2. >=Very Good 2. No 4.255
## 161300 2. Middle Atlantic 1. Industrial 1. <=Good 1. Yes 4.875
## 155159 2. Middle Atlantic 2. Information 2. >=Very Good 1. Yes 5.041
## 11443 2. Middle Atlantic 2. Information 1. <=Good 1. Yes 4.318
## 376662 2. Middle Atlantic 2. Information 2. >=Very Good 1. Yes 4.845
## wage
## 231655 75.04
## 86582 70.48
## 161300 130.98
## 155159 154.69
## 11443 75.04
## 376662 127.12
fit=lm(wage~poly(age,3),data=Wage)
coef(summary(fit))
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 111.7 0.7291 153.211 0.000e+00
## poly(age, 3)1 447.1 39.9335 11.195 1.571e-28
## poly(age, 3)2 -478.3 39.9335 -11.978 2.512e-32
## poly(age, 3)3 125.5 39.9335 3.143 1.687e-03
fit1=lm(wage~poly(age,4),data=Wage)
coef(summary(fit1))
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 111.70 0.7287 153.283 0.000e+00
## poly(age, 4)1 447.07 39.9148 11.201 1.485e-28
## poly(age, 4)2 -478.32 39.9148 -11.983 2.356e-32
## poly(age, 4)3 125.52 39.9148 3.145 1.679e-03
## poly(age, 4)4 -77.91 39.9148 -1.952 5.104e-02
fit2=lm(wage~poly(age,3,raw=T),data=Wage)
coef(summary(fit2))
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -7.524e+01 2.218e+01 -3.392 7.032e-04
## poly(age, 3, raw = T)1 1.019e+01 1.605e+00 6.348 2.511e-10
## poly(age, 3, raw = T)2 -1.680e-01 3.686e-02 -4.559 5.357e-06
## poly(age, 3, raw = T)3 8.495e-04 2.702e-04 3.143 1.687e-03
agelims=range(age)
age.grid=seq(from=agelims[1],to=agelims[2])
preds=predict(fit,newdata=list(age=age.grid),se=TRUE)
se.bands=cbind(preds$fit+2*preds$se.fit,preds$fit-2*preds$se.fit)
plot(age,wage,xlim=agelims,cex=.5,col="darkgrey")
title("Polinomio de grado 3")
lines(age.grid,preds$fit,lwd=2,col="blue")
matlines(age.grid,se.bands,lwd=1,col="blue",lty=3)
fit.1=lm(wage~age,data=Wage)
fit.2=lm(wage~poly(age,2),data=Wage)
fit.3=lm(wage~poly(age,3),data=Wage)
fit.4=lm(wage~poly(age,4),data=Wage)
fit.5=lm(wage~poly(age,5),data=Wage)
anova(fit.1,fit.2,fit.3,fit.4,fit.5)
## Analysis of Variance Table
##
## Model 1: wage ~ age
## Model 2: wage ~ poly(age, 2)
## Model 3: wage ~ poly(age, 3)
## Model 4: wage ~ poly(age, 4)
## Model 5: wage ~ poly(age, 5)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 2998 5022216
## 2 2997 4793430 1 228786 143.59 <2e-16 ***
## 3 2996 4777674 1 15756 9.89 0.0017 **
## 4 2995 4771604 1 6070 3.81 0.0510 .
## 5 2994 4770322 1 1283 0.80 0.3697
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
coef(summary(fit.5))
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 111.70 0.7288 153.2780 0.000e+00
## poly(age, 5)1 447.07 39.9161 11.2002 1.491e-28
## poly(age, 5)2 -478.32 39.9161 -11.9830 2.368e-32
## poly(age, 5)3 125.52 39.9161 3.1446 1.679e-03
## poly(age, 5)4 -77.91 39.9161 -1.9519 5.105e-02
## poly(age, 5)5 -35.81 39.9161 -0.8972 3.697e-01
fit.1=lm(wage~education+age,data=Wage)
fit.2=lm(wage~education+poly(age,2),data=Wage)
fit.3=lm(wage~education+poly(age,3),data=Wage)
anova(fit.1,fit.2,fit.3)
## Analysis of Variance Table
##
## Model 1: wage ~ education + age
## Model 2: wage ~ education + poly(age, 2)
## Model 3: wage ~ education + poly(age, 3)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 2994 3867992
## 2 2993 3725395 1 142597 114.70 <2e-16 ***
## 3 2992 3719809 1 5587 4.49 0.034 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
fit=glm(I(wage>250)~poly(age,4),data=Wage,family=binomial)
preds=predict(fit,newdata=list(age=age.grid),se=T)
pfit=exp(preds$fit)/(1+exp(preds$fit))
se.bands.logit = cbind(preds$fit+2*preds$se.fit, preds$fit-2*preds$se.fit)
se.bands = exp(se.bands.logit)/(1+exp(se.bands.logit))
preds=predict(fit,newdata=list(age=age.grid),type="response",se=T)
plot(age,I(wage>250),xlim=agelims,type="n",ylim=c(0,.2))
points(jitter(age), I((wage>250)/5),cex=.5,pch="|",col="darkgrey")
lines(age.grid,pfit,lwd=2, col="blue")
matlines(age.grid,se.bands,lwd=1,col="blue",lty=3)
Funciones paso
table(cut(age,4))
##
## (17.9,33.5] (33.5,49] (49,64.5] (64.5,80.1]
## 750 1399 779 72
fit=lm(wage~cut(age,4),data=Wage)
coef(summary(fit))
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 94.158 1.476 63.790 0.000e+00
## cut(age, 4)(33.5,49] 24.053 1.829 13.148 1.982e-38
## cut(age, 4)(49,64.5] 23.665 2.068 11.443 1.041e-29
## cut(age, 4)(64.5,80.1] 7.641 4.987 1.532 1.256e-01
Splines
library(splines)
# Spline funciones base
fit=lm(wage~bs(age,knots=c(25,40,60)),data=Wage)
coef(summary(fit))
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 60.49 9.460 6.3944 1.863e-10
## bs(age, knots = c(25, 40, 60))1 3.98 12.538 0.3175 7.509e-01
## bs(age, knots = c(25, 40, 60))2 44.63 9.626 4.6364 3.698e-06
## bs(age, knots = c(25, 40, 60))3 62.84 10.755 5.8426 5.691e-09
## bs(age, knots = c(25, 40, 60))4 55.99 10.706 5.2297 1.815e-07
## bs(age, knots = c(25, 40, 60))5 50.69 14.402 3.5196 4.387e-04
## bs(age, knots = c(25, 40, 60))6 16.61 19.126 0.8682 3.853e-01
pred=predict(fit,newdata=list(age=age.grid),se=T)
plot(age,wage,col="gray")
lines(age.grid,pred$fit,lwd=2)
lines(age.grid,pred$fit+2*pred$se,lty="dashed")
lines(age.grid,pred$fit-2*pred$se,lty="dashed")
# Spline natural
fit2=lm(wage~ns(age,df=4),data=Wage)
coef(summary(fit2))
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 58.556 4.235 13.827 3.394e-42
## ns(age, df = 4)1 60.462 4.190 14.430 1.108e-45
## ns(age, df = 4)2 41.963 4.372 9.597 1.669e-21
## ns(age, df = 4)3 97.020 10.386 9.341 1.800e-20
## ns(age, df = 4)4 9.773 8.657 1.129 2.590e-01
pred2=predict(fit2,newdata=list(age=age.grid),se=T)
lines(age.grid, pred2$fit,col="red",lwd=2)
# Suavizacion por Splines
plot(age,wage,xlim=agelims,cex=.5,col="darkgrey")
title("Suavizacion por Splines")
fit=smooth.spline(age,wage,df=16)
fit2=smooth.spline(age,wage,cv=TRUE)
## Warning: cross-validation with non-unique 'x' values seems doubtful
fit2$df
## [1] 6.795
lines(fit,col="red",lwd=2)
lines(fit2,col="blue",lwd=2)
legend("topright",legend=c("16 DF","6.8 DF"),col=c("red","blue"),lty=1,lwd=2,cex=.8)
Regresion local
plot(age,wage,xlim=agelims,cex=.5,col="darkgrey")
title("Regresion local")
fit=loess(wage~age,span=.2,data=Wage)
fit2=loess(wage~age,span=.5,data=Wage)
lines(age.grid,predict(fit,data.frame(age=age.grid)),col="red",lwd=2)
lines(age.grid,predict(fit2,data.frame(age=age.grid)),col="blue",lwd=2)
legend("topright",legend=c("Span = 0.2","Span = 0.5"),col=c("red","blue"),lty=1,lwd=2,cex=.8)
Modelos Aditivos Generalizados (GAM)
gam1=lm(wage~ns(year,4)+ns(age,5)+education,data=Wage)
library(gam)
## Warning: package 'gam' was built under R version 3.1.2
## Loaded gam 1.09.1
# s() Smoothing Spline
gam.m3=gam(wage~s(year,4)+s(age,5)+education,data=Wage)
par(mfrow=c(1,3))
plot(gam.m3, se=TRUE,col="blue")
plot.gam(gam1, se=TRUE, col="red")
gam.m1=gam(wage~s(age,5)+education,data=Wage)
gam.m2=gam(wage~year+s(age,5)+education,data=Wage)
anova(gam.m1,gam.m2,gam.m3,test="F")
## Analysis of Deviance Table
##
## Model 1: wage ~ s(age, 5) + education
## Model 2: wage ~ year + s(age, 5) + education
## Model 3: wage ~ s(year, 4) + s(age, 5) + education
## Resid. Df Resid. Dev Df Deviance F Pr(>F)
## 1 2990 3711731
## 2 2989 3693842 1 17889 14.5 0.00014 ***
## 3 2986 3689770 3 4071 1.1 0.34857
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
preds=predict(gam.m2,newdata=Wage)
gam.lo=gam(wage~s(year,df=4)+lo(age,span=0.7)+education,data=Wage)
plot.gam(gam.lo, se=TRUE, col="green")
