时空模型2

 

```{r}

library(INLA)

library(inlabru)

library(magrittr)

library(sf)

library(tidyverse)

sf_data <- readRDS('data/spatial_data.rds')

# 基于最早日期计算连续周数

st_covid <- readRDS('data/st_covid.rds') |>

  mutate(week = as.numeric(difftime(date, min(date, na.rm = TRUE), units = "weeks")) %/% 1 + 1) |>

  left_join(sf_data@data |> dplyr::select(MSOA11CD,LONG, LAT), by = "MSOA11CD") |>

  drop_na(cases, Population, IMD, perc.wb, perc.ba, age1, pm25) |>

  st_as_sf(coords=c("LONG", "LAT"))

```

```{r}

coordinates <- st_coordinates(st_covid)  # st_covid geometry 类型是 sfc_POINT

locs <- unique(coordinates)

# # 计算每个多边形的质心

# centroids <- st_centroid(sf_data)

# locs <- st_coordinates(centroids)

# 计算初始参数

initial_range <- diff(range(locs[, 1])) / 3

max_edge <- initial_range / 2

# 构建网格

mesh <- fm_mesh_2d_inla(

  st_covid,

  max.edge = c(1, 2) * max_edge,

  offset = c(initial_range, initial_range),

  cutoff = max_edge / 5

)

plot(mesh)

```

```{r}

# 3. 定义SPDE模型

spde <- inla.spde2.pcmatern(

  mesh = mesh,

  prior.range = c(initial_range, 0.5),

  prior.sigma = c(1, 0.01)

)

# 4. 设置时间趋势先验，

rw1_prior <- list(theta = list(prior = "pc.prec", param = c(1, 0.01)))

```

```{r}

form <- cases ~ 1 +

    offset(log(Population)) +

    IMD + perc.wb + perc.ba + age1 + pm25 +

    site(st_coordinates(geometry), model = spde) +

    trend(main = week, model = 'rw1', hyper = rw1_prior)

model1 <- bru(form,data = st_covid,family = "poisson",

    E = st_covid$Population,

    control.family = list(link = "log"),

    options = list(

    control.inla = list(

      reordering = "metis",

      int.strategy = "eb"

    ),

    verbose = TRUE,

    inla.mode = "experimental"

  ))

summary(model1)

```

```{r}

# 创建时空交互模型

form2 <- cases ~ 1 +

  offset(log(Population)) +

  IMD + perc.wb + perc.ba + age1 + pm25 +

  f(st_coordinates(geometry),model = spde, group = week,

              control.group = list(model = 'rw1', hyper = rw1_prior))

model2 <- bru(form2,data = st_covid,family = "poisson",

    E = st_covid$Population,

    control.family = list(link = "log"),

    options = list(

    control.inla = list(

      reordering = "metis",

      int.strategy = "eb"

    ),

    verbose = TRUE,

    inla.mode = "experimental"

  ))

summary(model2)

```

```{r}

saveRDS(model1,'/mnt/d/inlabru/data/model1_stf.rds')

saveRDS(model2,'/mnt/d/inlabru/data/model1_sth.rds')

```