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

Coeficientes de correlación

Así como vimos las pruebas de hipótesis para variables continuas y discretas, es necesario conocer el proceso para realizar pruebas de hipótesis entre variables continuas. Aunque los diagramas de dispersión proporcionan mucha información visual, cuando hay más de unas pocas variables, puede ser útil evaluar la relación entre cada par con un solo número.Una medida de la relación entre dos variables es la covarianza, que se puede calcular para dos variables cualesquiera usando la función cov():

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.
data <- na.omit(data)

cov(data$Finishing, data$Dribbling)
## [1] 306.9026

Si los valores x e y tienden a ir en la misma dirección (ambos más altos o ambos más bajos que sus respectivas medias) entre las observaciones, entonces tienen una covarianza positiva. Si cov (x, y) es cero, entonces no hay asociación (lineal) entre x e y. La covarianza negativa significa que las variables van en direcciones opuestas con respecto a sus medias: cuando x es menor, y tiende a ser mayor. El problema con esta medida es igual al de la varianza, no se encuentra en una medida adecuada para ser analizado. Por tanto, usaremos la función cor() para revisar las correlaciones de las variables.

cor(data$Finishing, data$Dribbling)
## [1] 0.8260533

¿Y cómo sabemos si es pequeña o grande una correlación? Sencillo, se llama débil entre 0.1 y 0.3, mediana entre 0.3 y 0.5 y fuerte por encima de 0.5. Debemos recordar que esto es una prueba de hipótesis, entonces debemos probar esta relación. Para ello, utilizaremos la función cor.test(), asumiendo que la hipótesis nula es que la correlación entre las variables es igual a cero.

cor.test(data$Finishing, data$Dribbling)
## 
##  Pearson's product-moment correlation
## 
## data:  data$Finishing and data$Dribbling
## t = 189.14, df = 16652, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.8211678 0.8308177
## sample estimates:
##       cor 
## 0.8260533

Modelos 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 esta sección veremos tres versiones de la regresión lineal

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

reg1 <- lm(data$Overall ~ data$Age)
reg1
## 
## Call:
## lm(formula = data$Overall ~ data$Age)
## 
## Coefficients:
## (Intercept)     data$Age  
##     48.9065       0.6837
summary(reg1)
## 
## Call:
## lm(formula = data$Overall ~ data$Age)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -30.9878  -4.2125  -0.4165   3.6855  26.1039 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 48.90651    0.26152  187.01   <2e-16 ***
## data$Age     0.68367    0.01019   67.12   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.22 on 16652 degrees of freedom
## Multiple R-squared:  0.213,  Adjusted R-squared:  0.2129 
## F-statistic:  4506 on 1 and 16652 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 
## -64.454  -5.006   2.994   9.108  39.440 
## 
## Coefficients:
##          Estimate Std. Error t value Pr(>|t|)    
## data$Age 2.555770   0.003305   773.3   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10.95 on 16653 degrees of freedom
## Multiple R-squared:  0.9729, Adjusted R-squared:  0.9729 
## F-statistic: 5.98e+05 on 1 and 16653 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  16652  644198                                 
## 2  16653 1997156 -1  -1352958 34973 < 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.6975  -3.9113  -0.3674   3.3880  28.6258 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 44.520190   0.251860  176.77   <2e-16 ***
## data$Age     0.603149   0.009439   63.90   <2e-16 ***
## data$Curve   0.136563   0.002418   56.48   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.698 on 16651 degrees of freedom
## Multiple R-squared:  0.3395, Adjusted R-squared:  0.3394 
## F-statistic:  4279 on 2 and 16651 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.6809918     0.6370678       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.0000000 0.3349401 0.4213158
## [2,] 0.3349401 1.0000000 0.7622167
## [3,] 0.4213158 0.7622167 1.0000000
vif(reg3)
##   data$Age data$Curve 
##   1.023346   1.023346

¿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)  
##             48.7676              -0.1416               4.7836
summary(reg4)
## 
## Call:
## lm(formula = data$Overall ~ log(data$Finishing) + log(data$Curve))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -19.9146  -4.2943  -0.2558   4.1786  29.3021 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          48.7676     0.3819 127.709   <2e-16 ***
## log(data$Finishing)  -0.1416     0.1630  -0.869    0.385    
## log(data$Curve)       4.7836     0.1903  25.136   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.605 on 16651 degrees of freedom
## Multiple R-squared:  0.1126, Adjusted R-squared:  0.1125 
## F-statistic:  1056 on 2 and 16651 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.10214                           0.08983  
## data$`International Reputation`2  data$`International Reputation`3  
##                          9.54532                          15.00010  
## data$`International Reputation`4  data$`International Reputation`5  
##                         19.29001                          22.72475
summary(reg5)
## 
## Call:
## lm(formula = data$Overall ~ data$Finishing + data$`International Reputation`)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -21.139  -3.671   0.215   3.934  20.910 
## 
## Coefficients:
##                                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                      61.102137   0.113157 539.974   <2e-16 ***
## data$Finishing                    0.089826   0.002331  38.534   <2e-16 ***
## data$`International Reputation`2  9.545319   0.177962  53.637   <2e-16 ***
## data$`International Reputation`3 15.000098   0.342417  43.806   <2e-16 ***
## data$`International Reputation`4 19.290007   0.828244  23.290   <2e-16 ***
## data$`International Reputation`5 22.724749   2.361179   9.624   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.779 on 16648 degrees of freedom
## Multiple R-squared:  0.3206, Adjusted R-squared:  0.3204 
## F-statistic:  1572 on 5 and 16648 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.