Elección y ajuste de modelos

La modelización presupone dar un resumen de baja dimensión del dataset. En este contexto lo que buscamos es partir los datos entre muestras que sean representantivas, y otra parte en residuos o ruido. Un patrón fuerte tiende esconder tendencias un poco más útiles, así que es muy habitual hacer uso de técnicas para ir “pelando” las capas de la información del dataset para obtener modelos relevantes, en este proceso de limpieza, modelización y tratamiento de datos de forma constante.

Estudiaremos los fundamentos de cómo funciona un modelo y las diferentes familias que hay:

El objetivo del modelo no es encontrar la verdad, es una aproximación de los datos lo suficientemente buena como para que sirva a nuestro objetivo.

En este modulo haremos uso de la librería modelr, donde trabajaremos con los datos simulados.

library(modelr)
library(tidyverse)
## -- Attaching packages ----------------------------------------------------------------------------------------------- tidyverse 1.2.1 --
## v ggplot2 3.1.0     v purrr   0.2.5
## v tibble  1.4.2     v dplyr   0.7.8
## v tidyr   0.8.2     v stringr 1.3.1
## v readr   1.1.1     v forcats 0.3.0
## -- Conflicts -------------------------------------------------------------------------------------------------- tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()

Introducción a los modelos lineales

Ejemplo de un modelo lineal de inicio a fin usando un dataset dentro del paquete de “modelr”:

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

Podemos ver un fuerte patrón en esta gráfica, parace ser que es un modelo que tendrá una relación lineal entre los patrones. Por lo tanto será:

`y = a~1~ x + a~0~             -> relacion lineal`

Ahora vamos a crear una tibble con muchos modelos posibles dentro de ella, para ver si alguna cuadra con el gráfico superior, si alguna muestra la supuesta relación lineal:

models <- tibble(
  a0 = runif(300,-20,40),
  a1 = runif(300, -5, 5)
)

Visualizamos los modelos creados encima de los valores de “sim1”:

sim1 %>%
  ggplot(aes(x,y))+
    geom_abline(aes(intercept = a0, slope = a1), data = models, alpha = 0.2)+
    geom_point()

Podemos ver que hay algunos modelos que pasan por algunos puntos, pero que en otros no pasan. También hay modelos que no se ajustan en absoluto a ninguno de ellos. Por lo tanto podriamos seguir probando de forma aleatoria si encontramos un modelo, sin embargo hay otros sistemas a seguir mucho mas óptimos. Debemos buscar un modelo que tenga la menor distancia entre la recta y los puntos, como se muestra en la siguiente imagen:

Evaluacion de un modelo lineal

Evaluacion de un modelo lineal

Por lo tanto, deberemos ver la diferencia desde la recta(predicted y) hasta cada uno de los puntos(y). Para solucionar este problema crearemos una función:

model1 <- function(a0, a1, data){
  a0 + data$x * a1
}
model1(3,  1.2, sim1)
##  [1]  4.2  4.2  4.2  5.4  5.4  5.4  6.6  6.6  6.6  7.8  7.8  7.8  9.0  9.0
## [15]  9.0 10.2 10.2 10.2 11.4 11.4 11.4 12.6 12.6 12.6 13.8 13.8 13.8 15.0
## [29] 15.0 15.0
# Buscamos la forma de evaluar si el modelo es bueno o no usando el error cuadrático medio
rmse <- function(mod, data){
  diff <- data$y - model1(mod[1], mod[2], data)
  sqrt(mean(diff^2))
}
# Ejemplo para los siguientes valores
rmse(c(3,1.2),sim1)
## [1] 6.735838

Ejemplo sobre modelo lineal

Calulamos cuál de los modelos es el mejor usando purrr:

sim1_rmse <- function(a0,a1){
  rmse(c(a0,a1),sim1)
}
# Creamos una variable que añada el rmse para cada uno de los modelos 
models <- models %>%
  mutate(rmse = purrr::map2_dbl(a0,a1,sim1_rmse))
# Muestra los 10 mejores modelos con mejor rmse respecto sim1
sim1 %>%
  ggplot(aes(x,y)) +
    geom_point(size = 2, color = "grey30") +
    geom_abline(aes(intercept = a0, slope = a1, color = rmse),
                data = filter(models, rank(rmse) <= 10))

¿Cómo encontrar el mejor modelo?

Tenemos que tomarnos como todos los modelos de la tibble que hemos creado antes, como observaciones. Para eso podemos crear un scatter plot, para comparar los distintos puntos:

models %>%
  ggplot(aes(a0,a1)) +
    geom_point(data = filter(models, rank(rmse) <= 10), size = 4, color = "red") +
    geom_point(aes(color = -rmse))

# con esta penúltima linea de código seleccionamos los 10 mejores modelos y los coloreamos de rojo para que nos lo muestre en el gráfico

Ahora vamos a crear una parrilla que cubra todo el plano y así detectar los mejores modelos, se llamará matriz de puntos equiespaciada.

grid <- expand.grid(
  a0 = seq(-5,20,length = 25),
  a1 = seq(0, 4, length = 25)
) %>%
  mutate(rmse = purrr::map2_dbl(a0, a1, sim1_rmse))
# Pintamos la parrilla encima del gráfico previo:
grid %>%
  ggplot(aes(a0,a1))+
    geom_point(data = filter(grid, rank(rmse) <= 10), size = 4, color = "red")+
    geom_point(aes(color = -rmse)) # lo ponemos en negativo para  que  me queden los mas intensons con el mínimo error cuadrático

Esta matriz de puntos equiespaciada nos muestran los buenos modelos muy concentrados. Ahora podríamos empezara a hacer iteraciones para ir acotando la información, para ir viendo cuáles son estos modelos:

sim1 %>%
  ggplot(aes(x,y))+
    geom_point(size = 2, color = "grey30") +
    geom_abline(
      aes(intercept = a0, slope = a1, color = -rmse),
      data = filter(grid, rank(rmse) <= 10)
    )

Ahora podemos ir mejorando los valores a0 y a1 para ir ajustando el modelo de forma manual:

grid <- expand.grid(
  a0 = seq(3,5,length = 25), # hemos cambiado los valores entre 3 y 5
  a1 = seq(1.9, 2.2, length = 25) # hemos  cambiado os valores entre 1.9 y 2.2
) %>%
  mutate(rmse = purrr::map2_dbl(a0, a1, sim1_rmse))
# Pintamos la parrilla encima del grafico previo:
grid %>%
  ggplot(aes(a0,a1))+
    geom_point(data = filter(grid, rank(rmse) <= 10), size = 4, color = "red")+
    geom_point(aes(color = -rmse)) # lo ponemos en negativo para  que  me queden los mas intensons con el minimo error cuadratico

Ahora representamos la nueva acotación sobre las variables que queremos predecir:

sim1 %>%
  ggplot(aes(x,y))+
    geom_point(size = 2, color = "grey30") +
    geom_abline(
      aes(intercept = a0, slope = a1, color = -rmse),
      data = filter(grid, rank(rmse) <= 10)
    )

Las diferencias de rmse se han visto reducido, por lo tanto hemos podido fabricar unos modelos que se adaptan mejor a a búsqueda de patrones en “sim1”.

La búsqueda de Newton-Raphson

La información para este método está muy completo en wikipedia, en el siguiente gráfico está la información resumida:

Metodo de Newton_Raphson

Metodo de Newton_Raphson

Usaremos la función optim para definir el mejor modelo:

best <- optim(c(0,0), rmse, data = sim1)
# hemos  elegido los valores 0 y 0 a ojo, después hemos elegido la función que queremos optimizar, y luego hemos suministrado los datos.
# Ahora ya hemos definido el mejor modelo. Miremos el parámetro par:
best$par
## [1] 4.222248 2.051204

Estos valore mas óptimos para a a0 y la a1. Hagamos el gráfico:

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

El algoritmo óptimo es capaz de encontrar el mejor modelo para “sim1”. Este sistema es ideal para trabajar con la familia de modelos lineales:

`y = a0 + a1*x1 + a2*x2 + a2*x3 + ... + an*xn`

Sin embargo, R tiene una herramienta que lo hace todo esto de forma automática, que se llama lm, y usa el formato formula, donde se puede definir la variable dependiente y la variable indpendiente.

(lm(y~x, data = sim1) -> sim1_mod)
## 
## Call:
## lm(formula = y ~ x, data = sim1)
## 
## Coefficients:
## (Intercept)            x  
##       4.221        2.052
# si queremos saber cuáles son los coeficientes:
coef(sim1_mod) # vemos que coinciden con los obtenidos per el método de newton-raphson
## (Intercept)           x 
##    4.220822    2.051533
#Usar la función summary para ver toda la información sobre los estadísticos del modelo:
summary(sim1_mod)
## 
## Call:
## lm(formula = y ~ x, data = sim1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.1469 -1.5197  0.1331  1.4670  4.6516 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   4.2208     0.8688   4.858 4.09e-05 ***
## x             2.0515     0.1400  14.651 1.17e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.203 on 28 degrees of freedom
## Multiple R-squared:  0.8846, Adjusted R-squared:  0.8805 
## F-statistic: 214.7 on 1 and 28 DF,  p-value: 1.173e-14

Esta función usa el modelo de la regresión lineal, que es mucho mas rápido y perfecto.

Visualización de predicciones y residuos

Para visualizar las predicciones a partir del modelo, generaremos una matriz de datos equiespaciados, que cubra toda la región que ocupan nuestros datos. Para ello tomaremos le primer argumento del dataset “sim1”, el de la variable indpendiente, la x:

# generamos una parrilla a partir de la columna x de sim1
grid <- sim1 %>%
  data_grid(x)# genera una parrilla de posibles vaores de x 
grid

Ahora haremos las predicciones usando la funcion de la libreria (modelr) llamada add_predictions(), que lo que hara es tomar el data frame y el modelo para hacer las predicciones en forma de una nueva columna que añadiremos al dataset de grid.

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

Ahora “pintamos” las predicciones:

sim1 %>%
  ggplot(aes(x))+
    geom_point(aes(y = y))+ # puntos de los datos originales
    geom_line(aes(y = pred), data = grid, color = "red", size = 1) # pintamos la predicción en rojo y linea

Las preddicciones nos dicen patrón que sigue el modelo, y que el modelo es capaz de explicar. Los residuos es la parte que el modelo es incapaz de explicar, son las diferencias entre el valor observado y la predicción.

Para añadir los residuos, podemos hacer uso de la función add_residuals() que añade los residuos, sobre el dataset original, nos calcula las diferencias sobre el dataset original y la predicción.

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

Para ver las diferencia entre estos valores, la forma mas fácil es hacer un polígono de frequencias o un histograma.

sim1 %>%
  ggplot(aes(resid)) +
  geom_freqpoly(binwidth = 0.5)

Esto es muy útil por que no ayuda a ver si la distribución de los errores sigue una normal o no. En general es lo que queremos, si los errores tubieran una forma exponencial o logarítmica querría decir que los errores tienen una información del modelo que aun no ha sido explicada, por lo tanto el modelo que habríamos creado no sería el correcto. Por lo tanto, lo más normal es que los errores sigan una distribución alrededor de 0.

Otra forma de estudiar los residuos es a través de la comparación del valores de estos respecto al valor de la predicción:

sim1 %>%
  ggplot(aes(x, resid)) + # representamos las variabels x frente a su valor de residuo
  geom_ref_line(h = 0) +# quiero centrarlo en 0
  geom_point()

Parece que estos errores son aleatorios, no siguen ningun patrón. Esto quiere decir que el modelo ha hecho un buen trabajo capturando los patrones del dataset original. Por lo tanto ha sido correcto elegir un modelo lineal para crear nuestro modelo.

Modelos según el tipo de variable

Tenemos que tener en cuenta siempre si la variable que se debe predecir es numérica o categórica. Tomemos un caso de una variable que represente el sexo de una persona:

`y ~ sex <-> y = a0 + a1 * sex`

Este modelo no tiene ningún sentido porque sex es una categoría, y no se puede multiplicar por nada. R convierte esa variable categórica en un valor. En este caso, para hombre el valor será 0, y para mujer, será 1.

`y ~ sex <-> y = a0 + a1 * sex_male (sex_male  = 0, 1)`

Creamos un ejemplo para poder darle sentido:

df <- tribble(
  ~sex, ~value, 
  "male",    1,
  "female",  5,
  "male",    1
)
model_matrix(df, value ~sex)

En este ejemplo R parece que no tiene en cuenta a las mujeres, pero la verdad es que no le hace falta crear una columna de sex female, dado que puede predecir que valores tendrá. Añadir esa columna dará redundancia a nuestros datos.

Vamos a trabajar con otro data frame:

sim2 %>%
  ggplot(aes(x,y)) +
    geom_point()

# este data set tiene cuato categorías

Vamos a crear un modelo en base a este dataset:

mod2 <- lm(y~x, data = sim2)
# Vamos a visuzalizar el modelo:
grid <- sim2 %>%
  data_grid(x) %>%
  add_predictions(mod2)
grid

Loa valores que predice es exactamente la media de cada una de las categorías, es quien minimiza el error cuadrático medio de cada categoría. Vamos a añadir el punto de predicción en las visualizión general:

sim2 %>%
  ggplot(aes(x)) +
    geom_point(aes(y = y)) +
    geom_point(data = grid, aes(y = pred), color = "red", size = 4)

¿Qué ocurre cuando juntamos una categoría y una variable continua en un dataset?

head(sim3)

Vamos a visulizar si información:

sim3 %>%
  ggplot(aes(x1, y)) +
    geom_point(aes(color = x2))

Ahora pasamos a la modelización, podemos hacer un modelo lineal con las suma de las dos variables, estimando que el efecto indpendiente de cada una de ellas por separado sobre el final:

mod1 <- lm(y ~x1 + x2, data = sim3)

Pero también se puede hacer otro modelo que no sea lineal, sino que sea el producto de la variable x1 por la variable x2 del dataset, con el objetivo de querer interrelacionarlas y ver como una afecta a la otra. Sería lo mismo que si hicieramos:

`y = a0 + (a1 * x1 + a2 * x2) + [a12 * x1 * x2]`
mod2 <- lm(y ~ x1 * x2, data = sim3)

Vamos a visualizar los dos modelos:

grid <- sim3 %>%
  data_grid(x1, x2) %>%
  gather_predictions(mod1, mod2) # añadiremos los modelos como una fila del nuevo dataset
head(grid)

Podemos también hacerlo con spread_prediction() para añadir los modelos en un dataframe:

sim3 %>%
  data_grid(x1, x2) %>%
  spread_predictions(mod1, mod2) %>%
  head()

Vamos a hacer las visualizaciones:

sim3 %>%
  ggplot(aes(x1, y, color = x2)) +
  geom_point() +
    geom_line(data = grid, aes(y = pred)) +
    facet_wrap(~model)

Con esta visualización podemos ver las diferencias para cada uno de los modelos, y las predicciones de cada uno según las variables categóricas y numéricas. Podemos ver que para el modelo1 tiene la misma pendiente pero diferenta “intercept”, mientras que el mod2 tiene en cuenta la categoría, cada pendiente es totalmente distinta y ajusta mucha mejor los puntos. El motivo es que el mod2, al tener en cuenta las variables categóricas, parece ser que es capaz de crear una recta totalmente independiente. ¿Cómo sabemos cuál es el mejor? Parece ser que sería el mejor, pero deberíamos comprobarlo con el gráfico de los “residuos”:

sim3 %>%
  gather_residuals(mod1, mod2)  %>% # añadimos los residuos en el dataset
  ggplot(aes(x1, resid, color = x2)) +
    geom_point() +
    facet_grid(model ~x2)

En este gráfico podemos ver que los residuos del modelo dos siguen patrones más caóticos, por lo tanto, al no aparecer patrones comunes, nos indica que sería mejor que el modelo 1. El modelo 1 sí que muestra residuos que siguen patrones, por ejemplo en el factor b no es un buen modelo.

Ahora miremos la interactuación entre dos variables continuas con el dataset que nos queda, “sim4”.

head(sim4)
mod1 <- lm(y ~ x1 + x2, data = sim4)
mod2 <- lm(y ~ x1 * x2, data = sim4)

Creamos la parrilla de datos:

grid <- sim4 %>%
  data_grid( # vamos a crear dos variables con el rango de las variables que queremos analizar, dividido por los trozos que queramos analizar. 
    x1 = seq_range(x1, 5),
    x2 = seq_range(x2, 5)
  ) %>%
  gather_predictions(mod1, mod2)
head(grid)
grid %>%
  ggplot(aes(x1, pred, color = x2)) +
    geom_line() +
    facet_wrap(~model)

Otra forma de visualizar los modelos más clara:

grid %>%
  ggplot(aes(x2, pred, color = x1, group = x1)) +
    geom_line() +
    facet_wrap(~model)

Tranformación de variables

¿Cómo podemos aplicar una transformación, que no es linal, dentro de un modelo? Ejemplo:

  `log(y) ~ sqrt(x1) + x2  <-> log(y) = a0 + a1 * sqrt(x1) + a2 * x2

Si nuestra tranformación involucra sumas, productos, elevados o menos, entre diferentes variables, incluso entre si misma.

log
## function (x, base = exp(1))  .Primitive("log")

Otros modelos no lineales

Los modelos lineales tienen una distribucion de