9 Regresión
R pre-introductorio cree en las regresiones. Es más, las adora.
9.1 Correlación
En el capítulo anterior estudiamos medidas descriptivas de una variable.
Aprendimos que la varianza es la desviación cuadrada promedio de los valores respecto a la media. Al ser “cuadrada”, la varianza tiene una desventaja: no la podemos interpretar en las unidades originales en las que está planteada la variable.
Una solución al problema de la naturaleza cuadrada de la varianza es -you guessed it- sacarle la raíz cuadrada, con lo cual la varianza se convierte en la desviación estandar y nuevamente tenemos una cifra que se puede interpretar en las unidades originales de la variable… Necesito que retengamos esta diferencia entre varianza y desviación estandar.
Trabajaremos ahora no con una sino con dos variables al mismo tiempo. Cuando analizamos dos variables nuestro interés es determinar si están correlacionadas: si una cambia en cuanto cambia la otra.
Hay tres formas de medir la correlación y cada una brinda información más sofisticada que la anterior.
La primera opción para medir la correlación entre dos variables consiste en determinar la dirección en la que se mueven. La covarianza sirve para saber si las variables cambian en la misma dirección, en direcciones opuestas, o si el cambio de una no tiene nada que ver con el cambio de la otra.
La covarianza de dos variables es igual al (1) promedio del (2) producto de las distancias respecto a la media.
El (2) producto de las distancias respecto a la media se calcula así:
## Compute
dis.products <- (iris$Petal.Length - mean(iris$Petal.Length)) * (iris$Petal.Width - mean(iris$Petal.Width))
## Print
head(dis.products)
#> [1] 2.356428 2.356428 2.456361 2.256495 2.356428 1.645028El (1) promedio del (2) producto de las distancias se computa como sigue:
Y ahora vamos a calcular la covarianza mediante la función cov(), incluyendo descomposición de código para corrobar que entendemos lo que estamos haciendo:
## Compute
cov(iris$Petal.Length, iris$Petal.Width)
#> [1] 1.295609
## Ok, let's prove it
## Multiply distances
dis.products <- (iris$Petal.Length - mean(iris$Petal.Length)) * (iris$Petal.Width - mean(iris$Petal.Width))
## Add distance products and divide by N - 1
sum(dis.products) / (dim(iris)[1] - 1)
#> [1] 1.295609La covarianza es importante, pero no es demasiado informativa. La covarianza sólo nos revela el signo de la correlación de las variables, siendo el signo lo mismo que la dirección: la covarianza delata si dos variables se mueven en la misma dirección (signo positivo) o en direcciones opuestas (signo negativo):
Si las dos variables tienden a estar a la vez por encima (o por debajo) de su media, la covarianza es positiva.
Si una tiende a estar por encima de su media mientras la otra está por abajo, la covarianza es negativa.
Si las variables no están relacionadas linealmente, la covarianza es cero.
Este problema de la covarianza -que la cifra en sí misma no nos dice nada, sólo nos sirve el signo- se puede solventar si la dividimos por el producto de las desviaciones estandar. De hecho, si efectuamos esta división lo que obtenemos es una segunda medida de correlación: el coeficiente de correlación.
R puede calcular el coeficiente de correlación con la función cor():
## Get the correlation
cor(iris$Petal.Length, iris$Petal.Width)
#> [1] 0.9628654
## Show your understanding
cov(iris$Petal.Length, iris$Petal.Width) /
(sd(iris$Petal.Length) * sd(iris$Petal.Width))
#> [1] 0.9628654El coeficiente de correlación es más informativo que la covarianza, un paso hacia adelante. La covarianza -recuerden- sólo nos indica la dirección de la correlación (la cifra en sí misma carece de significado). Al dividir la covarianza entre el producto de las desviaciones estandar lo que conseguimos es expresarla en un rango que va de -1 a 1; a esa covarianza “normalizada” se le conoce como coeficiente de correlación y es tremendo avance porque entonces nos informa sobre la dirección y la magnitud de la relación entre dos variables.
Comparen la covarianza y el coeficiente de correlación:
## Explore covariance
cov(na.omit(penguins[, 5:6]))
#> flipper_length_mm body_mass_g
#> flipper_length_mm 197.7318 9824.416
#> body_mass_g 9824.4161 643131.077
## Explore correlation
cor(na.omit(penguins[, 5:6]))
#> flipper_length_mm body_mass_g
#> flipper_length_mm 1.0000000 0.8712018
#> body_mass_g 0.8712018 1.0000000Observen cómo los coeficientes de correlación están normalizados: en tanto sólo pueden ir de -1 (los negativos más fuertes) a 1 (los positivos más fuertes), pasando por 0 (si no hay correlación), el coeficiente de correlación nos revela tanto la dirección (positiva o negativa) como la magnitud (fuerte o debil) de la correlación entre dos variables.
De ahí que a la gente le gusten estos gráficos:
## Prepare correlation matrix
cors <- cor(na.omit(penguins[, 3:6]))
## Plot
corrplot(cors, method = "number")
La verdad, todo esto está muy fácil…
La tercera forma de estimar la correlación empieza siendo muy visual. Vamos a usar el data set iris. Abajo, podemos apreciar que las dos variables que grafiqué presentan una relación lineal positiva. Lo sabemos porque los puntos van juntitos hacia arriba:
## Plot
plot(y = iris$Petal.Length, x = iris$Petal.Width) 
Y si quiero hasta le pongo colores (en función de la variable Species):
## Set colors
colors <- c("red", "gold", "blue")
## Plot
plot(y = iris$Petal.Length,
x = iris$Petal.Width,
col = colors[iris$Species]) 
Fasten your seatbelts, vienen temas serios.
El gráfico arriba es clarísimo sobre la correlación positiva que existe entre las dos variables, pero no nos da un número con el que la podamos expresar.
Para obtener esa medición primero tendremos que pedirle a R que dibuje una línea encima de los puntos. Pero no es cualquier línea la que nos interesa: queremos exactamente la línea que, en promedio, pasa más cerca de todos los puntos.
Efectivamente, R es capaz de encontrar la línea que mejor se ajusta a los puntos:
## Plot
plot(y = iris$Petal.Length,
x = iris$Petal.Width)
abline(lm(iris$Petal.Length ~ iris$Petal.Width), col = "blue")
La línea azul es la línea que, en promedio, pasa más cerca de todos los puntos. Esa línea tiene una pendiente, es decir, un grado de inclinación. Junto con la covarianza y el coeficiente de correlación, la pendiente de la línea azul es una medida de la correlación entre la variable en el eje X y la variable en el eje Y: la pendiente de la línea de regresión.
¿Cómo hace R para saber que la línea azul es la que mejor se ajusta, en promedio, a los datos? ¿Cómo la obtiene? Primero, es importante saber que esta línea azul resulta de un modelo de regresión lineal: es, en otras palabras, una línea de regresión.
Pronto vamos a explicar por qué la pendiente de la línea de regresión es mejor que cualquier otra medida de correlación.
9.2 Regresión lineal
El gráfico anterior expresa visualmente la corelación entre una variable \(Y\) y una variable \(X\). La mejor forma de describir conceptualmente la correlación entre \(Y\) y \(X\) es mediante un modelo de regresión lineal, que a su vez se expresa mediante la siguiente fórmula:
\[ Y = (constante + X * pendiente) + residuo \]
Es fácil si vamos por partes:
En la fórmula, \(Y\) es la variable dependiente y \(X\) es la variable independiente; justamente lo que queremos determinar es cómo el cambio en \(X\) se manifesta en \(Y\). Normalmente, como es nuestro caso ahora, uno conoce los valores de ambas variables pues están en el conjunto de datos con el que trabajamos.
En cambio, la \(constante\) y la \(pendiente\) no las conocemos. Es más, el fin del modelo de regresión es precisamente estimar estos dos coeficientes, que es justo lo que logra R con la función
lm(): R nos devuelve el coeficiente de la \(constante\) (intercept) y el coeficiente de la \(pendiente\) (este es el que más nos interesa).El \(residuo\) es… bueno, de él mejor hablemos en un rato. Una pista: de momento he dicho que el modelo de regresión lineal calcula los coeficientes de la \(constante\) y la \(pendiente\) pero no he dicho cómo… R escoge los coeficientes que garantizan el menor \(residuo\) posible.
En R, podemos correr un modelo de regresión lineal y extraer los coeficientes de la \(constante\) y la \(pendiente\) con impactante facilidad:
## Model
model <- lm(iris$Petal.Length ~ iris$Petal.Width) # Y ~ X
## Get coefficients
coef(model)
#> (Intercept) iris$Petal.Width
#> 1.083558 2.229940Arriba tenemos los coeficientes de la \(constante\) y la \(pendiente\) respectivamente.
Permítanme reacomodar la fórmula un poco:
\[ Y - residuo = (constante + X * pendiente) \]
Para demostrar cómo opera esta ecuación, tomemos solamente un par de filas del conjunto de datos iris:
## Subset
(iris_una_fila <- iris[1, 3:4])
#> Petal.Length Petal.Width
#> 1 1.4 0.2
(iris_otra_fila <- iris[130, 3:4])
#> Petal.Length Petal.Width
#> 130 5.8 1.6En términos de la fórmula del modelo de regresión lineal, Petal.Length es \(Y\) mientras Petal.Width es \(X\). O sea que ya de entrada tenemos \(Y\) y \(X\).
Ahora recortemos un pedacito de la ecuación. Trabajemos sólo con la expresión que está a la derecha, entre paréntesis:
\[ constante + X * pendiente \]
Si se fijan, ya contamos con todos los elementos para calcular la expresión \(constante + X * pendiente\), y pues eso mismo haremos:
## Remove names
coefs <- coef(model)
names(coefs) <- NULL
## Set terms
(Y1 = iris_una_fila[1] |> pull())
#> [1] 1.4
(X1 = iris_una_fila[2] |> pull())
#> [1] 0.2
(Y2 = iris_otra_fila[1] |> pull())
#> [1] 5.8
(X2 = iris_otra_fila[2] |> pull())
#> [1] 1.6
(constante = coefs[1])
#> [1] 1.083558
(pendiente = coefs[2])
#> [1] 2.22994
## Compute
constante + X1 * pendiente
#> [1] 1.529546
constante + X2 * pendiente
#> [1] 4.651463Los resultados que acabamos de obtener para la expresión \(constante + X * pendiente\) son los valores que el modelo de regresión lineal predice para (dos) valores específicos de \(Y\) dados (dos) valores específicos de \(X\).
Voy a agregar dichas predicciones de \(Y\) como una columna nueva en iris_una_fila y iris_otra_fila, a la que voy a llamar \(yPredicted\):
iris_una_fila$yPredicted <- constante + X1 * pendiente
iris_otra_fila$yPredicted <- constante + X2 * pendienteY ahora junto los data sets con la función rbind():
## Bind
(iris_dos_filas <- rbind(iris_una_fila, iris_otra_fila))
#> Petal.Length Petal.Width yPredicted
#> 1 1.4 0.2 1.529546
#> 130 5.8 1.6 4.651463Recordemos que Petal.Length es el valor real de \(Y\) y Petal.Width es el valor real de \(X\). Por su lado, \(yPredicted\) es el valor que el modelo de regresión lineal predice para \(Y\) en función de \(X\).
\(yPredicted\) resulta de haber llevado a término la expresión \(constante + X * pendiente\) que, a su vez, requirió obtener la \(constante\) y la \(pendiente\) mediante la función lm().
O sea:
\[ yPredicted = (constante + X * pendiente) \]
Por lo tanto, esta ecuación:
\[ Y - residuo = (constante + X * pendiente) \]
Es igual que esta otra:
\[ Y - residuo = yPredicted \]
Por cierto, no hace falta calcular \(yPredicted\) “a mano”; si lo hicimos así es estrictamente para entender bien el concepto. Del mismo modo que la función lm() produce los coeficientes \(constante\) y \(pendiente\), ella nos da también los valores \(yPredicted\), los cuales se pueden accesar así:
## Get fitted values
model$fitted.values |> head()
#> 1 2 3 4 5 6
#> 1.529546 1.529546 1.529546 1.529546 1.529546 1.975534Si estamos desconfiados de los valores yPredicted que tenemos en iris_dos_filas porque los calculamos “a mano”, los podemos comparar con los valores que arrojó el modelo de regresión lineal. Es tan fácil como incluirlos en el conjunto de datos original, iris:
## Create variable
iris$yPredicted <- model$fitted.values
## Subset
iris[1, c(3:4, 6)]
#> Petal.Length Petal.Width yPredicted
#> 1 1.4 0.2 1.529546
iris[130, c(3:4, 6)]
#> Petal.Length Petal.Width yPredicted
#> 130 5.8 1.6 4.651463
## Print
iris_dos_filas
#> Petal.Length Petal.Width yPredicted
#> 1 1.4 0.2 1.529546
#> 130 5.8 1.6 4.651463Así, hemos comprobado que el resultado que obtuvimos vía descomposición de código es consistente con el que R consiguió.
Intentaré destacar los dos puntos de iris_dos_filas en un scater plot para explicarlo todo de nuevo, pero más fácil esta vez:
## Plot
ggplot(iris, aes(x = Petal.Width, y = Petal.Length)) +
geom_point(color = "grey90") +
geom_point(aes(x = iris_dos_filas$Petal.Width[1],
y = iris_dos_filas$Petal.Length[1]),
color = "red") +
geom_point(aes(x = iris_dos_filas$Petal.Width[2],
y = iris_dos_filas$Petal.Length[2]),
color = "red") +
theme_light()
Marqué todos los datos de iris en un gris muy ténue y resalté en rojo sólo los dos puntos que corresponden a iris_dos_filas. Atención: los puntos se localizan en función de \(Y\), que está en el eje Y, y de \(X\), que está en el eje X. El valor en el eje Y que estamos observando es el valor \(Y\) real, esto es, el valor observado en el conjunto de datos.
Recordemos que los modelos de regresión lineal nos permiten predecir el valor de \(Y\) en función del valor \(X\). Dijimos antes que el resultado de esa predicción es \(yPredicted\).
Entonces voy a incluir \(yPredicted\) en la visualización:
## Plot
ggplot(iris, aes(x = Petal.Width, y = Petal.Length)) +
geom_point(color = "white") +
geom_point(aes(x = iris_dos_filas$Petal.Width[1],
y = iris_dos_filas$Petal.Length[1]),
color = "red") +
geom_point(aes(x = iris_dos_filas$Petal.Width[2],
y = iris_dos_filas$Petal.Length[2]),
color = "red") +
geom_point(aes(x = iris_dos_filas$Petal.Width[1],
y = iris_dos_filas$yPredicted[1]),
shape = 21,
colour = "black",
fill = "white",
size = 1.5,
stroke = 1) +
geom_point(aes(x = iris_dos_filas$Petal.Width[2],
y = iris_dos_filas$yPredicted[2]),
shape = 21,
colour = "black",
fill = "white",
size = 1.5,
stroke = 1) +
theme_light()
El valor que predice el modelo de regresión lineal (\(yPredicted\)), aunque intenta ser lo más próximo posible al valor real de \(Y\), casi nunca es igual. En los casos arriba, se observa que tenemos una predicción en la que el modelo sobreestimó el valor real -y casi lo pega- y otra en la que el modelo lo subestimó -y no estuvo tan cerca-.
Prometí que esta variedad de medición de la correlación entre dos variables sería muy visual. Había graficado los datos y luego le pedí a R que les encajara encima la línea que, en promedio, pasa más cerca de todos los puntos.
Voy a agregar esa línea, la línea de regresión, a este mismo gráfico:
## Plot
ggplot(iris, aes(x = Petal.Width, y = Petal.Length)) +
geom_point(color = "white") +
geom_point(aes(x = iris_dos_filas$Petal.Width[1],
y = iris_dos_filas$Petal.Length[1]),
color = "red") +
geom_point(aes(x = iris_dos_filas$Petal.Width[2],
y = iris_dos_filas$Petal.Length[2]),
color = "red") +
geom_point(aes(x = iris_dos_filas$Petal.Width[1],
y = iris_dos_filas$yPredicted[1]),
shape = 21,
colour = "black",
fill = "white",
size = 1.5,
stroke = 1) +
geom_point(aes(x = iris_dos_filas$Petal.Width[2],
y = iris_dos_filas$yPredicted[2]),
shape = 21,
colour = "black",
fill = "white",
size = 1.5,
stroke = 1) +
geom_smooth(method = "lm", se = F) +
theme_light()
#> `geom_smooth()` using formula = 'y ~ x'
Ok, acabamos de aprender algo re importante: la línea de regresión es, en realidad, la línea que pasa por los valores \(yPredicted\)… en otras palabras, los valores que el modelo predice para \(Y\) en función de los valores de \(X\) están todos reunidos en la línea de regresión.
Ya hemos traído a esta visualización todos los ingredientes de la fórmula… o casi todos: nos faltó el \(residuo\).
El \(residuo\) es justamente la distancia que hay entre el valor real de \(Y\) (el punto rojo) y el valor que predijo el modelo, \(yPredicted\) (el que está en la línea de regresión). Esto se aprecia claramente en la ecuación:
\[ Y - residuo = yPredicted \]
Si despejamos el \(residuo\), tenemos que:
\[ residuo = Y - yPredicted \]
El \(residuo\) es la distancia entre \(Y\) y \(yPredicted\). Voy a visualizarlos para salir de cualquier duda:
## Plot
ggplot(iris, aes(x = Petal.Width, y = Petal.Length)) +
geom_point(color = "white") +
geom_point(aes(x = iris_dos_filas$Petal.Width[1],
y = iris_dos_filas$Petal.Length[1]),
color = "red") +
geom_point(aes(x = iris_dos_filas$Petal.Width[2],
y = iris_dos_filas$Petal.Length[2]),
color = "red") +
geom_point(aes(x = iris_dos_filas$Petal.Width[1],
y = iris_dos_filas$yPredicted[1]),
shape = 21,
colour = "black",
fill = "white",
size = 1.5,
stroke = 1) +
geom_point(aes(x = iris_dos_filas$Petal.Width[2],
y = iris_dos_filas$yPredicted[2]),
shape = 21,
colour = "black",
fill = "white",
size = 1.5,
stroke = 1) +
geom_segment(data = iris |> filter(
Petal.Width == iris_dos_filas$Petal.Width & Petal.Length == iris_dos_filas$Petal.Length),
aes(xend = Petal.Width,
yend = yPredicted),
color = "red") +
geom_smooth(method = "lm", se = F) +
theme_light()
#> `geom_smooth()` using formula = 'y ~ x'
El \(residuo\) es lo que le faltó al valor que el modelo predice (\(yPredicted\)) para ser igual al valor real (\(Y\)). Por fin la ecuación cobra sentido.
Habíamos empezado con esta fórmula:
\[ Y = (constante + X * pendiente) + residuo \]
Aprendimos que \(Y\) y \(X\) las tenemos de entrada en nuestros datos, mientras R nos permite estimar la \(constante\) y la \(pendiente\).
Aprendimos que la expresión entre paréntesis es, de hecho, el valor que el modelo predice para \(Y\) dado \(X\):
\[ Y = yPredicted + residuo \]
Y aprendimos que el \(residuo\) es la diferencia entre \(Y\) y \(yPredicted\):
\[ Y - yPredicted = residuo \]
Este ejercicio, si lo repetimos para todas las observaciones, resulta en algo así:
## Create variable
iris$RESIDUALS <- iris$Petal.Length - iris$yPredicted
## Print
car::some(iris[, c(3:4, 6:7)])
#> Petal.Length Petal.Width yPredicted RESIDUALS
#> 49 1.5 0.2 1.529546 -0.02954613
#> 70 3.9 1.1 3.536493 0.36350742
#> 77 4.8 1.4 4.205475 0.59452527
#> 92 4.6 1.4 4.205475 0.39452527
#> 106 6.6 2.1 5.766433 0.83356693
#> 112 5.3 1.9 5.320445 -0.02044497
#> 116 5.3 2.3 6.212421 -0.91242117
#> 121 5.7 2.3 6.212421 -0.51242117
#> 133 5.6 2.2 5.989427 -0.38942712
#> 140 5.4 2.1 5.766433 -0.36643307Un resultado que puedo comprobar con el resulto que arrojó la función lm():
## Create variable
iris$residuals <- model$residuals
## Print
car::some(iris[, c(3:4, 7:8)])
#> Petal.Length Petal.Width RESIDUALS residuals
#> 31 1.6 0.2 0.07045387 0.07045387
#> 32 1.5 0.4 -0.47553423 -0.47553423
#> 33 1.5 0.1 0.19344792 0.19344792
#> 41 1.3 0.3 -0.45254018 -0.45254018
#> 54 4.0 1.3 0.01751932 0.01751932
#> 56 4.5 1.3 0.51751932 0.51751932
#> 57 4.7 1.6 0.04853717 0.04853717
#> 62 4.2 1.5 -0.22846878 -0.22846878
#> 63 4.0 1.0 0.68650147 0.68650147
#> 66 4.4 1.4 0.19452527 0.19452527Al final tenemos un \(residuo\) para cada observación. ¿Qué pasa si los sumo todos juntos?
No los puedo sumar directamente; el resultado sería -básicamente- cero ya que estas diferencias entre \(Y\) y \(yPredicted\) se cancelan entre sí.
Si los quiero sumar los debo elevar al cuadrado primero:
## Compute
sum(iris$RESIDUALS ^ 2)
#> [1] 33.84475El resultado anterior es la suma del cuadrado de los errores (SSE: Sum of Squares Error). La función lm() consigue los coeficientes de la \(constante\) y la \(pendiente\) que minimizan la suma del cuadrado de los errores, de ahí que la línea que finalmente dibujamos con base en esos coeficientes y los valores \(X\) y \(Y\) sea la línea que, en promedio, pasa más cerca de todos los datos.
En R, se calcularía así también:
Ahora que ya todo tiene sentido, quiero que conectemos la \(pendiente\) de la línea de regresión con la covarianza y el coeficiente de correlación que estudiamos atrás:
## Get the slope coefficient
coef(model)[2]
#> iris$Petal.Width
#> 2.22994El coeficiente de la \(pendiente\) es igual a la covarianza de \(X\) y \(Y\) entre la varianza de \(X\):
Que es lo mismo que la correlación entre \(X\) y \(Y\) dividida entre el producto de las desviaciones estandar de \(X\) y \(Y\):
## Compute
cor(iris$Petal.Width, iris$Petal.Length) /
sd(iris$Petal.Width) * sd(iris$Petal.Length)
#> [1] 2.22994La \(pendiente\) de la línea de regresión es la mejor de las tres mediciones de correlación por la siguiente razón:
La covarianza sólo nos da la dirección del cambio.
El coeficiente de correlación nos da la dirección y la magnitud del cambio.
La \(pendiente\) de la línea de regresión nos da la dirección, la magnitud y el tamaño del cambio.
## Get the slope coefficient
coef(model)[2]
#> iris$Petal.Width
#> 2.22994La \(pendiente\) de la línea de regresión nos revela el cambio en \(Y\) por una unidad de cambio en \(X\): es el tamaño exacto del cambio.
Esto ha sido muy fácil. Dos precauciones, eso sí:
Correlación no es causalidad: para poder darle una interpretación causal a los coeficientes de una regresión se han de cumplir supuestos demasiado exigentes, una materia lejos de los alcances de R pre-introductorio (pero sí le recomiendo unos buenos libros).
Causalidad no es correlación: las tres medidas de correlación que estudiamos hasta ahora (la covarianza, el coeficiente de correlación, la \(pendiente\) de la línea de regresión) sirven exclusivamente para determinar si existe una relación lineal entre dos variables. No obstante, hay variables que pueden tener otro tipo de relación funcional -cuadrática, por ejemplo- y estarían relacionadas entre sí aunque den negativo a un test de correlación.
9.3 Wrapping-up
Las temidas regresiones son muy fáciles en R.
Podemos empezar por visualizarlas con un scatter plot, como este que voy a hacer con el data frame penguins:
## Plot
ggplot(penguins, aes(x = flipper_length_mm, y = body_mass_g)) +
geom_point(aes(color = species), shape = 4) +
geom_smooth(method = "lm", se = FALSE, color = "black")
#> `geom_smooth()` using formula = 'y ~ x'
#> Warning: Removed 2 rows containing non-finite values
#> (`stat_smooth()`).
#> Warning: Removed 2 rows containing missing values
#> (`geom_point()`).
Esa línea negra es una línea producida por un modelo de regresión lineal, lo que quiere decir que es la línea que, en promedio, pasa más cerca de todos los puntos.
Ahora le cambio las variables al modelo y en este nuevo me encuentro una correlación negativa:
## Plot
ggplot(penguins, aes(x = bill_depth_mm, y = bill_length_mm)) +
geom_point(shape = 5) +
geom_smooth(method = "lm", se = FALSE, color = "purple")
#> `geom_smooth()` using formula = 'y ~ x'
#> Warning: Removed 2 rows containing non-finite values
#> (`stat_smooth()`).
#> Warning: Removed 2 rows containing missing values
#> (`geom_point()`).
En la etapa de análisis exploratorio de datos es aconsejable agruparlos de diversas formas para ver si encontramos patrones relevantes.
Miren qué fenómeno interesante sucede con estas mismas variables cuando agrupo los pinguinos por su especie y estimo líneas de regresión para cada grupo por separado:
## Plot
ggplot(penguins, aes(x = bill_depth_mm, y = bill_length_mm)) +
geom_point(aes(color = species), shape = 5) +
geom_smooth(method = "lm", se = FALSE, color = "purple") +
geom_smooth(data = penguins |> filter(species == "Adelie"),
method = "lm",
se = FALSE,
color = "black",
linetype = "dashed") +
geom_smooth(data = penguins |> filter(species == "Chinstrap"),
method = "lm",
se = FALSE,
color = "black",
linetype = "dashed") +
geom_smooth(data = penguins |> filter(species == "Gentoo"),
method = "lm",
se = FALSE,
color = "black",
linetype = "dashed")
#> `geom_smooth()` using formula = 'y ~ x'
#> Warning: Removed 2 rows containing non-finite values
#> (`stat_smooth()`).
#> `geom_smooth()` using formula = 'y ~ x'
#> Warning: Removed 1 rows containing non-finite values
#> (`stat_smooth()`).
#> `geom_smooth()` using formula = 'y ~ x'
#> `geom_smooth()` using formula = 'y ~ x'
#> Warning: Removed 1 rows containing non-finite values
#> (`stat_smooth()`).
#> Warning: Removed 2 rows containing missing values
#> (`geom_point()`).
Resultó que la relación dentro de cada grupo es exactamente la contraria a la que rige para todo el conjunto (si esto les causa inquietud, ver: Simpson’s Paradox). Porque el mundo del análisis de datos está lleno de sorpresas así, siempre hemos de ser exigentes y puntillosos en nuestros abordajes.
Recordemos que la línea de regresión es la línea que, en promedio, pasa más cerca de todos los puntos. ¿Pero qué es esa línea exactamente? Es la línea que reúne los valores que el modelo predice para \(Y\) en función de \(X\).
Normalmente el valor real \(Y\) y el valor que predice el modelo (\(yPredicted\)) no coinciden, lo que nos lleva al concepto de \(residuo\). La diferencia entre \(Y\) y \(yPredicted\) es el \(residuo\).
Con base en las variables \(X\) y \(Y\), el modelo de regresión lineal encuentra los coeficientes de la \(constante\) y la \(pendiente\) que aseguran que ese \(residuo\) será, en promedio, el menor para todos los puntos, y esa es la garantía de que la línea de regresión es la que mejor se ajusta a los datos.
Todo lo anterior obedece a la siguiente ecuación:
\[ Y = (constante + X * pendiente) + residuo \]
En el caso de las variables flipper_length_mm y body_mass_g, los coeficientes de la \(constante\) y la \(pendiente\) son los siguientes:
## Prepare data
df <- penguins |>
select(flipper_length_mm, body_mass_g) |>
rename(y = "body_mass_g",
x = "flipper_length_mm")
## Model
(model <- lm(y ~ x, df))
#>
#> Call:
#> lm(formula = y ~ x, data = df)
#>
#> Coefficients:
#> (Intercept) x
#> -5780.83 49.69Con los coeficientes de la \(constante\) y la \(pendiente\) podemos llevar a término la expresión \(constante + X * pendiente\) y obtendremos como resultado una predicción del valor de \(Y\); a esa predicción le llamamos \(yPredicted\).
En la siguiente visualización, los puntos negros son los valores reales, corresponden a la intersección de \(Y\) y \(X\)), y las predicciones \(yPredicted\) quedan todas en la línea azul:
Una vez más: las predicciones de \(Y\), que aquí hemos llamado \(yPredicted\), son los valores que quedan ubicados justo en la línea de regresión. Si nos lo llevamos todo a un data frame, se vería así:
## Decompose model
df <- df |>
add_predictions(model) |>
rename(
yPredicted = "pred")
## Explore data
tail(df)
#> # A tibble: 6 × 3
#> x y yPredicted
#> <int> <int> <dbl>
#> 1 195 3650 3908.
#> 2 207 4000 4504.
#> 3 202 3400 4256.
#> 4 193 3775 3808.
#> 5 210 4100 4653.
#> 6 198 3775 4057.Nótese, sin embargo, que el valor \(yPredicted\) -casi- nunca es el mismo que el valor \(Y\) real. A tal diferencia le llamamos \(residuo\). Esa diferencia es la que el modelo ha minimizado lo más que ha podido; esto quiere decir que no hay ninguna otra línea que, en promedio, pase más cerca de los puntos:
## Compute
df$residual <- df$y - df$yPredicted
## Inspect
tail(df)
#> # A tibble: 6 × 4
#> x y yPredicted residual
#> <int> <int> <dbl> <dbl>
#> 1 195 3650 3908. -258.
#> 2 207 4000 4504. -504.
#> 3 202 3400 4256. -856.
#> 4 193 3775 3808. -33.5
#> 5 210 4100 4653. -553.
#> 6 198 3775 4057. -282.El \(residuo\) se puede visualizar también:
## Plot
ggplot(penguins, aes(x = flipper_length_mm, y = body_mass_g)) +
geom_point(aes(color = species), shape = 15) +
geom_smooth(method = "lm", se = FALSE, color = "black") +
geom_segment(aes(xend = df$x,
yend = df$yPredicted,
color = species))
#> `geom_smooth()` using formula = 'y ~ x'
#> Warning: Removed 2 rows containing non-finite values
#> (`stat_smooth()`).
#> Warning: Removed 2 rows containing missing values
#> (`geom_point()`).
#> Warning: Removed 2 rows containing missing values
#> (`geom_segment()`).
El gráfico de arriba muestra todos los elementos:
Cada cuadrado representa un valor real: la intersección entre \(X\) y \(Y\).
La línea de regresión negra representa los valores que el modelo predice, \(yPredicted\).
Las líneas verticales de colores representan el \(residuo\), es decir, la diferencia entre el valor real, \(Y\), y el valor que el modelo predijo, \(yPredicted\). De hecho, el modelo de regresión lineal consiguió los coeficientes de la \(constante\) y la \(pendiente\) de modo tal que si uno eleva cada \(residuo\) al cuadrado y los suma todos, el resultado de esa suma es el mínimo posible (de ahí que la línea de regresión sea la que mejor se ajusta a los datos).
Es muy interesante observar que, de hecho, la línea de regresión pasa justo por la intersección de las medias de ambas variables. Produciré una visualización menos cargada para que esto se aprecie mejor:
## Plot x and y
plot(y = df$y, x = df$x)
abline(lm(df$y ~ df$x, data = df), col = "blue")
abline(v = mean(df$x, na.rm = TRUE), col = "red", lty = "dashed")
abline(h = mean(df$y, na.rm = TRUE), col = "red", lty = "dashed")
Como sea, lo más importante de los coeficientes de regresión es su interpretación:
## Model
lm(penguins$body_mass_g ~ penguins$flipper_length_mm) # Y ~ X
#>
#> Call:
#> lm(formula = penguins$body_mass_g ~ penguins$flipper_length_mm)
#>
#> Coefficients:
#> (Intercept) penguins$flipper_length_mm
#> -5780.83 49.69Si la aleta gana un milímetro de largo, los pinguinos ganan 50 gramos (\(pendiente\)) de masa corporal. El intercept (\(constante\)) es el valor que tendría \(Y\) si \(X\) fuese igual a cero (si el pingüino no tuviese aletas, pongamos); como en esta situación, a veces el intercept no tiene interpretación significativa (un peso negativo no es posible).
Esta indicación del tamaño aproximado del cambio de \(X\) en el cambio de \(Y\) no la proporcionan ni el coeficiente de correlación (que únicamente da la dirección y la magnitud del cambio, expresada en una escala de -1 a 1) ni la covarianza (que sólo da la dirección del cambio).
Las mediciones de correlación lo que ponen a prueba es la existencia de relación lineal entre las variables; las variables, sin embargo, pueden tener relaciones funcionales que no son lineales, con lo cual sería incorrecto afirmar que dos variables no están relacionadas entre sí si sólo tenemos evidencia de que no lo están linealmente.
Otro día explico más. Pasemos a otra cosa por ahora.