INLA mesh

 library(INLA)

library(sf)

library(gstat)

library(tidyverse)

data(sic97)

# 直接在一个流程中准备全部数据

df_rain <- sic_full %>% 

  st_as_sf() %>% 

  mutate(

    altitude = raster::extract(raster::raster(demstd), as.matrix(st_coordinates(.))),

    coords = st_coordinates(.)

  ) %>% 

  mutate(

    lat = coords[, 1] / 1000,

    lon = coords[, 2] / 1000

  ) %>% 

  dplyr::select(ID, rainfall, altitude, lat, lon)

# 检查空间依赖性

df_rain %>%  

  as_Spatial() %>%

  variogram(rainfall ~ 1, data = .) %>%

  graphics::plot()

# 创建网格

loc <- cbind(df_rain$lat, df_rain$lon)

Mesh <- inla.mesh.2d(

  loc.domain = loc,

  max.edge = c(20, 100),

  cutoff = 1

)

ggplot() +

  inlabru::gg(Mesh) +

  geom_point(data = data.frame(x = loc[,1], y = loc[,2]),

             aes(x = x, y = y), color = "red", shape = 1) +

  theme_minimal() +

  labs(title = "INLA Mesh with Data Points",

       x = "Latitude (km)", y = "Longitude (km)") +

  coord_equal()

# SPDE设置并直接创建堆栈

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

s.index <- inla.spde.make.index(name = "w", n.spde = spde$n.spde)

# 直接创建堆栈

stack_data <- inla.stack(

  tag = "data",

  data = list(y = df_rain$rainfall),

  A = list(1, 1, inla.spde.make.A(Mesh, loc = loc)),

  effects = list(

    Intercept = rep(1, nrow(df_rain)),

    X = model.matrix(~ -1 + altitude, data = df_rain) %>%  as.data.frame(),

    w = s.index

  )

)

# 拟合模型

fit <- inla(

  y ~ -1 + Intercept + altitude + f(w, model = spde),

  data = inla.stack.data(stack_data, spde = spde),

  family = "gaussian",

  control.compute = list(dic = TRUE, waic = TRUE),

  control.predictor = list(A = inla.stack.A(stack_data), compute = TRUE)

)

# 总结和可视化

summary(fit)

ggregplot::ggField(fit, Mesh, Groups = 1) + 

  scale_fill_brewer(palette = "Blues")