library(tidyverse)
library(openintro)
#install.packages("GGally")
library(GGally)
library(knitr)
library(kableExtra)
options(knitr.table.format = "html") 

En estas notas de clase veremos el concepto de correlación:

\[\rho_{x,y}=\frac{cov(x,y)}{\sigma_x \sigma_y}\]

mtcars

primero, veamos de qué se trata el dataset. Para esto, hacemos un head() de la tabla. Aprovechamos para usar la librería knitr, cuya función kable() permite realizar mejores presentaciones de resultados. Con kable_styling() podemos modificar algunas características de la tabla

mtcars %>% 
  head() %>% 
  kable() %>% 
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))
mpg cyl disp hp drat wt qsec vs am gear carb
Mazda RX4 21.0 6 160 110 3.90 2.620 16.46 0 1 4 4
Mazda RX4 Wag 21.0 6 160 110 3.90 2.875 17.02 0 1 4 4
Datsun 710 22.8 4 108 93 3.85 2.320 18.61 1 1 4 1
Hornet 4 Drive 21.4 6 258 110 3.08 3.215 19.44 1 0 3 1
Hornet Sportabout 18.7 8 360 175 3.15 3.440 17.02 0 0 3 2
Valiant 18.1 6 225 105 2.76 3.460 20.22 1 0 3 1

Con ggpairs(), podemos graficar todas las variables, y buscar las correlaciones

Coloreamos por:

-\(am\): Tipo de transmisión: automatico (am=0) o manual (am=1)

mtcars %>% 
  select(-carb,-vs) %>% 
  mutate(cyl = factor(cyl),
         am = factor(am)) %>% 
ggpairs(., 
        title = "Matriz de correlaciones",
        mapping = aes(colour= am))

Veamos la correlación entre:

Miramos el scatter plot y pareciera haber una relación negativa.

La mitad superior de la matriz muestra la estimación puntual de la correlación, para todos los datos y considerando cada conjunto por separado. Recordemos que la fórmula para calcular ese estimador es:

\[ r = \frac{\sum_{i=1}^n(x_i-\bar{x})(y_i-\bar{y})}{\sqrt{\sum_{i=1}^n(x_i-\bar{x})^2} \sqrt{\sum_{i=1}^n(y_i-\bar{y})^2}} \]

Si quisieramos testear la significatividad de este estimador:

\(H_0\) : ρ =0
\(H_1\) : ρ \(\neq\) 0

cor.test(mtcars$mpg,mtcars$hp)

    Pearson's product-moment correlation

data:  mtcars$mpg and mtcars$hp
t = -6.7424, df = 30, p-value = 1.788e-07
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.8852686 -0.5860994
sample estimates:
       cor 
-0.7761684 

Con este p-value rechazamos \(H_0\)

¿y si queremos comparar la relación entre \(drat\) y \(gear\)?

Con ggpairs() ya habíamos visto que la relación era diferente entre los automáticos y con transmisión manual. Sabiendo esto, volvamos a calcular los estimadores puntuales de cada grupo

mtcars %>% 
  group_by(am) %>% 
  summarise(cor = cor(drat, gear))

La correlación para los autos automáticos da súper alto! quedemosnos con ese grupo

mtcars2 <- mtcars %>% filter(am==0)
ggplot(mtcars2, aes(gear,drat, group=gear, fill = factor(gear)))+
  geom_boxplot(alpha= 0.75)

No parece muy correcto hacer un test de correlación de pearson, es decir buscar una relación lineal, con una variable que sólo toma dos valores.

Usemos el test de correlación de Spearman

cor.test(mtcars2$gear,mtcars2$drat, method = "pearson")

    Pearson's product-moment correlation

data:  mtcars2$gear and mtcars2$drat
t = 5.1262, df = 17, p-value = 8.421e-05
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.5030694 0.9110028
sample estimates:
      cor 
0.7792264 
cor.test(mtcars2$gear,mtcars2$drat, method = "spearman")
Cannot compute exact p-value with ties

    Spearman's rank correlation rho

data:  mtcars2$gear and mtcars2$drat
S = 383.98, p-value = 0.001968
alternative hypothesis: true rho is not equal to 0
sample estimates:
      rho 
0.6631736 

Noten que el test de Spearman ya no da tan significativo como el de Pearson

LS0tCnRpdGxlOiAiQ29ycmVsYWNpb24iCm91dHB1dDoKICBodG1sX25vdGVib29rOgogICAgdGhlbWU6IHNwYWNlbGFiCiAgICBkZl9wcmludDogcGFnZWQKLS0tCgpgYGB7ciBtZXNzYWdlPUZBTFNFfQpsaWJyYXJ5KHRpZHl2ZXJzZSkKbGlicmFyeShvcGVuaW50cm8pCiNpbnN0YWxsLnBhY2thZ2VzKCJHR2FsbHkiKQpsaWJyYXJ5KEdHYWxseSkKbGlicmFyeShrbml0cikKbGlicmFyeShrYWJsZUV4dHJhKQpvcHRpb25zKGtuaXRyLnRhYmxlLmZvcm1hdCA9ICJodG1sIikgCmBgYAoKCkVuIGVzdGFzIG5vdGFzIGRlIGNsYXNlIHZlcmVtb3MgZWwgY29uY2VwdG8gZGUgY29ycmVsYWNpw7NuOiAKCiQkXHJob197eCx5fT1cZnJhY3tjb3YoeCx5KX17XHNpZ21hX3ggXHNpZ21hX3l9JCQKCiMjIG10Y2FycwoKcHJpbWVybywgdmVhbW9zIGRlIHF1w6kgc2UgdHJhdGEgZWwgZGF0YXNldC4gUGFyYSBlc3RvLCBoYWNlbW9zIHVuIGBoZWFkKClgIGRlIGxhIHRhYmxhLiBBcHJvdmVjaGFtb3MgcGFyYSB1c2FyIGxhIGxpYnJlcsOtYSBga25pdHJgLCBjdXlhIGZ1bmNpw7NuIGBrYWJsZSgpYCBwZXJtaXRlIHJlYWxpemFyIG1lam9yZXMgcHJlc2VudGFjaW9uZXMgZGUgcmVzdWx0YWRvcy4gQ29uIGBrYWJsZV9zdHlsaW5nKClgIHBvZGVtb3MgbW9kaWZpY2FyIGFsZ3VuYXMgY2FyYWN0ZXLDrXN0aWNhcyBkZSBsYSB0YWJsYQoKCmBgYHtyfQptdGNhcnMgJT4lIAogIGhlYWQoKSAlPiUgCiAga2FibGUoKSAlPiUgCiAga2FibGVfc3R5bGluZyhib290c3RyYXBfb3B0aW9ucyA9IGMoInN0cmlwZWQiLCAiaG92ZXIiLCAiY29uZGVuc2VkIiwgInJlc3BvbnNpdmUiKSkKYGBgCgoKQ29uIGBnZ3BhaXJzKClgLCBwb2RlbW9zIGdyYWZpY2FyIHRvZGFzIGxhcyB2YXJpYWJsZXMsIHkgYnVzY2FyIGxhcyBjb3JyZWxhY2lvbmVzCgpDb2xvcmVhbW9zIHBvcjoKCi0kYW0kOiBUaXBvIGRlIHRyYW5zbWlzacOzbjogYXV0b21hdGljbyAoYW09MCkgbyBtYW51YWwgKGFtPTEpCgoKYGBge3IgbWVzc2FnZT1GQUxTRSwgd2FybmluZz1GQUxTRSwgZmlnLndpZHRoPTEwLCBmaWcuaGVpZ2h0PTEwfQptdGNhcnMgJT4lIAogIHNlbGVjdCgtY2FyYiwtdnMpICU+JSAKICBtdXRhdGUoY3lsID0gZmFjdG9yKGN5bCksCiAgICAgICAgIGFtID0gZmFjdG9yKGFtKSkgJT4lIApnZ3BhaXJzKC4sIAogICAgICAgIHRpdGxlID0gIk1hdHJpeiBkZSBjb3JyZWxhY2lvbmVzIiwKICAgICAgICBtYXBwaW5nID0gYWVzKGNvbG91cj0gYW0pKQpgYGAKCgpWZWFtb3MgbGEgY29ycmVsYWNpw7NuIGVudHJlOgoKLSAkbXBnJDogTWlsZXMvKFVTKSBnYWxsb24uIEVmaWNpZW5jaWEgZGUgY29tYnVzdGlibGUKLSAkaHAkOiBHcm9zcyBob3JzZXBvd2VyOiBQb3RlbmNpYSBkZWwgbW90b3IKCk1pcmFtb3MgZWwgc2NhdHRlciBwbG90IHkgcGFyZWNpZXJhIGhhYmVyIHVuYSByZWxhY2nDs24gbmVnYXRpdmEuIAoKTGEgbWl0YWQgc3VwZXJpb3IgZGUgbGEgbWF0cml6IG11ZXN0cmEgbGEgZXN0aW1hY2nDs24gcHVudHVhbCBkZSBsYSBjb3JyZWxhY2nDs24sIHBhcmEgdG9kb3MgbG9zIGRhdG9zIHkgY29uc2lkZXJhbmRvIGNhZGEgY29uanVudG8gcG9yIHNlcGFyYWRvLiBSZWNvcmRlbW9zIHF1ZSBsYSBmw7NybXVsYSBwYXJhIGNhbGN1bGFyIGVzZSBlc3RpbWFkb3IgZXM6CgoKJCQKciA9IFxmcmFje1xzdW1fe2k9MX1ebih4X2ktXGJhcnt4fSkoeV9pLVxiYXJ7eX0pfXtcc3FydHtcc3VtX3tpPTF9Xm4oeF9pLVxiYXJ7eH0pXjJ9IFxzcXJ0e1xzdW1fe2k9MX1ebih5X2ktXGJhcnt5fSleMn19CiQkCgoKClNpIHF1aXNpZXJhbW9zIHRlc3RlYXIgbGEgc2lnbmlmaWNhdGl2aWRhZCBkZSBlc3RlIGVzdGltYWRvcjoKCiRIXzAkIDogz4EgPTAgICAgICAgIAokSF8xJCA6IM+BICRcbmVxJCAwICAgICAgCgpgYGB7cn0KY29yLnRlc3QobXRjYXJzJG1wZyxtdGNhcnMkaHApCmBgYAoKQ29uIGVzdGUgcC12YWx1ZSByZWNoYXphbW9zICRIXzAkCgoKwr95IHNpIHF1ZXJlbW9zIGNvbXBhcmFyIGxhIHJlbGFjacOzbiBlbnRyZSAkZHJhdCQgeSAkZ2VhciQ/CgotICRkcmF0JDogIGxhIHJlbGFjacOzbiBkZSBlbmdyYW5hamUgZGVsIGVqZSB0cmFzZXJvOiBpbmRpY2EgZWwgbsO6bWVybyBkZSB2dWVsdGFzIGRlbCBlamUgZGUgdHJhbnNtaXNpw7NuIHBhcmEgY2FkYSByb3RhY2nDs24gZGVsIGVqZSBkZSBsYSBydWVkYS4gVW4gdmVow61jdWxvIGNvbiB1bmEgcmVsYWNpw7NuIGFsdGEgcHJvcG9yY2lvbmFyw61hIG3DoXMgcGFyIHksIHBvciBsbyB0YW50bywgbcOhcyBjYXBhY2lkYWQgZGUgcmVtb2xxdWUsIHBvciBlamVtcGxvCi0gJGdlYXIkOiBOw7ptZXJvIGRlIHZlbG9jaWRhZGVzIGhhY2lhIGFkZWxhbnRlCgpDb24gYGdncGFpcnMoKWAgeWEgaGFiw61hbW9zIHZpc3RvIHF1ZSBsYSByZWxhY2nDs24gZXJhIGRpZmVyZW50ZSBlbnRyZSBsb3MgYXV0b23DoXRpY29zIHkgY29uIHRyYW5zbWlzacOzbiBtYW51YWwuIFNhYmllbmRvIGVzdG8sIHZvbHZhbW9zIGEgY2FsY3VsYXIgbG9zIGVzdGltYWRvcmVzIHB1bnR1YWxlcyBkZSBjYWRhIGdydXBvIAoKYGBge3J9Cm10Y2FycyAlPiUgCiAgZ3JvdXBfYnkoYW0pICU+JSAKICBzdW1tYXJpc2UoY29yID0gY29yKGRyYXQsIGdlYXIpKQpgYGAKCkxhIGNvcnJlbGFjacOzbiBwYXJhIGxvcyBhdXRvcyBhdXRvbcOhdGljb3MgZGEgc8O6cGVyIGFsdG8hIHF1ZWRlbW9zbm9zIGNvbiBlc2UgZ3J1cG8KCmBgYHtyfQptdGNhcnMyIDwtIG10Y2FycyAlPiUgZmlsdGVyKGFtPT0wKQpnZ3Bsb3QobXRjYXJzMiwgYWVzKGdlYXIsZHJhdCwgZ3JvdXA9Z2VhciwgZmlsbCA9IGZhY3RvcihnZWFyKSkpKwogIGdlb21fYm94cGxvdChhbHBoYT0gMC43NSkKYGBgCgoKTm8gcGFyZWNlIG11eSBjb3JyZWN0byBoYWNlciB1biB0ZXN0IGRlIGNvcnJlbGFjacOzbiBkZSBwZWFyc29uLCBlcyBkZWNpciBidXNjYXIgdW5hIHJlbGFjacOzbiBsaW5lYWwsIGNvbiB1bmEgdmFyaWFibGUgcXVlIHPDs2xvIHRvbWEgZG9zIHZhbG9yZXMuCgpVc2Vtb3MgZWwgdGVzdCBkZSBjb3JyZWxhY2nDs24gZGUgU3BlYXJtYW4KCmBgYHtyfQoKY29yLnRlc3QobXRjYXJzMiRnZWFyLG10Y2FyczIkZHJhdCwgbWV0aG9kID0gInBlYXJzb24iKQpjb3IudGVzdChtdGNhcnMyJGdlYXIsbXRjYXJzMiRkcmF0LCBtZXRob2QgPSAic3BlYXJtYW4iKQpgYGAKCgpOb3RlbiBxdWUgZWwgdGVzdCBkZSBTcGVhcm1hbiB5YSBubyBkYSB0YW4gc2lnbmlmaWNhdGl2byBjb21vIGVsIGRlIFBlYXJzb24=