El material de este taller ha sido creado por R-luminatis para estudiantes de biología de la Universidad Nacional de Costa Rica (e interesados en aprender).

Si usted quiere ser partícipe de los talleres, o quiere organizar este mismo taller en su centro de estudio o grupo de trabajo, puede contactarnos a los correos:

Contacto: ronny.hernandez@gmail.com / alvarovh95@gmail.com

Taller de R Correlación

Correlación no implica causalidad.

¿Para qué usamos correlaciones?

-Se usa para indicar la relación entre variables cuantitativas.

-La magnitud va de 0 a 1 siendo 0 ninguna relación y 1 una relación perfectamente predictiva.

-Además del coeficiente existe la prueba de significancia (p-value).

Tipos de correlación:

  • Pearson: evalúa el grado de relación entre dos variables cuantitativas.
  • Spearman: evalúa el grado de relación entre dos variables de rango ordenado
  • Kendall: es una medida no paramétrica de un rango de correlación.

Supuestos:

  • Variables continuas
  • Distribución normal
  • Homocedasticidad

Paquetes a utilizar en la sesión:

install.packages(c("psych","ggplot2","car","Hmisc","corrplot"))

Para cargar los paquetes:

library("psych")
library("ggplot2")
library("car")
library("Hmisc")
library("corrplot")

Inicio de sesión:

Revisar directorio de trabajo

getwd() #Con este vemos el directorio de trabajo que tenemos

#Para cambarlo:
setwd("C:/Users/Mechas/Desktop/Taller_R")

Vamos a utilizar un set de datos que vienen incluidos en el paquete “car”:

head(Soils)
##   Group Contour Depth Gp Block   pH     N Dens   P    Ca   Mg    K   Na
## 1     1     Top  0-10 T0     1 5.40 0.188 0.92 215 16.35 7.65 0.72 1.14
## 2     1     Top  0-10 T0     2 5.65 0.165 1.04 208 12.25 5.15 0.71 0.94
## 3     1     Top  0-10 T0     3 5.14 0.260 0.95 300 13.02 5.68 0.68 0.60
## 4     1     Top  0-10 T0     4 5.14 0.169 1.10 248 11.92 7.88 1.09 1.01
## 5     2     Top 10-30 T1     1 5.14 0.164 1.12 174 14.17 8.12 0.70 2.17
## 6     2     Top 10-30 T1     2 5.10 0.094 1.22 129  8.55 6.92 0.81 2.67
##   Conduc
## 1   1.09
## 2   1.35
## 3   1.41
## 4   1.64
## 5   1.85
## 6   3.18
str(Soils)
## 'data.frame':    48 obs. of  14 variables:
##  $ Group  : Factor w/ 12 levels "1","2","3","4",..: 1 1 1 1 2 2 2 2 3 3 ...
##  $ Contour: Factor w/ 3 levels "Depression","Slope",..: 3 3 3 3 3 3 3 3 3 3 ...
##  $ Depth  : Factor w/ 4 levels "0-10","10-30",..: 1 1 1 1 2 2 2 2 3 3 ...
##  $ Gp     : Factor w/ 12 levels "D0","D1","D3",..: 9 9 9 9 10 10 10 10 11 11 ...
##  $ Block  : Factor w/ 4 levels "1","2","3","4": 1 2 3 4 1 2 3 4 1 2 ...
##  $ pH     : num  5.4 5.65 5.14 5.14 5.14 5.1 4.7 4.46 4.37 4.39 ...
##  $ N      : num  0.188 0.165 0.26 0.169 0.164 0.094 0.1 0.112 0.112 0.058 ...
##  $ Dens   : num  0.92 1.04 0.95 1.1 1.12 1.22 1.52 1.47 1.07 1.54 ...
##  $ P      : int  215 208 300 248 174 129 117 170 121 115 ...
##  $ Ca     : num  16.4 12.2 13 11.9 14.2 ...
##  $ Mg     : num  7.65 5.15 5.68 7.88 8.12 ...
##  $ K      : num  0.72 0.71 0.68 1.09 0.7 0.81 0.39 0.7 0.74 0.77 ...
##  $ Na     : num  1.14 0.94 0.6 1.01 2.17 2.67 3.32 3.76 5.74 5.85 ...
##  $ Conduc : num  1.09 1.35 1.41 1.64 1.85 3.18 4.16 5.14 5.73 6.45 ...

Si quieren saber de otros data sets para jugar R, pueden buscar acá.

Ya que hemos visto algunas características del data set, podríamos hacer algunas relaciones exploratorias:

cor(Soils[,6:14])
##                pH          N       Dens          P         Ca         Mg
## pH      1.0000000  0.6366540 -0.5890264  0.5910303  0.8086293 -0.3957821
## N       0.6366540  1.0000000 -0.8641559  0.8422007  0.8502155 -0.5215444
## Dens   -0.5890264 -0.8641559  1.0000000 -0.7936652 -0.7914376  0.4901171
## P       0.5910303  0.8422007 -0.7936652  1.0000000  0.6875874 -0.4889733
## Ca      0.8086293  0.8502155 -0.7914376  0.6875874  1.0000000 -0.4274958
## Mg     -0.3957821 -0.5215444  0.4901171 -0.4889733 -0.4274958  1.0000000
## K       0.5795727  0.6760033 -0.6670677  0.5557269  0.7209104 -0.3567182
## Na     -0.6932614 -0.8119353  0.7423018 -0.7728571 -0.7889082  0.5645363
## Conduc -0.7648104 -0.8037846  0.7625652 -0.7616939 -0.8320952  0.5082623
##                 K         Na     Conduc
## pH      0.5795727 -0.6932614 -0.7648104
## N       0.6760033 -0.8119353 -0.8037846
## Dens   -0.6670677  0.7423018  0.7625652
## P       0.5557269 -0.7728571 -0.7616939
## Ca      0.7209104 -0.7889082 -0.8320952
## Mg     -0.3567182  0.5645363  0.5082623
## K       1.0000000 -0.6932082 -0.7531033
## Na     -0.6932082  1.0000000  0.9724087
## Conduc -0.7531033  0.9724087  1.0000000

Esto nos regresa una matrix de correlación, donde se muestra el coeficiente de correlación de Pearson (que viene por default) entre todas las variables numéricas continuas.

Hay una manera de que R nos redondee los resultados, para así tener una mejor lectura de los coeficientes:

round(cor(Soils[,6:14]),2)
##           pH     N  Dens     P    Ca    Mg     K    Na Conduc
## pH      1.00  0.64 -0.59  0.59  0.81 -0.40  0.58 -0.69  -0.76
## N       0.64  1.00 -0.86  0.84  0.85 -0.52  0.68 -0.81  -0.80
## Dens   -0.59 -0.86  1.00 -0.79 -0.79  0.49 -0.67  0.74   0.76
## P       0.59  0.84 -0.79  1.00  0.69 -0.49  0.56 -0.77  -0.76
## Ca      0.81  0.85 -0.79  0.69  1.00 -0.43  0.72 -0.79  -0.83
## Mg     -0.40 -0.52  0.49 -0.49 -0.43  1.00 -0.36  0.56   0.51
## K       0.58  0.68 -0.67  0.56  0.72 -0.36  1.00 -0.69  -0.75
## Na     -0.69 -0.81  0.74 -0.77 -0.79  0.56 -0.69  1.00   0.97
## Conduc -0.76 -0.80  0.76 -0.76 -0.83  0.51 -0.75  0.97   1.00

Listo, con esto hemos indicado que queremos 2 decimales.

Aquí toca buscar las correlaciones más fuertes. Aún así, estos coeficientes (al ser de una muestra para hacer inferencias) no nos dicen si difieren significativamente de 0. (Es decir, si hay suficiente evidencia basada en la muestra para concluir que las correlaciones en la población difieren de cero)

Probando la significancia de la correlación

¿Cómo probamos la significancia estadística?

Recordar nuestras hipótesis: (siempre es el mismo juego de palabras)

  • Hipótesis nula/H0: No hay relación.
  • Hipótesis alternativa/H1: Hay relación.

La función cor.test (x, y, alternative = , method = ) nos ayuda a encontrar el nivel de significancia.

  • x = variable continua explicativa
  • y = variable continua respuesta
  • alternative = referencia a las colas: less, greater, two.side
  • method = tipo de correlación: pearson, spearman, kendall

Podemos entonces a proceder a realizar la prueba de significancia en una de las correlaciones que nos haya llamado la atención. (Por lo general son las que tienen coeficientes fuertes, es decir, lejos del cero)

Recordatorio de las colas:

¿Para qué son?

-Simplemente le indicamos a la prueba si queremos saber si la hipótesis de prueba es que la correlación de la población es menor a cero (usamos “less”). Si queremos saber si es mayor a cero usamos: (a que no adivina!) Y si solo queremos saber si es diferente de cero, usamos “two.side”

cor.test(Soils[,6],Soils[,10])
## 
##  Pearson's product-moment correlation
## 
## data:  Soils[, 6] and Soils[, 10]
## t = 9.3221, df = 46, p-value = 3.614e-12
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.6809493 0.8885997
## sample estimates:
##       cor 
## 0.8086293

¿Qué tenemos?

En la matriz anterior vimos que el pH estaba fuertemente y positivamente relacionado con el Ca por lo que la elegimos para encontrar si la correlación en la población es diferente de cero.

Al hacer al prueba, encontramos que la correlación es significativa, nos da el intervalo de confianza al 95%.

Y de esta manera seguiríamos con las demás correlaciones interesantes que hayamos visto.

Práctica:

Encontrar los niveles de significancia de todas las correlaciones que se encuentran en la matriz que obtuvimos.

….

Claramente sería un proceso sumamente tedioso hacer la función anterior una por una y aquí entra el trabajo comunal genial de R. Una función ya creada (del paquete psych) para obtener todos los p-values con una sola instrucción:

corr.test(Soils[,6:14],use="complete")
## Call:corr.test(x = Soils[, 6:14], use = "complete")
## Correlation matrix 
##           pH     N  Dens     P    Ca    Mg     K    Na Conduc
## pH      1.00  0.64 -0.59  0.59  0.81 -0.40  0.58 -0.69  -0.76
## N       0.64  1.00 -0.86  0.84  0.85 -0.52  0.68 -0.81  -0.80
## Dens   -0.59 -0.86  1.00 -0.79 -0.79  0.49 -0.67  0.74   0.76
## P       0.59  0.84 -0.79  1.00  0.69 -0.49  0.56 -0.77  -0.76
## Ca      0.81  0.85 -0.79  0.69  1.00 -0.43  0.72 -0.79  -0.83
## Mg     -0.40 -0.52  0.49 -0.49 -0.43  1.00 -0.36  0.56   0.51
## K       0.58  0.68 -0.67  0.56  0.72 -0.36  1.00 -0.69  -0.75
## Na     -0.69 -0.81  0.74 -0.77 -0.79  0.56 -0.69  1.00   0.97
## Conduc -0.76 -0.80  0.76 -0.76 -0.83  0.51 -0.75  0.97   1.00
## Sample Size 
## [1] 48
## Probability values (Entries above the diagonal are adjusted for multiple tests.) 
##          pH N Dens P Ca   Mg    K Na Conduc
## pH     0.00 0    0 0  0 0.01 0.00  0      0
## N      0.00 0    0 0  0 0.00 0.00  0      0
## Dens   0.00 0    0 0  0 0.00 0.00  0      0
## P      0.00 0    0 0  0 0.00 0.00  0      0
## Ca     0.00 0    0 0  0 0.01 0.00  0      0
## Mg     0.01 0    0 0  0 0.00 0.01  0      0
## K      0.00 0    0 0  0 0.01 0.00  0      0
## Na     0.00 0    0 0  0 0.00 0.00  0      0
## Conduc 0.00 0    0 0  0 0.00 0.00  0      0
## 
##  To see confidence intervals of the correlations, print with the short=FALSE option

Visualización de datos:

Lo anterior incluía matrices que en ocasiones, pueden hacer que se nos escapen cosas o que pasemos por alto algunas otras. Por lo qué, la visualización de los datos sería mejor.

Scatterplots básicos:

Este es uno de los scatterplots más básicos:

plot(Soils$pH~Soils$Ca)
abline(lm(Soils$pH~Soils$Ca,col="green"))
## Warning: In lm.fit(x, y, offset = offset, singular.ok = singular.ok, ...) :
##  extra argument 'col' will be disregarded

Pero ese representa solo una de las correlaciones que hicimos anteriormente. Podríamos generar un “scatterplot matrix”

pairs(~Soils$pH+Soils$N+Soils$Dens+Soils$P+Soils$Ca+Soils$Mg+Soils$K+Soils$Conduc,cex.labels=1.2)

E inclusive se pueden hacer en 3D. Vamos a elegir tres variables:

plot3d(Soils$pH,Soils$Ca,Soils$Na,main="Gráfico 3D")

Esos son básicos. Existen paquetes, funciones y todo un mundo para generar gráficos. Les recomendamos buscar más en google:

mcor <- cor(Soils[,6:14]) #Hacer selección de columnas

flattenCorrMatrix <- function(cormat, pmat) {
        ut <- upper.tri(cormat)
        data.frame(
                row = rownames(cormat)[row(cormat)[ut]],
                column = rownames(cormat)[col(cormat)[ut]],
                cor  =(cormat)[ut],
                p = pmat[ut]
        )
}

res<-rcorr(as.matrix(mcor))
flattenCorrMatrix(res$r, res$P)
##     row column        cor            p
## 1    pH      N  0.9636558 2.908136e-05
## 2    pH   Dens -0.9591039 4.375513e-05
## 3     N   Dens -0.9966360 7.245975e-09
## 4    pH      P  0.9536132 6.763450e-05
## 5     N      P  0.9926230 1.127174e-07
## 6  Dens      P -0.9896553 3.670119e-07
## 7    pH     Ca  0.9860148 1.050638e-06
## 8     N     Ca  0.9899583 3.308323e-07
## 9  Dens     Ca -0.9876766 6.759023e-07
## 10    P     Ca  0.9752456 7.670472e-06
## 11   pH     Mg -0.9063807 7.535475e-04
## 12    N     Mg -0.9312792 2.617978e-04
## 13 Dens     Mg  0.9191063 4.576975e-04
## 14    P     Mg -0.9279890 3.073469e-04
## 15   Ca     Mg -0.9145676 5.515442e-04
## 16   pH      K  0.9552624 5.968265e-05
## 17    N      K  0.9691990 1.638462e-05
## 18 Dens      K -0.9704107 1.425437e-05
## 19    P      K  0.9513868 7.951878e-05
## 20   Ca      K  0.9774658 5.532619e-06
## 21   Mg      K -0.8929425 1.188790e-03
## 22   pH     Na -0.9797221 3.832844e-06
## 23    N     Na -0.9912987 2.006347e-07
## 24 Dens     Na  0.9833387 1.934020e-06
## 25    P     Na -0.9857819 1.112918e-06
## 26   Ca     Na -0.9920846 1.441625e-07
## 27   Mg     Na  0.9270777 3.208857e-04
## 28    K     Na -0.9779212 5.153385e-06
## 29   pH Conduc -0.9842880 1.576509e-06
## 30    N Conduc -0.9888346 4.790387e-07
## 31 Dens Conduc  0.9827092 2.200715e-06
## 32    P Conduc -0.9815382 2.764885e-06
## 33   Ca Conduc -0.9942982 4.583147e-08
## 34   Mg Conduc  0.9174911 4.896974e-04
## 35    K Conduc -0.9820478 2.508100e-06
## 36   Na Conduc  0.9990872 7.558620e-11
corrplot(mcor, type="upper", order="hclust", tl.col="black", tl.srt=45)

  • Lo anterior es tomado del siguiente enlace

Supuestos:

A través de los scatterplots nos podemos dar cuenta si hemos cometido una violación. Así mismo, existe el “Anscombe´s quartet”, que explica de manera muy intuitiva lo que debemos de evitar en nuestras correlaciones:

Todos los datasets tienen las mismas propiedades estadísticas, pero al graficarlos, vemos cosas distintas.

  • El primer gráfico (dataset I) muestra lo que parece una relación lineal simple, correspondiente a dos variables correlacionadas cumpliendo con la suposición de normalidad.

  • El segundo gráfico (dataset II) no está distribuido normalmente, aunque se observa relación entre los datos, esta no es lineal y el coeficiente de correlación de Pearson no es importante.

  • En la tercera gráfica (dataset III) la distribución es lineal pero hay un outlier que influye como para alterar la línea de correlación de 0.816 hasta una perfecta de 1.

  • la cuarta gráfica (dataset IV) es un ejemplo de muestra en la que un valor atípico es suficiente para producir un coeficiente de correlación alto incluso cuando la relación entre las dos variables no es lineal.

Referencias:

Por si quieren verificar información (Y si encuentran algo erróneo avisarlo por favor)