Modelos con dependencia y heterogeneidad espacial

Modelos con dependencia y heterogeneidad espacial#

Una característica fundamental de los datos espaciales es la presencia simultánea de dependencia espacial y heterogeneidad espacial. La dependencia espacial se refiere al hecho de que las observaciones cercanas tienden a parecerse más entre sí que las distantes, debido a procesos subyacentes que operan en el espacio geográfico. Por su parte, la heterogeneidad espacial se manifiesta en diferencias sistemáticas entre regiones o ubicaciones, ya sea por condiciones ambientales, sociales o geológicas que varían en el territorio. Para abordar estas dos dimensiones, es necesario emplear modelos espaciales que puedan representarlas explícitamente en su estructura.

Los modelos autoregresivos condicionales (CAR), implementados en la librería INLA, permiten incorporar simultáneamente estos dos aspectos en un marco bayesiano. La dependencia espacial se modela mediante una matriz de vecindad (por ejemplo, tipo Queen o Rook), y la heterogeneidad se introduce agregando un efecto aleatorio no estructurado (iid), como en el modelo BYM (Besag-York-Mollié). Por otro lado, en la familia de modelos simultáneos autoregresivos (SAR) —implementados en R mediante la librería spatialreg— la dependencia espacial puede incorporarse a través de la variable dependiente (modelo SARlag), de los errores (SARerror), o de las covariables (modelo SLX). En estos modelos, la heterogeneidad espacial se representa explícitamente mediante interacciones entre las covariables continuas y un factor espacial categórico, como subcuencas o regiones, permitiendo que los efectos de las variables explicativas cambien según la unidad espacial. Este tipo de interacción permite capturar variaciones estructurales (regímenes espaciales) sin necesidad de introducir efectos aleatorios, y es una forma eficaz de modelar heterogeneidad espacial en un contexto frecuentista.

Marco conceptual: dependencia + heterogeneidad en un mismo modelo#

El capítulo anterior (11_CAR) introdujo cómo los modelos CAR capturan dependencia espacial mediante efectos aleatorios latentes (ICAR, BYM, Leroux). Este capítulo da un paso adelante: combinar esa dependencia espacial con heterogeneidad espacial, es decir, diferencias sistemáticas entre grupos o regiones que no se deben a la vecindad local sino a factores estructurales.

Dos fuentes de variación espacial#

En nuestro caso de estudio (subcuencas de los Andes colombianos distribuidas en tres cuencas principales: Atrato, Cauca y Magdalena), podemos identificar dos niveles de variación:

Nivel	Tipo	Estructura	Modelo
Entre subcuencas vecinas	Dependencia espacial local	Autocorrelación por vecindad	CAR (ICAR/BYM/Leroux)
Entre cuencas (Atrato, Cauca, Magdalena)	Heterogeneidad espacial	Diferencias estructurales por región	Efecto aleatorio i.i.d. por cuenca

Un modelo que capture ambas fuentes tiene la estructura general:

\[\log(\mu_i) = \mathbf{x}_i^T \boldsymbol{\beta} + u_{\text{cuenca}[i]} + \phi_i + \log(a_i)\]

donde:

$\mathbf{x}_i^T \boldsymbol{\beta}$: efectos fijos de las covariables (precipitación, elevación, pendiente)
$u_{\text{cuenca}[i]} \sim \mathcal{N}(0, \sigma_u^2)$: intercepto aleatorio por cuenca (heterogeneidad entre Atrato, Cauca, Magdalena)
$\phi_i$: efecto aleatorio espacial CAR (dependencia entre subcuencas vecinas)
$\log(a_i)$: offset para modelar tasas en lugar de conteos

Progresión de modelos en este capítulo#

Se presentan los modelos en orden de complejidad creciente, todos implementados con INLA:

m3_icar: Efectos fijos + intercepto aleatorio por cuenca (i.i.d.) + efecto espacial ICAR
m4_bym: Efectos fijos + intercepto aleatorio por cuenca (i.i.d.) + efecto espacial BYM (ICAR + i.i.d.)
m5_ler: Efectos fijos + pendiente aleatoria por cuenca + efecto espacial Leroux
m6_spatial_slope: Efectos fijos + intercepto aleatorio + efecto ICAR espacial + pendiente espacial ICAR
m7_basin_slope: Efectos fijos + intercepto aleatorio por cuenca + pendientes aleatorias por cuenca

Finalmente, se presentan los modelos SAR con regímenes espaciales mediante spatialreg, que modelan la heterogeneidad como coeficientes variables por cuenca en lugar de efectos aleatorios.

Objetivos de aprendizaje#

Al finalizar este notebook, el estudiante será capaz de:

Combinar dependencia espacial (efectos CAR) y heterogeneidad espacial (efectos aleatorios por cuenca) en un único modelo bayesiano.
Especificar interceptos y pendientes aleatorias simultáneamente con efectos ICAR/BYM/Leroux en INLA.
Interpretar los efectos aleatorios por cuenca como medida de heterogeneidad regional.
Modelar pendientes espacialmente variables (coeficientes que cambian en el espacio) con INLA.
Contrastar el enfoque bayesiano jerárquico (INLA) con el frecuentista de regímenes espaciales (spatialreg).

Requisitos previos: 11_CAR — modelos CAR con INLA; 12_Jerarquicos — modelos multinivel; R con spatialreg e INLA.

INLA#

Tradicionalmente, los modelos jerárquicos bayesianos han sido una herramienta poderosa para abordar estas características. No obstante, su implementación mediante métodos de muestreo como Markov Chain Monte Carlo (MCMC) puede resultar computacionalmente costosa, especialmente en contextos con grandes volúmenes de datos o estructuras espaciales complejas.

La librería INLA (Integrated Nested Laplace Approximation) para el lenguaje R se presenta como una alternativa eficiente y precisa para ajustar modelos bayesianos latentes gaussianos (LGMs). A través de aproximaciones deterministas, INLA permite estimar la distribución posterior de los parámetros sin recurrir a simulaciones estocásticas, lo que reduce drásticamente los tiempos de cómputo. INLA es especialmente útil para incorporar efectos espaciales estructurados, como los que se modelan mediante procesos autoregresivos condicionales (CAR), y también para representar campos espaciales continuos mediante el enfoque basado en ecuaciones diferenciales estocásticas (SPDE). Esta flexibilidad la convierte en una herramienta ideal para abordar tanto la dependencia espacial estructurada como la variabilidad no estructurada, facilitando la construcción de modelos robustos y adaptados a la complejidad geográfica de los fenómenos estudiados.

A continuación se presentan modelos de regresión de Poisson utilizando la librería INLA, incorporando dos tipos de efectos aleatorios: un intercepto aleatorio para cuenca y un efecto espacial CAR para las subcuencas.

Estructura interna de los modelos INLA con efectos mixtos#

La función f() en INLA es el mecanismo para especificar efectos aleatorios latentes. Su sintaxis general es:

f(indice, [covariable,] model = "tipo", graph = matriz_vecindad)

Los dos usos principales en este capítulo son:

Especificación	Código INLA	Significado
Intercepto aleatorio por cuenca	`f(cuenca, model = "iid")`	Cada cuenca tiene su propio nivel base de riesgo
Pendiente aleatoria por cuenca para `x`	`f(cuenca, x, model = "iid")`	El efecto de `x` varía por cuenca
Efecto espacial ICAR	`f(id, model = "besag", graph = aoi.mat)`	Suavizado espacial entre subcuencas vecinas
Efecto espacial BYM	`f(id, model = "bym", graph = aoi.mat)`	ICAR + i.i.d. por subcuenca
Efecto espacial Leroux	`f(id, model = "generic1", Cmatrix = Cmatrix)`	CAR propio con $\rho$ estimado

Nota sobre índices duplicados: INLA no permite usar el mismo índice dos veces en la misma fórmula. Por eso, en los modelos con intercepto aleatorio por cuenca (f(cuenca, model="iid")) y pendiente aleatoria por cuenca (f(cuenca, x, model="iid")), se debe crear un índice duplicado o usar índices distintos aunque refieran a la misma variable.

library(sf) # Simple features for spatial data
library(spdep) # Spatial dependence tools
library(spatialreg) # Spatial regression models
library(sf)
library(dplyr)
library(sp)
library(spdep)
library(INLA) 

aoi = st_read("https://github.com/edieraristizabal/ModeloMultinivel/blob/main/DATA/df_catchments_spatial.gpkg")
aoi <- aoi %>% mutate_at(c('elev_mean','slope_mean','RainfallDaysmean'), ~(scale(.) %>% as.vector))
aoi <- aoi %>% mutate(cuenca_num = case_when(cuenca == 'Atrato' ~ 1, cuenca == 'Cauca' ~ 2, cuenca == 'Magdalena' ~ 3))
aoi_sp <- as_Spatial(aoi)

#Spatial matrix
W.nb <- poly2nb(aoi)
W.mat <- nb2mat(W.nb, style="B")
W.list <- nb2listw(W.nb, style="B")


# intercepto aleatorio + ICAR model
m3_icar <- inla(formula = lands_rec ~ 1 + RainfallDaysmean + elev_mean + slope_mean + 
                  f(cuenca, model = "iid") +
                  f(id, model = "besag", graph = aoi.mat),
                  offset=log(area), data = as.data.frame(aoi), family ="poisson",
                  control.predictor = list(compute = TRUE),
                  control.compute = list(dic = TRUE, waic = TRUE))

summary(m3_icar)
m3_icar$summary.random$cuenca
aoi$ICAR <- m3_icar$summary.fitted[1:nrow(aoi), "mean"]
res_pearson_m3 <- (aoi$lands_rec - m3_icar$summary.fitted.values$mean) / sqrt(m3_icar$summary.fitted.values$mean)
aoi$res_pearson_m3 <- res_pearson_m3
moran_m3 <- moran.mc(x = res_pearson_m3, listw = aoi.listw,nsim = 999, alternative = "greater")

Interpretación del modelo m3_icar (intercepto aleatorio + ICAR)#

El modelo m3_icar combina dos fuentes de variación espacial:

1. Efectos fijos (summary(m3_icar)$fixed)

Los coeficientes de RainfallDaysmean, elev_mean y slope_mean representan el efecto promedio de cada variable sobre la tasa de deslizamientos, después de controlar tanto las diferencias entre cuencas como la autocorrelación espacial local. Al estar en escala log (modelo Poisson), se interpretan como: por cada unidad de aumento estandarizado en la variable $k$, la tasa esperada de deslizamientos se multiplica por $e^{\beta_k}$.

2. Efectos aleatorios de cuenca (m3_icar$summary.random$cuenca)

Este objeto devuelve una tabla con los interceptos estimados para cada una de las tres cuencas (Atrato=1, Cauca=2, Magdalena=3):

mean: intercepto aleatorio estimado $\hat{u}_{\text{cuenca}}$ — positivo indica que la cuenca tiene una tasa más alta que el promedio global; negativo, más baja.
0.025quant, 0.975quant: intervalo de credibilidad del 95% del intercepto.

Si los tres intervalos de credibilidad no se solapan, hay evidencia de heterogeneidad significativa entre cuencas. Si los intervalos se solapan fuertemente, la diferencia entre cuencas es pequeña y el intercepto aleatorio contribuye poco al ajuste.

3. Test de Moran sobre residuos (moran_m3)

Comparar con moran_m1 (modelo base): si el p-valor de Moran pasa de significativo a no significativo, el efecto ICAR ha absorbido la dependencia espacial residual del modelo base.

# BYM model
m4_bym<-inla(formula = lands_rec ~ 1 + RainfallDaysmean + elev_mean + slope_mean +
           f(cuenca, model = "iid") +
           f(id, model = "bym",graph = aoi.mat), 
           offset=log(area), data = as.data.frame(aoi), family = "poisson",
           control.compute = list(dic = TRUE, waic=T), 
           control.predictor = list(compute = TRUE))

summary(m4_bym)
m4_bym$summary.random$cuenca
aoi$BYM <- m4_bym$summary.fitted[1:nrow(aoi), "mean"]
aoi_sp$BYM <- m4_bym$summary.fitted.values[, "mean"]
res_pearson_m4 <- (aoi$lands_rec - m4_bym$summary.fitted.values$mean) / sqrt(m4_bym$summary.fitted.values$mean)
aoi$res_pearson_m4 <- res_pearson_m4
moran_m4 <- moran.mc(x = res_pearson_m4, listw = aoi.listw,nsim = 999, alternative = "greater")

Interpretación del modelo m4_bym (intercepto aleatorio + BYM)#

El modelo BYM añade sobre el ICAR un componente no estructurado a nivel de subcuenca (f(id, model="bym")), que puede capturar sobredispersión o variabilidad idiosincrásica que no sigue el patrón de vecindad.

¿BYM vs ICAR: cuándo aporta el BYM?

El BYM introduce dos hiperparámetros de precisión en lugar de uno:

Precision for ID (spatial) → controla la varianza del componente ICAR
Precision for ID (iid) → controla la varianza del componente no estructurado

Si la precisión del componente i.i.d. es muy alta (varianza cercana a cero), significa que el componente no estructurado es despreciable y el BYM colapsa al ICAR. En ese caso, el DIC/WAIC del BYM no debería ser sustancialmente mejor que el del ICAR.

Si la precisión del componente ICAR es muy alta, significa que el suavizado espacial es débil y la variación es mayormente no estructurada. En ese caso, un modelo i.i.d. puro podría ser suficiente.

Lectura práctica del summary(m4_bym):

Precision for ID (spatial): [valor]  → 1/σ²_espacial
Precision for ID (iid):     [valor]  → 1/σ²_no_estructurado

Comparar los inversos de estas precisiones (varianzas) da una idea de la proporción de la variabilidad espacial explicada por estructura vs. ruido.

Comparación DIC/WAIC entre modelos:

Modelo	DIC esperado	Interpretación
m1 (base, sin CAR)	Mayor	No captura dependencia espacial
m3_icar	Menor que m1	ICAR absorbe autocorrelación
m4_bym	Similar o menor que m3	BYM añade flexibilidad adicional
m5_ler	Similar a m4	Leroux estima $\rho$, potencialmente más parsimonioso

#Leroux et al. model
m5_leroux <- Diagonal(nrow(aoi.mat), apply(aoi.mat, 1, sum)) - aoi.mat
Cmatrix <- Diagonal(nrow(aoi), 1) -  m5_leroux
max(eigen(Cmatrix)$values) #just to check =1

m5_ler = inla(formula = lands_rec ~ 1 + RainfallDaysmean + elev_mean + slope_mean +
                f(cuenca, elev_mean, model = "iid") +
                f(id, model = "generic1", Cmatrix = Cmatrix),
                offset=log(area), data = as.data.frame(aoi), family ="poisson",
                control.predictor = list(compute = TRUE),
                control.compute = list(dic = TRUE, waic = TRUE))

summary(m5_ler)
m5_ler$summary.random$cuenca
aoi$LEROUX <- m5_ler$summary.fitted[1:nrow(aoi), "mean"]
res_pearson_m5 <- (aoi$lands_rec - m5_ler$summary.fitted.values$mean) / sqrt(m5_ler$summary.fitted.values$mean)
aoi$res_pearson_m5 <- res_pearson_m5
moran_m5 <- moran.mc(x = res_pearson_m5, listw = aoi.listw,nsim = 999, alternative = "greater")

Interpretación del modelo m5_ler (pendiente aleatoria por cuenca + Leroux)#

El modelo m5_ler introduce una variación importante respecto a los anteriores: en lugar de un intercepto aleatorio por cuenca, incluye una pendiente aleatoria para elev_mean por cuenca:

f(cuenca, elev_mean, model = "iid")

Esto significa que el efecto de la elevación media sobre la tasa de deslizamientos no es igual en las tres cuencas — puede ser positivo en una y negativo en otra, dependiendo de la geomorfología local.

Lectura de m5_ler$summary.random$cuenca:

Este objeto devuelve una tabla con tres filas (una por cuenca) y las siguientes columnas:

mean: pendiente aleatoria estimada $\hat{\gamma}_{\text{cuenca}}$ para elev_mean. El efecto total de la elevación en la cuenca $k$ es $\beta_{\text{elev}} + \hat{\gamma}_k$.
sd: incertidumbre de la estimación.
0.025quant, 0.975quant: intervalo de credibilidad del 95%.

Diferencia clave con el intercepto aleatorio:

Efecto aleatorio	Lo que varía por cuenca	Pregunta que responde
Intercepto (`f(cuenca, model="iid")`)	Nivel base de riesgo	¿Hay cuencas sistemáticamente más peligrosas?
Pendiente (`f(cuenca, elev_mean, model="iid")`)	Respuesta a la elevación	¿El efecto de la altitud es igual en todas las cuencas?

Si los intervalos de credibilidad de las tres pendientes son disjuntos, hay evidencia de que el efecto de la elevación es genuinamente distinto entre cuencas — lo que sugiere que procesos geomorfológicos o climáticos locales modulan cómo la altitud se relaciona con la ocurrencia de deslizamientos.

Modelos con pendientes espaciales variables (m6 y m7)#

Los modelos anteriores permitían que el intercepto (nivel base de riesgo) variara entre cuencas o entre subcuencas vecinas. Una extensión natural es permitir que también la pendiente de una covariable varíe espacialmente — es decir, que el efecto de la elevación o la pendiente no sea uniforme en el territorio, sino que cambie de forma continua a lo largo del espacio.

En INLA, esto se implementa añadiendo un segundo campo latente que pondera la covariable:

f(id_slope, slope_mean, model = "besag", graph = cl_graph)

Aquí id_slope es un índice de polígono (igual a id) y slope_mean es el modificador — el efecto latente $\gamma_i$ multiplica el valor de slope_mean en cada subcuenca, generando un coeficiente espacialmente variable $\beta_{\text{slope}} + \gamma_i$.

Intuición geomorfológica: En terrenos de alta elevación y pendiente pronunciada (cordillera), la pendiente puede ser un predictor más fuerte de deslizamientos que en zonas de piedemonte. Permitir que $\gamma_i$ varíe espacialmente captura esta heterogeneidad de manera continua y suavizada por la estructura de vecindad.

Diferencia entre m6 y m7:

Modelo	Variación de pendiente	Escala	Tipo de efecto
m6 (`id_slope ~ besag`)	Continua, polígono a polígono	Local (subcuenca)	Espacialmente estructurado (CAR)
m7 (`cuenca ~ iid`)	Discreta, por cuenca	Regional (Atrato/Cauca/Magdalena)	No estructurado (i.i.d.)

Para añadir coeficientes aleatorios espaciales (p.ej. para la pendiente) en un modelo CAR con INLA, lo más sencillo es definir un segundo campo latente que module justamente esa covariable. La idea es que en la fórmula se incluyan: (i) un campo latente para el intercepto espacial (como ya se hace con f(id, model="besag", graph=graph)), y (ii) un campo latente “pesado” por la covariable cuya pendiente se quiere aleatorizar (por ejemplo slope_mean).

# Primero, crea dos índices idénticos (uno para intercepto y otro para pendiente)
aoi$id_int   <- 1:nrow(aoi)
aoi$id_slope <- aoi$id_int

# Ahora la fórmula:
m6_spatial_slope <- inla(
  lands_rec ~
    1
    + RainfallDaysmean
    + elev_mean
    + slope_mean                  # efecto fijo de la pendiente
    + f(cuenca, model="iid")      # efecto aleatorio no espacial por cuenca
    + f(id_int,   model="besag", graph=cl_graph)              # intercepto espacial
    + f(id_slope, slope_mean,      model="besag", graph=cl_graph),  # pendiente espacial
  offset = log(area),
  family = "poisson",
  data = aoi,
  control.predictor = list(compute = TRUE),
  control.compute   = list(dic = TRUE, waic = TRUE)
)

Si el objetivo fuera que la pendiente de slope_mean variara entre cuencas (no dentro de ellas), basta con usar el índice cuenca en lugar de id_slope. Esto implica solo tres coeficientes de pendiente distintos (uno por cuenca), en contraste con el modelo m6 que estima un coeficiente por subcuenca. Esta es una versión más parsimoniosa, adecuada cuando se tienen pocas cuencas o cuando la heterogeneidad es principalmente regional y no local.

m7_basin_slope <- inla(
  lands_rec ~
    1 +
    RainfallDaysmean +
    elev_mean +
    slope_mean +
    f(cuenca,   model="iid") +                         # intercepto por cuenca
    f(cuenca, slope_mean, model="iid"),                # pendiente aleatoria por cuenca
  offset = log(area),
  family = "poisson",
  data = aoi,
  control.predictor = list(compute=TRUE),
  control.compute   = list(dic=TRUE, waic=TRUE)
)

Modelo complejo: intercepto + pendientes aleatorias + efecto CAR (m_complex)#

El modelo m_complex es la combinación más completa presentada en este capítulo. Incluye simultáneamente:

Efectos fijos globales ($\boldsymbol{\beta}$): pendientes medias de RainfallDaysmean, elev_mean y slope_mean.
Intercepto aleatorio por cuenca (f(basin_id, model="iid")): nivel base de riesgo distinto para cada cuenca.
Pendientes aleatorias por cuenca para cada covariable: el efecto de la precipitación, elevación y pendiente puede ser distinto en Atrato, Cauca y Magdalena.
Efecto espacial CAR a nivel de polígono (f(id, model="besag", graph=cl_graph)): suavizado espacial residual entre subcuencas vecinas, después de controlar todo lo anterior.

La ecuación completa del predictor lineal es:

\[\log(\mu_i) = \beta_0 + \beta_1 x_{1i} + \beta_2 x_{2i} + \beta_3 x_{3i} + u_{\text{c}[i]} + \gamma_{1,\text{c}[i]} x_{1i} + \gamma_{2,\text{c}[i]} x_{2i} + \gamma_{3,\text{c}[i]} x_{3i} + \phi_i + \log(a_i)\]

Cuándo usar este modelo:

Cuando se sospecha que tanto el nivel de riesgo como la sensibilidad a las covariables varía entre grandes unidades geográficas (cuencas).
Cuando aun así persiste autocorrelación espacial entre subcuencas vecinas no explicada por las covariables ni por las diferencias de cuenca.

Advertencia de sobreajuste: Cada nuevo efecto aleatorio añade parámetros y puede dificultar la convergencia de INLA o producir estimaciones inestables. Si el DIC/WAIC no mejora sustancialmente con la versión compleja, preferir el modelo más parsimonioso (m3_icar o m4_bym). El principio de parsimonia es fundamental en el modelado bayesiano jerárquico.

Para incluir a la vez intercepto aleatorio por cuenca, pendientes aleatorias por cuenca para cada covariable, y un efecto CAR espacial a nivel de ID (cada polígono), se puede usar una fórmula como la siguiente en R-INLA:

# índice numérico de cuenca (1,2,3)  
aoi$basin_id <- as.integer(aoi$cuenca)  
# índice numérico de polígono (ID)  
aoi$id        <- 1:nrow(aoi)

m_complex <- inla(
  lands_rec ~
    1
    + RainfallDaysmean          # efecto fijo
    + elev_mean                 # efecto fijo
    + slope_mean                # efecto fijo

    # 1) Intercepto aleatorio por cuenca
    + f(basin_id, model="iid")

    # 2) Pendientes aleatorias por cuenca
    + f(basin_id, RainfallDaysmean, model="iid")
    + f(basin_id, elev_mean,        model="iid")
    + f(basin_id, slope_mean,       model="iid")

    # 3) Efecto espacial CAR a nivel de polígono
    + f(id, model="besag", graph=cl_graph),

  data = aoi,
  family = "poisson",
  offset = log(area),
  control.predictor = list(compute = TRUE),
  control.compute   = list(dic = TRUE, waic = TRUE)
)

Modelos SAR con regímenes espaciales (spatialreg)#

Enfoque alternativo: heterogeneidad como regímenes, no como efectos aleatorios#

Los modelos INLA presentados arriba tratan la heterogeneidad espacial (diferencias entre cuencas) como efectos aleatorios bajo un marco bayesiano: los interceptos o pendientes por cuenca se asumen provenientes de una distribución normal común. Este enfoque es parsimonioso y permite la contracción (shrinkage): cuencas con pocos datos “toman prestada” información de las demás.

La alternativa frecuentista es modelar la heterogeneidad como regímenes espaciales (spatial regimes): se permiten coeficientes completamente distintos para cada grupo de unidades espaciales, sin ninguna restricción de distribución. En spatialreg, esto se hace interactuando las covariables con el factor de cuenca:

y ~ 0 + (x1 + x2 + x3):(cuenca)

El 0 elimina el intercepto global y permite que cada cuenca tenga su propio intercepto implícito a través de la interacción.

Diferencias clave entre ambos enfoques#

Característica	Efectos aleatorios (INLA)	Regímenes espaciales (spatialreg)
Marco estadístico	Bayesiano	Frecuentista
Supuesto sobre cuencas	Distribución normal común	Completamente libre (sin restricciones)
*Contracción (shrinkage)*	Sí — cuencas con pocos datos usan info de otras	No — cada cuenca se estima independientemente
Número de parámetros	Menor (más parsimonioso)	Mayor (un coeficiente por cuenca por variable)
Dependencia espacial	CAR (ICAR, BYM, Leroux)	SAR (lag, error, Mansky, SAC, SLX, SDM)
Interpretación	Distribución posterior de efectos	Coeficientes puntuales por régimen

Modelos SAR implementados con `spatialreg`#

A continuación se ajustan seis variantes de modelos SAR con regímenes por cuenca. Recordando las familias del capítulo 10:

SAR-Lag (lagsarlm): dependencia en $\mathbf{y}$ — el valor de una subcuenca depende del promedio de sus vecinas.
SEM (errorsarlm): dependencia en los errores — los residuos muestran autocorrelación espacial.
Mansky / SAC-mixed (sacsarlm, type="sacmixed"): combina lag en $\mathbf{y}$ y lag en las covariables.
SAC (sacsarlm, type="sac"): dependencia simultánea en $\mathbf{y}$ y en los errores.
SLX (lmSLX): lag solo en las covariables — el efecto de las vecinas actúa a través de $\mathbf{X}$.
SDM (lagsarlm, type="mixed"): lag en $\mathbf{y}$ y en las covariables (Spatial Durbin Model).

A continuación se presenta un modelo SAR que está intentando explicar la variación en la variable dependiente transformada (y_transf) en función de la elevación media, la pendiente media y la precipitación media anual, permitiendo que la relación entre estas variables predictoras y la variable dependiente sea diferente para cada cuenca. Además, el componente SAR del modelo tiene en cuenta la autocorrelación espacial en la variable dependiente, lo que significa que el valor de y_transf en una ubicación está influenciado por los valores de y_transf en sus ubicaciones vecinas (definidas por aoi.listw). La ausencia de un intercepto implica que el modelo está restringido a pasar por el origen para cada cuenca en el espacio de las variables predictoras.

# Regímenes - Modelo de Retardo Espacial (SAR) - cuenca
col.fit9 <- lagsarlm(y_transf ~ 0 + (elev_mean + slope_mean + RainfallDaysmean):(cuenca), data = aoi2,listw = aoi.listw)
summary(col.fit9,Nagelkerke=T)

# Regímenes - Modelo de Error Espacial (SEM) - cuenca
col.fit11 = errorsarlm(y_transf ~ 0 + (elev_mean + slope_mean + RainfallDaysmean):(cuenca), data = aoi2,listw = aoi.listw)
summary(col.fit11,Nagelkerke=T)

# Regímenes - Modelo Mansky (SARAR con restricciones) - cuenca
col.fit14 <- sacsarlm(y_transf ~ 0 + (elev_mean + slope_mean + RainfallDaysmean):(cuenca), data = aoi2,listw = aoi.listw, type="sacmixed")
summary(col.fit14,Nagelkerke=T)

# Regímenes - Modelo SAC (Simultaneous Autoregressive and Conditional Autoregressive) - cuenca
col.fit15 <- sacsarlm(y_transf ~ 0 + (elev_mean + slope_mean + RainfallDaysmean):(cuenca), data = aoi2,listw = aoi.listw, type="sac")
summary(col.fit15)

# Regímenes - Modelo SLX (Spatial Lag of X) - cuenca
col.fit17 = lmSLX(y_transf ~ 0 + (elev_mean + slope_mean + RainfallDaysmean):(cuenca), data = aoi2,listw = aoi.listw)
summary(col.fit17,Nagelkerke=T)
AIC(col.fit17)

# Regímenes - Modelo SDM (Spatial Durbin Model) - cuenca
col.fit19 = lagsarlm(y_transf ~ 0 + (elev_mean + slope_mean + RainfallDaysmean):(cuenca), data = aoi2,listw = aoi.listw,type = "mixed")
summary(col.fit19,Nagelkerke=T)

Interpretación de los modelos SAR con regímenes espaciales#

Lectura del `summary()` en `spatialreg`#

Cada modelo (col.fit9, col.fit11, col.fit14, etc.) produce un resumen que incluye:

Coeficientes por régimen: Con la especificación ~ 0 + (x1 + x2 + x3):(cuenca), el output mostrará coeficientes separados para cada combinación cuenca × variable:

elev_mean:cuencaAtrato     → efecto de la elevación en Atrato
elev_mean:cuencaCauca      → efecto de la elevación en Cauca
elev_mean:cuencaMagdalena  → efecto de la elevación en Magdalena

Para evaluar si el efecto de una variable es significativamente distinto entre cuencas, comparar los intervalos de confianza (p-valores en la tabla Coefficients). Si los tres coeficientes de elev_mean tienen el mismo signo y magnitud similar, la heterogeneidad para esa variable es pequeña.

Parámetro $\rho$ (o $\lambda$) de dependencia espacial:

En lagsarlm (SAR-Lag): Rho = [valor] — si es cercano a 0, la dependencia en $\mathbf{y}$ es débil.
En errorsarlm (SEM): Lambda = [valor] — si es cercano a 0, la autocorrelación en errores es débil.
En sacsarlm (SAC): reporta tanto Rho como Lambda.

Comparación de modelos con AIC / pseudo-R²#

spatialreg reporta AIC y, con summary(modelo, Nagelkerke=T), el pseudo-R² de Nagelkerke — una medida de ajuste global entre 0 y 1. Valores más altos son mejores.

Modelo	Dependencia capturada	Heterogeneidad
`col.fit9` (SAR-Lag)	En $\mathbf{y}$ (retardo)	Por cuenca (regímenes)
`col.fit11` (SEM)	En errores	Por cuenca (regímenes)
`col.fit14` (Mansky)	En $\mathbf{y}$ y en $\mathbf{X}$W	Por cuenca (regímenes)
`col.fit15` (SAC)	En $\mathbf{y}$ y en errores	Por cuenca (regímenes)
`col.fit17` (SLX)	Solo en $\mathbf{X}$W (lag de covariables)	Por cuenca (regímenes)
`col.fit19` (SDM)	En $\mathbf{y}$ y en $\mathbf{X}$W	Por cuenca (regímenes)

Síntesis: ¿INLA-CAR o spatialreg-SAR?#

Al finalizar este capítulo, es útil reflexionar sobre qué enfoque resulta más adecuado para el problema de los deslizamientos en cuencas andinas:

Si la variable respuesta es un conteo (número de deslizamientos) o una tasa (deslizamientos/área), los modelos INLA-Poisson son más apropiados — spatialreg asume distribución gaussiana.
Si hay claras diferencias estructurales entre cuencas y se quiere estimarlas sin restricciones distribucionales, los regímenes de spatialreg son más flexibles.
Si el tamaño de muestra por cuenca es pequeño (pocas subcuencas en alguna cuenca), los efectos aleatorios bayesianos de INLA son más robustos (contracción hacia la media global).
Los modelos INLA son más adecuados para predicción en nuevas ubicaciones (interpolación bayesiana), mientras que los modelos spatialreg son preferibles para inferencia causal sobre los coeficientes individuales.

Actividades#

Comparación de modelos: Construye una tabla resumen comparando los modelos m3_icar, m4_bym y m5_ler usando DIC, WAIC y el Moran I de los residuos de Pearson. ¿Añadir el intercepto aleatorio por cuenca mejora sustancialmente el ajuste respecto al modelo sin heterogeneidad (capítulo 11)?
Efectos aleatorios por cuenca: Del modelo m3_icar, extrae m3_icar$summary.random$cuenca. Grafica los intervalos de credibilidad al 95% de los interceptos aleatorios para las tres cuencas. ¿Son estadísticamente distintos entre sí? ¿Cuál cuenca tiene el mayor riesgo base de deslizamientos después de controlar las covariables?
Pendientes variables vs. fijas: Compara el modelo m3_icar (intercepto aleatorio por cuenca) con el modelo m5_ler (pendiente aleatoria de elev_mean por cuenca). ¿Mejora el DIC? ¿Es el efecto de la elevación realmente distinto entre cuencas? Interpreta los coeficientes de pendiente específicos de cuenca.
Regímenes SAR vs. efectos aleatorios INLA: Para los modelos de spatialreg, compara los coeficientes de slope_mean para cada cuenca (del modelo SAR-Lag con regímenes) con las pendientes aleatorias estimadas por INLA (del modelo m5_ler). ¿Son similares las estimaciones? Discute las diferencias conceptuales y prácticas entre ambos enfoques.

Recursos adicionales#

Documentación oficial#

Lecturas recomendadas#

Gómez-Rubio, V. (2020). Bayesian Inference with INLA. CRC Press.
Moraga, P. (2019). Geospatial Health Data: Modeling and Visualization with R-INLA and Shiny. CRC Press.
Rue, H., Martino, S., & Chopin, N. (2009). Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations. JRSS-B, 71(2), 319–392.
Riebler, A. et al. (2016). An intuitive Bayesian spatial model for disease mapping that accounts for scaling. Statistical Methods in Medical Research, 25(4), 1145–1165. (Modelo BYM2)

Modelos con dependencia y heterogeneidad espacial

Contents

Modelos con dependencia y heterogeneidad espacial#

Marco conceptual: dependencia + heterogeneidad en un mismo modelo#

Dos fuentes de variación espacial#

Progresión de modelos en este capítulo#

Objetivos de aprendizaje#

INLA#

Estructura interna de los modelos INLA con efectos mixtos#

Interpretación del modelo m3_icar (intercepto aleatorio + ICAR)#

Interpretación del modelo m4_bym (intercepto aleatorio + BYM)#

Interpretación del modelo m5_ler (pendiente aleatoria por cuenca + Leroux)#

Modelos con pendientes espaciales variables (m6 y m7)#

Modelo complejo: intercepto + pendientes aleatorias + efecto CAR (m_complex)#

Modelos SAR con regímenes espaciales (spatialreg)#

Enfoque alternativo: heterogeneidad como regímenes, no como efectos aleatorios#

Diferencias clave entre ambos enfoques#

Modelos SAR implementados con `spatialreg`#

Interpretación de los modelos SAR con regímenes espaciales#

Lectura del `summary()` en `spatialreg`#

Comparación de modelos con AIC / pseudo-R²#

Síntesis: ¿INLA-CAR o spatialreg-SAR?#

Actividades#

Recursos adicionales#

Documentación oficial#

Lecturas recomendadas#

Recursos interactivos#

Modelo	Dependencia capturada	Heterogeneidad
`col.fit9` (SAR-Lag)	En \(\mathbf{y}\) (retardo)	Por cuenca (regímenes)
`col.fit11` (SEM)	En errores	Por cuenca (regímenes)
`col.fit14` (Mansky)	En \(\mathbf{y}\) y en \(\mathbf{X}\)W	Por cuenca (regímenes)
`col.fit15` (SAC)	En \(\mathbf{y}\) y en errores	Por cuenca (regímenes)
`col.fit17` (SLX)	Solo en \(\mathbf{X}\)W (lag de covariables)	Por cuenca (regímenes)
`col.fit19` (SDM)	En \(\mathbf{y}\) y en \(\mathbf{X}\)W	Por cuenca (regímenes)

Modelos con dependencia y heterogeneidad espacial

Contents

Modelos con dependencia y heterogeneidad espacial#

Marco conceptual: dependencia + heterogeneidad en un mismo modelo#

Dos fuentes de variación espacial#

Progresión de modelos en este capítulo#

Objetivos de aprendizaje#

INLA#

Estructura interna de los modelos INLA con efectos mixtos#

Interpretación del modelo m3_icar (intercepto aleatorio + ICAR)#

Interpretación del modelo m4_bym (intercepto aleatorio + BYM)#

Interpretación del modelo m5_ler (pendiente aleatoria por cuenca + Leroux)#

Modelos con pendientes espaciales variables (m6 y m7)#

Modelo complejo: intercepto + pendientes aleatorias + efecto CAR (m_complex)#

Modelos SAR con regímenes espaciales (spatialreg)#

Enfoque alternativo: heterogeneidad como regímenes, no como efectos aleatorios#

Diferencias clave entre ambos enfoques#

Modelos SAR implementados con spatialreg#

Interpretación de los modelos SAR con regímenes espaciales#

Lectura del summary() en spatialreg#

Comparación de modelos con AIC / pseudo-R²#

Síntesis: ¿INLA-CAR o spatialreg-SAR?#

Actividades#

Recursos adicionales#

Documentación oficial#

Lecturas recomendadas#

Recursos interactivos#

Modelos SAR implementados con `spatialreg`#

Lectura del `summary()` en `spatialreg`#