PeMS 1, covariates and preprocessing

Go back to the Contents page.

---

Press Show to reveal the code chunks.

---

Go back to the About page.

This vignette processes half of the PeMS dataset, shows some illustrations, and produces pems_repl1_data.RData, which is a file that contains the graph with data and a covariate defined on the mesh. To see how the data is modeled, go to pems_repl2.html.

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

# Create a clipboard button on the rendered HTML page
source(here::here("clipboard.R")); clipboard

# Set seed for reproducibility
set.seed(1982) 
# 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 = FALSE,       
  # 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
)
# Start the figure counter
fig_count <- 0
# Define the captioner function
captioner <- function(caption) {
  fig_count <<- fig_count + 1
  paste0("Figure ", fig_count, ": ", caption)
}
# Define the function to truncate a number to two decimal places
truncate_to_two <- function(x) {
  truncated <- floor(x * 100) / 100
  sprintf("%.2f", truncated)
}

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

library(slackr)
source("keys.R")
slackr_setup(token = token) # token comes from keys.R

## [1] "Successfully connected to Slack"

capture.output(
  knitr::purl(here::here("functionality1.Rmd"), output = here::here("functionality1.R")),
  file = here::here("old/purl_log.txt")
)
source(here::here("functionality1.R"))

Below we define function gets_summary_parameters() to extract the summary of the parameters of the model.

gets_summary_parameters <- function(fit, model) {
  param_spde <- summary(rspde.result(fit, "field", model, parameterization = "spde"))
  param_matern <- summary(rspde.result(fit, "field", model, parameterization = "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(allparams)
}

---

Folder data.pems was downloaded from this GitHub repository.

Lines <- read_sf(here("data_files/data.pems/lines.shp"))
lines <- as_Spatial(Lines)

EtV <- read.csv(here("data_files/data.pems/E.csv"), header = T, row.names = NULL)
PtE <- read.csv(here("data_files/data.pems/PtE.csv"), header = T, row.names = NULL)
PtE[,1] <- PtE[,1] + 1
Y <- read.csv(here("data_files/data.pems/Y.csv"), header = T, row.names = NULL)
Y <- as.matrix(Y[,-1])
edge_length_m <- EtV[,4]
PtE[,2] = PtE[,2]/edge_length_m[PtE[,1]]

Matrix Y has dimension \(26\times325\). Each row of Y corresponds to a replicate (\(26\) replicates in total) and each column corresponds to a location (\(325\) locations in total) on the network.

We remove columns in Y, those with the same value in the first 13 replicates. We also remove the corresponding rows in PtE.

all_same <- function(col) {
  length(unique(col[1:13])) == 1
}
cols_to_keep <- apply(Y, 2, function(col) !all_same(col))

Y <- Y[, cols_to_keep]
PtE <- PtE[cols_to_keep,]

Half of the replicates (Y_for_summary) are used to compute Y_logstd and Y_mean, and the other half (Y again ) are used to fit a replicate model in pems_repl2.html.

sampled_numbers <- 1:13
not_sampled_numbers <- 14:26
Y_for_summary <- Y[sampled_numbers,]
Y_logstd <- apply(Y_for_summary, 2, sd) |> as.vector() |> log()
Y_mean <- apply(Y_for_summary, 2, mean) |> as.vector()
Y <- Y[not_sampled_numbers,]
save(Y_mean, file = here::here("data_files/Y_mean.RData"))

---

Observe that Y_logstd corresponds to \(\log\) of the standard deviation at each location.

Below we build the graph using the list element edges object from pems_repl.

# Build the graph
graph <- metric_graph$new(edges = pems_repl$edges, longlat = TRUE)
# Add the observations
graph$add_observations(data = data.frame(y = Y_logstd, 
                                         edge_number = PtE[,1], 
                                         distance_on_edge = PtE[,2]),
                       edge_number = "edge_number",
                       distance_on_edge = "distance_on_edge",
                       data_coords = "PtE",
                       normalized = TRUE, 
                       clear_obs = TRUE)
# Build the mesh
graph$build_mesh(h = 0.05) 

Below we fit \(\log(\mathrm{sd}(s_i))\sim N(\beta_0+u(s_i), \sigma_e^2)\).

# Build the model
rspde_model_stat_logstd <- rspde.metric_graph(graph,
                                       parameterization = "spde")
# Prepare the data for fitting
data_rspde_bru_stat <- graph_data_rspde(rspde_model_stat_logstd, bru = TRUE)
# Define the component
cmp_stat <- y ~ -1 +
  Intercept(1) +
  field(cbind(.edge_number, .distance_on_edge), model = rspde_model_stat_logstd)
# Fit the model
rspde_fit_stat_logstd <-
  bru(cmp_stat,
      data = data_rspde_bru_stat[["data"]],
      family = "gaussian",
      options = list(verbose = FALSE)
  )
save(rspde_fit_stat_logstd, rspde_model_stat_logstd, file = here::here("data_files/rspde_fit_stat_and_model_logstd.RData"))

load(here::here("data_files/rspde_fit_stat_and_model_logstd.RData"))
# Summarize the results
summary(rspde_fit_stat_logstd)

## inlabru version: 2.13.0
## INLA version: 25.11.22
## Components:
## Intercept: main = linear(1), group = exchangeable(1L), replicate = iid(1L), NULL
## field: main = cgeneric(cbind(.edge_number, .distance_on_edge)), group = exchangeable(1L), replicate = iid(1L), NULL
## Observation models:
##   Family: 'gaussian'
##     Tag: <No tag>
##     Data class: 'metric_graph_data', 'data.frame'
##     Response class: 'numeric'
##     Predictor: y ~ .
##     Additive/Linear: TRUE/TRUE
##     Used components: effects[Intercept, field], latent[]
## Time used:
##     Pre = 0.549, Running = 3.27, Post = 1.19, Total = 5.01 
## Fixed effects:
##            mean    sd 0.025quant 0.5quant 0.975quant  mode kld
## Intercept 1.756 0.133      1.491    1.756       2.02 1.756   0
## 
## Random effects:
##   Name     Model
##     field CGeneric
## 
## Model hyperparameters:
##                                           mean    sd 0.025quant 0.5quant
## Precision for the Gaussian observations  4.600 0.504      3.685    4.574
## Theta1 for field                         0.261 0.193     -0.120    0.261
## Theta2 for field                        -0.647 0.224     -1.084   -0.649
## Theta3 for field                         0.886 0.688     -0.432    0.873
##                                         0.975quant   mode
## Precision for the Gaussian observations      5.667  4.523
## Theta1 for field                             0.642  0.261
## Theta2 for field                            -0.201 -0.655
## Theta3 for field                             2.278  0.819
## 
## Marginal log-Likelihood:  -327.36 
##  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)')

gets_summary_parameters(rspde_fit_stat_logstd, rspde_model_stat_logstd)

We now compute the kriging predictor \(E(u(s_i)|\log(\mathrm{sd}(s_i)))\) (denoted as cov below) and standardize it.

# Prediction locations
data_prd_list_mesh <- graph$get_mesh_locations(bru = TRUE)
# Compute the kriging predictor
y_pred <- predict(rspde_fit_stat_logstd, newdata = data_prd_list_mesh, ~Intercept + field)
cov_almost <- y_pred$mean
# Standardize the kriging predictor of the log standard deviation
cov <- (cov_almost - mean(cov_almost))/sd(cov_almost)

Below we plot the results. We first generate two maps, one for zoom level 12 and the other for zoom level 13.

# 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]

# Set the style of the map
style_vector <- c(c(feature = "road", element = "geometry", visibility = "on"),
                  c("&style=", feature = "poi", element = "labels", visibility = "off"),
                  c("&style=", feature = "road", element = "labels", visibility = "off"),
                  c("&style=", feature = "administrative", element = "labels", visibility = "off"), #name of places
                  c("&style=", feature = "transit", element = "labels", visibility = "off"),
                  c("&style=", feature = "landscape", element = "geometry", visibility = "on"),
                  c("&style=", feature = "water", element = "geometry", visibility = "on")) 
register_stadiamaps(Sys.getenv("GGMAP_STADIAMAPS_API_KEY"))
# Obtain the maps
p12 <- ggmap(get_stadiamap(bbox = c(left = x_left, bottom = y_bottom, right = x_right, top = y_top),
              zoom = 12,
              maptype = "stamen_toner_lite", 
              scale =2,
              style = style_vector,
              color = "color")) + xlab(NULL) + ylab(NULL)
p13 <- ggmap(get_stadiamap(bbox = c(left = x_left, bottom = y_bottom, right = x_right, top = y_top),
              zoom = 13,
              maptype = "stamen_toner_lite", 
              scale =2,
              style = style_vector,
              color = "color")) + xlab(NULL) + ylab(NULL)

save(p12, p13, file = here::here("data_files/maps_zoom12and13from_stadia.RData"))

load(here::here("data_files/maps_zoom12and13from_stadia.RData"))

Below we plot the field on top of the map. We also add the point values of the standard deviation of the log standard deviation to the graph so we can plot it on top of the map and the field.

# Define the colors for the windows
lower_color <- "darkred"   
upper_color <- "darkblue"  
# 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
# Plot the field on top of the map
f12 <- graph$plot_function(X = cov, 
                          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 = cov, 
                          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)

# We add the point values of standard deviation of the log standard deviation to the graph so we can plot it
standardized_Y_logstd <- (Y_logstd - mean(Y_logstd))/sd(Y_logstd)
graph$add_observations(data = data.frame(y = standardized_Y_logstd, 
                                         edge_number = PtE[,1], 
                                         distance_on_edge = PtE[,2]), 
                       edge_number = "edge_number",
                       distance_on_edge = "distance_on_edge",
                       data_coords = "PtE",
                       normalized = TRUE, 
                       clear_obs = TRUE)
g12 <- graph$plot(data = "y", 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", 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/stand_log_sigma_3.png"), width = 11.2, height = 5.43, plot = final_plot, dpi = 500)
save(final_plot, file = here("data_files/stand_log_sigma_3.RData"))

# Print the combined plot
knitr::include_graphics(here::here("data_files/stand_log_sigma_3.png"))

Figure 1: Standardized logarithm of the standard deviation of the first part of the data. Points represent the point value at each available location, and edges represent the \(\mathrm{std.cov}(s)\) covariate.

We repeat the same process but now for Y_mean. We first add the observations to the graph.

graph$add_observations(data = data.frame(y = Y_mean, 
                                         edge_number = PtE[,1], 
                                         distance_on_edge = PtE[,2]),
                       edge_number = "edge_number",
                       distance_on_edge = "distance_on_edge",
                       data_coords = "PtE",
                       normalized = TRUE, 
                       clear_obs = TRUE)

We fit the model \(\mathrm{mean}(s_i)\sim N(\beta_0+u(s_i), \sigma_e^2)\).

# Build the model
rspde_model_stat_mean <- rspde.metric_graph(graph,
                                       parameterization = "spde")
# Prepare the data for fitting
data_rspde_bru_stat <- graph_data_rspde(rspde_model_stat_mean,
                                        bru = TRUE)
# Define the component
cmp_stat <- y ~ -1 +
  Intercept(1) +
  field(cbind(.edge_number, .distance_on_edge), model = rspde_model_stat_mean)
# Fit the model
rspde_fit_stat_mean <-
  bru(cmp_stat,
      data = data_rspde_bru_stat[["data"]],
      family = "gaussian",
      options = list(verbose = FALSE)
  )
save(rspde_fit_stat_mean, rspde_model_stat_mean, file = here::here("data_files/rspde_fit_stat_and_model_mean.RData"))

load(here::here("data_files/rspde_fit_stat_and_model_mean.RData"))
# Summarize the results
summary(rspde_fit_stat_mean)

## inlabru version: 2.13.0
## INLA version: 25.11.22
## Components:
## Intercept: main = linear(1), group = exchangeable(1L), replicate = iid(1L), NULL
## field: main = cgeneric(cbind(.edge_number, .distance_on_edge)), group = exchangeable(1L), replicate = iid(1L), NULL
## Observation models:
##   Family: 'gaussian'
##     Tag: <No tag>
##     Data class: 'metric_graph_data', 'data.frame'
##     Response class: 'numeric'
##     Predictor: y ~ .
##     Additive/Linear: TRUE/TRUE
##     Used components: effects[Intercept, field], latent[]
## Time used:
##     Pre = 0.538, Running = 8.84, Post = 0.626, Total = 10 
## Fixed effects:
##             mean    sd 0.025quant 0.5quant 0.975quant   mode kld
## Intercept 50.758 2.945     44.851   50.778     56.538 50.776   0
## 
## Random effects:
##   Name     Model
##     field CGeneric
## 
## Model hyperparameters:
##                                           mean    sd 0.025quant 0.5quant
## Precision for the Gaussian observations  0.020 0.002      0.016    0.019
## Theta1 for field                        -2.326 0.133     -2.603   -2.321
## Theta2 for field                        -0.753 0.261     -1.224   -0.766
## Theta3 for field                         1.209 0.542      0.256    1.175
##                                         0.975quant   mode
## Precision for the Gaussian observations      0.024  0.019
## Theta1 for field                            -2.081 -2.297
## Theta2 for field                            -0.199 -0.830
## Theta3 for field                             2.376  1.010
## 
## Marginal log-Likelihood:  -1193.09 
##  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)')

gets_summary_parameters(rspde_fit_stat_mean, rspde_model_stat_mean)

Below we compute the predictions for the mesh and the data location.

mean_on_mesh_pred <- predict(rspde_fit_stat_mean, newdata = data_prd_list_mesh, ~Intercept + field)
cov_for_mean_to_plot <- mean_on_mesh_pred$mean

mean_on_loc_pred <- predict(rspde_fit_stat_mean, newdata = rename(PtE, .edge_number = E, .distance_on_edge = length), ~Intercept + field)
cov_for_mean <- mean_on_loc_pred$mean
save(cov_for_mean_to_plot, cov_for_mean, file = here::here("data_files/cov_for_mean_to_plot_and_cov_for_mean.RData"))

Below we check that the predicted mean is close to the observed mean.

load(here::here("data_files/Y_mean.RData"))
load(here::here("data_files/cov_for_mean_to_plot_and_cov_for_mean.RData"))
plot(Y_mean, type = "l", col = "darkblue")
lines(cov_for_mean, col = "darkred")

Figure 2: Predicted mean and observed mean.

Below we plot the results as we did before.

# Plot the field on top of the map
f12 <- graph$plot_function(X = cov_for_mean_to_plot, 
                          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 = cov_for_mean_to_plot, 
                          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", 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", 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/mean_better_3.png"), width = 11.2, height = 5.43, plot = final_plot, dpi = 500)

knitr::include_graphics(here::here("data_files/mean_better_3.png"))

Figure 3: Average speeds of the first part of the data. Points represent the point value at each available location, and edges represent the \(\mathrm{mean.cov}(s)\) covariate.

We finally add the remaining half of the replicates to the graph and plot one of them.

df_rep <- lapply(1:nrow(Y), function(i){data.frame(y = Y[i,],
                                                   mean_value = cov_for_mean,
                                                   edge_number = PtE[,1],
                                                   distance_on_edge = PtE[,2],
                                                   repl = i)})
df_rep <- do.call(rbind, df_rep)

graph$add_observations(data = df_rep, 
                       edge_number = "edge_number",
                       distance_on_edge = "distance_on_edge",
                       data_coords = "PtE",
                       normalized = TRUE, 
                       clear_obs = TRUE, 
                       group = "repl")

Below we plot replicate 14.

g12 <- graph$plot(data = "y", 
                 p = p12,
                 group = 1, 
                 vertex_size = 0, 
                 data_size = 1) + 
  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) +
  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", 
                 p = p13,
                 group = 1, 
                 vertex_size = 0, 
                 data_size = 1) + 
  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 = "speed", x = "", y = "") +
  xlim(x_left, x_right) + 
  ylim(y_bottom, y_top) +
  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.png"), width = 11.2, height = 5.43, plot = final_plot, dpi = 500)

knitr::include_graphics(here::here("data_files/replicate14_3.png"))

Figure 4: 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.

save(graph, cov, cov_for_mean_to_plot, data_prd_list_mesh, file = here("data_files/pems_repl1_data.RData"))

1 References

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

We used R version 4.5.2 (R Core Team 2025) and the following R packages: cowplot v. 1.2.0 (Wilke 2025), ggmap v. 4.0.2 (Kahle and Wickham 2013), ggpubr v. 0.6.3 (Kassambara 2026), ggtext v. 0.1.2 (Wilke and Wiernik 2022), glue v. 1.8.0 (Hester and Bryan 2024), here v. 1.0.1 (Müller 2020), htmltools v. 0.5.8.1 (Cheng et al. 2024), INLA v. 25.11.22 (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.13.0 (Yuan et al. 2017; Bachl et al. 2019), knitr v. 1.50 (Xie 2014, 2015, 2025), latex2exp v. 0.9.8 (Meschiari 2026), Matrix v. 1.7.3 (Bates, Maechler, and Jagan 2025), MetricGraph v. 1.5.0.9000 (Bolin, Simas, and Wallin 2023a, 2023b, 2024, 2025; Bolin et al. 2024), OpenStreetMap v. 0.4.1 (Fellows and Stotz 2025), patchwork v. 1.3.1 (Pedersen 2025), plotly v. 4.11.0 (Sievert 2020), plotrix v. 3.8.14 (J 2006), renv v. 1.1.7 (Ushey and Wickham 2026), reshape2 v. 1.4.4 (Wickham 2007), reticulate v. 1.44.1 (Ushey, Allaire, and Tang 2025), rmarkdown v. 2.30 (Xie, Allaire, and Grolemund 2018; Xie, Dervieux, and Riederer 2020; Allaire et al. 2025), rSPDE v. 2.5.2.9000 (Bolin and Kirchner 2020; Bolin and Simas 2023; Bolin, Simas, and Xiong 2024), scales v. 1.4.0 (Wickham, Pedersen, and Seidel 2025), sf v. 1.1.0 (E. Pebesma 2018; E. Pebesma and Bivand 2023), slackr v. 3.4.0 (Kaye et al. 2025), sp v. 2.2.1 (E. J. Pebesma and Bivand 2005; Bivand, Pebesma, and Gomez-Rubio 2013), tidyverse v. 2.0.0 (Wickham et al. 2019), tikzDevice v. 0.12.6 (Sharpsteen and Bracken 2023), viridis v. 0.6.5 (Garnier et al. 2024), 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://doi.org/10.32614/CRAN.package.xaringanExtra.

Allaire, JJ, Yihui Xie, Christophe Dervieux, Jonathan McPherson, Javier Luraschi, Kevin Ushey, Aron Atkins, et al. 2025. 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. 2025. Matrix: Sparse and Dense Matrix Classes and Methods. https://doi.org/10.32614/CRAN.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. MetricGraph: Random Fields on Metric Graphs. https://CRAN.R-project.org/package=MetricGraph.

———. 2023b. “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.

———. 2025. “Markov Properties of Gaussian Random Fields on Compact Metric Graphs.” Bernoulli. https://doi.org/10.48550/arXiv.2304.03190.

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.

Cheng, Joe, Carson Sievert, Barret Schloerke, Winston Chang, Yihui Xie, and Jeff Allen. 2024. htmltools: Tools for HTML. https://github.com/rstudio/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 Jan-Peter Stotz. 2025. OpenStreetMap: Access to Open Street Map Raster Images. https://doi.org/10.32614/CRAN.package.OpenStreetMap.

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

Hester, Jim, and Jennifer Bryan. 2024. glue: Interpreted String Literals. https://glue.tidyverse.org/.

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. 2026. ggpubr: “ggplot2” Based Publication Ready Plots. https://doi.org/10.32614/CRAN.package.ggpubr.

Kaye, Matt, Bob Rudis, Andrie de Vries, and Jonathan Sidi. 2025. slackr: Send Messages, Images, r Objects and Files to “Slack” Channels/Users. https://github.com/mrkaye97/slackr.

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.

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. 2026. Latex2exp: Use LaTeX Expressions in Plots. https://doi.org/10.32614/CRAN.package.latex2exp.

Müller, Kirill. 2020. here: A Simpler Way to Find Your Files. https://doi.org/10.32614/CRAN.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. 2025. patchwork: The Composer of Plots. https://doi.org/10.32614/CRAN.package.patchwork.

R Core Team. 2025. 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.

Sharpsteen, Charlie, and Cameron Bracken. 2023. tikzDevice: R Graphics Output in LaTeX Format. https://doi.org/10.32614/CRAN.package.tikzDevice.

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

Ushey, Kevin, JJ Allaire, and Yuan Tang. 2025. reticulate: Interface to “Python”. https://doi.org/10.32614/CRAN.package.reticulate.

Ushey, Kevin, and Hadley Wickham. 2026. renv: Project Environments. https://doi.org/10.32614/CRAN.package.renv.

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. 2025. scales: Scale Functions for Visualization. https://scales.r-lib.org.

Wilke, Claus O. 2025. cowplot: Streamlined Plot Theme and Plot Annotations for “ggplot2”. https://doi.org/10.32614/CRAN.package.cowplot.

Wilke, Claus O., and Brenton M. Wiernik. 2022. ggtext: Improved Text Rendering Support for “ggplot2”. https://doi.org/10.32614/CRAN.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/.

———. 2025. 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.