Modelo Lineal simple¹

Modelo Básico- Clase 1

library(tidyverse)
library(modelr)
options(na.action = na.warn)

datos sim1

El paquete modlr viene con un set de datos de juguete llamado sim1

ggplot(sim1, aes(x, y)) + 
  geom_point()

Se puede ver un patrón fuerte en los datos. Pareciera que el modelo lineal y = a_0 + a_1 * x podría servir.

Modelos al azar

Para empezar, generemos aleatoriamente varios modelos lineales para ver qué pinta tienen. Para eso, podemos usar geom_abline () que toma una pendiente e intercepto como parámetros.

models <- tibble(
  a1 = runif(250, -20, 40),
  a2 = runif(250, -5, 5)
)
ggplot(sim1, aes(x, y)) + 
  geom_abline(aes(intercept = a1, slope = a2), data = models, alpha = 1/4) +
  geom_point() 

A simple vista podemos apreciar que algunos modelos son mejores que otros. Pero necesitamos una forma de cuantificar cuales son los mejores modelos.

distancias

Una forma de definir mejor es pensar en aquel modelo que minimiza la distancia vertical con cada punto:

Para eso, eligamos un modelo cualquiera:

\[ y= 7 + 1.5*x\]

(para que se vean mejor las distancias, corremos un poquito cada punto sobre el eje x)

dist1 <- sim1 %>% 
  mutate(
    dodge = rep(c(-1, 0, 1) / 20, 10),
    x1 = x + dodge,
    pred = 7 + x1 * 1.5
  )
ggplot(dist1, aes(x1, y)) + 
  geom_abline(intercept = 7, slope = 1.5, colour = "grey40") +
  geom_point(colour = "grey40") +
  geom_linerange(aes(ymin = y, ymax = pred), colour = "#3366FF") 

La distancia de cada punto a la recta es la diferencia entre lo que predice nuestro modelo y el valor real

Para computar la distancia, primero necesitamos una función que represente a nuestro modelo:

Ejercicio 1

Crear una función que reciba un vector con los parametros del modelo, y el set de datos, y genere la predicción:

model1 <- function(a, data) {
  a[1] + data$x * a[2]
}
model1(c(7, 1.5), sim1)

 [1]  8.5  8.5  8.5 10.0 10.0 10.0 11.5 11.5 11.5 13.0 13.0 13.0 14.5 14.5 14.5 16.0 16.0 16.0 17.5 17.5 17.5 19.0 19.0 19.0 20.5 20.5 20.5
[28] 22.0 22.0 22.0

Ahora, necesitamos una forma de calcular los residuos y agruparlos. Como se vió en la clase teórica, esto lo vamos a hacer con el error cuadrático medio

Ejercicio 2

Necesitamos una función, con el mismo input que la anterior, que calcule el promedio de los errores cuadráticos (ECM)

\[ECM = \sqrt\frac{\sum_i^n{(\hat{y_i} - y_i)^2}}{n}\]

measure_distance <- function(mod, data) {
  diff <- data$y - model1(mod, data)
  sqrt(mean(diff ^ 2))
}
measure_distance(c(7, 1.5), sim1)

[1] 2.665212

Evaluando los modelos aleatorios

Ahora podemos calcular el ECM para todos los modelos del dataframe models. Para eso utilizamos el paquete purrr, para ejecutar varias veces la misma función sobre varios elementos.

MAP²

La función map toma un input, una función para aplicar, y alguna otra cosa (por ejemplo parametros que necesite la función)

• map(.x, .f, …)
• map(VECTOR_OR_LIST_INPUT, FUNCTION_TO_APPLY, OPTIONAL_OTHER_STUFF)

Usamos map2 cuando tenemos que pasar dos input, que se aplican sobre una función:

• map2(.x, .y, .f, …)
• map2(INPUT_ONE, INPUT_TWO, FUNCTION_TO_APPLY, OPTIONAL_OTHER_STUFF)

Si tenemos más de dos…

• pmap(.l, .f, …)
• pmap(LIST_OF_INPUT_LISTS, FUNCTION_TO_APPLY, OPTIONAL_OTHER_STUFF)

Nosotros tenemos que pasar pasar los valores de a1 y a2 (dos parámetros –> map2), pero como nuestra función toma sólo uno (el vector a), nos armamos una función de ayuda para wrapear a1 y a2

sim1_dist <- function(a1, a2) {
  measure_distance(c(a1, a2), sim1)
}
models <- models %>% 
  mutate(dist = purrr::map2_dbl(a1, a2, sim1_dist))
models

A continuación, superpongamos los 10 mejores modelos a los datos. Coloreamos los modelos por -dist: esta es una manera fácil de asegurarse de que los mejores modelos (es decir, los que tienen la menor distancia) obtengan los colores más brillantes.

ggplot(sim1, aes(x, y)) + 
  geom_point(size = 2, colour = "grey30") + 
  geom_abline(
    aes(intercept = a1, slope = a2, colour = -dist), 
    data = filter(models, rank(dist) <= 10)
  )

También podemos pensar en estos modelos como observaciones y visualizar con un gráfico de dispersión de a1 vsa2, nuevamente coloreado por -dist. Ya no podemos ver directamente cómo se compara el modelo con los datos, pero podemos ver muchos modelos a la vez. Nuevamente, destaqué los 10 mejores modelos, esta vez dibujando círculos rojos debajo de ellos.

ggplot(models, aes(a1, a2)) +
  geom_point(data = filter(models, rank(dist) <= 10), size = 4, colour = "red") +
  geom_point(aes(colour = -dist))

Grid search

En lugar de probar muchos modelos aleatorios, podríamos ser más sistemáticos y generar una cuadrícula de puntos uniformemente espaciada (esto se denomina grid search). Elegimos los parámetros de la grilla aproximadamente mirando dónde estaban los mejores modelos en el gráfico anterior.

grid <- expand.grid(
  a1 = seq(-5, 20, length = 25),
  a2 = seq(1, 3, length = 25)
  ) %>% 
  mutate(dist = purrr::map2_dbl(a1, a2, sim1_dist))
grid %>% 
  ggplot(aes(a1, a2)) +
  geom_point(data = filter(grid, rank(dist) <= 10), size = 4, colour = "red") +
  geom_point(aes(colour = -dist)) 

Cuando superponemos los 10 mejores modelos en los datos originales, todos se ven bastante bien:

ggplot(sim1, aes(x, y)) + 
  geom_point(size = 2, colour = "grey30") + 
  geom_abline(
    aes(intercept = a1, slope = a2, colour = -dist), 
    data = filter(grid, rank(dist) <= 10)
  )

Modelo Básico- Clase 2

óptimo por métodos numéricos

Podríamos imaginar iterativamente haciendo la cuadrícula más fina y fina hasta que nos centramos en el mejor modelo. Pero hay una forma mejor de abordar ese problema: una herramienta de minimización numérica llamada búsqueda Newton-Raphson.

La intuición de Newton-Raphson es bastante simple: Se elige un punto de partida y se busca la pendiente más inclinada. Luego, desciende por esa pendiente un poco, y se repite una y otra vez, hasta que no se puede seguir bajando.

En R, podemos hacer eso con optim ():

• necesitamos pasarle un vector de puntos iniciales. Elegimos el origen por poner cualquier cosa
• le pasamos nuestra función de distancia, y los parámetros que nuestra función necesita (data)

optim(c(4,2), measure_distance, data = sim1)

$par
[1] 4.221029 2.051528

$value
[1] 2.128181

$counts
function gradient 
      49       NA 

$convergence
[1] 0

$message
NULL

best <- optim(c(0, 0), measure_distance, data = sim1)
best$par

[1] 4.222248 2.051204

ggplot(sim1, aes(x, y)) + 
  geom_point(size = 2, colour = "grey30") + 
  geom_abline(intercept = best$par[1], slope = best$par[2])

Óptimo para el modelo lineal

Este procedimiento es válido para muchas familias de modelos. Pero para el caso del modelo lineal, conocemos otras formas de resolverlo

Si nuestro modelo es

\[ y = a_1 + a_2x + \epsilon \]

La solución del óptima es

\[ \hat{a_1} = \bar{y} - \hat{a_2}\bar{x} \]

\[ \hat{a_2} = \frac{\sum_i^n (y_i -\bar{y})(x_i -\bar{x})}{\sum_i^n (x_i- \bar{x})} \]

R tiene una función específica para el modelo lineal lm(). Cómo esta función sirve tanto para regresiones lineales simples como múltiples, debemos especificar el modelo en las formulas: y ~ x

sim1_mod <- lm(y ~ x, data = sim1)
coef(sim1_mod)

(Intercept)           x 
   4.220822    2.051533 

Predicciones

Para visualizar las predicciones de un modelo, comenzamos por generar una cuadrícula de valores espaciados uniformemente que cubre la región donde se encuentran nuestros datos. La forma más fácil de hacerlo es usar modelr :: data_grid (). Su primer argumento es un dataframe, y para cada argumento posterior encuentra las variables únicas y luego genera todas las combinaciones:

grid <- sim1 %>% 
  data_grid(x) 
grid

(Cuando tengamos modelos de varias variables, esta función va a ser de mucha utilidad)

A continuación, agregamos predicciones. Usaremos modelr :: add_predictions () que toma un dataframe y un modelo. Agrega las predicciones del modelo a una nueva columna en el dataframe:

grid <- grid %>% 
  add_predictions(sim1_mod) 
grid

(También podríamos agregar predicciones al dataset original)

A continuación, trazamos las predicciones. La ventaja de este enfoque respecto de geom_abline () es que funcionará con cualquier modelo en R, desde el más simple hasta el más complejo.

ggplot(sim1, aes(x)) +
  geom_point(aes(y = y)) +
  geom_line(aes(y = pred), data = grid, colour = "red", size = 1)

Residuos del modelo

Agregamos residuos al DF con add_residuals (), que funciona de forma muy similar a add_predictions (). Sin embargo, tengan en cuenta que usamos el conjunto de datos original, no una cuadrícula fabricada. Esto se debe a que para calcular los residuos necesitamos valores y reales.

sim1 <- sim1 %>% 
  add_residuals(sim1_mod)
sim1

Hay varias formas diferentes de entender lo que nos dicen los residuos sobre el modelo. Una forma es simplemente dibujar un polígono de frecuencia que nos ayude a entender la dispersión de los residuos:

ggplot(sim1, aes(resid)) + 
  geom_freqpoly(binwidth = 0.5)

Noten que en promedio los residuos deben ser un valor muy cercano a 0

mean(sim1$resid)

[1] 5.121872e-15

Otra forma de graficar los residuos s

ggplot(sim1, aes(x, resid)) + 
  geom_ref_line(h = 0, size = 2,colour = "firebrick") +
  geom_point() 

Si los residuos no tienen estructura es que hicimos las cosas bien.

---

1. Estas notas estan basadas en Hadley Wickham and Garrett Grolemund. 2017. R for Data Science Import, Tidy, Transform, Visualize, and Model Data (1st ed.). O’Reilly Media, Inc. http://r4ds.had.co.nz.↩

2. basado en https://jennybc.github.io/purrr-tutorial/ls03_map-function-syntax.html, gracias Jenny!↩