Este texto se basa en los siguientes materiales:

> library(tidyverse)
> library(quantreg)
> library(gapminder)
> library(margins)

Introducción

La visualización de datos es más que generar cifras que muestran los números sin procesar de una tabla de datos. Desde el principio, implica resumir o transformar partes de los datos y luego graficar los resultados.

Los modelos estadísticos son una parte central de ese proceso.

Objetivos

  • Aprender a usar las bibliotecas broom y margins para extraer y trazar de manera ordenada las estimaciones de los modelos que ajustamos
  • Entender brevemente cómo ggplot puede usar varias técnicas de modelado directamente dentro de geoms.

Mirando adentro de los objetos

Aquí discutiremos algunas formas de tomar los modelos que estimamos y extraer información con la que es fácil trabajar en ggplot. Nuestro objetivo es pasar de la forma en que el objeto está almacenado a una tabla tidy de números que podamos plotear.

La mayoría de las clases de modelos estadísticos en R contendrán la información que necesitamos, o tendrán un conjunto especial de funciones o métodos diseñados para extraerla.

Podemos comenzar por aprender un poco más sobre cómo se almacena la salida de los modelos en R. Siempre estamos trabajando con objetos y los objetos tienen una estructura interna que consta de elementos con nombre. A veces, estos son números únicos, a veces vectores y, a veces, listas de cosas como vectores, matrices o fórmulas.

Hemos trabajado mucho con tibbles y dataframes. Estos almacenan tablas de datos con columnas con nombre, que quizás constan de diferentes clases de variables, como números enteros, caracteres, fechas o factores. Los objetos modelo vuelven a ser un poco más complicados.

Vamos a usar información de una organización sin fine de lucro llamada Gapminder. Fundada por Hans Rosling y otros en 2005 ofrece datos y visualizaciones sobre una amplia variedad de temas globales, principalmente relacionados con el desarrollo humano y el progreso social. Entre los temas de información disponible están:

  • Salud: indicadores de esperanza de vida, mortalidad infantil, mortalidad materna, acceso a servicios de salud, etc.
  • Economía: PIB per cápita, niveles de pobreza, desigualdad de ingresos, distribución de riqueza, etc.
  • Educación: tasas de alfabetización, acceso a la educación primaria, secundaria y superior.
  • Población: crecimiento poblacional, distribución por edades, migración, urbanización, etc.
  • Medio ambiente: emisiones de CO₂, uso de recursos naturales, acceso a agua potable, biodiversidad.
  • Desigualdades: desigualdades de género, brechas salariales, desigualdad entre países y dentro de ellos.
  • Tendencias históricas: cómo han cambiado estos indicadores a lo largo del tiempo, en diferentes regiones del mundo.

El paquete gapminder es una herramienta que facilita el acceso a un conjunto de datos basado en los indicadores globales de desarrollo recopilados por la organización Gapminder. Contiene datos de 195 países para los años 1952, 1957, 1962, y cada cinco años hasta 2007.

> gapminder <- gapminder
> str(gapminder)
## tibble [1,704 × 6] (S3: tbl_df/tbl/data.frame)
##  $ country  : Factor w/ 142 levels "Afghanistan",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ continent: Factor w/ 5 levels "Africa","Americas",..: 3 3 3 3 3 3 3 3 3 3 ...
##  $ year     : int [1:1704] 1952 1957 1962 1967 1972 1977 1982 1987 1992 1997 ...
##  $ lifeExp  : num [1:1704] 28.8 30.3 32 34 36.1 ...
##  $ pop      : int [1:1704] 8425333 9240934 10267083 11537966 13079460 14880372 12881816 13867957 16317921 22227415 ...
##  $ gdpPercap: num [1:1704] 779 821 853 836 740 ...

Siempre podemos usar str() para ver la estructura de los diferentes objetos con los que nos encontremos. Aquí hay mucha información sobre el objeto como un todo y cada variable que contiene. Vamos a repasar un poco lo que vimos en las primeras clases de este módulo.

De la misma forma, los modelos estadísticos en R tienen una estructura interna. Pero debido a que los modelos son entidades más complejas que las tablas de datos, su estructura es generalmente más complicada. Hay más elementos y más tipos de información. Toda esta información generalmente se almacena o se puede calcular a partir de partes de un objeto modelo.

Podemos crear un modelo lineal, una regresión MCO, utilizando los datos de gapminder.

Este conjunto de datos tiene una estructura país-año que hace que una especificación OLS como esta sea incorrecta. Pero eso no importa por ahora. Usamos la función lm() para ejecutar el modelo y lo almacenamos en un objeto llamado

> out <- lm(formula = lifeExp ~ gdpPercap + pop + continent, data = gapminder)

El primer argumento es la fórmula del modelo. lifeExp es la variable dependiente y el operador ~ se usa para designar los lados izquierdo y derecho de un modelo (incluso en casos, como vimos con facet_wrap() donde el modelo solo tiene un lado derecho).

Veamos los resultados pidiéndole a R que imprima un resumen del modelo.

El primer argumento es la fórmula del modelo. lifeExp es la variable dependiente y el operador ~ se usa para designar los lados izquierdo y derecho de un modelo (incluso en casos, como vimos con facet_wrap() donde el modelo solo tiene un lado derecho).

Veamos los resultados pidiéndole a R que imprima un resumen del modelo.

> summary(out)
## 
## Call:
## lm(formula = lifeExp ~ gdpPercap + pop + continent, data = gapminder)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -49.161  -4.486   0.297   5.110  25.175 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       4.781e+01  3.395e-01 140.819  < 2e-16 ***
## gdpPercap         4.495e-04  2.346e-05  19.158  < 2e-16 ***
## pop               6.570e-09  1.975e-09   3.326 0.000901 ***
## continentAmericas 1.348e+01  6.000e-01  22.458  < 2e-16 ***
## continentAsia     8.193e+00  5.712e-01  14.342  < 2e-16 ***
## continentEurope   1.747e+01  6.246e-01  27.973  < 2e-16 ***
## continentOceania  1.808e+01  1.782e+00  10.146  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.365 on 1697 degrees of freedom
## Multiple R-squared:  0.5821, Adjusted R-squared:  0.5806 
## F-statistic: 393.9 on 6 and 1697 DF,  p-value: < 2.2e-16

Cuando usamos la función summary() en out, no obtenemos un conjunto de estadísticas de resumen simples.

En este caso, lo que se imprime en la consola es en parte información que se almacena dentro del objeto del modelo y en parte información que la función summary() ha calculado y formateado para mostrarla en la pantalla.

Detrás de escena, summary() obtiene ayuda de otras funciones. Los objetos de diferentes clases tienen métodos predeterminados asociados con ellos, de modo que cuando la función de summary() genérica se aplica a un objeto de modelo lineal, la función sabe que debe pasar el trabajo a una función más especializada que hace un montón de cálculos y formatea adecuadamente a un objeto de modelo lineal.

Usamos la misma función summary() genérica en los dataframes: summary(gapminder). Para nosotros, lo que escribimos es lo mismo, pero en cada caso se aplica un método predeterminado diferente.

summary() proporciona un resumen del modelo, pero realmente no podemos hacer ningún análisis adicional con él directamente. Por ejemplo, ¿qué pasa si queremos plotear alguna información del modelo? La información necesaria para hacer gráficos está dentro del objeto, pero no es obvio cómo usarlo.

> str(out)
## List of 13
##  $ coefficients : Named num [1:7] 4.78e+01 4.50e-04 6.57e-09 1.35e+01 8.19 ...
##   ..- attr(*, "names")= chr [1:7] "(Intercept)" "gdpPercap" "pop" "continentAmericas" ...
##  $ residuals    : Named num [1:1704] -27.6 -26.1 -24.5 -22.4 -20.3 ...
##   ..- attr(*, "names")= chr [1:1704] "1" "2" "3" "4" ...
##  $ effects      : Named num [1:1704] -2455.1 311.1 42.6 101.1 -17.2 ...
##   ..- attr(*, "names")= chr [1:1704] "(Intercept)" "gdpPercap" "pop" "continentAmericas" ...
##  $ rank         : int 7
##  $ fitted.values: Named num [1:1704] 56.4 56.4 56.5 56.5 56.4 ...
##   ..- attr(*, "names")= chr [1:1704] "1" "2" "3" "4" ...
##  $ assign       : int [1:7] 0 1 2 3 3 3 3
##  $ qr           :List of 5
##   ..$ qr   : num [1:1704, 1:7] -41.2795 0.0242 0.0242 0.0242 0.0242 ...
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:1704] "1" "2" "3" "4" ...
##   .. .. ..$ : chr [1:7] "(Intercept)" "gdpPercap" "pop" "continentAmericas" ...
##   .. ..- attr(*, "assign")= int [1:7] 0 1 2 3 3 3 3
##   .. ..- attr(*, "contrasts")=List of 1
##   .. .. ..$ continent: chr "contr.treatment"
##   ..$ qraux: num [1:7] 1.02 1.02 1 1.01 1.04 ...
##   ..$ pivot: int [1:7] 1 2 3 4 5 6 7
##   ..$ tol  : num 1e-07
##   ..$ rank : int 7
##   ..- attr(*, "class")= chr "qr"
##  $ df.residual  : int 1697
##  $ contrasts    :List of 1
##   ..$ continent: chr "contr.treatment"
##  $ xlevels      :List of 1
##   ..$ continent: chr [1:5] "Africa" "Americas" "Asia" "Europe" ...
##  $ call         : language lm(formula = lifeExp ~ gdpPercap + pop + continent, data = gapminder)
##  $ terms        :Classes 'terms', 'formula'  language lifeExp ~ gdpPercap + pop + continent
##   .. ..- attr(*, "variables")= language list(lifeExp, gdpPercap, pop, continent)
##   .. ..- attr(*, "factors")= int [1:4, 1:3] 0 1 0 0 0 0 1 0 0 0 ...
##   .. .. ..- attr(*, "dimnames")=List of 2
##   .. .. .. ..$ : chr [1:4] "lifeExp" "gdpPercap" "pop" "continent"
##   .. .. .. ..$ : chr [1:3] "gdpPercap" "pop" "continent"
##   .. ..- attr(*, "term.labels")= chr [1:3] "gdpPercap" "pop" "continent"
##   .. ..- attr(*, "order")= int [1:3] 1 1 1
##   .. ..- attr(*, "intercept")= int 1
##   .. ..- attr(*, "response")= int 1
##   .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
##   .. ..- attr(*, "predvars")= language list(lifeExp, gdpPercap, pop, continent)
##   .. ..- attr(*, "dataClasses")= Named chr [1:4] "numeric" "numeric" "numeric" "factor"
##   .. .. ..- attr(*, "names")= chr [1:4] "lifeExp" "gdpPercap" "pop" "continent"
##  $ model        :'data.frame':   1704 obs. of  4 variables:
##   ..$ lifeExp  : num [1:1704] 28.8 30.3 32 34 36.1 ...
##   ..$ gdpPercap: num [1:1704] 779 821 853 836 740 ...
##   ..$ pop      : int [1:1704] 8425333 9240934 10267083 11537966 13079460 14880372 12881816 13867957 16317921 22227415 ...
##   ..$ continent: Factor w/ 5 levels "Africa","Americas",..: 3 3 3 3 3 3 3 3 3 3 ...
##   ..- attr(*, "terms")=Classes 'terms', 'formula'  language lifeExp ~ gdpPercap + pop + continent
##   .. .. ..- attr(*, "variables")= language list(lifeExp, gdpPercap, pop, continent)
##   .. .. ..- attr(*, "factors")= int [1:4, 1:3] 0 1 0 0 0 0 1 0 0 0 ...
##   .. .. .. ..- attr(*, "dimnames")=List of 2
##   .. .. .. .. ..$ : chr [1:4] "lifeExp" "gdpPercap" "pop" "continent"
##   .. .. .. .. ..$ : chr [1:3] "gdpPercap" "pop" "continent"
##   .. .. ..- attr(*, "term.labels")= chr [1:3] "gdpPercap" "pop" "continent"
##   .. .. ..- attr(*, "order")= int [1:3] 1 1 1
##   .. .. ..- attr(*, "intercept")= int 1
##   .. .. ..- attr(*, "response")= int 1
##   .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
##   .. .. ..- attr(*, "predvars")= language list(lifeExp, gdpPercap, pop, continent)
##   .. .. ..- attr(*, "dataClasses")= Named chr [1:4] "numeric" "numeric" "numeric" "factor"
##   .. .. .. ..- attr(*, "names")= chr [1:4] "lifeExp" "gdpPercap" "pop" "continent"
##  - attr(*, "class")= chr "lm"

Si echamos un vistazo a la estructura del objeto modelo con str(out), vemos que es un despiole. La mayoría de los objetos coplejos en R (como los outputs de un modelo) out está organizado como una lista con diferentes elementos.

Varios de estos elementos son en sí mismos listas. La figura siguiente ofrece un esquema una vista esquemática del contenido de un objeto de modelo lineal.

En esta lista hay vectores, dataframes, strings e incluso otras listas. Algunos son datos que le pasamos nosotros, otros son datos calculados por el modelo lineal.

Los objetos podrían considerarse organizados como un sistema de archivo: los gabinetes contienen cajones y el cajón puede contener páginas de información, documentos completos o grupos de carpetas con más documentos en su interior.

Podemos acceder a los elementos a través del nombre (y usando la sintaxis de R-base).

> out$coefficients
##       (Intercept)         gdpPercap               pop continentAmericas 
##      4.781408e+01      4.495058e-04      6.569823e-09      1.347594e+01 
##     continentAsia   continentEurope  continentOceania 
##      8.192632e+00      1.747269e+01      1.808330e+01

El resultado summary() se presenta de una manera compacta y eficiente en términos de transmitir información, pero también “no tidyficada”.

Hay una tabla de coeficientes, pero los nombres de las variables están en las filas. Los nombres de las columnas son incómodos y parte de la información (por ejemplo, en la parte inferior de la salida) se ha calculado e impreso, pero no se almacena en el objeto del modelo.

Presentando los modelos de forma correcta

Las cifras basadas en modelos estadísticos enfrentan todos los desafíos ordinarios de la visualización de datos efectiva, y más. Esto se debe a que los resultados del modelo suelen conllevar una considerable carga adicional de interpretación y conocimientos básicos necesarios.

Cuanto más complejo es el modelo, más complicado se vuelve transmitir esta información de manera eficaz y más fácil es llevar al público o a uno mismo al error.

Graficar un modelo está indisolublemente ligado a la estimación del propio modelo: un gráfico, por más bonito que sea, no es un sustituto de la comprensión del modelo en sí.

Algunas buenas prácticas

Veamos algunas ideas generales sobre los buenos gráficos basados en modelos y trabajar con algunos ejemplos de cómo ggplot y algunas bibliotecas adicionales pueden facilitar la obtención de buenos resultados.

1. Presentar los resultados en términos sustantivos:

Los gráficos útiles basados en modelos muestran los resultados de manera sustancialmente significativa y directamente interpretable con respecto a las preguntas que el análisis está tratando de responder. Esto significa mostrar resultados en un contexto en el que otras variables del análisis se mantienen en valores sensibles, como sus medias o medianas.

Con variables continuas, a menudo puede ser útil generar valores pronosticados que cubran algún movimiento sustancialmente significativo a través de la distribución, como del percentil 25 al 75, en lugar de un incremento de una sola unidad en la variable de interés.

Para las variables categóricas no ordenadas, los valores predichos pueden presentarse con respecto a la categoría modal en los datos, o para una categoría particular de interés teórico. Presentar hallazgos sustancialmente interpretables a menudo también significa usar (y a veces convertir a) una escala que los lectores puedan comprender fácilmente.

Si los informes de su modelo dan como resultado log-odds, por ejemplo, convertir las estimaciones en probabilidades predichas facilitará la interpretación. Todos estos consejos son bastante generales. Cada uno de estos puntos se aplica igualmente bien a la presentación de los resultados resumidos en una tabla que en un gráfico. No hay nada distintivamente gráfico en centrarse en el significado sustantivo de sus hallazgos.

2. “Blanquear” los niveles de confianza:

Lo mismo se aplica a la presentación del grado de incertidumbre o confianza que tiene en sus resultados. Las estimaciones del modelo vienen con varias medidas de precisión, confianza, credibilidad o importancia. Presentar e interpretar estas medidas es notoriamente propenso a malas interpretaciones o sobreinterpretaciones, ya que tanto los investigadores como el público exigen más de cosas como los intervalos de confianza y los valores p de lo que estas estadísticas pueden ofrecer.

Como mínimo, después de haber decidido una medida adecuada de ajuste del modelo o la evaluación correcta de la confianza, debe mostrar su rango cuando presente sus resultados. Una familia de geoms de ggplot relacionadas le permite mostrar un rango o intervalo definido por posición en el eje xy luego un rango ymin e ymax en el eje y. Estos geoms incluyen geom_pointrange() y geom_errorbar(), que veremos en acción en breve. Un geom relacionado, geom_ribbon() usa los mismos argumentos para dibujar áreas rellenas y es útil para trazar rangos de valores del eje y a lo largo de algunos ejes x que varían continuamente.

3. (Si es posible) mostrar los datos:

Trazar los resultados de un modelo multivariado generalmente significa una de dos cosas.

Primero, podemos mostrar lo que es en efecto una tabla de coeficientes con medidas de confianza asociadas, quizás organizando los coeficientes en grupos significativos, o por el tamaño de la asociación predicha, o ambos.

En segundo lugar, podemos mostrar los valores predichos de algunas variables (en lugar de solo los coeficientes de un modelo) en algún rango de interés. El último enfoque nos permite mostrar los puntos de datos originales si lo deseamos.

La forma en que ggplot crea gráficos capa por capa nos permite combinar fácilmente las estimaciones del modelo (por ejemplo, una línea de regresión y un rango asociado) y los datos subyacentes. En efecto, estas son versiones construidas manualmente de los gráficos generados automáticamente que hemos estado produciendo con geom_smooth() desde el comienzo de este libro.

Generar predicciones para graficar

Entonces, habiendo ajustado un modelo, es posible que deseemos obtener un gráfico de las estimaciones que produce en el rango de alguna variable en particular, manteniendo constantes otras covariables en algunos valores sensibles.

La función predict() es una forma genérica de utilizar objetos de modelo para producir este tipo de predicción. En R, las funciones “genéricas” toman sus entradas y las pasan a funciones más específicas detrás de escena, aquellas que son adecuadas para trabajar con el tipo particular de objeto modelo que tenemos.

Los detalles de obtener valores predichos de un modelo MCO, por ejemplo, serán algo diferentes de obtener predicciones de una regresión logística. Pero en cada caso podemos usar la misma predict(), teniendo cuidado de verificar la documentación para ver en qué forma se devuelven los resultados para el tipo de modelo con el que estamos trabajando. Muchas de las funciones más utilizadas en R son genéricas de esta manera. summary(), por ejemplo, trabaja en objetos de muchas clases diferentes, desde vectores hasta marcos de datos y modelos estadísticos, produciendo una salida apropiada en cada caso por medio de una función específica de clase en segundo plano.

Para que predict() calcule los nuevos valores por nosotros, necesita algunos datos nuevos para ajustar el modelo. Generaremos un nuevo data.frame cuyas columnas tienen los mismos nombres que las variables en los datos originales del modelo, pero donde las filas tienen nuevos valores.

Una función muy útil llamada expand.grid() nos ayudará a hacer esto. Le daremos una lista de variables, especificando el rango de valores que queremos que tome cada variable. Luego expand.grid() generará el multiplicará el rango completo de valores para todas las combinaciones de los valores que le damos, creando así un data.frame de datos con los nuevos datos que necesitamos.

En el siguiente fragmento de código, usamos min() y max() para obtener los valores mínimo y máximo del PIB per cápita, y luego creamos un vector con cien elementos espaciados uniformemente entre el mínimo y el máximo. Mantenemos la población constante en su mediana y dejamos que continente tome los cinco valores disponibles.

> min_gdp <- min(gapminder$gdpPercap)
> max_gdp <- max(gapminder$gdpPercap)
> med_pop <- median(gapminder$pop)
> 
> pred_df <- expand.grid(gdpPercap = (seq(from = min_gdp, to = max_gdp, length.out = 100)),
+     pop = med_pop, continent = c("Africa", "Americas", "Asia", "Europe", "Oceania"))
> head(pred_df)
##   gdpPercap     pop continent
## 1  241.1659 7023596    Africa
## 2 1385.4282 7023596    Africa
## 3 2529.6905 7023596    Africa
## 4 3673.9528 7023596    Africa
## 5 4818.2150 7023596    Africa
## 6 5962.4773 7023596    Africa

Ahora podemos usar predict(). Si le pasamos nuestros nuevos datos y el modelo, sin ningún otro argumento, va a calcular los valores predichos para cada fila en el data.frame.

Si especificamos interval = 'predict' como argumento, calculará intervalos de predicción del 95% además de la estimación puntual.

> pred_out <- predict(object = out, newdata = pred_df, interval = "predict")
> head(pred_out)
##        fit      lwr      upr
## 1 47.96863 31.54775 64.38951
## 2 48.48298 32.06231 64.90365
## 3 48.99733 32.57670 65.41797
## 4 49.51169 33.09092 65.93245
## 5 50.02604 33.60497 66.44711
## 6 50.54039 34.11885 66.96193

Como sabemos que, por diseño, los casos en pred_df y pred_out corresponden fila por fila, podemos unir los dos data.frame por columna. Este método de unir o fusionar tablas definitivamente no se recomienda cuando se trata de datos más complejos.

> pred_df <- cbind(pred_df, pred_out)
> head(pred_df)
##   gdpPercap     pop continent      fit      lwr      upr
## 1  241.1659 7023596    Africa 47.96863 31.54775 64.38951
## 2 1385.4282 7023596    Africa 48.48298 32.06231 64.90365
## 3 2529.6905 7023596    Africa 48.99733 32.57670 65.41797
## 4 3673.9528 7023596    Africa 49.51169 33.09092 65.93245
## 5 4818.2150 7023596    Africa 50.02604 33.60497 66.44711
## 6 5962.4773 7023596    Africa 50.54039 34.11885 66.96193

El resultado final es un data.frame ordenado, que contiene los valores predichos del modelo para el rango de valores que especificamos.

Ahora podemos graficar los resultados. Debido a que producimos una gama completa de valores predichos, podemos decidir si usarlos o no todos. Aquí subseteamos solo para Europa y África.

> p <- ggplot(data = pred_df %>%
+     filter(continent %in% c("Europe", "Africa")), aes(x = gdpPercap, y = fit, ymin = lwr,
+     ymax = upr, color = continent, fill = continent, group = continent))
> 
> p + geom_point(data = subset(gapminder, continent %in% c("Europe", "Africa")), aes(x = gdpPercap,
+     y = lifeExp, color = continent), alpha = 0.5, inherit.aes = FALSE) + geom_line() +
+     geom_ribbon(alpha = 0.2, color = FALSE) + scale_x_log10(labels = scales::dollar) +
+     theme_minimal()

Usamos una nueva geom aquí para dibujar el área cubierta por los intervalos de predicción: geom_ribbon(). Toma un argumento x como una línea, pero un argumento ymin e ymax como se especifica en el mapeo estético ggplot(). Esto define los límites superior e inferior del intervalo de predicción.

En la práctica, es posible que no utilicemos predict() directamente con tanta frecuencia. En su lugar, podemos escribir código utilizando bibliotecas adicionales que encapsulan el proceso de producir predicciones y diagramas a partir de modelos.

Estos son especialmente útiles cuando su modelo es un poco más complejo y la interpretación de los coeficientes se vuelve más complicada. Esto sucede, por ejemplo, cuando tiene una variable de resultado binaria y necesita convertir los resultados de una regresión logística en probabilidades predichas, o cuando tiene términos de interacción entre sus predicciones. Discutiremos algunas de estas bibliotecas auxiliares en las próximas secciones. Sin embargo, predict() y su capacidad para trabajar de forma segura con diferentes clases de modelos sustenta muchas de esas bibliotecas. Por lo tanto, es útil verlo en acción de primera mano para comprender lo que está haciendo..

Actividad

Generar el mismo gráfico pero para América Latina y Asia

> ###

Graficando efectos marginales

Hasta ahora, nuestro uso de predict() trataba de obtener estimaciones del efecto promedio de algún coeficiente, neto de los otros términos del modelo, como en cualquier regresión lineal. En un modelo lineal sin demasiados vericuetos \(y = \beta_{0} + \beta_{1}edad + \beta_{2}genero_{fem} + \epsilon\), resulta más o menos clara la interpretación de los coeficientes y sus efectos. Cuando tenemos una regresión lineal pero con términos más complicados \(y = \beta_{0} + \beta_{1}edad + \beta_{2}edad^2 + \beta_{3}genero_{fem} + \epsilon\). Si aplicamos derivadas la cosa:

\[ \frac{\partial E(y | edad, genero_{fem})}{\partial edad} = \beta_{1} + 2\beta_{2}edad \]

No hay un valor único del efecto. El efecto de la edad depende del valor de la edad: así hay un efecto a los 20, otro a los 50, etc.Algo similar ocurre con las interacciones. Pero en los modelos no lineales (como una regresión logística) la cosa es más complicada.

Recordemos una forma de expresar un modelo logístico

\[log(\frac{p}{1-p}) = \beta_{0} + \beta_{1}edad + \beta_{2}genero_{fem} + \epsilon\]

Los parámetros estimados están en la escala log-odds. Así, más que el signo, no tienen ninguna interpretación útil. En la ecuación anterior, \(\beta{1}\) es el efecto de la edad en el logartimo del odds ratio del resultado, no en la probabilidad, que es a menudo lo que importa estimar.

Una opción para etimar es presentar los odds-ratios. Pero estos pueden ser malinterpretados. Idealmente, queremos entender lo que dice el modelo en el escala de probabilidad y no en la escala de odds-ratio, mucho menos en la escala de estimación, el log-odds.

En la escala de probabilidad, todos los efectos son no lineales porque, condicional a los valores de las covariables, la probabilidad debe estar acotada entre 0 y 1. Es por eso que vamos a trabajar con otra forma: efectos marginales. O “post-estimación”.

Vamos a utilizar el modelo logístico para introducir efectos marginales. Pero los efectos marginales son aplicables a cualquier otro modelo. Se pueden usar para interpretar modelos lineales con formas funcionales más difíciles, con interacciones o con modelos de Poisson, GLM.

Vamos a estimar un modelo y vamos a utilizar predicciones para ayudarnos a interpretar el modelo. Durante la última década, estimar y graficar efectos parciales o marginales de un modelo se ha convertido en una forma cada vez más común de presentar predicciones precisas e interpretativamente útiles. Particularmente, en modelos de Machine Learning. Esto ha dado lugar a una área de investigación llamada Interpretable Machine Learning.

El interés en las gráficas de efectos marginales se estimuló al darse cuenta de que la interpretación de los términos en los modelos de regresión logística, en particular, era más complicada de lo que parecía, especialmente cuando había términos de interacción en el modelo. El paquete margins de Thomas Leeper puede ayudarnos mucho.

Para ver esto, vamos usar un dataset que contiene algunas características (género, edad, bloque y región) de la votación de la Ley de Interrupción Voluntaria del Embarzo que se produjo el año pasado. Este ejercicio está basado en este post que escribimos en ese momento.

Vamos a cargar los datos y vamos a trasnformar en factor la variable dependiente (votacion_final1).

> dip <- read_csv("https://github.com/gefero/idaes_viz/blob/main/data/proc/diputados_IVE.csv?raw=true")
> 
> dip <- dip %>%
+     mutate(votacion_final1 = as.factor(votacion_final1))

Estimemos, entonces, un modelo de regresión logística. En este modelo vamos a insertar todos las variables y dos bloques de interacción: entre genero y edad y entre bloque y provincia.

> out_bo <- glm(votacion_final1 ~ genero + edad + bloque + region + genero * edad +
+     region * bloque, family = "binomial", data = dip)
> summary(out_bo)
## 
## Call:
## glm(formula = votacion_final1 ~ genero + edad + bloque + region + 
##     genero * edad + region * bloque, family = "binomial", data = dip)
## 
## Coefficients:
##                                              Estimate Std. Error z value
## (Intercept)                                   4.91198    1.37045   3.584
## generoVarón                                  -0.65744    1.58008  -0.416
## edad                                         -0.05417    0.02350  -2.305
## bloqueJxC y aliados                          -2.39342    0.70821  -3.380
## bloqueresto y provinciales                   -0.83887    0.87014  -0.964
## regionCABA                                   -0.12381    1.24010  -0.100
## regionCentro                                 -0.32414    0.94937  -0.341
## regionCuyo                                   -2.68530    0.87814  -3.058
## regionNEA                                    -3.29236    0.88382  -3.725
## regionNOA                                    -1.89282    0.68506  -2.763
## regionPatagonia                               0.02818    0.94542   0.030
## generoVarón:edad                              0.01302    0.02961   0.440
## bloqueJxC y aliados:regionCABA               -0.40224    1.45467  -0.277
## bloqueresto y provinciales:regionCABA        17.40186 2869.10963   0.006
## bloqueJxC y aliados:regionCentro            -18.08596 1715.71204  -0.011
## bloqueresto y provinciales:regionCentro      -0.74205    1.26331  -0.587
## bloqueJxC y aliados:regionCuyo              -15.71397 4593.26973  -0.003
## bloqueresto y provinciales:regionCuyo         0.97106    1.29749   0.748
## bloqueJxC y aliados:regionNEA               -14.99910 3167.20326  -0.005
## bloqueresto y provinciales:regionNEA          0.01104    1.53766   0.007
## bloqueJxC y aliados:regionNOA               -15.91781 2648.57453  -0.006
## bloqueresto y provinciales:regionNOA         -1.35957    1.43911  -0.945
## bloqueJxC y aliados:regionPatagonia          -1.81654    1.53939  -1.180
## bloqueresto y provinciales:regionPatagonia   -1.58574    1.40059  -1.132
##                                            Pr(>|z|)    
## (Intercept)                                0.000338 ***
## generoVarón                                0.677349    
## edad                                       0.021156 *  
## bloqueJxC y aliados                        0.000726 ***
## bloqueresto y provinciales                 0.335018    
## regionCABA                                 0.920473    
## regionCentro                               0.732786    
## regionCuyo                                 0.002228 ** 
## regionNEA                                  0.000195 ***
## regionNOA                                  0.005728 ** 
## regionPatagonia                            0.976218    
## generoVarón:edad                           0.660011    
## bloqueJxC y aliados:regionCABA             0.782153    
## bloqueresto y provinciales:regionCABA      0.995161    
## bloqueJxC y aliados:regionCentro           0.991589    
## bloqueresto y provinciales:regionCentro    0.556946    
## bloqueJxC y aliados:regionCuyo             0.997270    
## bloqueresto y provinciales:regionCuyo      0.454210    
## bloqueJxC y aliados:regionNEA              0.996221    
## bloqueresto y provinciales:regionNEA       0.994269    
## bloqueJxC y aliados:regionNOA              0.995205    
## bloqueresto y provinciales:regionNOA       0.344794    
## bloqueJxC y aliados:regionPatagonia        0.237986    
## bloqueresto y provinciales:regionPatagonia 0.257555    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 354.75  on 255  degrees of freedom
## Residual deviance: 228.80  on 232  degrees of freedom
## AIC: 276.8
## 
## Number of Fisher Scoring iterations: 17

La función summary(), como hemos visto, reporta los coeficientes y alguna otra informacipon. Ahora bien, sobre esta base tenemos varias opciones para graficar información de este modelo.

Podmeos usar margins() para calcular los efectos marginales de cada una de las variables predictoras:

> library(margins)
> bo_m <- margins(out_bo)
> summary(bo_m)
##                      factor     AME     SE        z      p   lower   upper
##         bloqueJxC y aliados -0.5461 0.0486 -11.2388 0.0000 -0.6413 -0.4508
##  bloqueresto y provinciales -0.1878 0.0602  -3.1214 0.0018 -0.3056 -0.0699
##                        edad -0.0069 0.0020  -3.4196 0.0006 -0.0108 -0.0029
##                 generoVarón  0.0037 0.0509   0.0725 0.9422 -0.0960  0.1034
##                  regionCABA  0.0286 0.0832   0.3440 0.7309 -0.1345  0.1917
##                regionCentro -0.1940 0.0746  -2.6021 0.0093 -0.3402 -0.0479
##                  regionCuyo -0.4545 0.0942  -4.8245 0.0000 -0.6391 -0.2698
##                   regionNEA -0.5834 0.0808  -7.2244 0.0000 -0.7417 -0.4251
##                   regionNOA -0.4441 0.0745  -5.9579 0.0000 -0.5902 -0.2980
##             regionPatagonia -0.1745 0.0859  -2.0324 0.0421 -0.3428 -0.0062

La librería margins trae varios métodos de gráficación propios. Por ejemplo, podríamos probar lo siguiente:

> plot(bo_m)

Y vemos un gráfico de los efectos marginales promedio, inspirado en el diseño de stata. De hecho, si consultan el sitio del paquete van a poder ver que el autor intentó replicar la lógica del comando margins de Stata.

Otros métodos de visualización en la librería son

  • cplot(), que visualiza los efectos marginales condicionados a una segunda variable
  • image(), que muestra predicciones o efectos marginales como un mapa de calor relleno o un gráfico de contorno.

También, podemos tomar los resultados de margins y graficarlos nosotros mismos.

Primero, convertimos la salida del summary() en una tibble Luego, usamos mutate(factor=str_replace_all(...) para ordenar las etiquetas. Queremos eliminar los prefijos bloque, genero y region para ajustarlos y hacerlos más prolijos:

> bo_gg <- as_tibble(summary(bo_m))
> prefixes <- c("bloque", "genero", "region")
> replacement <- c("Bloque: ", "Genero: ", "Region: ")
> 
> bo_gg <- bo_gg %>%
+     mutate(factor = str_replace_all(factor, c(bloque = "Bloque: ", genero = "Género: ",
+         region = "Región: ", edad = "Edad")))
> 
> bo_gg %>%
+     select(factor, AME, lower, upper)
## # A tibble: 10 × 4
##    factor                            AME   lower    upper
##    <chr>                           <dbl>   <dbl>    <dbl>
##  1 Bloque: JxC y aliados        -0.546   -0.641  -0.451  
##  2 Bloque: resto y provinciales -0.188   -0.306  -0.0699 
##  3 Edad                         -0.00688 -0.0108 -0.00294
##  4 Género: Varón                 0.00369 -0.0960  0.103  
##  5 Región: CABA                  0.0286  -0.134   0.192  
##  6 Región: Centro               -0.194   -0.340  -0.0479 
##  7 Región: Cuyo                 -0.454   -0.639  -0.270  
##  8 Región: NEA                  -0.583   -0.742  -0.425  
##  9 Región: NOA                  -0.444   -0.590  -0.298  
## 10 Región: Patagonia            -0.174   -0.343  -0.00622
> p <- bo_gg %>%
+     ggplot(aes(x = reorder(factor, AME), y = AME, ymin = lower, ymax = upper))
> 
> p + geom_hline(yintercept = 0, color = "gray80") + geom_pointrange() + coord_flip() +
+     labs(x = NULL, y = "Efecto Marginal Promedio") + theme_minimal()

Si solo estamos interesados en obtener efectos condicionales para una variable en particular, entonces podemos pedir a los métodos de la librería que hagan el trabajo de calcular los efectos por nosotros, pero sin dibujar el gráfico.

Pueden devolver los resultados en un formato que podemos reusar fácilmente en ggplot, y con menos necesidad de limpieza, para la limpieza. Por ejemplo, con cplot():

> pv_cp <- cplot(out_bo, x = "bloque", draw = FALSE)
> pv_cp
##                  xvals     yvals     upper     lower
## 1                  FDT 0.8914640 1.0004901 0.7824378
## 2        JxC y aliados 0.4285822 0.6550378 0.2021266
## 3 resto y provinciales 0.7802142 1.0082583 0.5521702
> 
> p <- ggplot(data = pv_cp, aes(x = reorder(xvals, yvals), y = yvals, ymin = lower,
+     ymax = upper))
> 
> p + geom_hline(yintercept = 0, color = "gray80") + geom_pointrange() + coord_flip() +
+     labs(x = NULL, y = "Efecto condicional") + theme_minimal()

Para cerrar, vamos a llevar esta idea un poco más allá. Vamos a generar datasets “sintéticos” con diferentes valores en cada una de las variables independientes y vamos a hacer predicciones sobre esos datos. La idea es poder tratar de generar más información útil para poder interpretar los coeficientes del modelo.

Primero, generamos tres vectores con los valores que queremos evaluar en tres de las variables: genero, edad y bloque.

> gen <- unique(dip$genero)
> ed <- seq(20, 99, 2)
> bl <- c("JxC y aliados", "FDT", "resto y provinciales")

Segundo, llamamos a la función expand.grid() que va a calcular el producto cartesiano entre todos los vectores que le pasemos como argumentos; esto nos devuelve un data.frame con todas las combinaciones posibles

Tercero, llamamos a mutate y creamos una nueva columna llamada probs que contiene los resulados de hacer las predicciones de probabilidad para cada una de las filas.

Cuarto, graficamos.

> caba <- expand.grid(genero = gen, edad = ed, bloque = bl, region = "CABA") %>%
+     mutate(probs = predict(out_bo, ., type = "response"))
> 
> caba %>%
+     ggplot() + geom_line(aes(x = edad, y = probs, color = genero)) + ylim(0, 1) +
+     facet_wrap(~bloque, scales = "fixed") + theme_minimal() + labs(title = "CABA")

Regresiones, geometrías y “smoothers”

Los histogramas, diagramas de densidad, boxplots y otras geoms_ calculan números únicos o nuevas variables antes de trazarlos. Estos cálculos se realizan mediante funciones stat_, cada una de las cuales trabaja en paralelo a una o varias funciones geom_ predeterminadas, y viceversa.

Las funciones stat_ pueden realizar una buena cantidad de cálculos e incluso estimaciones de modelos sobre la marcha. La función geom_smooth() puede tomar una variedad de argumentos de método para ajustar LOESS, MCO, entre otros.

Tanto las funciones geom_smooth() como geom_quantile() también pueden recibir instrucciones para usar fórmulas diferentes para producir sus ajustes.

> p <- ggplot(data = gapminder, mapping = aes(x = log10(gdpPercap), y = lifeExp))
> 
> p + geom_point(alpha = 0.1) + geom_smooth(color = "tomato", fill = "tomato", method = MASS::rlm) +
+     geom_smooth(color = "steelblue", fill = "steelblue", method = "lm") + theme_minimal()

En el plot superior, accedemos a la función rlm de la biblioteca MASS para ajustar una línea de regresión robusta.

> p + geom_point(alpha = 0.1) + geom_smooth(color = "tomato", method = "lm", size = 1.2,
+     formula = y ~ splines::bs(x, 3), se = FALSE) + theme_minimal()

En el segundo panel, la función bs se invoca directamente desde la biblioteca de splines de la misma manera, para ajustar una curva polinominal a los datos.

Este es el mismo enfoque para acceder directamente a las funciones sin cargar una biblioteca completa que ya hemos usado varias veces cuando usamos funciones de la biblioteca de escalas.

> p + geom_point(alpha = 0.1) + geom_quantile(color = "tomato", size = 1.2, method = "rqss",
+     lambda = 1, quantiles = c(0.2, 0.5, 0.85)) + theme_minimal()

Mientras tanto, la función geom_quantile() es como una versión especializada de geom_smooth() que puede ajustarse a líneas de regresión de cuantiles utilizando una variedad de métodos.

El argumento quantiles toma un vector que especifica los cuantiles en los que se ajustan las líneas.

Mostrando ajustes diferentes de una sola vez

Como acabamos de ver en el primer plot, donde graficamos tanto un MCO como una línea de regresión robusta, podemos observar varios ajustes a la vez en el mismo gráfico colocando capas sobre nuevos suavizadores con geom_smooth().

Siempre que establezcamos el color y la estética del fill en diferentes valores para cada ajuste, podemos distinguirlos visualmente fácilmente. Sin embargo, ggplot no dibujará una leyenda que nos oriente sobre qué ajuste es cuál. Esto se debe a que los suavizadores no están conectados lógicamente entre sí. Existen como capas separadas. ¿Qué pasa si comparamos varios ajustes diferentes y queremos una leyenda que los describa?

Resulta que geom_smooth() puede hacer esto a través de un camino raro:

  1. mapear el color y el fill a una cadena de text que describe el modelo que estamos ajustando, y
  2. luego usar scale_color_manual() y scale_fill_manual() para crear la leyenda.

Primero usamosbrewer.pal () de la biblioteca RColorBrewer para extraer tres colores cualitativamente diferentes de una paleta más grande. Los colores se representan como valores hexadecimales. Como antes, usamos :: para usar la función sin cargar toda la biblioteca:

> model_colors <- RColorBrewer::brewer.pal(3, "Set1")
> model_colors
## [1] "#E41A1C" "#377EB8" "#4DAF4A"

Luego creamos un gráfico con tres smoothers diferentes, mapeando color y fill dentro de la función aes() como el nombre del smoother:

> p0 <- ggplot(data = gapminder, mapping = aes(x = log10(gdpPercap), y = lifeExp))
> 
> p1 <- p0 + geom_point(alpha = 0.2) + geom_smooth(method = "lm", aes(color = "OLS",
+     fill = "OLS")) + geom_smooth(method = "lm", formula = y ~ splines::bs(x, df = 3),
+     aes(color = "Cubic Spline", fill = "Cubic Spline")) + geom_smooth(method = "loess",
+     aes(color = "LOESS", fill = "LOESS"))
> 
> 
> p1 + scale_color_manual(name = "Models", values = model_colors) + scale_fill_manual(name = "Models",
+     values = model_colors) + theme(legend.position = "top") + theme_minimal()

De alguna manera le metimos un poco gato por liebre a ggplot aquí para que funcione.

Hasta ahora, siempre hemos mapeado las aes() a los nombres de las variables, no a strings como “OLS” o “Cubic Splines”. De hecho, en una de las primeras clases del curso vimos que pasaban cosas extrañas cuando hacíamos eso.

Aquí aprovechamos ese comportamiento, creando una nueva variable de valor único para el nombre de cada uno de nuestros modelos. ggplot construirá correctamente la leyenda si llamamos scale_color_manual() y scale_fill_manual(). Recuerden que tenemos que llamar a dos funciones de escala porque tenemos dos mapeos. El resultado es una única leyenda que contiene no solo nuestros tres suavizadores, sino también una leyenda apropiada para guiar al lector.

Estas características de ajuste de modelos hacen que ggplot sea muy útil para el trabajo exploratorio y facilitan la generación y comparación de tendencias basadas en modelos y otros resúmenes como parte del proceso de visualización de datos descriptivos. Las diversas funciones stat_ son una forma flexible de agregar estimaciones resumidas de varios tipos a los gráficos.

Pero también queremos más que esto, incluida la presentación de resultados de modelos que nos ajustamos a nosotros mismos.