10.3 Modeling!

\(Y_{it}\) : \(PM_{2.5}\) \(i\) are the locations at \(t\) times.

\[Y_{it} \sim N(\mu_{it}, \sigma^2_e )\]

\[\mu_{it} = \beta_0 + \xi(\boldsymbol{x}_i, t)\]

\(\sigma_2²\) = variance of the measurement error (not correlated in space and time)

\(\beta_0\) = intercept

\(\xi(\boldsymbol{x}_i, t)\): random effect changes in time first order autoregressive dynamics and spatially correlated innovations

\[\xi(\boldsymbol{x}_i, t) = a \xi(\boldsymbol{x}_i, t-1) + w(\boldsymbol{x}_i, t) \]

First part is first order autoregressive process (think \(Y_t = \beta_0 + \beta_1 Y_{t-1} + u_t\)) and second part is the classic :

follows a zero-mean Gaussian distribution temporally independent but spatially dependent at each time with Matérn covariance function.

This is a “simple”: we could imagine relation between time and space and had more covariate (temperature, precipitation, etc …)

10.3.1 Mesh construction

This is like the previous chapter:

coo <- cbind(d$x, d$y)
bnd <- inla.nonconvex.hull(st_coordinates(m)[, 1:2])
mesh <- inla.mesh.2d(
  loc = coo, boundary = bnd,
  max.edge = c(100000, 200000), # huge!
  cutoff = 1000
)

plot(mesh)
points(coo, col = "red")

10.3.2 SPDE model on the mesh with penalised complexity

Picking prior for range and sigma:

spde <- inla.spde2.pcmatern(
  mesh = mesh, alpha = 2, constr = TRUE, # same as last week
  prior.range = c(10000, 0.01), # P(range < 10000) = 0.01 / P(range < v_r) = p_r
  prior.sigma = c(3, 0.01) # P(sigma > 3) = 0.01 / P(sigma > v_s) = p_s
)

10.3.3 Index

timesn <- length(unique(d$year))
indexs <- inla.spde.make.index("s",
  n.spde = spde$n.spde,
  n.group = timesn
)
lengths(indexs) # pay attention to lengthS()

##       s s.group  s.repl 
##    2142    2142    2142

10.3.4 projection matrix:

group <- d$year - min(d$year) + 1 # we are just creating group similar to index
A <- inla.spde.make.A(mesh = mesh, loc = coo, 
                      group = group) # new stuff
# 299 x 2142 sparse Matrix of class "dgCMatrix"

10.3.5 Grid:

bb <- st_bbox(m)
x <- seq(bb$xmin - 1, bb$xmax + 1, length.out = 50) # do not know if i should convert first or not ?
y <- seq(bb$ymin - 1, bb$ymax + 1, length.out = 50)
dp <- as.matrix(expand.grid(x, y)) 

p <- st_as_sf(data.frame(x = dp[, 1], y = dp[, 2]),
  coords = c("x", "y")
) # this could be a function no ?

st_crs(p) <- st_crs(25830)
ind <- st_intersects(m, p)
dp <- dp[ind[[1]], ]
plot(dp, asp = 1)

We just repeat it 3 times (3 years of data)

dp <- rbind(cbind(dp, 1), cbind(dp, 2), cbind(dp, 3)) # 1,2,3 are recycled
dim(dp)

## [1] 3930    3

10.3.6 Ap + stack

coop <- dp[, 1:2]
groupp <- dp[, 3]
Ap <- inla.spde.make.A(mesh = mesh, loc = coop, group = groupp)

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

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

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

10.3.7 Model formula:

This will take a bit of time (around 1 min)!

rprior <- list(theta = list(prior = "pccor1", param = c(0, 0.9)))

formula <- y ~ 0 + b0 + f(s, #space time
  model = spde, 
  group = s.group, # years 
  control.group = list(model = "ar1", # auto regresive first order
                       hyper = rprior)
)


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

I have not reproduce the summary of res but we have a range around 17857 m (from 14510 to 21721).

list_marginals <- list(
"b0" = res$marginals.fixed$b0,
"precision Gaussian obs" =
res$marginals.hyperpar$"Precision for the Gaussian observations",
"range" = res$marginals.hyperpar$"Range for s", # yes this is how they are called inside res
"stdev" = res$marginals.hyperpar$"Stdev for s",
"rho" = res$marginals.hyperpar$"GroupRho for s"
)


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

If I am correct rho (and GroupRho) are correlation between group (ie years)