Un ejemplo simple de "ensamble learning" con modelos de regresión (Post I de varios...)

Ensamblando modelos

La idea de este post es introducir algunas nociones de ensamble de modelos y presentar un pequeño código para implementar uno para un modelo de regresión a título de ilustración. El post pretende ordenar (básicamente para mí mismo) algunas ideas al respecto e ilustrarlas con un código y un modelo muy simples.

En líneas generales los ensambles de modelos buscan aumentar la capacidad predictiva a partir de la generación de muchos "modelitos" diferentes sobre datasets diferentes. Obviamente, lo ideal sería tener muchas muestras independientes de la misma población para poder estimar modelos diferentes y promediarlos de alguna manera (de esta manera, se reduciría/suavizaría la variancia y se mantendría el sesgo relativamente constante). Como consecuencia, deberían mejorar las predicciones.

Ahora bien, como esto no suele ser posible (tenemos UNA sola muestra de la población) lo que puede utilizarse con diversas técnicas de remuestreo (por ejemplo, bootstrap) de una misma muestra. Sobre esas remuestras se estiman diversos modelos y se los agrega de diversas maneras. Hay muchos algoritmos de "ensamblado" (no sé si es la mejor traducción): los más conocidos son el bagging y el boosting.

El algoritmo bagging es una forma simple de incrementar el poder predictivo de un modelo. EL paper al respecto lo escribió el CAPO de Leo Breiman. Si bien es ampliamente utilizado en otras ramas de la ciencia... en sociales (¿cuándo no?) estamos atrasados en la aplicación de técnicas más sofisticadas en el análisis de datos. La intención de este post es presentarlo.

La idea general es generar un conjunto de modelos diferentes para agregarlos de alguna manera. Ahora, ¿cómo hacemos para generar esos diferentes modelos si, en general, tenemos un solo dataset? Básicamente, haciendo bootstrap. O sea, boostrapeamos, estimamos los modelos y agregramos los resultados (Bootstrap AGGregatING). Lo que se espera obtener es una reducción del sesgo de los modelos y su variancia "promediando" a través de muchos modelos posibles. En líneas generales, bagging funciona mejor con clasificadores "inestables" (por ejemplo, árboles de decisión, best-subset en regresiones lineales, regresiones con regularización, redes neuronales, etc.). Acá lo vamos a presentar con un modelo más bien simple, para no entrar en detalles, una regresión lineal.

Para ilustrar la cosa, vamos usar datos reales. El modelo va a ser bastante poco imaginativo: un modelo de regresión lineal con el ingreso total individual como variable dependiente; el nivel educativo, la cantidad de horas trabajadas, la calificación y la edad como predictoras (algo así como una regresión de Mincer ampliada). Usamos el dataset correspondiente al segundo trimestre del 2015 de la Encuesta Permanente de Hogares. Nos quedamos con las variables p47t (ingreso individual total), nivel_ed (nivel educativo), htot (horas trabajadas), calif (calfiicación) y ch06 (edad en años simples). Al ingreso lo tomamos en logaritmo 10.

Estimando el modelo lineal simple

En primera instancia estadarizamos las dos variables cuantitativas (cantidad de horas trabajadas por semana -htot-) y edad (ch06). Luego, generamos un índice para dividir nuestra base de datos en un test set y un training set (el modelo lo entrenamos, obviamente, en el training set).

library(foreach)
library(parallel)
library(doSNOW)
## Loading required package: iterators
## Loading required package: snow
## 
## Attaching package: 'snow'
## The following objects are masked from 'package:parallel':
## 
##     clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
##     clusterExport, clusterMap, clusterSplit, makeCluster,
##     parApply, parCapply, parLapply, parRapply, parSapply,
##     splitIndices, stopCluster
setwd("E:/PEN1/Blog datos y hechos/2016_Ensambles")
data_filt<-read.csv("data_filt.csv")

data_filt$ch06<-(data_filt$ch06-mean(data_filt$ch06))/sd(data_filt$ch06)
data_filt$htot<-(data_filt$htot-mean(data_filt$htot))/sd(data_filt$htot)

set.seed(1598)
tr_ind<-sample(1:nrow(data_filt),round(nrow(data_filt)*0.75,0),replace=FALSE)
train<-data_filt[tr_ind,]
test<-data_filt[-tr_ind,]

lmod<-lm(log10(p47t)~.,data=train)
lmod
## 
## Call:
## lm(formula = log10(p47t) ~ ., data = train)
## 
## Coefficients:
##               (Intercept)                       ch06  
##                   3.81032                    0.06335  
## nivel_ed2_Sec. comp y más                       htot  
##                   0.14487                    0.12183  
##      calif2_Op./No calif.  
##                  -0.14358
pr_simple<-predict(lmod,newdata=test)

Luego, estimamos (indexando con tr_ind) sobre el training set un modelo lineal y hacemos la predicción sobre el training set y calculamos el error cuadratico medio.

Vemos que el error es igual a

mean((pr_simple-log10(test$p47t))**2)
## [1] 0.08271727
Baggeando

Bien, ahora, tenemos que generar el ensamble a través del algoritmo "bagging". La idea es generar muchas muestras (en este caso, para simplificar, con muestreo aleatorio simple con reposición) y estimar para cada una de ellas un modelo de regresión lineal. Hay un montón de paquetes en R para hacer bagging pero está bueno programarlo desde cero para entender la lógica que lo sustenta.

Los pasos serian los siguientes

- generar una muestra boostrap del training set

- sobre ella estimamos una regresión

- con ese modelo predecimos la variable dependiente en el test set

Repetimos esto "rep" veces. (Está seteado 500 por defecto)

Para esto, tenemos que hacer un loop. Ahora, como los for loop en R son un tanto lentos, podemos usar el paquete foreach para paralelizar en cada uno de los núcleos disponibles y achicar los tiempos de procesamiento.

rep<-500
cl<-makeCluster(detectCores())
registerDoSNOW(cl)

ens<-foreach(m=1:rep,.combine="cbind") %dopar% {
        index<-sample(1:nrow(train),nrow(train),replace=TRUE)
        lmod<-lm(log10(p47t)~.,data=train[index,])
        predict(lmod,newdata=test)
}

pr_ens<-rowMeans(ens)

"rep" setea el total de repeticiones. "cl" y "registerDoSNOW" genera la cantidad de núcleos para paralelizar el cómputo (cuántos más núcleos, mejor).

El resto es un for loop pero usando el paquete "foreach", que permite paralelizarlo. Más allá de cuestiones referidas a la paralelización, puede verse que en cada iteración entre 1 y "rep" lo que hace la función es, básicamente, lo que habíamos planteado. La primera línea dentro del loop genera un índice para muestrear con reposición registros dentro del training set. La segunda línea estima el modelo lineal y la tercera realiza la predicción sobre el test set.

La predicción final es un agregado de todas las predicciones a lo largo de las "rep" muestras. En nuestro caso, usamos la media. Pero podría haber sido otra: mediana, "majority vote", medias ponderadas, etc.

En la corrida r=500, tenemos que el error es igual a

mean((pr_ens-log10(test$p47t))^2)
## [1] 0.08271705

Ahora, para distintos valores de "rep" tenemos que

- para rep=100 el error=0.08271209

- para rep=2000 el error=0.08271708

- para rep=5000 el error=0.08271716

- para rep=15000 el error=0.08271727

- para rep=20000 el error=0.08271711

Vemos, entonces, que aparece una ligera tendencia al descenso. O sea, "marginalmente" la cosa funciona. Un ejemplo similar (con variables simuladas) lo pueden ver acá.

A su vez, para dar una idea de los tiempos: paralelizando en cuatro núcleos (i5) correr 20000 modelos tarda aproximadamente 8 minutos.

Obviamente, este es un ejemplo general pero podría mejorarse notablemente. Por ejemplo, mejorando el muestreo tanto en la generación de la partición test-training como en la generación de las sucesivas muestras bootstrap. Podríamos usar, por ejemplo, las funciones que habíamos armado en los posts anteriores para muestrear de forma estratificada (un poco más parecida a lo que fue el diseño original de la encuesta): por ejemplo, por aglomerado.

A la vez, acá estamos usando un modelo lineal. Este suele ser bastante "estable" es por ello que la eficiencia del bagging parece marginal Podríamos usar algún otro (en eso estamos...). Lo dejamos para otros posts (y papers...).

Al mismo tiempo, estamos trabajando sobre una sola partición en test y training del dataset. Una opción (para más adelante) sería repetir muchas veces esa partición y agregar los resultados del error cuadrático medio (esto es lo que hace Breiman en su paper).

Y también dejamos el problema planteado de cómo lograr cuantificar la importancia o la magnitud del efecto de las diversas variables dependientes.

Finalmente: ¿qué pasaría si armáramos un ensamble de varios modelos diferentes? ¿Un modelo lineal y un random forest, por ejemplo?

PS: No sé si notaron el cambio en el diseño... Todo gracias a la magia de #RMarkdown.

Comentarios

Entradas populares de este blog

Sobre la multicausalidad

Número efectivo de partidos en elecciones presidenciales 1983-2011

El ruido de las capitales (vol. 1)