PeMS 2, modeling

Go back to the About page.

This vignette compares different models for PeMS data. It uses pems_repl1_data.RData, which is a file with a graph and data created in pems_repl1.html.

Let us set some global options for all code chunks in this document.

# Set seed for reproducibility
set.seed(1938) 
# Set global options for all code chunks
knitr::opts_chunk$set(
  # Disable messages printed by R code chunks
  message = FALSE,    
  # Disable warnings printed by R code chunks
  warning = FALSE,    
  # Show R code within code chunks in output
  echo = TRUE,        
  # Include both R code and its results in output
  include = TRUE,     
  # Evaluate R code chunks
  eval = TRUE,       
  # Enable caching of R code chunks for faster rendering
  cache = FALSE,      
  # Align figures in the center of the output
  fig.align = "center",
  # Enable retina display for high-resolution figures
  retina = 2,
  # Show errors in the output instead of stopping rendering
  error = TRUE,
  # Do not collapse code and output into a single block
  collapse = FALSE
)

Below we load the necessary libraries.

library(INLA)
library(inlabru)
library(rSPDE)
library(MetricGraph)

library(dplyr)
library(plotly)
library(scales)
library(patchwork)

library(ggplot2)
library(cowplot)
library(ggpubr) #annotate_figure()
library(grid) #textGrob()
library(ggmap)

library(viridis)
library(OpenStreetMap)


library(tidyr)
library(sf)

library(here)
library(rmarkdown)
library(grateful) # Cite all loaded packages

Below we define the function captioner() to generate captions for the figures and the function process_model_results() to extract the summary of the parameters of the model.

---

Press the Show button below to reveal the code.

---

# Start the figure counter
fig_count <- 0
# Define the captioner function
captioner <- function(caption) {
  fig_count <<- fig_count + 1
  paste0("Figure ", fig_count, ": ", caption)
}
process_model_results <- function(fit, model) {
  fit_spde <- rspde.result(fit, "field", model, parameterization = "spde")
  fit_matern <- rspde.result(fit, "field", model, parameterization = "matern")
  df_for_plot_spde <- gg_df(fit_spde)
  df_for_plot_matern <- gg_df(fit_matern)
  param_spde <- summary(fit_spde)
  param_matern <- summary(fit_matern)
  param_fixed <- fit$summary.fixed[,1:6]
  marginal.posterior.sigma_e = inla.tmarginal(
    fun = function(x) exp(-x/2), 
    marginal = fit[["internal.marginals.hyperpar"]][["Log precision for the Gaussian observations"]])
  quant.sigma_e <- capture.output({result_tmp <- inla.zmarginal(marginal.posterior.sigma_e)}, file = "/dev/null") 
  quant.sigma_e <- result_tmp
  statistics.sigma_e <- unlist(quant.sigma_e)[c(1,2,3,5,7)]
  mode.sigma_e <- inla.mmarginal(marginal.posterior.sigma_e)
  allparams <- rbind(param_fixed, param_spde, param_matern, c(statistics.sigma_e, mode.sigma_e))
  rownames(allparams)[nrow(allparams)] <- "sigma_e"
  return(list(allparams = allparams, df_for_plot_spde = df_for_plot_spde, df_for_plot_matern = df_for_plot_matern))
}

We first load the data in the file pems_repl1_data.RData and extract the data from the graph.

# Load the data
load(here("data_files/pems_repl1_data.RData"))
# Extract the data from the graph
data <- graph$get_data()

Below we extract the locations to compute the distance matrix. Using this matrix, we define the groups for cross-validation. Observe that we only compute the distance matrix for the first replicate and compute the groups for it. As all replicates share the same locations, we can use the groups structure from the first replicate for all replicates.

---

Press the Show button below to reveal the code.

---

# Define aux data frame to compute the distance matrix
aux <- data |> filter(repl == 1) |>
  rename(distance_on_edge = .distance_on_edge, edge_number = .edge_number) |> # Rename the variables (because graph$compute_geodist_PtE() requires so)
  as.data.frame() |> # Transform to a data frame (i.e., remove the metric_graph class)
  dplyr::select(edge_number, distance_on_edge)

# Compute the distance matrix
distmatrix <- graph$compute_geodist_PtE(PtE = aux,
                                             normalized = TRUE,
                                             include_vertices = FALSE)
# Define the distance vector
distance = seq(from = 0, to = 10, by = 0.1)
# Compute the groups for one replicate
GROUPS <- list()
for (j in 1:length(distance)) {
  GROUPS[[j]] = list()
  for (i in 1:nrow(aux)) {
    GROUPS[[j]][[i]] <- which(as.vector(distmatrix[i, ]) <= distance[j])
  }
}
# Compute the groups for all replicates, based on the groups of the first replicate
nrowY <- length(unique(data$repl))
ncolY <- nrow(filter(data, repl == 1))
NEW_GROUPS <- list()
for (j in 1:length(distance)) {
  my_list <- GROUPS[[j]]
  aux_list <- list()
  for (i in 0:(nrowY - 1)) {
  added_vectors <- lapply(my_list, function(vec) vec + i*ncolY)
  aux_list <- c(aux_list, added_vectors)
  }
  NEW_GROUPS[[j]] <- aux_list
}

GROUPS <- NEW_GROUPS

Below we plot to check that the groups are correctly defined.

point_of_interest <- 3 # Any number between 1 and nrow(data)
small_neighborhood <- GROUPS[[20]][[point_of_interest]]
large_neighborhood <- GROUPS[[50]][[point_of_interest]]
p <- graph$plot(vertex_size = 0) +
  geom_point(data = data, aes(x = .coord_x, y = .coord_y), color = "darkviolet", size = 2) +
  geom_point(data = data[large_neighborhood, ], aes(x = .coord_x, y = .coord_y), color = "green", size = 1.5) +
  geom_point(data = data[small_neighborhood, ], aes(x = .coord_x, y = .coord_y), color = "blue", size = 1) +
  geom_point(data = data[point_of_interest, ], aes(x = .coord_x, y = .coord_y), color = "red", size = 0.5) +
  ggtitle("Groups") + 
  theme_minimal() + 
  theme(text = element_text(family = "Palatino")) +
  coord_fixed()
ggplotly(p)

Figure 1: Illustrations of groups for cross-validation based on the distance matrix.

Below we define the non-stationary parameters.

# Non-stationary parameters
B.tau = cbind(0, 1, 0, cov, 0)
B.kappa = cbind(0, 0, 1, 0, cov)

We now model the speed records \(y_i\) as 13 independent replicates satisfying \[\begin{equation} \label{applimodel} y_i|u(\cdot)\sim N(\beta_0 + \beta_1\text{mean.cov}(s_i) + u(s_i),\sigma_\epsilon^2),\;i = 1,\dots, 314, \end{equation}\] where \(u(\cdot)\) is a Gaussian process on the highway network. We consider stationary models with \(\kappa,\tau>0\) and non-stationary models where \(\tau\) and \(\kappa\) are given by \[\begin{equation} \label{logregressions} \begin{aligned} \log(\tau(s)) &= \theta_1 + \theta_3 \text{std.cov}(s),\\ \log(\kappa(s)) &= \theta_2 + \theta_4 \text{std.cov}(s). \end{aligned} \end{equation}\]

For each of the two classes of models, we consider three cases: when (1) \(\nu\) is fixed to 0.5 or (2) 1.5, and (3) \(\nu\) is estimated from the data.

Below cov refers to \(\text{std.cov}(s)\) and mean_value refers to \(\text{mean.cov}(s)\).

Case \(\nu = 0.5\)

We first consider the stationary model.

---

Press the Show button below to reveal the code.

---

# Build the model
rspde_model_stat <- rspde.metric_graph(graph,
                                       parameterization = "spde",
                                       nu = 0.5)
# Prepare the data for fitting
data_rspde_bru_stat <- graph_data_rspde(rspde_model_stat,
                                        repl = ".all",
                                        bru = TRUE,
                                        repl_col = "repl")
# Define the component
cmp_stat <- y ~ -1 +
  Intercept(1) +
  mean_value +
  field(cbind(.edge_number, .distance_on_edge), 
        model = rspde_model_stat,
        replicate = repl)
# Fit the model
rspde_fit_stat <-
  bru(cmp_stat,
      data = data_rspde_bru_stat[["data"]],
      family = "gaussian",
      options = list(verbose = FALSE)
  )

output_from_models <-process_model_results(rspde_fit_stat, rspde_model_stat)
parameters_statistics <- output_from_models$allparams
mean_and_mode_params_statnu0.5 <- parameters_statistics[, c(1,6)]
rspde_fit_statnu0.5 <- rspde_fit_stat
# Summarize the results
summary(rspde_fit_stat)

## inlabru version: 2.12.0.9002
## INLA version: 24.12.11
## Components:
## Intercept: main = linear(1), group = exchangeable(1L), replicate = iid(1L), NULL
## mean_value: main = linear(mean_value), group = exchangeable(1L), replicate = iid(1L), NULL
## field: main = cgeneric(cbind(.edge_number, .distance_on_edge)), group = exchangeable(1L), replicate = iid(repl), NULL
## Likelihoods:
##   Family: 'gaussian'
##     Tag: ''
##     Data class: 'metric_graph_data', 'data.frame'
##     Response class: 'numeric'
##     Predictor: y ~ .
##     Used components: effects[Intercept, mean_value, field], latent[]
## Time used:
##     Pre = 0.221, Running = 12.9, Post = 5.87, Total = 19 
## Fixed effects:
##               mean    sd 0.025quant 0.5quant 0.975quant    mode kld
## Intercept  -14.463 1.556    -17.532  -14.461    -11.403 -14.460   0
## mean_value   1.276 0.021      1.236    1.276      1.317   1.276   0
## 
## Random effects:
##   Name     Model
##     field CGeneric
## 
## Model hyperparameters:
##                                           mean    sd 0.025quant 0.5quant
## Precision for the Gaussian observations  0.014 0.000      0.013    0.014
## Theta1 for field                        -1.497 0.048     -1.590   -1.498
## Theta2 for field                        -2.983 0.197     -3.398   -2.974
##                                         0.975quant   mode
## Precision for the Gaussian observations      0.014  0.014
## Theta1 for field                            -1.399 -1.503
## Theta2 for field                            -2.625 -2.931
## 
## Deviance Information Criterion (DIC) ...............: 29843.76
## Deviance Information Criterion (DIC, saturated) ....: 4794.42
## Effective number of parameters .....................: 717.50
## 
## Watanabe-Akaike information criterion (WAIC) ...: 29813.53
## Effective number of parameters .................: 602.21
## 
## Marginal log-Likelihood:  -15148.04 
##  is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')

parameters_statistics

---

Press the Show button below to reveal the code.

---

# Plot the estimates of the parameters
ggplot(output_from_models$df_for_plot_spde) + geom_line(aes(x = x, y = y)) + 
  facet_wrap(~parameter, scales = "free") + labs(y = "Density") +
  theme(text = element_text(family = "Palatino"))

Figure 2: Posterior distributions of the spde parameters.

---

Press the Show button below to reveal the code.

---

ggplot(output_from_models$df_for_plot_matern) + geom_line(aes(x = x, y = y)) + 
  facet_wrap(~parameter, scales = "free") + labs(y = "Density") + 
  theme(text = element_text(family = "Palatino"))

Figure 3: Posterior distributions of the matern parameters.

We now fit the non-stationary model.

---

Press the Show button below to reveal the code.

---

# Build the model
rspde_model_nonstat <- rspde.metric_graph(graph,
                                          B.tau = B.tau,
                                          B.kappa =  B.kappa,
                                          parameterization = "spde",
                                          nu = 0.5)
# Prepare the data for fitting
data_rspde_bru_nonstat <- graph_data_rspde(rspde_model_nonstat,
                                        repl = ".all",
                                        bru = TRUE,
                                        repl_col = "repl")
# Define the component
cmp_nonstat <- y ~ -1 +
  Intercept(1) +
  mean_value +
  field(cbind(.edge_number, .distance_on_edge), 
        model = rspde_model_nonstat,
        replicate = repl)
# Fit the model
rspde_fit_nonstat <-
  bru(cmp_nonstat,
      data = data_rspde_bru_nonstat[["data"]],
      family = "gaussian",
      options = list(verbose = FALSE)
  )

output_from_models <- process_model_results(rspde_fit_nonstat, rspde_model_nonstat)
parameters_statistics <- output_from_models$allparams
mean_and_mode_params_nonstatnu0.5 <- parameters_statistics[, c(1,6)]
rspde_fit_nonstatnu0.5 <- rspde_fit_nonstat
# Summarize the results
summary(rspde_fit_nonstat)

## inlabru version: 2.12.0.9002
## INLA version: 24.12.11
## Components:
## Intercept: main = linear(1), group = exchangeable(1L), replicate = iid(1L), NULL
## mean_value: main = linear(mean_value), group = exchangeable(1L), replicate = iid(1L), NULL
## field: main = cgeneric(cbind(.edge_number, .distance_on_edge)), group = exchangeable(1L), replicate = iid(repl), NULL
## Likelihoods:
##   Family: 'gaussian'
##     Tag: ''
##     Data class: 'metric_graph_data', 'data.frame'
##     Response class: 'numeric'
##     Predictor: y ~ .
##     Used components: effects[Intercept, mean_value, field], latent[]
## Time used:
##     Pre = 0.173, Running = 32.5, Post = 5.28, Total = 38 
## Fixed effects:
##              mean    sd 0.025quant 0.5quant 0.975quant   mode kld
## Intercept  -4.528 1.994     -8.425   -4.533     -0.604 -4.532   0
## mean_value  1.126 0.030      1.068    1.126      1.184  1.126   0
## 
## Random effects:
##   Name     Model
##     field CGeneric
## 
## Model hyperparameters:
##                                           mean    sd 0.025quant 0.5quant
## Precision for the Gaussian observations  0.015 0.000      0.014    0.015
## Theta1 for field                        -1.184 0.062     -1.303   -1.185
## Theta2 for field                        -6.287 2.100    -11.243   -5.956
## Theta3 for field                        -0.614 0.035     -0.686   -0.613
## Theta4 for field                         3.002 1.675      0.696    2.740
##                                         0.975quant   mode
## Precision for the Gaussian observations      0.016  0.015
## Theta1 for field                            -1.061 -1.187
## Theta2 for field                            -3.401 -4.458
## Theta3 for field                            -0.547 -0.610
## Theta4 for field                             6.953  1.556
## 
## Deviance Information Criterion (DIC) ...............: 29452.98
## Deviance Information Criterion (DIC, saturated) ....: 4717.64
## Effective number of parameters .....................: 638.98
## 
## Watanabe-Akaike information criterion (WAIC) ...: 29491.49
## Effective number of parameters .................: 580.43
## 
## Marginal log-Likelihood:  -14963.96 
##  is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')

parameters_statistics

---

Press the Show button below to reveal the code.

---

# Plot the estimates of the parameters
ggplot(output_from_models$df_for_plot_spde) + geom_line(aes(x = x, y = y)) + 
  facet_wrap(~parameter, scales = "free") + labs(y = "Density") + 
  theme(text = element_text(family = "Palatino"))

Figure 4: Posterior distributions of the spde parameters.

---

Press the Show button below to reveal the code.

---

ggplot(output_from_models$df_for_plot_matern) + geom_line(aes(x = x, y = y)) + 
  facet_wrap(~parameter, scales = "free") + labs(y = "Density") + 
  theme(text = element_text(family = "Palatino"))

Figure 5: Posterior distributions of the matern parameters.

Below we consider the prediction of replicate 14.

---

Press the Show button below to reveal the code.

---

# Load the maps p12 and p13 from pems_repl1 vignette
load(here("data_files/maps_zoom12and13from_stadia.RData"))
# We consider replicate 14
replicate.number <- 1
# Prepare the data for prediction
data_prd_list_for_rep <- data_prd_list_mesh
data_prd_list_for_rep[["mean_value"]] <- cov_for_mean_to_plot
data_prd_list_for_rep[["repl"]] <- rep(replicate.number, nrow(data_prd_list_mesh))
# Perform the prediction
repl1_pred_full <- predict(rspde_fit_nonstat, newdata = data_prd_list_for_rep, ~Intercept + mean_value + field_eval(cbind(.edge_number, .distance_on_edge), replicate = repl))
repl1_pred_mean <- repl1_pred_full$mean
# Extract the Euclidean coordinates of the mesh points
xypoints <- graph$mesh$V
# Extract the range of the coordinates 
x_left <- range(xypoints[,1])[1]
x_right <- range(xypoints[,1])[2]
y_bottom <- range(xypoints[,2])[1]
y_top <- range(xypoints[,2])[2]
# Define coordinates for small windows
coordx_lwr1 <- -121.878
coordx_upr1 <- -121.828
coordy_lwr1 <- 37.315
coordy_upr1 <- 37.365

coordx_lwr2<- -122.075
coordx_upr2 <- -122.025
coordy_lwr2 <- 37.365
coordy_upr2 <- 37.415
# Define the colors for the windows
lower_color <- "darkred"   # Dark purple
upper_color <- "darkblue"  # Yellow
# Plot the field on top of the map
f12 <- graph$plot_function(X = repl1_pred_mean, 
                          vertex_size = 0, 
                          p = p12,
                          edge_width = 0.5) + 
  theme_minimal() + 
  theme(text = element_text(family = "Palatino"), 
        axis.text = element_text(size = 8),
        legend.text = element_text(size = 8),
        plot.margin = unit(-0.4*c(1,0,1,1), "cm")
        ) +
  labs(color = "", x = "", y = "") +
  xlim(x_left, x_right) + 
  ylim(y_bottom, y_top)
# Plot the field on top of the map
f13 <- graph$plot_function(X = repl1_pred_mean, 
                          vertex_size = 0, 
                          p = p13,
                          edge_width = 0.5) + 
  theme_minimal() + 
  theme(text = element_text(family = "Palatino"), 
        axis.text = element_text(size = 8),
        legend.text = element_text(size = 8),
        plot.margin = unit(-0.4*c(1,0,1,1), "cm")
        ) +
  labs(color = "", x = "", y = "") +
  xlim(x_left, x_right) + 
  ylim(y_bottom, y_top)

g12 <- graph$plot(data = "y", group = 1, vertex_size = 0, p = f12, edge_width = 0, data_size = 1) + 
  labs(color = "", x = "", y = "") + 
  annotate("segment", x = coordx_lwr1, y = coordy_lwr1, xend = coordx_upr1, yend = coordy_lwr1, 
           linewidth = 0.4, color = upper_color) +  # Bottom line
  annotate("segment", x = coordx_lwr1, y = coordy_upr1, xend = coordx_upr1, yend = coordy_upr1, 
           linewidth = 0.4, color = upper_color) +  # Top line
  annotate("segment", x = coordx_lwr1, y = coordy_lwr1, xend = coordx_lwr1, yend = coordy_upr1, 
           linewidth = 0.4, color = upper_color) +  # Left line
  annotate("segment", x = coordx_upr1, y = coordy_lwr1, xend = coordx_upr1, yend = coordy_upr1, 
           linewidth = 0.4, color = upper_color) +  # Right line
  annotate("segment", x = coordx_lwr2, y = coordy_lwr2, xend = coordx_upr2, yend = coordy_lwr2, 
           linewidth = 0.4, color = lower_color) +  
  annotate("segment", x = coordx_lwr2, y = coordy_upr2, xend = coordx_upr2, yend = coordy_upr2, 
           linewidth = 0.4, color = lower_color) +  
  annotate("segment", x = coordx_lwr2, y = coordy_lwr2, xend = coordx_lwr2, yend = coordy_upr2, 
           linewidth = 0.4, color = lower_color) +  
  annotate("segment", x = coordx_upr2, y = coordy_lwr2, xend = coordx_upr2, yend = coordy_upr2, 
           linewidth = 0.4, color = lower_color)
g13 <- graph$plot(data = "y", group = 1, vertex_size = 0, p = f13, edge_width = 0, data_size = 1) + 
  labs(color = "", x = "", y = "") + 
  annotate("segment", x = coordx_lwr1, y = coordy_lwr1, xend = coordx_upr1, yend = coordy_lwr1, 
           linewidth = 0.4, color = upper_color) +  # Bottom line
  annotate("segment", x = coordx_lwr1, y = coordy_upr1, xend = coordx_upr1, yend = coordy_upr1, 
           linewidth = 0.4, color = upper_color) +  # Top line
  annotate("segment", x = coordx_lwr1, y = coordy_lwr1, xend = coordx_lwr1, yend = coordy_upr1, 
           linewidth = 0.4, color = upper_color) +  # Left line
  annotate("segment", x = coordx_upr1, y = coordy_lwr1, xend = coordx_upr1, yend = coordy_upr1, 
           linewidth = 0.4, color = upper_color) +  # Right line
  annotate("segment", x = coordx_lwr2, y = coordy_lwr2, xend = coordx_upr2, yend = coordy_lwr2, 
           linewidth = 0.4, color = lower_color) +  
  annotate("segment", x = coordx_lwr2, y = coordy_upr2, xend = coordx_upr2, yend = coordy_upr2, 
           linewidth = 0.4, color = lower_color) +  
  annotate("segment", x = coordx_lwr2, y = coordy_lwr2, xend = coordx_lwr2, yend = coordy_upr2, 
           linewidth = 0.4, color = lower_color) +  
  annotate("segment", x = coordx_upr2, y = coordy_lwr2, xend = coordx_upr2, yend = coordy_upr2, 
           linewidth = 0.4, color = lower_color) 

r1 <- g13 + xlim(coordx_lwr1, coordx_upr1) + 
                    ylim(coordy_lwr1, coordy_upr1) + 
  theme(legend.position = "none", 
        plot.margin = unit(-0.2*c(1,1,1,1), "cm"))

r2 <- g13 + xlim(coordx_lwr2, coordx_upr2) + 
                    ylim(coordy_lwr2, coordy_upr2) + 
  theme(legend.position = "none", 
        plot.margin = unit(-0.2*c(1,1,1,1), "cm"))

# Arrange p2 and p3 horizontally
left_col <- plot_grid(r2, r1, labels = NULL, ncol = 1, nrow = 2, rel_heights = c(1,1))

# Combine the top row with p1 in a grid
combined_plot <- plot_grid(left_col, g12, labels = NULL, ncol = 2, rel_widths = c(1,2)) 
final_plot <- annotate_figure(combined_plot, left = textGrob("Latitude", rot = 90, vjust = 1, gp = gpar(cex = 0.8)),
                              bottom = textGrob("Longitude", vjust = -0.5, gp = gpar(cex = 0.8)))
ggsave(here("data_files/replicate14_3_with_prediction.png"), width = 11.2, height = 5.43, plot = final_plot, dpi = 500)
# Print the combined plot
print(final_plot)

Figure 6: Speed observations (in mph) on the highway network of the city of San Jose in California, recorded on April 3, 2017. The left panels are zoomed-in areas of the panel to the right.

Case \(\nu = 1.5\)

We first consider the stationary model.

---

Press the Show button below to reveal the code.

---

# Build the model
rspde_model_stat <- rspde.metric_graph(graph,
                                       parameterization = "spde",
                                       nu = 1.5)
# Prepare the data for fitting
data_rspde_bru_stat <- graph_data_rspde(rspde_model_stat,
                                        repl = ".all",
                                        bru = TRUE,
                                        repl_col = "repl")
# Define the component
cmp_stat <- y ~ -1 +
  Intercept(1) +
  mean_value +
  field(cbind(.edge_number, .distance_on_edge), 
        model = rspde_model_stat,
        replicate = repl)
# Fit the model
rspde_fit_stat <-
  bru(cmp_stat,
      data = data_rspde_bru_stat[["data"]],
      family = "gaussian",
      options = list(verbose = FALSE)
  )

output_from_models <-process_model_results(rspde_fit_stat, rspde_model_stat)
parameters_statistics <- output_from_models$allparams
mean_and_mode_params_statnu1.5 <- parameters_statistics[, c(1,6)]
rspde_fit_statnu1.5 <- rspde_fit_stat
# Summarize the results
summary(rspde_fit_stat)

## inlabru version: 2.12.0.9002
## INLA version: 24.12.11
## Components:
## Intercept: main = linear(1), group = exchangeable(1L), replicate = iid(1L), NULL
## mean_value: main = linear(mean_value), group = exchangeable(1L), replicate = iid(1L), NULL
## field: main = cgeneric(cbind(.edge_number, .distance_on_edge)), group = exchangeable(1L), replicate = iid(repl), NULL
## Likelihoods:
##   Family: 'gaussian'
##     Tag: ''
##     Data class: 'metric_graph_data', 'data.frame'
##     Response class: 'numeric'
##     Predictor: y ~ .
##     Used components: effects[Intercept, mean_value, field], latent[]
## Time used:
##     Pre = 0.18, Running = 11.9, Post = 4.55, Total = 16.7 
## Fixed effects:
##               mean    sd 0.025quant 0.5quant 0.975quant    mode kld
## Intercept  -16.096 1.308    -18.686  -16.088    -13.555 -16.088   0
## mean_value   1.308 0.023      1.264    1.307      1.352   1.307   0
## 
## Random effects:
##   Name     Model
##     field CGeneric
## 
## Model hyperparameters:
##                                           mean    sd 0.025quant 0.5quant
## Precision for the Gaussian observations  0.013 0.000      0.012    0.013
## Theta1 for field                        -1.104 0.105     -1.306   -1.106
## Theta2 for field                        -1.370 0.084     -1.540   -1.368
##                                         0.975quant   mode
## Precision for the Gaussian observations      0.014  0.013
## Theta1 for field                            -0.893 -1.113
## Theta2 for field                            -1.207 -1.363
## 
## Deviance Information Criterion (DIC) ...............: 29832.43
## Deviance Information Criterion (DIC, saturated) ....: 4580.50
## Effective number of parameters .....................: 506.73
## 
## Watanabe-Akaike information criterion (WAIC) ...: 29789.61
## Effective number of parameters .................: 418.30
## 
## Marginal log-Likelihood:  -15154.67 
##  is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')

parameters_statistics

---

Press the Show button below to reveal the code.

---

# Plot the estimates of the parameters
ggplot(output_from_models$df_for_plot_spde) + geom_line(aes(x = x, y = y)) + 
  facet_wrap(~parameter, scales = "free") + labs(y = "Density") + 
  theme(text = element_text(family = "Palatino"))

Figure 7: Posterior distributions of the spde parameters.

---

Press the Show button below to reveal the code.

---

ggplot(output_from_models$df_for_plot_matern) + geom_line(aes(x = x, y = y)) + 
  facet_wrap(~parameter, scales = "free") + labs(y = "Density") + 
  theme(text = element_text(family = "Palatino"))

Figure 8: Posterior distributions of the matern parameters.

We now fit the non-stationary model.

---

Press the Show button below to reveal the code.

---

# Build the model
rspde_model_nonstat <- rspde.metric_graph(graph,
                                          B.tau = B.tau,
                                          B.kappa =  B.kappa,
                                          parameterization = "spde",
                                          nu = 1.5)
# Prepare the data for fitting
data_rspde_bru_nonstat <- graph_data_rspde(rspde_model_nonstat,
                                        repl = ".all",
                                        bru = TRUE,
                                        repl_col = "repl")
# Define the component
cmp_nonstat <- y ~ -1 +
  Intercept(1) +
  mean_value +
  field(cbind(.edge_number, .distance_on_edge), 
        model = rspde_model_nonstat,
        replicate = repl)
# Fit the model
rspde_fit_nonstat <-
  bru(cmp_nonstat,
      data = data_rspde_bru_nonstat[["data"]],
      family = "gaussian",
      options = list(verbose = FALSE)
  )

output_from_models <- process_model_results(rspde_fit_nonstat, rspde_model_nonstat)
parameters_statistics <- output_from_models$allparams
mean_and_mode_params_nonstatnu1.5 <- parameters_statistics[, c(1,6)]
rspde_fit_nonstatnu1.5 <- rspde_fit_nonstat
# Summarize the results
summary(rspde_fit_nonstat)

## inlabru version: 2.12.0.9002
## INLA version: 24.12.11
## Components:
## Intercept: main = linear(1), group = exchangeable(1L), replicate = iid(1L), NULL
## mean_value: main = linear(mean_value), group = exchangeable(1L), replicate = iid(1L), NULL
## field: main = cgeneric(cbind(.edge_number, .distance_on_edge)), group = exchangeable(1L), replicate = iid(repl), NULL
## Likelihoods:
##   Family: 'gaussian'
##     Tag: ''
##     Data class: 'metric_graph_data', 'data.frame'
##     Response class: 'numeric'
##     Predictor: y ~ .
##     Used components: effects[Intercept, mean_value, field], latent[]
## Time used:
##     Pre = 0.189, Running = 24.9, Post = 5.04, Total = 30.1 
## Fixed effects:
##              mean    sd 0.025quant 0.5quant 0.975quant   mode kld
## Intercept  -4.456 1.572     -7.533   -4.458     -1.365 -4.458   0
## mean_value  1.124 0.025      1.076    1.124      1.173  1.124   0
## 
## Random effects:
##   Name     Model
##     field CGeneric
## 
## Model hyperparameters:
##                                           mean    sd 0.025quant 0.5quant
## Precision for the Gaussian observations  0.014 0.000      0.013    0.014
## Theta1 for field                        -0.337 0.112     -0.549   -0.340
## Theta2 for field                        -2.015 0.141     -2.307   -2.010
## Theta3 for field                        -0.921 0.057     -1.034   -0.921
## Theta4 for field                         0.689 0.097      0.504    0.687
##                                         0.975quant   mode
## Precision for the Gaussian observations      0.014  0.014
## Theta1 for field                            -0.109 -0.353
## Theta2 for field                            -1.755 -1.987
## Theta3 for field                            -0.811 -0.919
## Theta4 for field                             0.887  0.678
## 
## Deviance Information Criterion (DIC) ...............: 29494.44
## Deviance Information Criterion (DIC, saturated) ....: 4469.92
## Effective number of parameters .....................: 395.27
## 
## Watanabe-Akaike information criterion (WAIC) ...: 29537.62
## Effective number of parameters .................: 386.62
## 
## Marginal log-Likelihood:  -14973.80 
##  is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')

parameters_statistics

---

Press the Show button below to reveal the code.

---

# Plot the estimates of the parameters
ggplot(output_from_models$df_for_plot_spde) + geom_line(aes(x = x, y = y)) + 
  facet_wrap(~parameter, scales = "free") + labs(y = "Density") + 
  theme(text = element_text(family = "Palatino"))

Figure 9: Posterior distributions of the spde parameters.

---

Press the Show button below to reveal the code.

---

ggplot(output_from_models$df_for_plot_matern) + geom_line(aes(x = x, y = y)) + 
  facet_wrap(~parameter, scales = "free") + labs(y = "Density") + 
  theme(text = element_text(family = "Palatino"))

Figure 10: Posterior distributions of the matern parameters.

Case \(\nu\) estimated

We first consider the stationary model.

---

Press the Show button below to reveal the code.

---

# Build the model
rspde_model_stat <- rspde.metric_graph(graph,
                                       parameterization = "spde")
# Prepare the data for fitting
data_rspde_bru_stat <- graph_data_rspde(rspde_model_stat,
                                        repl = ".all",
                                        bru = TRUE,
                                        repl_col = "repl")
# Define the component
cmp_stat <- y ~ -1 +
  Intercept(1) +
  mean_value +
  field(cbind(.edge_number, .distance_on_edge), 
        model = rspde_model_stat,
        replicate = repl)
# Fit the model
rspde_fit_stat <-
  bru(cmp_stat,
      data = data_rspde_bru_stat[["data"]],
      family = "gaussian",
      options = list(verbose = FALSE)
  )

output_from_models <-process_model_results(rspde_fit_stat, rspde_model_stat)
parameters_statistics <- output_from_models$allparams
mean_and_mode_params_statnuest <- parameters_statistics[, c(1,6)]
rspde_fit_statnuest <- rspde_fit_stat
# Summarize the results
summary(rspde_fit_stat)

## inlabru version: 2.12.0.9002
## INLA version: 24.12.11
## Components:
## Intercept: main = linear(1), group = exchangeable(1L), replicate = iid(1L), NULL
## mean_value: main = linear(mean_value), group = exchangeable(1L), replicate = iid(1L), NULL
## field: main = cgeneric(cbind(.edge_number, .distance_on_edge)), group = exchangeable(1L), replicate = iid(repl), NULL
## Likelihoods:
##   Family: 'gaussian'
##     Tag: ''
##     Data class: 'metric_graph_data', 'data.frame'
##     Response class: 'numeric'
##     Predictor: y ~ .
##     Used components: effects[Intercept, mean_value, field], latent[]
## Time used:
##     Pre = 0.222, Running = 50.5, Post = 13.6, Total = 64.3 
## Fixed effects:
##               mean    sd 0.025quant 0.5quant 0.975quant    mode kld
## Intercept  -14.452 1.532    -17.476  -14.448    -11.457 -14.448   0
## mean_value   1.276 0.020      1.236    1.276      1.316   1.276   0
## 
## Random effects:
##   Name     Model
##     field CGeneric
## 
## Model hyperparameters:
##                                           mean    sd 0.025quant 0.5quant
## Precision for the Gaussian observations  0.014 0.000      0.013    0.014
## Theta1 for field                        -1.496 0.042     -1.575   -1.497
## Theta2 for field                        -2.971 0.078     -3.137   -2.967
## Theta3 for field                        -1.091 0.037     -1.159   -1.092
##                                         0.975quant   mode
## Precision for the Gaussian observations      0.014  0.014
## Theta1 for field                            -1.410 -1.503
## Theta2 for field                            -2.832 -2.947
## Theta3 for field                            -1.013 -1.100
## 
## Deviance Information Criterion (DIC) ...............: 29846.42
## Deviance Information Criterion (DIC, saturated) ....: 4792.55
## Effective number of parameters .....................: 713.61
## 
## Watanabe-Akaike information criterion (WAIC) ...: 29814.15
## Effective number of parameters .................: 597.47
## 
## Marginal log-Likelihood:  -15152.96 
##  is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')

parameters_statistics

---

Press the Show button below to reveal the code.

---

# Plot the estimates of the parameters
ggplot(output_from_models$df_for_plot_spde) + geom_line(aes(x = x, y = y)) + 
  facet_wrap(~parameter, scales = "free") + labs(y = "Density") + 
  theme(text = element_text(family = "Palatino"))

Figure 11: Posterior distributions of the spde parameters.

---

Press the Show button below to reveal the code.

---

ggplot(output_from_models$df_for_plot_matern) + geom_line(aes(x = x, y = y)) + 
  facet_wrap(~parameter, scales = "free") + labs(y = "Density") + 
  theme(text = element_text(family = "Palatino"))

Figure 12: Posterior distributions of the matern parameters.

We now fit the non-stationary model.

---

Press the Show button below to reveal the code.

---

# Build the model
rspde_model_nonstat <- rspde.metric_graph(graph,
                                          B.tau = B.tau,
                                          B.kappa =  B.kappa,
                                          parameterization = "spde")
# Prepare the data for fitting
data_rspde_bru_nonstat <- graph_data_rspde(rspde_model_nonstat,
                                        repl = ".all",
                                        bru = TRUE,
                                        repl_col = "repl")
# Define the component
cmp_nonstat <- y ~ -1 +
  Intercept(1) +
  mean_value +
  field(cbind(.edge_number, .distance_on_edge), 
        model = rspde_model_nonstat,
        replicate = repl)
# Fit the model
rspde_fit_nonstat <-
  bru(cmp_nonstat,
      data = data_rspde_bru_nonstat[["data"]],
      family = "gaussian",
      options = list(verbose = FALSE)
  )

output_from_models <- process_model_results(rspde_fit_nonstat, rspde_model_nonstat)
parameters_statistics <- output_from_models$allparams
mean_and_mode_params_nonstatnuest <- parameters_statistics[, c(1,6)]
rspde_fit_nonstatnuest <- rspde_fit_nonstat
# Summarize the results
summary(rspde_fit_nonstat)

## inlabru version: 2.12.0.9002
## INLA version: 24.12.11
## Components:
## Intercept: main = linear(1), group = exchangeable(1L), replicate = iid(1L), NULL
## mean_value: main = linear(mean_value), group = exchangeable(1L), replicate = iid(1L), NULL
## field: main = cgeneric(cbind(.edge_number, .distance_on_edge)), group = exchangeable(1L), replicate = iid(repl), NULL
## Likelihoods:
##   Family: 'gaussian'
##     Tag: ''
##     Data class: 'metric_graph_data', 'data.frame'
##     Response class: 'numeric'
##     Predictor: y ~ .
##     Used components: effects[Intercept, mean_value, field], latent[]
## Time used:
##     Pre = 0.216, Running = 409, Post = 12.2, Total = 421 
## Fixed effects:
##              mean    sd 0.025quant 0.5quant 0.975quant   mode kld
## Intercept  -4.176 1.718     -7.541   -4.178     -0.800 -4.178   0
## mean_value  1.122 0.026      1.070    1.122      1.174  1.122   0
## 
## Random effects:
##   Name     Model
##     field CGeneric
## 
## Model hyperparameters:
##                                           mean    sd 0.025quant 0.5quant
## Precision for the Gaussian observations  0.015 0.000      0.014    0.015
## Theta1 for field                         0.136 0.070     -0.006    0.138
## Theta2 for field                        -2.423 0.066     -2.547   -2.424
## Theta3 for field                        -0.734 0.044     -0.821   -0.734
## Theta4 for field                         0.664 0.060      0.545    0.664
## Theta5 for field                         0.944 0.065      0.817    0.944
##                                         0.975quant   mode
## Precision for the Gaussian observations      0.015  0.015
## Theta1 for field                             0.270  0.145
## Theta2 for field                            -2.289 -2.431
## Theta3 for field                            -0.649 -0.733
## Theta4 for field                             0.783  0.665
## Theta5 for field                             1.072  0.945
## 
## Deviance Information Criterion (DIC) ...............: 29418.94
## Deviance Information Criterion (DIC, saturated) ....: 4672.34
## Effective number of parameters .....................: 588.94
## 
## Watanabe-Akaike information criterion (WAIC) ...: 29453.51
## Effective number of parameters .................: 534.73
## 
## Marginal log-Likelihood:  -14959.25 
##  is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')

parameters_statistics

---

Press the Show button below to reveal the code.

---

# Plot the estimates of the parameters
ggplot(output_from_models$df_for_plot_spde) + geom_line(aes(x = x, y = y)) + 
  facet_wrap(~parameter, scales = "free") + labs(y = "Density") + 
  theme(text = element_text(family = "Palatino"))

Figure 13: Posterior distributions of the spde parameters.

---

Press the Show button below to reveal the code.

---

ggplot(output_from_models$df_for_plot_matern) + geom_line(aes(x = x, y = y)) + 
  facet_wrap(~parameter, scales = "free") + labs(y = "Density") + 
  theme(text = element_text(family = "Palatino"))

Figure 14: Posterior distributions of the matern parameters.

Below we perform leave-group-out pseudo cross-validation (Liu and Rue 2022) following the strategy from (Bolin, Simas, and Xiong 2023).

---

Press the Show button below to reveal the code.

---

mse.statnu0.5 <- mse.nonstatnu0.5 <- ls.statnu0.5 <- ls.nonstatnu0.5 <- rep(0,length(distance))
mse.statnu1.5 <- mse.nonstatnu1.5 <- ls.statnu1.5 <- ls.nonstatnu1.5 <- rep(0,length(distance))
mse.statnuest <- mse.nonstatnuest <- ls.statnuest <- ls.nonstatnuest <- rep(0,length(distance))

# cross-validation for-loop
for (j in 1:length(distance)) {
  print(j)
  # cross-validation of the stationary model
  cv.statnu0.5 <- inla.group.cv(rspde_fit_statnu0.5, groups = GROUPS[[j]])
  cv.statnu1.5 <- inla.group.cv(rspde_fit_statnu1.5, groups = GROUPS[[j]])
  cv.statnuest <- inla.group.cv(rspde_fit_statnuest, groups = GROUPS[[j]])
  # cross-validation of the nonstationary model
  cv.nonstatnu0.5 <- inla.group.cv(rspde_fit_nonstatnu0.5, groups = GROUPS[[j]])
  cv.nonstatnu1.5 <- inla.group.cv(rspde_fit_nonstatnu1.5, groups = GROUPS[[j]])
  cv.nonstatnuest <- inla.group.cv(rspde_fit_nonstatnuest, groups = GROUPS[[j]])
  # obtain MSE and LS
  mse.statnu0.5[j] <- mean((cv.statnu0.5$mean - data$y)^2)
  mse.statnu1.5[j] <- mean((cv.statnu1.5$mean - data$y)^2)
  mse.statnuest[j] <- mean((cv.statnuest$mean - data$y)^2)
  
  
  mse.nonstatnu0.5[j] <- mean((cv.nonstatnu0.5$mean - data$y)^2)
  mse.nonstatnu1.5[j] <- mean((cv.nonstatnu1.5$mean - data$y)^2)
  mse.nonstatnuest[j] <- mean((cv.nonstatnuest$mean - data$y)^2)
  
  
  ls.statnu0.5[j] <- mean(log(cv.statnu0.5$cv))
  ls.statnu1.5[j] <- mean(log(cv.statnu1.5$cv))
  ls.statnuest[j] <- mean(log(cv.statnuest$cv))
  
  ls.nonstatnu0.5[j] <- mean(log(cv.nonstatnu0.5$cv))
  ls.nonstatnu1.5[j] <- mean(log(cv.nonstatnu1.5$cv))
  ls.nonstatnuest[j] <- mean(log(cv.nonstatnuest$cv))
}
## [1] 1
## [1] 2
## [1] 3
## [1] 4
## [1] 5
## [1] 6
## [1] 7
## [1] 8
## [1] 9
## [1] 10
## [1] 11
## [1] 12
## [1] 13
## [1] 14
## [1] 15
## [1] 16
## [1] 17
## [1] 18
## [1] 19
## [1] 20
## [1] 21
## [1] 22
## [1] 23
## [1] 24
## [1] 25
## [1] 26
## [1] 27
## [1] 28
## [1] 29
## [1] 30
## [1] 31
## [1] 32
## [1] 33
## [1] 34
## [1] 35
## [1] 36
## [1] 37
## [1] 38
## [1] 39
## [1] 40
## [1] 41
## [1] 42
## [1] 43
## [1] 44
## [1] 45
## [1] 46
## [1] 47
## [1] 48
## [1] 49
## [1] 50
## [1] 51
## [1] 52
## [1] 53
## [1] 54
## [1] 55
## [1] 56
## [1] 57
## [1] 58
## [1] 59
## [1] 60
## [1] 61
## [1] 62
## [1] 63
## [1] 64
## [1] 65
## [1] 66
## [1] 67
## [1] 68
## [1] 69
## [1] 70
## [1] 71
## [1] 72
## [1] 73
## [1] 74
## [1] 75
## [1] 76
## [1] 77
## [1] 78
## [1] 79
## [1] 80
## [1] 81
## [1] 82
## [1] 83
## [1] 84
## [1] 85
## [1] 86
## [1] 87
## [1] 88
## [1] 89
## [1] 90
## [1] 91
## [1] 92
## [1] 93
## [1] 94
## [1] 95
## [1] 96
## [1] 97
## [1] 98
## [1] 99
## [1] 100
## [1] 101

# Create data frames
mse_df <- data.frame(
  distance,
  Statnu0.5 = mse.statnu0.5,
  Nonstatnu0.5 = mse.nonstatnu0.5,
  Statnu1.5 = mse.statnu1.5,
  Nonstatnu1.5 = mse.nonstatnu1.5,
  Statnuest = mse.statnuest,
  Nonstatnuest = mse.nonstatnuest
)

ls_df <- data.frame(
  distance,
  Statnu0.5 = -ls.statnu0.5,
  Nonstatnu0.5 = -ls.nonstatnu0.5,
  Statnu1.5 = -ls.statnu1.5,
  Nonstatnu1.5 = -ls.nonstatnu1.5,
  Statnuest = -ls.statnuest,
  Nonstatnuest = -ls.nonstatnuest
)

Below we plot the cross-validation results.

---

Press the Show button below to reveal the code.

---

choose_index <- seq(2, nrow(mse_df), by = 3)
mse_df_red <- mse_df[choose_index,]
ls_df_red <- ls_df[choose_index,]
# Convert to long format
mse_long <- mse_df_red %>%
  pivot_longer(cols = -distance, names_to = "nu", values_to = "MSE")

ls_long <- ls_df_red %>%
  pivot_longer(cols = -distance, names_to = "nu", values_to = "LogScore")


# Update the label mappings with the new legend title
label_mapping <- c(
  "Statnu0.5" = "1", 
  "Nonstatnu0.5" = "1", 
  "Statnu1.5" = "2", 
  "Nonstatnu1.5" = "2", 
  "Statnuest" = paste(round(mean_and_mode_params_statnuest[5,1]+0.5, 3), "(est)"), 
  "Nonstatnuest" = paste(round(mean_and_mode_params_nonstatnuest[7,1]+0.5, 3), "(est)")
)

# Define color and linetype mapping
color_mapping <- c(
  "Statnu0.5" = "blue", 
  "Nonstatnu0.5" = "blue", 
  "Statnu1.5" = "black", 
  "Nonstatnu1.5" = "black", 
  "Statnuest" = "red", 
  "Nonstatnuest" = "red"
)

linetype_mapping <- c(
  "Statnu0.5" = "dotdash", 
  "Nonstatnu0.5" = "solid", 
  "Statnu1.5" = "dotdash", 
  "Nonstatnu1.5" = "solid", 
  "Statnuest" = "dotdash", 
  "Nonstatnuest" = "solid"
)

# Plot MSE
mse_plot <- ggplot(mse_long, aes(x = distance, y = MSE, color = nu, linetype = nu)) +
  geom_line(linewidth = 1) +
  labs(y = "MSE", x = "Distance in km") +
  scale_color_manual(values = color_mapping, labels = label_mapping, name = expression(alpha)) +
  scale_linetype_manual(values = linetype_mapping, labels = label_mapping, name = expression(alpha)) +
  theme_minimal() +
  theme(text = element_text(family = "Palatino"))

# Plot negative log-score
ls_plot <- ggplot(ls_long, aes(x = distance, y = LogScore, color = nu, linetype = nu)) +
  geom_line(linewidth = 1) +
  labs(y = "Negative Log-Score", x = "Distance in km") +
  scale_color_manual(values = color_mapping, labels = label_mapping, name = expression(alpha)) +
  scale_linetype_manual(values = linetype_mapping, labels = label_mapping, name = expression(alpha)) +
  theme_minimal() +
  theme(text = element_text(family = "Palatino"))

# Combine plots with a shared legend at the top in a single line
combined_plot <- mse_plot + ls_plot + 
  plot_layout(guides = 'collect') & 
  theme(legend.position = 'right') & 
  guides(color = guide_legend(ncol = 1), linetype = guide_legend(nrow = 1))

# Save combined plot
ggsave(here("data_files/crossval_pems.png"), plot = combined_plot, width = 9.22, height = 4.01, dpi = 500)
# Display combined plot
print(combined_plot)

Figure 15: MSE and negative Log-Score as functions of distance (in km) for the stationary (dotdash line, \(\boldsymbol{\cdot-\cdot}\)) and non-stationary (solid line, \(\boldsymbol{-\!\!\!-\!\!\!-}\)) cases with \(\nu = 0.5\), \(\nu = 1.5\), and \(\nu\) estimated (est).

Save some of the objects to be used in the next vignette.

# Save the results
list_to_save <- list(mean_and_mode_params_statnu0.5 = mean_and_mode_params_statnu0.5,
                     mean_and_mode_params_nonstatnu0.5 = mean_and_mode_params_nonstatnu0.5,
                     mean_and_mode_params_statnu1.5 = mean_and_mode_params_statnu1.5,
                     mean_and_mode_params_nonstatnu1.5 = mean_and_mode_params_nonstatnu1.5, 
                     mean_and_mode_params_statnuest = mean_and_mode_params_statnuest,
                     mean_and_mode_params_nonstatnuest = mean_and_mode_params_nonstatnuest, 
                     mse_df = mse_df, 
                     ls_df = ls_df, 
                     B.tau = B.tau, 
                     B.kappa = B.kappa, 
                     graph = graph)
save(list_to_save, file = here("data_files/pems_repl2_results.RData"))

References

cite_packages(output = "paragraph", out.dir = ".")

We used R version 4.4.1 (R Core Team 2024) and the following R packages: cowplot v. 1.1.3 (Wilke 2024), ggmap v. 4.0.0.900 (Kahle and Wickham 2013), ggpubr v. 0.6.0 (Kassambara 2023), ggtext v. 0.1.2 (Wilke and Wiernik 2022), here v. 1.0.1 (Müller 2020), htmltools v. 0.5.8.1 (Cheng et al. 2024), INLA v. 24.12.11 (Rue, Martino, and Chopin 2009; Lindgren, Rue, and Lindström 2011; Martins et al. 2013; Lindgren and Rue 2015; De Coninck et al. 2016; Rue et al. 2017; Verbosio et al. 2017; Bakka et al. 2018; Kourounis, Fuchs, and Schenk 2018), inlabru v. 2.12.0.9002 (Yuan et al. 2017; Bachl et al. 2019), knitr v. 1.48 (Xie 2014, 2015, 2024), latex2exp v. 0.9.6 (Meschiari 2022), Matrix v. 1.6.5 (Bates, Maechler, and Jagan 2024), MetricGraph v. 1.4.0.9000 (Bolin, Simas, and Wallin 2023b, 2023a, 2023c, 2024; Bolin et al. 2024), OpenStreetMap v. 0.4.0 (Fellows and JMapViewer library by Jan Peter Stotz 2023), osmdata v. 0.2.5 (Mark Padgham et al. 2017), patchwork v. 1.2.0 (Pedersen 2024), plotly v. 4.10.4 (Sievert 2020), plotrix v. 3.8.4 (J 2006), reshape2 v. 1.4.4 (Wickham 2007), rmarkdown v. 2.28 (Xie, Allaire, and Grolemund 2018; Xie, Dervieux, and Riederer 2020; Allaire et al. 2024), rSPDE v. 2.4.0.9000 (Bolin and Kirchner 2020; Bolin and Simas 2023; Bolin, Simas, and Xiong 2024), scales v. 1.3.0 (Wickham, Pedersen, and Seidel 2023), sf v. 1.0.19 (E. Pebesma 2018; E. Pebesma and Bivand 2023), sp v. 2.1.4 (E. J. Pebesma and Bivand 2005; Bivand, Pebesma, and Gomez-Rubio 2013), tidyverse v. 2.0.0 (Wickham et al. 2019), viridis v. 0.6.4 (Garnier et al. 2023), xaringanExtra v. 0.8.0 (Aden-Buie and Warkentin 2024).

Aden-Buie, Garrick, and Matthew T. Warkentin. 2024. xaringanExtra: Extras and Extensions for “xaringan” Slides. https://CRAN.R-project.org/package=xaringanExtra.

Allaire, JJ, Yihui Xie, Christophe Dervieux, Jonathan McPherson, Javier Luraschi, Kevin Ushey, Aron Atkins, et al. 2024. rmarkdown: Dynamic Documents for r. https://github.com/rstudio/rmarkdown.

Bachl, Fabian E., Finn Lindgren, David L. Borchers, and Janine B. Illian. 2019. “inlabru: An R Package for Bayesian Spatial Modelling from Ecological Survey Data.” Methods in Ecology and Evolution 10: 760–66. https://doi.org/10.1111/2041-210X.13168.

Bakka, Haakon, Håvard Rue, Geir-Arne Fuglstad, Andrea I. Riebler, David Bolin, Janine Illian, Elias Krainski, Daniel P. Simpson, and Finn K. Lindgren. 2018. “Spatial Modelling with INLA: A Review.” WIRES (Invited Extended Review) xx (Feb): xx–. http://arxiv.org/abs/1802.06350.

Bates, Douglas, Martin Maechler, and Mikael Jagan. 2024. Matrix: Sparse and Dense Matrix Classes and Methods. https://CRAN.R-project.org/package=Matrix.

Bivand, Roger S., Edzer Pebesma, and Virgilio Gomez-Rubio. 2013. Applied Spatial Data Analysis with R, Second Edition. Springer, NY. https://asdar-book.org/.

Bolin, David, and Kristin Kirchner. 2020. “The Rational SPDE Approach for Gaussian Random Fields with General Smoothness.” Journal of Computational and Graphical Statistics 29 (2): 274–85. https://doi.org/10.1080/10618600.2019.1665537.

Bolin, David, Mihály Kovács, Vivek Kumar, and Alexandre B. Simas. 2024. “Regularity and Numerical Approximation of Fractional Elliptic Differential Equations on Compact Metric Graphs.” Mathematics of Computation 93 (349): 2439–72. https://doi.org/10.1090/mcom/3929.

Bolin, David, and Alexandre B. Simas. 2023. rSPDE: Rational Approximations of Fractional Stochastic Partial Differential Equations. https://CRAN.R-project.org/package=rSPDE.

Bolin, David, Alexandre B. Simas, and Jonas Wallin. 2023a. “Markov Properties of Gaussian Random Fields on Compact Metric Graphs.” arXiv Preprint arXiv:2304.03190. https://doi.org/10.48550/arXiv.2304.03190.

———. 2023b. MetricGraph: Random Fields on Metric Graphs. https://CRAN.R-project.org/package=MetricGraph.

———. 2023c. “Statistical Inference for Gaussian Whittle-Matérn Fields on Metric Graphs.” arXiv Preprint arXiv:2304.10372. https://doi.org/10.48550/arXiv.2304.10372.

———. 2024. “Gaussian Whittle-Matérn Fields on Metric Graphs.” Bernoulli 30 (2): 1611–39. https://doi.org/10.3150/23-BEJ1647.

Bolin, David, Alexandre B. Simas, and Zhen Xiong. 2024. “Covariance-Based Rational Approximations of Fractional SPDEs for Computationally Efficient Bayesian Inference.” Journal of Computational and Graphical Statistics 33 (1): 64–74. https://doi.org/10.1080/10618600.2023.2231051.

Bolin, David, Alexandre B Simas, and Zhen Xiong. 2023. “Covariance–Based Rational Approximations of Fractional SPDEs for Computationally Efficient Bayesian Inference.” Journal of Computational and Graphical Statistics, 1–11.

Cheng, Joe, Carson Sievert, Barret Schloerke, Winston Chang, Yihui Xie, and Jeff Allen. 2024. htmltools: Tools for HTML. https://CRAN.R-project.org/package=htmltools.

De Coninck, Arne, Bernard De Baets, Drosos Kourounis, Fabio Verbosio, Olaf Schenk, Steven Maenhout, and Jan Fostier. 2016. “Needles: Toward Large-Scale Genomic Prediction with Marker-by-Environment Interaction.” Genetics 203 (1): 543–55. https://doi.org/10.1534/genetics.115.179887.

Fellows, Ian, and using the JMapViewer library by Jan Peter Stotz. 2023. OpenStreetMap: Access to Open Street Map Raster Images. https://CRAN.R-project.org/package=OpenStreetMap.

Garnier, Simon, Ross, Noam, Rudis, Robert, Camargo, et al. 2023. viridis(Lite) - Colorblind-Friendly Color Maps for r. https://doi.org/10.5281/zenodo.4679423.

J, Lemon. 2006. “Plotrix: A Package in the Red Light District of r.” R-News 6 (4): 8–12.

Kahle, David, and Hadley Wickham. 2013. “ggmap: Spatial Visualization with Ggplot2.” The R Journal 5 (1): 144–61. https://journal.r-project.org/archive/2013-1/kahle-wickham.pdf.

Kassambara, Alboukadel. 2023. ggpubr: “ggplot2” Based Publication Ready Plots. https://CRAN.R-project.org/package=ggpubr.

Kourounis, D., A. Fuchs, and O. Schenk. 2018. “Towards the Next Generation of Multiperiod Optimal Power Flow Solvers.” IEEE Transactions on Power Systems PP (99): 1–10. https://doi.org/10.1109/TPWRS.2017.2789187.

Lindgren, Finn, and Håvard Rue. 2015. “Bayesian Spatial Modelling with R-INLA.” Journal of Statistical Software 63 (19): 1–25. http://www.jstatsoft.org/v63/i19/.

Lindgren, Finn, Håvard Rue, and Johan Lindström. 2011. “An Explicit Link Between Gaussian Fields and Gaussian Markov Random Fields: The Stochastic Partial Differential Equation Approach (with Discussion).” Journal of the Royal Statistical Society B 73 (4): 423–98.

Liu, Zhedong, and Haavard Rue. 2022. “Leave-Group-Out Cross-Validation for Latent Gaussian Models.” arXiv Preprint arXiv:2210.04482.

Mark Padgham, Bob Rudis, Robin Lovelace, and Maëlle Salmon. 2017. “Osmdata.” Journal of Open Source Software 2 (14): 305. https://doi.org/10.21105/joss.00305.

Martins, Thiago G., Daniel Simpson, Finn Lindgren, and Håvard Rue. 2013. “Bayesian Computing with INLA: New Features.” Computational Statistics and Data Analysis 67: 68–83.

Meschiari, Stefano. 2022. Latex2exp: Use LaTeX Expressions in Plots. https://CRAN.R-project.org/package=latex2exp.

Müller, Kirill. 2020. here: A Simpler Way to Find Your Files. https://CRAN.R-project.org/package=here.

Pebesma, Edzer. 2018. “Simple Features for R: Standardized Support for Spatial Vector Data.” The R Journal 10 (1): 439–46. https://doi.org/10.32614/RJ-2018-009.

Pebesma, Edzer J., and Roger Bivand. 2005. “Classes and Methods for Spatial Data in R.” R News 5 (2): 9–13. https://CRAN.R-project.org/doc/Rnews/.

Pebesma, Edzer, and Roger Bivand. 2023. Spatial Data Science: With applications in R. Chapman and Hall/CRC. https://doi.org/10.1201/9780429459016.

Pedersen, Thomas Lin. 2024. patchwork: The Composer of Plots. https://CRAN.R-project.org/package=patchwork.

R Core Team. 2024. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. https://www.R-project.org/.

Rue, Håvard, Sara Martino, and Nicholas Chopin. 2009. “Approximate Bayesian Inference for Latent Gaussian Models Using Integrated Nested Laplace Approximations (with Discussion).” Journal of the Royal Statistical Society B 71: 319–92.

Rue, Håvard, Andrea I. Riebler, Sigrunn H. Sørbye, Janine B. Illian, Daniel P. Simpson, and Finn K. Lindgren. 2017. “Bayesian Computing with INLA: A Review.” Annual Reviews of Statistics and Its Applications 4 (March): 395–421. http://arxiv.org/abs/1604.00860.

Sievert, Carson. 2020. Interactive Web-Based Data Visualization with r, Plotly, and Shiny. Chapman; Hall/CRC. https://plotly-r.com.

Verbosio, Fabio, Arne De Coninck, Drosos Kourounis, and Olaf Schenk. 2017. “Enhancing the Scalability of Selected Inversion Factorization Algorithms in Genomic Prediction.” Journal of Computational Science 22 (Supplement C): 99–108. https://doi.org/10.1016/j.jocs.2017.08.013.

Wickham, Hadley. 2007. “Reshaping Data with the reshape Package.” Journal of Statistical Software 21 (12): 1–20. http://www.jstatsoft.org/v21/i12/.

Wickham, Hadley, Mara Averick, Jennifer Bryan, Winston Chang, Lucy D’Agostino McGowan, Romain François, Garrett Grolemund, et al. 2019. “Welcome to the tidyverse.” Journal of Open Source Software 4 (43): 1686. https://doi.org/10.21105/joss.01686.

Wickham, Hadley, Thomas Lin Pedersen, and Dana Seidel. 2023. scales: Scale Functions for Visualization. https://CRAN.R-project.org/package=scales.

Wilke, Claus O. 2024. cowplot: Streamlined Plot Theme and Plot Annotations for “ggplot2”. https://CRAN.R-project.org/package=cowplot.

Wilke, Claus O., and Brenton M. Wiernik. 2022. ggtext: Improved Text Rendering Support for “ggplot2”. https://CRAN.R-project.org/package=ggtext.

Xie, Yihui. 2014. “knitr: A Comprehensive Tool for Reproducible Research in R.” In Implementing Reproducible Computational Research, edited by Victoria Stodden, Friedrich Leisch, and Roger D. Peng. Chapman; Hall/CRC.

———. 2015. Dynamic Documents with R and Knitr. 2nd ed. Boca Raton, Florida: Chapman; Hall/CRC. https://yihui.org/knitr/.

———. 2024. knitr: A General-Purpose Package for Dynamic Report Generation in r. https://yihui.org/knitr/.

Xie, Yihui, J. J. Allaire, and Garrett Grolemund. 2018. R Markdown: The Definitive Guide. Boca Raton, Florida: Chapman; Hall/CRC. https://bookdown.org/yihui/rmarkdown.

Xie, Yihui, Christophe Dervieux, and Emily Riederer. 2020. R Markdown Cookbook. Boca Raton, Florida: Chapman; Hall/CRC. https://bookdown.org/yihui/rmarkdown-cookbook.

Yuan, Yuan, Bachl, Fabian E., Lindgren, Finn, Borchers, et al. 2017. “Point Process Models for Spatio-Temporal Distance Sampling Data from a Large-Scale Survey of Blue Whales.” Ann. Appl. Stat. 11 (4): 2270–97. https://doi.org/10.1214/17-AOAS1078.