a) Regresión Lineal Clásica

# Lm
ajus.lm <- lm(medv ~  ., data = BostonHousing)

tabla <- kable(as.data.frame(summary(ajus.lm)$coefficients))


for(i in 1:nrow(summary(ajus.lm)$coefficients)){
  if(as.data.frame(summary(ajus.lm)$coefficients)$`Pr(>|t|)`[i]>.05){
    tabla <- row_spec(tabla,i,background = "#e76f51", bold = T)
  }
}

tabla

# indus y age no significativas

indus y age no resultaron significativas en este modelo, por lo que se decide sacarlas y ajustar un nuevo modelo de regresión lineal múltiple.

ajus.lm <- lm(medv ~  crim+zn+chas+nox+rm+dis+rad+tax+ptratio+b+lstat, data = BostonHousing)

tabla <- kable(as.data.frame(summary(ajus.lm)$coefficients))

for(i in 1:nrow(summary(ajus.lm)$coefficients)){
  if(as.data.frame(summary(ajus.lm)$coefficients)$`Pr(>|t|)`[i]>.05){
    tabla <- row_spec(tabla,i,background = "#e76f51", bold = T)
  }
}

tabla

# # - - - - - - - ECM - - - - - - - - 
# ec <- c()
# for(i in 1:nrow(BostonHousing)){
#   ajus.lm_CV <- lm(medv ~  crim+zn+chas+nox+rm+dis+rad+tax+ptratio+b+lstat, data = BostonHousing)
#   ec[i] <- (predict(ajus.lm_CV, BostonHousing[i,])-BostonHousing$medv[i])^2
#   
#   # Mostrar el proceso
#   cat("Proceso: ", i, "/", nrow(BostonHousing), " (", round(100*i/nrow(BostonHousing),2), "%)\n")
# }
# 
# # Calculo de error cuadratico medio
# ecm_cv_LM <- mean(ec)
# ecm_cv_LM



mae<-c()
for(i in 1:nrow(BostonHousing))
{
  ajus.lm <- lm(medv ~  crim+zn+chas+nox+rm+dis+rad+tax+ptratio+b+lstat, data = BostonHousing[-i,])
  mae[i]<-abs(predict(ajus.lm,BostonHousing[i,])-BostonHousing$medv[i])
}

MAE_lm <- mean(mae)










# Plotear la capacidad predictiva del modelo
graf <- data.frame(Observados=BostonHousing$medv)
graf$Predichos <- predict(ajus.lm, newdata=BostonHousing)

ggplot(graf, aes(x=Observados,y=Predichos))+
  geom_point(alpha=.2)+
  geom_abline(intercept = 0,
              slope = 1,
              color='#4cadb2',
              lwd=1.5)+
  ggtitle("Valores de medv vs. predichos por modelo lineal múltiple")+
  theme_minimal()

A través de la técnica de validación cruzada (LOOCV), el Error Medio Absoluto estimado fue 3.37.

b) Regresión Lineal Robusta

Con las variables indus y age el modelo no converge, por lo que se ajustó un modelo robusto directamente sin esas dos variables.

# Robustlm
ajus.lmrob <- lmrob(medv ~  ., data = BostonHousing)
# summary(ajus.lmrob)
# Con indus y age no converge
ajus.lmrob <- lmrob(medv ~  crim+zn+chas+nox+rm+dis+rad+tax+ptratio+b+lstat, data = BostonHousing)
summary(ajus.lmrob)

## 
## Call:
## lmrob(formula = medv ~ crim + zn + chas + nox + rm + dis + rad + tax + ptratio + 
##     b + lstat, data = BostonHousing)
##  \--> method = "MM"
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -14.80831  -1.86264  -0.04985   2.31354  33.13749 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  13.515876   9.064954   1.491 0.136599    
## crim         -0.109638   0.034343  -3.192 0.001501 ** 
## zn            0.033619   0.014536   2.313 0.021145 *  
## chas1         1.263976   0.646051   1.956 0.050973 .  
## nox         -10.044991   2.639202  -3.806 0.000159 ***
## rm            5.429161   1.224246   4.435 1.14e-05 ***
## dis          -0.830972   0.169640  -4.898 1.31e-06 ***
## rad           0.178821   0.054318   3.292 0.001066 ** 
## tax          -0.011710   0.002519  -4.649 4.28e-06 ***
## ptratio      -0.755307   0.089350  -8.453 3.20e-16 ***
## b             0.011174   0.002222   5.030 6.89e-07 ***
## lstat        -0.332847   0.056731  -5.867 8.13e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Robust residual standard error: 3.003 
## Multiple R-squared:  0.8305, Adjusted R-squared:  0.8267 
## Convergence in 44 IRWLS iterations
## 
## Robustness weights: 
##  11 observations c(162,365,366,368,369,370,371,372,373,375,413)
##   are outliers with |weight| = 0 ( < 0.0002); 
##  56 weights are ~= 1. The remaining 439 ones are summarized as
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## 0.001491 0.838600 0.945300 0.867600 0.985200 0.999000 
## Algorithmic parameters: 
##        tuning.chi                bb        tuning.psi        refine.tol 
##         1.548e+00         5.000e-01         4.685e+00         1.000e-07 
##           rel.tol         scale.tol         solve.tol       eps.outlier 
##         1.000e-07         1.000e-10         1.000e-07         1.976e-04 
##             eps.x warn.limit.reject warn.limit.meanrw 
##         1.293e-09         5.000e-01         5.000e-01 
##      nResample         max.it       best.r.s       k.fast.s          k.max 
##            500             50              2              1            200 
##    maxit.scale      trace.lev            mts     compute.rd fast.s.large.n 
##            200              0           1000              0           2000 
##                   psi           subsampling                   cov 
##            "bisquare"         "nonsingular"         ".vcov.avar1" 
## compute.outlier.stats 
##                  "SM" 
## seed : int(0)

tabla <- kable(as.data.frame(summary(ajus.lmrob)$coefficients))

for(i in 1:nrow(summary(ajus.lmrob)$coefficients)){
  if(as.data.frame(summary(ajus.lmrob)$coefficients)$`Pr(>|t|)`[i]>.05){
    tabla <- row_spec(tabla,i,background = "#e76f51", bold = T)
  }
}

tabla

Si bien la variable chas no resultó significativa, el p-valor asociado a su coeficiente es apenas superior a 0.05, por lo que se decidió mantener la variable en el modelo

mae<-c()
for(i in 1:nrow(BostonHousing))
{
  ajus.lmrob <- lmrob(medv ~  crim+zn+chas+nox+rm+dis+rad+tax+ptratio+b+lstat, data = BostonHousing[-i,])
  mae[i]<-abs(predict(ajus.lmrob,BostonHousing[i,])-BostonHousing$medv[i])
}

MAE_lmrob <- mean(mae)



# Plotear la capacidad predictiva del modelo
graf <- data.frame(Observados=BostonHousing$medv)
graf$Predichos <- predict(ajus.lmrob, newdata=BostonHousing)

ggplot(graf, aes(x=Observados,y=Predichos))+
  geom_point(alpha=.2)+
  geom_abline(intercept = 0,
              slope = 1,
              color='#4cadb2',
              lwd=1.5)+
  ggtitle("Valores de medv vs. predichos por modelo lineal múltiple robusto")+
  theme_minimal()

A través de la técnica de validación cruzada (LOOCV), el Error Medio Absoluto estimado fue 3.22.

c) Regresión Ridge

# Penalización por Ridge

set.seed(788)
lambda_vals <-rev(seq(0, 6, by = 0.01))
modelo_boston <- glmnet(x=BostonHousing[,c(1:13)],y=BostonHousing[,14],alpha=0,lambda=lambda_vals)
plot(modelo_boston,xvar="lambda",label=TRUE)

cv_bostonridge<-cv.glmnet(x=data.matrix(BostonHousing[,1:13]),y=BostonHousing[,14],alpha=0,type.measure = "mse",lambda=lambda_vals)
lambda_minimo <- cv_bostonridge$lambda.min
plot(cv_bostonridge)

print(cv_bostonridge)

## 
## Call:  cv.glmnet(x = data.matrix(BostonHousing[, 1:13]), y = BostonHousing[,      14], lambda = lambda_vals, type.measure = "mse", alpha = 0) 
## 
## Measure: Mean-Squared Error 
## 
##     Lambda Index Measure    SE Nonzero
## min   0.09   592   23.40 2.468      13
## 1se   3.03   298   25.86 2.633      13

#lambda.min=0.78



# Plotear la capacidad predictiva del modelo
graf <- data.frame(Observados=BostonHousing$medv)
graf$Predichos <- predict(modelo_boston, newx = data.matrix(BostonHousing[,1:13]), s=lambda_minimo)

ggplot(graf, aes(x=Observados,y=Predichos))+
  geom_point(alpha=.2)+
  geom_abline(intercept = 0,
              slope = 1,
              color='#4cadb2',
              lwd=1.5)+
  ggtitle("Valores de medv vs. predichos por modelo lineal múltiple robusto")+
  theme_minimal()

mae <- c()
for(i in 1:nrow(BostonHousing))
{
  modelo_ridge <- glmnet(x=data.matrix(BostonHousing[-i,1:13]),y=BostonHousing[-i,14], alpha=0, lambda=lambda_minimo)
  mae[i]<-abs(predict(modelo_ridge,newx=data.matrix(BostonHousing[i,1:13]))-BostonHousing$medv[i])
}

MAE_ridge <- mean(mae)

A través de la técnica de validación cruzada (LOOCV), el MAE estimado fue 3.365.

d) Regresión LASSO

# Penalización por Lasso

modelo_bostonlasso <- glmnet(x=BostonHousing[,c(1:13)],y=BostonHousing[,14],alpha=1,lambda=lambda_vals)
plot(modelo_bostonlasso, label=TRUE)

cv_boston <- cv.glmnet(x=data.matrix(BostonHousing[,c(1:13)]),y=BostonHousing[,14],alpha=1,type.measure = "mae", lambda=lambda_vals)
lambda_minimo_lasso <- cv_boston$lambda.min
lambda_minimo_lasso

## [1] 0.08

## Lambda minimo = 0.08

plot(cv_boston)

coef_blasso <- coef(modelo_bostonlasso, s=lambda_minimo_lasso)
coef_blasso

## 14 x 1 sparse Matrix of class "dgCMatrix"
##                        s1
## (Intercept)  28.340058177
## crim         -0.080287391
## zn            0.033454814
## indus         .          
## chas          2.618556883
## nox         -14.344875200
## rm            3.983208293
## age           .          
## dis          -1.219107254
## rad           0.168628170
## tax          -0.006327870
## ptratio      -0.899967301
## b             0.008539745
## lstat        -0.522399241

mae<-c()
for(i in 1:nrow(BostonHousing))
{
  modelo_lasso <- glmnet(x=data.matrix(BostonHousing[-i,1:13]),y=BostonHousing[-i,14], alpha=1, lambda=lambda_minimo_lasso)
  mae[i]<-abs(predict(modelo_lasso,newx=data.matrix(BostonHousing[i,1:13]))-BostonHousing$medv[i])
}

MAE_lasso <- mean(mae)

A través de la técnica de validación cruzada (LOOCV), el MAE estimado fue 3.365.

e) Regresión GAM (Generalized Additive Model)

# GAM
#  (Deberíamos meter solo las variables elegidas con lasso y evaluar interacciones ¿con PPR?)
ajus.gam <- gam::gam(medv ~  s(crim) + s(zn) + s(indus) + 
                             s(nox) + s(rm) + s(age) + s(dis) + 
                             s(rad) + s(tax) + s(ptratio) + s(b) +
                             s(lstat), data = BostonHousing)
# la variable chas no entra por ser factor


lineales <- as.data.frame(summary(ajus.gam)$parametric.anova)
tabla <- kable(lineales)

for(i in 1:nrow(lineales)){
  if(!is.na(lineales$`Pr(>F)`[i])&lineales$`Pr(>F)`[i]>.05){
    tabla <- row_spec(tabla,i,background = "#e76f51", bold = T)

  }
}

tabla %>%
   add_header_above(c("Efectos lineales" = 6))

Todos los efectos lineales resultaron significativos.

nolineales <- as.data.frame(summary(ajus.gam)$anova)
tabla <- kable(nolineales)

for(i in 1:nrow(nolineales)){
  if(!is.na(nolineales$`Pr(F)`[i])&nolineales$`Pr(F)`[i]>.05){
    tabla <- row_spec(tabla,i,background = "#e76f51", bold = T)

  }
}

tabla %>%
   add_header_above(c("Efectos no lineales" = 4))

Los efectos no lineales de gran parte de las variables resultaron significativos.

plot(ajus.gam)

mae<-c()
for(i in 1:nrow(BostonHousing))
{
  modelo_gam <- gam::gam(medv ~  s(crim) + s(zn) + s(indus) + 
                             s(nox) + s(rm) + s(age) + s(dis) + 
                             s(rad) + s(tax) + s(ptratio) + s(b) +
                             s(lstat), data = BostonHousing[-i,])
  mae[i]<-abs(predict(modelo_gam,BostonHousing[i,])-BostonHousing$medv[i])
}

MAE_gam <- mean(mae)



# Plotear la capacidad predictiva del modelo
graf <- data.frame(Observados=BostonHousing$medv)
graf$Predichos <- ajus.gam$fitted.values

ggplot(graf, aes(x=Observados,y=Predichos))+
  geom_point(alpha=.2)+
  geom_abline(intercept = 0,
              slope = 1,
              color='#4cadb2',
              lwd=1.5)+
  ggtitle("Valores de medv vs. predichos por modelo GAM")+
  theme_minimal()

A través de la técnica de validación cruzada (LOOCV), el MAE estimado fue 2.54.

f) Regresión PPR (Projection Pursuit Regression)

# PPR
ajus.ppr <- ppr(medv ~ crim+zn+indus+nox+rm+age+dis+rad+tax+ptratio+b+lstat, data=BostonHousing, 
                nterms = 3, 
                max.terms = 4,
                trace=F)



tabla <- kable(as.data.frame(ajus.ppr$alpha))

tabla %>%
   add_header_above(c("Vectores de dirección de proyección" = 4))

tabla <- kable(as.data.frame(ajus.ppr$beta))

tabla %>%
   add_header_above(c("Coeficientes de los términos de ridge" = 2))

plot(ajus.ppr)

A partir de los gráficos (particularmente el segundo), se infiere una interacción entre las variables de mayor peso en el segundo término (nox, rm y lstat).

mae<-c()
for(i in 1:nrow(BostonHousing))
{
  ajus.ppr <- ppr(medv ~ crim+zn+indus+nox+rm+age+dis+rad+tax+ptratio+b+lstat, data=BostonHousing[-i,], 
                nterms = 3, 
                max.terms = 4,
                trace=F)
  mae[i]<-abs(predict(ajus.ppr,BostonHousing[i,])-BostonHousing$medv[i])
}

MAE_ppr <- mean(mae)




# Plotear la capacidad predictiva del modelo
graf <- data.frame(Observados=BostonHousing$medv)
graf$Predichos <- predict(ajus.ppr, newdata=BostonHousing)

ggplot(graf, aes(x=Observados,y=Predichos))+
  geom_point(alpha=.2)+
  geom_abline(intercept = 0,
              slope = 1,
              color='#4cadb2',
              lwd=1.5)+
  ggtitle("Valores de medv vs. predichos por modelo GAM")+
  theme_minimal()

A través de la técnica de validación cruzada (LOOCV), el MAE estimado fue 2.612.

g) Regresión mediante ANN (Artificial Neural Networks)

# ANN
set.seed(195)
ajus.nnet <- nnet(medv ~ ., # Formula
                  data=BostonHousing, # Datos 
                  size=5, # Cantidad de neuronas 
                  decay=0.5, 
                  maxit=3000, # Epocas 
                  linout=1, 
                  trace=F) # No muestra por donde va
plotnet(ajus.nnet)

mae<-c()
for(i in 1:nrow(BostonHousing))
{
  
  set.seed(195)
  ajus.nnet <- nnet(medv ~ ., # Formula
                  data=BostonHousing[-i,], # Datos 
                  size=5, # Cantidad de neuronas 
                  decay=0.5, 
                  maxit=3000, # Epocas 
                  linout=1, 
                  trace=F) # No muestra por donde va
  mae[i]<-abs(predict(ajus.nnet,BostonHousing[i,])-BostonHousing$medv[i])

}

MAE_nnet <- mean(mae)



# Plotear la capacidad predictiva del modelo
graf <- data.frame(Observados=BostonHousing$medv)
graf$Predichos <- predict(ajus.nnet, newdata=BostonHousing)

ggplot(graf, aes(x=Observados,y=Predichos))+
  geom_point(alpha=.2)+
  geom_abline(intercept = 0,
              slope = 1,
              color='#4cadb2',
              lwd=1.5)+
  ggtitle("Valores de medv vs. predichos por modelo GAM")+
  theme_minimal()

A través de la técnica de validación cruzada (LOOCV), el MAE estimado fue 2.683.