时空数据建模

 library(bmstdr)

m <- getData(name = "GADM", country = "Spain", level = 0) %>%

  st_as_sf() %>%

  st_cast("POLYGON") %>%

  mutate(area = st_area(.)) %>%

  arrange(desc(area)) %>%

  slice(1) %>%  st_transform(25830)

d <- read.csv("./dataPM25.csv") %>% 

  dplyr::select(year='ReportingYear',id='StationLocalId',

                long='SamplingPoint_Longitude',lat='SamplingPoint_Latitude',value='AQValue') %>% 

  st_as_sf(coords=c('long','lat')) %>% 

  st_set_crs(4326) %>% 

  st_transform(25830) %>% 

  filter(row_number() %in% st_intersects(m,.)[[1]]) %>% 

  group_by(id) %>% 

  filter(n()==3) %>% 

  ungroup() %>% 

  mutate(yr = year - min(year) + 1,

         s.index = as.integer(as.factor(id))) %>%

  arrange(s.index, yr) %>% 

  mutate(Longitude = st_coordinates(.)[,1], Latitude = st_coordinates(.)[,2])

N <- 45

burn.in <- 5

n.report <- 2

ggplot(m) + 

  geom_sf() + 

  coord_sf(datum = NA) +

  geom_sf(

    data = d, aes(color = value),

    size = 2

  ) +

  labs(x = "", y = "") +

  scale_color_viridis() +

  facet_wrap(~year) +

  theme_bw()

# scale.transform 响应变量的转换。它可以采用三个值：SQRT、LOG 或 NONE。默认值为“NONE”。

# mchoice - 如果通过输入参数 mchoice=TRUE 请求了计算模型选择统计信息

# plotit逻辑标量值：是否根据观测值绘制预测值。

m1 <- Bsptime(model = "AR",formula=value~yr,data=d,coordtype = 'utm',coords=7:8,package = 'inla', 

             mchoice = TRUE,n.report = n.report, N = N, burn.in = burn.in,plotit=T)

summary(m1)

names(m1)

m1$fit$summary.fixed

m1$fit$summary.random

m1$fit$summary.hyperpar

m1$fit$summary.fitted.values

list_marginals <- list(

  "b0" = m1$fit$marginals.fixed$X....1,

  "precision Gaussian obs" =

    m1$fit$marginals.hyperpar$"Precision for the Gaussian observations",

  "range" = m1$fit$marginals.hyperpar$"Range for spatial.field",

  "stdev" = m1$fit$marginals.hyperpar$"Stdev for spatial.field",

  "rho" = m1$fit$marginals.hyperpar$"GroupRho for spatial.field"

)

marginals <- data.frame(do.call(rbind, list_marginals))

marginals$parameter <- rep(names(list_marginals),

                           times = sapply(list_marginals, nrow)

)

ggplot(marginals, aes(x = x, y = y)) + 

  geom_line() +

  facet_wrap(~parameter, scales = "free") +

  labs(x = "", y = "Density") + theme_bw()