Modelo GLM de Poisson no-homogéneo

Modelo GLM de Poisson no-homogéneo#

Los modelos GLM de Poisson no-homogeneo, también llamados Log-Gaussianos Cox (LGCPs), se refieren a que la intensidad $\lambda(s)$ varía en el espacio. Un LGCP es un proceso de Poisson con una intensidad variable. Dada un area $A$, la probabilidad de observar un cierto número de puntos en dicha área sigue una distribución de Poisson con intensidad (valor esperado):

\[ \lambda_A=\int{\lambda(s)ds} \]

La parte gaussiana de LGCP viene de la modelación $log(\lambda(s))$ como un Gaussiano latente (condicionado a un conjunto de parámetros), en un marco típico de GLM or GAM.

El proceso de Cox log-Gaussiano (LGCP) es un modelo probabilístico de patrones de puntos que se observan típicamente en el espacio o el tiempo. Tiene dos componentes principales. Primero, un campo de intensidad subyacente $\lambda(s)$ de valores reales positivos se modela sobre todo el dominio $X$ utilizando un proceso Gaussiano transformado exponencialmente, lo que obliga a que $\lambda$ sea positivo. Luego, este campo de intensidad se utiliza para parametrizar un proceso de puntos de Poisson, que representa un mecanismo estocástico para colocar puntos en el espacio. Algunos fenómenos que pueden representarse con este modelo incluyen la incidencia de casos de cáncer en un condado o las ubicaciones espaciotemporales de eventos delictivos en una ciudad. Tanto las dimensiones espaciales como temporales se pueden manejar de manera equivalente dentro de este marco, aunque este tutorial solo aborda datos en dos dimensiones espaciales.

En términos más formales, si tenemos un espacio $X$ y $A\subseteq X$, la distribución sobre el número de puntos $Y_A$ que ocurren dentro del subconjunto $A$ se define como:

\[Y_A \sim Poisson\left(\int_A \lambda(s) ds\right)\]

y el campo de intensidad se define como

\[\log \lambda(s) \sim GP(\mu(s), K(s,s'))\]

donde $GP(\mu(s), K(s,s'))$ denota un proceso Gaussiano con función de media $\mu(s)$ y núcleo de covarianza $K(s,s')$ para una ubicación $s \in X$.

Un LGCP es un proceso estocástico que tambien puede expresarce de la forma:

$\lambda_s = e^{Z(s)}$,

donde Z = {Z(s) : s ∈ R²} es un proceso Gaussiano. Entonces, condicionado a $\lambda(s)$, el proceso de puntos es un proceso de Poisson no-homogéneo con intensidad $\lambda(s)$. Esto implica que el número de eventos en cualquier región A se distribuye según una Poisson con media $\int_A \lambda(s)$ ds, y las ubicaciones de los eventos son una muestra aleatoria independiente de la distribución en A con una densidad de probabilidad proporcional a $\lambda(s)$. Un modelo LGCP también puede incluir variables explicativas espaciales, proporcionando un enfoque flexible para describir y predecir una amplia gama de fenómenos espaciales.

Ejemplo 1#

En este notebook, asumimos que hemos observado un patrón espacial de puntos de ubicaciones de eventos {x_i : i = 1, …, n} que ha sido generado como una realización de un LGCP, y mostramos cómo ajustar un modelo LGCP a los datos utilizando los enfoques INLA y SPDE. El Capítulo 15 introdujo el enfoque SPDE y describió su implementación en el contexto de la geoestadística basada en modelos utilizando un ejemplo de predicción de contaminación del aire. Aquí, describimos cómo usar SPDE para ajustar un modelo LGCP a un patrón de puntos de especies de plantas para estimar la intensidad del proceso.

library("sf")
library("spocc")
library("ggplot2")
library(sf)
library(terra)
library(rnaturalearth)
library(INLA)
library(rgeos)

devtools::install_github("cran/rgeos")

En este ejemplo, estimamos la intensidad de las especies de plantas de Solanum en Bolivia desde enero de 2015 hasta diciembre de 2022, obtenidas de la base de datos Global Biodiversity Information Facility (GBIF) con el paquete spocc. Recuperamos los datos utilizando la función occ(), especificando el nombre científico de la especie de planta, la fuente de datos, las fechas y el código del país. También especificamos has_coords = TRUE para devolver solo las ocurrencias que tienen coordenadas, y limit = 1000 para establecer el límite en el número de registros. Aquí, mostramos cómo formular y ajustar un modelo LGCP para las especies de plantas de Solanum en Bolivia utilizando un campo aleatorio Gaussiano continuo con INLA y SPDE. El modelo nos permite estimar la intensidad del proceso que genera las ubicaciones.

df <- occ(query = "solanum", from = "gbif",
          date = c("2015-01-01", "2022-12-31"),
          gbifopts = list(country = "BO"),
          has_coords = TRUE, limit = 1000)
d <- occ2df(df)

Utilizamos la función st_as_sf() para crear un objeto sf llamado d que contiene las nrow(d) = 263 ubicaciones recuperadas. Establecemos el sistema de referencia de coordenadas (CRS) al código EPSG 4326 ya que las coordenadas de las ubicaciones están dadas por valores de longitud y latitud.

d <- st_as_sf(d[, 2:3], coords = c("longitude", "latitude"))
st_crs(d) <- "EPSG:4326"

Para trabajar con kilómetros en lugar de grados, proyectamos los datos a UTM 19S correspondiente al código EPSG 5356 con kilómetros como unidades. Para hacerlo, obtenemos st_crs("EPSG:5356")$proj4string y cambiamos +units=m por +units=km.

st_crs("EPSG:5356")$proj4string
projUTM <- "+proj=utm +zone=19 +south +ellps=GRS80
+towgs84=0,0,0,0,0,0,0 +units=km +no_defs"
d <- st_transform(d, crs = projUTM)

También obtenemos el mapa de Bolivia con el paquete rnaturalearth y lo proyectamos a UTM 19S con kilómetros como unidades.

map <- ne_countries(type = "countries", country = "Bolivia",
                    scale = "medium", returnclass = "sf")
map <- st_transform(map, crs = projUTM)

ggplot() + geom_sf(data = map) +
  geom_sf(data = d) + coord_sf(datum = projUTM)

Finalmente, creamos un data frame llamado coo con las ubicaciones de los eventos.

coo <- st_coordinates(d)

Ahora, construimos una matriz con las ubicaciones coop donde queremos predecir la intensidad del proceso puntual. Para hacerlo, primero creamos un raster que cubra el mapa utilizando la función rast() del paquete terra. Luego, recuperamos las coordenadas de las celdas con la función xyFromCell() del paquete terra.

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

Creamos un objeto sf llamado dp con las ubicaciones de predicción utilizando st_as_sf(), y usamos st_filter() para conservar las ubicaciones de predicción que se encuentran dentro del mapa. También recuperamos los índices de los puntos dentro del mapa utilizando st_intersects(), configurando sparse = FALSE.

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

ggplot() + geom_sf(data = map) +
  geom_sf(data = dp) + coord_sf(datum = projUTM)

coop <- st_coordinates(dp)

Utilizamos un LGCP para modelar el patrón de puntos de las especies de plantas. Así, asumimos que el proceso que origina las ubicaciones de las especies de plantas es un proceso de Poisson con una intensidad variable expresada como:

$log(Λ(s)) = β₀ + Z(s),$

donde β₀ es el intercepto y Z(⋅) es un proceso espacial gaussiano de media cero con función de covarianza Matérn.

Para ajustar el modelo utilizando INLA y SPDE, primero construimos una malla. En el análisis de patrones de puntos, generalmente no empleamos las ubicaciones como vértices de la malla. Construimos una malla que cubre la región de estudio utilizando la función inla.mesh.2d(), estableciendo loc.domain igual a una matriz con las ubicaciones de los puntos en el límite de la región. Los otros argumentos son los siguientes: max.edge denota las longitudes máximas permitidas de los bordes de los triángulos en la región interior y la extensión. offset especifica el tamaño de las extensiones interior y exterior alrededor de las ubicaciones de los datos. cutoff es la distancia mínima permitida entre puntos que utilizamos para evitar construir muchos triángulos pequeños alrededor de ubicaciones agrupadas.


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

loc.d <- cbind(st_coordinates(map)[, 1], st_coordinates(map)[, 2])

mesh <- inla.mesh.2d(loc.domain = loc.d, max.edge = c(50, 100), offset = c(50, 100), cutoff = 1)

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

También creamos las variables nv con el número de vértices de la malla y la variable n con el número de eventos del patrón de puntos. Más adelante, utilizaremos estas variables para construir las pilas de datos.

(nv <- mesh$n)

(n <- nrow(coo))

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

Utilizamos la función inla.spde2.matern() para construir el modelo SPDE en la malla.

book.mesh.dual <- function(mesh) {
    if (mesh$manifold=='R2') {
        ce <- t(sapply(1:nrow(mesh$graph$tv), function(i)
            colMeans(mesh$loc[mesh$graph$tv[i, ], 1:2])))
        library(parallel)
        pls <- mclapply(1:mesh$n, function(i) {
            p <- unique(Reduce('rbind', lapply(1:3, function(k) {
            j <- which(mesh$graph$tv[,k]==i)
            if (length(j)>0) 
            return(rbind(ce[j, , drop=FALSE],
            cbind(mesh$loc[mesh$graph$tv[j, k], 1] +
            mesh$loc[mesh$graph$tv[j, c(2:3,1)[k]], 1], 
            mesh$loc[mesh$graph$tv[j, k], 2] +
            mesh$loc[mesh$graph$tv[j, c(2:3,1)[k]], 2])/2))
            else return(ce[j, , drop=FALSE])
            })))
            j1 <- which(mesh$segm$bnd$idx[,1]==i)
            j2 <- which(mesh$segm$bnd$idx[,2]==i)
            if ((length(j1)>0) | (length(j2)>0)) {
            p <- unique(rbind(mesh$loc[i, 1:2], p,
            mesh$loc[mesh$segm$bnd$idx[j1, 1], 1:2]/2 +
            mesh$loc[mesh$segm$bnd$idx[j1, 2], 1:2]/2, 
            mesh$loc[mesh$segm$bnd$idx[j2, 1], 1:2]/2 +
            mesh$loc[mesh$segm$bnd$idx[j2, 2], 1:2]/2))
            yy <- p[,2]-mean(p[,2])/2-mesh$loc[i, 2]/2
            xx <- p[,1]-mean(p[,1])/2-mesh$loc[i, 1]/2
            }
            else {
            yy <- p[,2]-mesh$loc[i, 2]
            xx <- p[,1]-mesh$loc[i, 1]
            }
            Polygon(p[order(atan2(yy,xx)), ])
        })
        return(SpatialPolygons(lapply(1:mesh$n, function(i)
            Polygons(list(pls[[i]]), i))))
    }
    else stop("It only works for R2!")
}

Aquí, creamos los vectores con el número observado y esperado de eventos. Primero, creamos la malla dual que consiste en un conjunto de polígonos alrededor de cada vértice de la malla original. Podemos crear la malla dual utilizando la función book.mesh.dual() que se proporciona en Krainski et al. (2019).

dmesh <- book.mesh.dual(mesh)
plot(dmesh)
axis(1)
axis(2)

Para ajustar el LGCP, los vértices de la malla se consideran como puntos de integración. Los valores esperados correspondientes a los vértices de la malla son proporcionales a las áreas alrededor de los vértices de la malla, es decir, las áreas de los polígonos de la malla dual. Calculamos un vector de pesos llamado w con las áreas de la intersección entre cada polígono de la malla dual y la región de estudio utilizando el siguiente código.

# Domain polygon is converted into a SpatialPolygons
domain.polys <- Polygons(list(Polygon(loc.d)), '0')
domainSP <- SpatialPolygons(list(domain.polys))

# Because the mesh is larger than the study area, we need to
# compute the intersection between each polygon
# in the dual mesh and the study area

w <- sapply(1:length(dmesh), function(i) {
  if (gIntersects(dmesh[i, ], domainSP))
    return(gArea(gIntersection(dmesh[i, ], domainSP)))
  else return(0)
})

sum(w) # sum weights

st_area(map) # area of the study region

plot(mesh)
points(mesh$loc[which(w > 0), 1:2], col = "black", pch = 20)
points(mesh$loc[which(w == 0), 1:2], col = "red", pch = 20)

Luego, creamos vectores de los conjuntos de datos aumentados con los valores observados y esperados. Los valores observados se especificarán en la fórmula del modelo como respuesta. Los valores esperados se especificarán en la fórmula del modelo como el componente E de la media para la verosimilitud de Poisson definida como

$E_i = \exp(\eta_i)$

donde $\eta_i$ es el predictor lineal.

El vector y.pp contiene la variable de respuesta. Los primeros nv elementos son ceros correspondientes a los vértices de la malla. Los últimos n elementos son unos correspondientes a los eventos observados.

El vector e.pp contiene los valores esperados. Los primeros nv elementos son los pesos w que representan la intersección entre las áreas alrededor de cada uno de los vértices de la malla y la región de estudio. Los siguientes n elementos son ceros correspondientes a los eventos puntuales.

y.pp <- rep(0:1, c(nv, n))
e.pp <- c(w, rep(0, n))

head(cbind(y.pp, e.pp))

tail(cbind(y.pp, e.pp))

Construimos la matriz de proyección A.pp para proyectar el campo aleatorio Gaussiano desde las observaciones a los vértices de la triangulación. Esta matriz se construye utilizando la matriz de proyección para los vértices de la malla, que es una matriz diagonal (A.int), y la matriz de proyección para las ubicaciones de eventos (A.y).

# Projection matrix for the integration points (mesh vertices)
A.int <- Diagonal(nv, rep(1, nv))
# Projection matrix for observed points (event locations)
A.y <- inla.spde.make.A(mesh = mesh, loc = coo)
# Projection matrix for mesh vertices and event locations
A.pp <- rbind(A.int, A.y)

También creamos la matriz de proyección Ap.pp para las ubicaciones de predicción.

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

Ahora usamos la función inla.stack() para construir pilas para estimación y predicción que organizan los datos, efectos y matrices de proyección. En los argumentos de inla.stack(), data es una lista con los valores observados (y) y esperados (e). El argumento A contiene las matrices de proyección, y el argumento effects es una lista con los efectos fijos y aleatorios. Luego, las pilas de estimación y predicción se combinan en una pila completa.

# stack for estimation
stk.e.pp <- inla.stack(tag = "est.pp",
data = list(y = y.pp, e = e.pp), 
A = list(1, A.pp),
effects = list(list(b0 = rep(1, nv + n)), list(s = 1:nv)))

# stack for prediction stk.p
stk.p.pp <- inla.stack(tag = "pred.pp",
data = list(y = rep(NA, nrow(coop)), e = rep(0, nrow(coop))),
A = list(1, Ap.pp),
effects = list(data.frame(b0 = rep(1, nrow(coop))),
               list(s = 1:nv)))

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

La fórmula se especifica incluyendo la variable de respuesta en el lado izquierdo y los efectos aleatorios en el lado derecho.

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

Ajustamos el modelo llamando a inla(). En la función, especificamos link = 1 para calcular los valores ajustados que se encuentran en res$summary.fitted.values y res$marginals.fitted.values con la misma función de enlace que la familia especificada en el modelo.

res <- inla(formula,  family = 'poisson',
data = inla.stack.data(stk.full.pp),
control.predictor = list(compute = TRUE, link = 1,
                         A = inla.stack.A(stk.full.pp)),
E = inla.stack.data(stk.full.pp)$e)

Un resumen de los resultados se puede inspeccionar escribiendo summary(res). El marco de datos res$summary.fitted.values contiene los valores ajustados. Los índices de las filas correspondientes a las predicciones se pueden obtener con inla.stack.index() especificando la etiqueta “pred.pp” del stack de predicción.

index <- inla.stack.index(stk.full.pp, tag = "pred.pp")$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"]

Luego, añadimos capas a la cuadrícula raster con la media posterior y los percentiles 2.5 y 97.5 en las celdas que están dentro del mapa.

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

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

Finalmente, creamos mapas de la media posterior y los límites inferior y superior de los intervalos creíbles del 95% de la intensidad del proceso puntual de especies en Bolivia (Figura 23.4). Para ello, utilizamos la función levelplot() del paquete rasterVis, especificando names.attr con el nombre de cada panel y layout con el número de columnas y filas.

library(rasterVis)
levelplot(raster::brick(grid), layout = c(3, 1),
names.attr = c("Mean", "2.5 percentile", "97.5 percentile"))

En general, observamos una baja intensidad de especies, con mayor intensidad en la parte centro-occidental de Bolivia. Cabe destacar que hemos modelado los datos de ocurrencia de especies recuperados de GBIF asumiendo que el patrón puntual espacial observado es una realización del proceso subyacente que genera las ubicaciones de las especies. En aplicaciones reales, es importante entender cómo se recolectaron los datos y evaluar posibles sesgos en los datos, como la sobrerepresentación de ciertas áreas, que pueden invalidar los resultados. Además, es importante incorporar el conocimiento de expertos para crear modelos que incluyan covariables relevantes y efectos aleatorios para tener en cuenta diversos tipos de variabilidad, lo que permite una comprensión más completa de la variable en investigación.

Ejemplo 2#

Un patrón de puntos registra la ocurrencia de eventos en una región de estudio. Ejemplos típicos incluyen la ubicación de árboles en un bosque o las coordenadas GPS de casos de enfermedades en una región. Las ubicaciones de los eventos observados dependen de un proceso espacial subyacente, que a menudo se modela utilizando una función de intensidad $\lambda(s)$. La función de intensidad mide el número promedio de eventos por unidad de espacio, y puede ser modelada para depender de covariables y otros efectos.

Bajo el supuesto del modelo de proceso de puntos log-Cox, modelamos la log-intensidad del proceso de Cox con un predictor lineal Gaussiano. En este caso, el proceso log-Cox se conoce como un proceso Cox log-Gaussiano (LGCP, Møller, Syversveen y Waagepetersen 1998), y la inferencia se puede realizar utilizando INLA. Un proceso de Cox es simplemente un nombre para un proceso de Poisson con intensidad variable; por lo tanto, utilizamos la verosimilitud de Poisson. El enfoque original utilizado para ajustar estos modelos en INLA (y en otro software) divide la región de estudio en celdas, que forman una cuadrícula, y cuenta el número de puntos en cada una (Møller y Waagepetersen 2003). Estos conteos pueden modelarse utilizando una verosimilitud de Poisson condicionada a un predictor lineal Gaussiano, y se puede usar INLA para ajustar el modelo (Illian, Sørbye y Rue 2012).

En este capítulo nos enfocamos en un nuevo enfoque que considera modelos Stochastic Partial Differential Equation (SPDE) directamente, desarrollado en Simpson et al. (2016). Este enfoque tiene una buena justificación teórica y considera una aproximación directa de la verosimilitud del modelo de proceso de puntos log-Cox. Las observaciones se modelan considerando su ubicación exacta en lugar de agruparlas en celdas. Junto con la flexibilidad para definir una malla, este enfoque puede manejar áreas no rectangulares sin desperdiciar esfuerzo computacional en un área rectangular grande.

library(spatstat)
library(RandomFields)
library(rgeos)

La función rLGCP utiliza la función GaussRF() del paquete RandomFields (Schlather et al. 2015) para simular desde el campo espacial sobre una cuadrícula. Hay un parámetro interno para controlar la resolución de la cuadrícula, que especificamos para dar 300 píxeles en cada dirección:

win <- owin(c(0, 3), c(0, 3))

Modelamos la intensidad como

\[ \log(\lambda(s)) = \beta_0 + S(s), \]

donde $\beta_0$ es un valor fijo y $S(s)$ es un proceso espacial Gaussiano con covarianza de Matérn y media cero. El parámetro $\beta_0$ puede considerarse como un nivel medio global para la log-intensidad; es decir, la log-intensidad fluctúa alrededor de este valor según el proceso espacial $S(s)$.

Si no hay campo espacial, el número esperado de puntos es $e^{\beta_0}$ multiplicado por el área de la ventana. Esto significa que el número esperado de puntos es:

beta0 <- 3
exp(beta0) * diff(range(win$x)) * diff(range(win$y))

Por lo tanto, este valor de $\beta_0$ producirá un número razonable de puntos en las siguientes simulaciones. Si establecemos $\beta_0$ demasiado bajo, obtendremos casi ningún punto, y no seremos capaces de producir resultados razonables. También es posible utilizar una función en varias covariables, por ejemplo, un GLM.

En este taller, utilizamos una función de covarianza de Matérn con $\nu = 1$. Los otros parámetros son la varianza y la escala. Los siguientes valores para estos parámetros producirán una intensidad suave del proceso puntual:

sigma2x <- 0.2
range <- 1.2
nu <- 1

El valor de $\sigma^2_x$ se establece para hacer que la log-intensidad varíe un poco alrededor de la media, pero siempre dentro de un rango razonable de valores. Además, con estos parámetros $\nu = 1$ y el rango del proceso espacial $S(s)$ es (aproximadamente) 2, lo que produce cambios suaves en la ventana de estudio actual. Valores más pequeños del rango práctico producirán un proceso espacial $S(s)$ (y, a su vez, la intensidad del proceso espacial) que cambia rápidamente en la ventana de estudio. De manera similar, valores muy grandes del rango práctico producirán un proceso espacial casi constante $S(s)$, de modo que la log-intensidad estará muy cerca de $\beta_0$ en todos los puntos de la ventana de estudio.

Los puntos del proceso puntual se simulan de la siguiente manera:

set.seed(1)
lg.s <- rLGCP('matern', beta0, var = sigma2x,scale = range / sqrt(8), nu = nu, win = win)

Tanto el campo espacial como el patrón de puntos son devueltos. Las coordenadas de los eventos observados del patrón de puntos se pueden obtener de la siguiente manera:

xy <- cbind(lg.s$x, lg.s$y)[, 2:1]

El numero de puntos simulados es:

(n <- nrow(xy))

La exponencial de los valores simulados del campo espacial se devuelve como el atributo Lambda del objeto. A continuación, extraemos los valores de λ(s) y resumimos el log(λ(s)).

Lam <- attr(lg.s, 'Lambda')
rf.s <- log(Lam$v)
summary(as.vector(rf.s))

Siguiendo a Simpson et al. (2016), los parámetros del modelo de proceso puntual de Cox log-Gaussiano pueden estimarse con INLA.

INLA es un paquete que permite ajustar una amplia gama de modelos. Utiliza la aproximación de Laplace para ajustar modelos bayesianos de manera mucho más rápida que algoritmos como MCMC (Cadena de Markov Monte Carlo). INLA permite ajustar modelos geoestadísticos a través de ecuaciones diferenciales parciales estocásticas (EDPE). Puedes encontrar más información sobre esto en estos dos Gitbooks:

spde-gitbook [https://becarioprecario.bitbucket.io/spde-gitbook]
inla-gitbook [https://becarioprecario.bitbucket.io/inla-gitbook]

Ajustar un modelo espacial en INLA requiere un conjunto específico de pasos:

Crear una malla: Se crea una malla para aproximar el efecto espacial. Esta malla discretiza el espacio de estudio en pequeñas unidades (píxeles o celdas) que permiten modelar la variación espacial continua.
Crear una matriz de proyección: Esta matriz vincula las observaciones puntuales a la malla creada. Básicamente, indica a qué celda de la malla pertenece cada observación.
Definir la ecuación diferencial parcial estocástica (EDPE): La EDPE describe la relación espacial entre las variables del modelo. Diferentes tipos de EDPE capturan diferentes patrones de dependencia espacial.
Especificar opcionalmente un conjunto de datos para predicciones: Si deseas realizar predicciones en ubicaciones no observadas, puedes especificar un conjunto de datos adicional que contenga las coordenadas de los puntos de predicción.
Apilar los objetos en un objeto stack: INLA trabaja con un objeto especial llamado “stack” que combina la información de la malla, la matriz de proyección, la EDPE y cualquier otro dato relevante para el modelo.
Ajustar el modelo: Una vez creado el objeto stack, se utiliza la función inla() del paquete INLA para ajustar el modelo y obtener los resultados de la inferencia bayesiana.

library("devtools")
devtools::install_github(repo = "https://github.com/hrue/r-inla", ref = "stable", subdir = "rinla", build = FALSE)

inla.upgrade() # for the stable version

library(INLA)

En términos simplificados, construiremos un conjunto de datos ampliado y ejecutaremos una regresión de Poisson con INLA. El conjunto de datos ampliado se compone de una respuesta binaria, con 1 para los puntos observados y 0 para algunas observaciones ficticias. Tanto las observaciones reales como las ficticias tendrán valores “esperados” o pesos asociados que se incluirán en la regresión de Poisson. Esto se explicará paso a paso en las siguientes secciones.

Para una inferencia adecuada con el LGCP, debemos tener cuidado al construir la malla. En el caso particular del análisis de patrones de puntos, generalmente no usamos los puntos de ubicación como nodos de la malla. Necesitamos una malla que cubra la región de estudio; para esto utilizamos loc.domain para construir la malla. Además, solo usamos una pequeña primera extensión exterior, pero no una segunda extensión exterior.

loc.d <- 3 * cbind(c(0, 1, 1, 0, 0), c(0, 0, 1, 1, 0))
mesh <- inla.mesh.2d(loc.domain = loc.d, offset = c(0.3, 1), 
  max.edge = c(0.3, 0.7), cutoff = 0.05)
nv <- mesh$n

El modelo SPDE se definirá considerando los PC-priors derivados en Fuglstad et al. (2018) para los parámetros del modelo: rango y desviación estándar marginal. Estos se definen de la siguiente manera:

spde <- inla.spde2.pcmatern(mesh = mesh,
  # PC-prior on range: P(practic.range < 0.05) = 0.01
  prior.range = c(0.05, 0.01),
  # PC-prior on sigma: P(sigma > 1) = 0.01
  prior.sigma = c(1, 0.01)) 

El enfoque SPDE para el análisis de patrones puntuales define el modelo en los nodos de la malla. Para ajustar el modelo de proceso puntual de log-Cox, estos puntos se consideran como puntos de integración. El método en Simpson et al. (2016) define que el número esperado de eventos es proporcional al área alrededor del nodo (las áreas de los polígonos en la malla dual). Esto significa que en los nodos de la malla con triángulos más grandes, también hay valores esperados más grandes. El comando inla.mesh.fem(mesh)$va proporciona este valor para cada nodo de la malla. Estos valores para los nodos en el dominio interno se pueden usar para calcular la intersección entre los polígonos de la malla dual y el polígono del dominio de estudio. Para ello, utilizamos la función book.mesh.dual():

dmesh <- book.mesh.dual(mesh)

Esta función está disponible en el archivo spde-book-functions.R y devuelve la malla dual en un objeto de la clase SpatialPolygons. El polígono del dominio se puede convertir en una clase SpatialPolygons de la siguiente manera:

domain.polys <- Polygons(list(Polygon(loc.d)), '0')
domainSP <- SpatialPolygons(list(domain.polys))

El vector de pesos que hemos calculado es exactamente lo que necesitamos usar como exposición (E) en la verosimilitud de Poisson en INLA (con la pequeña modificación de que log(E) se define como cero si E=0. Aumentamos el vector de unos para las observaciones (que representan los puntos) con una secuencia de ceros (que representan los nodos de la malla):

w <- sapply(1:length(dmesh), function(i) {
  if (gIntersects(dmesh[i, ], domainSP))
    return(gArea(gIntersection(dmesh[i, ], domainSP)))
  else return(0)
})

El vector de exposición y la matriz de proyección se definen. Para los puntos de integración, esta es simplemente una matriz diagonal porque estos lugares son solo los vértices de la malla. La matriz de proyección completa es:

y.pp <- rep(0:1, c(nv, n))
e.pp <- c(w, rep(0, n)) 
imat <- Diagonal(nv, rep(1, nv))
lmat <- inla.spde.make.A(mesh, xy)
A.pp <- rbind(imat, lmat)

Configuramos el conjunto de datos de la siguiente manera:

stk.pp <- inla.stack(
  data = list(y = y.pp, e = e.pp), 
  A = list(1, A.pp),
  effects = list(list(b0 = rep(1, nv + n)), list(i = 1:nv)),
  tag = 'pp')

Las marginales posteriores para todos los parámetros del modelo se obtienen ajustando el modelo con INLA:

pp.res <- inla(y ~ 0 + b0 + f(i, model = spde), 
  family = 'poisson', data = inla.stack.data(stk.pp), 
  control.predictor = list(A = inla.stack.A(stk.pp)), 
  E = inla.stack.data(stk.pp)$e)

pp.res$summary.hyperpar

Modelo GLM de Poisson no-homogéneo

Contents

Modelo GLM de Poisson no-homogéneo#

Ejemplo 1#

Ejemplo 2#