8.3 Spatial modeling of rainfall in Paraná, Brazil

8.3.1 Data

Average rainfall recorded over different year over May-June at 143 stations

data.frame(cbind(parana$coords, Rainfall = parana$data)) |>
ggplot()+
  geom_point(aes(east, north, color = Rainfall), size = 2) +
  coord_fixed(ratio = 1) +
  scale_color_gradient(low = "blue", high = "orange") +
  geom_path(data = data.frame(parana$border), aes(east, north)) +
  theme_bw()

8.3.2 Model

\[ Y_i \sim N (\mu_i, \sigma^2), i = 1,2, ..., n\]

\[\mu_i = \beta_0 + Z(si) \]

Rainfall as a normal distribution with a mean that follow a zero-mean Gaussian process with a Matèrn covariance function.

8.3.2.1 Mesh construction

coo <- parana$coords # I think they are in degree
summary(dist(coo))

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     1.0   143.5   231.2   244.1   333.6   619.5

mesh <- inla.mesh.2d(
  loc = coo, # our coords
  offset = c(50, 100), # inner border and outer border
  cutoff = 1, # avoid building small and plenty triangle if point are close
  max.edge = c(30, 60) # small triangle (30) and big one outside (60)
)

# we have a warning of a deprecated function
plot(mesh)
points(coo, col = "red")

str(mesh, max.level = 1)

## List of 8
##  $ meta    :List of 5
##  $ manifold: chr "R2"
##  $ n       : int 1189
##  $ loc     : num [1:1189, 1:3] 219 679 819 819 614 ...
##  $ graph   :List of 5
##  $ segm    :List of 2
##  $ idx     :List of 2
##  $ crs     : NULL
##  - attr(*, "class")= chr "inla.mesh"

On this one we did not use boundary:

bnd <- inla.nonconvex.hull(coo)
meshb <- inla.mesh.2d(
  boundary = bnd, offset = c(50, 100),
  cutoff = 1, max.edge = c(30, 60)
)
plot(meshb)
points(coo, col = "red")

8.3.2.2 SPDE model on the mesh

\[ \alpha = v + d/2\]

d = 2 in spatial 2D case

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

indexs <- inla.spde.make.index("s", spde$n.spde)

8.3.2.3 Projection matrix:

We need to build a matrix that contains n rows (n = number of locations) and m column (m = vertex of triangle).

A <- inla.spde.make.A(mesh = mesh, loc = coo)
rowSums(A) |> sum()

## [1] 143

8.3.2.4 Prediction data

We nee to build a grid: 50 x 50

bb <- bbox(parana$border) # bounding box
x <- seq(bb[1, "min"] - 1, bb[1, "max"] + 1, length.out = 50) # increase a bit the size 
y <- seq(bb[2, "min"] - 1, bb[2, "max"] + 1, length.out = 50)
coop <- as.matrix(expand.grid(x, y))

Then you keep point only inside the border:

ind <- point.in.polygon(
  coop[, 1], coop[, 2],
  parana$border[, 1], parana$border[, 2]
)
coop <- coop[which(ind == 1), ]
plot(coop, asp = 1)

And create a projection matrix with it:

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

## [1] 1533 1189

8.3.2.5 Building estimation and prediction “stack”

# stack for estimation stk.e
stk.e <- inla.stack(
  tag = "est",
  data = list(y = parana$data),
  A = list(1, A),
  effects = list(data.frame(b0 = rep(1, nrow(coo))), s = indexs) # our models b0 and spatial effect
)

# stack for prediction stk.p
stk.p <- inla.stack(
  tag = "pred",
  data = list(y = NA), # to be filled
  A = list(1, Ap),
  effects = list(data.frame(b0 = rep(1, nrow(coop))), s = indexs)
)

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

8.3.2.6 Formula

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

8.3.2.7 Inla() call

This is not intuitive:

res <- inla(formula,
  data = inla.stack.data(stk.full), #it has both data and grid
  control.predictor = list(
    compute = TRUE,
    A = inla.stack.A(stk.full)
  )
)

8.3.2.8 result!

We will get the index of the rows with prediction and use it to extract the correct mean and quantiles:

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

We build a tidy data frame with it:

dpm <- rbind(
  data.frame(
    east = coop[, 1], north = coop[, 2],
    value = pred_mean, variable = "pred_mean"
  ),
  data.frame(
    east = coop[, 1], north = coop[, 2],
    value = pred_ll, variable = "pred_ll"
  ),
  data.frame(
    east = coop[, 1], north = coop[, 2],
    value = pred_ul, variable = "pred_ul"
  )
)
dpm$variable <- as.factor(dpm$variable)

And map it:

ggplot(dpm) + geom_tile(aes(east, north, fill = value)) +
  facet_wrap(~variable, nrow = 1) +
  coord_fixed(ratio = 1) +
  scale_fill_gradient(
    name = "Rainfall",
    low = "blue", high = "orange"
  ) +
  theme_bw()

8.3.2.9 just the spatial field:

newloc <- cbind(c(219, 678, 818), c(20, 20, 160))
Aproj <- inla.spde.make.A(mesh, loc = newloc)
Aproj %*% res$summary.random$s$mean

## 3 x 1 Matrix of class "dgeMatrix"
##            [,1]
## [1,] -3.1661566
## [2,]  0.8556036
## [3,] -1.2270720

rang <- apply(mesh$loc[, c(1, 2)], 2, range)
proj <- inla.mesh.projector(mesh,
  xlim = rang[, 1], ylim = rang[, 2],
  dims = c(300, 300)
)

mean_s <- inla.mesh.project(proj, res$summary.random$s$mean)
sd_s <- inla.mesh.project(proj, res$summary.random$s$sd)

df <- expand.grid(x = proj$x, y = proj$y)
df$mean_s <- as.vector(mean_s)
df$sd_s <- as.vector(sd_s)

library(viridis)
library(cowplot)

gmean <- ggplot(df, aes(x = x, y = y, fill = mean_s)) +
  geom_raster() +
  scale_fill_viridis(na.value = "transparent") +
  coord_fixed(ratio = 1) + theme_bw()

gsd <- ggplot(df, aes(x = x, y = y, fill = sd_s)) +
  geom_raster() +
  scale_fill_viridis(na.value = "transparent") +
  coord_fixed(ratio = 1) + theme_bw()

plot_grid(gmean, gsd)