Skip to contents

Spatio-temporal models

David Bolin and Alexandre B. Simas

Created: 2024-10-23. Last modified: 2024-12-20.

Source: vignettes/spacetime.Rmd
spacetime.Rmd

Introduction

The rSPDE package implements the following spatio-temporal model du+γ(κ2+κd/2ρ⋅∇−Δ)αu=dWQ,on T×D d u + \gamma(\kappa^2 + \kappa^{d/2}\rho\cdot \nabla - \Delta)^{\alpha} u = dW_Q, \quad \text{on } T\times D where TT is a temporal interval and DD is a spatial domain which can be an interval, a bounded subset of ℝ2\mathbb{R}^2 or a metric graph. Here κ>0\kappa>0 is a spatial range parameter, ρ\rho is a drift parameter which is in ℝ\mathbb{R} for spatial domains that are intervals or metric graphs, and in ℝ2\mathbb{R}^2 for spatial domains which are bounded subsets of ℝ2\mathbb{R}^2. Further, WQW_Q is a QQ-Wiener process with spatial covariance operator σ2(κ2−Δ)−β\sigma^2(\kappa^2 - \Delta)^{-\beta}, where σ2\sigma^2 is a variance parameter. Thus, the model has two smoothness parameters α\alpha and β\beta which are assumed to be integers. The model is therefore a generalization of the spatio-temporal models introduced in Lindgren et al. (2024), where the generalization is to allow for drift and to allow for metric graphs as spatial domains. The model is implemented using a finite element discretization of the corresponding precision operator σ−2(d+γ(κ2+κd/2ρ⋅∇−Δ)α)(κ2−Δ)β((d+γ(κ2−κd/2ρ⋅∇−Δ)α)) \sigma^{-2}(d + \gamma(\kappa^2 + \kappa^{d/2}\rho\cdot \nabla - \Delta)^{\alpha})(\kappa^2 - \Delta)^{\beta} ((d + \gamma(\kappa^2 - \kappa^{d/2}\rho\cdot \nabla - \Delta)^{\alpha})) in both space and time, similarly to the discretization introduced in Lindgren et al. (2024). This parameterization of the drift term, using ρκd/2\rho \kappa^{d/2} in place of the standard ρ\rho, is chosen to simplify the enforcement of theoretical bounds on the range of ρ\rho, ensuring that the equation remains well-posed and also providing numerical stability for finite-dimensional approximations.

Implementation details

Let us begin by loading some packages needed for making the plots

The function spacetime.operators() can be used to define the model. The function requires specifying the two smoothness parameters, and the discretization points for the spatial and temporal discretizations. The spatial discretization can be specified through a mesh object from the fmesher package, as a graph from the MetricGraph package, or as the mesh nodes for models on intervals. The temporal discretization can be specified either by specifying the mesh nodes or by providing a mesh object.

Assume that we want to define a model on the spatial interval [0,20][0,20] and the temporal domain [0,10][0,10]. We can then simply specify the mesh nodes as

s <- seq(from = 0, to = 20, length.out = 101)
t <- seq(from = 0, to = 10, length.out = 21)

We can now use spacetime.operators() to construct the model

kappa <- 5
sigma <- 10
gamma <- 1/20
rho_1d <- 1
alpha <- 1
beta <- 1

op <- spacetime.operators(space_loc = s, time_loc = t,
                          kappa = kappa, sigma = sigma, alpha = alpha,
                          beta = beta, rho = rho_1d, gamma = gamma)

The spacetime.operators object has a plot_covariances method which for univariate spatial domains simply plots the covariance C(u(s,t),u(s0,t0))C(u(s,t), u(s_0, t_0)) for a fixed spatio-temporal location (s0,t0)(s_0, t_0) specified by the indices in the spatial and temporal discretizations. For example:

op$plot_covariances(t.ind = 15, s.ind = 50)

The object op contains the matrices needed for evaluating the model, and we have here initialized it by providing values for all parameters.

We can simulate from the model using simulate():

u <- simulate(op)

There is also built-in support for kriging prediction. To illustrate this, we use the simulation to create some noisy observations of the process. For this, we first construct the observation matrix linking the FEM basis functions to the locations where we want to simulate. We first randomly generate some observation locations and then construct the matrix.

n.obs <- 500

obs.loc <- data.frame(x = max(s)*runif(n.obs), t = max(t)*runif(n.obs))
A <- op$make_A(obs.loc$x, obs.loc$t)

We now generate the observations as Yi=u(si)+εiY_i = u(s_i) + \varepsilon_i, where εi∼N(0,σe2)\varepsilon_i \sim N(0,\sigma_e^2) is Gaussian measurement noise.

x <- simulate(op, nsim = 1) 
sigma.e <- 0.01
Y <- as.vector(A%*%x + sigma.e*rnorm(n.obs))

Finally, we compute the kriging prediction of the process uu at the locations in s based on these observations. To specify which locations that should be predicted, the argument Aprd is used. This argument should be an observation matrix that links the mesh locations to the prediction locations.

Aprd <- op$make_A(rep(s, length(t)), rep(t, each = length(s)))
u.krig <- predict(op, A = A, Aprd = Aprd, Y = Y, sigma.e = sigma.e)

The process simulation, and the kriging prediction are shown in the following figure.

data.df <- data.frame(space = obs.loc$x, time = obs.loc$t, field = Y, type = "Data")
krig.df <- data.frame(space = rep(s, length(t)), time = rep(t, each = length(s)),
                      field = as.vector(u.krig$mean), type = "Prediction")
df_plot <- rbind(data.df, krig.df)
ggplot(df_plot) + aes(x = space, y = time, fill = field) +
        facet_wrap(~type) + geom_raster(data = krig.df) +
        geom_point(data = data.df, aes(colour = field),
        show.legend = FALSE) + 
        scale_fill_viridis() + scale_colour_viridis()

Parameter estimation

To estimate the model parameters based on this data, we can use the our rspde_lme function or our inlabru implementation. For this, we collect the data in a data frame, that also contains the spatial locations and the time points.

df_1d <- data.frame(y = as.matrix(Y), space = obs.loc$x, time = obs.loc$t)

rspde_lme implementation

We now fit the model:

res <- rspde_lme(y ~ 1, loc = "space", loc_time = "time", data = df_1d, model = op, parallel = TRUE)

In the call, y~1 indicates that we also want to estimate a mean value of the model, and the arguments loc and loc_time provide the names of the spatial and temporal coordinates in the data frame. Let us see a summary of the fitted model:

summary(res)
#> 
#> Latent model - Spatio-temporal with alpha =  1 , beta =  1
#> 
#> Call:
#> rspde_lme(formula = y ~ 1, loc = "space", loc_time = "time", 
#>     data = df_1d, model = op, parallel = TRUE)
#> 
#> Fixed effects:
#>             Estimate Std.error z-value Pr(>|z|)
#> (Intercept)   0.0419    0.1137   0.369    0.712
#> 
#> Random effects:
#>       Estimate Std.error z-value
#> kappa 4.678250  0.388209  12.051
#> sigma 9.117339  0.875866  10.410
#> gamma 0.054523  0.009569   5.698
#> rho   1.164626  0.042404  27.465
#> 
#> Measurement error:
#>          Estimate Std.error z-value
#> std. dev  0.01704   0.01306   1.304
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
#> 
#> Log-Likelihood:  -568.7374 
#> Number of function calls by 'optim' = 39
#> Optimization method used in 'optim' = L-BFGS-B
#> 
#> Time used to:     fit the model =  18.2336 secs 
#>   set up the parallelization = 4.70659 secs

Let us compare the estimated results with the true values:

results <- data.frame(kappa = c(kappa, res$coeff$random_effects[1]), 
                      sigma = c(sigma, res$coeff$random_effects[2]),
                      gamma = c(gamma, res$coeff$random_effects[3]), 
                      rho = c(rho_1d, res$coeff$random_effects[4]),
                      sigma.e = c(sigma.e, res$coeff$measurement_error),
                      intercept = c(0, res$coeff$fixed_effects),
                      row.names = c("True", "Estimate"))
                      
print(results)
#>            kappa     sigma     gamma      rho    sigma.e  intercept
#> True     5.00000 10.000000 0.0500000 1.000000 0.01000000 0.00000000
#> Estimate 4.67825  9.117339 0.0545233 1.164626 0.01703642 0.04189725

Finally, we can also do prediction based on the fitted model as

pred.data <- data.frame(x = rep(s,length(t)), t = rep(t,each=length(s)))
pred <- predict(res, newdata = pred.data, loc = "x", time = "t")
data.df <- data.frame(space = obs.loc$x, time = obs.loc$t, field = Y, type = "Data")
krig.df <- data.frame(space = rep(s, length(t)), time = rep(t, each = length(s)),
                      field = as.vector(pred$mean), type = "Prediction")
df_plot <- rbind(data.df, krig.df)
ggplot(df_plot) + aes(x = space, y = time, fill = field) +
        facet_wrap(~type) + geom_raster(data = krig.df) +
        geom_point(data = data.df, aes(colour = field),
        show.legend = FALSE) + 
        scale_fill_viridis() + scale_colour_viridis()

inlabru Implementation

Let us now fit the model using our inlabru implementation. We start by creating the model object with the rspde.spacetime() function:

library(inlabru)
#> Loading required package: fmesher
st_bru <- rspde.spacetime(space_loc = s, time_loc = t, alpha = 1, beta=1)

We now create the model component, which requires the user to pass the index as a list containing the elements space with the spatial indices and time with the temporal indices:

cmp <- y ~ -1 + Intercept(1) + field(list(space=space, time = time), model = st_bru)

We are now in a position to fit the model:

bru_fit <- bru(cmp, data = df_1d, options = list(num.threads = "1:1"))

Let us now compare the estimated results (the means of the parameters) with the true values:

param_st_1d <- transform_parameters_spacetime(
                        bru_fit$summary.hyperpar$mean[2:5], 
                                            st_bru)

results <- data.frame(kappa = c(kappa,  param_st_1d[[1]]), 
                      sigma = c(sigma,  param_st_1d[[2]]),
                      gamma = c(gamma,  param_st_1d[[3]]), 
                      rho = c(rho_1d, param_st_1d[[4]]),
                      sigma.e = c(sigma.e, sqrt(1/bru_fit$summary.hyperpar$mean[1])),
                      intercept = c(0, bru_fit$summary.fixed$mean),
                      row.names = c("True", "Estimate"))
                      
print(results)
#>             kappa     sigma      gamma      rho     sigma.e  intercept
#> True     5.000000 10.000000 0.05000000 1.000000 0.010000000 0.00000000
#> Estimate 4.688129  9.109715 0.05301967 1.184015 0.006417642 0.01322333

A spatial example

Let us now illustrate how to implement a spatial version. We start by creating a region of interest and a spatial mesh using the fmesher package:

library(fmesher)
n_loc <- 1000
loc_2d_mesh <- matrix(runif(n_loc * 2), n_loc, 2)
mesh_2d <- fm_mesh_2d(
  loc = loc_2d_mesh,
  cutoff = 0.07,
  max.edge = c(0.2, 0.6)
)
plot(mesh_2d, main = "")

We now proceed as previously by defining a temporal region and the model

t <- seq(from = 0, to = 10, length.out = 11)
kappa <- 10
sigma <- 10
gamma <- 1/200
rho_2d <- 0.3*c(1,1)
alpha <- 1
beta <- 1

op <- spacetime.operators(mesh_space = mesh_2d, time_loc = t,
                          kappa = kappa, sigma = sigma, alpha = alpha,
                          beta = beta, rho = rho_2d, gamma = gamma)
op$plot_covariances(s.ind = 100, t.ind = 5, t.shift = c(-2,0,2))

#> TableGrob (2 x 1) "arrange": 2 grobs
#>   z     cells    name           grob
#> 1 1 (1-1,1-1) arrange gtable[layout]
#> 2 2 (2-2,1-1) arrange gtable[layout]
#> TableGrob (2 x 1) "arrange": 2 grobs
#>   z     cells    name           grob
#> 1 1 (1-1,1-1) arrange gtable[layout]
#> 2 2 (2-2,1-1) arrange gtable[layout]

The spacetime.operators object has a plot_covariances method which can be used to vizualise marginal spatial and temporal covariances. The function takes as input t.ind which is the index of the location in the time discretization to plot the marginal spatial covariance for, and an input s.ind which is the index of the location in the space discretization to show the marginal temporal covariance for. It further takes an input t.shift which can be used to plot covariances C(u(⋅,ti),u(⋅,tj))C(u(\cdot, t_i), u(\cdot, t_j)), where tit_i is t[t.ind] and tjt_j is t[t.ind + t.shift]. For example

We can simulate from the model using simulate():

u <- simulate(op)

Let us plot the simulation for a few time points

proj <- fm_evaluator(mesh_2d, dims = c(100, 100))
U <- matrix(u, nrow = mesh_2d$n, ncol = length(t))
field1 <- fm_evaluate(proj, field = as.vector(U[,2]))
field2 <- fm_evaluate(proj, field = as.vector(U[,3]))
field3 <- fm_evaluate(proj, field = as.vector(U[,4]))

field1.df <- data.frame(x1 = proj$lattice$loc[,1], x2 = proj$lattice$loc[,2], 
                       u = as.vector(field1), type = "u1")
field2.df <- data.frame(x1 = proj$lattice$loc[,1], x2 = proj$lattice$loc[,2], 
                       u = as.vector(field2), type = "u2")
field3.df <- data.frame(x1 = proj$lattice$loc[,1], x2 = proj$lattice$loc[,2], 
                       u = as.vector(field3), type = "u3")
field.df <- rbind(field1.df, field2.df, field3.df)
ggplot(field.df) + aes(x = x1, y = x2, fill = u) + facet_wrap(~type) + geom_raster() + xlim(0,1) + ylim(0,1) + scale_fill_viridis()
#> Warning: Removed 18468 rows containing missing values or values outside the scale range
#> (`geom_raster()`).

We now generate the observations as Yi=u(si)+εiY_i = u(s_i) + \varepsilon_i, where εi∼N(0,σe2)\varepsilon_i \sim N(0,\sigma_e^2) is Gaussian measurement noise.

n.obs <- 500
obs.loc <- data.frame(x = runif(n.obs),y = runif(n.obs), t = max(t)*runif(n.obs))
A <- op$make_A(cbind(obs.loc$x,obs.loc$y), obs.loc$t)
sigma.e <- 0.01
Y <- as.vector(A%*%u + sigma.e*rnorm(n.obs))

To estimate the model parameters based on this data, we can use the rspde_lme function. For this, we collect the data in a data frame, that also contanis the spatial locations, and we fit the model:

df_2d <- data.frame(Y = as.matrix(Y), x = obs.loc$x, y = obs.loc$y, t = obs.loc$t)
res <- rspde_lme(Y ~ 1, loc = c("x", "y"), loc_time = "t", data = df_2d, model = op, parallel = TRUE)
#> Warning in rspde_lme(Y ~ 1, loc = c("x", "y"), loc_time = "t", data = df_2d, :
#> The optimization failed to provide a numerically positive-definite Hessian. You
#> can try to obtain a positive-definite Hessian by setting 'improve_hessian' to
#> TRUE or by setting 'parallel' to FALSE, which allows other optimization methods
#> to be used.
#> Warning in sqrt(diag(inv_fisher)): NaNs produced

Let us see a summary of the fitted model:

summary(res)
#> 
#> Latent model - Spatio-temporal with alpha =  1 , beta =  1
#> 
#> Call:
#> rspde_lme(formula = Y ~ 1, loc = c("x", "y"), loc_time = "t", 
#>     data = df_2d, model = op, parallel = TRUE)
#> 
#> Fixed effects:
#>             Estimate Std.error z-value Pr(>|z|)
#> (Intercept)   0.1502    0.3159   0.476    0.634
#> 
#> Random effects:
#>        Estimate Std.error z-value
#> kappa 11.178757  1.059719  10.549
#> sigma 11.419295  1.408109   8.110
#> gamma  0.005183  0.001050   4.934
#> rho    0.164721  0.949677   0.173
#> NA     0.190030  1.140014   0.167
#> 
#> Measurement error:
#>           Estimate Std.error z-value
#> std. dev 0.0002992       NaN     NaN
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
#> 
#> Log-Likelihood:  -644.4856 
#> Number of function calls by 'optim' = 78
#> Optimization method used in 'optim' = L-BFGS-B
#> 
#> Time used to:     fit the model =  2.94761 mins 
#>   set up the parallelization = 5.71778 secs

Let us compare the estimated results with the true values:

results <- data.frame(kappa = c(kappa, res$coeff$random_effects[1]), 
                      sigma = c(sigma, res$coeff$random_effects[2]),
                      gamma = c(gamma, res$coeff$random_effects[3]), 
                      rho_1 = c(rho_2d[1], res$coeff$random_effects[4]),
                      rho_2 = c(rho_2d[2], res$coeff$random_effects[5]),
                      sigma.e = c(sigma.e, res$coeff$measurement_error),
                      intercept = c(0, res$coeff$fixed_effects),
                      row.names = c("True", "Estimate"))
                      
print(results)
#>             kappa    sigma       gamma     rho_1   rho_2      sigma.e intercept
#> True     10.00000 10.00000 0.005000000 0.3000000 0.30000 0.0100000000 0.0000000
#> Estimate 11.17876 11.41929 0.005182506 0.1647205 0.19003 0.0002992413 0.1502212

Let us now use our inlabru implementation. We first define the model. To this end we need to create a mesh for the time object:

mesh_time <- fmesher::fm_mesh_1d(t)
st_bru_field <- rspde.spacetime(mesh_space = mesh_2d,
        mesh_time = mesh_time, alpha = alpha, beta = beta)

Let us now prepare the data.frame object:

df_bru <- data.frame(y = Y, coord_x = obs.loc$x, coord_y = obs.loc$y, time = obs.loc$t)

Now, the component:

cmp <- y ~ -1 + Intercept(1) + field(list(space=cbind(coord_x, coord_y), 
                time = time), model = st_bru_field)

We are now in a position to fit the model:

bru_fit_field <- bru(cmp, data = df_bru, options = list(num.threads = "1:1"))

Let us now compare the estimated results (the means of the parameters) with the true values:

param_st <- transform_parameters_spacetime(bru_fit_field$summary.hyperpar$mean[2:6], 
                                            st_bru_field)

results <- data.frame(kappa = c(kappa,  param_st[[1]]), 
                      sigma = c(sigma,  param_st[[2]]),
                      gamma = c(gamma,  param_st[[3]]), 
                      rho_1 = c(rho_2d[1], param_st[[4]]),
                      rho_2 = c(rho_2d[2], param_st[[5]]),
                      sigma.e = c(sigma.e, sqrt(1/bru_fit_field$summary.hyperpar$mean[1])),
                      intercept = c(0, bru_fit_field$summary.fixed$mean),
                      row.names = c("True", "Estimate"))
                      
print(results)
#>             kappa    sigma       gamma     rho_1     rho_2    sigma.e intercept
#> True     10.00000 10.00000 0.005000000 0.3000000 0.3000000 0.01000000 0.0000000
#> Estimate 11.31637 11.59464 0.005199523 0.2462587 0.2686612 0.01128791 0.6125935

Fit with bounded_rho = FALSE

In cases where the estimated value of ρ\rho approaches the upper bound (available in st_bru$bound_rho), the model can be re-fitted with bounded_rho = FALSE. This removes the bounding constraint on ρ\rho, which can lead to a better fit but may introduce numerical instability. Below, we demonstrate this for both the 1D and spatial examples.

1D Example with bounded_rho = FALSE

Let us first start with the rspde_lme() example:

op_1d_unbounded <- spacetime.operators(space_loc = s, time_loc = t,
                          kappa = kappa, sigma = sigma, alpha = alpha,
                          beta = beta, rho = rho_1d, gamma = gamma,
                          bounded_rho = FALSE)

res_unbounded <- rspde_lme(y ~ 1, loc = "space", loc_time = "time", 
                      data = df_1d, model = op_1d_unbounded, parallel = TRUE)

results <- data.frame(kappa = c(kappa, res_unbounded$coeff$random_effects[1]), 
                      sigma = c(sigma, res_unbounded$coeff$random_effects[2]),
                      gamma = c(gamma, res_unbounded$coeff$random_effects[3]), 
                      rho = c(rho_1d, res_unbounded$coeff$random_effects[4]),
                      sigma.e = c(sigma.e, res_unbounded$coeff$measurement_error),
                      intercept = c(0, res_unbounded$coeff$fixed_effects),
                      row.names = c("True", "Estimate"))
                      
print(results)
#>              kappa    sigma      gamma      rho   sigma.e    intercept
#> True     10.000000 10.00000 0.00500000 1.000000 0.0100000  0.000000000
#> Estimate  5.939571 13.06381 0.04256791 1.521502 0.3809751 -0.005708659

Now, the inlabru implementation:

### Fit with bounded_rho = FALSE
st_bru_unbounded <- rspde.spacetime(space_loc = s, time_loc = t, alpha = 1, 
            beta = 1, bounded_rho = FALSE)

### Fitting the Model
cmp_unbounded <- y ~ -1 + Intercept(1) + field(list(space = space, time = time), 
                                          model = st_bru_unbounded)

bru_fit_unbounded <- bru(cmp_unbounded, data = df_1d, options = list(num.threads = "1:1", verbose=TRUE))

### Extract and Compare Results
param_unbounded <- transform_parameters_spacetime(
                bru_fit_unbounded$summary.hyperpar$mean[2:5], 
                st_bru_unbounded)

results_unbounded <- data.frame(
  kappa = c(kappa, param_unbounded$kappa),
  sigma = c(sigma, param_unbounded$sigma),
  gamma = c(gamma, param_unbounded$gamma),
  rho = c(rho_1d, param_unbounded$rho),
  sigma.e = c(sigma.e, sqrt(1 / bru_fit_unbounded$summary.hyperpar$mean[1])),
  intercept = c(0, bru_fit_unbounded$summary.fixed$mean),
  row.names = c("True", "Estimate")
)

print(results_unbounded)
#>              kappa   sigma      gamma      rho   sigma.e    intercept
#> True     10.000000 10.0000 0.00500000 1.000000 0.0100000  0.000000000
#> Estimate  6.007034 13.1311 0.03956561 1.743888 0.3744443 -0.005926273

References

Lindgren, Finn, Haakon Bakka, David Bolin, Elias Krainski, and Håvard Rue. 2024. “A Diffusion-Based Spatio-Temporal Extension of Gaussian Matérn Fields.” SORT Statistics and Operations Research Transactions 48 (1): 3–66.