Procesos Gaussianos con R

Créditos: El contenido de este cuaderno ha sido tomado de varias fuentes, pero especialmente de Lionel Hertzog, Spatial Statistics for Data Science: Theory and Practice with R. El compilador se disculpa por cualquier omisión involuntaria y estaría encantado de agregar un reconocimiento.

Procesos Gaussianos con R#

La geoestadística basada en modelos puede utilizarse para analizar datos espaciales relacionados con un fenómeno subyacente continuo en el espacio que se ha recolectado en un conjunto finito de ubicaciones. La geoestadística basada en modelos emplea modelos estadísticos para capturar la estructura de correlación espacial en los datos, lo que permite realizar inferencias estadísticas rigurosas y facilita la producción de predicciones espaciales junto con medidas de incertidumbre del fenómeno de interés (Diggle, Tawn y Moyeed 1998).

Asumiendo datos gaussianos observados en un conjunto de \(n\) ubicaciones, \({Y₁, …, Yₙ}\), podemos considerar el siguiente modelo para obtener predicciones en ubicaciones no muestreadas:

\(Yᵢ | S(sᵢ) ∼ N(μ + S(sᵢ), τ²), i = 1, …, n.\)

Aquí, \(μ\) es un efecto medio constante, y \(S(⋅)\) es un campo gaussiano espacial de media cero. Este modelo puede extenderse a situaciones en las que la variación estocástica en los datos no sea gaussiana, así como para incluir covariables y otros efectos aleatorios para tener en cuenta otros tipos de variabilidad.

La inferencia en la geoestadística basada en modelos puede realizarse utilizando los enfoques INLA y la ecuación diferencial parcial estocástica (SPDE), los cuales proporcionan una alternativa computacionalmente eficiente a los métodos MCMC (Lindgren y Rue 2015). En resumen, esto implica resolver una SPDE en una malla discreta de puntos e interpolar para obtener una solución continua a lo largo del dominio espacial (Krainski et al. 2019), que se calcula utilizando INLA (Rue, Martino y Chopin 2009).

library(geoR)
library(ggplot2)
library(tidyverse)
library(INLA)

# the ca20 dataset
# load the example dataset,
# calcium content in soil samples in Brazil
data(ca20)
# put this in a data frame
dat <- data.frame(x = ca20$coords[,1],y = ca20$coords[,2],calcium = ca20$data, 
                  elevation = ca20$covariate[,1], 
                  region = factor(ca20$covariate[,2]))

# meshes in 2D space can be created as follow:
mesh <- inla.mesh.2d(loc = dat[,c("x", "y")], max.edge = c(50, 5000))

La función inla.mesh.2d del paquete INLA se utiliza para crear una malla para el modelado espacial. Esta malla aproxima el efecto espacial en sus datos y juega un papel crucial en la captura de la dependencia espacial.

Aquí un desglose de los puntos clave:

  • Función: inla.mesh.2d

  • Propósito: Crea una malla bidimensional para el análisis espacial en modelos INLA.

  • Argumentos:

    • loc: Esto especifica la ubicación de sus muestras. Suele ser un marco de datos con columnas para latitud y longitud (u otras coordenadas espaciales).

    • max.edge: Controla la precisión de la malla. Define la distancia máxima permitida entre dos nodos conectados en la malla. Un valor más bajo (por ejemplo, 50 metros) crea una malla más fina con más detalles, pero requiere más recursos computacionales. Por el contrario, un valor más alto (por ejemplo, 5000 metros) crea una malla más gruesa que es computacionalmente más rápida pero puede pasar por alto variaciones espaciales más finas.

Compromiso entre precisión y cómputo:

Existe una compensación entre la precisión de la malla y el tiempo de cómputo. Una malla más precisa (max.edge más pequeño) conduce a una predicción más suave que captura patrones espaciales sutiles, pero tarda más en calcularse. Una malla más gruesa (max.edge más grande) es computacionalmente más rápida pero puede pasar por alto detalles espaciales importantes.

Exploración interactiva con meshbuilder:

El paquete INLA también proporciona la función meshbuilder para la exploración interactiva de la creación de mallas. Esto le permite visualizar diferentes opciones de malla basadas en su configuración max.edge y elegir la que mejor equilibre la precisión y la viabilidad computacional para su análisis específico.

En resumen, inla.mesh.2d es una herramienta esencial para crear mallas en el modelado espacial de INLA. Al considerar cuidadosamente la compensación entre precisión y cómputo, y potencialmente usar meshbuilder para la exploración, puede crear una malla adecuada para capturar los efectos espaciales en sus datos.

Amat <- inla.spde.make.A(mesh, loc = as.matrix(dat[,c("x", "y")]))

Efecto espacial en INLA con ecuaciones diferenciales parciales estocásticas (EDPE)#

INLA estima el efecto espacial utilizando una herramienta matemática compleja llamada ecuación diferencial parcial estocástica (EDPE). La idea básica es que podemos estimar un efecto espacial continuo utilizando un conjunto de puntos discretos (los nodos definidos en la malla) y funciones base, similar a las splines de regresión. Esto facilita mucho la estimación de los campos espaciales.

El enfoque de la ecuación diferencial parcial estocástica (SPDE) implementado en el paquete R-INLA proporciona una forma flexible y computacionalmente eficiente de modelar datos geoestadísticos y realizar predicciones en ubicaciones no muestreadas (Lindgren y Rue 2015). Suponemos que, subyacente a los datos observados, hay una variable continua en el espacio que puede modelarse utilizando un campo aleatorio gaussiano (GRF) con una función de covarianza de Matérn, que se define como:

\(Cov(x(sᵢ), x(sⱼ)) = (σ² / 2^(ν - 1)Γ(ν))(κ ||sᵢ - sⱼ||)^ν Kν(κ ||sᵢ - sⱼ||).\)

Aquí, \(σ²\) denota la varianza marginal del campo espacial. \(Kν(⋅)\) se refiere a la función de Bessel modificada de segunda especie y orden \(ν > 0\). El valor entero de \(ν\) determina la suavidad del campo y típicamente se fija, ya que es difícil de estimar en aplicaciones. \(κ > 0\) está relacionado con el rango \(ρ\), que representa la distancia a la cual la correlación entre dos puntos se vuelve aproximadamente 0. Específicamente, \(ρ = √(8ν) / κ\), y a esta distancia la correlación espacial es cercana a 0.1 (Cameletti et al. 2013).

Como se muestra en Whittle (1963), un GRF con una matriz de covarianza de Matérn puede representarse como la solución de la siguiente SPDE en dominio continuo:

\((κ² − Δ)^(α/2)(τx(s)) = W(s).\)

Aquí, \(x(s)\) representa un GRF, y \(W(s)\) es un proceso de ruido blanco espacial gaussiano. El parámetro \(α\) controla la suavidad exhibida por el GRF, \(τ\) controla su varianza, y \(κ > 0\) es un parámetro de escala. El laplaciano \(Δ\) se define como \(∑(dᵢ = 1 ∂² / ∂x²ᵢ)\), donde \(d\) es la dimensión del dominio espacial.

Los parámetros de la función de covarianza de Matérn y la SPDE están relacionados de la siguiente manera. El parámetro de suavidad \(ν\) de la función de covarianza de Matérn se expresa como \(ν = α − d/2\), y la varianza marginal \(σ²\) está relacionada con la SPDE mediante:

\(σ² = (Γ(ν) / Γ(α))(4π)^(d/2)κ^(−2ν)τ².\)

En el caso donde \(d = 2\) y \(ν = 1/2\), lo que corresponde a la función de covarianza exponencial, el parámetro \(α = ν + d/2 = 1/2 + 1 = 3/2\). En el paquete R-INLA, el valor predeterminado es \(α = 2\), aunque también están disponibles opciones dentro del rango \(0 ≤ α < 2\).

El método de elementos finitos puede usarse para encontrar una solución aproximada a la SPDE. Este método implica dividir el dominio espacial en un conjunto de triángulos no superpuestos, creando una malla triangulada con \(n\) nodos y \(n\) funciones base. Las funciones base, denotadas como \(ψₖ(⋅\))$, son funciones lineales por partes en cada triángulo. Toman el valor de 1 en el vértice \(k\), y 0 en todos los demás vértices.

Luego, el campo gaussiano indexado de manera continua \(x\) se representa como un campo aleatorio de Markov gaussiano indexado de manera discreta (GMRF) mediante una suma de funciones base definidas en la malla triangulada:

\(x(s) = ∑(nₖ = 1) ψₖ(s)xₖ,\)

donde \(n\) es el número de vértices de la triangulación, \(ψₖ(⋅)\) representa las funciones base lineales por partes, y \({xₖ}\) denotan pesos distribuidos de forma gaussiana con media cero.

La distribución conjunta del vector de pesos se asigna a una distribución gaussiana representada como \(x = (x₁, …, xₙ) ∼ N(0, Q^−1(τ, κ))\)

Esta distribución aproxima la solución \(x(s)\) de la SPDE en los nodos de la malla. Las funciones base transforman la aproximación \(x(s)\) desde los nodos de la malla a las otras ubicaciones espaciales de interés.

Ahora bien, la parte complicada de configurar la EDPE es definir los previos para:

  • Alcance del efecto espacial (parámetro κ): Este parámetro representa la distancia a partir de la cual dos puntos se pueden considerar espacialmente independientes. En otras palabras, ¿a qué distancia deben estar dos ubicaciones para que su efecto espacial sea prácticamente nulo? Un κ alto indica un efecto espacial de largo alcance (las ubicaciones distantes se influencian entre sí), mientras que un κ bajo indica un efecto de corto alcance (la influencia espacial se limita a las ubicaciones cercanas).

  • Variación del efecto espacial (parámetro δ): Este parámetro representa la variabilidad del campo espacial de un punto a otro. Un δ alto indica un campo espacial con mucha variación (los valores del efecto espacial cambian drásticamente a lo largo del espacio), mientras que un δ bajo indica un campo espacial con poca variación (los valores del efecto espacial son relativamente similares en todo el espacio).

Encontrar los valores óptimos para κ y δ es crucial para capturar adecuadamente el efecto espacial en su modelo. INLA permite especificar valores iniciales para estos parámetros, y luego estima sus valores finales durante el proceso de ajuste del modelo.

spde <- inla.spde2.pcmatern(mesh, 
                            prior.range = c(500, 0.5),
                            prior.sigma = c(2, 0.05))

Definiendo Priores para el Efecto Espacial en INLA: Alcance y Variación#

Establecer valores iniciales (priores) para el alcance y la variación del efecto espacial en INLA (representados por los parámetros κ y δ respectivamente) puede ser una tarea delicada, ya que estos valores tienen un gran impacto en el modelo ajustado. Afortunadamente, como el ajuste de modelos en INLA es relativamente rápido, es fácil realizar un análisis de sensibilidad para comprender qué configuraciones funcionan mejor.

Alcance del Efecto Espacial (parámetro κ)#

El prior para el alcance del efecto espacial corresponde con la siguiente fórmula:

P(κ < κ₀) = p₀

  • κ: Representa el alcance del efecto espacial (distancia a partir de la cual se considera que dos puntos son independientes).

  • κ₀: Valor inicial específico del alcance.

  • p₀: Probabilidad asociada al valor inicial κ₀.

Por ejemplo, si definimos:

P(κ < 500) = 0.5

Significa que estamos estableciendo una probabilidad del 50% de que el alcance real del efecto espacial sea menor a 500 metros. En otras palabras, hay una posibilidad del 50% de que dos ubicaciones separadas por más de 500 metros tengan efectos espaciales independientes. A mayor valor de κ, mayor alcance del efecto espacial (las ubicaciones distantes se influencian entre sí).

Variación del Efecto Espacial (parámetro δ)#

De manera similar, el prior para la variación del efecto espacial se establece mediante la siguiente fórmula:

P(δ > δ₀) = p₀

  • δ: Representa la variación del efecto espacial (cuánto cambia el efecto de un punto a otro).

  • δ₀: Valor inicial específico de la variación.

  • p₀: Probabilidad asociada al valor inicial δ₀.

Por ejemplo:

P(δ > 2) = 0.05

Indica que asignamos una probabilidad del 5% a la posibilidad de que la variación del efecto espacial sea mayor a 2. En otras palabras, hay una probabilidad del 95% de que la variación sea menor o igual a 2, lo que sugiere un campo espacial con una variación relativamente baja (los valores del efecto cambian poco en el espacio).

Recomendaciones:

  • Experimentar con diferentes valores de κ₀, δ₀ y p₀ para observar cómo afectan el modelo ajustado.

  • Consultar fuentes adicionales para encontrar configuraciones apropiadas según el tipo de análisis que se esté realizando.

  • Se recomienda establecer valores iniciales relativamente fuertes para δ (p₀ cercano a 0 y alejado de 0.5). Priores demasiado vagos (p₀ cercano a 0.5) pueden causar problemas en el modelo, especialmente en modelos complejos.

Tenga en cuenta que INLA permite especificar estos valores iniciales y luego estima los valores finales durante el proceso de ajuste del modelo. El análisis de sensibilidad y la referencia a otras fuentes te ayudarán a encontrar una buena configuración inicial para tus modelos espaciales en INLA.

# create the data stack
dat_stack <- inla.stack(data = list(calcium = dat$calcium), # the response variable
                        A = list(Amat, 1, 1, 1), # the projection matrix
                        effects = list(i = 1:spde$n.spde, # the spatial effect
                                       Intercept = rep(1, nrow(dat)), 
                                       elevation = dat$elevation,
                                       region = factor(dat$region)))

La clave está en la proyección de los efectos#

El punto clave aquí es el argumento A donde especificamos la proyección de los diferentes efectos. El efecto espacial se denomina i (aunque podemos nombrarlo como queramos) y está indexado por el número de nodos de la malla. Recordemos que cuanto más fina sea la malla, más precisa será la estimación del efecto espacial. Este efecto espacial i está vinculado a los datos a través de la matriz de proyección A_mat.

Los demás efectos se vinculan directamente a los datos, por lo que no necesitan matrices de proyección.

Aquí desglosamos un poco más la terminología:

  • Efecto espacial (i): Representa la variación espacial que se estima a través de la EDPE. Se modela utilizando los nodos de la malla y las funciones base.

  • Matriz de proyección (A_mat): Esta matriz vincula el efecto espacial estimado en cada nodo de la malla con las observaciones puntuales. Básicamente, indica cómo contribuye el efecto espacial en cada nodo al valor de la variable respuesta en cada ubicación observada.

  • Efectos directos: Estos efectos son covariables que se incluyen en el modelo además del efecto espacial. A diferencia del efecto espacial, se relacionan directamente con las observaciones puntuales sin necesidad de una matriz de proyección.

Predicción en Modelos Espaciales INLA#

En INLA, por lo general, es más sencillo obtener predicciones del modelo pasando directamente los nuevos datos que se quieren utilizar para predecir al ajuste del modelo. En otras palabras, necesitamos definir estos nuevos datos antes de ajustar el modelo.

Aquí vamos a ver cómo predecir el efecto de la elevación y la región sobre la variable respuesta, teniendo en cuenta el efecto espacial:

  1. Datos para la predicción:

    • Comenzaremos por definir un nuevo conjunto de datos que solo contenga las variables de elevación y región para las ubicaciones donde queremos predecir la variable respuesta.

    • Este nuevo conjunto de datos debe tener el mismo formato que las variables de elevación y región utilizadas en el modelo original.

    • Es importante asegurarse de que las ubicaciones para las que se desea predecir estén dentro del rango del área de estudio cubierta por el modelo original.

  2. Ajuste del modelo:

    • Durante el ajuste del modelo en INLA, se incluirá este nuevo conjunto de datos junto con los datos originales utilizados para entrenar el modelo.

    • INLA utilizará los valores de elevación y región en los nuevos datos para predecir los valores de la variable respuesta en esas ubicaciones específicas.

Al incluir los nuevos datos en el ajuste del modelo, INLA tendrá en cuenta tanto el efecto directo de la elevación y la región como el efecto espacial estimado en el modelo original. Esto permite obtener predicciones más precisas que consideren la dependencia espacial y las relaciones con las covariables de interés.

# a newdata to get the predictions
modmat <- expand.grid(elevation = seq(min(dat$elevation), 
                                      max(dat$elevation), 
                                      length.out = 10),
                      region = unique(dat$region))

# the stack for these predictions
pred_stack_fixef <- inla.stack(data = list(calcium = NA),
                               A = list(1, 1, 1),
                               effects = list(Intercept = rep(1, nrow(modmat)),
                                              elevation = modmat$elevation,
                                              region = factor(modmat$region)),
                               tag = "prd_fixef")

Puntos clave de la predicción en INLA#

Aquí hay un desglose de los puntos clave de la predicción en el modelo INLA que se acaba de describir:

  • Datos para la predicción:

    • Se define un nuevo conjunto de datos que solo incluye las variables de elevación y región para las ubicaciones donde se desean las predicciones.

    • El modelo utilizará estos valores para predecir la variable respuesta en esas ubicaciones específicas.

    • Un aspecto crucial es que se establece calcium=NA en este conjunto de datos. Esto indica a INLA que estime los valores de calcio (calcium) basándose en los efectos (incluido el espacial) y los parámetros del modelo.

    • La etiqueta prd_fixef en el stack permite posteriormente extraer fácilmente los valores predichos.

  • Predicción espacial:

    • Debido al modelado espacial, las predicciones también se pueden realizar a lo largo del espacio.

    • Una opción sería predecir solo en base al campo espacial, pero lo más interesante es tener en cuenta también las covariables (elevación y región) para obtener predicciones espaciales más precisas.

    • Obtener el stack de predicción para este escenario es un poco más complejo. Se necesitan valores de elevación y región no solo en las ubicaciones observadas, sino en todo el espacio de estudio.

  • Pasos adicionales:

    • Se requieren algunos pasos previos a INLA que pueden parecer complicados.

    • El objetivo final es crear rásteres con información de elevación y región a partir de los datos disponibles. Estos rásteres se utilizarán luego para predecir la variable respuesta en todo el espacio de estudio.

En resumen, la predicción en INLA para modelos espaciales implica definir datos para predicciones puntuales y, si se desea, generar predicciones espaciales teniendo en cuenta los efectos espaciales y las covariables. Los pasos adicionales previos a INLA pueden incluir la creación de rásteres a partir de los datos para cubrir todo el espacio de estudio.

library(raster)
library(fields) # for Tps

## first we define an empty raster to hold the coordinates of the predictions
r <- raster(xmn = min(dat$x), xmx = max(dat$x),
            ymn = min(dat$y), ymx = max(dat$y),
            resolution = 25)

## the we use thin-plate spline to derive elevation across the data
elev_m <- Tps(dat[,c("x","y")], dat$elevation)

## put this into a raster
elev <- interpolate(r, elev_m)

## for the region info we create a SpatialPolygons 
## based on the coordinates given in the ca20 object
pp <- SpatialPolygons(list(Polygons(list(Polygon(ca20[[5]])), ID = "reg1"),
                           Polygons(list(Polygon(ca20[[6]])), ID = "reg2"),
                           Polygons(list(Polygon(ca20[[7]])), ID = "reg3")))
# turn the SpatialPolygon into a raster object
region <- rasterize(pp, r)

# the new data frame with coordinates from the raster
# plus elevation and region information
newdat <- as.data.frame(xyFromCell(r, cell = 1:ncell(r)))
newdat$elevation <- values(elev)
newdat$region <- factor(values(region))
# remove NAs
newdat <- na.omit(newdat)

# create a new projection matrix for the points
Apred <- inla.spde.make.A(mesh,
                          loc = as.matrix(newdat[,c("x", "y")]))

# put this in a new stack
pred_stack_alleff <- inla.stack(data = list(calcium = NA),
                               A = list(Apred, 1, 1, 1),
                               effects = list(i = 1:spde$n.spde,
                                              Intercept = rep(1, nrow(newdat)),
                                              elevation = newdat$elevation,
                                              region = factor(newdat$region)),
                               tag = "prd_alleff")

Obteniendo predicciones espaciales con covariables en INLA#

Aquí detallamos el proceso para obtener predicciones espaciales teniendo en cuenta las covariables (elevación y región) en un modelo INLA:

  1. Información espacial de rásteres:

    • Se parte de un ráster que representa la región de interés. Este ráster debe contener información sobre las variables de elevación y región.

  2. Datos para la predicción (newdat):

    • Se crea un nuevo objeto de datos (newdat) que contiene los valores de elevación y región extraídos del ráster para las ubicaciones donde se desean las predicciones.

    • Es importante asegurarse de que la extensión del ráster cubra el área de interés para las predicciones.

    • Al igual que en el caso anterior, se establece calcium=NA en este conjunto de datos para indicar que INLA debe estimar los valores de la variable respuesta (calcio) a partir de los efectos y parámetros del modelo.

  3. Matriz de proyección:

    • Se define una nueva matriz de proyección que refleje cómo contribuyen el efecto espacial estimado (campo espacial), la elevación y la región a las predicciones en cada ubicación.

    • Esta matriz de proyección es más compleja que la utilizada para predicciones puntuales porque tiene que tener en cuenta todos los efectos que se utilizan para la predicción.

  4. Nueva pila (stack) de predicción:

    • Se crea una nueva pila de predicción en INLA que incluya:

      • El objeto newdat que contiene los valores de elevación y región para las predicciones.

      • La matriz de proyección recién definida.

      • Una etiqueta única para identificar esta pila de predicción (por ejemplo, “prediccion_espacial”).

# put all the stacks together
all_stack <- inla.stack(dat_stack, pred_stack_fixef,
                      pred_stack_alleff)

# fit the model
m_inla <- inla(calcium ~ -1 + Intercept + elevation + region + f(i, model = spde),
            data = inla.stack.data(all_stack),
            control.predictor = list(A = inla.stack.A(all_stack), compute = TRUE),
            quantiles = NULL)

Ejecutando el modelo INLA para predicción espacial#

Se estima que la ejecución del modelo dure alrededor de 30 segundos. Aquí desglosamos los pasos principales del código que acabaste de mencionar:

1. Ajuste del modelo:

  • fórmula del modelo: El primer argumento define la fórmula del modelo. Se utiliza -1 para eliminar la intercepción interna y ajustarla por separado.

  • efecto espacial aleatorio (f()): La función f() se utiliza para especificar un efecto aleatorio (i) que sigue el modelo de EDPE definido anteriormente.

  • datos y matriz de proyección: Se pasan los datos de predicción (newdat) y la matriz de proyección recién creada.

  • estimación de valores (compute=TRUE): Se configura compute=TRUE para indicar a INLA que estime los valores de la variable respuesta (calcio) que se proporcionaron como NA en newdat.

2. Resumen del modelo:

Una vez ajustado el modelo, se puede obtener su resumen utilizando las funciones de INLA para ver información como:

  • Efectos fijos (coeficientes estimados para elevación y región)

  • Parámetros de correlación espacial (nu y rho)

  • Valores p y otros indicadores de significación estadística

summary(m_inla)

Examinando las predicciones del modelo INLA#

Una vez ajustado el modelo INLA para las predicciones espaciales, podemos analizar los resultados de diferentes maneras:

1. Predicciones basadas solo en efectos fijos:

Ahora que lo mencionas, podemos comenzar por observar las predicciones que solo tienen en cuenta los efectos fijos (elevación y región). Esto implica promediar las variaciones espaciales para ver el impacto general de estas covariables en la variable respuesta:

  • INLA proporciona funciones para extraer las predicciones basadas únicamente en los efectos fijos.

  • Al analizar estas predicciones, podemos evaluar cómo cambian los valores predichos en función de la elevación y la región, promediando el efecto espacial.

2. Predicciones espaciales completas:

Además de los efectos fijos, las predicciones espaciales completas también incorporan el efecto espacial estimado en el modelo.

  • INLA permite extraer las predicciones espaciales totales, que reflejan la variación de la variable respuesta en todo el espacio de estudio.

  • Estas predicciones considerarán tanto la influencia de las covariables (elevación y región) como la dependencia espacial.

3. Visualización de las predicciones:

Una vez extraídas las predicciones, ya sean de efectos fijos o espaciales completas, podemos visualizarlas para comprender mejor los patrones espaciales y la relación con las covariables:

  • Podemos utilizar herramientas de representación geográfica (por ejemplo, mapas) para mostrar las predicciones espaciales en el contexto del área de estudio.

  • La visualización de las predicciones en relación con los valores observados de la variable respuesta permite evaluar la precisión del modelo.

## first we create an index to easily find these 
## prediction within the fitted model
id_fixef <- inla.stack.index(all_stack, "prd_fixef")$data

## add to modmat the prediction and their sd
modmat$calcium <- m_inla$summary.fitted.values[id_fixef, "mean"]
modmat$sd <- m_inla$summary.fitted.values[id_fixef, "sd"]

## a plot with the original data
ggplot(dat, aes(x = elevation, y = calcium)) +
  geom_ribbon(data = modmat, aes(ymin = calcium - 2 * sd,
                                 ymax = calcium + 2 * sd,
                                 fill = region),
              alpha = 0.2) +
  geom_line(data = modmat, aes(color = region)) +
  geom_point(aes(color = region))

Visualizando las predicciones de efectos fijos en INLA#

Ahora que hemos visto la teoría, veamos cómo visualizar las predicciones basadas únicamente en los efectos fijos (elevación y región) en tu modelo INLA:

1. Extraer predicciones de efectos fijos:

Como mencionaste anteriormente, utilizaste una etiqueta específica (por ejemplo, “prediccion_efectos_fijos”) al definir la pila de predicción. Esto nos permite extraer fácilmente las predicciones relevantes del objeto del modelo.

En INLA, el resumen del modelo (summary) contiene un marco de datos llamado summary.fitted.values. Este marco de datos almacena información sobre las predicciones, y podemos filtrarlo utilizando la etiqueta definida anteriormente.

2. Obtener media y desviación estándar:

Una vez filtradas las predicciones de efectos fijos, podemos calcular la media y la desviación estándar. Estas métricas resumen el impacto general de las covariables (elevación y región) en la variable respuesta, promediando la variación espacial.

3. Representación gráfica:

Finalmente, podemos representar gráficamente la media y la desviación estándar junto con los datos originales. Esto te permitirá:

  • Visualizar cómo cambian los valores predichos en función de la elevación y la región.

  • Comparar las predicciones de efectos fijos con los datos observados para evaluar la capacidad del modelo para capturar las tendencias generales.

4. El mapa genial:

¡Y ahora viene la parte emocionante: el mapa! Una vez que tengas las predicciones espaciales completas (que incluyen el efecto espacial), podrás visualizarlas en un mapa. Esto te permitirá ver la variación de la variable respuesta en todo el espacio de estudio, teniendo en cuenta tanto las covariables como la dependencia espacial.

Herramientas de visualización:

Existen varias herramientas de software y bibliotecas que puedes utilizar para crear mapas a partir de tus predicciones espaciales INLA. Algunos ejemplos populares incluyen:

  • ggplot2 (en R)

  • folium (en Python)

  • ArcGIS Pro

  • QGIS

# again get the correct indices
id_alleff <- inla.stack.index(all_stack, "prd_alleff")$data

# now add the model predictions
newdat$pred <- m_inla$summary.fitted.values[id_alleff, "mean"]
newdat$sd <- m_inla$summary.fitted.values[id_alleff, "sd"]
# get lower and upper confidence interval
newdat$lower_ci <- with(newdat, pred - 2 * sd)
newdat$upper_ci <- with(newdat, pred + 2 * sd)

# some data wraggling
nn <- pivot_longer(newdat, cols = c("pred", "lower_ci", "upper_ci"))

ggplot(nn, aes(x=x, y=y, fill=value)) +
  geom_raster() +
  facet_wrap(~name) +
  scale_fill_continuous(type = "viridis")

Ejemplo con INLA#

library(sf)
library(rnaturalearth)
library(sf)
library(terra)
library(geodata)
library(INLA)
library(rasterVis)

Los promedios anuales de los niveles de concentración de PM2.5 registrados en 1429 estaciones de monitoreo de la Agencia de Protección Ambiental de los Estados Unidos en 2022 se encuentran en el archivo PM25USA2022.csv, que se puede descargar desde este sitio web. Usamos la función read.csv() para leer los datos, que contienen los valores de longitud y latitud de las estaciones de monitoreo, y los valores registrados de PM2.5 en microgramos por metro cúbico. Luego, usamos la función st_as_sf() para transformar el data.frame obtenido en un objeto sf con CRS geográfico dado por el código EPSG 4326.


f <- file.path("https://www.paulamoraga.com/book-spatial/", "data/PM25USA2022.csv")
d <- read.csv(f)
d <- st_as_sf(d, coords = c("longitude", "latitude"))
st_crs(d) <- "EPSG:4326"

Luego, obtenemos el mapa de los EE. UU. con la función ne_countries() de rnaturalearth. Usamos st_crop() para eliminar Alaska y otras áreas que están fuera de la región comprendida por los valores de longitud (–130, 60) y los valores de latitud (18, 72).

map <- ne_countries(type = "countries",
                    country = "United States of America",
                    scale = "medium", returnclass = "sf")
map <- st_crop(map, xmin = -130, xmax = -60, ymin = 18, ymax = 72)

d <- st_filter(d, map)
nrow(d)

library(ggplot2)
library(viridis)
ggplot() + geom_sf(data = map) +
  geom_sf(data = d, aes(col = value)) +
  scale_color_viridis()

Aquí, construimos una matriz coop con las ubicaciones donde se predecirán los niveles de contaminación del aire. Primero, creamos una cuadrícula ráster de 100 × 100 celdas que cubre el mapa usando la función rast() de terra. Luego, obtenemos las coordenadas de las celdas con la función xyfromCell() de terra.

# raster grid covering map
grid <- terra::rast(map, nrows = 100, ncols = 100)
# coordinates of all cells
xy <- terra::xyFromCell(grid, 1:ncell(grid))

Luego, usamos la función st_as_sf() para crear un objeto sf con las coordenadas de las ubicaciones de predicción, especificando las coordenadas como un data frame, el nombre de las coordenadas y el CRS. Obtenemos los índices de las coordenadas de los puntos que están dentro del mapa con st_intersects() configurando sparse = FALSE. Luego usaremos estos índices para identificar las ubicaciones de predicción. También obtenemos las coordenadas de los puntos que están dentro del mapa con sf_filter(). La Figura muestra las ubicaciones de predicción.

# transform points to a sf object
dp <- st_as_sf(as.data.frame(xy), coords = c("x", "y"),
                 crs = st_crs(map))

# indices points within the map
indicespointswithin <- which(st_intersects(dp, map,
                                           sparse = FALSE))

# points within the map
dp <- st_filter(dp, map)

# plot
ggplot() + geom_sf(data = map) +
  geom_sf(data = dp)

En nuestro modelo, utilizamos la temperatura promedio y la precipitación como covariables. Los valores mensuales de estas variables a nivel global pueden obtenerse con la función worldclim_global() de geodata.

covtemp <- worldclim_global(var = "tavg", res = 10,
                            path = tempdir())
covprec <- worldclim_global(var = "prec", res = 10,
                            path = tempdir())

Después de descargar los datos, calculamos los promedios mensuales y extraemos los valores en las ubicaciones de observación y predicción con la función extract() de terra.

# Extract at observed locations
d$covtemp <- extract(mean(covtemp), st_coordinates(d))[, 1]
d$covprec <- extract(mean(covprec), st_coordinates(d))[, 1]
# Extract at prediction locations
dp$covtemp <- extract(mean(covtemp), st_coordinates(dp))[, 1]
dp$covprec <- extract(mean(covprec), st_coordinates(dp))[, 1]

library("patchwork")
p1 <- ggplot() + geom_sf(data = map) +
  geom_sf(data = d, aes(col = covtemp)) +
  scale_color_viridis()
p2 <- ggplot() + geom_sf(data = map) +
  geom_sf(data = d, aes(col = covprec)) +
  scale_color_viridis()
p1/p2

Los datos con los que estamos trabajando tienen un CRS geográfico que referencia ubicaciones utilizando valores de longitud y latitud. Para trabajar con kilómetros en lugar de grados, usamos st_transform() para transformar el CRS de los objetos sf con los datos correspondientes a las ubicaciones observadas (d) y las ubicaciones de predicción (dp) de geográfico a un CRS proyectado. Específicamente, utilizamos la proyección de Mercator dada por el código EPSG 3857 y usamos kilómetros como unidades. Para ello, usamos la proyección dada por st_crs("EPSG:3857")$proj4string, reemplazando +units=m por +units=km.

st_crs("EPSG:3857")$proj4string
projMercator<-"+proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0
+x_0=0 +y_0=0 +k=1 +units=km +nadgrids=@null +wktext +no_defs"
d <- st_transform(d, crs = projMercator)
dp <- st_transform(dp, crs = projMercator)

# Observed coordinates
coo <- st_coordinates(d)

# Predicted coordinates
coop <- st_coordinates(dp)

Ahora especificamos el modelo que utilizamos para predecir los valores de PM2.5 en ubicaciones no muestreadas. Suponemos que \(Yᵢ\), los valores de PM2.5 medidos en las ubicaciones \(i = 1, …, n\), pueden modelarse como

\(Yᵢ ∼ N(μᵢ, σ²)\),

\(μᵢ = β₀ + β₁ × tempᵢ + β₂ × precᵢ + S(x\)ᵢ),

donde \(β₀\) es la intersección, y \(β₁\) \(β₂\) son, respectivamente, los coeficientes de temperatura y precipitación. \(S(⋅)\) es un efecto aleatorio espacial que se modela como un proceso gaussiano de media cero con función de covarianza de Matérn.

summary(dist(coo)) # summary of distances between locations

Para ajustar el modelo utilizando el enfoque SPDE, primero creamos una malla triangulada que cubre la región de estudio donde aproximamos el campo aleatorio gaussiano como un campo aleatorio de Markov gaussiano. INLA produce buenas aproximaciones utilizando una malla fina compuesta por triángulos muy pequeños y con una gran distancia de separación entre las ubicaciones y el borde de la malla para evitar efectos de borde que aumentan la varianza cerca del límite. En algunas aplicaciones, el uso de una malla tan fina podría ser computacionalmente intensivo, y normalmente trabajamos con mallas que aún producen buenas aproximaciones consistentes en una región interna con triángulos pequeños donde se necesita precisión, y una extensión externa con triángulos más grandes donde no se necesitan aproximaciones precisas.

Aquí, creamos la malla con la función inla.mesh.2d() de R-INLA. Pasamos como argumentos loc = coo con las coordenadas de ubicación, y max.edge = c(200, 500) con las longitudes máximas permitidas de los bordes de los triángulos en la región y la extensión para tener triángulos más pequeños dentro de la región que en la extensión. También especificamos cutoff = 1 con la distancia mínima permitida entre puntos para evitar la construcción de muchos triángulos pequeños en áreas donde las ubicaciones están cerca unas de otras (Figura 15.4). El número de vértices de la malla se puede obtener con mesh$n, y la malla se puede graficar de la siguiente manera.

mesh <- inla.mesh.2d(loc = coo, max.edge = c(200, 500),
                     cutoff = 1)
mesh$n

plot(mesh)
points(coo, col = "red")
axis(1)
axis(2)

Luego, usamos la función inla.spde2.matern() para construir el modelo SPDE. Esta función tiene parámetros mesh con la malla triangulada construida y constr = TRUE para imponer una restricción de integración a cero. Además, establecemos el parámetro de suavidad \(ν\) igual a 1. En el caso espacial \(d = 2\) y \(α = ν + d / 2 = 2\).

spde <- inla.spde2.matern(mesh = mesh, alpha = 2, constr = TRUE)

Luego, creamos un conjunto de índices para el modelo SPDE usando la función inla.spde.make.index(), donde proporcionamos el nombre del efecto (s) y el número de vértices en el modelo SPDE (spde$n.spde). Esta función genera una lista con el vector s que va del 1 al spde$n.spde. Además, crea dos vectores, s.group y s.repl, que contienen todos los elementos establecidos en 1 y longitudes iguales al número de vértices de la malla.

indexs <- inla.spde.make.index("s", spde$n.spde)
lengths(indexs)

Usamos la función inla.spde.make.A() de R-INLA, pasando la malla (mesh) y las coordenadas (coo) para construir fácilmente una matriz de proyección A que proyecta el campo aleatorio gaussiano espacialmente continuo desde las observaciones hasta los nodos de la malla.

A <- inla.spde.make.A(mesh = mesh, loc = coo)

Podemos ver que la matriz de proyección \(A\) tiene un número de filas igual al número de observaciones y un número de columnas igual al número de vértices de la malla. También observamos que los elementos de cada fila de \(A\) suman 1.

# dimension of the projection matrix
dim(A)

También creamos una matriz de proyección para las ubicaciones de predicción.

Ap <- inla.spde.make.A(mesh = mesh, loc = coop)

Ahora creamos un stack con los datos para estimación y predicción que organiza datos, efectos y matrices de proyección. Creamos stacks para estimación (stk.e) y predicción (stk.p) usando tag para identificar el tipo de datos, data con la lista de vectores de datos, A con las matrices de proyección, y effects con una lista de efectos fijos y aleatorios. Primero, creamos un stack llamado stk.e que contiene los datos para estimación, el cual está etiquetado con la cadena “est”. En data, especificamos el vector de respuesta con los valores observados de PM2.5. La matriz de proyección se da en el argumento A, que es una lista donde el segundo elemento es la matriz de proyección para los efectos aleatorios (A) y el primer elemento se establece en 1 para indicar que los efectos fijos se mapean directamente uno a uno a la respuesta. Para definir los efectos, pasamos una lista que contiene los efectos fijos y aleatorios. Los efectos fijos son un data.frame que consta de una intersección (b0) y covariables de temperatura (covtemp) y precipitación (covprec). El efecto aleatorio está representado por el campo aleatorio gaussiano espacial s que contiene una lista con los índices del objeto SPDE (indexs). Además, construimos otro stack llamado stk.p para predicción, que está etiquetado con la etiqueta “pred”. Los datos, la matriz de proyección y los efectos se especifican para las ubicaciones de predicción. El vector de respuesta en el argumento data de este stack se establece en una lista con NA porque estos son los valores que queremos predecir. Finalmente, combinamos stk.e y stk.p en un solo stack completo llamado stk.full.

# stack for estimation stk.e
stk.e <- inla.stack(tag = "est",
data = list(y = d$value), A = list(1, A),
effects = list(data.frame(b0 = rep(1, nrow(A)),
covtemp = d$covtemp, covprec = d$covprec),
s = indexs))

# stack for prediction stk.p
stk.p <- inla.stack(tag = "pred",
data = list(y = NA), A = list(1, Ap),
effects = list(data.frame(b0 = rep(1, nrow(Ap)),
covtemp = dp$covtemp, covprec = dp$covprec),
s = indexs))

# stk.full has stk.e and stk.p
stk.full <- inla.stack(stk.e, stk.p)

Luego, especificamos la fórmula incluyendo la variable de respuesta, el símbolo ~ y los efectos fijos y aleatorios. En la fórmula, eliminamos la intersección añadiendo 0 e incluimos la intersección como un término de covariable añadiendo b0. Este paso asegura que todos los términos de covariable se capturen adecuadamente dentro de la matriz de proyección.

formula <- y ~ 0 + b0 + covtemp + covprec + f(s, model = spde)

Finalmente, llamamos a inla() especificando la fórmula, la familia, el stack con los datos y las opciones. Establecemos control.predictor = list(compute = TRUE) y control.compute = list(return.marginals.predictor = TRUE) para calcular y devolver las marginales para el predictor lineal.

res <- inla(formula, family = "gaussian",
       data = inla.stack.data(stk.full),
       control.predictor = list(compute = TRUE,
                                A = inla.stack.A(stk.full)),
       control.compute = list(return.marginals.predictor = TRUE))

res$summary.fixed

Observamos que el coeficiente de temperatura es \(^β₁ = 0.239\) con un intervalo creíble del 95% igual a (0.201, 0.28). El coeficiente de precipitación es \(^β₂ = 0.003\) con un intervalo creíble del 95% igual a (–0.003, 0.009). Así, la temperatura está significativamente asociada con PM2.5, mientras que la precipitación no es significativa.

El objeto res$summary.fitted.values contiene la media posterior y los cuantiles de los valores ajustados. Podemos obtener los índices correspondientes a las ubicaciones de predicción utilizando la función inla.stack.index() pasando el stack completo y tag = "pred". Luego, recuperamos la columna “mean” con la media posterior, y las columnas “0.025quant” y “0.975quant” con los límites inferior y superior de los intervalos creíbles del 95% que denotan la incertidumbre de las predicciones.

index <- inla.stack.index(stack = stk.full, tag = "pred")$data
pred_mean <- res$summary.fitted.values[index, "mean"]
pred_ll <- res$summary.fitted.values[index, "0.025quant"]
pred_ul <- res$summary.fitted.values[index, "0.975quant"]

Asignamos los valores predichos a sus celdas correspondientes dentro del mapa que están en el objeto grid que contiene las ubicaciones de predicción.

grid$mean <- NA
grid$ll <- NA
grid$ul <- NA

grid$mean[indicespointswithin] <- pred_mean
grid$ll[indicespointswithin] <- pred_ll
grid$ul[indicespointswithin] <- pred_ul

summary(grid) # negative values for the lower limit

Luego, graficamos la media posterior y los intervalos creíbles del 95% de los valores predichos de PM2.5 con la función levelplot() del paquete rasterVis. La Figura 15.5 muestra mapas con el patrón espacial de los niveles predichos de PM2.5 así como su incertidumbre asociada.


levelplot(grid, layout = c(1, 3),
names.attr = c("Mean", "2.5 percentile", "97.5 percentile"))

También podemos obtener probabilidades de que PM2.5 supere un umbral específico con la función inla.pmarginal(). Específicamente, calculamos las probabilidades de que los niveles de PM2.5 superen los 10 microgramos por metro cúbico. Es decir, P(PM2.5 > 10) = 1 – P(PM2.5 ≤ 10).

excprob <- sapply(res$marginals.fitted.values[index],
FUN = function(marg){1-inla.pmarginal(q = 10, marginal = marg)})

Luego, añadimos las probabilidades de excedencia como una capa en grid, y las graficamos con levelplot(). En levelplot(), establecemos margin = FALSE para ocultar los gráficos marginales de los resúmenes de columnas y filas del objeto raster. La Figura muestra las probabilidades de que los niveles de PM2.5 superen los 10 microgramos por metro cúbico. Observamos altas probabilidades en la costa oeste y en la parte sur del país.

grid$excprob <- NA
grid$excprob[indicespointswithin] <- excprob

levelplot(grid$excprob, margin = FALSE)
Contents