library(mgcv) # Fit and interrogate GAMs
library(tidyverse) # Tidy and flexible data manipulation
library(marginaleffects) # Compute conditional and marginal effects
library(ggplot2) # Flexible plotting
library(patchwork) # Combining ggplot objectsGams
Fuente: https://ecogambler.netlify.app/blog/interpreting-gams/
Un GAM (Generalized Additive Model) es un modelo flexible que permite modelar relaciones no lineales entre las variables independientes y la variable dependiente utilizando funciones suavizadas, manteniendo la estructura aditiva. Los GAMs extienden los modelos lineales generalizados (GLM) al permitir que las relaciones entre las variables independientes y la respuesta no sean estrictamente lineales.
Cargamos los datos de CO2 de los conjuntos de datos incorporados que se instalan con R. Estos datos contienen mediciones de un experimento sobre la tolerancia al frío de la especie de hierba Echinochloa crus-galli, sobre el que puede leer más utilizando ?CO2
data(CO2, package = "datasets")
glimpse(CO2)Rows: 84
Columns: 5
$ Plant <ord> Qn1, Qn1, Qn1, Qn1, Qn1, Qn1, Qn1, Qn2, Qn2, Qn2, Qn2, Qn2, …
$ Type <fct> Quebec, Quebec, Quebec, Quebec, Quebec, Quebec, Quebec, Queb…
$ Treatment <fct> nonchilled, nonchilled, nonchilled, nonchilled, nonchilled, …
$ conc <dbl> 95, 175, 250, 350, 500, 675, 1000, 95, 175, 250, 350, 500, 6…
$ uptake <dbl> 16.0, 30.4, 34.8, 37.2, 35.3, 39.2, 39.7, 13.6, 27.3, 37.1, …plant <- CO2 %>%
as_tibble() %>%
rename(plant = Plant, type = Type, treatment = Treatment) %>%
mutate(plant = factor(plant, ordered = FALSE))Representamos los datos:
plant %>%
ggplot(aes(x=conc,y=uptake,color=treatment)) +
geom_line(aes(group=plant)) +
geom_point() +
facet_wrap(~type) 
plant %>%
ggplot(aes(x=treatment,fill=type,y=uptake)) +
geom_boxplot()
plant %>%
ggplot(aes(x=conc,color=type,shape=treatment,y=uptake)) +
geom_point()
plant %>%
ggplot(aes(x=conc,y=uptake,color=treatment)) +
geom_line(aes(group=plant),alpha=0.4) +
geom_point() +
facet_wrap(~type) +
geom_smooth(linewidth=1.1,se=FALSE)
Comenzamos ajustando un modelo aditivo generalizado a estos datos, especificando una interacción paramétrica entre tratamiento y tipo, así como ordenadas aleatorias para cada planta y utilizando la absorción de dióxido de carbono (uptake) como respuesta no negativa de valor real. La parte no lineal procede de una interacción suavizada de factores que permite que un efecto no lineal de la concentración de dióxido de carbono (conc) varíe para diferentes niveles de la variable de tratamiento del frío (treatment). El modelo es de la forma:
\[ \mu = tipo*tratamiento + f_1(concentracion|chilles) + f_2(concentracion|Nonchilled) \]
model_0 <- gam(uptake ~ treatment * type +
s(plant, bs = "re") +
s(conc, by = treatment, k = 7),
data = plant,
method = "REML")summary(model_0)
Family: gaussian
Link function: identity
Formula:
uptake ~ treatment * type + s(plant, bs = "re") + s(conc, by = treatment,
k = 7)
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 35.333 1.298 27.232 < 2e-16 ***
treatmentchilled -3.581 1.835 -1.952 0.0553 .
typeMississippi -9.381 1.835 -5.112 2.96e-06 ***
treatmentchilled:typeMississippi -6.557 2.595 -2.527 0.0139 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(plant) 5.523 8.000 2.23 0.00369 **
s(conc):treatmentnonchilled 4.565 5.180 46.77 < 2e-16 ***
s(conc):treatmentchilled 4.134 4.782 27.35 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
R-sq.(adj) = 0.906 Deviance explained = 92.6%
-REML = 228.64 Scale est. = 10.946 n = 84plot(model_0,pages=1)
Sin embargo, dada la naturaleza de la respuesta, una distribución Gamma() con un enlace logarítmico es un modelo de observación más apropiado. Ello es porque en tal caso estamos asumiendo que:
- La variable dependiente \(Y\) sigue una distribución Gamma. Esto significa que:
\(Y\) es continua y __toma solo valores positivos._. Por ello la distribución Gamma se usa comúnmente para modelar cantidades positivas (como concentraciones, tiempos de espera, etc.).
Su varianza es proporcional al cuadrado de la media. Esto significa que cuando la media \(\mu\) aumenta, la dispersión de \(Y\) también lo hace, cosa que suele ser común en muchas variables (cuanto mayor sea la concentración, o el tiempo de espera, más variabilidad cabe esperar)
La distribución Gamma tiene dos parámetros: una media \(\mu\) y un parámetro de forma \(k\) (que generalmente se ajusta con el modelo).
- El enlace logarítmico implica que la relación entre la media \(\mu\) de la variable dependiente y las variables independientes se da mediante la función logarítmica. Es decir: \[ \log(\mu) = f_1(x_1) + f_2(x_2) + \dots + f_p(x_p) \] donde \(f_1, f_2, \dots, f_p\) son funciones suavizadas (en general, funciones no lineales) que relacionan cada variable independiente \(x_1, x_2, \dots, x_p\) con la media de \(Y\). Esto permite capturar relaciones más complejas que en un modelo lineal simple. El anterior enlace logarítmico se traduce en que: \[ \mu = \exp(f_1(x_1) + f_2(x_2) + \dots + f_p(x_p)) \] Esto garantiza que la media \(\mu\) de \(Y\) sea siempre positiva, lo cual es coherente con la naturaleza de la distribución Gamma.
El ajuste al modelo gam con la distribución gamma y el enlace logarítmico sería:
model_1 <- gam(uptake ~ treatment * type +
s(plant, bs = "re") +
s(conc, by = treatment, k = 7),
data = plant,
method = "REML",
family = Gamma(link = "log"))summary(model_1)
Family: Gamma
Link function: log
Formula:
uptake ~ treatment * type + s(plant, bs = "re") + s(conc, by = treatment,
k = 7)
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.51609 0.06376 55.143 < 2e-16 ***
treatmentchilled -0.11321 0.09018 -1.255 0.213987
typeMississippi -0.31196 0.09018 -3.459 0.000978 ***
treatmentchilled:typeMississippi -0.36041 0.12753 -2.826 0.006311 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(plant) 7.032 8.00 8.022 <2e-16 ***
s(conc):treatmentnonchilled 5.193 5.68 83.888 <2e-16 ***
s(conc):treatmentchilled 5.013 5.55 58.970 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
R-sq.(adj) = 0.957 Deviance explained = 95.5%
-REML = 239.08 Scale est. = 0.010327 n = 84Parece que hay muchos efectos «significativos», pero ¿cómo interpretar estos términos suavizados? El resumen no da ninguna indicación de la magnitud, dirección o forma general de estos efectos suaves. Ni tampoco los coeficientes:
coef(model_1) (Intercept) treatmentchilled
3.516091294 -0.113205794
typeMississippi treatmentchilled:typeMississippi
-0.311956077 -0.360409395
s(plant).1 s(plant).2
-0.041226486 -0.015294663
s(plant).3 s(plant).4
0.056521148 -0.040103411
s(plant).5 s(plant).6
0.026931174 0.013172237
s(plant).7 s(plant).8
-0.055925539 0.049887951
s(plant).9 s(plant).10
0.006037588 -0.215582180
s(plant).11 s(plant).12
0.097118479 0.118463702
s(conc):treatmentnonchilled.1 s(conc):treatmentnonchilled.2
5.232322869 1.133239546
s(conc):treatmentnonchilled.3 s(conc):treatmentnonchilled.4
-1.904830565 0.129187847
s(conc):treatmentnonchilled.5 s(conc):treatmentnonchilled.6
0.236060128 1.202602610
s(conc):treatmentchilled.1 s(conc):treatmentchilled.2
3.690591324 2.117499952
s(conc):treatmentchilled.3 s(conc):treatmentchilled.4
-2.367697749 0.063467378
s(conc):treatmentchilled.5 s(conc):treatmentchilled.6
0.236541597 0.967747793 Los coeficientes etiquetados s(conc):treatmentnonchilled y s(conc):treatmentchilled son los pesos de la función base que controlan la forma de estos splines. Pero sin tener una comprensión completa de como son estas funciones de base y cómo actúan colectivamente para formar el spline, estamos bastante perdidos para la interpretación.
Recorremos 3 pasos que se pueden utilizar para ayudar con la interpretación y presentación de informes de los efectos no lineales de GAM. Estos pasos no son en absoluto exhaustivos.
Paso 1: Visualización de los resultados de los GAM
La visualización es, con diferencia, la mejor herramienta para entender los GAM. Esta es la razón por la que el paquete {mgcv} tiene algunas rutinas de trazado simples pero potentes que podemos aplicar directamente a los objetos gam ajustados.
Gráficos de efectos GAM en la escala de enlace
La visualización más común cuando trabajamos con GAMs son los gráficos de efectos parciales en la escala de enlace. Básicamente, estos gráficos muestran el efecto componente individual de una función suave en la escala de enlace, condicional a que todos los demás términos del modelo se fijen en cero. Dado que la mayoría de los términos suaves en {mgcv} están centrados en cero para mantener la identificabilidad, es conveniente observarlos porque podemos hacer inmediatamente algunas interpretaciones. Por ejemplo, podemos ver el efecto parcial de la concentración para el tratamiento no refrigerado seleccionando el segundo término suave en el modelo:
plot(model_1, select = 2, shade = TRUE)
abline(h = 0, lty = 'dashed')
Este gráfico muestra que, para el tratamiento sin refrigeración, el efecto de la concentración de CO2 sobre el predictor lineal es mayoritariamente negativo (por debajo de cero) para valores pequeños de conc, pero rápidamente se vuelve positivo (por encima de cero) para valores intermedios de conc. A continuación, empieza a estabilizarse en valores mayores de conc. También podemos observar el efecto para el otro nivel de tratamiento:
plot(model_1, select = 3, shade = TRUE)
abline(h = 0, lty = 'dashed')
Estos gráficos son un buen primer paso para interrogar a los GAM, pero tienen varias limitaciones:
Sólo muestran lo que el modelo espera ver si se eliminan todos los demás términos. No son útiles si queremos agregar múltiples suavizados (es decir, ¿cuál es el efecto suavizado medio de conc?) o en situaciones en las que las covariables aparecen en múltiples términos. Por ejemplo, una fórmula como: y ~ s(rain) + s(temp) + s(year) + ti(rain, temp) + ti(rain, year) + ti(temp, year) está perfectamente permitida en gam() (y es muy útil en muchos contextos). Pero en este caso los efectos de nuestros predictores de interés se extienden a través de múltiples funciones suaves y no podemos interpretarlas individualmente.
Estos gráficos se muestran en la escala de enlace, que a menudo puede traducirse de forma no lineal en efectos en la escala de interés (la escala de respuesta. Por ejemplo, tenemos una función de enlace logarítmica en este modelo de observación Gamma(), por lo que los efectos parciales no son los mismos que esperaríamos en la escala real de uptake.
No podemos hacer preguntas más específicas utilizando estos gráficos, tales como ¿Cómo difiere el efecto de la concentración entre los grupos de tratamiento? o ¿Cuánta más absorción esperaríamos si la concentración aumentara de 150 a 300 mL/L?
Para el resto de este artículo, pasaré a hacer gráficos y resúmenes utilizando el paquete {marginaleffects}. Resulta que las funciones disponibles en este paquete hacen que responder a todas las preguntas anteriores sea fácil y sin problemas. Además, podríamos ajustar fácilmente diferentes tipos de modelos (es decir, GLM, modelos de escala de localización, modelos de supervivencia, etc.) a estos mismos datos y utilizar exactamente las mismas funciones para realizar comparaciones de modelos. Para empezar a ilustrarlo, abordaré el primer punto anterior calculando y trazando el efecto suavizado medio de conc en todos los niveles de tratamiento. De nuevo, este gráfico se realiza en la escala de enlace:
plot_predictions(model_1, condition = 'conc',
type = 'link') +
labs(y = "Linear predictor (response scale)",
title = "Average smooth effect of concentration",
subtitle = "Aggregated across treatments and types")
Dividimos estas predicciones según tratamiento y tipo:
plot_predictions(model_1, condition = c('conc', 'treatment', 'type'),
type = 'link') +
labs(y = "Linear predictor (link scale)",
title = "Average smooth effect of concentration",
subtitle = "Per treatment, per type")
Hay muchas más posibilidades para hacer gráficos simples utilizando plot_predictions(). Pero podemos hacer mucho más que trazar predicciones condicionales utilizando este marco. Por ejemplo, a menudo queremos entender cómo cambia la pendiente de una función ajustada a través del rango de una covariable. Podemos calcular las pendientes (es decir, las primeras derivadas) de nuestras funciones suaves utilizando la función plot_slopes(), que tiene argumentos y estructura muy similares a plot_predictions():
plot_slopes(model_1, variables = 'conc',
condition = c('conc', 'treatment'),
type = 'link') +
geom_hline(yintercept = 0, linetype = 'dashed') +
labs(y = "1st derivative of the linear predictor",
title = "Conditional slopes of the concentration effect",
subtitle = "Per treatment, per type")
Este gráfico deja más claro que la pendiente es positiva hasta que alcanzamos valores conc cercanos a 260, más allá de los cuales la función se aplana. Ahora abordaremos el segundo punto anterior, que se refiere probablemente a la segunda pregunta más común de los usuarios cuando trabajan con GAMs en {mgcv}:
¿Cómo puedo trazar efectos suaves en la escala de respuesta?
Gráficos de efectos GAM en la escala de respuesta
De hecho, esto es muy sencillo con {marginaleffects}. Simplemente cambiamos type = ‘response’ y ¡listo! Por supuesto, ahora tenemos opciones para mostrar los valores observados en estos gráficos, que es algo que a menudo queremos hacer:
plot_predictions(model_1, condition = 'conc',
type = 'response', points = 0.5,
rug = TRUE) +
labs(y = "Expected response",
title = "Average smooth effect of concentration",
subtitle = "Aggregated across treatments and types")
Aquí se muestra de nuevo la media suave de conc, pero esta vez en la escala de la variable de respuesta. Los puntos muestran los valores observados, lo que deja claro que un único efecto de conc no bastaría para captar la heterogeneidad de los datos. La división por tratamiento y tipo vuelve a dejar las cosas más claras
plot_predictions(model_1, condition = c('conc', 'treatment', 'type'),
type = 'response', points = 0.6,
rug = TRUE) +
labs(y = "Expected response",
title = "Average smooth effect of concentration",
subtitle = "Per treatment, per type")
También puede trazar pendientes en la escala de la respuesta:
plot_slopes(model_1, variables = 'conc',
condition = c('conc', 'treatment', 'type'),
type = 'response') +
geom_hline(yintercept = 0, linetype = 'dashed') +
labs(y = "1st derivative of the expected response",
title = "Conditional slopes of the concentration effect",
subtitle = "Per treatment, per type")
Estos tipos de gráficos nos dan más capacidad para interrogar y criticar nuestros modelos ajustados. Algunas de las preguntas que puede empezar a plantearse cuando intentamos interpretar o informar sobre funciones a partir de GAM incluyen:
¿Alcanza la función una asíntota en algún punto de su rango?
¿Es la función muy ondulada o parece suave?
Parece que la función tiene múltiples modas, ¿tiene esto sentido?
Existen grandes regiones en las que hay pocos puntos de datos. ¿La incertidumbre de la función crece adecuadamente?
¿Hay puntos que parezcan atípicos, a los que la función pueda estar respondiendo fuertemente?
Diferencias entre suavizados: ¿cómo varía un efecto suavizado entre grupos?
Podemos calcular la diferencia media esperada en los valores de consumo entre tratamientos agregando el tipo y la concentración. Esto nos da el efecto «global» del tratamiento, independientemente del tipo o de los valores de concentración, y puede ser una medida de resumen útil:
Diferencia entre tratamientos según región:
conc_seq <- seq(from = min(plant$conc), max(plant$conc), length.out = 500)
avg_comparisons(model_1, newdata = datagrid(conc = conc_seq,
treatment = unique,
type = unique),
variables = "treatment",
by = "type")
Term type Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
treatment Quebec -4.76 3.04 -1.57 0.117 3.1 -10.7 1.19
treatment Mississippi -10.71 2.20 -4.87 <0.001 19.8 -15.0 -6.40
Type: response
Comparison: mean(chilled) - mean(nonchilled)
Columns: term, contrast, type, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high, predicted_lo, predicted_hi, predicted Diferencia entre regiones según tratamiento
avg_comparisons(model_1, newdata = datagrid(conc = conc_seq,
treatment = unique,
type = unique),
variables = "type",
by = "treatment")
Term treatment Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
type nonchilled -10 2.55 -3.93 <0.001 13.5 -15 -5.02
type chilled -16 2.58 -6.20 <0.001 30.8 -21 -10.93
Type: response
Comparison: mean(Mississippi) - mean(Quebec)
Columns: term, contrast, treatment, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high, predicted_lo, predicted_hi, predicted Diferencia entre tratamientos para valores concretos de concentración:
¿cuál es la diferencia media entre los tratamientos a lo largo de una secuencia de valores conc (es decir, cuál es la diferencia entre los suavizados)?
avg_comparisons(model_1, newdata = datagrid(conc = conc_seq,
treatment = unique,
type = unique),
variables = "treatment",
by = c("conc", "type"))
Term conc type Estimate Std. Error z Pr(>|z|) S 2.5 %
treatment 95.0 Quebec 0.413 1.60 0.258 0.79644 0.3 -2.73
treatment 95.0 Mississippi -3.114 1.02 -3.057 0.00224 8.8 -5.11
treatment 96.8 Quebec 0.373 1.61 0.232 0.81662 0.3 -2.78
treatment 96.8 Mississippi -3.183 1.03 -3.102 0.00192 9.0 -5.19
treatment 98.6 Quebec 0.332 1.62 0.205 0.83731 0.3 -2.84
97.5 %
3.55
-1.12
3.53
-1.17
3.51
--- 990 rows omitted. See ?print.marginaleffects ---
treatment 996.4 Mississippi -11.426 2.75 -4.156 < 0.001 14.9 -16.81
treatment 998.2 Quebec -4.436 3.97 -1.119 0.26334 1.9 -12.21
treatment 998.2 Mississippi -11.427 2.76 -4.147 < 0.001 14.9 -16.83
treatment 1000.0 Quebec -4.435 3.98 -1.115 0.26467 1.9 -12.23
treatment 1000.0 Mississippi -11.428 2.76 -4.138 < 0.001 14.8 -16.84
97.5 %
-6.04
3.34
-6.03
3.36
-6.02
Type: response
Comparison: mean(chilled) - mean(nonchilled)
Columns: term, contrast, conc, type, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high, predicted_lo, predicted_hi, predicted Esto nos da la diferencia esperada entre los dos suavizados para una amplia secuencia de valores de conc. Es extremadamente útil porque ya tiene en cuenta cualquier variación en los interceptos u otros efectos donde el tratamiento o la conc podrían aparecer en el modelo. Por supuesto, mirar a un tibble devuelto con 1.000 filas no es muy útil, así que podemos dibujar estas diferencias en su lugar:
plot_comparisons(model_1, newdata = datagrid(conc = conc_seq,
treatment = unique,
type = unique),
variables = "treatment",
by = c("conc", "type"),
type = 'link') +
geom_hline(yintercept = 0, linetype = 'dashed') +
labs(y = "Estimated difference",
title = "Difference between treatment levels",
subtitle = "Chilled - nonchilled, per type")