Cálculo de la covarianza en R

En R, utilizamos la función cov para calcular la covarianza entre dos vectores o columnas de un data.frame.

library(tidyr)


iris %>% 
  summarise(
    cov_petal_length__petal_width = cov(petal_length, petal_width),
    cov_petal_length__sepal_width = cov(petal_length, sepal_width)
  )

## # A tibble: 1 x 2
##   cov_petal_length__petal_width cov_petal_length__sepal_width
##                           <dbl>                         <dbl>
## 1                          1.30                        -0.330

Propiedades de la covarianza

Dada una muestra \((x_1, y_1), (x_2, y_2), ..., (x_n, y_n)\), se cumplen las siguientes propiedades relacionadas con la covarianza \(S_{XY}\):

• Si transformamos linealmente las variables originales \(X'=a+bX\), \(Y'=c+dY\), la covarianza \(S_{X'Y'}\) entre \(X'\), \(Y'\) es la covarianza original multiplicada por \(bd\). Las constantes que se suman no alteran el resultado. En fórmula, tenemos

\[ S_{X'Y'} = bdS_{XY} \]

iris %>% 
  summarise(
    cov_1 = cov(petal_length, petal_width),
    cov_2 = cov(5 + 2 * petal_length, 3 - 4 * petal_width),
    cov_3 = 2 * (-4) * cov(petal_length, petal_width)
  )

## # A tibble: 1 x 3
##   cov_1 cov_2 cov_3
##   <dbl> <dbl> <dbl>
## 1  1.30 -10.4 -10.4

• La covarianza de una variable consigo misma es la varianza de la variable

\[ S_{XX}=S^2_X \]

iris %>% 
  summarise(
    cov = cov(petal_length, petal_length),
    var = var(petal_length)
  )

## # A tibble: 1 x 2
##     cov   var
##   <dbl> <dbl>
## 1  3.12  3.12

• La covarianza entre \(X\) e \(Y\) es igual a la covarianza entre \(Y\) y \(X\)

\[ S_{XY} = S_{YX} \]

iris %>% 
  summarise(
    cov_1 = cov(petal_length, petal_width),
    cov_2 = cov(petal_width, petal_length),
  )

## # A tibble: 1 x 2
##   cov_1 cov_2
##   <dbl> <dbl>
## 1  1.30  1.30

• La covarianza puede calcularse también de la siguiente manera

\[ S_{XY} = \frac{\sum_{i=1}^{n}x_iy_i}{n-1} - \frac{n}{n-1}\bar{X}\bar{Y} \]

cov_alt <- function(.x, .y) {
  n <- length(.x)
  (sum(.x * .y) / (n - 1)) - (n / (n - 1)) * mean(.x) * mean(.y)
}

iris %>% 
  summarise(
    cov_1 = cov(petal_length, petal_width),
    cov_2 = cov_alt(petal_length, petal_width),
  )

## # A tibble: 1 x 2
##   cov_1 cov_2
##   <dbl> <dbl>
## 1  1.30  1.30

• Dada la variable \(Z = X + Y\), podemos calcular la varianza de \(Z\) en función de las varianzas de \(X\), \(Y\) y de su covarianza

\[ S^2_{Z}=S^2_X + 2S_{XY} + S^2_Y \]

iris %>% 
  summarise(
    var_1 = var(petal_length + petal_width),
    var_2 = var(petal_length) + 2*cov(petal_length, petal_width) + var(petal_width),
  )

## # A tibble: 1 x 2
##   var_1 var_2
##   <dbl> <dbl>
## 1  6.29  6.29

• Dada la variable \(Z = X - Y\), podemos calcular la varianza de \(Z\) en función de las varianzas de \(X\), \(Y\) y de su covarianza

\[ S^2_{Z}=S^2_X - 2S_{XY} + S^2_Y \]

• La covarianza de una variable \(X\) con una constante \(c\) es 0

iris %>% 
  mutate(
    constant = 10
  ) %>% 
  summarise(
    cov = cov(petal_length, constant),
  )

## # A tibble: 1 x 1
##     cov
##   <dbl>
## 1     0

Cómo se ve la cóvarianza

Si dibujamos el mismo gráfico agregando líneas para \(\bar{X}\) y \(\bar{Y}\), tendríamos algo como esto:

Para \(X\) comparamos con la línea vertical:
- \((x_i - \bar{X})\) es positivo a la derecha
- \((x_i - \bar{X})\) es negativo a la izquierda

Para \(Y\) comparamos con la línea horizontal:
- \((y_i - \bar{Y})\) es positivo arriba
- \((y_i - \bar{Y})\) es negativo abajo

La mayoría de los puntos están en los cuadrantes “++” y “–”, y en estos cuadrantes \((x_i - \bar{X})(y_i - \bar{Y})\) es positivo; por eso la covarianza es positiva.

Veamos otro caso:

Ahora la mayoría de los puntos están en los cuadrantes “+-” y “-+”, y en esos cuadrantes \((x_i - \bar{X})(y_i - \bar{Y})\) es negativo, así que tendremos una covarianza negativa.

Paradoja de Simpson

Es posible que este último sea un ejemplo de la paradoja de Simpson, en la que la relación entre variables cambia drásticamente cuando se observa dentro de distintos grupos en lugar de observarla en la muestra completa.

La covarianza entre petal_length y sepal_width es -0.3296564 cuando la calculamos usando todos los datos. Veamos qué sucede cuando agrupamos por species.

iris %>% 
  group_by(species) %>% 
  summarise(
    cov_petal_length__sepal_width = cov(petal_length, sepal_width)
  )

## # A tibble: 3 x 2
##   species    cov_petal_length__sepal_width
##   <fct>                              <dbl>
## 1 setosa                            0.0117
## 2 versicolor                        0.0827
## 3 virginica                         0.0714

Y ahora todas las covarianzas de petal_length con sepal_width dentro de cada grupo de species son positivas. Por esta razón, es importante visualizar los datos antes de interpretar un estadístico. Otros ejemplos de por qué es sugerible visualizar los datos lo puedes encontrar en el famoso cuarteto de Ansconbe.