Este texto se basa en los siguientes materiales:

Introducción

La regresión lineal es una técnica estadística muy poderosa. Muches estamos algo familiarizades con los modelos de regresión solo por leer las noticias. Hemos visto alguna vez un gráfico de dispersión. Los modelos lineales se pueden utilizar para varias tareas. Probablemente, las dos más importantes sean

  • evaluar si existe una relación lineal entre una variable numérica en el eje horizontal y el promedio de la variable numérica en el eje vertical
  • hacer predicciones (seguramente vamos a discutir bastante alrededor de este término) de la variable dependiente

De ajustes, residuos y correlacion

Al considerar la regresión lineal, es útil pensar en el proceso de ajuste de línea. ¿Qué signifca “ajustar” una recta a una nube de puntos? Vamos a definir la forma de un modelo lineal. Luego, vamos a tratar de explorar los criterios para definir que un ajuste es “bueno”. Y vamos a introducir una nueva métrica llamda “correlación lineal”.

Ajustando una línea a los datos

El gráfico siguiente muestra dos variables. Podemos modelar esa relación entre ambas de forma perfecta con una línea recta. La ecuación de esa recta es

\[y = 5 + 64.96 X\]

¿Qué quiere decir esa idea de relación perfecta? Básicamente que podemos predecir sin error el valor de \(y\) conociendo el valor de \(X\). Obviamente, esto es un ejemplo de juguete y estas relaciones perfectas no existen práctiamente en ninguna situación de la vida real. Por ejemplo, si tomamos el ingreso familiar de una persona (\(X\)) este valor seguramente proporcione una buena cantidad de información útil sobre el monto que un banco le ofrecerá a esta persona si quisiera acceder a cierto tipo de instrumentos crediticios. Sin embargo, si quisiéramos hacer una predicción del monto prestado en base \(X\) esa predicción no va a ser perfecta.

Pensemos… ¿por qué?

La regresión lineal es un método estadístico que se usa para ajustar una línea a los datos donde la relación entre dos variables, \(X\) e \(y\) puede ser modelada como si fuera una línea recta tolerando algún grado de error.

\[ y = \beta{0} + \beta{1}X + \epsilon\] - \(\beta{0}\) representa el intercepto o la ordenada al origen. ¿Cómo podemos interpretar este valor? Como el valor que asume \(y\) cuando \(X = 0\).

  • \(\beta_{1}\) es la pendiente de la recta, es decir, cuánto aumenta \(y\) cuando \(X\) aumenta una unidad

  • \(\epsilon\) representa el error está representado por

En general, estos valores se calculan en base a los datos. Si los datos observados son una muestra aleatoria de una población objetivo sobre la que estamos interesados en hacer inferencias, estos valores se consideran estimaciones puntuales para los parámetros de la población. Más adelante vamos a trabajar sobre cómo hacer inferencias sobre los parámetros de un modelo lineal.

Cuando usamos \(X\) para predecir \(y\), solemos llamarla variable predictora y a \(y\) variable dependiente, output o variable resultado. Es raro que todos los datos caigan perfectamente en una línea recta. En cambio, es más común que los datos aparezcan como una nube de puntos, como los ejemplos que se muestran en el gráfico siguiente.

En cada caso, los datos caen alrededor de una línea recta, incluso si ninguna de las observaciones cae exactamente sobre la línea. La primera gráfica muestra una relación lineal descendente relativamente fuerte, donde la variabilidad restante en los datos alrededor de la línea es menor en relación con la fuerza de la relación entre \(X\) y \(y\).

La segunda gráfica muestra una relación hacia arriba (positiva) que, aunque evidente, no es tan fuerte como la primera.

El último gráfico muestra una relación negativa (a la baja) muy débil en los datos, tan leve que apenas podemos notarlo.

En cada gráfico podríamos preguntarnos, ¿deberíamos mover la línea un poco hacia arriba o hacia abajo? ¿Deberíamos inclinarla más o menos? La primera pregunta nos remite a \(\beta_{0}\). La segunda, a \(\beta_{1}\). Dado que cada una controla un aspecto de una recta, entonces, estas cuestiones nos plantean la existencia de una cierta incertidumbre con respecto a nuestras estimaciones de los parámetros del modelo. Vamos a ir avanzando sobre esto.

También hay casos en los que ajustar una línea recta a los datos, incluso si existe una relación clara entre las variables, no es útil. No funciona. Uno de estos casos se muestra en la figura siguiente.

Se vé una relación perfecta entre las variables, aunque la misma no es lineal. ¿Qué tipo de relación se les ocurre que podría funcionar aquí? Por ahora vamos a centrarnos en los modelos lineales.

Usando regresión lineal para predecir el peso de una población

Vamos a trabajar con un dataset que proviene del siguiente repositorio. Son datos censales parciales de un área de Botswana y, específicamente, de la etnia !Kung San. Son datos compilados a partir de e entrevistas realizadas por Nancy Howell a fines de la década de 1960. Los datasets completos (para les antropólogues que les interese) pueden ser recuperados desde acá.

Vamos a trabajar con un recorte de estos datos, particularmente, con datos de edad, sexo, altura y peso.

> df <- read_delim('https://raw.githubusercontent.com/rmcelreath/rethinking/master/data/Howell1.csv', delim=";")
> 
> df <- df %>%
+   mutate(male = as.factor(case_when(
+           male == 0 ~ 'No',
+           TRUE ~ 'Yes'
+   )))

Hagamos un breve análisis descriptivo de estas cuatro variables.

> summary(df)
##      height           weight            age         male    
##  Min.   : 53.98   Min.   : 4.252   Min.   : 0.00   No :287  
##  1st Qu.:125.09   1st Qu.:22.008   1st Qu.:12.00   Yes:257  
##  Median :148.59   Median :40.058   Median :27.00            
##  Mean   :138.26   Mean   :35.611   Mean   :29.34            
##  3rd Qu.:157.48   3rd Qu.:47.209   3rd Qu.:43.00            
##  Max.   :179.07   Max.   :62.993   Max.   :88.00

Tenemos 544 filas, correspondientes a 544 individios. De cada individuo tenemos el registro de su altura (height), peso (weight), edad (age) y si es masculino o no (male). Hagamos un gráfico de dispersión del peso y la altura.

> df %>%
+         ggplot(aes(x=height, y=weight)) + 
+                 geom_point() + 
+                 theme_minimal()

Pregunta: ¿qué pueden decir de esta relación?

Ahora bien, por ahora, vamos a restringirnos a la población adulta (mayores de 18 años). Básicamente, porque la edad y la altura están muy correlacionadas entre si hasta finalizada la adolescencia.

> df %>%
+         filter(age >= 18) %>%
+         ggplot(aes(x=height, y=weight)) + 
+                 geom_point() + 
+                 theme_minimal()

Acá vemos que, si bien la relación entre ambas no es perfectamente lineal, la información sobre el peso de una persona ayuda a predecir la altura de la misma. Queremos, entonces, de alguna forma explicar el peso utilizando la altura. Para eso, vamos a realizar una regresión lineal.

> df_mayores <- df %>%
+         filter(age >= 18)
> 
> 
> df_mayores %>%
+   ggplot(aes(x=height, y=weight)) + 
+                 geom_point() + 
+                 geom_smooth(method='lm', color = 'red', se = FALSE) +
+                 theme_minimal()

La ecuación para esta línea es:

\[\hat{wei} = -52.3162 + 0.6294 \times height\]

El sombrerito arriba de \(wei\) significa que es una estimación y no el valor observado. Podemos usar la ecuación de esta recta para hacer predicciones. Por ejemplo, una persona que mide 1.80 mts (180 cm) va a pesar…

\[\hat{wei} = -52.32 + 0.63 \times height = -52.3162 + 0.6294 \times 180 = 60.98\]

La estimación puede verse como un promedio: la ecuación predice que las personas con una altutra total de 180 cm tendrán un peso de unos 61 kilogramos. A falta de más información sobre una persona de 180 cm, la predicción del peso en función de la altura promedio parece razonable.

Puede haber otras variables que podrían ayudarnos a predecir el peso de una persona de la etnia !Kung San. Tal vez la relación sea diferente entre hombres y no hombres; o quizás sea diferente para diferentes áreas geográficas. El gráfico siguiente muestra la relación entre la altura y el peso de las personas !Kung San separando por la variable male.

> df_mayores %>%
+   ggplot(aes(x=height, y=weight, color=male)) + 
+                 geom_point() + 
+                 theme_minimal()

Los hombres (representados por los puntitos azules) parecen ser más altos y más “pesados”.

El gráfico siguiente muestra la misma relación teniendo en cuenta su edad. Es más difícil saber si la edad cambia la relación entre peso y altura.

> df_mayores %>%
+   ggplot(aes(x=height, y=weight, color=age)) + 
+                 geom_point() + 
+                 scale_color_viridis_c() +
+                 theme_minimal()

Traten de replicar este mismo gráfico con les menores de 18 años.

> ###

Residuos

Los residuos son la variación restante en los datos después de tener en cuenta el ajuste del modelo:

\[ Datos\ observados = Modelo\ ajustado + Residuos \] Vamos a poder calcular para cada observación un residuo. En el gráfico siguiente vemos tres que calculamos para la recta de regresión de nuestro dataset. Si una observación está por encima de la línea de regresión, entonces su residuo (la distancia vertical desde la observación hasta la línea) es positiva. Las observaciones debajo de la línea tienen residuos negativos. Uno de los objetivos al momento de estimar un modelo lineal es que estos residuos sean lo más pequeños posible.

> df_mayores %>%
+   ggplot(aes(x=height, y=weight)) + 
+                 geom_smooth(method = "lm", se = FALSE, color = "lightgrey") +
+                 geom_segment(aes(xend = height, yend = predicted), alpha = .2) +  # alpha to fade lines
+                 geom_point(aes(color = (resid))) +
+                 scale_color_viridis_c() +
+                 theme_minimal()

En este gráfico tenemos mucha información:

  • los datos observados
  • la recta de regresión
  • los residuos y
  • el valor de los residuos.

Cuanto más amarillo es el color de los puntitos, más positiva será la diferencia entre el valor observado y el valor predicho por el modelo. Es decir, el modelo “sub-estima” la observación. A la inversa, cuanto más azul sea el color del puntito, más negativa será la diferencia entre ambos. Es decir, el modelo “sobre-estima”.

Pregunta: ¿Cuál es el punto más sobre-estimado?

Residuos: la diferencia entre lo observado y lo predicho

El residuo de la i-ésima observación \((x_{i}, y_{i})\) es la diferencia entre el valor observado (\(y_{i}\)) y el valor que predice el modelo (\(\hat{y_{i}}\)):

\[\epsilon_{i} = y_{i} - \hat{y_{i}}\]

Veamos el punto de máximo residuo: es el que tiene height=139.7

¿Cómo lo calculamos?

Primero, calculamos el valor predicho: \[\hat{wei} = 52.3162 + 0.6294 \times height = 52.3162 + 0.6294 \times 139.7 = 35.61\]

Luego, computamos la diferencia entre el valor observado y el predicho: \[\epsilon = wei - \hat{wei}= 50.3 - 35.61 = 14.7\]

Pregunta: ¿Cómo interpretamos ese valor?

Los residuos son útiles para evaluar qué tan bien se ajusta (“fitea”) un modelo lineal a un conjunto de datos. Es común representarlos en un gráfico como el que sigue.

> df_mayores %>%
+   ggplot(aes(x=predicted, y=resid)) + # Puede verse cómo calcular la variable resid más abajo...
+     geom_point() +
+     geom_hline(yintercept=0, linetype='dashed') +
+     theme_minimal() +
+     labs(x='Valores predichos de height',
+          y="Residuos")

A menudo los mostramos en un diagrama de dispersión como el anterior. Los residuos en el eje y se grafican con el valor predicho por el modelo en el eje horizontal. Crear un residual es algo así como volcar el gráfico de dispersión para que la línea de regresión sea horizontal, como lo indica la línea discontinua. ¿Cuál sería el valor de máximo residuo?

Uno de los propósitos de las residual plots identificar características o patrones que aún aparecen en los datos después de ajustar un modelo. La siguiente muestra tres diagramas de dispersión con modelos lineales en la primera fila y diagramas residuales en la segunda fila. ¿Pueden identificar algún patrón que quede en los residuos?

En el primer conjunto de datos (primera columna), los residuos no muestran patrones obvios. Los residuos parecen estar dispersos al azar alrededor de la línea discontinua que representa 0.

El segundo conjunto de datos muestra un patrón bien claro en los residuos. Hay cierta curvatura en el gráfico de dispersión, que es más evidente en el gráfico residual. No deberíamos usar una línea recta para modelar estos datos. En su lugar, se debe usar una técnica más avanzada para modelar la relación curva, como cierto tipo de transformaciones de variables que vamos a trabajar más adelante.

El último gráfico muestra muy poca tendencia al alza y los residuos tampoco muestran patrones evidentes. Es razonable tratar de ajustar un modelo lineal a los datos. Sin embargo, no está claro si hay evidencia de que el parámetro de pendiente sea diferente de cero. La estimación puntual del parámetro de pendiente \(\beta_{1}\) no es cero, pero podríamos preguntarnos si esto podría deberse al azar. Esto nos lleva a la cuestión de la inferencia sobre los parámetros de un modelo lineal. Volveremos sobre esto en unas 5 o 6 clases.

Describiendo relaciones lineales con correlaciones

Hemos visto gráficos con relaciones lineales fuertes y otras con relaciones lineales muy débiles. Sería útil si pudiéramos cuantificar la fuerza de estas relaciones lineales con alguna métrica o estadístico.

Correlación: fuerza de una relación lineal.

El coefiiente de correlación lineal (R de Pearson) siempre toma valores entre -1 y 1 y describe la fuerza y dirección de la relación lineal entre dos variables. Denotamos la correlación por \(r\).

El valor \(r\) no tiene unidades y no se verá afectado por un cambio lineal en las unidades (por ejemplo, pasar de pulgadas a centímetros).

Podemos calcular la correlación usando una fórmula, tal como lo hicimos con la media muestral y la desviación estándar.

\[ r_{xy} = \frac{\sum_{i=1}^n (x_{i} - \overline{x}_{i})(y_{i} - \overline{y}_{i})}{\sqrt{\sum_{i=1}^n (x_{i} - \overline{x}_{i})^2} \sqrt{\sum_{i=1}^n (y_{i} - \overline{y}_{i})^2}}\]

La fórmula del R de Pearson es es bastante compleja, como podrán ver. Al igual que con otras estadísticas, generalmente realizamos los cálculos en una computadora o calculadora.

La anterior muestra ocho parcelas y sus correlaciones correspondientes. Solo cuando la relación es perfectamente lineal, la correlación es -1 o 1. Si la relación es fuerte y positiva, la correlación será cercana a +1. Si es fuerte y negativo, estará cerca de -1. Si no existe una relación lineal aparente entre las variables, entonces la correlación será cercana a cero.

La correlación pretende cuantificar la fuerza de una tendencia lineal. Las tendencias no lineales, incluso cuando son fuertes, a veces producen correlaciones que no reflejan la fuerza de la relación. Como se ven en la figura anterior.

Los siguientes diagramas de dispersión muestran las relaciones entre los rendimientos de varios cultivos en los países. En las gráficas, cada punto representa un país diferente. Las variables x e y representan la proporción del rendimiento total en los últimos 50 años que se debe a ese tipo de cultivo.

Ordenen los seis diagramas de dispersión de la relación lineal negativa más fuerte a la positiva más fuerte.

Un aspecto importante del coeficiente \(r\) es que no tiene unidades. Es decir, a diferencia de la pendiente de una línea (volveremos sobre esto la semana que viene) que nos dice cuánto aumenta \(y\) cuando \(X\) aumenta una unidad, en \(r\) no hay unidades asociadas con la correlación de \(X\) e \(y\). La Figura siguiente muestra la relación entre pesos y alturas de 507 individuos físicamente activos.

En el gráfico A, el peso se mide en kilogramos (kg) y la altura en centímetros (cm). En el gráfico B, el peso se ha convertido a libras (lbs) y la altura a pulgadas (in). El coeficiente de correlación (\(r=0.72\) ) es el mismo en ambos gráficos. Tampoco ha cambiado la forma de la relación.El único cambio visual en el gráfico es el etiquetado del eje de los puntos.

Implementando una regresión lineal simple

Hasta aquí la aproximación conceputal a una regresión lineal simple. Vamos a ver ahora el código para implementarla en R. La función que vamos a usar se llama lm() por linear models. Pueden consultar la ayuda de la función para más detalles pero aquí por ahora vamos a continuar con nuestro ejemplo de pesos y alturas.

Vamos a usar height como variable predictora (\(X\)) y weight como variable dependiente (\(y\)). La sintaxis básica de lm() es la siguiente: lm(y~x, data), donde x es la variable independiente e y la dependiente y data la tibble o dataframe con datos.

En nuestro caso, vamos a guardar los resultados en un objeto llamado lm_1.

> lm_1 <- lm(weight~height, data=df_mayores)

Podemos ver que lm_1 es una lista. Por ahora, pidamos un resumen del objeto:

> summary(lm_1)
## 
## Call:
## lm(formula = weight ~ height, data = df_mayores)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -11.8022  -3.0183  -0.2293   2.8117  14.7348 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -52.31618    4.52650  -11.56   <2e-16 ***
## height        0.62942    0.02924   21.52   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.242 on 350 degrees of freedom
## Multiple R-squared:  0.5696, Adjusted R-squared:  0.5684 
## F-statistic: 463.3 on 1 and 350 DF,  p-value: < 2.2e-16

También acá hay bastante información por ahora centremos la mirada en la columna de Estimate. Puede verse que los valores son los mismos que mencionamos al inicio de la clase. Intecept se corresponde con \(\beta_{0} = 52.31\) y height, con \(\beta_{1} = 0.63\).

Podemos usar la función names() para averiguar qué otra información se almacenan en lm_1.

> names(lm_1)
##  [1] "coefficients"  "residuals"     "effects"       "rank"         
##  [5] "fitted.values" "assign"        "qr"            "df.residual"  
##  [9] "xlevels"       "call"          "terms"         "model"

Vemos que están los coeficientes de modelo, los residuos, etc… Contiene mucha información que vamos a ir trabajando a lo largo de las clases. Aunque podemos extraer estas objetos por su nombre lm_1$coficientses más seguro usar el extractor funciones como coef() para acceder a ellas.

> coef(lm_1)
## (Intercept)      height 
##  -52.316178    0.629421

La función predict() puede ser usada para hacer predicciones sobre nuevos datos a partir de la ecuación de la recta que se genera en lm_1. Así, si yo quisiera predecir el valor del peso para una persona de 180 cm. (como hicimos más arriba a mano) podríamos escribir:

> predict(lm_1, data.frame(height = 180)) 
##        1 
## 60.97961

Como pudemos ver, lm_1 tiene un vector llamado fitted.values: son los valores predichos por el modelo. También podemos generarlos mediante predict:

> predict(lm_1, df_mayores) %>%
+   head()
##        1        2        3        4        5        6 
## 43.20791 35.61394 33.61553 46.40537 39.21108 50.80187

De esta forma, podemos generar una columna en el dataframe original y calcular a mano los residuos:

> df_mayores <- df_mayores %>%
+         mutate(weight_predicted = predict(lm_1, df_mayores)) %>%
+         mutate(weight_resid = weight - weight_predicted)
> 
> head(df_mayores)
## # A tibble: 6 × 6
##   height weight   age male  weight_predicted weight_resid
##    <dbl>  <dbl> <dbl> <fct>            <dbl>        <dbl>
## 1   152.   47.8    63 Yes               43.2        4.62 
## 2   140.   36.5    63 No                35.6        0.872
## 3   137.   31.9    65 No                33.6       -1.75 
## 4   157.   53.0    41 Yes               46.4        6.64 
## 5   145.   41.3    51 No                39.2        2.07 
## 6   164.   63.0    35 Yes               50.8       12.2

Por último, podemos calcular los coeficientes de correlación de Pearson \(r\) de la siguiente forma:

> df_mayores %>%
+   select(height, weight) %>%
+   cor(method="pearson")
##           height    weight
## height 1.0000000 0.7547479
## weight 0.7547479 1.0000000

Esto nos devuelve la matriz de correlaciones. Es una matriz cuadrada (tiene la misma cantidad de filas que de columnas). En este caso como tiene solamente dos variables, es una matriz 2x2. En cada celda tenemos la correlación de pearson (\(r\)) para la variable en fila y la variable en la columna. De esta forma, para los mayores de 18 años \(r=0.754\).

Vamos a la práctica…