8.6 Ejemplo correlación de Pearson

Se dispone de un data set con información sobre diferentes coches. Se quiere estudiar si existe una correlación entre el peso de un vehículo (Weight) y la potencia de su motor (Horsepower).

R contiene funciones que permiten calcular los diferentes tipos de correlaciones y sus niveles de significancia: cor() y cor.test(). La segunda función es más completa ya que además de calcular el coeficiente de correlación indica su significancia (p-value) e intervalo de confianza.

library(MASS)
library(ggplot2)
data("Cars93")

8.6.1 ¿Es lineal?

En primer lugar se representan las dos variables mediante un diagrama de dispersión para intuir si existe relación lineal o monotónica. Si no la hay, no tiene sentido calcular este tipo de correlaciones.

ggplot(data = Cars93, aes(x = Weight, y = Horsepower)) + 
  geom_point(colour = "red4") +
  ggtitle("Diagrama de dispersión") +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5))

El diagrama de dispersión parece indicar una posible relación lineal positiva entre ambas variables.

Para poder elegir el coeficiente de correlación adecuado, se tiene que analizar el tipo de variables y la distribución que presentan. En este caso, ambas variables son cuantitativas continuas y pueden transformarse en rangos para ordenarlas, por lo que a priori los tres coeficientes podrían aplicarse. La elección se hará en función de la distribución que presenten las observaciones.

8.6.2 ¿Es normal?

par(mfrow = c(1, 2))
hist(Cars93$Weight, breaks = 10, main = "", xlab = "Weight", border = "darkred")

hist(Cars93$Horsepower, breaks = 10, main = "", xlab = "Horsepower",
     border = "blue")

qqnorm(Cars93$Weight, main = "Weight", col = "darkred")
qqline(Cars93$Weight)

qqnorm(Cars93$Horsepower, main = "Horsepower", col = "blue")
qqline(Cars93$Horsepower)

par(mfrow = c(1,1))

Existe un test que nos permite analizar si la distribución de nuestra variable es normal: Shapiro. Se calcula en R muy fácilmente:

shapiro.test(Cars93$Weight)

## 
##  Shapiro-Wilk normality test
## 
## data:  Cars93$Weight
## W = 0.97432, p-value = 0.06337

Ahora lo calculamos para la segunda variable:

shapiro.test(Cars93$Horsepower)

## 
##  Shapiro-Wilk normality test
## 
## data:  Cars93$Horsepower
## W = 0.93581, p-value = 0.0001916

El análisis gráfico y el contraste de normalidad muestran que para la variable Horsepower no se puede asumir normalidad y que la variable Weight está en el límite. Siendo estrictos, este hecho excluye la posibilidad de utilizar el coeficiente de Pearson, dejando como alternativas el de Spearman o Kendall. Sin embargo, dado que la distribución no se aleja mucho de la normalidad y de que el coeficiente de Pearson tiene cierta robustez, a fines prácticos sí que se podría utilizar siempre y cuando se tenga en cuenta este hecho en los resultados. Otra posibilidad es tratar de transformar las variables para mejorar su distribución.

# Representación gráfica
par(mfrow = c(1, 2))
hist(log10(Cars93$Horsepower), breaks = 10, main = "", xlab = "Log10(Horsepower)",
     border = "blue")
     
qqnorm(log10(Cars93$Horsepower), main = "", col = "blue")
qqline(log10(Cars93$Horsepower))

par(mfrow = c(1, 1))
shapiro.test(log10(Cars93$Horsepower))

## 
##  Shapiro-Wilk normality test
## 
## data:  log10(Cars93$Horsepower)
## W = 0.98761, p-value = 0.5333

La trasformación logarítmica de la variable Horsepower consigue una distribución de tipo normal.

8.6.3 ¿Es homocedástica?

La homocedasticidad implica que la varianza se mantenga constante. Puede analizarse de forma gráfica representando las observaciones en un diagrama de dispersión y viendo si mantiene una homogeneidad en su dispersión a lo largo del eje X. Una forma cónica es un claro indicativo de homocedasticidad. En algunos libros se menciona el test de Goldfeld-Quandt o el de Breusch-Pagan como test de hipótesis para la homocedasticidad en correlación y regresión.

Tal como muestra el diagrama de dispersión generado al inicio del ejercicio, sí hay un patrón cónico. Esto debe de tenerse en cuenta si se utiliza Pearson puesto que viola una de sus condiciones.

ggplot(data = Cars93, aes(x = Weight, y = Horsepower)) + 
  geom_point(colour = "red4") +
  geom_segment(aes(x = 1690, y = 70, xend = 3100, yend = 300),linetype="dashed") +
  
    geom_segment(aes(x = 1690, y = 45, xend = 4100, yend = 100),linetype="dashed") +
  
    ggtitle("Diagrama de dispersión") +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5))

8.6.4 ¡Ahora sí! Correlación

Aplicamos una correlación de Pearson:

cor(x = Cars93$Weight, y = log10(Cars93$Horsepower), method = "pearson")

## [1] 0.809672

Para analizar el resultado de nuestra correlación de Pearson debemos tener en cuenta que:

t es el valor del estadístico t-test
df son los grados de libertad
p-value es el p valor que comprueba el nivel de significancia del t-test
conf-int es el intervalo de confianza de a un nivel de confianza de 95%
sample estimates es el coeficiente de correlación

Ahora probemos con el coeficiente de Spearman:

cor(x = Cars93$Weight, y = log10(Cars93$Horsepower), method = "spearman")

## [1] 0.8042527

Ambos test muestran una correlación alta (>0.8). Sin embargo para poder considerar que existe realmente correlación entre las dos variables es necesario calcular su significancia, de lo contrario podría deberse al azar.

8.6.5 ¿Es significativa?

Aunque nuestros coeficientes de correlación sean altos, debemos chequear que la relación sea significativa. De lo contrario, la relación se debe considerar inexistente.

Para saber si la relación analizada es significante o nó, primero analizamos el p-valor. Para ello debemos comparar el resultado del p-valor con el nivel de alfa, que es 0.05 (dado nuestro nivel de confianza de 95%).

Alfa es el grado de significancia, que es la probabilidad de rechazar nuestra hipótesis nula cuando es comprobada. Por ejemplo, el nivel de significancia de 0.05 indica que existe un 5% de riesgo de concluir que existe una diferencia cuando no hay una diferencia realmente.

En este caso, dado que el p-valor es menor al grado de significancia (alpha=0.05) podemos concluir que wt y mpg están correlacionados de manera significante con un coeficiente de -0.87 y un p-valor de 1.29410^{-10}.

Otra manera de comprobar la significancia es mediante una función aparte:

cor.test(x = Cars93$Weight,
         y = log10(Cars93$Horsepower), 
         alternative = "two.sided",
         conf.level  = 0.95,
         method      = "pearson")

## 
##  Pearson's product-moment correlation
## 
## data:  Cars93$Weight and log10(Cars93$Horsepower)
## t = 13.161, df = 91, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.7256502 0.8699014
## sample estimates:
##      cor 
## 0.809672

Interpretamos: “La relación es…”

Para extraer el p-valor directamente:

res$p.value

Y para extraer el coeficiente de correlación:

res$estimate

8.6.6 ¿Cuán “grande” es?

Ahora veremos cuál es el “tamaño” del efecto, mediante el análisis del R2 (R cuadrado).

R2_pearson <- cor(x = Cars93$Weight,
                  y = log10(Cars93$Horsepower),
                  method = "pearson")
R2_pearson <- R2_pearson^2

R2_pearson

## [1] 0.6555688

8.6.7 Interpretación

Existe evidencia de que hay una correlación significativa entre el peso del vehículo y la potencia de su motor (r=0.8, p-value < 2.2e-16), con un tamaño de efecto medio-alto (R2=0.66).

8.6.8 Conclusión

Funciones importantes:

visualizar datos ggplot()
correlación cor.test(x,y) o bien cor(x,y, method=“pearson”)
shapiro.test() o bien qqnorm()