Modelos de regresión para dependencia espacial tipo CAR

Modelos de regresión para dependencia espacial tipo CAR#

Los modelos CAR con Campos Aleatorios de Markov (MRF) se utilizan ampliamente en estadísticas espaciales para modelar datos observados, así como variables latentes y efectos aleatorios que varían espacialmente. Las estructuras CAR se implementan en modelos jerárquicos con efectos aleatorios latentes (Schmidt 2014). CAR es un tipo de Campo Aleatorio de Markov (MRF) (Besag 1991), que es el enfoque dominante para analizar datos espaciales discretos dentro de un marco jerárquico. Un proceso espacial se considera que tiene propiedades de Markov si el estado futuro esperado de un parámetro depende únicamente del estado adyacente inmediato, haciendo que el estado futuro sea condicionalmente independiente. Extender esto a múltiples efectos aleatorios da como resultado un MRF (Clifford 1990, Rue 2005). La forma más común de MRF representa una conectividad de vecindad de primer orden en términos de contigüidad, donde las observaciones que comparten un límite se consideran vecinas.

Los efectos aleatorios latentes (o simplemente efectos aleatorios espaciales) ofrecen una forma diferente y muy flexible de capturar la dependencia espacial, a menudo representando la variabilidad espacial no observada o no explicada por las covariables fijas en tu modelo.

Los efectos aleatorios latentes pueden ser de dos tipos principales en el contexto espacial:

Efecto Aleatorio Espacial Estructurado (Structured Spatial Random Effect): Este componente captura la dependencia espacial “suave” o de “vecindad”. Modela la correlación entre observaciones cercanas, asumiendo que las ubicaciones más próximas son más similares.
Efecto Aleatorio Espacial No Estructurado (Unstructured Spatial Random Effect): Este es un componente aleatorio independiente e idénticamente distribuido (i.i.d.) para cada ubicación. Captura la variabilidad “ruidosa” que no tiene una estructura espacial discernible. Podría pensarse como una variabilidad aleatoria intrínseca a cada observación o ubicación.

La idea es que la combinación de estos efectos aleatorios (junto con los efectos fijos de las covariables) puede explicar completamente la variación en tu variable dependiente, incluyendo la dependencia espacial.

Varios modelos han sido propuestos dentro de esta clase general de estructuras CAR, incluidos los modelos intrínsecos y de convolución (Besag 1991), así como alternativas como la propuesta por Leroux 2000. El modelo CAR intrínseco (ICAR) es el CAR más simple, que asume un efecto aleatorio autorregresivo espacial para abordar las asociaciones entre regiones vecinas (Besag 1991). La expectativa condicional es igual a la media de los efectos aleatorios en las unidades de mapeo de terreno vecinas, mientras que la varianza condicional es inversamente proporcional al número de sus vecinos. Este modelo representa estructuras de correlación espacial fuertes y puede no ser adecuado si los datos están débilmente correlacionados.

Besag 1991 propuso el modelo de convolución, o modelo Besag-York-Mollié (BYM), que combina dos efectos aleatorios latentes: un efecto latente ICAR ($\rho$) y un efecto latente Gaussiano i.i.d. ($\sigma$). Sin embargo, los dos componentes de efectos aleatorios separados no se pueden identificar individualmente, y solo se puede identificar su suma.

Leroux 2000 propuso utilizar un único efecto aleatorio en su lugar, que incluye un parámetro ($\rho$) para medir el nivel de correlación espacial entre las unidades de mapeo de terreno, tomando valores en el intervalo unitario [0-1]. Este modelo es una variación de los modelos BYM e ICAR. Al igual que el modelo ICAR, tiene un parámetro de efecto aleatorio espacial para cada observacion. Sin embargo, la distribución condicional incorpora características tanto de efectos aleatorios espaciales estructurados como no estructurados (del modelo BYM) en un solo parámetro ($\rho$). El modelo de Leroux generaliza el modelo ICAR y el modelo independiente (es decir, un modelo sin ningún efecto aleatorio espacial estructurado). Cuando $\rho=1$, se recupera el modelo ICAR; cuando $\rho=0$, se recupera el modelo independiente. Así, el modelo de Leroux busca equilibrar estos dos modelos estimando el valor de $\rho$.

Objetivos de aprendizaje#

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

Diferenciar conceptualmente entre los enfoques SAR (especificación conjunta) y CAR (especificación condicional, propiedad de Markov).
Comprender la estructura de los modelos ICAR, BYM y Leroux como efectos aleatorios latentes en un marco bayesiano.
Ajustar modelos CAR con CARBayes en R y evaluar la convergencia de las cadenas MCMC.
Ajustar modelos ICAR, BYM y Leroux con INLA en R aprovechando la aproximación de Laplace anidada.
Comparar modelos mediante DIC/WAIC e interpretar los hiperparámetros de precisión.

Requisitos previos: 10_SAR — modelos SAR; nociones de estadística bayesiana (distribuciones prior/posterior); R básico.

Diferencia conceptual entre SAR y CAR#

Antes de adentrarnos en los modelos CAR, es fundamental entender en qué se diferencian de los modelos SAR (Simultaneous Autoregressive) presentados en el capítulo anterior. Aunque ambos modelan dependencia espacial, lo hacen de formas matemáticamente distintas.

Especificación: global (SAR) vs. condicional (CAR)#

Característica	SAR	CAR
Especificación	Distribución conjunta	Distribución condicional
Perspectiva	Global: todo el sistema simultáneamente	Local: unidad por unidad dado su entorno
Propiedad de Markov	No	Sí — depende solo de vecinas directas
Parámetro de dependencia	$\rho \in (-1, 1)$	$\lambda$ (precisión condicional)
Enfoque estadístico	Frecuentista (máx. verosimilitud)	Bayesiano jerárquico (MCMC o INLA)
Implementación en R	`spatialreg` (`lagsarlm`, `errorsarlm`)	`CARBayes`, `INLA`

SAR: sistema simultáneo#

Un modelo SAR-lag define la variable dependiente como función de todo el sistema a la vez:

\[\mathbf{y} = \rho \mathbf{W}\mathbf{y} + \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon}\]

Aquí $\rho$ captura cuánto del valor de $y_i$ se explica por el promedio ponderado de sus vecinas $\sum_j w_{ij} y_j$. La solución requiere invertir la matriz $(I - \rho W)$, lo que implica que cada unidad influye —aunque indirectamente— en todas las demás.

CAR: propiedad de Markov espacial#

Los modelos CAR especifican la distribución de cada efecto aleatorio $\phi_i$ condicionada en todas las demás unidades:

\[[\phi_i \mid \phi_{j \neq i}] \sim \mathcal{N}\!\left(\frac{\sum_{j \in \mathcal{N}(i)} \phi_j}{n_i},\ \frac{\sigma^2}{n_i}\right)\]

donde $\mathcal{N}(i)$ es el conjunto de vecinas directas de $i$ y $n_i$ su cardinalidad. Esta es la propiedad de Markov espacial: el estado esperado de $\phi_i$ es simplemente el promedio de sus vecinas inmediatas, con una varianza inversamente proporcional al número de vecinos (áreas con más vecinos tienen estimaciones más precisas). No hay influencia de unidades lejanas más allá de la cadena de vecindad.

¿Cuándo preferir uno u otro?#

SAR es adecuado cuando se quiere modelar la dependencia espacial directamente en la variable observada, en un marco de estimación frecuentista con interpretación directa de coeficientes.
CAR es preferible en modelos bayesianos jerárquicos, donde la dependencia espacial se introduce como un efecto aleatorio latente $\phi_i$ que suaviza la variabilidad residual entre unidades vecinas, capturando variación espacial no explicada por las covariables fijas. Este es el uso que desarrollaremos a lo largo de este capítulo.

Modelos con la librería CARBayes en R#

Existen diferentes librerías en R para trabajar modelos espaciales. A continuación se presenta una de ellas.

library(spBayes)
library(maps)
library(RANN)
library(gjam)
library(CARBayes)
library(CARBayesdata)
library(mgcv)

#### Set up a square lattice region
m <- 12
xEast  <- 1:m
xNorth <- 1:m
grid   <- expand.grid(xEast, xNorth)
n      <- nrow(grid)
plot( NULL, xlim = c(0, m), ylim = c(0, m), xlab='East', ylab='North' )
abline(v=grid[,1], h=grid[,2])
text(grid[,1] - .5, grid[,2] - .5, 1:n, cex=.8)

Se construye una matriz de vecindad (W).

D <- W <- as.matrix(dist(grid))
W[W != 1] <- 0 

Q <- 3
x <- matrix( rnorm(Q*n), n, Q )
x[,1] <- 1
x2    <- x[,2]
x3    <- x[,3]
beta  <- matrix( rnorm(Q), Q, 1)
sigma <- .1

form <- as.formula(y ~ x2 + x3)

gaussianModel <- S.CARleroux(formula = form, family  = 'gaussian', W = W, 
                             burnin = 20000, n.sample = 100000, thin = 10, verbose = F)
gaussianModel

#random effect
fv <- gaussianModel$fitted.values
mf <- min(fv)
cc <- fv - mf
ss <- seq(0, max(cc), length.out=10)
cc <- findInterval(cc, ss)

colM <- colorRampPalette( c("red","orange","blue"))
colm <- colM(10)

symbols(x=grid[,1], y=grid[,2], squares = cc*0+1, bg=colm[cc],
        fg=colm[cc],inches=F, xlab='East', ylab='North')

lambda <- exp(x%*%beta + phi[,2] + rnorm(n, 0, sigma))
y <- rpois(n, lambda)

poissonModel <- S.CARbym(formula=form, family="poisson",
                         W=W, burnin=20000, n.sample=100000, thin=10, verbose=F)
poissonModel

Modelos con la librería INLA en R#

La librería INLA (Integrated Nested Laplace Approximation) es una herramienta poderosa para ajustar modelos bayesianos latentes gaussianos (LGMs), ampliamente utilizados en estadística espacial. A diferencia de los métodos MCMC tradicionales, INLA ofrece una aproximación determinista, mucho más rápida y eficiente, especialmente útil en contextos con grandes volúmenes de datos o estructuras espaciales complejas. INLA no elimina la incertidumbre bayesiana. Más bien, lo que hace es utilizar aproximaciones tipo Laplace para calcular directamente la media, desviación estándar y cuantiles de las distribuciones posteriores. Estas aproximaciones son deterministas, por lo tanto reproducibles. Los resultados siguen siendo probabilísticos: por ejemplo, se puede obtener una distribución posterior para cada parámetro o predicción, pero no mediante muestreo aleatorio.

Para ilustrar cómo se ajustan los modelos espaciales con INLA, se utilizará el conjunto de datos de leucemia de Nueva York. Éste ha sido ampliamente analizado en la literatura (ver, por ejemplo, Waller y Gotway, 2004) y está disponible en el paquete DClusterm. El conjunto de datos registra una serie de casos de leucemia en el norte del estado de Nueva York a nivel de distrito censal. Algunas de las variables en el conjunto de datos son:

Casos: Número de casos de leucemia en el período 1978-1982.
POP8: Población en 1980.
PCTOWNHOME: Proporción de personas que son propietarias de su vivienda.
PCTAGE65P: Proporción de personas de 65 años o más.
AVGIDIST: Distancia inversa promedio al sitio de tricloroetileno (TCE) más cercano.

Tenga en cuenta que el interés se centra en la exposición al TCE, utilizando AVGIDIST como un indicador de exposición. Las variables PCTOWNHOME y PCTAGE65P actuarán como posibles factores de confusión que deben incluirse en el modelo. Sin embargo, no lo haremos aquí porque queremos probar cómo los efectos espaciales latentes capturan la variación espacial residual.

El conjunto de datos se puede cargar de la siguiente manera:

library(spdep)
library(DClusterm)
data(NY8)

Dado que el interés se centra en estudiar el riesgo de leucemia en el norte del estado de Nueva York, primero se calcula el número esperado de casos. Esto se hace calculando la tasa de mortalidad general (total de casos dividido por la población total) y multiplicándola por la población:

rate <- sum(NY8$Cases) / sum(NY8$POP8)
NY8$Expected <- NY8$POP8 * rate

Una vez que se obtiene el número esperado de casos, se puede obtener una estimación aproximada del riesgo con la razón de mortalidad estandarizada (SMR), que se calcula como el número de casos observados dividido por el número de casos esperados.

NY8$SMR <- NY8$Cases / NY8$Expected

En epidemiología, es importante producir mapas para mostrar la distribución espacial del riesgo relativo. En este ejemplo, nos centraremos en la ciudad de Syracuse para reducir el tiempo de cómputo necesario para producir el mapa. Por lo tanto, creamos un índice con las áreas de la ciudad de Syracuse:

# Subset Syracuse city
syracuse <- which(NY8$AREANAME == "Syracuse city")

Un mapa de enfermedades se puede crear simplemente con la función spplot (del paquete sp):

library(viridis)
spplot(NY8[syracuse, ], "SMR", #at = c(0.6, 0.9801, 1.055, 1.087, 1.125, 13),
   col.regions = rev(magma(16))) #gray.colors(16, 0.9, 0.4))

Modelo gaussiano clásico#

El primer modelo que consideraremos es un modelo Gausiano sin efectos aleatorios latentes, ya que proporcionará una línea base para comparar con otros modelos. Para ajustar el modelo con INLA, se utiliza la función inla:

library(INLA)

m1 <- inla(Cases ~ 1 + AVGIDIST,
   data = as.data.frame(NY8),
   family = "Gaussian",verbose=T,
   E = NY8$Expected, control.predictor = list(compute = TRUE),
   control.compute = list(dic = TRUE, waic = TRUE))

summary(m1)
NY8$FIXED.EFF <- m1$summary.fitted[, "mean"]

Modelo Gaussiano con efectos aleatorios NO estructurados#

Se pueden agregar efectos aleatorios latentes al modelo para tener en cuenta la sobredispersión incluyendo efectos aleatorios Gaussianos i.i.d. en el predictor lineal. Para ajustar el modelo con INLA, primero se crea un índice para identificar los efectos aleatorios (ID). Los efectos aleatorios latentes se especifican con la función f en INLA:

NY8$ID <- 1:nrow(NY8)
m2 <- inla(Cases ~ 1 + AVGIDIST + f(ID, model = "iid"),
  data = as.data.frame(NY8), family = "Gaussian", 
  E = NY8$Expected,
  control.predictor = list(compute = TRUE),
  control.compute = list(dic = TRUE, waic = TRUE))

summary(m2)
NY8$IID.EFF <- m2$summary.fitted[, "mean"]

spplot(NY8[syracuse, ], c("SMR", "FIXED.EFF", "IID.EFF"),
  col.regions = rev(magma(16)))

Modelo Gaussiano con efectos aleatorios estructurados#

Tenemos observaciones $y = {y_i}_{i=1}^n$ donde n es el número de áreas. A y se le asigna una distribución multivariante que tiene en cuenta la dependencia espacial. Una forma común de describir la proximidad espacial en datos discretos es mediante una matriz de adyacencia W. El elemento $W_{i,j}$ es distinto de cero si las áreas i y j son vecinas. Por lo general, dos áreas son vecinas si comparten un límite común.

La matriz de adyacencia se puede calcular utilizando la función poly2nb en el paquete spdep. Esta función considerará dos áreas como vecinas si sus bordes se tocan al menos en un punto (es decir, adyacencia de “Queen”):

NY8.nb <- poly2nb(NY8)

Esto devolverá un objeto nb con la definición de la estructura del vecindario.

NY8.nb

Además, estos objetos creados con la función nb se pueden visualizar gráficamente cuando se conoce información adicional sobre las áreas, como la ubicación de su centro.

plot(NY8) 
plot(NY8.nb, coordinates(NY8), add = TRUE, pch = ".", col = "gray")

Modelo ICAR (Intrinsic CAR)#

El modelo ICAR (Intrinsic Conditional Autoregressive) es la forma más simple de estructura CAR. Su distribución conjunta es impropia (no integra a 1), lo que significa que no puede usarse como una distribución de probabilidad absoluta, sino como una prior impropia en un modelo bayesiano. A pesar de ello, el modelo posterior sí es propio siempre que los datos sean informativos.

Intuición: El ICAR impone que el efecto aleatorio $\phi_i$ tenga como valor esperado el promedio simple de sus vecinas:

\[[\phi_i \mid \boldsymbol{\phi}_{-i}] \sim \mathcal{N}\!\left(\bar{\phi}_{\mathcal{N}(i)},\ \frac{\sigma^2}{n_i}\right)\]

Esto genera un suavizado espacial fuerte: áreas vecinas tenderán a tener efectos similares. El modelo no tiene intercepto propio (la suma de todos los efectos debe ser cero, restricción de identificabilidad), por lo que captura solo variación espacial relativa, no niveles absolutos.

Cuándo usarlo: Es adecuado cuando se espera que los datos presenten una autocorrelación espacial fuerte y se quiere capturar patrones de variación espacial suave sin modelar explícitamente su intensidad.

El primer ejemplo se basará en la especificación ICAR. Tenga en cuenta que el efecto latente espacial se define usando la función f. Esto requerirá un índice para identificar los efectos aleatorios en cada área, el tipo de modelo y la matriz de adyacencia. Para ello, se utilizará una matriz dispersa.

# Create sparse adjacency matrix
NY8.mat <- as(nb2mat(NY8.nb, style = "B"), "Matrix")

# Fit model
m.icar <- inla(Cases ~ 1 +  AVGIDIST + 
    f(ID, model = "besag", graph = NY8.mat),
  data = as.data.frame(NY8), E = NY8$Expected, family ="Gaussian",
  control.predictor = list(compute = TRUE),
  control.compute = list(dic = TRUE, waic = TRUE))

summary(m.icar)
NY8$ICAR <- m.icar$summary.fitted.values[, "mean"]

Modelo BYM (Besag, York y Mollié)#

El modelo BYM fue propuesto por Besag, York y Mollié en 1991 y se convirtió en el estándar de facto para el análisis de enfermedades en áreas pequeñas (small area estimation). Su innovación central es descomponer la variabilidad espacial en dos fuentes independientes:

\[\phi_i = u_i + v_i\]

donde:

$u_i$ es un efecto espacialmente estructurado (ICAR): captura la variación que comparten unidades vecinas (autocorrelación espacial).
$v_i$ es un efecto no estructurado i.i.d.: captura variabilidad intrínseca de cada unidad que no sigue ningún patrón espacial (heterogeneidad pura o sobredispersión).

Intuición: Imaginemos que la tasa de deslizamientos en una subcuenca es más alta de lo esperado. Hay dos explicaciones posibles: (a) sus vecinas también tienen tasas altas, sugiriendo un proceso geográfico compartido ($u_i$); o (b) algo particular de esa subcuenca no relacionado con la vecindad la hace distinta ($v_i$). El BYM modela ambas explicaciones simultáneamente.

Limitación: Los componentes $u_i$ y $v_i$ no son individualmente identificables — solo su suma $\phi_i = u_i + v_i$ puede estimarse de forma robusta. Esto motivó el posterior modelo BYM2 (Riebler et al. 2016), que reparametriza la descomposición para hacer los componentes identificables.

El modelo de Besag, York y Mollié incluye dos efectos aleatorios latentes: un efecto latente ICAR y un efecto latente Gaussiano i.i.d. No es necesario definir estos dos efectos latentes por separado si se establece el argumento model como "bym" cuando se define el efecto aleatorio latente con la función f.

m.bym = inla(Cases ~ 1 +  AVGIDIST + 
    f(ID, model = "bym", graph = NY8.mat),
  data = as.data.frame(NY8), E = NY8$Expected, family ="Gaussian",
  control.predictor = list(compute = TRUE),
  control.compute = list(dic = TRUE, waic = TRUE))

summary(m.bym)
NY8$BYM <- m.bym$summary.fitted.values[, "mean"]

Modelo de Leroux et al. (CAR propio)#

El modelo de Leroux et al. (2000) surge como una solución elegante al problema de la distribución impropia del ICAR. En lugar de separar efectos estructurados y no estructurados como hace el BYM, propone un único efecto aleatorio que interpola entre los dos extremos mediante un solo parámetro $\rho \in [0, 1]$:

\[[\phi_i \mid \boldsymbol{\phi}_{-i}] \sim \mathcal{N}\!\left(\frac{\rho \sum_{j \in \mathcal{N}(i)} \phi_j}{\rho n_i + (1 - \rho)},\ \frac{\sigma^2}{\rho n_i + (1 - \rho)}\right)\]

Interpretación de $\rho$:

$\rho = 0$: el efecto es puramente i.i.d. (sin estructura espacial) — equivale a un modelo independiente.
$\rho = 1$: el efecto es puramente ICAR (estructura espacial máxima) — recupera el modelo intrínseco.
$0 < \rho < 1$: el modelo interpola entre ambos extremos; $\rho$ se estima a partir de los datos.

Ventaja clave: La distribución conjunta del modelo de Leroux es propia para cualquier $\rho \in [0, 1)$, lo que lo hace más manejable computacionalmente que el ICAR. Además, el valor estimado de $\rho$ proporciona una medida directa de la fuerza de la dependencia espacial en los datos — algo que ni el ICAR ni el BYM ofrecen directamente.

Implementación en INLA: Como INLA no implementa directamente el modelo de Leroux, se aproxima mediante el efecto latente generic1, que utiliza la siguiente matriz de precisión:

\[\Sigma^{-1} = \frac{1}{\tau} \left(I_n - \rho \cdot \lambda_{\max} \cdot C\right), \quad \rho \in [0, 1)\]

donde $C = I_n - M$ y $M = D - W$ es la matriz laplaciana del grafo de vecindad ($D$ es la matriz diagonal de grados). Se puede verificar que $\lambda_{\max}(C) = 1$, por lo que la construcción se simplifica.

Este modelo se define utilizando una “mezcla” de matrices (modelo de Leroux et al.) para definir la matriz de precisión del efecto latente. Para ajustar el modelo, el primer paso es crear la matriz $M$.

ICARmatrix <- Diagonal(nrow(NY8.mat), apply(NY8.mat, 1, sum)) - NY8.mat
Cmatrix <- Diagonal(nrow(NY8), 1) -  ICARmatrix

El modelo se ajusta como de costumbre con la función inla. Tenga en cuenta que la matriz C se pasa a la función f usando el argumento Cmatrix:

m.ler = inla(Cases ~ 1 +  AVGIDIST +
    f(ID, model = "generic1", Cmatrix = Cmatrix),
  data = as.data.frame(NY8), E = NY8$Expected, family ="Gaussian",
  control.predictor = list(compute = TRUE),
  control.compute = list(dic = TRUE, waic = TRUE))
summary(m.ler)
NY8$LEROUX <- m.ler$summary.fitted.values[, "mean"]

Ejemplo con datos de cuencas de la implementación de modelos CAR en R en CARBayes#

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)

Inicialmente se crea una matriz de pesos espaciales binaria basada en la contigüidad (tipo “queen”). tyle="B" especifica que el estilo de los pesos es binario. Similar a la matriz W.mat, el objeto W.list representará las relaciones de vecindad con pesos de 1 para los vecinos y 0 para los no vecinos. La diferencia principal con W.mat es la estructura de datos subyacente; W.list es generalmente más eficiente para cálculos con matrices de pesos espaciales dispersas (donde la mayoría de los elementos son cero). nb2listw() convierte la lista de vecinos (W.nb) en un objeto de lista de pesos espaciales (listw), que es un formato más flexible que a menudo se utiliza en las funciones de modelado espacial en spdep y spatialreg.

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

Inicialmente se presenta un modelo de regresión de Poisson bayesiano SIN dependencia espacial y evalúa la autocorrelación espacial restante en los residuos.

#Standard Poisson model in CARBayes
m1_carbayes <- S.glm(formula=lands_rec ~ 1 + RainfallDaysmean + elev_mean + slope_mean + offset(log(area)),
                    data=aoi, family="poisson", 
                    burnin=100000, n.sample=300000, thin=100, n.chains=3, n.cores=3)


m1_carbayes
m1_carbayes$modelfit
aoi_sp$carm1 <- m1_carbayes$fitted.values
moran.mc(x=residuals(m1_carbayes), listw=W.list, nsim=1000)

A continuación se presenta un modelo de regresión de Poisson bayesiano CON dependencia espacial tipo Leroux propio con un nivel moderado de dependencia espacial inicial (rho=0.5), examina la convergencia y los resultados del modelo, obtiene los valores ajustados y los residuos, y evalúa la autocorrelación espacial restante en los residuos.

el término offset(log(areaa)) especifica un término de offset en la regresión de Poisson. En modelos de conteo como la regresión de Poisson, un offset permite modelar tasas en lugar de solo conteos. Por ejemplo, si estás modelando el número de eventos en diferentes áreas, el área en sí misma podría influir en el número de eventos. Al incluir log(area) como un offset, esencialmente estás modelando la tasa de deslizamientos por unidad de área.

rho=0.5 establece el valor inicial para el parámetro de dependencia espacial (rho) en el modelo. Rho típicamente varía entre 0 y 1 e indica la fuerza de la autocorrelación espacial. Un valor de 0 significa que no hay autocorrelación espacial, mientras que un valor más cercano a 1 indica una fuerte autocorrelación espacial. Cuando rho se establece en 1, el modelo CAR Leroux se convierte en un modelo Intrinsic Conditional Autoregressive (ICAR).

plot( cariid$samples$rho, bty = 'n' ) grafica las muestras MCMC para el parámetro de dependencia espacial rho. Examinar el gráfico de traza de las muestras MCMC para rho ayuda a evaluar si la cadena de Markov ha convergido.

#Leroux=proper
cariid <- S.CARleroux(formula=lands_rec ~ 1 + RainfallDaysmean + elev_mean + slope_mean + offset(log(area)),
                      data=aoi, family="poisson", W=W.mat,rho=0.5,
                      burnin=100000, n.sample=300000, thin=100, n.chains=3, n.cores=3)
cariid
cariid$modelfit
summary(cariid$samples)
summary.beta <- summary(cariid$samples$beta, quantiles=c(0.025, 0.975))
beta.mean <- summary.beta$statistics[ ,"Mean"]
beta.ci <- summary.beta$quantiles
beta.results <- cbind(beta.mean, beta.ci)
rownames(beta.results) <- colnames(cariid$X)
round(beta.results, 5)
aoi_sp$cariid <- cariid$fitted.values
aoi_sp$cariid_res <- cariid$residuals
moran.mc(x=residuals(cariid), listw=W.list, nsim=1000)

A continuación se presenta un modelo ICAR. A diferencia de los modelos CAR “propios” (donde rho < 1), la distribución conjunta de un modelo ICAR es “impropia”, lo que significa que no se integra a 1. Sin embargo, sigue siendo útil para modelar la dependencia espacial, especialmente en contextos donde solo se conocen las diferencias relativas entre las áreas. Los modelos ICAR son particularmente útiles para modelar el suavizado espacial. Imponen una estructura donde las unidades vecinas tienen distribuciones posteriores que tienden a estar cerca unas de otras. Esto es útil cuando se espera que las áreas geográficamente cercanas tengan tasas o niveles de la variable de respuesta similares debido a factores espaciales no medidos. Los efectos espaciales en un modelo ICAR implican una restricción de suma a cero. Esto significa que la suma de los efectos espaciales a través de todas las unidades espaciales es cero.

Cuando decimos que la distribución conjunta de un modelo ICAR es “impropia” y “no se integra a 1”, significa que la función matemática que describe la probabilidad de todas las posibles configuraciones de los efectos espaciales no suma (o integra) a uno sobre todo el espacio de posibles valores. Por lo que técnicamente, no podemos interpretar directamente la función como una distribución de probabilidad sobre la que se puedan calcular probabilidades absolutas. Los modelos ICAR se centran en modelar las diferencias en el efecto espacial entre las áreas vecinas. La distribución impropia refleja esta idea de que solo las diferencias relativas están bien definidas, no los valores absolutos de los efectos espaciales.

La restricción de suma a cero es necesaria para la identificabilidad del modelo. Sin ella, habría una infinidad de posibles conjuntos de efectos espaciales que producirían el mismo patrón de autocorrelación espacial. Al imponer la restricción de suma a cero, se fija un punto de referencia y se asegura que haya una única solución para los efectos espaciales que mejor se ajuste a los datos. La restricción de suma a cero centra los efectos espaciales alrededor de cero. Los valores positivos de $ϕ_i$ indican que la unidad espacial i tiene un efecto mayor que el promedio, mientras que los valores negativos indican un efecto menor que el promedio, después de tener en cuenta las covariables. Al centrar los efectos espaciales, los coeficientes de las otras covariables en el modelo pueden interpretarse como el efecto de esas variables promediando la influencia espacial. El efecto espacial representa la desviación de este patrón promedio debido a la ubicación geográfica.

#Leroux=ICAR
caricar <- S.CARleroux(formula=lands_rec ~ 1 + RainfallDaysmean + elev_mean + slope_mean + offset(log(area)),
                       data=aoi, family="poisson", W=W.mat,rho=1,
                       burnin=100000, n.sample=300000, thin=100, n.chains=3, n.cores=3)
caricar
caricar$modelfit
summary(caricar$samples)
summary.beta <- summary(caricar$samples$beta, quantiles=c(0.025, 0.975))
beta.mean <- summary.beta$statistics[ ,"Mean"]
beta.ci <- summary.beta$quantiles
beta.results <- cbind(beta.mean, beta.ci)
rownames(beta.results) <- colnames(caricar$X)
round(beta.results, 5)
aoi_sp$caricar <- caricar$fitted.values
aoi_sp$caricar_res <- caricar$residuals
moran.mc(x=residuals(caricar), listw=W.list, nsim=1000)

En el siguiente caso el parámetro rho no está especificado directamente en la función S.CARleroux. Cuando rho no se especifica, la función estima como parte del proceso de muestreo MCMC. En el caso de S.CARleroux sin especificar rho, el modelo permite que la fuerza de la dependencia espacial (autocorrelación) sea estimada a partir de los datos dentro del rango (0, 1).

```r
carleroux <- S.CARleroux(formula=lands_rec ~ 1 + RainfallDaysmean + elev_mean + slope_mean + offset(log(area)),
                  data=aoi, family="poisson", W=W.mat,
                  burnin=100000, n.sample=300000, thin=100, n.chains=3, n.cores=3)
carleroux
carleroux$modelfit
plot( carleroux$samples$rho, bty = 'n' )
summary(carleroux$samples)
summary.beta <- summary(carleroux$samples$beta, quantiles=c(0.025, 0.975))
beta.mean <- summary.beta$statistics[ ,"Mean"]
beta.ci <- summary.beta$quantiles
beta.results <- cbind(beta.mean, beta.ci)
rownames(beta.results) <- colnames(carleroux$X)
round(beta.results, 5)
aoi_sp$carleroux <- carleroux$fitted.values
aoi$carleroux_res <- carleroux$residuals
moran.mc(x=residuals(carleroux), listw=W.list, nsim=1000)
```

A continuación se implementa el modelo BYM (Besag, York, y Mollié). La principal diferencia con los anteriores modelos radica en la estructura de los efectos aleatorios espaciales. El modelo BYM descompone la variación espacial en dos componentes:

Un componente espacialmente estructurado (o correlacionado): Este componente modela la dependencia espacial, es decir, la tendencia de las unidades vecinas a tener valores similares. Al igual que el modelo CAR Leroux, utiliza la matriz de pesos espaciales (W=W.mat) para definir las relaciones de vecindad.
Un componente espacialmente no estructurado (o no correlacionado): Este componente modela la variabilidad espacial aleatoria o idiosincrásica que no está explicada por la estructura espacial definida por la matriz de pesos. Se asume que este componente sigue una distribución normal independiente para cada unidad espacial.

El componente no estructurado del BYM puede ayudar a modelar la sobredispersión en los datos que no está relacionada con la autocorrelación espacial. La sobredispersión ocurre cuando la varianza de la variable de respuesta es mayor que la media. El modelo BYM ofrece mayor flexibilidad al permitir la existencia de variación espacial que no sigue el patrón de vecindad definido por la matriz W.

#BYM
carbym <- S.CARbym(formula=lands_rec ~ 1 + RainfallDaysmean + elev_mean + slope_mean + offset(log(area)),
                   data=aoi, family="poisson", W=W.mat,
                   burnin=100000, n.sample=300000, thin=100, n.chains=3, n.cores=3)
carbym
carbym$modelfit
plot( carbym$samples$rho, bty = 'n' )
summary(carbym$samples)
summary.beta <- summary(carbym$samples$beta, quantiles=c(0.025, 0.975))
beta.mean <- summary.beta$statistics[ ,"Mean"]
beta.ci <- summary.beta$quantiles
beta.results <- cbind(beta.mean, beta.ci)
rownames(beta.results) <- colnames(carbym$X)
round(beta.results, 5)
aoi$carbym <- carbym$fitted.values
aoi$carbym_res <- carbym$residuals
moran.mc(x=residuals(carbym), listw=W.list, nsim=1000)

Interpretación y comparación de los modelos CARBayes#

Una vez ajustados los modelos (m1_carbayes, cariid, caricar, carleroux, carbym), la comparación y evaluación debe hacerse en varios niveles:

1. Diagnóstico de autocorrelación espacial residual#

El test de Moran sobre los residuos (moran.mc(residuals(modelo), ...)) es el primer indicador de si el modelo ha logrado capturar la dependencia espacial. Se espera que:

El modelo base (m1_carbayes, sin estructura CAR) muestre un Moran I significativo → confirma la presencia de dependencia espacial.
Los modelos CAR (Leroux, ICAR, BYM) reduzcan el Moran I de los residuos hasta que no sea estadísticamente significativo (p > 0.05).

Si los residuos siguen mostrando autocorrelación significativa después de incluir el efecto CAR, puede indicar que la matriz de vecindad $W$ no captura adecuadamente la estructura espacial relevante, o que la dependencia es de orden mayor (vecinas de vecinas).

2. Ajuste del modelo: `modelfit`#

El objeto $modelfit de cada modelo CARBayes reporta:

DIC (Deviance Information Criterion): penaliza la complejidad. Valores más bajos indican mejor ajuste (similar al AIC).
WAIC (Widely Applicable AIC): alternativa bayesiana al DIC, más robusta.
Deviance (mean): desviación promedio del modelo.

Para comparar modelos: un modelo con DIC/WAIC menor es preferible, pero una diferencia < 5 unidades generalmente no es sustancial.

3. Convergencia MCMC: cadenas y trazas#

CARBayes utiliza MCMC para la estimación. Antes de interpretar coeficientes, se debe verificar convergencia:

plot(carleroux$samples$rho) — la traza de $\rho$ debe verse como “ruido blanco” (sin tendencia ni ciclos), indicando que la cadena ha convergido.
summary(modelo$samples) — el diagnóstico de Geweke (proporcionado automáticamente) compara la primera y última parte de la cadena; un z-score $|z| < 2$ sugiere convergencia.
Si hay problemas de convergencia: aumentar burnin y n.sample, o reducir thin.

4. Interpretación de los coeficientes $\boldsymbol{\beta}$#

Los coeficientes se extraen de $samples$beta. El output de round(beta.results, 5) muestra:

Mean: el valor esperado a posteriori del coeficiente — análogo al coeficiente en regresión frecuentista.
2.5% y 97.5%: intervalo de credibilidad del 95%. Si no incluye el cero, el efecto es estadísticamente significativo.

En modelos de Poisson con offset log(area), los coeficientes se interpretan como log-razones de incidencia (log Incidence Rate Ratios): $e^{\beta_k}$ es el factor multiplicativo en la tasa de deslizamientos por unidad de cambio en la covariable $k$.

5. Interpretación del parámetro $\rho$ en el modelo de Leroux#

El parámetro $\rho$ estimado en carleroux (el modelo sin $\rho$ fijado) indica la fuerza de la dependencia espacial en los datos:

$\hat{\rho}$ cercano a 1 → la dependencia espacial es fuerte; el modelo se acerca al ICAR.
$\hat{\rho}$ cercano a 0 → poca estructura espacial; un modelo i.i.d. sería suficiente.
$\hat{\rho}$ intermedio → el modelo Leroux aporta genuina flexibilidad frente al ICAR fijo.

La traza de plot(carleroux$samples$rho) también permite evaluar si $\rho$ está bien identificado (distribución posterior concentrada) o si los datos no son suficientemente informativos para estimarlo (distribución posterior plana).

ggplot() + geom_sf(data=aoi,aes(fill=carbym),color = "black") +
  annotation_scale(location="br",style = "ticks") +
  annotation_north_arrow(location = "tr",which_north = "true", height = unit(1, "cm"), width = unit(1, "cm"),) +
  scale_fill_gradientn(colors=c("white","red"),name = "Landslides") +
  theme_bw() +
  theme(
    panel.grid.major = element_blank(), 
    panel.grid.minor = element_blank(),
    panel.border = element_rect(color = "black",fill = NA,size = 1),
    axis.ticks.length=unit(-0.1, "cm"),
    axis.text.x = element_text(size = 10, margin = unit(c(t = 1, r = 0, b = 0, l = 0), "mm")),
    axis.text.y = element_text(size = 10, margin = unit(c(t = 0, r = 1, b = 0, l = 0), "mm")),
    legend.text = element_text(size=12),
    legend.title.align = 0,
    legend.position = c(0.34,0.9), 
    legend.key.size = unit(0.5, 'cm'),
    legend.justification = "center",
    legend.direction = "horizontal",
    legend.title = element_text(size=14, vjust = .8, hjust = .5))

Ejemplo con datos de cuencas de la implementación de modelos CAR en R en INLA#

#Spatial matrix
aoi.nb <- poly2nb(aoi) #Queen matrix
aoi.mat <- as(nb2mat(aoi.nb, style = "B"), "Matrix")
aoi.listw = nb2listw(aoi.nb)
colnames(aoi.mat) <- rownames(aoi.mat) 
mat <- as.matrix(aoi.mat[1:dim(aoi.mat)[1], 1:dim(aoi.mat)[1]])

A continuación se representa un modelo de regresión de Poisson estándar donde el conteo de deslizamientos se modela directamente por las variables de precipitación, elevación y pendiente, ajustando por el área como un offset.

m1_inla <- inla(lands_rec ~ 1 + RainfallDaysmean + elev_mean + slope_mean,
           offset=log(area), data = as.data.frame(aoi),family = "poisson",
           control.predictor = list(compute = TRUE),
           control.compute = list(dic = TRUE, waic = TRUE))
summary(m1_inla)
aoi$m1_inla <- m1_inla$summary.fitted[1:nrow(aoi), "mean"]
res_pearson_m1 <- (aoi$lands_rec - m1_inla$summary.fitted.values$mean) / sqrt(m1_inla$summary.fitted.values$mean)
aoi$res_pearson_m1 <- res_pearson_m1
moran_m1 <- moran.mc(x = res_pearson_m1, listw = aoi.listw,nsim = 999, alternative = "greater")

A continuación el modelo m3_icar utiliza INLA para ajustar un modelo de regresión de Poisson que incorpora explícitamente la dependencia espacial entre las unidades geográficas a través de un efecto aleatorio ICAR, utilizando la matriz de vecindad definida.

# ICAR model
m3_icar <- inla(formula = lands_rec ~ 1 + RainfallDaysmean + elev_mean + slope_mean + 
                  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)
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")

Este código implementa el modelo BYM (Besag, York, y Mollié) utilizando la librería INLA.

# BYM model
m4_bym<-inla(formula = lands_rec ~ 1 + RainfallDaysmean + elev_mean + slope_mean +
           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)
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")

A continuación se implementa un modelo CAR Leroux utilizando una especificación más directa de la matriz de precisión dentro de INLA, empleando model = “generic1”. Las primeras líneas construyen la matriz de precisión para el modelo Leroux. apply(aoi.mat, 1, sum) calcula el el número de vecinos de cada unidad espacial sumando las filas de la matriz de adyacencia aoi.mat. Posteriormente se crea una matriz diagonal con el grado de cada nodo en la diagonal principal. A esta matriz diagonal se le resta la matriz de adyacencia aoi.mat. El resultado, m5_leroux, es similar a la matriz Laplaciana no normalizada del grafo de vecindad. Cmatrix <- Diagonal(nrow(aoi), 1) - m5_leroux:

Diagonal(nrow(aoi), 1) crea una matriz identidad de dimensión igual al número de unidades espaciales. A esta matriz identidad se le resta m5_leroux. La matriz resultante, Cmatrix, representa la matriz de precisión (inversa de la matriz de covarianza) del modelo CAR Leroux sin el parámetro de dependencia espacial (ρ). Finalmente se calcula los autovalores de la matriz Cmatrix y luego encuentra el máximo de estos autovalores. INLA internamente estimará el parámetro de precisión (inverso de la varianza) asociado con este efecto.

#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(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)
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")

cuenca_colors <- c("Atrato" = "green", "Cauca" = "red", "Magdalena" = "blue")

p1 <- ggplot() + 
  geom_sf(data = aoi, aes(fill = res_pearson_m1, color = cuenca), size = 0.7) +
  scale_color_manual(values = cuenca_colors, guide = "none") +
  scale_fill_gradient2(
    low = "yellow", mid = "white", high = "purple", midpoint = 0,
    name = "Res. Pearson"
  ) +
  annotation_scale(
    location = "br", style = "ticks",
    text_cex = 1.0  # Tamaño de fuente en escala
  ) +
  annotation_north_arrow(
    location = "tr", which_north = "true",
    height = unit(1, "cm"), width = unit(1, "cm"),
    pad_x = unit(0.1, "cm"), pad_y = unit(0.1, "cm")
  ) +
  theme_bw() +
  theme(
    legend.position = c(0.3, 0.9),
    legend.direction = "horizontal",
    legend.title = element_text(size = 12),
    legend.text = element_text(size = 11),
    legend.key.height = unit(0.4, "cm"),
    legend.key.width = unit(0.5, "cm"),
    axis.text = element_text(size = 11),          # Coordenadas del mapa
    axis.title = element_text(size = 13),         # Si decides añadir títulos
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    plot.title = element_blank()
  )

p2 <- ggplot(aoi, aes(x = cuenca, y = res_pearson_m1, fill = cuenca)) +
  geom_boxplot(notch = TRUE, outlier.shape = 1) +
  scale_fill_manual(values = cuenca_colors) +
  theme_bw() +
  labs(
    y = "Res. Pearson", x = ""  # Etiqueta eje Y
  ) +
  theme(
    legend.position = "none",
    plot.title = element_blank(),
    axis.text.x = element_text(size = 12),         # Etiquetas cuenca
    axis.text.y = element_text(size = 10),         # Ticks eje Y
    axis.title.y = element_text(size = 12),  # Título eje Y
    aspect.ratio = 1.95
  )

combined_plot <- plot_grid(p1, p2, ncol = 2, rel_widths = c(1.2, 0.8))

final_plot <- ggdraw(combined_plot) +
  draw_plot_label(
    label = c("A", "B"),
    x = c(0.1, 0.7), y = c(0.92, 0.92),
    hjust = 0, vjust = 1,
    fontface = "bold", size = 14)

print(final_plot)

Interpretación del output de INLA y comparación de modelos#

Lectura del `summary()` en INLA#

El objeto devuelto por inla() contiene toda la inferencia bayesiana posterior. Los elementos clave son:

Efectos fijos (summary(modelo)$fixed):

Campo	Significado
`mean`	Media posterior del coeficiente (interpretación análoga a MLE)
`sd`	Desviación estándar posterior
`0.025quant`, `0.975quant`	Intervalo de credibilidad del 95%
`mode`	Moda de la distribución posterior

Si el intervalo [0.025quant, 0.975quant] no contiene el cero, el efecto es estadísticamente relevante.

Hiperparámetros (summary(modelo)$hyperpar):

Para modelos con efectos aleatorios, INLA estima la precisión $\tau = 1/\sigma^2$. Un valor alto de $\tau$ indica efectos aleatorios pequeños (poca varianza no explicada).
Para el modelo Leroux (generic1), el hiperparámetro Beta for ID corresponde a $\rho$ — el nivel de dependencia espacial estimado.

Selección de modelos con DIC y WAIC#

INLA calcula automáticamente DIC y WAIC cuando se especifica control.compute = list(dic=TRUE, waic=TRUE). La lógica de comparación:

Criterio	Regla	Umbral práctico
DIC	Menor = mejor	$\Delta$DIC > 10 es sustancial
WAIC	Menor = mejor	$\Delta$WAIC > 5 sugiere diferencia real

Se espera la siguiente jerarquía en términos de ajuste: Modelo base (sin CAR) < ICAR < BYM ≈ Leroux. Sin embargo, el modelo más complejo no siempre gana — un modelo más simple puede ser preferible si la diferencia de DIC/WAIC es pequeña.

Diagnóstico de residuos: Moran I secuencial#

La secuencia de tests de Moran (moran_m1, moran_m3, moran_m4, moran_m5) permite evaluar si la estructura espacial explícita reduce la dependencia residual:

Modelo 1 (sin CAR) → Moran I significativo (p < 0.05) ✓ Hay dependencia espacial
Modelo 3 (ICAR)   → Moran I no significativo           ✓ Dependencia capturada
Modelo 4 (BYM)    → Moran I no significativo           ✓ Dependencia capturada
Modelo 5 (Leroux) → Moran I no significativo           ✓ Dependencia capturada

Si los modelos CAR todavía muestran Moran I significativo, considerar: (a) redefinir la matriz de vecindad, (b) incluir covariables adicionales, o (c) usar modelos de mayor complejidad (BYM jerárquico, ver siguiente capítulo).

Visualización de los efectos espaciales#

Los valores ajustados m3_icar$summary.fitted[, "mean"] representan la tasa esperada de deslizamientos por subcuenca, suavizada por la estructura espacial. Al mapear estos valores junto con los residuos de Pearson:

Valores ajustados: muestran el patrón espacial capturado por el modelo — zonas de alto/bajo riesgo.
Residuos de Pearson positivos (observado > esperado): el modelo subestima el riesgo en esa subcuenca.
Residuos de Pearson negativos (observado < esperado): el modelo sobreestima el riesgo.

La comparación de los mapas de residuos entre el modelo base y el modelo ICAR/BYM/Leroux debe mostrar una reducción en la estructura espacial de los residuos (menos clústers visibles), confirmando que la dependencia espacial fue absorbida por el efecto aleatorio latente.

Actividades#

Comparación ICAR, BYM y Leroux: Ajusta los tres modelos CAR (ICAR, BYM y Leroux) para el dataset de cuencas usando CARBayes. Construye una tabla comparativa con: DIC, WAIC, Moran I de residuos (p-valor), y $\hat{ ho}$ (para Leroux). ¿Qué modelo recomendarías y por qué?
Convergencia MCMC: Para el modelo Leroux, grafica las trazas MCMC de los tres coeficientes $eta$ y del parámetro $ ho$. Evalúa la convergencia usando el diagnóstico de Geweke. Si alguna cadena no converge, ¿qué ajustes en los parámetros de muestreo (burnin, n.sample, thin) harías?
Interpretación geomorfológica de coeficientes: Para el modelo BYM ajustado al dataset de cuencas, calcula los Incidence Rate Ratios ($e^{\hat{eta}_k}$) para cada covariable. Interpreta cada uno en términos geomorfológicos: ¿qué significa un IRR de 1.5 para la pendiente media?
Efectos aleatorios espaciales: Del modelo ICAR ajustado con INLA, extrae los efectos aleatorios $\hat{\phi}_i$ (m3_icar$summary.random$ID$mean). Mapea estos efectos. Las subcuencas con $\hat{\phi}_i > 0$ tienen una tasa de deslizamientos más alta de lo esperado por las covariables. ¿Se agrupan geográficamente estas subcuencas? ¿Qué variables no medidas podrían explicarlas (litología, cobertura vegetal, etc.)?

Recursos adicionales#

Documentación oficial#

Lecturas recomendadas#

Besag, J., York, J., & Mollié, A. (1991). Bayesian image restoration, with two applications in spatial statistics. Annals of the Institute of Statistical Mathematics, 43(1), 1–20.
Leroux, B. G., Lei, X., & Breslow, N. (2000). Estimation of disease rates in small areas: A new mixed model for spatial dependence. Statistical Models in Epidemiology, the Environment, and Clinical Trials, 179–191.
Lee, D. (2013). CARBayes: An R package for Bayesian spatial modeling with conditional autoregressive priors. Journal of Statistical Software, 55(13), 1–24.
Rue, H. & Held, L. (2005). Gaussian Markov Random Fields: Theory and Applications. CRC Press.
Gómez-Rubio, V. (2020). Bayesian Inference with INLA. CRC Press. Disponible online

Modelos de regresión para dependencia espacial tipo CAR

Contents

Modelos de regresión para dependencia espacial tipo CAR#

Objetivos de aprendizaje#

Diferencia conceptual entre SAR y CAR#

Especificación: global (SAR) vs. condicional (CAR)#

SAR: sistema simultáneo#

CAR: propiedad de Markov espacial#

¿Cuándo preferir uno u otro?#

Modelos con la librería CARBayes en R#

Modelos con la librería INLA en R#

Modelo gaussiano clásico#

Modelo Gaussiano con efectos aleatorios NO estructurados#

Modelo Gaussiano con efectos aleatorios estructurados#

Modelo ICAR (Intrinsic CAR)#

Modelo BYM (Besag, York y Mollié)#

Modelo de Leroux et al. (CAR propio)#

Ejemplo con datos de cuencas de la implementación de modelos CAR en R en CARBayes#

Interpretación y comparación de los modelos CARBayes#

1. Diagnóstico de autocorrelación espacial residual#

2. Ajuste del modelo: `modelfit`#

3. Convergencia MCMC: cadenas y trazas#

4. Interpretación de los coeficientes \(\boldsymbol{\beta}\)#

5. Interpretación del parámetro \(\rho\) en el modelo de Leroux#

Ejemplo con datos de cuencas de la implementación de modelos CAR en R en INLA#

Interpretación del output de INLA y comparación de modelos#

Lectura del `summary()` en INLA#

Selección de modelos con DIC y WAIC#

Diagnóstico de residuos: Moran I secuencial#

Visualización de los efectos espaciales#

Actividades#

Recursos adicionales#

Documentación oficial#

Lecturas recomendadas#

Recursos interactivos#

Criterio	Regla	Umbral práctico
DIC	Menor = mejor	\(\Delta\)DIC > 10 es sustancial
WAIC	Menor = mejor	\(\Delta\)WAIC > 5 sugiere diferencia real

Modelos de regresión para dependencia espacial tipo CAR

Contents

Modelos de regresión para dependencia espacial tipo CAR#

Objetivos de aprendizaje#

Diferencia conceptual entre SAR y CAR#

Especificación: global (SAR) vs. condicional (CAR)#

SAR: sistema simultáneo#

CAR: propiedad de Markov espacial#

¿Cuándo preferir uno u otro?#

Modelos con la librería CARBayes en R#

Modelos con la librería INLA en R#

Modelo gaussiano clásico#

Modelo Gaussiano con efectos aleatorios NO estructurados#

Modelo Gaussiano con efectos aleatorios estructurados#

Modelo ICAR (Intrinsic CAR)#

Modelo BYM (Besag, York y Mollié)#

Modelo de Leroux et al. (CAR propio)#

Ejemplo con datos de cuencas de la implementación de modelos CAR en R en CARBayes#

Interpretación y comparación de los modelos CARBayes#

1. Diagnóstico de autocorrelación espacial residual#

2. Ajuste del modelo: modelfit#

3. Convergencia MCMC: cadenas y trazas#

4. Interpretación de los coeficientes \(\boldsymbol{\beta}\)#

5. Interpretación del parámetro \(\rho\) en el modelo de Leroux#

Ejemplo con datos de cuencas de la implementación de modelos CAR en R en INLA#

Interpretación del output de INLA y comparación de modelos#

Lectura del summary() en INLA#

Selección de modelos con DIC y WAIC#

Diagnóstico de residuos: Moran I secuencial#

Visualización de los efectos espaciales#

Actividades#

Recursos adicionales#

Documentación oficial#

Lecturas recomendadas#

Recursos interactivos#

2. Ajuste del modelo: `modelfit`#

Lectura del `summary()` en INLA#