Ejercicio 5

Objetivo: Determinar si para el conjunto de datos de los salarios presentados en el Tema 1 es más conveniente utilizar un método muy flexible o poco flexible.

library(ISLR2)
data(Wage)
is(Wage)
## [1] "data.frame" "list"       "oldClass"   "vector"
dim(Wage)
## [1] 3000   11
head(Wage)
##        year age           maritl     race       education             region
## 231655 2006  18 1. Never Married 1. White    1. < HS Grad 2. Middle Atlantic
## 86582  2004  24 1. Never Married 1. White 4. College Grad 2. Middle Atlantic
## 161300 2003  45       2. Married 1. White 3. Some College 2. Middle Atlantic
## 155159 2003  43       2. Married 3. Asian 4. College Grad 2. Middle Atlantic
## 11443  2005  50      4. Divorced 1. White      2. HS Grad 2. Middle Atlantic
## 376662 2008  54       2. Married 1. White 4. College Grad 2. Middle Atlantic
##              jobclass         health health_ins  logwage      wage
## 231655  1. Industrial      1. <=Good      2. No 4.318063  75.04315
## 86582  2. Information 2. >=Very Good      2. No 4.255273  70.47602
## 161300  1. Industrial      1. <=Good     1. Yes 4.875061 130.98218
## 155159 2. Information 2. >=Very Good     1. Yes 5.041393 154.68529
## 11443  2. Information      1. <=Good     1. Yes 4.318063  75.04315
## 376662 2. Information 2. >=Very Good     1. Yes 4.845098 127.11574
table(Wage$region)
## 
##        1. New England    2. Middle Atlantic 3. East North Central 
##                     0                  3000                     0 
## 4. West North Central     5. South Atlantic 6. East South Central 
##                     0                     0                     0 
## 7. West South Central           8. Mountain            9. Pacific 
##                     0                     0                     0

Para empezar, recordemos las variables que previamente concluimos que serían predictores apropiados para el modelo.

  1. year, es el año en el que se registró la información salarial (entre 2003 y 2009).

  2. age, es la edad del trabajador (entre 18 y 80 años).

  3. maritl, es un factor indicando el estado civil del trabajador con niveles 1. Nunca casado, 2. Casado, 3. Viudo, 4. Divorciado y 5. Separado.

  4. race, es un factor indicando la raza del trabajador con niveles 1. Blanco, 2. Negro, 3. Asiático y 4. Otro.

  5. education, es un factor indicando el nivel de educación del trabajador con niveles 1. Menos que secundaria, 2. Secundaria, 3. Comenzó grado, 4. Grado y 5. Postgrado.

  6. jobclass, es un factor indicando el tipo de labor realizada por el trabajador con niveles 1. Industrial y 2. Información.

  7. health, es un factor indicando el nivel de salud del trabajador con niveles 1. Bueno o peor que bueno, 2. Muy bueno o mejor que muy bueno.

  8. health_ins, es un factor indicando si el trabajador tiene un seguro de salud con niveles 1. Si y 2. No.

  9. logwage, es el logaritmo del salario del trabajador (en logaritmos de miles de dolares).

  10. wage, es el salario del trabajador (en miles de dolares).

Notemos que hemos quitado de la lista la variable location porqué esta tiene un unico valor (todos los encuestados son de la misma ciudad). Por otra parte, hemos mencionado que el analisis se hará en base a logwage y no wage por tener una distribución más simétrica.

No podemos usar todos los predictores para una correlación ya que algunas son variables cualitativas nominales (sin orden) o no binarias.

#preparación de librerias y colore para graficar
library(ggplot2)
color_1 <- "turquoise2"
color_2 <- "salmon2"
library(gridExtra)
install.packages("plotly")
## 
## The downloaded binary packages are in
##  /var/folders/jd/zn7_v9ls5n35srxc2v42_ymh0000gn/T//RtmpXB6D7k/downloaded_packages
library(plotly)

Lo que hemos visto anteriormente y se puede reutilizar

En clase hemos observado como diferentes regresiones se ajustan a un conjunto de datos similar al de este problema (Ambos data frames giran en torno a la predicción o inferencia del ingreso de una población), más especificamente, hemos usado una regresion lineal, y dos B-splines para mostrar la correlación entre Education y Seniority con Income.

Con base en eso podemos partir nuestro análisis con una hipotesis: Usar un modelo lineal hará que perdamos mucha información de la correlación, y un método muy flexible causará sobre ajuste; por lo cual, lo mejor será usar un modelo poco flexible.

El análisis

¿Que variables nos importan/sirven? En clase hemos discutido esto previamente. Sin embargo, vamos a retomar la discusión. year no parecia tener alguna correlación relevante, jobclass y health no tienen una importancia obvia (puede ser que haya correlación pero no es tan evidente) y race no puede ser usada para una regresión ya que no podemos convertirla a una variable binaria. Esto nos deja con una variable cuantitativa: age, una cualitativa ordinal: education, y una cualitativa binaria: health_ins

Empecemos por el scatter plot de age

 plot_ly(Wage,x=~age,y=~logwage, type="scatter",
  symbols=20,mode="markers",color=I(color_1))

Ahora, hagamos los box plots para nuestra variable cualitativa ordinal education

plot_ly(Wage,x=~education,y=~logwage, type="box",
       symbols=20,mode="markers",color=I(color_1))
## Warning: 'box' objects don't have these attributes: 'mode'
## Valid attributes include:
## 'alignmentgroup', 'boxmean', 'boxpoints', 'customdata', 'customdatasrc', 'dx', 'dy', 'fillcolor', 'hoverinfo', 'hoverinfosrc', 'hoverlabel', 'hoveron', 'hovertemplate', 'hovertemplatesrc', 'hovertext', 'hovertextsrc', 'ids', 'idssrc', 'jitter', 'legendgroup', 'legendgrouptitle', 'legendrank', 'line', 'lowerfence', 'lowerfencesrc', 'marker', 'mean', 'meansrc', 'median', 'mediansrc', 'meta', 'metasrc', 'name', 'notched', 'notchspan', 'notchspansrc', 'notchwidth', 'offsetgroup', 'opacity', 'orientation', 'pointpos', 'q1', 'q1src', 'q3', 'q3src', 'quartilemethod', 'sd', 'sdsrc', 'selected', 'selectedpoints', 'showlegend', 'stream', 'text', 'textsrc', 'transforms', 'type', 'uid', 'uirevision', 'unselected', 'upperfence', 'upperfencesrc', 'visible', 'whiskerwidth', 'width', 'x', 'x0', 'xaxis', 'xcalendar', 'xhoverformat', 'xperiod', 'xperiod0', 'xperiodalignment', 'xsrc', 'y', 'y0', 'yaxis', 'ycalendar', 'yhoverformat', 'yperiod', 'yperiod0', 'yperiodalignment', 'ysrc', 'key', 'set', 'frame', 'transforms', '_isNestedKey', '_isSimpleKey', '_isGraticule', '_bbox'

y para health_ins

plot_ly(Wage,x=~health_ins,y=~logwage, type="box",
       symbols=20,mode="markers",color=I(color_1))
## Warning: 'box' objects don't have these attributes: 'mode'
## Valid attributes include:
## 'alignmentgroup', 'boxmean', 'boxpoints', 'customdata', 'customdatasrc', 'dx', 'dy', 'fillcolor', 'hoverinfo', 'hoverinfosrc', 'hoverlabel', 'hoveron', 'hovertemplate', 'hovertemplatesrc', 'hovertext', 'hovertextsrc', 'ids', 'idssrc', 'jitter', 'legendgroup', 'legendgrouptitle', 'legendrank', 'line', 'lowerfence', 'lowerfencesrc', 'marker', 'mean', 'meansrc', 'median', 'mediansrc', 'meta', 'metasrc', 'name', 'notched', 'notchspan', 'notchspansrc', 'notchwidth', 'offsetgroup', 'opacity', 'orientation', 'pointpos', 'q1', 'q1src', 'q3', 'q3src', 'quartilemethod', 'sd', 'sdsrc', 'selected', 'selectedpoints', 'showlegend', 'stream', 'text', 'textsrc', 'transforms', 'type', 'uid', 'uirevision', 'unselected', 'upperfence', 'upperfencesrc', 'visible', 'whiskerwidth', 'width', 'x', 'x0', 'xaxis', 'xcalendar', 'xhoverformat', 'xperiod', 'xperiod0', 'xperiodalignment', 'xsrc', 'y', 'y0', 'yaxis', 'ycalendar', 'yhoverformat', 'yperiod', 'yperiod0', 'yperiodalignment', 'ysrc', 'key', 'set', 'frame', 'transforms', '_isNestedKey', '_isSimpleKey', '_isGraticule', '_bbox'

Regresiones lineales

Empecemos por age

# Modelo
linear_fit_age <- lm(logwage~age,data=Wage)

#grafica
plot(Wage$age, Wage$logwage, main = "Regresión Lineal Age", xlab = "Age", ylab = "log wage", pch = 19, col = color_1)
abline(linear_fit_age, col=I(color_2), lwd = 2)

#ECM
predicted_age <- predict(linear_fit_age)
residuals_age <- Wage$logwage - predicted_age
ECM_lineal_age <- mean(residuals_age^2)


cat("ECM para age:", ECM_lineal_age, "\n")
## ECM para age: 0.1178164

ECM para age con la regresion lineal: 0.1178164

Para education,

#Modelo
linear_fit_education <- lm(logwage ~ education, data = Wage)

# Grafica
plot(Wage$education, Wage$logwage, main = "Regresión Lineal Education", xlab = "Education", ylab = "Log Wage", col = color_1, pch = 19)
predicted_values <- predict(linear_fit_education, newdata = data.frame(education = levels(Wage$education)))
lines(1:5, predicted_values, col = color_2, lwd = 2)  

#ECM
residuals_education <- Wage$logwage - predicted_values
ECM_lineal_education <- mean(residuals_education^2)

cat("ECM para education:", ECM_lineal_education, "\n")
## ECM para education: 0.1651646

ECM para education con la regresion lineal: 0.1651646

Para ambas variables se puede observar que la regresión lineal no captura en su toalidad la relación entre las variables, y tiene un error cuadratico medio enorme.

y para health_ins solo ponemos el ECM de la regresión, ya que, al no tener un orden, no hace sentido la gráfica.

#Modelo
linear_fit_health_ins <- lm(logwage ~ health_ins, data = Wage)
#ECM
predicciones <- predict(linear_fit_health_ins)  
residuos <- Wage$logwage - predicciones 
ECM_lineal_health_ins <- mean(residuos^2)
cat ("ECM para health ins:",ECM_lineal_health_ins)
## ECM para health ins: 0.1067801

ECM para health_ins con la regresion lineal: 0.1067801

Ajustes poco flexibles

Ahora probemos con un modelo más flexible que el lineal pero de grado bajo.

library(splines)

#Modelo
modelo_bspline_age <- lm(logwage ~ bs(age, degree = 4), data = Wage)

#grafica
plot(Wage$age, Wage$logwage, main = "B-spline de Grado 4 para Age", xlab = "Age", ylab = "Log Wage", pch = 19, col = color_1)
age_seq <- seq(min(Wage$age), max(Wage$age), length.out = 100) #curva
predicted_bspline_age <- predict(modelo_bspline_age, newdata = data.frame(age = age_seq))
lines(age_seq, predicted_bspline_age, col = color_2, lwd = 2)

#ECM
residuos <- Wage$logwage - predicted_bspline_age
ECM_age_bspline_bajo <- mean(residuos^2)
cat ("ECM para age:",ECM_age_bspline_bajo)
## ECM para age: 0.1560217
# Modelo
modelo_bspline_education <- lm(logwage ~ bs(as.numeric(education), degree = 4), data = Wage)

#Grafica
plot(Wage$education, Wage$logwage, main = "B-spline de Grado 4 para Education", xlab = "Education", ylab = "Log Wage", col = color_1, pch = 19)
education_seq <- seq(min(as.numeric(Wage$education)), max(as.numeric(Wage$education)), length.out = 100)
predicted_bspline_education <- predict(modelo_bspline_education, newdata = data.frame(education = education_seq))
lines(education_seq, predicted_bspline_education, col = color_2, lwd = 2)

#ECM
residuos <- Wage$logwage - predicted_bspline_education
ECM_education_bspline_bajo <- mean(residuos^2)
cat ("ECM para education:",ECM_education_bspline_bajo)
## ECM para education: 0.1481386
#Modelo
modelo_bspline_health_ins <- lm(logwage ~ bs(as.numeric(health_ins), degree = 4), data = Wage)

#ECM
predicciones <- predict(modelo_bspline_health_ins)  
residuos <- Wage$logwage - predicciones 
ECM_healt_ins_bspline_bajo <- mean(residuos^2)
cat("ECM para health_ins", ECM_healt_ins_bspline_bajo, "\n")
## ECM para health_ins 0.1067801

Ajustes muy flexibles

Ajustes más flexibles

Ahora probemos con un modelo más flexible que el lineal pero de grado bajo.

library(splines)

#Modelo
modelo_bspline_age <- lm(logwage ~ bs(age, degree = 8), data = Wage)

#grafica
plot(Wage$age, Wage$logwage, main = "B-spline de Grado 8 para Age", xlab = "Age", ylab = "Log Wage", pch = 19, col = color_1)
age_seq <- seq(min(Wage$age), max(Wage$age), length.out = 100) #curva
predicted_bspline_age <- predict(modelo_bspline_age, newdata = data.frame(age = age_seq))
lines(age_seq, predicted_bspline_age, col = color_2, lwd = 2)

#ECM
residuos <- Wage$logwage - predicted_bspline_age
ECM_age_bspline_alto <- mean(residuos^2)
cat ("ECM para age:",ECM_age_bspline_alto)
## ECM para age: 0.1575472
# Modelo
modelo_bspline_education <- lm(logwage ~ bs(as.numeric(education), degree = 8), data = Wage)

#Grafica
plot(Wage$education, Wage$logwage, main = "B-spline de Grado 8 para Education", xlab = "Education", ylab = "Log Wage", col = color_1, pch = 19)
education_seq <- seq(min(as.numeric(Wage$education)), max(as.numeric(Wage$education)), length.out = 100)
predicted_bspline_education <- predict(modelo_bspline_education, newdata = data.frame(education = education_seq))
## Warning in predict.lm(modelo_bspline_education, newdata = data.frame(education
## = education_seq)): prediction from rank-deficient fit; attr(*, "non-estim") has
## doubtful cases
lines(education_seq, predicted_bspline_education, col = color_2, lwd = 2)

#ECM
residuos <- Wage$logwage - predicted_bspline_education
ECM_education_bspline_alto <- mean(residuos^2)
cat ("ECM para education:",ECM_education_bspline_alto)
## ECM para education: 15.18815
#Modelo
modelo_bspline_health_ins <- lm(logwage ~ bs(as.numeric(health_ins), degree = 8), data = Wage)

#ECM
predicciones <- predict(modelo_bspline_health_ins)  
residuos <- Wage$logwage - predicciones 
ECM_healt_ins_bspline_alto <- mean(residuos^2)
cat("ECM para health_ins", ECM_healt_ins_bspline_alto, "\n")
## ECM para health_ins 0.1067801

Conclusiones

Empecemos por los más relevante. Debemos notar los errores cuadraticos medios para cada uno de los modelos

Para el modelo lineal los valores son:

ECM_lineal_age : 0.1178 ECM_lineal_education : 0.1652 ECM_lineal_health_ins : 0.1068

Para el modelo poco flexible los valores son:

ECM_age_bspline_bajo : 0.156O ECM_education_bspline_bajo : 0.1481 ECM_healt_ins_bspline_bajo : 0.1068

Para el modelo muy flexible los valores son:

ECM_age_bspline_alto : 0.1575 ECM_education_bspline_alto : 15.18815 ECM_healt_ins_bspline_alto : 0.1068

Podemos concluir que necesitamos un modelo poco flexible para evitar el sobre ajuste. La felxibilidad del modelo afecta fuertemente las predicciones que pueden ser generadas en base a los valores cuantitativos, y por ende es mejor mantenerla baja.

Restricciónes

  1. No he considerado todas las variables que podrian llegar a ser predictores aptos
  2. Solo he observado la correlación de las variables individualmente y no como un conjunto