Práctico 3. Modelos Aditivos Generalizados

Caso 1.

En un estudio del bentos marino se muestrearon 45 puntos a lo largo de la costa de Holanda, registrando en cada punto la riqueza de especies. Se registraron además una serie de variables predictoras: pendiente, exposición, salinidad, temperatura, altura respecto al nivel medio de mareas (NAP), penetrabilidad y tamaño medio de grano del sustrato. En este caso, realizaremos un modelo aditivo generalizado utilizando la altura (NAP) y el factor exposición (exposure) como variables explicativas. Se plantea un GAM ya que los gráficos de dispersión sugieren una relación no lineal entre estas variables y la riqueza de especies, y debido a que al aplicar modelos lineales o modelos lineales generalizados se continua observando patrones entre residuos y predichos. Veremos también la forma de introducir términos aleatorios.

library(mgcv) 
library(gamm4)
library(gratia)
library(tidygam)
library(ggplot2)
library(patchwork)

RIKZ <- read.table("RIKZ.txt", header = TRUE)
#para calcular la riqueza sumamos las filas 2 a la 76
RIKZ$Richness <- rowSums(RIKZ[, 2:76] > 0)
#convertir "exposure" en factor
RIKZ$exposure <- as.factor(RIKZ$exposure)

#exploración
g1 <- ggplot(RIKZ, aes(x= NAP, y = Richness, color = exposure)) + 
  geom_point() + scale_fill_viridis_d()
g2 <- ggplot(RIKZ, aes(x= Richness, y = exposure, fill = exposure)) +
  geom_violin() + scale_fill_viridis_d()
g1 | g2

# GAM
# (no se muestra quasipoisson, tiene scale = 1.14)
# recordar siempre revisar la sobredispersión
fit1 <- gam(Richness ~ exposure + s(NAP), data = RIKZ, 
            family = poisson, method = "REML")
summary(fit1)

# prueba de supuestos del modelo
layout(matrix(1:4,2,2)) 
gam.check(fit1)
layout(1)

# visualización del spline
plot(fit1)
plot(fit1, seWithMean = TRUE, shade = TRUE, 
     shade.col = "palegreen", trans = exp)

# idem en gratia
sm1 <- smooth_estimates(fit1) |> add_confint() |> transform_fun(fun = exp)
ggplot(sm1, aes(NAP, est)) + geom_line() +
  geom_ribbon(aes(ymin = lower_ci, ymax = upper_ci), 
              alpha = 0.2, fill = "forestgreen")

# idem en tidygam
tp1 <- predict_gam(fit1, tran_fun = exp, length_out = 50)
plot(tp1, series = "NAP") 
ggplot(tp1, aes(NAP, Richness, color = exposure)) + geom_line() +
  geom_ribbon(aes(ymin = lower_ci, ymax = upper_ci, 
                  fill = exposure), alpha = 0.2)

# ¿Interacción con un factor?
# smooths separados para cada nivel
fit2 <- gam(Richness ~ exposure +  s(NAP, by = exposure),
            data = RIKZ, family = poisson, method = "REML")
summary(fit2)

layout(matrix(1:4,2,2))
gam.check(fit2)
layout(1)

layout(matrix(1:3, 1, 3)) # se esperan 3 splines
plot(fit2, seWithMean = TRUE, shade = TRUE, 
     shade.col = "palegreen", trans = exp)
layout(1)

sm2 <- smooth_estimates(fit2)
sm2 <- sm2 |> add_confint() |> transform_fun(fun = exp)

g2 <- ggplot(sm2, aes(x = NAP, y = est)) +  
  geom_line(colour = "forestgreen", linewidth = 1.5) +
  geom_ribbon(aes(ymin = lower_ci, ymax = upper_ci), 
              alpha = 0.2, fill = "forestgreen") +
  facet_wrap( . ~ exposure)
g2 #gratia

# smooths con misma rugosidad
# (más adecuado para muchos niveles poco interesantes) 
fit3 <- gam(Richness ~ s(NAP, exposure, bs = "fs"),
            data = RIKZ, family = poisson, method = "REML")
summary(fit3)

layout(matrix(1:4,2,2))
gam.check(fit3)
layout(1)

sm3 <- smooth_estimates(fit3)
sm3 <- sm3 |> add_confint() |> transform_fun(fun = exp)
g3 <- ggplot(sm3, aes(x = NAP, y = est)) +  
  geom_line(colour = "forestgreen", linewidth = 1.5) +
  geom_ribbon(aes(ymin = lower_ci, ymax = upper_ci), 
              alpha = 0.2, fill = "forestgreen") +
  facet_wrap( . ~ exposure)
g3 # gratia

# no hay método en tidygam para bs = fs
# el plot.gam de mgcv es bastante malo

# comparación de modelos
AIC(fit1, fit2, fit3)
anova(fit1, fit2, test = "Chisq")

# puedo poner un efecto random sencillo tipo (1 | Beach)
# entonces puedo usar todas las familias de gam!
fit_r <- gam(Richness ~ exposure + s(NAP) + s(Beach, bs = "re"), data = RIKZ, 
             family = poisson, method = "REML")

AIC(fit1, fit_r)
summary(fit_r)

# puedo hacer lo mismo en gamm4 (menos familias)
# y puedo poner efectos random más complejos
rm1 <- gamm4(Richness ~ exposure + s(NAP), random = ~ (1|Beach), data = RIKZ,
             family = poisson, REML = TRUE)
rm2 <- gamm4(Richness ~ exposure + s(NAP), 
             random = ~ (1|Beach) + (0+NAP|Beach), data = RIKZ,
             family = poisson, REML = TRUE)

AIC(rm1$mer, rm2$mer)
summary(rm1$gam)

sm4 <- smooth_estimates(rm1$gam)
sm4 <- sm4 |> add_confint() |> transform_fun(fun = exp)
g4 <- ggplot(sm4, aes(x = NAP, y = est)) +  
  geom_line(colour = "forestgreen", linewidth = 1.5) +
  geom_ribbon(aes(ymin = lower_ci, ymax = upper_ci), 
              alpha = 0.2, fill = "forestgreen")
g4 # gratia

Caso 2.

En un trabajo de de selección fenotípica sobre rasgos de la orquídea Cyclopogon elatus se intenta conocer si el éxito reproductivo de estas plantas (medido en número de polinios exportados) es influenciado por rasgos como el número de flores y la profundidad de sus nectarios.
Utilizaremos estos datos como ejemplo para obtener medidas de concurvidad, y explorar la posible presencia de interacciones entre variables continuas. En particular, existen dos posibles maneras de analizar la interacción entre variables continuas en un gam.
Los modelos isotrópicos se caracterizan por ser invariantes a la rotación de los ejes, pero no al cambio en la escala. Especialmente buenos para modelos donde ambas variables tienen la misma escala (por ejemplo datos espaciales).
Los modelos invariantes a la escala poseen las propiedades inversas: son afectados por la rotación de los ejes (no-isotrópicos).

cyc <- read.table("cyclop.txt", header = TRUE)
cyc<-na.omit(cyc) # mgcv no admite datos faltantes
colnames(cyc)<-c("nec", "flo", "lar", "fru", "PF", "pol", "PP")

# modelo isotrópico puramente aditivo
# la base por defecto de s() es "tp" thin plate spline
# selec = TRUE para DOBLE PENALIDAD (selección de variables)
m1 <- gam(pol ~ s(nec, k = 10, bs = "tp") + s(flo, k = 10, bs = "tp"), 
          data=cyc, family = poisson, method = "REML", select = TRUE)
summary(m1)

# examen del modelo
concurvity(m1)
layout(matrix(1:4,2,2))
gam.check(m1)
layout(1)

# gráfico para cada spline
layout(matrix(1:2,1,2)) # se esperan 2 splines
plot(m1, seWithMean = TRUE, shade = TRUE, 
     shade.col = "pink", trans = exp)
layout(1)

# gráfico para superficies via mgcv
vis.gam(m1, view = c("nec", "flo"), type = "response",
        plot.type = "contour", color = "cm")

# interacción. Isotrópico, sensible a escala 
m1 <- gam(pol ~ s(nec, k = 10, bs = "tp") + s(flo, k = 10, bs = "tp"), 
          data=cyc, family = poisson, method = "REML", 
          select = TRUE)
m2 <- gam(pol ~ s(nec, flo, k = 100, bs = "tp"), data=cyc, 
          family = poisson, method = "REML", 
          select = TRUE)
AIC(m1, m2)

summary(m2)
vis.gam(m2, view = c("nec", "flo"), type = "response", 
        plot.type = "contour", color = "cm") # ajá

# interacción. No isotrópico, invariantes a la escala
# Modelo aditivo, cambiamos la base a "cr" cubic spline
m3 <- gam(pol ~ s(nec, bs = "cr") + s(flo, bs = "cr"), 
          data=cyc, family = poisson, method = "REML", 
          select = TRUE)
m4 <- gam(pol ~ te(nec, flo), data = cyc, 
          family = poisson, method = "REML",
          select = TRUE)
m5 <- gam(pol ~ s(nec, bs = "cr") + s(flo, bs = "cr") + ti(nec, flo),
          data = cyc, family = poisson, method = "REML", 
          select = TRUE)
AIC(m3, m4, m5)

summary(m5)
vis.gam(m5, view = c("nec", "flo"), type = "response", 
        plot.type = "contour", color = "cm")

# gráfico para superficies via tidygam
m5_p <- predict_gam(m5, tran_fun = exp, length_out = 50)
ggplot(m5_p, aes(nec, flo, z = pol)) + 
  geom_raster(aes(fill = pol)) + 
  geom_contour(colour = "white") +
  scale_fill_viridis_c() +
  theme_minimal()

Ejercicios

El archivo elecciones2023.txt contiene los resultados de las elecciones PASO 2023 para los circuitos electorales de la ciudad de Córdoba, junto a la geolocalización aproximada del centro de estos circuitos. Elija un partido político y construya un mapa de la estructura espacial de sus votos utilizando gam. El set de datos incluye el total de electores (padrón) pero no el número de votos válidos (suma de todos los partidos), ni blancos, nulos o impugnados.

Tip: base para descargar un mapa y usarlo como base en ggplot2

library(OpenStreetMap)
library(ggplot2)
cba_map <- openmap(upperLeft = c(-31.303436, -64.320779), 
                  lowerRight = c(-31.534660, -64.052269), zoom = 12,
                  type = "osm", mergeTiles = TRUE)
sa_map2 <- openproj(sa_map)

cba_plot <- autoplot.OpenStreetMap(sa_map2) +
  geom_point(...)

El archivo wesdr.txt (originalmente del paquete gss) contiene datos sobre la incidencia de retinoplastía diabética (ret, codificado en 0 o 1) y se intenta modelar esta enfermedad en función del tiempo desde la aparición de la diabetes (dur), el índice de masa corporal (bmi) y el porcentaje de hemoglobina glicosilada en sangre (gly). Construya un modelo y utilice matrices de confusión y ROC para revisar su adecuación.
¿Tiene usted datos donde espera asociaciones no lineales? (n > 50). Discuta con sus compañeres de mesa, prestando especial atención a: a) el proceso ue originaron los datos, b) la presencia de estructuras espaciales o temporales c) Presente sus datos con gráficos y análisis.

Práctico 3. Modelos Aditivos Generalizados

Santiago Benitez-Vieyra

Caso 1.

Caso 2.

Ejercicios