set.seed(12345)
x <- rnorm(n = 1e6, mean = 0, sd = 0.5)
y <- rnorm(n = 1e6, mean = 0, sd = 0.5)
cov(x, y)[1] -0.0004273584Correlaciones
Cuadernos prácticos de Aprendizaje Supervisado I del Grado en Ciencia de Datos Aplicada (curso 2023-2024)
1 Intro a la modelización lineal
El objetivo de la asignatura es ser capaces de construir modelos predictivos para predecir una variable \(Y\) numérica y continua (conocida como variable objetivo o variable dependiente), haciendo uso de un conjunto de predictoras \(X_1,\ldots,X_p\) (conocidas como covariables o variables independientes) tal que
\[Y = f\left( X_1,\ldots,X_p \right) + \varepsilon\]
Concretamente, dado que la asignatura gira en torno a modelos lineales, esa función \(f()\) que relacione las predictoras con la variable objetivo tendrá que ser una función lineal. En matemáticas decimos que una función \(f(x)\) es lineal cuando se cumple:
Propiedad aditiva: \(f(x + y) = f(x) + f(y)\)
Propiedad homogénea: \(f(k*x) = k*f(x)\) (donde \(k\) es una constante en \(\mathbb{R}\)).
Ambas se pueden resumir en \(f(a*x + b*y) = a*f(x) + b*f(y)\)
Ejemplos lineales: \(y = 2*x_1 - 3\) o \(y = 4 - \frac{x_1}{2} + 3*x_2\)
Ejemplos no lineales: \(y = 2*\frac{1}{x_1}\) o \(y = 4 - x_{1}^{2} - x_2\) o \(y = ln(x_1) + cos(x_2)\)
Así, nuestra regresión lineal será:
\[Y = \beta_0 + \beta_1 X_1 + \beta_p X_p + \varepsilon, \quad E[\varepsilon] = 0\] lo que traduce en
\[E[Y|X = x] = \beta_0 + \beta_1 x_1 + \beta_p x_p\]
La idea es encontrar el «mejor» modelo posible atendiendo a diferentes aspectos
¿Cuánto se equivoca? (conocemos el error ya que es aprendizaje supervisado: de mi conjunto con el que entreno el modelo, conozco los valores reales de \(Y\). 📚 Ver «The elements of Statistical Learning» Hastie et al., 2008)
¿Cómo de complicado es el modelo? ¿Existe alguno más sencillo (equivocándose similar)?
Amén de lo que se equivoca, ¿cómo varían dichas equivocaciones? ¿Qué varianza tienen mis estimaciones?
¿Qué % de información estoy capturando?
Ese modelo lo podemos reducir a encontrar los «mejores» estimadores \(\left( \hat{\beta}_0, \ldots, \hat{\beta}_p \right)\)
\[\widehat{E[Y|X = x]} = \hat{\beta}_0 + \hat{\beta}_1 x_1 + \ldots + \hat{\beta}_p x_p\]
Así, nuestro error estimado se podrá calcular como
\[\hat{\varepsilon}= Y - \hat{Y}, \quad E[\hat{\varepsilon}] = 0\]
1.1 Caso particular: regresión univariante
El caso univariante es el más sencillo y es aquel en el que solo tenemos una variable predictora \(X_1 := X\), tal que
\[Y = \beta_0 + \beta_1 X_1 + \varepsilon := \beta_0 + \beta_1 X + \varepsilon\]
El objetivo en este primera aproximación es, dado un dataset con una variable objetivo y múltiples predictoras, ser capaces de decidir la predictora cuyo efecto lineal en \(Y\) sea más relevante. Y para ello debemos saber un concepto fundamental estadística: la covarianza y la correlación.
2 Covarianza
2.1 Repaso teórico
Es habitual en estadística hablar de medidas de centralización (media, moda, mediana) y de medidas de dispersión, que cuantifican cómo varían los datos y en torno a qué valor lo hacen.
La medida de dispersión más famosa es la quizás la varianza. ¿Qué es en realidad la varianza?
\[s_{x}^{2} = \frac{1}{n} \sum_{i=1}^{n} \left(x_i - \overline{x} \right)^2 = \overline{x^2} - \overline{x}^2 \geq 0\]
Si atendemos simplemente a la fórmula y su definición, la varianza es el promedio de las desviaciones al cuadrado (respecto a la media), lo que nos permite cuantificar de alguna manera la relación de una variable CONSIGO MISMA. Si te fijas en la fórmula, podríamos separar los cuadrados de la siguiente manera:
\[s_{x}^{2} = \frac{1}{n} \sum_{i=1}^{n} \left(x_i - \overline{x} \right)\left(x_i - \overline{x} \right) = \overline{x*x} - \overline{x}*\overline{x} \geq 0\]
¿Y si quiésemos medir la relación de una variable X respecto a otra variable (en lugar de consigo misma)?
Esta es la base del concepto de la covarianza: sustituir en la fórmula anterior una de las desviaciones de x por las desviaciones respecto a OTRA variable
\[s_{xy} = \frac{1}{n} \sum_{i=1}^{n} \sum_{j=1}^{n} \left(x_i - \overline{x} \right)\left(y_j - \overline{y} \right) = \overline{x*y} - \overline{x}*\overline{y}\]
2.1.1 Propiedades de la covarianza
Es importante entender algunas propiedades de la covarianza
Signo: la covarianza puede ser tanto positiva como negativa como 0: al eliminar el cuadrado de la varianza, ya no es necesario que sea positiva
¿Qué cuantifica? La covarianza SOLO mide la relación LINEAL (en torno a una recta) entre dos variables, pero no captura relaciones no lineales
¿Qué dice su signo? El signo de la covarianza nos indicará la dirección de la dependencia lineal: si es positiva, la relación será creciente (cuando X crece, Y crece); si es negativa, la relación será decreciente (cuando X crece, Y decrece)
¿Qué dice su magnitud? Al igual que pasa con la varianza, la covarianza depende de las unidades y magnitudes de los datos, así que no podemos comparar covarianzas de variables y datasets distintos: si me dicen que la covarianza es de, por ejemplo, 5.5, no sabemos si es mucho o poco.
2.2 Herramientas en R
Para calcular la covarianza entre dos variables en R la forma más sencilla es usar la función cov().
Importante
Recuerda que los softwares estadísticos nos devuelven siempre la cuasi covarianza, dividido entre \(n-1\) y no entre \(n\). La cuasivarianza y la cuasicovarianza son los mejores estimadores muestrales (insesgados) de los respectivos parámetros poblaciones
2.2.1 Simulación: variables sin relación
Por ejemplo, vamos a simular dos variables \(X \sim N(0, 0.5)\) e \(Y \sim N(0, 0.5)\) que no tienen ningún tipo de relación (ni lineal ni no lineal entre ellas)
Vemos como efectivamente la covarianza es prácticamente cero, ¿pero cómo cercano a 0 sería?
2.2.2 Simulación: variables sin relación lineal
¿Y si simulamos una variable \(X\) de manera aleatoria y otra variable \(Y\) tal que \(Y = X^2\)?
set.seed(12345)
x <- rnorm(n = 1e6, mean = 0, sd = 0.5)
y <- x^2
cov(x, y)[1] 0.0002907904Aunque haya una relación de dependencia perfecta entre X e Y (una función que las relaciona), ¡la covarianza puede que sea incluso más pequeña que antes (en valor absoluto)!.
Esto se debe a que la covarianza/correlación no captura relaciones no lineales, así que para saber si son independientes realmente deberíamos de pasar un contraste de independencia. Para ello podemos hacer uso de chisq.test() (contraste no parámetrico).
Esto lo podemos ilustrar con ejemplo más sencillo:
- \(X = -1, 0, 1\) cuya media es \(\overline{x} = 0\)
- \(Y = X^2 = 1, 0, 1\) cuya media es \(\overline{y} = \frac{2}{3}\)
- \(XY = X*X^2 = X^3 = -1, 0, 1\) cuya media es \(\overline{xy} = 0\)
Si hacemos la covarianza
\[s_{xy} = \overline{xy} - \overline{x}{y} = 0 - 0 \frac{2}{3} = 0\]
Moraleja: covarianza nula NO IMPLICA independencia (independencia si implica, en particular, ausencia de relación lineal)
2.2.3 Simulación: variables con relación lineal
¿Y si simulamos dos variables con relación lineal?
set.seed(12345)
x <- rnorm(n = 1e6, mean = 0, sd = 0.5)
y <- 2*x
cov(x, y)[1] 0.5013159La covarianza en este caos es 0.50132 y la pregunta es…¿es mucho o poco? Eso no podemos resolverlo hasta hablar de correlación: de momento solo podemos decir que, en caso de existir, la relación lineal será positiva. Esto cambiaría si la hubiéramos definido \(Y\) de otra manera, por ejemplo y <- -2*x
y <- -2*x
cov(x, y)[1] -0.50131592.3 Matriz de covarianzas
Antes de hablar de correlación, vamos a introducir una generalización de una simple covarianza, y es la conocida como matriz de covarianzas, y que jugará un papel fundamental en estadística ya que resume un conjunto de datos multivariantes.
La idea de dicha matriz que llamaremos \(\Sigma\) es, que si en lugar de tener una sola predictora tenemos \(\left(X_1, \ldots, X_p \right)\) covariables, cada elemento \(\Sigma_{i,j}\) se defina como
\[\Sigma = \left(\Sigma_{i,j} \right)_{i=1,\ldots,p}^{j=1,\ldots,p} = \begin{pmatrix} s_{x_1 x_1} = s_{x_1}^2 & s_{x_1 x_2} & \ldots & s_{x_1 x_p} \\ s_{x_2 x_1} & s_{x_2 x_2} = s_{x_2}^2 & \ldots & s_{x_2 x_p} \\ \vdots & \vdots & \ddots & \vdots \\ s_{x_p x_1} & s_{x_p x_2} & \ldots & s_{x_p x_p} = s_{x_p}^2 \end{pmatrix}\]
En el futuro veremos como los autovalores y autovectores de esta matriz serán fundamentales. Para calcularla en R basta con pasarle cov() a un dataset (¡importante, solo numéricas!)
library(tidyverse)
starwars |>
select(where(is.numeric)) |> # importante: solo numéricas
drop_na() |> # importante: covarianza de ausentes --> ausente
cov() height mass birth_year
height 957.3802 676.8867 -2164.10
mass 676.8867 46313.2034 17402.55
birth_year -2164.1002 17402.5466 28603.042.4 Datos agrupados
Un pequeño matiz en caso de que tuviésemos que hacer el cálculo a mano: si tenemos una base de datos agrupados (no tenemos el dato en bruto sino los pares de datos únicos \((x_i, y_j)\) y su frecuencia \(n_{ij}\)), la matriz anterior se calcula como
\[\Sigma = \left(\Sigma_{i,j} \right)_{i=1,\ldots,p}^{j=1,\ldots,p}, \quad \Sigma_{i,j} = s_{x_i x_j} = \frac{1}{n} \sum_{i=1}^{k_i} \sum_{j=1}^{k_j} n_{ij} \left(x_i - \overline{x} \right)\left(y_j - \overline{y} \right)\]
3 Correlación
3.1 Repaso teórico
Como decíamos, al igual que pasaba con la varianza, la covarianza depende de las unidades y magnitudes de los datos, así que lo que haremos será estandarizar la covarianza.
Definiremos el coeficiente correlación lineal (de Pearson) como la covarianza dividida entre el producto de las desviaciones típicas (algo adimensional)
\[r_{xy} = \rho_{xy} = \frac{s_{xy}}{s_x s_y} = \frac{\overline{x*y} - \overline{x} * \overline{y}}{\sqrt{\overline{x^2} - \overline{x}^2}\sqrt{\overline{y^2} - \overline{y}^2}}\]
3.1.1 Propiedades de la correlación
Es importante entender algunas propiedades de la correlación
Signo: mismo signo que la covarianza ya que el denominador es siempre positivo.
¿Qué cuantifica?: sus valores siempre están entre -1 y 1 por lo que nos sirve como una escala de la fortaleza de una relación lineal. Al igual qe con la covarianza no captura relaciones no lineales
¿Qué dice su signo? El signo de la correlación nos indicará la dirección de la dependencia lineal: si es positiva, la relación será creciente (cuando X crece, Y crece); si es negativa, la relación será decreciente (cuando X crece, Y decrece)
¿Qué dice su magnitud? Más cerca de -1 o 1 implica (a priori…) una relación más fuerte, más cerca de 0 ausencia (a priori…) de relación lineal
3.2 Herramientas en R
Para calcular la correlación entre dos variables en R la forma más sencilla es usar la función cor().
Importante
Aunque los softwares estadísticos nos devuelven siempre la cuasi covarianza, dividido entre \(n-1\) y no entre \(n\), esto no afecta la correlación (ya que los denominadores se cancelan)
Como hemos hecho antes con la covarianza, podemos tener una generalización de una simple correlación, y es la conocida como matriz de correlaciones, y que jugará también un papel fundamental en estadística ya que resume un conjunto de datos multivariantes (estandrizados).
La idea de dicha matriz que llamaremos \(R\) es, que si en lugar de tener una sola predictora tenemos \(\left(X_1, \ldots, X_p \right)\) covariables, cada elemento \(r_{i,j}\) se defina como
\[R = \left(r_{i,j} \right)_{i=1,\ldots,p}^{j=1,\ldots,p} = \begin{pmatrix} 1 & r_{x_1 x_2} & \ldots & r_{x_1 x_p} \\ r_{x_2 x_1} & 1 & \ldots & r_{x_2 x_p} \\ \vdots & \vdots & \ddots & \vdots \\ r_{x_p x_1} & r_{x_p x_2} & \ldots & 1 \end{pmatrix}\]
Fíjate que ahora, además de ser simétrica como la matriz de covarianzas, la diagonal es siempre 1 (la correlación entre una variable consigo misma es la máxima posible). Para calcularla en R basta con pasarle cor() a un dataset (¡importante, solo numéricas!)
library(tidyverse)
starwars |>
select(where(is.numeric)) |> # importante: solo numéricas
drop_na() |> # importante: correlación de ausentes --> ausente
cor() height mass birth_year
height 1.0000000 0.1016533 -0.4135510
mass 0.1016533 1.0000000 0.4781391
birth_year -0.4135510 0.4781391 1.0000000
Una propiedad importante que relaciona ambas matrices. Si tenemos los datos estandarizados por media y varianza, la matriz de covarianzas coincide con la de correlaciones
\[R\left(X_1, \ldots, X_p \right) = \Sigma \left(\frac{X_1 - \overline{X_1}}{s_{X_1}}, \ldots, \frac{X_p - \overline{X_p}}{s_{X_p}} \right)\]
La matriz de correlaciones no es más que la matriz de covarianzas cuando los datos están estandarizados
library(tidyverse)
# Correlación
starwars |>
select(where(is.numeric)) |>
drop_na() |>
cor() height mass birth_year
height 1.0000000 0.1016533 -0.4135510
mass 0.1016533 1.0000000 0.4781391
birth_year -0.4135510 0.4781391 1.0000000# Covarianza con datos estandarizados
starwars |>
select(where(is.numeric)) |>
drop_na() |>
# Estandarizamos
mutate(across(everything(), function(x) { (x - mean(x))/sd(x) })) |>
cov() height mass birth_year
height 1.0000000 0.1016533 -0.4135510
mass 0.1016533 1.0000000 0.4781391
birth_year -0.4135510 0.4781391 1.0000000Efectivamente es la misma matriz.
3.2.1 Herramientas extras: paquete corrr y paquete corrplot
Para poder hacer un mejor uso de las correlaciones contamos con el paquete {corrr} que nos proporciona la matriz de correlaciones en un tibble propio, pudiendo realizar una profundo análisis exploratorio de la misma. La función más importante es correlate() (el equivalente a cor()), que vamos a aplicar al conjunto mtcars del paquete {datasets}
library(corrr)
mtcars |>
correlate()Correlation computed with
• Method: 'pearson'
• Missing treated using: 'pairwise.complete.obs'# A tibble: 11 × 12
term mpg cyl disp hp drat wt qsec vs am
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 mpg NA -0.852 -0.848 -0.776 0.681 -0.868 0.419 0.664 0.600
2 cyl -0.852 NA 0.902 0.832 -0.700 0.782 -0.591 -0.811 -0.523
3 disp -0.848 0.902 NA 0.791 -0.710 0.888 -0.434 -0.710 -0.591
4 hp -0.776 0.832 0.791 NA -0.449 0.659 -0.708 -0.723 -0.243
5 drat 0.681 -0.700 -0.710 -0.449 NA -0.712 0.0912 0.440 0.713
6 wt -0.868 0.782 0.888 0.659 -0.712 NA -0.175 -0.555 -0.692
7 qsec 0.419 -0.591 -0.434 -0.708 0.0912 -0.175 NA 0.745 -0.230
8 vs 0.664 -0.811 -0.710 -0.723 0.440 -0.555 0.745 NA 0.168
9 am 0.600 -0.523 -0.591 -0.243 0.713 -0.692 -0.230 0.168 NA
10 gear 0.480 -0.493 -0.556 -0.126 0.700 -0.583 -0.213 0.206 0.794
11 carb -0.551 0.527 0.395 0.750 -0.0908 0.428 -0.656 -0.570 0.0575
# ℹ 2 more variables: gear <dbl>, carb <dbl>¿Qué diferencias hay?
- La salida es a su vez un tibble por lo que podemos filtrar de manera muy sencilla quedándonos solo con la información que nos interesa
# Filtramos solo las correlaciones respecto a 3 variables
mtcars |>
correlate() |>
filter(term %in% c("mpg", "disp", "gear"))Correlation computed with
• Method: 'pearson'
• Missing treated using: 'pairwise.complete.obs'# A tibble: 3 × 12
term mpg cyl disp hp drat wt qsec vs am gear
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 mpg NA -0.852 -0.848 -0.776 0.681 -0.868 0.419 0.664 0.600 0.480
2 disp -0.848 0.902 NA 0.791 -0.710 0.888 -0.434 -0.710 -0.591 -0.556
3 gear 0.480 -0.493 -0.556 -0.126 0.700 -0.583 -0.213 0.206 0.794 NA
# ℹ 1 more variable: carb <dbl># Si queremos predecir mpg, filtramos aquellas variables que tienen una correlación (en valor absoluto) superior a 0.7
mtcars |>
correlate() |>
filter(abs(mpg) > 0.7)Correlation computed with
• Method: 'pearson'
• Missing treated using: 'pairwise.complete.obs'# A tibble: 4 × 12
term mpg cyl disp hp drat wt qsec vs am gear
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 cyl -0.852 NA 0.902 0.832 -0.700 0.782 -0.591 -0.811 -0.523 -0.493
2 disp -0.848 0.902 NA 0.791 -0.710 0.888 -0.434 -0.710 -0.591 -0.556
3 hp -0.776 0.832 0.791 NA -0.449 0.659 -0.708 -0.723 -0.243 -0.126
4 wt -0.868 0.782 0.888 0.659 -0.712 NA -0.175 -0.555 -0.692 -0.583
# ℹ 1 more variable: carb <dbl>Los valores negativos vienen marcados en rojo
Podemos especificar que valor queremos en la diagonal
mtcars |>
correlate(diag = 1)Correlation computed with
• Method: 'pearson'
• Missing treated using: 'pairwise.complete.obs'# A tibble: 11 × 12
term mpg cyl disp hp drat wt qsec vs am
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 mpg 1 -0.852 -0.848 -0.776 0.681 -0.868 0.419 0.664 0.600
2 cyl -0.852 1 0.902 0.832 -0.700 0.782 -0.591 -0.811 -0.523
3 disp -0.848 0.902 1 0.791 -0.710 0.888 -0.434 -0.710 -0.591
4 hp -0.776 0.832 0.791 1 -0.449 0.659 -0.708 -0.723 -0.243
5 drat 0.681 -0.700 -0.710 -0.449 1 -0.712 0.0912 0.440 0.713
6 wt -0.868 0.782 0.888 0.659 -0.712 1 -0.175 -0.555 -0.692
7 qsec 0.419 -0.591 -0.434 -0.708 0.0912 -0.175 1 0.745 -0.230
8 vs 0.664 -0.811 -0.710 -0.723 0.440 -0.555 0.745 1 0.168
9 am 0.600 -0.523 -0.591 -0.243 0.713 -0.692 -0.230 0.168 1
10 gear 0.480 -0.493 -0.556 -0.126 0.700 -0.583 -0.213 0.206 0.794
11 carb -0.551 0.527 0.395 0.750 -0.0908 0.428 -0.656 -0.570 0.0575
# ℹ 2 more variables: gear <dbl>, carb <dbl>Además contamos con algunas funciones que nos pueden ayudar:
focus(): realiza un filtro específico de las variables de interés
mtcars |>
correlate() |>
focus(mpg, cyl)Correlation computed with
• Method: 'pearson'
• Missing treated using: 'pairwise.complete.obs'# A tibble: 9 × 3
term mpg cyl
<chr> <dbl> <dbl>
1 disp -0.848 0.902
2 hp -0.776 0.832
3 drat 0.681 -0.700
4 wt -0.868 0.782
5 qsec 0.419 -0.591
6 vs 0.664 -0.811
7 am 0.600 -0.523
8 gear 0.480 -0.493
9 carb -0.551 0.527Podemos indicar que variables queremos quitar con un signo negativo sobre ellas y mirror = TRUE
mtcars |>
correlate() |>
focus(-(mpg:drat), mirror = TRUE)Correlation computed with
• Method: 'pearson'
• Missing treated using: 'pairwise.complete.obs'# A tibble: 6 × 7
term wt qsec vs am gear carb
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 wt NA -0.175 -0.555 -0.692 -0.583 0.428
2 qsec -0.175 NA 0.745 -0.230 -0.213 -0.656
3 vs -0.555 0.745 NA 0.168 0.206 -0.570
4 am -0.692 -0.230 0.168 NA 0.794 0.0575
5 gear -0.583 -0.213 0.206 0.794 NA 0.274
6 carb 0.428 -0.656 -0.570 0.0575 0.274 NA rearrange(): nos permite ordenar por correlación, de manera que nos agrupa las variables más altamente correladas juntas
mtcars |>
correlate() |>
rearrange(absolute = TRUE)Correlation computed with
• Method: 'pearson'
• Missing treated using: 'pairwise.complete.obs'# A tibble: 11 × 12
term mpg vs drat am gear qsec carb hp wt
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 mpg NA 0.664 0.681 0.600 0.480 0.419 -0.551 -0.776 -0.868
2 vs 0.664 NA 0.440 0.168 0.206 0.745 -0.570 -0.723 -0.555
3 drat 0.681 0.440 NA 0.713 0.700 0.0912 -0.0908 -0.449 -0.712
4 am 0.600 0.168 0.713 NA 0.794 -0.230 0.0575 -0.243 -0.692
5 gear 0.480 0.206 0.700 0.794 NA -0.213 0.274 -0.126 -0.583
6 qsec 0.419 0.745 0.0912 -0.230 -0.213 NA -0.656 -0.708 -0.175
7 carb -0.551 -0.570 -0.0908 0.0575 0.274 -0.656 NA 0.750 0.428
8 hp -0.776 -0.723 -0.449 -0.243 -0.126 -0.708 0.750 NA 0.659
9 wt -0.868 -0.555 -0.712 -0.692 -0.583 -0.175 0.428 0.659 NA
10 disp -0.848 -0.710 -0.710 -0.591 -0.556 -0.434 0.395 0.791 0.888
11 cyl -0.852 -0.811 -0.700 -0.523 -0.493 -0.591 0.527 0.832 0.782
# ℹ 2 more variables: disp <dbl>, cyl <dbl>shave(): nos muestra solo la diagonal inferior
mtcars |>
correlate() |>
rearrange(absolute = TRUE) |>
shave()Correlation computed with
• Method: 'pearson'
• Missing treated using: 'pairwise.complete.obs'# A tibble: 11 × 12
term mpg vs drat am gear qsec carb hp wt disp
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 mpg NA NA NA NA NA NA NA NA NA NA
2 vs 0.664 NA NA NA NA NA NA NA NA NA
3 drat 0.681 0.440 NA NA NA NA NA NA NA NA
4 am 0.600 0.168 0.713 NA NA NA NA NA NA NA
5 gear 0.480 0.206 0.700 0.794 NA NA NA NA NA NA
6 qsec 0.419 0.745 0.0912 -0.230 -0.213 NA NA NA NA NA
7 carb -0.551 -0.570 -0.0908 0.0575 0.274 -0.656 NA NA NA NA
8 hp -0.776 -0.723 -0.449 -0.243 -0.126 -0.708 0.750 NA NA NA
9 wt -0.868 -0.555 -0.712 -0.692 -0.583 -0.175 0.428 0.659 NA NA
10 disp -0.848 -0.710 -0.710 -0.591 -0.556 -0.434 0.395 0.791 0.888 NA
11 cyl -0.852 -0.811 -0.700 -0.523 -0.493 -0.591 0.527 0.832 0.782 0.902
# ℹ 1 more variable: cyl <dbl>fashion(): nos ofrece un resumen más limpio y resumido
mtcars |>
correlate() |>
shave() |>
fashion()Correlation computed with
• Method: 'pearson'
• Missing treated using: 'pairwise.complete.obs' term mpg cyl disp hp drat wt qsec vs am gear carb
1 mpg
2 cyl -.85
3 disp -.85 .90
4 hp -.78 .83 .79
5 drat .68 -.70 -.71 -.45
6 wt -.87 .78 .89 .66 -.71
7 qsec .42 -.59 -.43 -.71 .09 -.17
8 vs .66 -.81 -.71 -.72 .44 -.55 .74
9 am .60 -.52 -.59 -.24 .71 -.69 -.23 .17
10 gear .48 -.49 -.56 -.13 .70 -.58 -.21 .21 .79
11 carb -.55 .53 .39 .75 -.09 .43 -.66 -.57 .06 .27 network_plot(): nos permite visualizar las relaciones lineales mediante un grafo (podemos decidir el umbral de correlación conmin_cor = ...)
mtcars |>
correlate() |>
network_plot(min_cor = 0.7)Correlation computed with
• Method: 'pearson'
• Missing treated using: 'pairwise.complete.obs'
Además del paquete {corrr} podemos hacer uso de {corrrplot} para una mejor visualización (importante: necesita la matriz de correlaciones sin formato)
library(corrplot)corrplot 0.92 loadedmtcars |>
cor() |>
corrplot()
Dando valores distintos al argumento method = ... podemos obtener visualizaciones diferentes.
mtcars |>
cor() |>
corrplot(method = "color")
mtcars |>
cor() |>
corrplot(method = "number")
mtcars |>
cor() |>
corrplot(method = "ellipse")
Podemos hacer uso también del argumento type = ... para visualizar solo una mitad de la matriz (simétrica)
mtcars |>
cor() |>
corrplot(method = "ellipse")
Incluso podemos elegir la paleta de colores con col = ... y las opciones del paquete {RColorBrewer} para construir paletas
library(RColorBrewer)
mtcars |>
cor() |>
corrplot(type="upper", col = brewer.pal(n = 8, name = "PuOr"))
Podemos incluso realizar un contraste de hipótesis cuya hipótesis nula es que la correlación es nula
cor_test <-
mtcars |>
cor.mtest()
cor_test$p # matriz de p-valores mpg cyl disp hp drat
mpg 0.000000e+00 6.112687e-10 9.380327e-10 1.787835e-07 1.776240e-05
cyl 6.112687e-10 0.000000e+00 1.802838e-12 3.477861e-09 8.244636e-06
disp 9.380327e-10 1.802838e-12 0.000000e+00 7.142679e-08 5.282022e-06
hp 1.787835e-07 3.477861e-09 7.142679e-08 0.000000e+00 9.988772e-03
drat 1.776240e-05 8.244636e-06 5.282022e-06 9.988772e-03 0.000000e+00
wt 1.293959e-10 1.217567e-07 1.222320e-11 4.145827e-05 4.784260e-06
qsec 1.708199e-02 3.660533e-04 1.314404e-02 5.766253e-06 6.195826e-01
vs 3.415937e-05 1.843018e-08 5.235012e-06 2.940896e-06 1.167553e-02
am 2.850207e-04 2.151207e-03 3.662114e-04 1.798309e-01 4.726790e-06
gear 5.400948e-03 4.173297e-03 9.635921e-04 4.930119e-01 8.360110e-06
carb 1.084446e-03 1.942340e-03 2.526789e-02 7.827810e-07 6.211834e-01
wt qsec vs am gear
mpg 1.293959e-10 1.708199e-02 3.415937e-05 2.850207e-04 5.400948e-03
cyl 1.217567e-07 3.660533e-04 1.843018e-08 2.151207e-03 4.173297e-03
disp 1.222320e-11 1.314404e-02 5.235012e-06 3.662114e-04 9.635921e-04
hp 4.145827e-05 5.766253e-06 2.940896e-06 1.798309e-01 4.930119e-01
drat 4.784260e-06 6.195826e-01 1.167553e-02 4.726790e-06 8.360110e-06
wt 0.000000e+00 3.388683e-01 9.798492e-04 1.125440e-05 4.586601e-04
qsec 3.388683e-01 0.000000e+00 1.029669e-06 2.056621e-01 2.425344e-01
vs 9.798492e-04 1.029669e-06 0.000000e+00 3.570439e-01 2.579439e-01
am 1.125440e-05 2.056621e-01 3.570439e-01 0.000000e+00 5.834043e-08
gear 4.586601e-04 2.425344e-01 2.579439e-01 5.834043e-08 0.000000e+00
carb 1.463861e-02 4.536949e-05 6.670496e-04 7.544526e-01 1.290291e-01
carb
mpg 1.084446e-03
cyl 1.942340e-03
disp 2.526789e-02
hp 7.827810e-07
drat 6.211834e-01
wt 1.463861e-02
qsec 4.536949e-05
vs 6.670496e-04
am 7.544526e-01
gear 1.290291e-01
carb 0.000000e+00Esa matriz de p-valores podemos incluirla en corrplot() para que directamente nos marque las variables cuya hipótesis nula de incorrelación no puede rechazarse
mtcars |>
cor() |>
corrplot(method = "color", p.mat = cor_test$p, sig.level = 0.05) 
3.2.2 Simulación: variables sin relación
Vamos a volver a simular las dos variables \(X \sim N(0, 0.5)\) e \(Y \sim N(0, 0.5)\) que no tienen ningún tipo de relación (ni lineal ni no lineal entre ellas)
set.seed(12345)
x <- rnorm(n = 1e6, mean = 0, sd = 0.5)
y <- rnorm(n = 1e6, mean = 0, sd = 0.5)
cor(x, y)[1] -0.001707521Ahora sí que pdoemos decir que, efectivamente, la relación lineal es inexistente ya que la correlación es prácticamente nula
3.2.3 Simulación: variables sin relación lineal
Simulamos de nuevo una variable \(X\) de manera aleatoria y otra variable \(Y\) tal que \(Y = X^2\)
set.seed(12345)
x <- rnorm(n = 1e6, mean = 0, sd = 0.5)
y <- x^2
cor(x, y)[1] 0.001639215Aunque haya una relación de dependencia perfecta entre X e Y (una función que las relaciona), la correlación vuelve a ser prácticamente nula.
3.2.4 Simulación: variables con relación lineal
Y vamos a volver a simular dos variables con relación lineal
set.seed(12345)
x <- rnorm(n = 1e6, mean = 0, sd = 0.5)
y <- 2*x
cor(x, y)[1] 1La covarianza en este caso es 1, ya que tenemos una relación lineal perfecta dado que ambas variables vienen relacionadas por una fórmula
3.2.5 Simulación con ruido
¿Y si a la anterior relación lineal le incorporamos ruido, con una varianza cada vez más grande?
set.seed(12345)
x <- rnorm(n = 1e5, mean = 0, sd = 0.5)
sigmas <- seq(0.0001, 10, l = 100)
correlations <- tibble("sigmas" = sigmas, "cor" = NA)
for (i in 1:length(sigmas)) {
# vamos añadiendo ruido y calculamos la correlación
y <- 2*x + rnorm(n = 1e5, mean = 0, sd = sigmas[i])
correlations$cor[i] <- cor(x, y)
}
ggplot(correlations, aes(x = sigmas, y = cor)) +
geom_point(aes(color = sigmas), size = 3, alpha = 0.7) +
labs(title = "Correlación vs ruido",
x = "sigma", y = "correlación",
color = "sd ruido") +
theme_minimal()
La correlación captura la relación lineal pero si entre ambas variables hay demás algo aleatorio que no podemos capturar, según aumente la varianza del ruido, la correlación será más débil ya que perturbará la dependencia entre ellas.
4 Correlación vs causalidad
Por último, vamos a mencionar un matiz importante: la idea de causalidad
Imagina que \(X\) e \(Y\) fuesen variables dependientes (lineales o no) como en la imagen…¿implicaría que el ascenso/descenso de una PROVOCA el ascenso/descenso de la otra?
En el caso de la imagen el nivel de bronceado (X) y el consumo de helados (Y). parece que son dependientes? Pero…¿una causa la otra?
Diremos que dos variables tienen una relación causal o de causalidad cuando una variable es la CAUSA del comportamiento de la otra, algo que no sabremos solo con estadística (sino con conocimiento experto, en este caso de nutricionistas y médicos). Y es que, mal que nos pese a veces a los matemáticos/estadísticos/científicos de datos, correlación NO IMPLICA causalidad (al menos no tiene por qué)
Otro ejemplo podemos verlo en la siguiente gráfica que, aparentemente, nos dice que comer chocolate nos hará ganar el Nobel. Ojalá :P
Este fenómeno es conocido como correlaciones espúreas, y aparece cuando dos variables presentan una relación matemática pero sin ningún tipo de relación causal o lógica. Puedes ver más en https://www.tylervigen.com/spurious-correlations
Dicho patrón matemático puede deberse al mero azar o la existencia de lo que se conoce como variables de confusión (en el caso del helado, el verano es el que causa ambas).
La rama que se dedica al estudio de la causalidad, combinando las ramas de la estadística, filosofía, sociología y psicología se le conoce como inferencia causal y es fundamental en un análisis más profundo de las relaciones entre las variables (sobre todo en campos como la economía, la psicología o la sociología).
Recuerda
Así que recuerda: un buen estadístico/a o científico/a de datos debe ser humilde en sus atribuciones. Nuestro trabajo llega hasta encontrar patrones. Deja a otros expertos que sean los que, con su conocimiento de la matería, atribuyan causas y consecuencias.
5 Caso práctico I: anscombe
En el paquete {datasets} se encuentra el dataset conocido como cuarteto de Anscombe, un dataset que cuenta en realidad con 4 conjuntos de datos.
library(tidyverse)
anscombe_tb <- as_tibble(datasets::anscombe)
anscombe_tb# A tibble: 11 × 8
x1 x2 x3 x4 y1 y2 y3 y4
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 10 10 10 8 8.04 9.14 7.46 6.58
2 8 8 8 8 6.95 8.14 6.77 5.76
3 13 13 13 8 7.58 8.74 12.7 7.71
4 9 9 9 8 8.81 8.77 7.11 8.84
5 11 11 11 8 8.33 9.26 7.81 8.47
6 14 14 14 8 9.96 8.1 8.84 7.04
7 6 6 6 8 7.24 6.13 6.08 5.25
8 4 4 4 19 4.26 3.1 5.39 12.5
9 12 12 12 8 10.8 9.13 8.15 5.56
10 7 7 7 8 4.82 7.26 6.42 7.91
11 5 5 5 8 5.68 4.74 5.73 6.89Para su correcto manejo vamos a convertirlo a tidydata
anscombe_x <-
anscombe_tb |>
pivot_longer(cols = x1:x4, names_to = "dataset", values_to = "x") |>
select(-contains("y")) |>
mutate(dataset = str_remove_all(dataset, "x"))
anscombe_y <-
anscombe_tb |>
pivot_longer(cols = y1:y4, names_to = "dataset", values_to = "y") |>
select(-contains("x"), -dataset)
anscombe_tidy <-
anscombe_x |> mutate("y" = anscombe_y$y)
anscombe_tidy# A tibble: 44 × 3
dataset x y
<chr> <dbl> <dbl>
1 1 10 8.04
2 2 10 9.14
3 3 10 7.46
4 4 8 6.58
5 1 8 6.95
6 2 8 8.14
7 3 8 6.77
8 4 8 5.76
9 1 13 7.58
10 2 13 8.74
# ℹ 34 more rows¿Qué características muestrales tienen? Para conocerlos un poco mejor vamos a calcular la media, varianza, desv. típica, covarianza y correlación en cada dataset
anscombe_tidy |>
summarise(mean_x = mean(x), mean_y = mean(y),
var_x = var(x), var_y = var(y),
sd_x = sd(x), sd_y = sd(y),
cov = cov(x, y), cor = cor(x, y), .by = dataset)# A tibble: 4 × 9
dataset mean_x mean_y var_x var_y sd_x sd_y cov cor
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 1 9 7.50 11 4.13 3.32 2.03 5.50 0.816
2 2 9 7.50 11 4.13 3.32 2.03 5.5 0.816
3 3 9 7.5 11 4.12 3.32 2.03 5.50 0.816
4 4 9 7.50 11 4.12 3.32 2.03 5.50 0.817Si te fijas todos los datasets tienen los mismos momentos muestrales, y en concreto tienen la misma correlación y misma covarianza.
¿Serán el mismo dataset desordenado?
ggplot(anscombe_tidy) +
geom_point(aes(x = x, y = y, color = dataset), size = 3, alpha = 0.7) +
ggthemes::scale_color_colorblind() +
facet_wrap(~dataset) +
theme_minimal()
Por suerte o por desgracia no todo son matemáticas: antes de pensar que modelo es mejor para nuestros datos, es importantísimo realizar un análisis exploratorio de los datos (incluyendo, por supuesto, una correcta visualización).
6 Caso práctico II: datasauRus
Podemos visualizarlo de manera aún más extrema con el dataset datasaurus_dozen del paquete {datasauRus} (ver más en https://www.research.autodesk.com/publications/same-stats-different-graphs/)
library(datasauRus)
datasaurus_dozen |>
summarise(mean_x = mean(x), mean_y = mean(y),
var_x = var(x), var_y = var(y),
sd_x = sd(x), sd_y = sd(y),
cov = cov(x, y), cor = cor(x, y), .by = dataset)# A tibble: 13 × 9
dataset mean_x mean_y var_x var_y sd_x sd_y cov cor
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 dino 54.3 47.8 281. 726. 16.8 26.9 -29.1 -0.0645
2 away 54.3 47.8 281. 726. 16.8 26.9 -29.0 -0.0641
3 h_lines 54.3 47.8 281. 726. 16.8 26.9 -27.9 -0.0617
4 v_lines 54.3 47.8 281. 726. 16.8 26.9 -31.4 -0.0694
5 x_shape 54.3 47.8 281. 725. 16.8 26.9 -29.6 -0.0656
6 star 54.3 47.8 281. 725. 16.8 26.9 -28.4 -0.0630
7 high_lines 54.3 47.8 281. 726. 16.8 26.9 -30.9 -0.0685
8 dots 54.3 47.8 281. 725. 16.8 26.9 -27.2 -0.0603
9 circle 54.3 47.8 281. 725. 16.8 26.9 -30.8 -0.0683
10 bullseye 54.3 47.8 281. 726. 16.8 26.9 -31.0 -0.0686
11 slant_up 54.3 47.8 281. 726. 16.8 26.9 -31.0 -0.0686
12 slant_down 54.3 47.8 281. 726. 16.8 26.9 -31.2 -0.0690
13 wide_lines 54.3 47.8 281. 726. 16.8 26.9 -30.1 -0.0666De nuevo tenemos la misma correlación pero…¿tendrán la misma relación lineal de verdad?
ggplot(datasaurus_dozen |>
filter(dataset != "wide_lines"),
aes(x = x, y = y, color = dataset)) +
geom_point(size = 2.5, alpha = 0.75) +
geom_smooth(method = "lm", se = FALSE) +
facet_wrap(~dataset, ncol = 3) +
theme_minimal()`geom_smooth()` using formula = 'y ~ x'
Moraleja
Los resúmenes numéricos son eso, resúmenes, y como veremos más adelante, pueden existir distintos factores que enmascaran la verdadera relación de las variables. Por eso saber visualizar de una manera correcta es vital.
7 Caso práctico III: wine.csv
Vamos a poner en práctica lo aprendido con el dataset wine.csv disponible en el campus.
datos <- read_csv(file = "./datos/wine.csv")Rows: 27 Columns: 7
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
dbl (7): Year, Price, WinterRain, AGST, HarvestRain, Age, FrancePop
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.datos# A tibble: 27 × 7
Year Price WinterRain AGST HarvestRain Age FrancePop
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 1952 7.50 600 17.1 160 31 43184.
2 1953 8.04 690 16.7 80 30 43495.
3 1955 7.69 502 17.2 130 28 44218.
4 1957 6.98 420 16.1 110 26 45152.
5 1958 6.78 582 16.4 187 25 45654.
6 1959 8.08 485 17.5 187 24 46129.
7 1960 6.52 763 16.4 290 23 46584.
8 1961 8.49 830 17.3 38 22 47128.
9 1962 7.39 697 16.3 52 21 48089.
10 1963 6.71 608 15.7 155 20 48799.
# ℹ 17 more rowsEl conjunto de datos está formado por 27 observaciones (cosechas de vino rojo Burdeos) y 7 variables
Year,Age: año de la cosecha y número de años en barrica.Price: precio en 1990-1991 de cada cosecha (en escala log)WinterRain: lluvia (en mm) que cayó ese año en invierno.AGST: crecimiento medio de la temperatura (en grados Celsius) durante la temporada.HarvestRain: lluvia (en mm) que cayó ese año durante la cosecha.FrancePop: población (miles de habitantes) de Francia.
Puedes ver más en detalles en https://doi.org/10.1080/09332480.1995.10542468
Si el objetivo fuese predecir el precio, \(Y = precio\), con una sola predictora \(X\), ¿cuál tendríamos que elegir? Si el modelo lineal, parece ya claro que deberemos elegir a priori aquella con mayor correlación lineal (en valor absoluto). Pero como hemos visto con los ejemplos, a veces los resúmenes numéricos pueden engañarnos así que lo primero que haremos será realizar un análisis exploratorio.
7.1 Análisis exploratorio
Algunas de las primeras preguntas que podemos hacernos son:
- ¿Las variables son numéricas (continuas)?
- ¿Tienen problemas de rango (por ejemplo, pesos negativos)? ¿Tienen datos ausentes?
Para responder a dichas preguntas lo más sencillos es hacer uso de la función skim() del paquete {skimr}
library(skimr)
Attaching package: 'skimr'The following object is masked from 'package:corrr':
focusdatos |> skim()| Name | datos |
| Number of rows | 27 |
| Number of columns | 7 |
| _______________________ | |
| Column type frequency: | |
| numeric | 7 |
| ________________________ | |
| Group variables | None |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| Year | 0 | 1 | 1966.81 | 8.25 | 1952.00 | 1960.50 | 1967.00 | 1973.50 | 1980.00 | ▅▇▇▇▇ |
| Price | 0 | 1 | 7.04 | 0.63 | 6.20 | 6.51 | 6.98 | 7.44 | 8.49 | ▇▅▆▂▂ |
| WinterRain | 0 | 1 | 608.41 | 129.03 | 376.00 | 543.50 | 600.00 | 705.50 | 830.00 | ▃▃▇▃▃ |
| AGST | 0 | 1 | 16.48 | 0.66 | 14.98 | 16.15 | 16.42 | 17.01 | 17.65 | ▂▃▇▃▅ |
| HarvestRain | 0 | 1 | 144.81 | 73.07 | 38.00 | 88.00 | 123.00 | 185.50 | 292.00 | ▇▇▇▁▅ |
| Age | 0 | 1 | 16.19 | 8.25 | 3.00 | 9.50 | 16.00 | 22.50 | 31.00 | ▇▇▇▇▅ |
| FrancePop | 0 | 1 | 50085.44 | 3793.00 | 43183.57 | 46856.00 | 50650.41 | 53511.21 | 55110.24 | ▃▃▅▅▇ |
En este caso no tenemos ausentes ni problemas de codificación
Profundizando un poco podemos preguntarnos…
- ¿Cómo se distribuyen las variables?
- ¿Hay datos atípicos?
Para ello podemos empezar visualizando la distribución de cada variable haciendo uso de densidades, histogramas y/o boxplots (por ejemplo), pero si queremos ahorrar tiempo en visualizar convertimos nuestro dataset en tidydata
datos_tidy <-
datos |>
pivot_longer(cols = everything(), names_to = "variable", values_to = "values")
ggplot(datos_tidy, aes(x = values)) +
geom_density(aes(color = variable, fill = variable),
alpha = 0.75) +
ggthemes::scale_color_colorblind() +
ggthemes::scale_fill_colorblind() +
facet_wrap(~variable, scales = "free") +
theme_minimal()
ggplot(datos_tidy, aes(x = values)) +
geom_histogram(aes(color = variable, fill = variable),
alpha = 0.75, bins = 10) +
ggthemes::scale_color_colorblind() +
ggthemes::scale_fill_colorblind() +
facet_wrap(~variable, scales = "free") +
theme_minimal()
ggplot(datos_tidy, aes(y = values, x = variable)) +
geom_boxplot(aes(color = variable, fill = variable),
alpha = 0.75) +
ggthemes::scale_color_colorblind() +
ggthemes::scale_fill_colorblind() +
facet_wrap(~variable, scales = "free") +
theme_minimal()
No se aprecian valores atípicos (al menos respecto a los percentiles)
Podemos incluso gráficar todos los scatter plot sin transformar los datos con facet_matrix() del paquete {ggforce}
Code
library(ggforce)
ggplot(datos) +
geom_point(aes(x = .panel_x, y = .panel_y)) +
facet_matrix(vars(everything())) +
ggthemes::scale_color_colorblind() +
ggthemes::scale_fill_colorblind() +
theme_minimal()
Podemos también visualizar con un scatter plot todas las variables y además sus correlaciones, con el paquete {GGally} y la función `ggpairs(), sin necesidad de convertir a tidydata.
Code
library(GGally)Registered S3 method overwritten by 'GGally':
method from
+.gg ggplot2Code
ggpairs(datos) +
theme_minimal()
7.2 Análisis de dependencia
Tras una primera visualización de los datos, de cada variable por separado, nuestro objetivo real sería responder a
- ¿Qué predictora está más correlacionada (linealmente) con la variable objetivo a predecir? ¿Existe otro tipo de dependencia ?
- ¿Cómo se relacionan las predictoras entre sí? ¿Están correlacionadas? Este último paso será crucial en el contexto multivariante pero en este caso simplemente vamos a ver como se relacionan linealmente las predictoras entre sí (por saber como se haría), y cuál de ellas es la más adecuada para predecir linealmente
precio
El primer paso es la matriz de correlaciones con la función correlate() del paquete {corrr}
datos |>
correlate() |>
fashion()Correlation computed with
• Method: 'pearson'
• Missing treated using: 'pairwise.complete.obs' term Year Price WinterRain AGST HarvestRain Age FrancePop
1 Year -.46 .05 -.29 -.06 -1.00 .99
2 Price -.46 .13 .67 -.51 .46 -.48
3 WinterRain .05 .13 -.32 -.27 -.05 .03
4 AGST -.29 .67 -.32 -.03 .29 -.30
5 HarvestRain -.06 -.51 -.27 -.03 .06 -.03
6 Age -1.00 .46 -.05 .29 .06 -.99
7 FrancePop .99 -.48 .03 -.30 -.03 -.99 Podemos visualizar dichas correlaciones con un grafo con network_plot()
datos |>
correlate() |>
network_plot(min_cor = 0.5)Correlation computed with
• Method: 'pearson'
• Missing treated using: 'pairwise.complete.obs'
En cuando a la dependencia entre predictoras, las variables Age, Year y FrancePop presentan la misma información. Si queremos centrarnos solo en las variables vs la predictora, podemos hacer focus() en ella
datos |>
correlate() |>
corrr::focus(Price)Correlation computed with
• Method: 'pearson'
• Missing treated using: 'pairwise.complete.obs'# A tibble: 6 × 2
term Price
<chr> <dbl>
1 Year -0.460
2 WinterRain 0.135
3 AGST 0.668
4 HarvestRain -0.507
5 Age 0.460
6 FrancePop -0.481Respecto a Y, las predictoras con mayor correlación lineal son AGST (más calor, menos cosechas, sube el precio) y HarvestRain (más lluvias, más cosechas, baja el precio, ¡el signo importa!)
También podemos usar corrplot() del paquete {corrplot}, al que le pasamos una matriz de correlaciones clásica y nos la visualiza.
library(corrplot)
datos |>
cor() |>
corrplot(method = "ellipse")
Otra opción es visualizar con un scatter plots todas las predictoras vs Y, pivotando antes nuestro dataset (solo pivotamos las predictoras)
datos_tidy <-
datos |> pivot_longer(cols = -Price, names_to = "variable",
values_to = "values")
ggplot(datos_tidy, aes(x = values, y = Price)) +
geom_point(aes(color = variable), alpha = 0.7) +
geom_smooth(method = "lm") +
ggthemes::scale_color_colorblind() +
facet_wrap(~variable, scales = "free_x") +
theme_minimal()`geom_smooth()` using formula = 'y ~ x'
No solo comprobamos que las rectas con más pendiente son AGST y HarvestRain, además los puntos parecen poder ajustarse a una recta sin otro patrón identificable. Esto es importante hacerlo ya que debemos descartar posibles correlaciones espúreas (ver ejemplo datasaurus)
Por lo tanto parece claro: la mejor covariable para predecir linealmente el precio sería AGST (en sentido positivo) y HarvestRain (en sentido negativo)