Profesor del programa de Ingeniería Industrial, Universidad Sergio Arboleda. Candidato a doctor en Ingeniería de la Universidad de los Andes.

Regresiones lineales

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():

library(readr)
## Warning: package 'readr' was built under R version 3.6.3
data <- read_csv("Clase 2.csv")
## Warning: Missing column names filled in: 'X1' [1]
## Parsed with column specification:
## cols(
##   .default = col_double(),
##   Name = col_character(),
##   Nationality = col_character(),
##   Club = col_character(),
##   `Preferred Foot` = col_character(),
##   `Work Rate` = col_character(),
##   `Body Type` = col_character(),
##   Position = col_character(),
##   Joined = col_character(),
##   `Contract Valid Until` = col_character(),
##   Height = col_character(),
##   Weight = col_character()
## )
## See spec(...) for full column specifications.
na.omit(data)
## # A tibble: 16,654 x 49
##       X1     ID Name    Age Nationality Overall Potential Club  `Preferred Foot`
##    <dbl>  <dbl> <chr> <dbl> <chr>         <dbl>     <dbl> <chr> <chr>           
##  1     0 158023 "L. ~    31 Argentina        94        94 "FC ~ Left            
##  2     1  20801 "Cri~    33 Portugal         94        94 "Juv~ Right           
##  3     2 190871 "Ney~    26 Brazil           92        93 "Par~ Right           
##  4     3 193080 "De ~    27 Spain            91        93 "Man~ Right           
##  5     4 192985 "K. ~    27 Belgium          91        92 "Man~ Right           
##  6     5 183277 "E. ~    27 Belgium          91        91 "Che~ Right           
##  7     6 177003 "L. ~    32 Croatia          91        91 "Rea~ Right           
##  8     7 176580 "L. ~    31 Uruguay          91        91 "FC ~ Right           
##  9     8 155862 "Ser~    32 Spain            91        91 "Rea~ Right           
## 10     9 200389 "J. ~    25 Slovenia         90        93 "Atl~ Right           
## # ... with 16,644 more rows, and 40 more variables: `International
## #   Reputation` <dbl>, `Weak Foot` <dbl>, `Skill Moves` <dbl>, `Work
## #   Rate` <chr>, `Body Type` <chr>, Position <chr>, `Jersey Number` <dbl>,
## #   Joined <chr>, `Contract Valid Until` <chr>, Height <chr>, Weight <chr>,
## #   Crossing <dbl>, Finishing <dbl>, HeadingAccuracy <dbl>, ShortPassing <dbl>,
## #   Volleys <dbl>, Dribbling <dbl>, Curve <dbl>, FKAccuracy <dbl>,
## #   LongPassing <dbl>, BallControl <dbl>, Acceleration <dbl>,
## #   SprintSpeed <dbl>, Agility <dbl>, Reactions <dbl>, Balance <dbl>,
## #   ShotPower <dbl>, Jumping <dbl>, Stamina <dbl>, Strength <dbl>,
## #   LongShots <dbl>, Aggression <dbl>, Interceptions <dbl>, Positioning <dbl>,
## #   Vision <dbl>, Penalties <dbl>, Composure <dbl>, Marking <dbl>,
## #   StandingTackle <dbl>, SlidingTackle <dbl>
reg1 <- lm(data$Overall ~ data$Age)
reg1
## 
## Call:
## lm(formula = data$Overall ~ data$Age)
## 
## Coefficients:
## (Intercept)     data$Age  
##     49.4262       0.6692
summary(reg1)
## 
## Call:
## lm(formula = data$Overall ~ data$Age)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -30.872  -4.157  -0.480   3.835  25.858 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 49.426245   0.249876  197.80   <2e-16 ***
## data$Age     0.669227   0.009779   68.44   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.162 on 18205 degrees of freedom
## Multiple R-squared:  0.2046, Adjusted R-squared:  0.2046 
## F-statistic:  4683 on 1 and 18205 DF,  p-value: < 2.2e-16

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$Overall ~ data$Age)
abline(reg1)

Si queremos una regresión que pase por el origen (sin intercepto) deberemos usar la siguiente expresión:

reg2 <- lm(data$Overall ~ data$Age - 1)
summary(reg2)
## 
## Call:
## lm(formula = data$Overall ~ data$Age - 1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -65.122  -4.987   3.010   9.152  39.152 
## 
## Coefficients:
##          Estimate Std. Error t value Pr(>|t|)    
## data$Age 2.570949   0.003171   810.7   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10.93 on 18206 degrees of freedom
## Multiple R-squared:  0.973,  Adjusted R-squared:  0.973 
## F-statistic: 6.572e+05 on 1 and 18206 DF,  p-value: < 2.2e-16
plot(data$Overall ~ data$Age)
abline(reg2)

Y, ¿cómo sabemos que modelo es mejor? Es sencillo, podemos comparar el R cuadrado y el R cuadrado ajustado. Adicionalmente, podemos usar un anova entre los modelos para ver si hay diferencias significativas entre los modelos.

anova(reg1, reg2)
## Analysis of Variance Table
## 
## Model 1: data$Overall ~ data$Age
## Model 2: data$Overall ~ data$Age - 1
##   Res.Df     RSS Df Sum of Sq     F    Pr(>F)    
## 1  18205  691211                                 
## 2  18206 2176752 -1  -1485541 39126 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Para la regresión lineal múltiple, a las variables predictoras las agregaremos con el comando +:

reg3 <- lm(data$Overall ~ data$Age + data$Curve)
summary(reg3)
## 
## Call:
## lm(formula = data$Overall ~ data$Age + data$Curve)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -24.6766  -3.8795  -0.3114   3.3989  28.4727 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 44.909174   0.241432  186.01   <2e-16 ***
## data$Age     0.593966   0.009063   65.54   <2e-16 ***
## data$Curve   0.136076   0.002301   59.14   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.645 on 18156 degrees of freedom
##   (48 observations deleted due to missingness)
## Multiple R-squared:  0.3336, Adjusted R-squared:  0.3336 
## F-statistic:  4545 on 2 and 18156 DF,  p-value: < 2.2e-16
plot(reg3)

Además, podemos revisar los coeficientes y graficarlos para saber la importancia de cada uno de los predictores. Para esto utilizaremos la funcion coefplot():

library(coefplot)
## Warning: package 'coefplot' was built under R version 3.6.3
## Loading required package: ggplot2
## Warning: package 'ggplot2' was built under R version 3.6.3
coefplot(reg3 , intercept=FALSE, lwdOuter =1.5,)

Sin embargo, hay una trinidad que debemos revisar para asegurar que se cumplan tres condiciones: multicolinealidad, heterocedasticidad y independencia de los residuos.

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

plot(reg3$residuals)

qqnorm(reg3$residuals)
qqline(reg3$residuals)

Como podemos ver, los datos no parecen comportarse adecuadamente. Para esto veremos varias variaciones que nos pueden solucionar el problema más adelante.

Además, podemos usar el test de Durbin Watson para saber si se comportan normalmente. La función que nos permite hacer eso es durbinWatsonTest() de la librería car. La hipótesis nula de esta prueba es que

library(car)
## Warning: package 'car' was built under R version 3.6.3
## Loading required package: carData
## Warning: package 'carData' was built under R version 3.6.3
durbinWatsonTest(reg3)
##  lag Autocorrelation D-W Statistic p-value
##    1       0.6838447     0.6314172       0
##  Alternative hypothesis: rho != 0

Este valor tiene un valor de y podemos ver que no tenemos problemas de autocorrelación de los residuos.

Finalmente, tenemos que probar que nuestros datos no estén correlacionados y eliminamos la multicolinealidad. Para usaremos las correlaciones y el índice de inflación de la varianza.

cor(cbind(data$Overall, data$Finishing, data$Curve))
##      [,1] [,2] [,3]
## [1,]    1   NA   NA
## [2,]   NA    1   NA
## [3,]   NA   NA    1
vif(reg3)
##   data$Age data$Curve 
##   1.020958   1.020958

¿Cómo 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. En nuestro caso podemos ver que no tenemos multicolinealidad.

Entonces, ¿cuál es el problema? Parece que la forma funcional de los datos no nos ayuda. Para eso, la solución más fácil es linealizar las variables. Para eso utilizaremos la función log().

reg4 <- lm(data$Overall ~ log(data$Finishing) + log(data$Curve))
reg4
## 
## Call:
## lm(formula = data$Overall ~ log(data$Finishing) + log(data$Curve))
## 
## Coefficients:
##         (Intercept)  log(data$Finishing)      log(data$Curve)  
##             49.0287              -0.1434               4.7333
summary(reg4)
## 
## Call:
## lm(formula = data$Overall ~ log(data$Finishing) + log(data$Curve))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -19.9838  -4.2207  -0.1949   4.1577  29.1744 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          49.0287     0.3642  134.61   <2e-16 ***
## log(data$Finishing)  -0.1434     0.1542   -0.93    0.353    
## log(data$Curve)       4.7333     0.1800   26.29   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.518 on 18156 degrees of freedom
##   (48 observations deleted due to missingness)
## Multiple R-squared:  0.1114, Adjusted R-squared:  0.1113 
## F-statistic:  1138 on 2 and 18156 DF,  p-value: < 2.2e-16
plot(reg4)

¿Que pasó? Hay que explorar la forma funcional de las variables para solucionar el problema.

plot(data$Overall ~ data$Finishing)

plot(data$Overall ~ log(data$Finishing))

El verdadero problema es que esta variable no sirve para explicar la otra. En este caso el problema ya no es estadístico, es de especificación de las variables. Ahí es donde el conocimiento específico de su área puede ser muy útil.

Por eso, hemos decidido incluir una variable international reputation que es un factor y eliminar la variable curve:

data$`International Reputation` <- as.factor(data$`International Reputation`)
reg5 <- lm(data$Overall ~ data$Finishing + data$`International Reputation` )
reg5
## 
## Call:
## lm(formula = data$Overall ~ data$Finishing + data$`International Reputation`)
## 
## Coefficients:
##                      (Intercept)                    data$Finishing  
##                         61.26025                           0.08828  
## data$`International Reputation`2  data$`International Reputation`3  
##                          9.45386                          14.71772  
## data$`International Reputation`4  data$`International Reputation`5  
##                         19.22243                          22.68725
summary(reg5)
## 
## Call:
## lm(formula = data$Overall ~ data$Finishing + data$`International Reputation`)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -21.1281  -3.6453   0.2101   3.9133  20.7687 
## 
## Coefficients:
##                                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                      61.26025    0.10789 567.807   <2e-16 ***
## data$Finishing                    0.08828    0.00221  39.949   <2e-16 ***
## data$`International Reputation`2  9.45386    0.16866  56.053   <2e-16 ***
## data$`International Reputation`3 14.71772    0.33029  44.560   <2e-16 ***
## data$`International Reputation`4 19.22243    0.80328  23.930   <2e-16 ***
## data$`International Reputation`5 22.68725    2.33622   9.711   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.719 on 18153 degrees of freedom
##   (48 observations deleted due to missingness)
## Multiple R-squared:  0.3162, Adjusted R-squared:  0.316 
## F-statistic:  1679 on 5 and 18153 DF,  p-value: < 2.2e-16

Como pueden notar en el resultado, este modelo hace más cosas que los anteriores al incluir cada una de las categorías por separado, las cuales parecen ser significativas. En este caso solo restaría hacer de nuevo las pruebas y revisar que se cumplan los supuestos y ya, tendremos un modelo completamente funcional.

Regresiones logísticas

Ahora, ,miramemos un modelo más de regresión. la regresión logística es un tipo de análisis de regresión utilizado para predecir el resultado de una variable categórica en función de las variables independientes o predictoras. Es útil para modelar la probabilidad de un evento ocurriendo como función de otros factores.

La característica central de un modelo logístico es la siguiente: relaciona la probabilidad de un resultado con una función exponencial de una variable predictora. Al modelar la probabilidad de un resultado, un modelo logístico logra dos cosas. Primero, modela más directamente lo que nos interesa, que es una probabilidad o proporción, como la probabilidad de que un cliente determinado compre un producto o la proporción esperada de un segmento que responderá a una promoción. En segundo lugar, limita el modelo al rango apropiado para una proporción, que es [0, 1]. Un modelo lineal básico como se genera con lm () no tiene tal límite y podría estimar una probabilidad sin sentido como 1.05 o -0.04. La ecuación de la función logística es: Ecuación 1: Ecuación Logística Veamos un ejemplo de lo que hace la función aplicando la probabilidad logística con plogis:

plogis( exp (0) / (exp (0) + 1))
## [1] 0.6224593
plogis(-Inf)
## [1] 0
plogis ( -0.2)
## [1] 0.450166

El análisis de regresión logística se enmarca en el conjunto de Modelos Lineales Generalizados (GLM por sus siglas en inglés). Las probabilidades que describen el posible resultado de un único ensayo se modelan, como una función de variables explicativas, utilizando una función logística. Por lo tanto, los modelos lineales generalizados se pueden utilizar para modelar recuentos de datos (como el número de compras) o intervalos de tiempo (como el tiempo que se pasa en un sitio web) o variables binarias (por ejemplo, compró / no compró). La característica común de todos los modelos GLM es que relacionan predictores distribuidos normalmente con un resultado no normal mediante una función conocida como enlace. Esto significa que pueden adaptarse a modelos para muchas distribuciones diferentes utilizando un marco único y coherente.

Para construir nuestra función, tenemos que asegurarnos que la variable que vamos a predecir esté compuesta por factores. En este caso, lo haremos con la variable bebedor (o Drinker). Además, como haremos uso de modelos lineales generalizados (como vimos arriba), la función a utilizar es glm().

data1 <- read_csv("Clase 3.csv")
## Parsed with column specification:
## cols(
##   ID = col_double(),
##   City = col_character(),
##   Gender = col_character(),
##   Age = col_double(),
##   Income = col_double(),
##   Illness = col_character(),
##   Height1 = col_double(),
##   Height2 = col_double(),
##   Height3 = col_double(),
##   Weight1 = col_double(),
##   Weight2 = col_double(),
##   Weight3 = col_double(),
##   Smoker = col_double(),
##   Drinker = col_double(),
##   Cocaine = col_double(),
##   ExerciseTime = col_double()
## )
table(data1$Drinker)
## 
##     0     1 
## 90000 60000
rl1 <- glm(data1$Drinker ~ data1$Height1 + data1$Weight1, family = binomial)
rl1
## 
## Call:  glm(formula = data1$Drinker ~ data1$Height1 + data1$Weight1, 
##     family = binomial)
## 
## Coefficients:
##   (Intercept)  data1$Height1  data1$Weight1  
##    -3.630e-01     -2.857e-04      9.884e-05  
## 
## Degrees of Freedom: 149999 Total (i.e. Null);  149997 Residual
## Null Deviance:       201900 
## Residual Deviance: 201900    AIC: 201900
summary(rl1)
## 
## Call:
## glm(formula = data1$Drinker ~ data1$Height1 + data1$Weight1, 
##     family = binomial)
## 
## Deviance Residuals: 
##    Min      1Q  Median      3Q     Max  
## -1.023  -1.011  -1.008   1.353   1.368  
## 
## Coefficients:
##                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   -3.630e-01  4.268e-02  -8.505   <2e-16 ***
## data1$Height1 -2.857e-04  2.105e-04  -1.358    0.175    
## data1$Weight1  9.884e-05  2.115e-04   0.467    0.640    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 201904  on 149999  degrees of freedom
## Residual deviance: 201901  on 149997  degrees of freedom
## AIC: 201907
## 
## Number of Fisher Scoring iterations: 4

¿Cómo sabemos si le pegamos a los resultados? Creamos una matriz de confusión con los datos. Esta no es más que una matriz de contigencia para modelos de clasificación.

predicciones <- ifelse(test = rl1$fitted.values > 0.4, yes = 1, no = 0)
t1 <- table(data1$Drinker, predicciones)

Esto es solamente un primer acercamiento a los modelos lineales. Queda mucho por explorar y por conocer ¿Hay algún otro tema que sea de su interés? Los invito a mirar las librerias y paquetes que hay en R para resolver muchos de los problemas que tenemos en nuestro campo profesional.