Modelización: modelos de regresión

La raíz del error cuadrático medio

Haremos uso de medidas para saber la eficiencia de nuestro modelo. Nos mide cuanto nos hemos equivocado en nuestra predicción:

dat <- read.csv("../DataSets/rmse.csv")
head(dat)

Ahora queremos saber si los datos que hemos predicho son buenos o no lo són. Para evitar que los valores negativos contrarresten las positivo, elevaremos todo al cuadrado, y luego dividiremos entre el número de muestras. Por último, haremos la raíz cuadrada para volver a las dimensiones reales:

rmse <- sqrt(mean((dat$price - dat$pred)^2))
rmse

## [1] 2.934995

Este valor puede ser mucho en relación al valor absoluto de la realidad (no es lo mismo equivocarse de 3€ sobre un producto que vale 900€ que sobre uno que vale 3€).

plot(dat$price, dat$pred, xlab = "Actual", ylab = "Predicho")
abline(0,1)

Podemos usar una función para calcular automaticamente el error cuadrático medio con raíz:

rmse <- function(actual, predicted){
  return(sqrt(mean((actual - predicted)^2)))
}
rmse(dat$price, dat$pred)

## [1] 2.934995

K nearest neighbors

library(FNN)
library(scales)
library(caret)

## Loading required package: lattice

## Loading required package: ggplot2

library(dummies)

## dummies-1.5.6 provided by Decision Patterns

Usaremos un dataset sobre escuelas:

edu <- read.csv("../DataSets/education.csv")
head(edu)

La variables state y region son variables categoricas y las debemos transformar a numéricas con dummify. Queremos preveer el gasto del año que viene (“expense”) según las otras variables:

dms <- dummy(edu$region, sep = "_")
edu <- cbind(edu, dms)
head(edu,6)

Ahora vamos a estandarizar (rescalar) las variables numéricas al rango 0/1:

edu$urban.s <- rescale(edu$urban)
edu$income.s <- rescale(edu$income)
edu$under18.s <- rescale(edu$under18)

Ahora vamos a crear el caso:

set.seed(2018)
t.ids <- createDataPartition(edu$expense,
                             p = 0.6, 
                             list = F)
tr <- edu[t.ids,]
temp <- edu[-t.ids,]
v.ids <- createDataPartition(temp$expense, p = 0.5, list = F)
val <- temp[v.ids,]
test <- temp[-v.ids,]

Ahora podemos usar el RMSE para evaluar que modelo es mejor:

reg1 <- knn.reg(tr[,7:12], # predictores que usaremos para predecir 
                val[,7:12], # validadores, los predictores que queremos usar para validar
                tr$expense, # salida que tiene que dar en el conjunto de entrenamiento
                k = 1, # número de vecinos
                algorithm = "brute") # algoritmos que se pueden usar
rmse1 <- sqrt(mean((reg1$pred - val$expense)^2))
rmse1

## [1] 57.5465

Provamos otra k:

reg2 <- knn.reg(tr[,7:12], val[,7:12], tr$expense, k = 2, algorithm = "brute")
rmse2 <- rmse(val$expense, reg2$pred)
rmse2

## [1] 45.4978

Podemos ver que la root mean squared ha bajado, por lo tanto el hecho de considerar más k, más vecinos, aumenta por el momento la rmse. Seguimos probando más k:

reg3 <- knn.reg(tr[,7:12], val[,7:12], tr$expense, k = 3, algorithm = "brute")
rmse3 <- rmse(val$expense, reg3$pred)
rmse3

## [1] 42.46672

Por último, un vecino más:

reg4 <- knn.reg(tr[,7:12], val[,7:12], tr$expense, k = 4, algorithm = "brute")
rmse4 <- rmse(val$expense, reg4$pred)
rmse4

## [1] 51.62521

Hacemos una vuisualización de los errores:

errors <- c(rmse1, rmse2, rmse3, rmse4)
plot(errors, type = "o", xlab = "k", ylab = "RMSE")

Podemos como conluir que los 3 vecinos son los que nos dan un mejor resultado, con el RMSE más bajo.

Ahora que ya hemos visto el mejor número de vecinos, vamos a comprovar el resultado de k para el conjunto de testing, no el de validación:

reg.test <- knn.reg(tr[,7:12], test[,7:12], tr$expense, k = 3, algorithm = "brute")
rmse.test <- rmse(test$expense, reg.test$pred)
rmse.test

## [1] 45.1279

El error ha subido, pero se mantiene bajo. Trabajar con un conjunto de validación y uno de testing permite hacer un double check a la k encontrada, o también se puede usar técnicas de validación cruzada como el train control del paquete caret, y así evitar que la k se ajuste solo al conjunte de entrenamiento.

K nearest neighbors sin partición de validación

En vez de usar un conjunto de validación para el modelo, vamos a usar la validación cruzada:

t.ids <- createDataPartition(edu$expense, p = 0.7, list = F)
tr <- edu[t.ids,]
val <- edu[-t.ids,]

reg <- knn.reg(tr[,7:12], test = NULL, y = tr$expense, k = 3, algorithm = "brute") 

El modelo que hemos creado ya contiene los residuso, ahora debemos elevarlos al cuadrado:

rmse.reg <- sqrt(mean((reg$residuals)^2))
rmse.reg

## [1] 54.07932

Podemos ver que el error es superior a antes, sin margen de mejora.

Tomamos otra vez el conjunto de entrenamiento, validación y test:

set.seed(2018)
t.ids <- createDataPartition(edu$expense, p = 0.6, list = F)
tr <- edu[t.ids,]
temp <- edu[-t.ids,]
v.ids <- createDataPartition(temp$expense, p = 0.5, list = F)
val <- temp[v.ids,]
test <- temp[-v.ids,]

Vamos a ver como podemos hacer todo el proceso anterior con una función:

rdacb.knn.reg <- function(tr_predictors, val_predictors, 
                          tr_target, val_target, k){
  library(FNN)
  res <- knn.reg(tr_predictors, val_predictors,
                 tr_target, k, algorithm = "brute")
  rmserror <- sqrt(mean((val_target - res$pred)^2))
  cat(paste("RMSE para k = ", toString(k), ": ", rmserror, "\n", sep = ""))
  rmserror
}

Ahora aplicamos la función creada para diferentes k:

• k = 1

rdacb.knn.reg(tr[,7:12], val[,7:12],
              tr$expense, val$expense, 1)

## RMSE para k = 1: 57.5465029345833

## [1] 57.5465

• k = 2

rdacb.knn.reg(tr[,7:12], val[,7:12],
              tr$expense, val$expense, 2)

## RMSE para k = 2: 45.497802144719

## [1] 45.4978

• k = 3

rdacb.knn.reg(tr[,7:12], val[,7:12],
              tr$expense, val$expense, 3)

## RMSE para k = 3: 42.4667189952582

## [1] 42.46672

• k = 4

rdacb.knn.reg(tr[,7:12], val[,7:12],
              tr$expense, val$expense, 4)

## RMSE para k = 4: 51.625211863972

## [1] 51.62521

Podemos ver que el mejor resultado es para k = 3.

Ahora vamos a caluclar una función donde podamos introducir un rank de k en vez de ir comproband cada k:

radcb.knn.reg.multi <- function(tr_predictors, val_predictors,
                                tr_target, val_target, start_k, end_k){
  rms_errors <- vector()
  for(k in start_k:end_k){
    rms_error <- rdacb.knn.reg(tr_predictors, val_predictors, 
                               tr_target, val_target, k)
    rms_errors <- c(rms_errors, rms_error)
  }
  plot(rms_errors, type = "o", xlab = "k", ylab = "RMSE")
}

radcb.knn.reg.multi(tr[,7:12], val[,7:12], tr$expense, val$expense, 1, 10)                          

## RMSE para k = 1: 57.5465029345833
## RMSE para k = 2: 45.497802144719
## RMSE para k = 3: 42.4667189952582
## RMSE para k = 4: 51.625211863972
## RMSE para k = 5: 53.1055552649626
## RMSE para k = 6: 53.1944284885384
## RMSE para k = 7: 52.5103876992594
## RMSE para k = 8: 50.7557262779285
## RMSE para k = 9: 49.4976742189612
## RMSE para k = 10: 50.9256516894973

Aquí podemos ver claramente todas las k y el plot para ver que el vecino 3 es la mejor.

La regresión lineal

Vamos a intentar predecir el consumo de un coche a partir del siguiente dataset:

library(caret)
auto <- read.csv("../DataSets/auto-mpg.csv")
head(auto)

Convertimos a factor:

auto$cylinders <- factor(auto$cylinders, 
                         levels = c(3,4,5,6,8),
                         labels = c("3c","4c","5c","6c","8c"))

Hacemos la partición de los datos:

set.seed(2018)
t.id <- createDataPartition(auto$mpg, p = 0.7, list = F)

Creamos el modelo lineal, excluyendo algunas variables que no nos estan añadiendo ningún valor: (P.D.: el paquete mod ya se encarga de crear varaibles dummies para las categorias que introduzcamos):

mod <- lm(mpg ~., data = auto[t.id, -c(1,8,9)])
mod

## 
## Call:
## lm(formula = mpg ~ ., data = auto[t.id, -c(1, 8, 9)])
## 
## Coefficients:
##  (Intercept)   cylinders4c   cylinders5c   cylinders6c   cylinders8c  
##    38.607312      7.212652      5.610350      3.307172      6.211343  
## displacement    horsepower        weight  acceleration  
##     0.006878     -0.072209     -0.005156      0.024852

Podemos ver que ha encontrado una relación lineal entre mpg y el resto de variables. POdemos ver que no aparecen “cylinders3c”, es decir que no hay ninguno en nuestro dataset de training, porque deben haber muy pocos en el dataset total.

La información que nos está diciendo este modelo es:

mpg = 38.607312 + 7.212652*4c + 5.610350*5c + 3.307172*6c + 6.211343*8c + 0.006878*displacement + horsepower*-0.072209 + weight*-0.00515 + acceleration*0.024852

¿Cómo podemos saber si este modelo predice bien o mal?

summary(mod)

## 
## Call:
## lm(formula = mpg ~ ., data = auto[t.id, -c(1, 8, 9)])
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -10.0708  -2.5374  -0.6035   1.8329  15.5441 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  38.6073120  3.9287480   9.827  < 2e-16 ***
## cylinders4c   7.2126522  2.9190400   2.471   0.0141 *  
## cylinders5c   5.6103498  4.1073757   1.366   0.1731    
## cylinders6c   3.3071717  3.1514235   1.049   0.2949    
## cylinders8c   6.2113432  3.5388015   1.755   0.0804 .  
## displacement  0.0068775  0.0100493   0.684   0.4943    
## horsepower   -0.0722087  0.0177733  -4.063 6.36e-05 ***
## weight       -0.0051560  0.0009051  -5.697 3.18e-08 ***
## acceleration  0.0248515  0.1424361   0.174   0.8616    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.997 on 271 degrees of freedom
## Multiple R-squared:  0.7538, Adjusted R-squared:  0.7466 
## F-statistic: 103.7 on 8 and 271 DF,  p-value: < 2.2e-16

Vemos información sobre el residuo para ver el error que cometemos. Podemos acceder a muchas varaibles del modelo:

boxplot(mod$residuals)

La mayoría de resduos estan en torno al 0, así que parece correcto.

Vemos otras métricas que son relevantes para el modelo:

• p-vales: nos interesan que sean muy bajos, ya que un p value alto implica que ese factor no tiene nada que ver con el modelo, y su valor podría ser 0. Con este valor podríamos elegir que variables o facotores descartar del modelo para así mejorarlo. Las estrellas de la derecha indican los términos que más aportan, relacionadas con el p-value, los que tiene más significatividad.

• residual standard error: podríamos calcularlo nosotros mismos aunque en summary está ajustado a los grados de libertad the chi-squared:

sqrt(mean((mod$fitted.values - auto[t.id,]$mpg)^2))

## [1] 3.932162

• el r-squared y el adjusted r-squared nos dicen que porcentage de la variación está explicada por parte de nuestro modelo.

Si queremos evaluar la predicción con el conunto de validación, hacemos uso de la función predict():

pred <- predict(mod, auto[-t.id, -c(1,8,9)])
sqrt(mean((pred - auto[-t.id,]$mpg)^2))

## [1] 4.17699

Conclusión: nuestro error cuadrático es muy bueno, dado que podemos predecir +/- 4 millas del mpg real.

Ahora vamos a interpretar los parametros estimados respecto el dataset orignal:

par(mfrow = c(2,2))
plot(mod)

Entendiendo los gráficos de los residuos en un modelo lineal

El p vaor nos permite identificar cuando los resultados son significativos y cuando no. Aún así, podemos obtener más información leiendo los gráficos que tenemos encima. Estos nos expicarán si el modelo trabaja bien con la información. Endender los residuos nos permitirá idenfiticar patrones que puede que hayamos tenido en cuenta, afctados por errores de los datos, outliers, …

Residuals vs fitted values

Sirven para detectar varianza no constante en la distribución de los errores. Nos muestra si los de entrada respecto los valores de salida tienen un patrón no lineal.

Si nos econtramos un gráfico como el de la izquierda, podemos ver los residuos más o menos repartidos en torno al eje horizontal, tendiendo a un error 0, sin encontrar un patrón relamente claro. Esto es un buen indicativo ya que los errores no siguen un patrón ni ninguna tendencia lineal.

A la derecha vemos que los residuos sí que siguen un patrón, y que el modelo no ha explicado de forma correcta los errores que tienen. En este caso los errores tienen forma de parábola, por lo tanto el modelo debería buscar una relación cuadrática entre las variables, y no una relación lineal.

Gráficos cuantil-cuantil

Se intenta representar si los errores se distribuyen según una normal estándar o no. Si los puntos se acumulan cerca de larecta y = x, como en el gráfico inferior a la izquierda, tienen una distribución normal.

En el segundo gráfico, a la derecha, valores parecen alejarse más de la recta central, por lo tanto no siguen una distribución normal.

Gráficos Scale-Location

La escala y localización de los residuos, este gráfico nos intenta explicar los rangos que toman los datos. Son los errores estandarizados y con la raíz cuadrada. Con este gráfico se quiere comprovar si los datos tienen todos la misma variancia.

Si tenemos todos los datos que no siguen un patrón, sino que están centrados y no hay ninguna tendencia existente, como en el caso 1, entonces quiere decir que el modelo es bueno. Si exite un patrón, con resiudos con diferentes varianzas, puede que el modelo no esté explicando una parte de los datos.

Residuos vs apalancamiento

Este gráfico nos permite encontrar posibles sujetos influyentes dentro de nuestro modelo. Podrían existir datos como outliers que fueran muy influentes a la hora de crear la recta de regresión lineal.

En el caso 1 no hay ningun error que influencie o cause distorsión en el modelo lineal.

Los casos que estan fuera de la distancia de cook, como en el caso 2, son grandes influenciadores del modelo, y deberían ser analizados con más detalle.

Entendiendo nuestro modelo lineal

par(mfrow = c(2,2))
plot(mod)

El modelo esplic bastante bien todos los residuos:

• Residual vs fitted: no parece que haya ningún patrón que el modelo no explique, con algunos outliers

• Normal Q-Q, se ajusta mucho a la recta x = y, por lo tanto los datos estan normalizados. Aparecen los mismos outliers en que en residuals vs fitted.

• Scale-Location: El gráfico es el mismo que el residual vs fitted pero reducido con la raíz cuadrada con los residuos estandarizados. No hay mucha tendencia por lo tanto no podíamos decir que existe un patrón.

• Residual vs laverage (gráfico de apalancamiento): hay objetos muy lejos pero no sobrepasan la distancia de cook

Vamos a dare un vistazo a les errores más destacados:

• El el coche nº 300 pesa muchos y es un 6 cilindros, de los que no había mucho en la muestra. Puede que nos hiciese falta más datos de cochos con 6 cilindros para obtener un mejor modelo.

• EL coche nº 299 es un 4 cilindros que consume muy poco, siendo un caso atípico.

Opciones para las fórmulas de un modelo lineal

Se puede forzar que e linear model use un factor específico. En el anterior caso ha cogido como referencia el primer de los factores; hay 5 factores, 3 cilindros, 4 cilindros…, y ha cogido el de 3 porque era el que estaba antes. Si queremos ue coja el d 4 cilindros dado que es el que tiene más número de coches, lo haríamos de la siguiente forma:

auto <- within(auto, 
               cylinders <- relevel(cylinders, ref = "4c")) # hemos puesto que el 4c se rejustará sobre 4c

Vamos a ver si el modelo cambia:

mod <- lm(mpg ~., data = auto[t.id, -c(1,8,9)])
mod

## 
## Call:
## lm(formula = mpg ~ ., data = auto[t.id, -c(1, 8, 9)])
## 
## Coefficients:
##  (Intercept)   cylinders3c   cylinders5c   cylinders6c   cylinders8c  
##    45.819964     -7.212652     -1.602302     -3.905480     -1.001309  
## displacement    horsepower        weight  acceleration  
##     0.006878     -0.072209     -0.005156      0.024852

Ahora sí que aparece el 3cilindros, antes se tomaba como valor basico. Ahora se considera que lo estándard es el modelo de 4cilindros, un coche por defecto es de 4 cilindros. Recomendación: elegir aquella categoría que tenga más elementos, así a la hora de hacer computos exactos se pueden llevar a cabo menos cálculos.

Usamos este modelo para hacer la nueva predicción:

pred <- predict(mod, auto[-t.id, -c(1,8,9)])
sqrt(mean((pred - auto[-t.id,]$mpg)^2))

## [1] 4.17699

plot(mod)

postResample(pred, auto[-t.id,]$mpg)

##      RMSE  Rsquared       MAE 
## 4.1769896 0.6971008 3.1935939

¿Qué ocurre cuando nuestro modelo no es lineal?

POr ejemplo: cuando observamos cierta tendencia en los residuos que nos indica que los valores siguan un patrón cuadrático. La tabla a continuación podría defeinier la mayoría de funciones LM:

Función step para simplificar un modelo lineal

¿Debemos introducir todas las variables en el modelo, o puede que alguna sobre?

library(MASS)

Hacemos la selección de las variables relevantes para el modelo:

step.model <- stepAIC(mod, direction = "backward") # step Acaique Information Criterium, si ponemos forwards, irá añadiendo variables, si hacemos backward irá quitando variables a partir de un modeo completo

## Start:  AIC=784.75
## mpg ~ cylinders + displacement + horsepower + weight + acceleration
## 
##                Df Sum of Sq    RSS    AIC
## - acceleration  1      0.49 4329.8 782.78
## - displacement  1      7.48 4336.8 783.23
## <none>                      4329.3 784.75
## - horsepower    1    263.69 4593.0 799.30
## - cylinders     4    600.26 4929.6 813.10
## - weight        1    518.45 4847.8 814.42
## 
## Step:  AIC=782.78
## mpg ~ cylinders + displacement + horsepower + weight
## 
##                Df Sum of Sq    RSS    AIC
## - displacement  1      7.03 4336.8 781.23
## <none>                      4329.8 782.78
## - horsepower    1    402.85 4732.7 805.69
## - cylinders     4    600.05 4929.9 811.12
## - weight        1    633.52 4963.3 819.01
## 
## Step:  AIC=781.23
## mpg ~ cylinders + horsepower + weight
## 
##              Df Sum of Sq    RSS    AIC
## <none>                    4336.8 781.23
## - horsepower  1    420.53 4757.4 805.15
## - cylinders   4    593.99 4930.8 809.17
## - weight      1    713.02 5049.9 821.85

Podemos ver como la función step ha ido eliminando paso a paso toda la información, hasta quedarse solo con horsepower, cylinders y weight. Este proceso lo ha seguido al añadir el parámetro “backward”, donde ha ido evaluando cada variables y quitando las no relevantes. Decide las variables no relevantes cuando estan por debajo del AIC global (none).

Vamos a ver la inforrmación con mas detalle:

summary(step.model)

## 
## Call:
## lm(formula = mpg ~ cylinders + horsepower + weight, data = auto[t.id, 
##     -c(1, 8, 9)])
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -9.9463 -2.5897 -0.6286  1.8488 15.9613 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 46.134348   1.540252  29.952  < 2e-16 ***
## cylinders3c -7.597071   2.850293  -2.665  0.00815 ** 
## cylinders5c -1.524612   2.897902  -0.526  0.59924    
## cylinders6c -3.488165   0.846567  -4.120 5.02e-05 ***
## cylinders8c -0.180428   1.373267  -0.131  0.89557    
## horsepower  -0.070401   0.013683  -5.145 5.12e-07 ***
## weight      -0.004850   0.000724  -6.700 1.19e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.986 on 273 degrees of freedom
## Multiple R-squared:  0.7534, Adjusted R-squared:  0.748 
## F-statistic:   139 on 6 and 273 DF,  p-value: < 2.2e-16

Si lo hacemos al reves (el AIC):

step.model2 <- stepAIC(mod, direction = "forward")

## Start:  AIC=784.75
## mpg ~ cylinders + displacement + horsepower + weight + acceleration

Los árboles de regresión

Modelos en forma de árbol para llevar acabo estimación de valores.

library(rpart)
library(rpart.plot)
library(caret)
bh <- read.csv("../DataSets/BostonHousing.csv")

Vamos a predecir la variables de habitabilidad de la vivienda.

Creamos una partición para hacer el modelo:

set.seed(2018)
t.id <- createDataPartition(bh$MEDV, p = 0.7, list = F)
bfit <- rpart(MEDV ~., data = bh[t.id,])
bfit

## n= 356 
## 
## node), split, n, deviance, yval
##       * denotes terminal node
## 
##  1) root 356 30671.9100 22.69101  
##    2) RM< 6.9715 301 11935.0000 19.93090  
##      4) LSTAT>=14.405 122  2049.6700 14.97787  
##        8) CRIM>=6.99237 53   608.5698 12.10189 *
##        9) CRIM< 6.99237 69   665.9983 17.18696 *
##      5) LSTAT< 14.405 179  4852.4720 23.30670  
##       10) LSTAT>=5.145 160  2537.4140 22.40625  
##         20) LSTAT>=10.14 73   452.9088 20.60411 *
##         21) LSTAT< 10.14 87  1648.4910 23.91839 *
##       11) LSTAT< 5.145 19  1092.8580 30.88947 *
##    3) RM>=6.9715 55  3894.3790 37.79636  
##      6) RM< 7.435 32   933.8472 32.70938 *
##      7) RM>=7.435 23   980.3443 44.87391 *

Entendiendo la estructura de est árbol:

• root: primer nodo y a partir de aquí las condiciones para dividir el árbol

• n: el número de casos que caen en el nodo (total de casos son 356)

• deviance: la suma del cuadrado de los errores en cada nodo

• yval: valor promedio de las variables de salida de todos los casos del nodo, valor estimado cuando alguien cae en ese nodo

• el * indica el nodo terminal

Visualización del modelo:

prp(bfit, type = 2, nn = T, fallen.leaves = T, faclen = 4,
    varlen = 8, shadow.col = "gray", roundint = F)

El número final nos muestra el valor estimado en los nodos hoja.

bfit$cptable

##           CP nsplit rel error    xerror       xstd
## 1 0.48391276      0 1.0000000 1.0077885 0.09969609
## 2 0.16408695      1 0.5160872 0.5592724 0.06462367
## 3 0.06456030      2 0.3520003 0.3935890 0.05349306
## 4 0.03984754      3 0.2874400 0.3211048 0.05065157
## 5 0.02527075      4 0.2475924 0.3063209 0.05227560
## 6 0.01421543      5 0.2223217 0.2792451 0.04993005
## 7 0.01000000      6 0.2081063 0.2856330 0.05017668

Esta tabla nos muestra valores para cada uno de los nodos:

• CP: número de complejidad del árbol

• nsplit: divisiones que se han hecho en el árbol para llegar a este nivel de complejidad

• rel error: para el mejor árbol con el número especifico de dvisiones, se calcula el error de clasificación total con raíz cuadrada.

• xerro: error de validación cruzada promedio, suma de los errores dividido por el número de árboles llevado a cabo.

• xstd: desviación estándard con la validación cruzada, aplicada sobre el mejor árbol que ha salido con el número específico de los splits (divisiones).

Para elegir el mejor número de splits, podemos elegir el árbol que tenga menor número de corss validation (xerror), o utilizar la regla de una desviación típica y elegir el árbol que se encuentra dentro de una desviación típica del error mínimo que tenga menos nodos.

Este último enfoque nos permite elegir árboles con menos divisiones que el original, y que sea más sencillo que el orginal.

0.2856330 + 0.05017668 # sumamos el error de validación cruzada  más pequeño con su desviación standard para establecer el error mínimo

## [1] 0.3358097

El primer árbol que se encuentra por debajo de este margen que hemos establecido, sería el de CP igual a 4, y podríamos tomarlo como bueno.

plotcp(bfit)

Con esta representación podemos ver que árbol tomar, y así escoger el modelo más óptimo. Podemos ver que el 4 es el mejor. ¿Cómo lo podríamos calcular el árbol “pruned”?

  bfitpruned <- prune(bfit, cp = 0.03984754) # ponemos el valor de complejidad para el árbol 4

prp(bfitpruned, type = 2, nn = T, fallen.leaves = T, faclen = 4,
    varlen = 8, shadow.col = "gray", roundint = F)

Aquí podemos ver que hemos hecho una poda de árbol simplificando el árbol anterior. Ahora que tenemos el modelo, lo aplicamos:

• Validación en el conjunto de entrenamiento con poda

preds <- predict(bfitpruned, bh[t.id,])
sqrt(mean((preds - bh[t.id,]$MEDV)^2))

## [1] 4.61864

postResample(preds, bh[t.id,]$MEDV)

##      RMSE  Rsquared       MAE 
## 4.6186400 0.7524076 3.2394648

• Validación en el test con poda

preds <- predict(bfitpruned, bh[-t.id,])
sqrt(mean((preds - bh[-t.id,]$MEDV)^2))

## [1] 4.915902

postResample(preds, bh[-t.id,]$MEDV)

##      RMSE  Rsquared       MAE 
## 4.9159021 0.7000611 3.4991285

Ahora comparamos los resultados sin la “poda” del árbol

• Validación en el conjunto de entrenamiento

preds <- predict(bfit, bh[t.id,])
sqrt(mean((preds - bh[t.id,]$MEDV)^2))

## [1] 4.234362

postResample(preds, bh[t.id,]$MEDV)

##      RMSE  Rsquared       MAE 
## 4.2343621 0.7918937 2.8517371

• Validación en el test

preds <- predict(bfit, bh[-t.id,])
sqrt(mean((preds - bh[-t.id,]$MEDV)^2))

## [1] 4.469432

postResample(preds, bh[-t.id,]$MEDV)

##      RMSE  Rsquared       MAE 
## 4.4694321 0.7520928 3.0418549

Los reslutados son mejores, pero con el árbol podado optimizaremos mucho más el sistema, tampoco son mucho peores y evitaremos overfitting.

Las técnicas de Bagging and Boosting

¿Qué ocurriría si alguna de las variables predictoras fuera categórica?

Queremos calcular el gasto en educación para el siguiente dataset:

ed <- read.csv("../DataSets/education.csv")
ed$region <- factor(ed$region)

t.id <- createDataPartition(ed$expense, p = 0.7, list = F)

fit <- rpart(expense ~ region + urban + income + under18, data = ed[t.id,])

prp(fit, type = 2, nn = T, fallen.leaves = T, faclen = 4, varlen = 8, 
    shadow.col = "grey", roundint = F)

Se pueden generar árboles de regresión usando el método de ensamblaje:

• bagging

• boosting

Bagging

El bagging o del bootstrap agggregation, es un método de ensamblaje que combina de forma conjunta las predicciones de diferentes árboles por separado, para dar un resultado mejor.

library(ipred)

Usaremos el mismo dataset que antes con la misma partición sobre el conjunto de entrenamiento:

bagging.fit <- bagging(MEDV ~., data = bh[t.id, ])
prediction.t <- predict(bagging.fit, bh[t.id,])
sqrt(mean((prediction.t - bh[t.id,]$MEDV)^2))

## [1] 2.64692

caret::postResample(prediction.t, bh[t.id,]$MEDV)

##      RMSE  Rsquared       MAE 
## 2.6469200 0.8399787 1.8710443

I sobre el conjunto de validación:

prediction.t <- predict(bagging.fit, bh[-t.id,])
sqrt(mean((prediction.t - bh[-t.id,]$MEDV)^2))

## [1] 7.386543

caret::postResample(prediction.t, bh[-t.id,]$MEDV)

##      RMSE  Rsquared       MAE 
## 7.3865434 0.4473396 5.0234741

Boosting

Este método nos permite reducir el desvío cuando los nuevos modelos se construyen para aprender de una mala clasificación de modelos prévios. Es decir, cada vez que se hacen nuevos modelos, se aprenden de sus errores, y el nuevo modelo “boosted” lo hace mejor.

library(gbm)

## Loaded gbm 2.1.4

gbm.fit <- gbm( MEDV ~., data = bh, distribution = "gaussian") # la muestra es demasiado pequeña y nos vemos obligados  a coger el conjunto del dataset

Bosques aleatorios para regresión

library(randomForest)

## randomForest 4.6-14

## Type rfNews() to see new features/changes/bug fixes.

## 
## Attaching package: 'randomForest'

## The following object is masked from 'package:ggplot2':
## 
##     margin

library(caret)

bh <- read.csv("../DataSets/BostonHousing.csv")

set.seed(2018)
t.id <- createDataPartition(bh$MEDV, p = 0.7, list = F)

mod <- randomForest(x = bh[t.id, 1:13], y = bh[t.id, 14],
                    ntree = 1000, 
                    xtest = bh[-t.id, 1:13], ytest = bh[-t.id,14],
                    importance = T, keep.forest = T) # importance nos dice si tiene que computar o no las puntuaciones de cada variable  predictora
mod

## 
## Call:
##  randomForest(x = bh[t.id, 1:13], y = bh[t.id, 14], xtest = bh[-t.id,      1:13], ytest = bh[-t.id, 14], ntree = 1000, importance = T,      keep.forest = T) 
##                Type of random forest: regression
##                      Number of trees: 1000
## No. of variables tried at each split: 4
## 
##           Mean of squared residuals: 11.55283
##                     % Var explained: 86.59
##                        Test set MSE: 10.61
##                     % Var explained: 86.81

Datos que podemos ver:

• 1000 árboles creados

• Número de variables elegidos en cada división han llegado a ser hasta 4

• El primer % Var explained es sobre el conjunto de entrenamiento, y el otro % Var explained es sobre el conjunto de testing.

mod$importance

##           %IncMSE IncNodePurity
## CRIM     8.377029     1999.1072
## ZN       1.021550      324.6305
## INDUS    7.666375     1950.1366
## CHAS     1.017831      201.5456
## NOX      7.240169     1660.1109
## RM      36.749369     9121.9544
## AGE      3.437088      877.4536
## DIS      6.333173     1660.5814
## RAD      1.213486      249.2855
## TAX      3.142002      763.1997
## PTRATIO  7.935050     2043.6182
## B        1.547519      622.8623
## LSTAT   55.177598     8606.9462

La componente de importancia del modelo, nos dice de cada uno de sus predictores la importancia.

Vamos a ver las predicciones sobre el conjunto de entrenamiento:

plot(bh[t.id,]$MEDV, predict(mod, newdata = bh[t.id,]),
     xlab = "Actual", ylab = "Predichos")
abline(0,1)

Las predicciones sobre el conjunto de entrenamiento parecen ser muy buenas. Vemos ahora sobre el conjunto de testing:

plot(bh[-t.id,]$MEDV, predict(mod, newdata = bh[-t.id,]),
     xlab = "Actual", ylab = "Predichos")
abline(0,1)

Vemos que también se ajusta bastante bien. Podemos ver que hay algunos problemas a la derecha del gráfico.

Podemos usar parámetros para mejorar nuestro randomForest:

• mtry: cuantos predictores vamos a tomar en cada una de las divisiones del bosque aleatorio. Por defecto, este parámetro está configurado como “m/3”, m = nº predictores seleccionados aleatoriamente.

• nodsize: mínimo tamaño que necesitan tener los nodos terminales para ser considerados. Por defecto siempre vale 5, como mínimo necesitamos que 5 elementos caigan dentro del nodo para ser considerado como terminal.

• maxnodes: número máximo de nodes terminales. Por defecto buscan el tamáño máximo posible dentro del nodsize.

Redes neuronales para regresión

library(nnet)
library(caret)
library(devtools)

bh <- read.csv("../DataSets/BostonHousing.csv")
set.seed(2018)
t.id <- createDataPartition(bh$MEDV, p= 0.7, list = F)

Vamos a encontrar el rango de variabilidad de la respuesta para poder escalarlo a 0/1.

summary(bh$MEDV) # podemos ver que el máximo es 50, pero a la hora de construir la red neuronal voy a decirle ue me prediga MEDV pero dividido entre 50 para que los valores esten escalados entre 0 i 1.

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    5.00   17.02   21.20   22.53   25.00   50.00

Ahora escalaremos los resultados de MEDV cuando creemos el modelo:

fit <- nnet(MEDV/50 ~., data=bh[t.id, ], # escalamos a 50 el MEDV
            size = 6, # número de máximo de nodos en la capa oculta, conectando los  datos de entrada por los datos de salida
            decay = 0.1, # decaimiento que tiene que sufrir cada nodo para entrar en uno y salir en otro
            maxit = 1000, # nº máx iteraciones de la red neuronal
            linout=T) # nos  especifica si queremos una salida lineal en función de los parámetros de entrada, y no logística

## # weights:  91
## initial  value 35.628223 
## iter  10 value 10.965976
## iter  20 value 8.483621
## iter  30 value 6.294559
## iter  40 value 4.852913
## iter  50 value 2.799252
## iter  60 value 2.588793
## iter  70 value 2.201279
## iter  80 value 1.823838
## iter  90 value 1.657157
## iter 100 value 1.615447
## iter 110 value 1.525545
## iter 120 value 1.510856
## iter 130 value 1.477779
## iter 140 value 1.408658
## iter 150 value 1.383164
## iter 160 value 1.332337
## iter 170 value 1.288604
## iter 180 value 1.256027
## iter 190 value 1.243384
## iter 200 value 1.242302
## iter 210 value 1.242262
## iter 220 value 1.242236
## iter 230 value 1.242077
## iter 240 value 1.241900
## iter 250 value 1.241742
## iter 260 value 1.239741
## iter 270 value 1.237501
## iter 280 value 1.236978
## final  value 1.236975 
## converged

Vamos a visualizar la red neuronal usando un repositior de git:

source_url("https://gist.githubusercontent.com/fawda123/7471137/raw/466c1474d0a505ff044412703516c34f1a4684a5/nnet_plot_update.r")

## SHA-1 hash of file is 74c80bd5ddbc17ab3ae5ece9c0ed9beb612e87ef

plot(fit, max.sp = T) # max.sp nos intenta dar el max space para  visualizar los nodos

## Loading required package: reshape

Podemos ver la 13 variabes de entrada, en la capa ocuta máximo los 6 nodos.

El resultado que nos da esta red neuronal esta escalado entre 0 i 1, por lo tanto deberemos transformar a datos lógicos:

# error cuadrático medio sobre el valor de la casa
sqrt(mean((fit$fitted.values*50-bh[t.id,"MEDV"])^2)) # nos da un error muy bajo

## [1] 2.453378

# comprobamos los resultados sobre el conjunto de validación
pred <- predict(fit, bh[-t.id,])
sqrt(mean((pred*50 -  bh[-t.id,"MEDV"])^2))

## [1] 3.479434

Implementando una k-fold cross validation en R

Algunas de las técnicas que hemos visto para classificación o árboles de regresión, utilizan directamente la validación cruzada automaticamente para ayudar a la selección del modelo y evitar el overfitting. Por otro lado, hay algunas técnicas que no lo hacen de serie.

Vamos a ver como programar el metodo de k-fold cross validation para la regresión lineal:

#k-fold cross validation
bh <- read.csv("../DataSets/BostonHousing.csv")

# Vamos a crear dos funciones:

# 1.k-fold cross validation
kfold.crossval.reg <- function(df, nfolds){
  fold <- sample(1:nfolds, nrow(df), replace = T)
  cat(fold)
  
  mean.sqr.errs <- sapply(1:nfolds, 
                          kfold.cval.reg.iter,
                          df, fold)
  
  list("MSE "= mean.sqr.errs,
       "Overall_Mean_Sqr_Error" = mean(mean.sqr.errs),
       "Std_Mean_Sqr_Error" = sd(mean.sqr.errs))
}

kfold.cval.reg.iter <- function(k, df, fold){
  
  tr.ids <- !fold %in% c(k)
  test.ids <- fold %in% c(k)
  mod <- lm(MEDV ~., data = df[tr.ids,])
  pred <- predict(mod, df[test.ids,])
  sqr.errs <- (pred - df[test.ids,"MEDV"])^2
  mean(sqr.errs)
}

res <- kfold.crossval.reg(bh, 5)

## 5 5 1 5 3 3 3 3 5 4 3 3 5 5 2 1 5 2 1 3 5 5 3 3 1 4 4 4 3 5 4 1 4 1 2 1 2 1 5 5 1 3 4 1 5 5 4 2 5 5 4 1 5 5 5 5 1 5 1 1 5 4 2 5 2 3 3 4 1 1 2 4 3 1 5 3 2 4 1 5 2 1 4 1 5 4 3 5 5 1 3 1 4 5 2 1 4 1 3 3 4 4 2 1 5 3 2 5 1 4 1 1 1 1 5 1 1 4 2 5 2 2 4 3 4 5 3 5 2 1 4 4 2 2 4 3 3 1 3 1 2 1 5 1 5 4 1 3 4 2 5 4 5 1 5 3 1 4 3 3 4 2 1 5 4 3 2 3 2 3 2 2 4 5 2 5 3 4 3 2 5 2 5 2 5 3 1 2 5 5 5 2 4 5 4 1 4 4 5 5 2 2 2 3 1 5 1 3 4 3 5 1 2 2 1 3 2 5 1 3 4 1 5 5 3 5 1 2 3 5 2 4 1 1 4 3 4 1 2 4 3 3 5 5 3 4 5 3 5 1 5 1 3 3 5 3 4 4 1 5 1 3 1 2 4 2 2 5 4 4 3 2 1 1 2 2 1 1 3 3 3 3 4 1 1 5 5 2 4 2 2 3 4 2 5 3 4 3 5 3 1 3 2 1 4 5 2 2 2 3 2 1 2 3 2 2 2 2 1 3 4 4 1 2 2 5 1 5 3 5 4 5 1 5 1 3 4 1 2 1 5 2 5 4 4 3 5 1 1 3 4 3 1 3 1 2 1 3 3 2 4 4 1 1 3 1 2 4 4 1 1 1 4 1 4 5 5 4 4 1 2 5 1 2 3 3 2 2 5 3 2 4 1 1 4 1 3 3 5 1 1 1 1 4 4 5 2 4 2 2 1 1 4 5 4 1 2 3 3 5 2 5 3 4 4 5 5 1 2 1 3 1 4 4 2 3 1 5 3 4 4 3 3 4 2 1 3 4 3 1 1 1 5 1 5 1 4 4 4 3 1 1 5 3 5 2 5 5 5 3 4 5 3 5 1 3 3 1 4 5 1 1 5 4 5 3 2 3 3 3 4 3 3 2 3 2 4 3 4 2 5 1 3 5 4 1

loocv.reg <- function(df){
  mean.sqr.errors <- sapply(1:nrow(df), 
                            loocv.reg.iter, df)
  list("MSE"=mean.sqr.errors,
       "overall_mean_sqr_errors" = mean(mean.sqr.errors),
       "sd_mean_sqr_errors" = sd(mean.sqr.errors))
}

loocv.reg.iter <- function(k, df){
  mod <- lm(MEDV~., data = df[-k,])
  pred <- predict(mod, df[k,])
  sqr.error <- (pred - df[k,"MEDV"])^2
  sqr.error
}

res <- loocv.reg(bh)