Regresión Lineal, un ejemplo sencillo.

Un poco de la regresión

La regresión lineal es quizás el método más conocido para “predecir” el comportamiento de los datos o intentar hacerlo. Es como el caballito de batalla en los métodos predictivos, en  Machine Learning (ML) es considerada como una de las técnica centrales del aprendizaje supervisado.

Siguiendo las ideas de ML siempre se tienen datos de entrada X y de salida Y, el método que se usa para analizar los datos, es la regresión lineal cuando los datos de salida son “números reales“, como 1.2, 3.1415, 2.768, etc.

Pero la técnica no solo se estudia en ML, es de los principales temas que se enseñan en estadística, econometría o en ciencias sociales y políticas.

La idea de ajustar los datos experimentales a una curva, no se limita a una recta,  en el peor de los casos puede ser una curva más elaborada, por ejemplo polinomios. Para la gente que no tiene conocimientos de lo que es un polinomio, basta decir que es una cosa muy poco lineal y el costo de ajustar por un polinomio es sobre estimar los datos.

Desde los primeros años, entre la secundaria y la preparatoria, nos enseñan a graficar parejas de números, pensando en que se tienen valor x y la función f(x). En el problema de la regresión y en general de ML,  solo se tienen las parejas de datos (x,y) y la cuestión es definir  una función f(x) donde los valores sean muy cercanos a y.

En el caso de la regresión lineal, la función f(x)=ax+b es la que se considera “mejor” para describir la relación y el siguiente problema que surge es cómo elegir “b” y “a” de modo que entre todas nuestras posibles rectar sea la “mejor” para el modelo.

Entonces el problema se reduce a elegir un modelo y estimar sus parámetros, pero teniendo una estrategia para discernir entre los que serán los mejores parámetros.

Hago un ejemplo fácil con datos de la temperatura global. Los datos se pueden descargar desde la página de Robert H. Shumway.

#Cargamos los datos
temperaturas.globales<-scan("data/globtemp.dat")
x=temperaturas.global[45:142]
#Asignamos los tiempos de nuestros datos
t=1900:1997
#Estimamos la recta que describe "mejor" nuestros datos
fit=lm(x~t)
plot(t,x,type="o",col="2",xlab="Años",ylab="Temperaturas globales",main="Ejemplo de Regresión Lineal")
abline(fit)

Regresion_temperaturas_globales

La imagen muestra como las temperaturas  tiene “cierto” comportamiento lineal con respecto al tiempo, es decir, siguen una tendencia. Se observa que en el código anterior el comando lm es el que estima los valores de la recta, más adelante explico qué es lo que informa este comando y cómo interpretar esa información.

Regresando el modelo lineal, una vez elegida la recta, es valido preguntar si el ajuste es “bueno” o qué tanto explica la recta el comportamiento de los datos.

Hago otro ejemplo, considero datos de fumadores y no fumadores y sus respectivas edades de muerte. Como buena práctica primero realizo una exploración y después estimo el modelo lineal.

#Cargamos los datos
edades<-read.csv("data/longevity.csv")
head(edades)
#Les cambiamos el nombre a las columnas
names(edades)<-c("Fumadores","Edad_de_Muerte")
#Graficamos explorando el comportamiento de la densidad entre hombres y mujeres.
ggplot(edades2, aes(x = Edad_de_Muerte, fill = factor(Fumadores)))+geom_density() + facet_grid(Fumadores ~ .)

Comparacion_entre_fumadores

Lo que se observar y  se concluye a primera vista, es que parece que la gente que no fuma vive más años, pero la cuestión es cómo justificar la predicción o el descubrimiento al explorar los datos.

Avanzando un poco a justificar lo que se ve en la gráfica. Calculo la media de la edad de muerte y se tiene que es 72.72 años (la media se puede calcular con el comando mean). Se puede pensar que la media de la edad es aproximadamente  73 años y  se usa este número como la edad pronosticada.

Lo que trato de hacer para justificar la elección de la predicción, es calcular lo que se llama MSE (mean squared error), error cuadrático medio, que surge de comparar la predicción contra los valores de los datos, datos menos predicción. A estos nuevos datos (los errores) se les calcula su media.

El MSE tienen una explicación estadística que en otro momento complementaré, en breve se puede decir que es el segundo momento de los errores.

#Calculamos MSE
predicion<-73
with(edades,mean(Edad_de_Muerte-predicion)^2))

Lo que obtengo es un valor de 32.99, esto por si solo no dice nada respecto a la propuesta de pronóstico de edad de muerte. Lo que busco es que la edad que propongo minimice la diferencia entre los datos y el pronóstico, entonces pensando en que la edad está entre los 60 y los 85 años comparo si alguna otra edad minimiza mejor el MSE contra la propuesta de 73 años.

#Calculamo el MSE para varios pronósticos propuestos
presicion.pronostico<-data.frame()
for(predicion in seq(60,85,by=1){
error.prediccion<-with(edades,mean(Edad_de_Muerte-predicion)^2))
presicion.pronostico<-rbind(presicion.pronostico,data.frame(Pronóstico=pronostico,Error=error.prediccion))}
ggplot(presicion.pronostico,aes(x=Pronóstico,y=Error))+geom_point()+geom_line()

Comparacion_de_predicion

Se ve en la gráfica que la predicción 73 años, no es tan mala como valor medio de edades de muerte, de hecho está casi en la parte más baja de la curva, lo cual indica que es un óptimo.

Además  se observa que el numero 32.99 es complicado de interpretarlo solo, haciendo lo anterior toma sentido decir que la predicción no es tan mala y que ese 32.99 se refiere al valor donde se minimiza el error.

Ahora analizo de nuevo la información  pero teniendo en cuenta las personas  que son fumadoras.

Si uno calcula la regresión, el modelo que se obtiene es Edad_de_Muerte=75.254-5.062*Fumadores+ a, entonces se puede considerar que la variable Fumadores es categórica y “a” es una constante. También uno puede pensar en qué sucede si se calcula el promedio solo de los fumadores y el de los no fumadores, para comparar hago una gráfica de las densidades.

#Comparación RMSE
media<-with(edades,mean(Edad_de_Muerte))
error.cuadratico_medio<-with(edades, sqrt(mean((Edad_de_Muerte-media)^2)))
#Consideramos las medias de fumadores y no fumadores
media.fumador<-with(subset(edades,Fumador==1),mean(Edad_de_muerte))
media.nofumadores<-with(subset(edades, Fumador==0),mean(Edad_de_Muerte))

#Ahora calculamos el RMSE por categoría

edades<-transform(edades,NuevaPrediccion=Ifelse(Fumador==0,media.nofumador,media.fumador))

error.cuadratico_categorias<-with(edades,sqrt(mean((Edad_de_Muerte-NuevaPrediccion)^2)))

Si se compara la medición RMSE, se nota que considerar la información de los fumadores, modifica la calidad de la predicción, casi en un 10%.

Aun que no se dije quien es RMSE ( root mean squared error), basta decir que es la raíz cuadrada de MSE y que funciona como la desviación estándar del error de los valores predichos con respecto a los datos, es una medida de dispersión del error entre el modelo y los datos.

El ejemplo anterior parece confuso a primera vista, lo que hice fue como primer paso explorar las densidades de los datos por categorías, lo segundo fue calcular la media de todos los datos, después validar de manera gráfica que la media era una buena predicción de la edad y posteriormente calcular la regresión de las Edades contra ser Fumador o no. Esta ultima variable solo tienen valores, 0 y 1, el modelo lineal da valores aproximados a las medias de las edades de cada grupo, fumadores y no fumadores. Como ultimo paso uso RMSE para compara cuanto cambia la edad predicha entre categorías.

Antes de hacer un ejemplo más interesante, regreso a la pregunta ,¿cómo saber si el modelo es bueno ?, eso es sutil y las herramientas estadísticas por si solas no pueden confirmarlo en todos los casos.

La ventaja es que en la regresión lineal existe el coeficiente  R cuadrada, también es conocido como coeficiente de determinación, el cual es la proporción de la variación  total de nuestras Y explicadas por la regresión con X.

Este índice se encuentra entre 0 y 1, es una herramienta rápida y la regla es que un valor cercano a 1 es un buen ajuste y cercano a 0, es malo el ajuste. Este índice es solo un estadístico descriptivo, por ello solo sirve para poder evaluar de manera rápida los modelos.

Un valor cercano a 0  puede deberse a diversos factores, por lo cual es importante explorar la información en busca de posibles anomalías y probar otro tipo de modelos que no sean lineales.

Hasta ahora lo dicho se resume en lo siguiente:

  1. Se puede construir un modelo que relacione de manera lineal los datos de entrada con los datos de salida.
  2. Proponer la media de los datos observados, para la cual se calculan los errores entre la media propuesta y los datos observados.
  3. Teniendo los errores observados se calcula su media (MSE) y la desviación estándar de los errores (RMSE). Se puede considerar que la media propuesta es mejor si el RMSE es menor.
  4. Se tiene un índice que de manera descriptiva indica si el modelo lineal es bueno o malo, el cual ayuda a validar la calidad del modelo.

Predicción del tráfico en la Web

En los ejemplos anterior solo puse la relación entre dos variables, pero la regresión se puede generalizar a más de una variable para explicar el comportamiento de los datos de salida.  Ha este tipo de modelos se les llama, modelos de regresión lineal múltiple.

En mi opinión son más natural los modelos de regresión lineal múltiple, es decir; Y=a+bX1+cx2+dX3….+zXn. Principalmente creo que es complicado que una sola variable explique el comportamiento de otra, por lo cual la regresión múltiple creo que es más natural.

Sigo el ejemplo de Neil Kodner, el cual revisa los 1000 sitios web más visitados en el 2011, este ejemplo es discutido en la referencia [1]. Para replicar el código los datos se pueden descargados desde GitHub.

#Cargamos los datos
data.file<-file.path("data","top_1000_sites.tsv")
#Pedimos que no lea las cadenas como Factores
top.1000.sites<-read.csv(data.file,sep="/t",stringsAsFactor=FALSE)
#Visualizamos los primeros 30 sites
head(top.1000.sites,n=30)

Hago dos observaciones de la carga de los datos, la primera es que se pide que sean separados por tabuladores y la segunda que las cadenas no sean tomadas como “factores”. Un factor es un vector usado por R que clasifica de manera discreta los datos de otro vector de la misma longitud, se puede revisar más información en el manual de R. Ejemplo, si cargamos un dato con los valores (1,2,3,1,2,3,3,4), si los leemos como factores lo que dará R  es un factor con esos valores y con los niveles (1,2,3,4) lo cual dice los manejará con un orden.

Se observa, mediante el comando head que los datos tienen ciertas cabeceras en las columnas que describen en general la información y se observa que no toda la información de las columnas es numérica. Las únicas columnas que tienen interés ser estudiadas son: PageViews, UniqueVisitors, HasAdvestising y IsEnglish.

Primero se explora la relación entre UniqueVisitors y PageViews, lo que observa al graficar las dos variables es que no muestra alguna relación lineal, se ven solo puntos cercanos al origen de la gráfica. Lo que hago como segundo paso es transformar los datos por medio del logaritmo, es decir se hace una gráfica de log(PageViews) vs log(UniqueVisitors).

Lo anterior cambia la escala y permite ver los datos como si uno se acercara con una lupa, con este cambio se observa ahora una relación lineal entre las dos variables.

#Graficamos las variables
ggplot(top.1000.site,aes(x=PageViews,y=UniqueVisitors))+geom_point(colour="2")

#Transformamos los valores de las variables
ggplot(top.1000.site,aes(x=log(PageViews),y=log(UniqueVisitors))+geom_point(colour="2",size=3)

#Ahora agregamos la recta de la regresión
ggplot(top.1000.site,aes(x=log(PageViews),y=log(UniqueVisitors)))+geom_point(colour="2",size=3)+geom_smooth(method='lm',size=3)

Una justificación de por qué transformar los datos, es decir; considerar log(PageViews) en lugar de solo PageViews, tiene que ver con el rango y basta ver que la densidad no se comportan como una distribución normal en la primera de las siguientes gráficas, hasta complica un poco su visualización. Ese tipo de comportamientos implican que puede ser útil transformar los datos a otra escala (segunda gráfica).

Denisidad_paginasvistas

Denisdad_log_paginasvistas

Paginasvistas_vs_visitas

log_Paginasvistas_vs_log_visitas

Tendencia

El análisis gráfico insinúa que existe una relación entre las dos variables. Realizando una regresión se encuentra hasta qué punto una variable explica la variabilidad de la otra.

Para el ejemplo que muestro no explico formalmente qué pasa con el cálculo de la regresión, ni tampoco entro en detalles formales de qué significan cada estadístico, solo explico la idea intuitiva de qué se obtiene al aplicar la regresión en R project. El siguiente código es la respuesta de R con el comando summary.

#Regresión lineal    
summary(lm.fit)

Call:
lm(formula = log(UniqueVisitors) ~ log(PageViews), data = top.1000.site)

Residuals:
 Min 1Q Median 3Q Max 
-1.95645 -0.35798 -0.04481 0.27530 2.34629 

Coefficients:
             Estimate Std. Error t-value Pr(>|t|) 
(Intercept)    9.83288 0.22670   43.37    <2e-16 ***
log(PageViews) 0.34544 0.01181   29.25    <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.551 on 998 degrees of freedom
Multiple R-squared: 0.4616, Adjusted R-squared: 0.4611 
F-statistic: 855.6 on 1 and 998 DF, p-value: < 2.2e-16

Hay dos aspectos para revisar de la regresión observada en el código anterior. El primero es el significa de lo que se obtiene del software  y el segundo, que uno está obligado aprender un poco, es el análisis de incertidumbre.

El último es fundamental, pero requiere entrar en detalles teóricos, por lo cual recomiendo acudir a un buen texto para consultarlo, recomiendo el libro de A.Gelman [2], pero también en la referencia [3] se cubre todo lo necesario desde una perspectiva de teoría del aprendizaje.

En breve, Estimate da el valor de los coeficiente en la regresión y Std.Error da el error estándar del coeficiente calculado. Lo que se denota como t-valor tienen que ver con la proporción o el comportamiento de coeficiente estimado y rango de error que se tienen, es decir; si se calcula 1.33628/0.04568=29.251 es igual al t-valor. Ese t-valor está relacionado con las desviaciones estándar del coeficiente estimado, tener un un t-valor pequeño habla de que posiblemente el coeficiente no es relevante  y que es viable considerar el valor del coeficiente como cero .

Como prueba estadística se considera que si Pr(>|t|) es grande uno debe de rechazar la relación entre las variables, es decir; el coeficiente puede ser tomado como cero.

En el ejemplo uno puede estar tranquilo, existen pistas estadísticas para pensar en que la relación lineal no es un mal ajuste, ya que el Pr(>|t|) es muy pequeño, así que los coeficientes “a” y “b” de la regresión no son cero.

Por último  Multiple R-squared y Adjusted R-squared,  hablan de la varianza explicada de los valores de  Y, lo cual de cierto modo indica que tan bueno es el modelo.

Anteriormente expliqué lo que es la R cuadrada, que el ejemplo su valor es de 0.4616. Esto indica que la variable log(PageViews) explica el 46% la variabilidad de la variable log(UniqueVisitors).

El aspecto que me falta mencionar es el estudio de los residuos, todos los estadístico comentados anteriormente están relacionados con ellos, pero uno puede analizarlos y lo deseable teóricamente es que tengan distribución normal para el modelo lineal estimado. En general no pasa eso y para ello existen pruebas estadísticas o herramientas gráficas que permiten confirmar cómo se distribuyen los residuos.

Una gráfica sencilla es revisar los cauntiles con respecto al comportamiento de los cuatiles de una distribución normal. El ejemplo gráfico es el siguiente:

Qqplot

Se ve que los residuos no tienen distribución normal, sobre todo en las “colas” de las distribuciones se desvían mucho del comportamiento deseable, el comportamiento deseable se muestra  en la gráfica con la línea roja.

Se vio que la regresión anterior solo explica el 46%, lo cual no es malo. Ahora hago una regresión multilineal considerando además de la variable log(UniqueVisitors), también las variables HasAdvertising y InEnglish.

#Regresión multilineal
lm.fit=lm(log(UniqueVisitors)~log(PageViews)+InEnglish+HasAdvertising,data=top.1000.site)    

#Revisamos los estadísticos de nuestra regresión multilineal

summary(lm.fit)

#Observamos los valores

Call:
lm(formula = log(UniqueVisitors) ~ log(PageViews) + InEnglish + 
 HasAdvertising, data = top.1000.site)

Residuals:
 Min        1Q      Median    3Q     Max 
-1.22011 -0.30299 -0.06927 0.26348 1.42979 

Coefficients:
                Estimate  Std.Error t-value Pr(>|t|) 
(Intercept)       12.63974 0.20283 62.318 <2e-16 ***
log(PageViews)    0.19314 0.01077 17.936 <2e-16 ***
InEnglishNo       1.15772 0.07351 15.750 <2e-16 ***
InEnglishYes      1.57684 0.06224 25.335 <2e-16 ***
HasAdvertisingYes -0.03810 0.03601 -1.058 0.29 
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.4168 on 995 degrees of freedom
Multiple R-squared: 0.6928, Adjusted R-squared: 0.6915 
F-statistic: 560.9 on 4 and 995 DF, p-value: < 2.2e-16

Se observa que agregando esas dos nuevas variables solo contribuyen aproximadamente con 1% para explicar la variabilidad de los datos. Por lo cual uno puede observar la poca relevancia que tienen las nuevas variables al ver su t-valores y sus p-valores, Pr(>|t|). En general el modelo mejora un poco con respecto al de una sola variable. Uno también podría revisar los residuos y analizar su distribución o aplicar alguna prueba estadística para ver que tanto se puede considerar que siguen una distribución normal, lo cual no hago en este ejemplo.

Depende del investigador y de las preguntas que se busque responder para considerar otras variables o construir alguna que pueda tener relevancia. Se suelen agregar variables Dummy para conocer si impacta alguna categorización de los datos o no, en el ejemplo anterior la variable InEnglish puede verse como variable Dummy.

Como siempre, no existe un solo modelo posible, ni tampoco basta un solo modelo para explicar cualitatívamente lo que sucede, tratar de entender un fenómeno es algo más delicado y no solo requiere aplicar alguna técnica sobre la información que se tiene, se requiere profundizar más en el tema o fenómeno.

Varios puntos faltan mencionar en la entrada, entre ellos los aspectos a revisar para como fallas comunes en los modelos de regresión. Esto implica el comportamiento de los outlier como impactan en el modelo, la multi-colinealidad, la homocedasticidad. Otro aspecto es que uno puede agregar y agregar variables predictoras para modelar la variable de respuesta, pero uno puede usar algunas técnicas para probar con varios modelos y “reducir dimensiones”; es decir, usar menos variables predictoras que satisfagan condiciones estadísticas necesarias para tener un buen modelo lineal.

Correlación

Me faltó mencionar un indicador o medida de linealidad entre variables, la correlación; de la correlación surge el coeficiente de correlación y es una medida de qué tan alta es la relación lineal entre dos variables. En R, es muy simple calcularla, pero antes de mostrar un ejemplo cabe resaltar que si las variables son “independientes” su correlación es cero y que la existencia de correlación no implica la existencia de causalidad.

Esto último es sumamente importante, por que un puede ponerse hacer cálculos de la correlación sin tener ninguna relevancia y creer que se encontró una explicación a tal o cual fenómeno, lo cual en la mayoría de los casos es falso.

Una divertida página que muestra cosas absurdas al evaluar la correlación es la siguiente: spurious correlations. Es altamente recomendada para reírse entre las relaciones encontradas.

Hago un ejemplo de como calcular la correlación en R project.

#Ejemplo de cálculo de la correlación
x<-rnorm(200)
y<-x+rnorm(200)
ggplot(data.frame(X=x,Y=y),aes(x=X,y=Y))+geom_point()+geom_smooth()+ggtitle("Ejemplo de Correlación")+theme(plot.title = element_text(lineheight=.8, face="bold"))
cor(x,y)

El valor debe de ser mayor a 0.5 lo cual indica que se tienen una correlación positiva y la gráfica muestra que la línea tienen una pendiente positiva, lo cual va acorde con la correlación entre las variables.

Ejemplo_Correlación

Espero que estos pocos ejemplos y explicaciones inciten a buscar otras fuentes o hacer otra serie de ejemplos sobre la regresión lineal.

En lo comentado no hice mención de temas y métodos actuales de la regresión lineal (Lasso,Rigde, elastic Net), estos se pueden leer en la referencia [3,4,5].

Referencias:

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

2.-http://www.stat.columbia.edu/~gelman/arm/

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

4.-http://statweb.stanford.edu/~owen/courses/305/Rudyregularization.pdf

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

Anuncios

3 comentarios sobre “Regresión Lineal, un ejemplo sencillo.

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