Clase 7 - Modelos de regresión

Modelos y datos

A partir de la clase de hoy, nos adentraremos en la metodología central de la ciencia de datos: la creación de un modelo estadístico. Un modelo es una representación compacta de los datos y, por lo tanto, de los procesos subyacentes que producen los datos. Así como los arquitectos crean modelos de edificios para comprender mejor cómo se distribuyen los diversos componentes de un edificio, los científicos de datos crean modelos estadísticos para obtener una comprensión análoga de los datos.

Hay dos tipos de tareas de modelado estadístico: descriptivas y predictivas. En el modelado descriptivo, el objetivo es identificar las relaciones entre las variables para aprender más sobre la estructura de los datos. En el modelado predictivo, el objetivo es predecir el valor de una variable dado el resto de las variables. La variable que se va a predecir se denomina variable dependiente o objetivo, y las variables utilizadas para la predicción se denominan variables independientes o predictoras. Las tareas de modelado descriptivo y predictivo a menudo van de la mano: podemos utilizar la misma técnica de modelado para ambos tipos de análisis.

Sin embargo, todos los modelos no son equivalentes; al ser una representación compacta, un modelo es casi siempre una representación simplificada de los datos. Para cada modelo, existe un supuesto de modelado subyacente, por ejemplo, en el caso de modelos lineales, asumimos una relación lineal entre variables. Estos supuestos imponen limitaciones al modelo en términos de los casos en los que podemos aplicarlo y cómo interpretar sus resultados. Es importante conocer las suposiciones del modelo, al igual que conocer las advertencias de efectos secundarios en la etiqueta de un medicamento médico.

En la clase del día de hoy veremos uno de los enfoques de modelado estadístico más utilizados: la regresión. En la siguiente semana veremos la clasificación. En la regresión, la variable objetivo es numérica, por ejemplo, pulgadas de lluvia durante el día, mientras que en la clasificación, la variable objetivo es categórica, por ejemplo, pronóstico del tiempo para el día (soleado, nublado, lluvia, nieve). Existe una larga lista de enfoques de modelado para la regresión; veremos algunos de los aspectos fundamentales y de uso frecuente.

Análisis de regresión

El análisis de regresión implica identificar relaciones entre variables. Dado un conjunto de variables independientes o predictoras, nuestro objetivo es predecir una variable dependiente o objetivo. En el análisis de regresión, la variable objetivo es numérica y las variables predictoras pueden ser numéricas, categóricas u ordinales.

En última instancia, la regresión implica aprender una función que calcula el valor predicho de la variable objetivo dadas las variables predictoras. Convencionalmente, usamos los símbolos x1,. . . , xm para denotar las variables predictoras e y para denotar la variable objetivo. Suponiendo que el conjunto de datos contiene n puntos de datos, estas variables serían n vectores dimensionales: y = (y1,…, Yn).

La mayoría de las técnicas de regresión, y la mayoría de las técnicas de modelado en general, implican un proceso de dos pasos:

Ajuste del modelo: para un conjunto de datos dado con variables x1,. . . , xm e y calculamos los parámetros del modelo w que mejor cumplen con ciertos criterios.
Modelo de predicción: dado un conjunto de datos con las variables predictoras x1,. . . , xm y los parámetros del modelo w, calculamos el valor de la variable objetivo y. El conjunto de datos podría ser el mismo que usamos para el ajuste del modelo o también podría ser un nuevo conjunto de datos.

Dependiendo de la elección de la función f(), existen cientos de técnicas de regresión, cada una con un método para el ajuste y la predicción del modelo. La familia de técnicas más comúnmente utilizada es la regresión lineal, donde asumimos que f() es una función lineal. En esta sección, primero examinamos algunas de las técnicas de regresión lineal y luego pasaremos a los métodos de regresión no paramétrica. Una ventaja de usar R es que la mayoría de las funciones de regresión tienen interfaces similares, por lo que podemos alternar entre técnicas sin demasiado esfuerzo. Esto se debe a que la mayoría de las funciones de regresión tienen interfaces similares en términos de sus entradas.

Además, el día de hoy empezaremos a trabajar con el paquete caret() para realizar regresiones lineales. Además, este paquete nos ayudará durante las próximas semanas.

Cargue de datos

Lo primero que hacemos cada una de las clases es revisar el cargue de nuestros datos;

data <- read.csv("Clase 7.csv", sep=";")
str(data)

## 'data.frame':    352 obs. of  9 variables:
##  $ Marca     : chr  "A Speyside Distillery" "Ainslie's" "Angel's Envy" "Asyla" ...
##  $ Pais      : chr  "Scotland" "Scotland" "United States" "Scotland" ...
##  $ Whiskies  : int  11 32 23 11 12 364 11 12 1071 62 ...
##  $ Votos     : int  10 12 13 7 4 249 11 3 908 26 ...
##  $ Puntuacion: num  86.7 82.3 83.4 81.8 81 ...
##  $ Categoria : chr  "C" "D" "D" "D" ...
##  $ Sabor     : int  8 5 3 3 9 6 3 7 8 7 ...
##  $ Color     : int  6 6 6 2 10 2 7 9 5 9 ...
##  $ Maltas    : int  8 8 1 8 8 1 9 1 1 9 ...

El día de hoy trabajaremos con una base de datos muy sencilla, la cual nos muestra las diferentes características de un grupo de whiskies de diferentes países del mundo. La variable que le interesa estimar es el puntaje del whisky. El link para acceder a los datos es el siguiente: https://www.dropbox.com/s/2c9tr2jcf1wp550/Clase%207.csv?dl=0.

Regresiones lineales simples

Empezaremos con el modelo de regresión más sencillo, la regresión lineal simple. La regresión lineal simple es un caso especial de regresión lineal con una sola variable predictiva x y una variable objetivo y. Para realizar una regresión usaremos el comando lm(). Esta es una función que nos permitirá hacer uso de diferentes variaciones de modelos de regresión lineal (Nota: solo se encuentra disponible la versión con mínimos cuadrados ordinarios). La forma de utilizarlo es el siguiente:

reg1 <- lm(data$Puntuacion ~ data$Maltas)
reg1

## 
## Call:
## lm(formula = data$Puntuacion ~ data$Maltas)
## 
## Coefficients:
## (Intercept)  data$Maltas  
##   83.023531    -0.001233

Vamos a empezar viendo el detalle general del primer modelo de regresión. Como podrán ver, el intercepto corta en el punto 83 y el beta de la variable es -0.001, indicando que entre más maltas menor es el puntaje. Pero esto no nos da mucho detalle sobre el modelo en general, por lo que usaremos la función summary()para revisarlo:

summary(reg1)

## 
## Call:
## lm(formula = data$Puntuacion ~ data$Maltas)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -80.022  -2.462   1.082   3.183  10.228 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 83.023531   0.584675 141.999   <2e-16 ***
## data$Maltas -0.001233   0.099337  -0.012     0.99    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.162 on 350 degrees of freedom
## Multiple R-squared:  4.4e-07,    Adjusted R-squared:  -0.002857 
## F-statistic: 0.000154 on 1 and 350 DF,  p-value: 0.9901

Lo primero que podremos notar es que el p-valor de los betas es no significativo, lo que nos indica que el modelo no es funcional en su esencia. Además, podemos notar que el R cuadrado es bajo y la regresión es no significativa. Si queremos ver que tan bien nos fue con los datos, usando un gráfico de dispersión de nuestras variables podemos añadir nuestra línea de regresión con la función abline().

plot(data$Puntuacion ~ data$Maltas)
abline(reg1)

Parece que el modelo no nos funciona tan bien como esperábamos, por lo que empezaremos a generar cambios sobre el modelo o sobre la función para probar diferentes escenarios.Si queremos una regresión que pase por el origen (sin intercepto) deberemos usar la siguiente expresión:

reg2 <- lm(data$Puntuacion ~ data$Maltas - 1)
summary(reg2)

## 
## Call:
## lm(formula = data$Puntuacion ~ data$Maltas - 1)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -43.44 -12.76  18.92  71.16  79.49 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## data$Maltas  11.6686     0.4266   27.35   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 47.11 on 351 degrees of freedom
## Multiple R-squared:  0.6807, Adjusted R-squared:  0.6798 
## F-statistic: 748.1 on 1 and 351 DF,  p-value: < 2.2e-16

plot(data$Puntuacion ~ data$Maltas)
abline(reg2)

Este modelo parece funcionar mejor que el anterior, pero miraremos nuevas variaciones ya que el arreglo es artificial y realmente no nos ayuda a predecir. Si intentamos incluir una variable discreta como predictor podemos hacerlo directamente, y no necesitamos crear variables indicadoras individuales:

reg2 <- lm(data$Puntuacion ~ data$Categoria)
summary(reg2)

## 
## Call:
## lm(formula = data$Puntuacion ~ data$Categoria)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -52.728  -0.900  -0.003   1.105  14.192 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       91.448      1.388  65.880  < 2e-16 ***
## data$CategoriaC   -4.685      1.418  -3.304  0.00105 ** 
## data$CategoriaD   -8.611      1.418  -6.074  3.3e-09 ***
## data$CategoriaE  -13.350      1.465  -9.115  < 2e-16 ***
## data$CategoriaF  -18.817      1.726 -10.905  < 2e-16 ***
## data$CategoriaG  -35.720      1.963 -18.196  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.4 on 346 degrees of freedom
## Multiple R-squared:  0.699,  Adjusted R-squared:  0.6947 
## F-statistic: 160.7 on 5 and 346 DF,  p-value: < 2.2e-16

Como podemos ver, parece que la variable categoría es un buen predictor del puntaje de un Whisky. Acabamos de crear k-1 variables que corresponden a las diferentes categorías, en la cual todos los niveles fueron significativos.

Regresiones lineales múltiples

El modelo estadístico en regresión lineal múltiple es una generalización del regresión lineal simple para k covariables. Para la regresión lineal múltiple, a las variables predictoras las agregaremos con el comando +:

reg3 <- lm(data$Puntuacion ~ data$Categoria+ data$Votos)
summary(reg3)

## 
## Call:
## lm(formula = data$Puntuacion ~ data$Categoria + data$Votos)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -52.724  -0.943   0.047   1.103  14.194 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      9.144e+01  1.388e+00  65.870  < 2e-16 ***
## data$CategoriaC -4.850e+00  1.428e+00  -3.395 0.000765 ***
## data$CategoriaD -8.655e+00  1.419e+00  -6.101 2.83e-09 ***
## data$CategoriaE -1.336e+01  1.465e+00  -9.122  < 2e-16 ***
## data$CategoriaF -1.882e+01  1.726e+00 -10.904  < 2e-16 ***
## data$CategoriaG -3.572e+01  1.963e+00 -18.194  < 2e-16 ***
## data$Votos       6.893e-04  7.098e-04   0.971 0.332110    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.4 on 345 degrees of freedom
## Multiple R-squared:  0.6998, Adjusted R-squared:  0.6946 
## F-statistic: 134.1 on 6 and 345 DF,  p-value: < 2.2e-16

Comparación de modelos

Ahora tenemos varios modelos y queremos saber si una nueva variable mejora el modelo. Para eso usaremos la función anova() del paquete MASS:

library(MASS)
anova(reg2,reg3)

## Analysis of Variance Table
## 
## Model 1: data$Puntuacion ~ data$Categoria
## Model 2: data$Puntuacion ~ data$Categoria + data$Votos
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1    346 4000.1                           
## 2    345 3989.2  1    10.907 0.9433 0.3321

Modelos de regresión Log-lineal

Un modelo de regresión lineal hace un supuesto de modelado sobre la relación lineal entre variables. Esta suposición no es un problema cuando las variables tienen una relación lineal o casi lineal. Si los datos tienen una tendencia no lineal, un modelo de regresión lineal no puede ajustar los datos con precisión. Un enfoque para modelar datos no lineales es ajustar un modelo para una variable transformada. Una transformación de uso común es usar la función log() en las variables.

Los modelos de regresión log-lineal no son fundamentalmente diferentes de sus contrapartes lineales. Usamos el logaritmo de nuestra variable como variable predictora y luego realizamos una regresión lineal. Reutilizamos la función lm(), ya que el cálculo del coeficiente del modelo tampoco es diferente para los modelos lineales y log-lineales. La transformación se puede aplicar a la variable de destino o una o más de las variables predictoras.

reg4 <- lm(data$Puntuacion ~ data$Categoria+ log(data$Votos))
summary(reg4)

## 
## Call:
## lm(formula = data$Puntuacion ~ data$Categoria + log(data$Votos))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -52.120  -1.024   0.126   1.021  14.293 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      90.7324     1.3955  65.020  < 2e-16 ***
## data$CategoriaC  -5.5512     1.4346  -3.869  0.00013 ***
## data$CategoriaD  -9.0940     1.4127  -6.438 4.07e-10 ***
## data$CategoriaE -13.5552     1.4509  -9.343  < 2e-16 ***
## data$CategoriaF -18.8560     1.7075 -11.043  < 2e-16 ***
## data$CategoriaG -35.6127     1.9427 -18.331  < 2e-16 ***
## log(data$Votos)   0.3661     0.1262   2.901  0.00396 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.364 on 345 degrees of freedom
## Multiple R-squared:  0.7062, Adjusted R-squared:  0.7011 
## F-statistic: 138.2 on 6 and 345 DF,  p-value: < 2.2e-16

Parece que la forma de nuestros datos se vio favorecido con este cambio, ya que los problemas de la variable parece que se daban por la forma de la variables

Revisión de los supuestos

Los supuestos en un modelo de regresión se pueden escribir de la siguiente forma:

Los errores tienen distribución normal.
Los errores tienen media cero.
Los errores tiene varianza constante.
Los errores no están correlacionados.

Para determinar la media de los residuales solo necesitamos aplicar la función mean() a los residuales obtenidos en el modelo:

mean(reg4$residuals)

## [1] 8.07361e-18

Además, hay una trinidad que debemos revisar para asegurar que se cumplan las otras tres condiciones: multicolinealidad, heterocedasticidad y independencia de los residuos.

Normalidad de los residuales

Para determinar si los residuos son independientes, usaremos un gráfico QQ. Usaremos las funciones qqnorm() y qqline():

plot(reg4$residuals)

qqnorm(reg4$residuals)
qqline(reg4$residuals)

Como podemos ver, los datos parecen comportarse adecuadamente.

Independencia de los errores

Para probar si los errores se comportan normalmente podemos usar el test de Durbin Watson. La función que nos permite hacer eso es dwtest() de la librería lmtest.

library(lmtest)

## Loading required package: zoo

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

dwtest(reg4)

## 
##  Durbin-Watson test
## 
## data:  reg4
## DW = 1.9537, p-value = 0.3322
## alternative hypothesis: true autocorrelation is greater than 0

El p-valor es mayor a 0.05, por lo que no tenemos suficiente evidencia para rechazar la hipótesis nula y podemos ver que no tenemos problemas de autocorrelación de los residuos. Ahora, si la correlación es para más de un periodo podemos usar la prueba de Breusch-Godfrey con la función bgtest():

bgtest(reg4)

## 
##  Breusch-Godfrey test for serial correlation of order up to 1
## 
## data:  reg4
## LM test = 0.18592, df = 1, p-value = 0.6663

Por lo tanto, estas dos pruebas nos mostraron que no tenemos autocorrelación.

Multicolinealidad

Finalmente, tenemos que probar que nuestros datos no estén correlacionados y eliminamos la multicolinealidad. Para usaremos las correlaciones entre las variables:

reg5 <- lm(data$Puntuacion ~ data$Sabor + data$Color + data$Maltas)
cor(cbind(data$Sabor, data$Color, data$Maltas))

##              [,1]         [,2]       [,3]
## [1,] 1.0000000000 0.0006716295 0.05516003
## [2,] 0.0006716295 1.0000000000 0.03475009
## [3,] 0.0551600290 0.0347500870 1.00000000

¿Cómo más sabemos si hay problemas de multicolinealidad? Simple, si tenemos unos altos niveles de correlación entre las variables o un vif mayor a 10 (aunque por lo general se dice que 5) podemos afirmar que nuestro modelo tiene problemas. Para hallar el índice de inflación de la varianza usaremos la función vif():

library(car)

## Loading required package: carData

vif(reg5)

##  data$Sabor  data$Color data$Maltas 
##    1.003053    1.001211    1.004266

En nuestro caso podemos ver que no tenemos multicolinealidad. Pero acá no hemos acabado, todavía podemos usar una versión más avanzada.

Paquete CARET

El paquete caret (classification and regression training, Kuhn (2016)) incluye una serie de funciones que facilitan el uso de decenas de métodos complejos de clasificación y regresión. Utilizar este paquete en lugar de las funciones originales de los métodos presenta dos ventajas:

Permite utilizar un código unificado para aplicar reglas de clasificación muy distintas, implementadas en diferentes paquetes.

*Es más fácil poner en práctica algunos procedimientos usuales en problemas de clasificación. Por ejemplo, hay funciones específicas para dividir la muestra en datos de entrenamiento y datos de test o para ajustar parámetros mediante validación cruzada.

Con este paquete nos centraremos en el uso para problemas de regresión y clasificación limitadas, pero son más de 200 modelos diferentes que tendrán a su disposición.

Para ello lo primero que haremos será realizar una partición en los datos en dos grupos, train y test. Esto nos permite entrenar el modelo y, además, probar si funciona. La proporción no es exacta y su elección no suele ser aleatoria, por lo que se aconseja una proporción de 80/20. La función que usaremos es createDataPartition():

library(caret)

## Loading required package: lattice

## Loading required package: ggplot2

trainIndex <- createDataPartition(data$Puntuacion, p=0.8, list=F)

Ahora que tenemos los índices podemos crear nuestros nuevos data frames al usar este índice para hacer la división:

train <- data[trainIndex,]
test <- data[-trainIndex,]

Ahora deberemos definir el resampling para la evaluación del modelo. La finalidad última de un modelo es predecir la variable respuesta en observaciones futuras o en observaciones que el modelo no ha “visto” antes. Para conseguir una estimación más certera, y antes de recurrir al conjunto de test, se pueden emplear estrategias de validación basadas en resampling. caret incorpora los métodos boot, boot632, optimism_boot, boot_all, cv, repeatedcv, LOOCV, LGOCV. Cada uno funciona internamente de forma distinta, pero todos ellos se basan en la idea de ajustar y evaluar el modelo múltiples veces con distintos subconjuntos creados a partir de los datos de entrenamiento, obteniendo en cada repetición una estimación del error. El promedio de todas las estimaciones tiende a converger en el valor real del error de test.

Para especificar el tipo de validación, así como el número de repeticiones, se crea un control de entrenamiento mediante la función trainControl(), que se pasa al argumento trControl de la función train(). A efectos prácticos, cuando se aplican métodos de resampling hay que tener en cuenta dos cosas: el coste computacional que implica ajustar múltiples veces un modelo, cada vez con un subconjunto de datos distinto, y la reproducibilidad en la creación de las particiones. caret permite paralelizar el proceso para que sea más rápido y establecer semillas para asegurar que cada resampling o partición pueda crearse de nuevo con exactamente las mismas observaciones.

fitControl <- trainControl( method = "cv",  number = 10, savePredictions = 'final', classProbs = FALSE)

Lo siguiente será definir las variables que utilizaremos como predictoras y la variable a predecir:

predictors<-c("Sabor", "Color","Maltas")
outcomeName<-'Puntuacion'

Ahora solo nos falta correr el modelo utilizando la función train(). Lo que definirá qué modelo utilizaremos será el parámetro de method ya que la implementación es igual para todos los modelos en caret:

model1 <- train(train[,predictors], train[,outcomeName], method='lm', trControl = fitControl)
summary(model1)

## 
## Call:
## lm(formula = .outcome ~ ., data = dat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -79.618  -2.318   1.376   3.161   9.303 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 82.74896    1.50314  55.051   <2e-16 ***
## Sabor       -0.08965    0.16069  -0.558    0.577    
## Color        0.09183    0.15735   0.584    0.560    
## Maltas       0.03756    0.11471   0.327    0.744    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.397 on 280 degrees of freedom
## Multiple R-squared:  0.00262,    Adjusted R-squared:  -0.008066 
## F-statistic: 0.2452 on 3 and 280 DF,  p-value: 0.8647

Ahora, lo interesante que tiene caret es que tiene mayor variedad, incluyendo backward and forward estimation, e incluso variaciones más robustas. En este link podrán encontrar todas las posibles variaciones: https://topepo.github.io/caret/train-models-by-tag.html#linear-regression