Coeficientes de Regresión

\[y = \alpha + \beta x\] Repasemos brevemente cómo interpretar los coeficientes para un modelo de regresión lineal simple. Recuerda que esto implica calcular dos términos de coeficientes que relacionan la variable dependiente con la variable independiente \(x\).

Para nuestros modelos, la variable:

• \(y\): Es la esperanza de vida en relación con el año (variable \(x\))
• \(\alpha\): El coeficiente del intercepto nos dice la esperanza en el año 0. Esto no es significativo para nuestros datos, por lo que no nos centraremos en este término,
• \(\beta\): En cambio, estamos interesados en la estimación del coeficiente del año que, para un modelo lineal simple, corresponde directamente a la pendiente.

Usando la función tidy() en el primer modelo, aprendemos que con cada año que pasa, la esperanza de vida promedio de la población de este país aumenta aproximadamente 0,23 años. Este enfoque puede brindarle información sobre el crecimiento o la falta de crecimiento en la esperanza de vida a lo largo del tiempo para los países que está modelando.

tidy(gap_models$model[[1]])

## # A tibble: 2 × 5
##   term        estimate std.error statistic  p.value
##   <chr>          <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept) -397.     12.4         -32.1 2.48e-37
## 2 year           0.236   0.00622      38.0 3.72e-41

Coeficientes de modelos múltiples

Puede generar estos coeficientes mapeando la función tidy() para cada uno de sus modelos y luego simplificando el nuevo marco de datos usando la función unnest(). Esto da como resultado un tibble que contiene la estimación para cada coeficiente de cada modelo de país.

gap_models %>% 
mutate(coef = map(model,~tidy(.x))) %>% 
  unnest(coef)

## # A tibble: 370 × 8
## # Groups:   country [185]
##    country             data     model  term    estimate std.e…¹ stati…²  p.value
##    <fct>               <list>   <list> <chr>      <dbl>   <dbl>   <dbl>    <dbl>
##  1 Albania             <tibble> <lm>   (Inter… -3.97e+2 1.24e+1   -32.1 2.48e-37
##  2 Albania             <tibble> <lm>   year     2.36e-1 6.22e-3    38.0 3.72e-41
##  3 Algeria             <tibble> <lm>   (Inter… -1.10e+3 4.05e+1   -27.2 1.51e-33
##  4 Algeria             <tibble> <lm>   year     5.86e-1 2.04e-2    28.8 7.84e-35
##  5 Angola              <tibble> <lm>   (Inter… -7.48e+2 1.12e+1   -67.0 2.03e-54
##  6 Angola              <tibble> <lm>   year     4.01e-1 5.62e-3    71.4 6.69e-56
##  7 Antigua and Barbuda <tibble> <lm>   (Inter… -3.79e+2 1.56e+1   -24.2 5.19e-31
##  8 Antigua and Barbuda <tibble> <lm>   year     2.26e-1 7.87e-3    28.8 7.64e-35
##  9 Argentina           <tibble> <lm>   (Inter… -3.56e+2 7.67e+0   -46.4 8.83e-46
## 10 Argentina           <tibble> <lm>   year     2.15e-1 3.86e-3    55.7 4.58e-50
## # … with 360 more rows, and abbreviated variable names ¹​std.error, ²​statistic

Práctica 1 - Ordena los coeficientes de tus modelos

En este ejercicio, aprovechará el flujo de trabajo de la columna de lista junto con la función tidy() de broom para extraer y explorar los coeficientes de los modelos que creo.

Recuerde que el marco de datos gap_models contiene un modelo que predice la esperanza de vida por año para los paises del marco de datos.

# Extraemos las estadísticas de coeficientes de cada modelo del marcos de datos anidados
coef_mod_anidado <- gap_models %>% 
    mutate(coef = map(model, ~tidy(.x)))

# Simplificamos los marcos de datos coef para cada modelo
coef_mod <- coef_mod_anidado %>%
    unnest(coef)

# Trazamos un histograma de las estimaciones de coeficientes para el año
coef_mod %>% 
  filter(term == "year") %>% 
  ggplot(aes(x = estimate)) +
  geom_histogram()

Muy bien, ahora que tenemos la pendiente de cada modelo, exploremos su distribución.

Echemos un vistazo a nuestros modelos

Para hacer esto, usamos map() y glance() para crear un marco de datos de estadísticas de resumen para cada modelo almacenado como la columna coef. Luego puede simplificar estos marcos de datos usando la función unnest(). Esto da como resultado un tibble que contiene las estadísticas del modelo para cada modelo de país.

#

model_perf <- gap_models %>% 
  mutate(coef = map(model,~glance(.x))) %>% 
  unnest(coef)

model_perf

## # A tibble: 185 × 15
## # Groups:   country [185]
##    country    data     model r.squ…¹ adj.r…² sigma stati…³  p.value    df logLik
##    <fct>      <list>   <lis>   <dbl>   <dbl> <dbl>   <dbl>    <dbl> <dbl>  <dbl>
##  1 Albania    <tibble> <lm>    0.963   0.963 0.772  1443.  3.72e-41     1  -65.1
##  2 Algeria    <tibble> <lm>    0.938   0.937 2.53    828.  7.84e-35     1 -133. 
##  3 Angola     <tibble> <lm>    0.989   0.989 0.698  5091.  6.69e-56     1  -59.3
##  4 Antigua a… <tibble> <lm>    0.938   0.937 0.977   828.  7.64e-35     1  -78.6
##  5 Argentina  <tibble> <lm>    0.983   0.982 0.479  3103.  4.58e-50     1  -37.9
##  6 Armenia    <tibble> <lm>    0.288   0.275 1.57     22.2 1.70e- 5     1 -106. 
##  7 Aruba      <tibble> <lm>    0.882   0.880 0.964   412.  3.28e-27     1  -77.8
##  8 Australia  <tibble> <lm>    0.983   0.983 0.540  3240.  1.42e-50     1  -44.7
##  9 Austria    <tibble> <lm>    0.989   0.989 0.430  4949.  1.45e-55     1  -31.7
## 10 Azerbaijan <tibble> <lm>    0.679   0.673 1.54    116.  3.48e-15     1 -105. 
## # … with 175 more rows, 5 more variables: AIC <dbl>, BIC <dbl>, deviance <dbl>,
## #   df.residual <int>, nobs <int>, and abbreviated variable names ¹​r.squared,
## #   ²​adj.r.squared, ³​statistic

Mirando los valores \(R^2\) de los primeros 4 modelos tienen un \(R^2\) alto, lo que sugiere que se han ajustado bien a los datos de ese país.

Modelos con mejor y peor ajuste

Ahora puede explorarf el ajuste de los modelos. Por ejemplo, puede usar la función slice_max() para encontrar los modelos que mejor se ajustan. Asimismo, pueden encontrar los modelos con el peor ajuste utilizando la función slice_min(). Hechemos un vistazo al código y los resultados generados,

mejores_modelos <- model_perf %>%  filter(r.squared >0.995)
mejores_modelos

## # A tibble: 3 × 15
## # Groups:   country [3]
##   country     data     model r.squ…¹ adj.r…² sigma stati…³  p.value    df logLik
##   <fct>       <list>   <lis>   <dbl>   <dbl> <dbl>   <dbl>    <dbl> <dbl>  <dbl>
## 1 Bahamas     <tibble> <lm>    0.995   0.995 0.235  11428. 1.74e-65     1  2.67 
## 2 Israel      <tibble> <lm>    0.996   0.996 0.250  15626. 3.29e-69     1 -0.826
## 3 Switzerland <tibble> <lm>    0.996   0.995 0.244  12349. 2.08e-66     1  0.508
## # … with 5 more variables: AIC <dbl>, BIC <dbl>, deviance <dbl>,
## #   df.residual <int>, nobs <int>, and abbreviated variable names ¹​r.squared,
## #   ²​adj.r.squared, ³​statistic

peores_modelos <- model_perf %>% filter(r.squared < 0.02)
peores_modelos

## # A tibble: 2 × 15
## # Groups:   country [2]
##   country   data     model  r.squared adj.r…¹ sigma stati…² p.value    df logLik
##   <fct>     <list>   <list>     <dbl>   <dbl> <dbl>   <dbl>   <dbl> <dbl>  <dbl>
## 1 Botswana  <tibble> <lm>   0.000821  -0.0173  4.98 0.0452    0.832     1  -171.
## 2 Swaziland <tibble> <lm>   0.0000692 -0.0181  5.83 0.00380   0.951     1  -180.
## # … with 5 more variables: AIC <dbl>, BIC <dbl>, deviance <dbl>,
## #   df.residual <int>, nobs <int>, and abbreviated variable names
## #   ¹​adj.r.squared, ²​statistic

Práctica 2 - Eche un vistazo al ajuste de sus modelos

En esta práctica utilizará glance() para calcular qué tan bien se ajustan los modelos lineales a los datos de cada país.

#  Extraigamos las estadísticas de ajuste de cada modelo en marcos de datos
model_perf_anidado <- gap_models %>% 
    mutate(fit = map(model, ~glance(.x)))

# Simplifiquemos los marcos de datos de ajuste para cada modelo
model_perf <- model_perf_anidado %>% 
    unnest(fit)

# Miremos las primeras seis filas de model_perf
head(model_perf)

## # A tibble: 6 × 15
## # Groups:   country [6]
##   country     data     model r.squ…¹ adj.r…² sigma stati…³  p.value    df logLik
##   <fct>       <list>   <lis>   <dbl>   <dbl> <dbl>   <dbl>    <dbl> <dbl>  <dbl>
## 1 Albania     <tibble> <lm>    0.963   0.963 0.772  1443.  3.72e-41     1  -65.1
## 2 Algeria     <tibble> <lm>    0.938   0.937 2.53    828.  7.84e-35     1 -133. 
## 3 Angola      <tibble> <lm>    0.989   0.989 0.698  5091.  6.69e-56     1  -59.3
## 4 Antigua an… <tibble> <lm>    0.938   0.937 0.977   828.  7.64e-35     1  -78.6
## 5 Argentina   <tibble> <lm>    0.983   0.982 0.479  3103.  4.58e-50     1  -37.9
## 6 Armenia     <tibble> <lm>    0.288   0.275 1.57     22.2 1.70e- 5     1 -106. 
## # … with 5 more variables: AIC <dbl>, BIC <dbl>, deviance <dbl>,
## #   df.residual <int>, nobs <int>, and abbreviated variable names ¹​r.squared,
## #   ²​adj.r.squared, ³​statistic

Ha calculado con éxito las estadísticas de ajuste para los modelos de cada país. A continuación, exploraremos estos resultados.

En este ejercio responderemos las siguientes preguntas:

• En general ¿qué tan bien se ajustan sus modelos a sus datos?
• ¿Cuáles son los modelos que mejor se ajustan?
• ¿Qué modelos no se ajustan bien a los datos?

library(dplyr)
# Trace un histograma de rsquared para los modelos
model_perf %>% 
  ggplot(aes(x = r.squared)) + 
  geom_histogram()

# Extraiga los 4 modelos que mejor se ajustan
best_fit <- model_perf %>% 
  slice_max(r.squared, n = 4)

# Extraiga los 4 modelos con el peor ajuste
worst_fit <- model_perf %>% 
  slice_min(r.squared, n = 4)

Ahora ha preparado dos marcos de datos, uno que contiene los cuatro modelos de mejor ajuste y otro los cuatro modelos de peor ajuste. En la siguiente sección, usaremos la función Augment() para explorar estos ajustes visualmente.

Creación de marco de datos aumentados

Para hacer esto, primero debe crear un marco de datos que contenga tanto los valores predichos como los originales. Esto requiere primero usar map() y Augment() para trabajar en la columna de la lista que contiene los modelos para crear marcos de datos anidados que contengan tanto los valores originales como los predichos. Luego puede usar unnest() en esta nueva columna para simplificar estos marcos de datos y permitir una mayor exploración.

augment_models <- gap_models %>% 
  mutate(augmented = map(model,~augment(.x))) %>% 
  unnest(augmented)

augment_models

## # A tibble: 10,545 × 11
## # Groups:   country [185]
##    country data     model  life_exp…¹  year .fitted .resid   .hat .sigma .cooksd
##    <fct>   <list>   <list>      <dbl> <int>   <dbl>  <dbl>  <dbl>  <dbl>   <dbl>
##  1 Albania <tibble> <lm>         62.9  1960    65.7 -2.80  0.0684  0.672 0.517  
##  2 Albania <tibble> <lm>         63.9  1961    65.9 -1.99  0.0648  0.728 0.245  
##  3 Albania <tibble> <lm>         64.8  1962    66.1 -1.30  0.0614  0.758 0.0989 
##  4 Albania <tibble> <lm>         65.6  1963    66.4 -0.778 0.0581  0.772 0.0332 
##  5 Albania <tibble> <lm>         66.2  1964    66.6 -0.434 0.0549  0.777 0.00970
##  6 Albania <tibble> <lm>         66.6  1965    66.9 -0.260 0.0518  0.779 0.00327
##  7 Albania <tibble> <lm>         66.9  1966    67.1 -0.206 0.0489  0.779 0.00193
##  8 Albania <tibble> <lm>         67.1  1967    67.3 -0.213 0.0461  0.779 0.00192
##  9 Albania <tibble> <lm>         67.3  1968    67.6 -0.239 0.0435  0.779 0.00227
## 10 Albania <tibble> <lm>         67.6  1969    67.8 -0.245 0.0409  0.779 0.00224
## # … with 10,535 more rows, 1 more variable: .std.resid <dbl>, and abbreviated
## #   variable name ¹​life_expectancy

Vamos a visualizar algunos de estos modelos.

Modelo para Bahamas \(R^2: 0.995\)

Note que dado su \(R^2\) es bastante alto, por lo tanto podemos asumir que el modelo lineal se ajusta bien a los datos. Puede comparar el ajuste del modelo con los datos originales trazando ambos en el mismo gráfico. En este ejemplo, usé ggplot2 para trazar los valores originales de la esperanza de vida como un diagrama de dispersión usando la capa geom_point() y agregué el ajuste del modelo lineal como una línea roja usando la capa geom_line().

augment_models %>% filter(country == "Bahamas") %>% 
  ggplot(aes(x = year, y = life_expectancy)) + 
  geom_point() + 
  geom_line(aes(y = .fitted),color = "red") +
  labs(title = "Modelo Regresión Lineal para Bahamas",
           x = "Año",
           y = "Experanza de vida") +
  theme_classic()

Usando esta gráfica, podemos pensar que un modelo de regresión lineal se va ajustando bien a los datos de este país en particular

Modelo para Ukraine \(R^2: 0.17\)

A continuación veamos el modelo de Ukraine, que tiene un valor de r cuadrado más bajo que el modelo de Bahamas. Al graficar estos datos, puede ver que un modelo lineal hace un trabajo deficiente.

augment_models %>% filter(country == "Ukraine") %>% 
  ggplot(aes(x = year, y = life_expectancy)) + 
  geom_point() + 
  geom_line(aes(y = .fitted),color = "red") +
  labs(title = "Modelo Regresión Lineal para Ukraine",
           x = "Año",
           y = "Experanza de vida") +
  theme_classic()

Modelo para Armenia \(R^2 = 0.28\)

Finalmente veamos un ejemplo, donde el r cuadrado continua siendo bajo aunque un poco mejor que el caso previo,. Al graficar estos datos, puede ver que un modelo lineal continua no siendo el mejor.

augment_models %>% filter(country == "Armenia") %>% 
  ggplot(aes(x = year, y = life_expectancy)) + 
  geom_point() + 
  geom_line(aes(y = .fitted),color = "red") +
  labs(title = "Modelo Regresión Lineal para Armenia",
           x = "Año",
           y = "Experanza de vida") +
  theme_classic()

Como pudo ver en estos tres ejemplos, Augment() y ggplot() facilitan la exploración visual del ajuste de un modelo.

Práctica 3 - Aumentar los valores ajustados de cada modelo

En este ejercicio, prepararemos los cuatro modelos de mejor y peor ajuiste para una mayor exploración al aumentar los datos de su modelo con augment().

best_augmented <- mejores_modelos %>% 
  # Creamos el marco de datos aumentado para cada modelo de país
  mutate(augmented = map(model, ~augment(.x))) %>% 
  # Expandimos los marcos de datos aumentados
  unnest(augmented)


worst_augmented <- peores_modelos %>% 
  # Creeamos el marco de datos aumentado para cada modelo de país
  mutate(augmented = map(model, ~augment(.x))) %>% 
  # Expandimos los marcos de datos aumentados
  unnest(augmented)

Ahora, exploremos los mejores y peores modelos de ajuste comparando las lineas ajustadas con los valores reales.

# Mejores 3 modelos 
best_augmented %>% 
  ggplot(aes(x = year)) +
  geom_point(aes(y = life_expectancy)) + 
  geom_line(aes(y = .fitted), color = "red") +
  facet_wrap(~country, scales = "free_y")+
  theme_minimal()

# Peores 2 modelos 
worst_augmented  %>% 
  ggplot(aes(x = year)) +
  geom_point(aes(y = life_expectancy)) + 
  geom_line(aes(y = .fitted), color = "red") +
  facet_wrap(~country, scales = "free_y")+
  theme_minimal()

Parcelas geniales! Puede ver que un modelo lineal hace un gran trabajo para los mejores 3 modelos de ajuste, pero los peores 2 modelos de ajuste no parecen tener una relación lineal. Trabajará para mejorar este ajuste en la próxima serie de ejercicios mediante la incorporación de funciones adicionales.

Modelo de Regresión Lineal Múltiple

\[y = \alpha + \beta_{1}x_{1} + \beta_{2}x_{2} + . . . \]

Este modelo es una extensión natural del modelo de regresión lineal simple. La diferencia clave es que usa más variables explicativas para explicar el resultado, lo que significa que, en lugar de ajustar una linea de mejor ajuste, estamos ajustando un plano multidimensional. En el conjunto de datos de gapminder, podemos usar características adicionals de nuestras observaciones para modelar la esperanza de vida. Entonces, vamos a usarlo.

La elección de que características usar se puede controlar en el campo de fórmula de la función \(lm()\). Recuerda que para un modelo simple usaste la fórmula de la esperanza de vida explicada por año. De manera similar, para un modelo de regresión lineal múltiple, puede definir explícitamente la fórmula incluyendo el nombre de cada característica separada por un signo más o si sabe que desea incluir todas las características, puede capturarlas usando un punto, como veremos posterioremente.

Uso de broom con modelos de regresión lineal múltiple

El comportamiento de las funciones de broom sigue siendo el mismo. tidy() devuelve las estimaciones de los coeficientes de los modelos, esto ahora incluye estimaciones para las cuatro características adicionales. Lo mismo ocurre con augment(), además de los valores ajustados para cada observación, se devuelven los valores de cuatro características nuevas. Y aunque la salida esperada de glance() sigue siendo la misma, tenemos que cambiar nuestro enfoque del valor de r cuadrado al valor de r cuadrado ajustado al evaluar el ajuste de nuestros modelos o comparar modelos de regresión lineal simple y multiple.

\(R^2\) ajustado

Recuerda que R-cuadrado mide la variación explicada por el modelo. Agregar cualquier característica nueva a un modelo, independientemente de su relación con la variable dependiente, siempre aumentará el valor r cuadrado del modelo. Esto se vuelve problemático cuando se compara el ajuste de modelos con diferente número de características explicativas utilizadas. Para compensar esto, en su lugar, utilizará el valor R-cuadrado ajustado, esta es una métrica rcuadrada modificada cuyo cálculo tiene en cuenta la cantidad de características utilizadas en el modelo. La interpretación del valor R-cuadrado ajustado es muy similar al R-cuadrado y lo usará para evaluar el ajuste de sus nuevos modelos y compararlos con los modelos lineales simples creados anteriormente.

Práctica 4: Construir Mejores Modelos

Anteriormente, creó una colección de modelos simples para ajustarse a la expectativa de vida usando la característica de año. Su análisis anterior mostró que algunos de estos modelos no encajaban muy bien.

En este ejercicio, construirá modelos de regresión multiple para cada país utilizando todas las funciones disponibles. Puede que le interese comparar el rendimiento de los dos modelos con el peor ajuste para que su ajuste \(R^2\) se proporcionan a continuación:

# Cree un modelo lineal para cada país usando algunas caracteristicas
gap_fullmodel <- gap_anidado %>% 
  mutate(model = map(data, ~lm(life_expectancy~year+population+fertility+gdp, data = .x)))


fullmodel_perf <- gap_fullmodel %>% 
  # Extraigaimos las estadísticas de ajuste de cada modelo en marcos de datos
  mutate(fit = map(model, ~glance(.x))) %>% 
  # Simplifiquemos los marcos de datos de ajuste para cada modelo
  unnest(fit)

# Vea el rendimiento de los 2 países con el peor ajuste, es decir, 
# los dos modelos simples que viste antes
fullmodel_perf %>% 
  filter(country %in% peores_modelos$country) %>% 
  select(country, adj.r.squared)

## # A tibble: 2 × 2
## # Groups:   country [2]
##   country   adj.r.squared
##   <fct>             <dbl>
## 1 Botswana          0.736
## 2 Swaziland         0.905

Si bien el \(R^2\)-ajustado nos dice qué tan bien se ajusta el modelo a nuestros datos, no da ninguna indicación sobre cómo se desempeñaría con nuevos datos. En el próximo capítulo, aprenderá cómo estimar el rendimiento del modelo utilizando los datos retenidos de la construcción del modelo.