library(kableExtra)
library(ggplot2)

Coeficiente de Correlación

Vamos ahora a considerar que existe una relación lineal entre dos variables, pero sin asumir que una es funcionalmente dependiente de la otra. En este caso estaremos realizando un análisis de correlación, y usaremos el coeficiente de correlación ( r ), como el estadístico para medir dicha relación.

El coeficiente de correlación simple o de Pearson - entre dos variables, x y y, sin asumir una relación funcional de dependencia entre ellas - se calcula de la siguiente manera:
\[r = \frac{\sum{xy}}{\sqrt{\sum x^2 \sum y^2}}\] Donde: \[\sum{xy} = \sum{X_i Y_i} - \frac{\sum X_i \sum Y_i}{n}\] \[\sum x^2 = \sum X{_i}^2 - \frac{(\sum X_i)^2}{n}\] \[\sum y^2 = \sum Y{_i}^2 - \frac{(\sum Y_i)^2}{n}\]

Aunque el denominador de la primera ecuación ( r ) es siempre positivo, el numerador puede ser positivo, cero o negativo. El valor absoluto del numerador tampoco puede ser mayor que uno. Por lo tanto:
\[-1 \leq r \leq 1\]

Ejemplo 1. Vamos a calcular el coeficiente de correlación de Pearson, para la posible relación lineal entre la longitud del ala (cm) de un especie de ave, y el largo de su cola (cm).

ala <- c(10.4,10.8,11.1,10.2,10.3,10.2,
         10.7,10.5,10.8,11.2,10.6,11.4)
cola <- c(7.4,7.6,7.9,7.2,7.4,7.1,
          7.4,7.2,7.8,7.7,7.8,8.3)
scatter <- data.frame(ala,cola)
ggplot(scatter, aes(ala, cola)) +
  geom_point() +
  stat_ellipse()

Cálculo manual de r

#datos
ala <- c(10.4,10.8,11.1,10.2,10.3,10.2,
         10.7,10.5,10.8,11.2,10.6,11.4)
cola <- c(7.4,7.6,7.9,7.2,7.4,7.1,
          7.4,7.2,7.8,7.7,7.8,8.3)
#cálculos de sumatorias
n <- length(ala)
##X
sumXi <- sum(ala)
sumXi2 <- sum(ala^2)
sumx2 <- sumXi2 - ((sumXi)^2)/n
##Y
sumYi <- sum(cola)
sumYi2 <- sum(cola^2)
sumy2 <- sumYi2 - ((sumYi)^2)/n
#XY
sumXiYi <- sum(ala*cola)
sumxy <- sumXiYi - (sumXi*sumYi)/n
#tablas de resultados parciales
SumX <- data.frame(sumXi,sumYi2,sumx2)
SumY <- data.frame(sumYi,sumYi2,sumy2)
SumXY <- data.frame(sumXiYi,sumxy)
kable(SumX, format = "markdown", col.names = c("Sumatoria  X","Sumatoria X^2", "Suma Cuadrados X"))
Sumatoria X Sumatoria X^2 Suma Cuadrados X
128.2 688.4 1.716667 
kable(SumY, format = "markdown", col.names = c("Sumatoria  Y","Sumatoria Y^2", "Suma Cuadrados Y"))
Sumatoria Y Sumatoria Y^2 Suma Cuadrados Y
90.8 688.4 1.346667 
kable(SumXY, format = "markdown", col.names = c("Sumatoria XY", "Suma Productos XY"))
Sumatoria XY Suma Productos XY
971.37 1.323333 
#Coeficiente de Correlación
coefcorr <- sumxy/sqrt(sumx2*sumy2)
sprintf("r = %.3f", coefcorr)
## [1] "r = 0.870"

También podemos calcular el coeficiente dd determinación, como un indicador de la variabilidad en una variable, tomando en cuenta su correlación con la otra. Matemáticamente es simplemente el cuadrado de r. Por otra parte, el error estándar del coeficiente de correlación se calcula de la siguiente manera: \[s_r = \sqrt \frac{1-r^2}{n-2}\]

#coeficiente de determinación
coefdet <- coefcorr^2
#errorestr
sr <- sqrt((1-coefdet)/(n-2))
coeferr <- data.frame(coefdet,sr)
kable(coeferr, format = "markdown", col.names = c("Coef. Determinación", "Error Estándar r"))
Coef. Determinación Error Estándar r
0.7575171 0.1557186

En el análisis de correlación debemos asumir que ambas variables provienen de una distribución normal, y que los errores conjuntos de desviación de la linealidad, son aleatorios. La normalidad se puede probar mediante una gráfica de q-q y la prueba de Shapito-Wilk.

Prueba de hipótesis para r

El coeficiente de correlación, r, es un estimador del parámetro poblacional \(\rho\). Podemos preguntarnos si ciertamente existe una correlación entre los valores X y los valores Y de la población, y probar \(H_0: \rho=0\).

Una manera directa de hacerlo es utilizando una tabla de valores críticos de r, con n - 2 grados de libertad. Estos valores, en realidad, se basan en el estadístico t, mediante la fórmula: \[r_{\alpha,\nu} = \sqrt \frac{t_{\alpha,\nu}^2}{t_{\alpha,\nu}^2 + \nu}\]

A continuación los cálculos para el ejemplo anterior y la hipótesis: \[H_0: \rho = 0 \quad H_A: \rho \neq 0\]

#valor crítico para r usando una función
#basada en la fórmula r_alfa,nu
critical.r <- function( n, alpha) {
  df <- n - 2
  critical.t <- qt(alpha/2, df, lower.tail = F)
  critical.r <- sqrt( (critical.t^2) / ( (critical.t^2) + df ) )
  return(critical.r)
}
#en nuestro caso n = 12, alfa = 0.05
rcrit <- critical.r(12,0.05)
sprintf("valor crítico de r(0.05,10) = %.3f", rcrit)
## [1] "valor crítico de r(0.05,10) = 0.576"

(Verificar con la tabla)
El valor de \(r\) es 0.870, por lo tanto rechazamos la \(H_0\) y podemos decir que la correlación entre la longitud del ala y la cola es significativa.

Podemos hacer la prueba con el estadístico t, usando la siguiente ecuación para el estadístico de la muestra: \[t = \frac{r}{s_r}\] con n - 2 grados de libertad.

#calculamos el estadístico t
tcalc <- coefcorr/sr
#y el valor crítico de t(alfa,nu)
alfa = 0.05
tcrit <- qt(alfa/2,n-2, lower.tail = F)
tcomp <- data.frame(tcalc,tcrit)
kable(tcomp, format = "markdown", col.names = c("t calculado", "t crítico"))
t calculado t crítico
5.589277 2.228139

Igualmente, el valor calculado es mayor que el crítico, obtenido de la tabla o mediante función en R, por lo tanto podemos rechazar con seguridad la hipótesis nula, y aceptar que hay una correlación entre las variables.

Análisis de Correlación con R

Mediante la función cor.test() de R, podemos realizar el análisis de correlación simple paramétrico (Pearson) y no-paramétricos (Kendall, Spearman). La correlación paramétrica provee la prueba de hipótesis \(H_0: \rho = 0\) mediante el estadístico t, el intervalo de confianza (95 %) para r, y el valor de r. Utilizando cor.test también podemos realizar las pruebas no-paramétricas, o correlaciones de rangos de Kendall y Spearman. La interpretación se hace examinando el valor de p, la probabilidad de cometer error Tipo I.

#prueba normalidad
shapiro.test(ala)
## 
##  Shapiro-Wilk normality test
## 
## data:  ala
## W = 0.94168, p-value = 0.5201
shapiro.test(cola)
## 
##  Shapiro-Wilk normality test
## 
## data:  cola
## W = 0.94449, p-value = 0.5583
#analisis de correlacion paramétricos y no-paramétricos
corranalp <- cor.test(ala,cola, method = "pearson")
corranalp
## 
##  Pearson's product-moment correlation
## 
## data:  ala and cola
## t = 5.5893, df = 10, p-value = 0.0002311
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.5923111 0.9631599
## sample estimates:
##       cor 
## 0.8703546
corranalk <- cor.test(ala,cola, method = "kendall")
corranalk
## 
##  Kendall's rank correlation tau
## 
## data:  ala and cola
## z = 3.1418, p-value = 0.001679
## alternative hypothesis: true tau is not equal to 0
## sample estimates:
##       tau 
## 0.7202074
corranals <- cor.test(ala,cola, method = "spearman")
corranals
## 
##  Spearman's rank correlation rho
## 
## data:  ala and cola
## S = 42.59, p-value = 0.0004467
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.8510852

Referencias

Kabacoff, R. 2015. R in Action. Data analysis and graphics with R, Second Ed. ed. Manning, Shelter Island.

Rosner, B. 2016. Fundamentals of biostatistics, 8th edition. ed. Cengage Learning, Boston, MA.

Zar, J.H. 2014. Biostatistical analysis. Pearson India, Noida.