Regresión no lineal, Cross-Validation y Regularization.

En esta entrada sigo la postura del texto Machine Learning for Hackers[1], la cual es explicar los tres temas juntos para mostrar como se conectan. Se pueden revisar por separado de manera amplia y desde una perspectiva estadística en la referencia [2].

En breve, se puede mencionar que Cross-Validation y Regularization funcionan bien usándolos para determinar un buen modelo o determinar cuando se está “sobre estimando” o “sobre ajustado” (overfitting).

La regresión no-lineal, como su nombre lo indica; es además de una generalización de la regresión lineal, también otra opción para construir un modelo. Al ser más general, no es extraño que se pueda caer en la idea de pensar que cuanto más complejo el modelo es mejor, lo cual no es necesariamente cierto.

Algo que suele ocurrir es que uno busca por medio de un modelo lineal estudiar un fenómenos no lineales, o en otro caso implementar alguna técnica que funciona bien con fenómenos no lineales en fenómenos que pueden ser “bien” estudiados con modelos lineales. Por lo cual es importante explorar la información y tratar de construir más de un modelo, lineal o no-lineal y hacer uso de herramientas para determinar cual es un mejor modelo.

Sobre la idea de complejidad de un modelo, se pude hacer un análisis bastante serio al respecto, pero uno se acercaría al concepto de complejidad  y este puede tener varios puntos de vista los cuales dejo sin detalle en esta entrada, basta mencionar que la complejidad tienen varias aristas y puntos de vista, que van desde temas de matemática puras hasta temas de sistemas complejos y teoría de computación.

Hago un  primer  un ejemplo de ajuste lineal y no-lineal.

#Ejemplo no lineal
x<-seq(-10,10,by=0.01)
y<-1-x^2+rnorm(length(x),0,7)
ggplot(data.frame(X=x,Y=y),aes(x=X,y=Y))+geom_point(col=2)+geom_smooth(size=3,method='lm',se=FALSE)
ggplot(data.frame(X=x,Y=y),aes(x=X,y=Y))+geom_point(col=2)+geom_smooth(size=3,se=FALSE)

Ajuste_LinealAjuste_No-Lineal

La gráfica muestra que al calcular de manera lineal la regresión, la recta estimada corta los datos y al calcular la regresión no-lineal se estima una curva, la cual parece un mejor ajuste. Pero si pensamos en transformar los datos de entrada, elevándolos al cuadrado, se obtiene otro tipo de comportamiento.

#Ejemplo de ajusto con datos transformados
x.cuadrados<-x^2
ggplot(data.frame(X.Cuadrados=x.cuadrados,Y=y),aes(x=X.Cuadrados,y=Y))+geom_point(col=2)+geom_smooth(size=3,se=FALSE)+ ggtitle("Ajuste de datos transformados")+theme(plot.title = element_text(lineheight=.8, face="bold"))

 

Ajuste_No_lineal_con_datos_transformados

La gráfica anterior muestra casi una relación lineal con el ajuste no-lineal, se puede revisar los residuos o la r cuadrada para estimar cuanto de la variabilidad de Y es explicada con los ajustes. Lo hago en R  por medio del comando summary y el resultado es el siguiente:

#Revisamos las estimaciones
> summary(lm(y~x))$r.squared
[1] 0.0001486085
> summary(lm(y~x.cuadrada))$r.squared
[1] 0.9444159

En el ejemplo se observa como una transformación en los datos hace que se pase de un ajuste que explica prácticamente el 0%, a un ajuste que explica casi el 95% de la variabilidad de los datos de salida.

Del ejemplo se puede concluir, que una transformación  hace una gran diferencia y en este caso beneficia al modelo, pero lo que no se visualiza es el costo de transformar los datos, nada limita a pensar que se  puede considerar cualquier potencia de la información y revisar el ajuste. Con potencia me refiero a elevar a la tercera, cuarta o quinta potencia nuestros datos de entrada.

Antes de abordar el problema de determinar cuando un ajuste no-lineal es mejor que otro, hago otro ejemplo para mostrar los ajustes no-paramétricos. Para este ejemplo uso datos que pueden descargarse desde la página de Robert H.Shumway.

#Serie de tiempo de la mortalidad cardiovascular
mortalidad<-scan("cmort.dat")
t=1:length(mortalidad)
lm(mortalidad~t)
t2=t^2
t3=t^3
c=cos(2*pi*t/52)
s=sin(2*pi*t/52)
fit=lm(mortalidad~t+t2+t3+c+s)
#Graficamos el ajuste
plot(t,mortalidad,type="p",main="Datos de Mortalidad Cardiovascular",xlab="Tiempo",ylab="Cantidad de muertes",lwd="2")
lines(fit$fit,col="2")
abline(a=96.65,b=-0.0312,col="3")
lines(lowess(t,mortalidad),col="4")
lines(ksmooth(t,mortalidad,"normal",bandwidth=104),col="5")
lines(ksmooth(t,mortalidad,"normal",bandwidth=5),col="6")

Regresiones_lineales_y_nolineales

La gráfica del ejemplo anterior muestra varios ajustes, el primero es el lineal, el segundo el polinomial y el tercero y cuarto son métodos no paramétricos que pueden ser considerados como extensiones de la regresión polinomial (regresión no-lineal). Cada uno muestra cierto comportamiento, y puede servir para explorar la información y posteriormente pensar en el tipo de modelo que se puede construir. Con la gráfica anterior solo trato mostrar que las regresiones no se limitan solo a lo lineal y no-lineal, uno puede hacer uso de técnicas no-paramétrica para estimar otro tipo de regresiones. Para mayor detalles se puede revisar la presentación de Andrew W. Moore [3].

Lo que hago ahora es enfocarme a la regresión polinomial, construyo un ejemplo por medio del cual comparo el ajuste lineal y el polinomial.

#Ejemplo no lineal
x<-seq(0,1,by=0.01)
y<-sin(2*pi*x)+rnorm(length(x),0,0.1)
df<-data.frame(X=x,Y=y)
ggplot(df,aes(x=X,y=Y))+geom_point()+ggtitle("Relación no lineal")+theme(plot.title = element_text(lineheight=.8, face="bold"))

Relación_No_Lineal

Se observa en la gráfica que los datos parecen como si estuvieran en una montaña rusa, tienen partes altas y bajas. Pareciera que siguen una trayectoria curva y queda pensar en aplicar dos regresiones, una lineal y otra polinomial.

#Ajuste lineal y polinomial
summary(lm(Y~X,data=df))
Call:
lm(formula = Y ~ X, data = df)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.91351 -0.39942  0.02522  0.39833  0.98606 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.94573    0.09068   10.43   |t|)    
(Intercept)  -0.13568    0.04621  -2.936  0.00415 ** 
X            11.50556    0.40215  28.610  < 2e-16 ***
X2          -33.75407    0.93713 -36.018  < 2e-16 ***
X3           22.55399    0.61585  36.622  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.1204 on 97 degrees of freedom
Multiple R-squared:  0.9726,    Adjusted R-squared:  0.9717 
F-statistic:  1147 on 3 and 97 DF,  p-value: < 2.2e-16

De los ajuste se  puede observar que el lineal explica aproximadamente el 59% y el polinomial con un polinomio de grado 3, explica aproximadamente el 97%. Pero nada limita a considerar polinomios de grado mayor, hago la estimación con un polinomio de grado 14 y  comparo como se comporta el ajuste, para esto usamos la función poly de R.

#Ajuste lineal y polinomial
summary(lm(Y~poly(X,degree=14),data=df))
Call:
lm(formula = Y ~ poly(X, degree = 14), data = df)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.243497 -0.066481 -0.001563  0.066558  0.240361 

Coefficients:
                        Estimate Std. Error t value Pr(>|t|)    
(Intercept)             0.004371   0.009707   0.450    0.654    
poly(X, degree = 14)1  -5.516374   0.097551 -56.549  < 2e-16 ***
poly(X, degree = 14)2   0.058760   0.097551   0.602    0.549    
poly(X, degree = 14)3   4.410330   0.097551  45.211  < 2e-16 ***
poly(X, degree = 14)4  -0.034386   0.097551  -0.352    0.725    
poly(X, degree = 14)5  -0.708767   0.097551  -7.266 1.58e-10 ***
poly(X, degree = 14)6  -0.005587   0.097551  -0.057    0.954    
poly(X, degree = 14)7   0.155736   0.097551   1.596    0.114    
poly(X, degree = 14)8   0.147905   0.097551   1.516    0.133    
poly(X, degree = 14)9   0.022788   0.097551   0.234    0.816    
poly(X, degree = 14)10  0.055761   0.097551   0.572    0.569    
poly(X, degree = 14)11 -0.087618   0.097551  -0.898    0.372    
poly(X, degree = 14)12  0.003417   0.097551   0.035    0.972    
poly(X, degree = 14)13  0.147284   0.097551   1.510    0.135    
poly(X, degree = 14)14 -0.075299   0.097551  -0.772    0.442    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.09755 on 86 degrees of freedom
Multiple R-squared:  0.984,     Adjusted R-squared:  0.9814 
F-statistic: 378.9 on 14 and 86 DF,  p-value: < 2.2e-16

Se observa que el ajuste explicado por el polinomio de grado 14 es del 98% de variabilidad en los datos, por cual puede hacer pensar que es mejor modelos, comparando los 3 que realicé.

Pero nada hace especial el polinomio de grado 14 con respecto a cualquier otro, bien podrías preferir uno de grado 20, 30 o 50.

Si reviso la R cuadrada de varios modelos polinomiales de grado mayor que 14 obtendré que explican del 98% en adelante y claro sin llegar  explicar el 100%. Pero el problema que surge tienen que ver con la siguiente pregunta,¿cómo elegir el mejor modelo de entre un conjunto de modelos?

En Machine Learning la idea no es encontrar un solo modelo, es construir varios, ya que más de una técnica puede funcionar para el mismo tipo de situación. Pero resulta crucial tener un modo o método para  elegir entre todos los modelos posibles.

Cross-Validation y Regularization

“All models are wrong, but  some models are userful”.

George Box

Este tema puedo ser consultado para más detalles y breve en las famosas nota de Andrew Ng, pero también de modo breve puede ser revisado en las notas de Andrew W. Moore. Recomiendo las notas del último, son más amigables ya que las de Andrew Ng requieren formación en matemáticas (probabilidad, álgebra lineal y cálculo de varias variables). Para ver estos dos temas de manera amplia en la referencia[2]

La idea de Cross-Validation es la siguiente:

  1. Elegir dos conjuntos, el de entrenamiento (training set) y el de prueba (test set).
  2.  Estimar el modelo en los datos de entrenamiento.
  3. Calcular MSE o RMSE de los datos de prueba.
  4. Comparar entre los Modelos con datos del paso 3.

En general se elige de todos los datos un 80-20 o 70-30,  es decir; 80% de los datos para entrenamiento y 20% para prueba. Pero se pueden tener algunos problemas con esto, el principal es que si el modelo requiere muchos datos por ser no-paramétrico es posible que la estrategia falle si no se cuenta con suficientes datos. Entonces se hace uso de una técnica conocida como LOOCV (Leave-one-out Cross-Validation), que en lugar de fijar el 20% para prueba solo no se considera uno de los datos, se calcula el modelo y se mide el error del dato no considerado y se hace ese proceso para toda la información, por último se calcula MSE o RMSE con la información de los errores.

Puede parecer que LOOCV cumple bien ante las situaciones donde la cantidad es importante, pero se puede usar otra alternativa, en lugar de hacer solo para una muestra de 20% fija, se divide la muestra en una cantidad finita de subconjuntos. Suponemos en K subconjuntos disjuntos, entonces se hace la estimación del modelo con cada uno de los subconjuntos tomados como el conjunto de prueba y se procede a calcular MSE o RMSE. Esta técnicas se llama K-Fold Cross Validation

La regularización está íntimamente relacionada con el concepto de complejidad en un modelo, se puede pensar que tener un modelo de un polinomio de grado 3 es más complejo que un polinomio de grado 2. Los fundamentos teóricos están relacionados con la noción de un “espacio normado“. No explico a que se refiere el último concepto, lo que interesa en esta entrada es entender cómo validar que un modelo es más complejo que otro.

Como la idea no es desarrollar la teoría, entonces lo que hago es usar una librería en R que realiza la estimación de los modelos lineales generalizados (GLM). Parte de sus algoritmos de estimación se realiza por regularización, así se puede usar uno de sus parámetros de estimación como medida de regularización.

Antes de hacer el ejemplo en R, lo que hago es mostrar como estimamos un modelo de regresión polinomial con respecto a unos datos y analizar el parámetro de regularización por medio de cross-validation para determinar cual es el mejor modelo. Considero que desea usar como modelo la regresión de un polinomio de grado 10, así se elige un conjunto de entrenamiento (training set) y uno de prueba (test set).

#Ejemplo del uso de regularización y cross-validation
x<-seq(0,1,by=0.01)
y<-sin(2*pi*x)+rnorm(length(x),0,0.1)
n<-lenght(x)
indices<-sort(sample(1:n,round(0.5*n))
training.x<-x[indices]
training.y<-y[indices]
test.x<-x[-indices]
test.y<-y[-indices]
df<-data.frame(X=X,Y=y)
training.df<-data.frame(X=training.x,Y=training.y)
test.df<-data.frame(X=test.x,Y=test.y)
#Usamos el conjunto de entrenamiento
glmnet.fit<-with(training.df,glmnet(poly(X,degree=10),Y)
lambdas<-glmnet.fit$lambda
modelo<-dta.frame()
for( lambda in lambdas){
modelo<-rbind(modelo,data.frame(Lambda=lambda,RMSE=rsme(test.y,with(test.df,predict(glmnet.fit,poly(X,degree=10),s=lambda)))))
}
#Visualizamos la relación entre Lambdas vs RMSE
ggplot(modelo,aes(X=Lambda,Y=RMSE))+geom_point()+geom_line()
#Elegimos Lambda que minimiza RMSE
best.lambda<-with(modelo, Lambda[which(RMSE==min(RMSE))])
glmnet.fit<-with(df,glmnet(poly(X,degree=10),Y))
#Visualizamos los coeficientes estimados del "mejor" modelos
coef(glmnet.fit,s=best.lambda

Del código anterior y del proceso que se realiza, se puede concluir que entre las dos herramientas Regularización y Cross-Validation se puede elegir el modelo que tienen menos “desviación estándar” en la distribución de los errores.

Se aprecia en la gráfica de abajo que el mínimo de la curva  ayuda a determinar el mejor modelo.

Lambda_vs_RMSE

 Caso de estudio

El siguiente ejemplo forma parte de los ejemplos presentados en el texto de Machine Learning for Hackers[1]. Lo que se hace en el texto es tomar datos de los 100 libros más populares de la editorial O’Reilly y tratar de definir partiendo de los textos de las contra portadas el posible éxito o alta venta.

Para desarrollar el ejemplo se usan las librerías tm y glmnet. Por medio de la primera se construye una matriz mediante los textos de la contra portada y ya teniendo la matriz se usa la regresión lineal general contra una variable para tratar de definir un rango de ventas partiendo de la lista de palabras que  contiene.

Mediante de la regularización (lambda) y la cross-validation, se trata de elegir el modelo que minimice RMSE con respecto al parámetro lambda del modelo.

#Análisis de contra-portadas
rango<-read.csv(data.file,stringsAsFactor=FALSE)
documentos<-data.frame(Texto=rango$Long.Desc.)
row.names(documentos)<-1:nrow(documentos)
library(tm)
corpus<-Corpus(DataframeSource(documentos))
corpus<-tm_map(corpus, content_transformer(tolower))
corpus<-tm_map(corpus,stripWhitespaces)
corpus<-tm_map(corpus,removeWords,stopwords('english'))
dtm<-DocumentTermMatrix(corpus)

rendimiento<-data.frame()

for(lambda in c(0.1,0.25,0.5,1,2,2.5))
{ 
  for ( i in 1:50){

    indices<-sample(1:100,80)
    training.x<-x[indices,]
    training.y<-y[indices]
    test.x<-x[-indices,]
    test.y<-y[-indices]
    glmnet.fit<-glmnet(training.x,training.y)
    prediccion.y<-predict(glm.fit,test.x,s=lambda)
    rmse<-sqrt(mean((prediccion.y-test.y)^2))
    rendimiento<-rbind(rendimiento,data.frame(Lambda=lambda,Interaccion=i,RMSE=rmse))
    }
}

El código anterior hace los siguiente:

  1. Se define el rango del parámetro de regularización (lambda).
  2. Se eligen subconjunto, 80% para entrenamiento y 20% para probar.
  3. Se estima la regresión y se compara contra el conjunto de prueba, con esta comparación se calcula RMSE para cierto lambda.

Análisis_Text_RMSE vs Lambda

Análisis_Text_RMSE vs Lambda2

 

La primera gráfica muestra que no se logra definir un buen modelo. En general RMSE vs Lambda nos indican que no se puede obtener mucho de los datos analizados, es decir; no resulta evidente como poder distinguir entre los 50 libros más vendidos con ese modelo.

Por lo cual se hace uso de una variable dummy para pasar de una regresión a una clasificación, por medio de la cual se puede responder cuando un libro es considerado uno de los 50 más vendidos solo partiendo de su contra portada. En este ejemplo no se trata solo de hacer el modelo, sino de verificar mediante regularización y cross-validation como es posible elegir un modelo con menor complejidad o uno mejor entre los posibles

Las observaciones y pasos a seguir son los siguientes; al usar la variable dummy y aplicar una regresión logística, se eligen subconjuntos de entrenamiento y prueba. Después de varias interacciones se puede analizar el comportamiento del parámetro lambda y  en lugar de RMSE  se usa la tasa de errores entre las predicciones y el subconjunto de prueba.

#Análisis de contra-portadas
y <- rep(c(1, 0), each = 50)
regularizacion.fit <- glmnet(x, y, family = 'binomial')
predict(regularizacion.fit, newx = x, s = 0.001)
ifelse(predict(regularizacion.fit, newx = x, s = 0.001) > 0, 1, 0)

performance <- data.frame()

for (i in 1:250)
{
 indices <- sample(1:100, 80)
 
 training.x <- x[indices, ]
 training.y <- y[indices]
 
 test.x <- x[-indices, ]
 test.y <- y[-indices]
 
 for (lambda in c(0.0001, 0.001, 0.0025, 0.005, 0.01, 0.025, 0.5, 0.1))
 {
 glm.fit <- glmnet(training.x, training.y, family = 'binomial')
 
 predicted.y <- ifelse(predict(glm.fit, test.x, s = lambda) > 0, 1, 0)
 
 tasa.error <- mean(predicted.y != test.y)

 performance <- rbind(performance,data.frame(Lambda = lambda,Iteration = i,TasaError = tasa.error))
 }
}

ggplot(performance, aes(x = Lambda, y = ErrorRate)) +
 stat_summary(fun.data = 'mean_cl_boot', geom = 'errorbar') +
 stat_summary(fun.data = 'mean_cl_boot', geom = 'point') +
 scale_x_log10()

 Análisis_Text_Tasa_Error vs log_ Lambda

En conclusión, se observa que al realizar el proceso con una regresión logística, se puede analizar el comportamiento de lambda y la tasa de error. Por lo tanto se puede clasificar los textos y elegir un modelo con la tasa de error correspondiente a lambda menor (como en la gráfica anterior).

¿Por qué el último modelo permite hacer una elección entre ellos y el anterior no?

El último modelo muestra constancia en la tasa de error, entonces se puede considerar que la menor Lambda hace considerar al modelo como el mejor.

En el modelo previo no se veía un patrón definido, elegir el que tiene menor RMSE o la menor Lambda no asegura que se este haciendo una buena elección.

Lo anterior puede parecer superficial, pero pensándolo un poco no es  una mala estrategia proceder con una variable dummy para modificar el modelo y tomar las tasas de error en lugar de RMSE.

Espero la entrada de una idea general de como se usa la regularización y la cross-validation, también como el modelo no-lineal puede beneficiarnos en algunos casos. En otra ocasión espero extender los ejemplo y hacer otro que puedan ayudar a clarificar los conceptos.

Algunas cosas que desearía fortalecer de esta entrada más adelante es resaltar con otro ejemplo como la regresión logística toma un rol de método de clasificación lineal, el rol de la regularización y la validación-cruzada o cross-validation para elegir modelos. Un par de buenos artículos sobre el tema se pueden ver en las referencias [2,5,6,7].

 Referencias:

1.-http://shop.oreilly.com/product/0636920018483.do

2.-http://statweb.stanford.edu/~tibs/ElemStatLearn/

3.-http://openclassroom.stanford.edu/MainFolder/DocumentPage.php?course=MachineLearning&doc=exercises/ex5/ex5.html

4.-http://www.amazon.com/Machine-Learning-Probabilistic-Perspective-Computation/dp/0262018020

5.-http://robjhyndman.com/hyndsight/crossvalidation/

6.-http://stats.stackexchange.com/

7.-http://projecteuclid.org/download/pdfview_1/euclid.ssu/1268143839

Anuncios

3 comentarios sobre “Regresión no lineal, Cross-Validation y Regularization.

Responder

Por favor, inicia sesión con uno de estos métodos para publicar tu comentario:

Logo de WordPress.com

Estás comentando usando tu cuenta de WordPress.com. Cerrar sesión / Cambiar )

Imagen de Twitter

Estás comentando usando tu cuenta de Twitter. Cerrar sesión / Cambiar )

Foto de Facebook

Estás comentando usando tu cuenta de Facebook. Cerrar sesión / Cambiar )

Google+ photo

Estás comentando usando tu cuenta de Google+. Cerrar sesión / Cambiar )

Conectando a %s