Whittle--Matérn fields with general smoothness
David Bolin, Alexandre B. Simas
Created: 2022-11-23. Last modified: 2024-10-22.
Source:vignettes/fem_models.Rmd
fem_models.Rmd
Introduction
In this vignette we will introduce how to fit Whittle–Matérn fields
with general smoothness based on finite element and rational
approximations. The theory for this approach is provided in Bolin et
al. (2023) and Bolin, Simas, and Xiong (2023). For the
implementation, we make use of the rSPDE
package for the rational approximations.
These models are thus implemented using finite element approximations. Such approximations are not needed for integer smoothness parameters, and for the details about the exact models we refer to the vignettes
For details on the construction of metric graphs, see Working with metric graphs
For further details on data manipulation on metric graphs, see Data manipulation on metric graphs
Constructing the graph and the mesh
We begin by loading the rSPDE
and
MetricGraph
packages:
As an example, we consider the following metric graph
edge1 <- rbind(c(0,0),c(1,0))
edge2 <- rbind(c(0,0),c(0,1))
edge3 <- rbind(c(0,1),c(-1,1))
theta <- seq(from=pi,to=3*pi/2,length.out = 20)
edge4 <- cbind(sin(theta),1+ cos(theta))
edges = list(edge1, edge2, edge3, edge4)
graph <- metric_graph$new(edges = edges)
graph$plot()
To construct a FEM approximation of a Whittle–Matérn field with general smoothness, we must first construct a mesh on the graph.
graph$build_mesh(h = 0.01)
graph$plot(mesh=TRUE)
In the command build_mesh
, the argument h
decides the largest spacing between nodes in the mesh. As can be seen in
the plot, the mesh is very coarse, so let’s reduce the value of
h
and rebuild the mesh:
graph$build_mesh(h = 0.01)
We are now ready to specify the model
for the Whittle–Matérn field
.
For this, we use the matern.operators
function from the
rSPDE
package:
sigma <- 1.3
range <- 0.15
nu <- 0.8
rspde.order <- 2
op <- matern.operators(nu = nu, range = range, sigma = sigma,
parameterization = "matern",
m = rspde.order, graph = graph)
As can be seen in the code, we specify
via the practical correlation range
.
Also, the model is not parametrized by
but instead by
.
Here, sigma
denotes the standard deviation of the field and
nu
is the smoothness parameter, which is related to
via the relation
.
The object op
contains the matrices needed for evaluating
the distribution of the stochastic weights in the FEM approximation.
Let us simulate the field at the mesh locations and plot the result:
u <- simulate(op)
graph$plot_function(X = u, plotly = TRUE)
## Warning: The `plotly` argument of `plot()` is deprecated as of MetricGraph 1.3.0.9000.
## ℹ Please use the `type` argument instead.
## ℹ The argument `plotly` was deprecated in favor of the argument `type`.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
If we want to evaluate
at some locations
,
we need to multiply the weights with the FEM basis functions
evaluated at the locations. For this, we can construct the observation
matrix
,
with elements
,
which links the FEM basis functions to the locations. This can be done
by the function fem_basis
in the metric graph object. To
illustrate this, let us simulate some observation locations on the graph
and construct the matrix:
obs.per.edge <- 100
obs.loc <- NULL
for(i in 1:graph$nE) {
obs.loc <- rbind(obs.loc,
cbind(rep(i,obs.per.edge), runif(obs.per.edge)))
}
n.obs <- obs.per.edge*graph$nE
A <- graph$fem_basis(obs.loc)
In the code, we generate observation locations per edge in the graph, drawn at random. It can be noted that we assume that the observation locations are given in the format where denotes the edge of the observation and is the position on the edge, i.e., the relative distance from the first vertex of the edge.
To compute the precision matrix from the covariance-based rational
approximation one can use the precision()
method on object
returned by the matern.operators()
function:
Q <- precision(op)
As an illustration of the model, let us compute the covariance
function between the process at
,
that is, the point at edge 2 and distance on edge 0.1, and all the other
mesh points. To this end, we can use the helper function
cov_function_mesh
that is contained in the op
object:
Using the model for inference
There is built-in support for computing log-likelihood functions and
performing kriging prediction in the rSPDE
package which we
can use for the graph model. To illustrate this, we use the simulation
to create some noisy observations of the process. We generate the
observations as
,
where
is Gaussian measurement noise,
and
are covariates generated the relative positions of the observations on
the graph.
sigma.e <- 0.1
x1 <- obs.loc[,1]
x2 <- obs.loc[,2]
Y <- 1 + 2*x1 - 3*x2 + as.vector(A %*% u + sigma.e * rnorm(n.obs))
Let us now fit the model. To this end we will use the
graph_lme()
function (that, for the finite element models,
acts as a wrapper for the rspde_lme()
function from the
rSPDE
package). To this end, let us now assemble the
data.frame()
with the observations, the observation
locations and the covariates:
df_data <- data.frame(y = Y, edge_number = obs.loc[,1],
distance_on_edge = obs.loc[,2],
x1 = x1, x2 = x2)
Let us now add the data to the graph object and plot it:
graph$add_observations(data = df_data, normalized = TRUE)
## Adding observations...
## list()
graph$plot(data = "y")
We can now fit the model. To this end, we use the
graph_lme()
function and set the model to
'WM
’.
fit <- graph_lme(y ~ x1 + x2, graph = graph, model = "WM")
Let us obtain a summary of the model:
summary(fit)
##
## Latent model - Whittle-Matern
##
## Call:
## graph_lme(formula = y ~ x1 + x2, graph = graph, model = "WM")
##
## Fixed effects:
## Estimate Std.error z-value Pr(>|z|)
## (Intercept) 0.4326 0.7208 0.600 0.54843
## x1 2.1130 0.2021 10.456 < 2e-16 ***
## x2 -2.3626 0.7190 -3.286 0.00102 **
##
## Random effects:
## Estimate Std.error z-value
## alpha 1.254255 0.023625 53.091
## tau 0.057719 0.006288 9.179
## kappa 15.560492 2.454061 6.341
##
## Random effects (Matern parameterization):
## Estimate Std.error z-value
## nu 0.75426 0.02362 31.927
## sigma 1.34735 0.13608 9.901
## range 0.15786 0.02365 6.676
##
## Measurement error:
## Estimate Std.error z-value
## std. dev 0.097240 0.006306 15.42
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Log-Likelihood: -126.9842
## Number of function calls by 'optim' = 183
## Optimization method used in 'optim' = L-BFGS-B
##
## Time used to: fit the model = 33.21914 secs
We can also obtain additional information by using the function
glance()
:
glance(fit)
## # A tibble: 1 × 9
## nobs sigma logLik AIC BIC deviance df.residual model alpha
## <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <chr> <dbl>
## 1 400 0.0972 -127. 268. 296. 254. 393 Covariance-Based M… 1.25
Let us compare the values of the parameters of the latent model with the true ones:
print(data.frame(sigma = c(sigma, fit$matern_coeff$random_effects[2]),
range = c(range, fit$matern_coeff$random_effects[3]),
nu = c(nu, fit$matern_coeff$random_effects[1]),
row.names = c("Truth", "Estimates")))
## sigma range nu
## Truth 1.300000 0.1500000 0.8000000
## Estimates 1.347354 0.1578632 0.7542551
Kriging
Given that we have estimated the parameters, let us compute the kriging predictor of the field given the observations at the mesh nodes.
We will perform kriging with the predict()
method. To
this end, we need to provide a data.frame
containing the
prediction locations, as well as the values of the covariates at the
prediction locations.
df_pred <- data.frame(edge_number = graph$mesh$VtE[,1],
distance_on_edge = graph$mesh$VtE[,2],
x1 = graph$mesh$VtE[,1],
x2 = graph$mesh$VtE[,2])
u.krig <- predict(fit, newdata = df_pred, normalized = TRUE)
The estimate is shown in the following figure
graph$plot_function(as.vector(u.krig$mean))
We can also use the augment()
function to easily plot
the predictions. Let us a build a 3d plot now and add the observed
values on top of the predictions:
Fitting a model with replicates
Let us now illustrate how to simulate a data set with replicates and
then fit a model to such data. To simulate a latent model with
replicates, all we do is set the nsim
argument to the
number of replicates.
n.rep <- 30
u.rep <- simulate(op, nsim = n.rep)
Now, let us generate the observed values :
Note that
is a matrix with 20 columns, each column containing one replicate. We
need to turn y
into a vector and create an auxiliary vector
repl
indexing the replicates of y
:
y_vec <- as.vector(Y.rep)
repl <- rep(1:n.rep, each = n.obs)
df_data_repl <- data.frame(y = y_vec,
edge_number = rep(obs.loc[,1], n.rep),
distance_on_edge = rep(obs.loc[,2], n.rep),
repl = repl)
Let us clear the previous observations and add the new data to the graph:
graph$add_observations(data = df_data_repl, normalized = TRUE,
group = "repl", clear_obs = TRUE)
## Adding observations...
## list()
We can now fit the model in the same way as before by using the
rspde_lme()
function. Note that we can optimize in parallel
by setting parallel
to TRUE
. If we do not
specify which replicate to consider, in the which_repl
argument, all replicates will be considered.
fit_repl <- graph_lme(y ~ -1, graph = graph, model = "WM", parallel = TRUE)
Observe that we have received a warning saying that the Hessian was
not positive-definite, which ended up creating NaN
s for the
standard errors. Indeed, let us see a summary of the fit:
summary(fit_repl)
##
## Latent model - Whittle-Matern
##
## Call:
## graph_lme(formula = y ~ -1, graph = graph, model = "WM", parallel = TRUE)
##
## No fixed effects.
##
## Random effects:
## Estimate Std.error z-value
## alpha 1.28076 NaN NaN
## tau 0.05429 NaN NaN
## kappa 15.47022 0.43245 35.77
##
## Random effects (Matern parameterization):
## Estimate Std.error z-value
## nu 0.780759 NaN NaN
## sigma 1.323295 0.025305 52.29
## range 0.161550 0.004889 33.05
##
## Measurement error:
## Estimate Std.error z-value
## std. dev 0.301126 0.002945 102.2
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Log-Likelihood: -9741.752
## Number of function calls by 'optim' = 86
## Optimization method used in 'optim' = L-BFGS-B
##
## Time used to: fit the model = 47.96698 secs
## set up the parallelization = 2.69565 secs
Let us, then, follow the suggestion from the warning and refit the
model setting improve_hessian
to TRUE
. This
will obtain a more precise estimate of the Hessian, which can possibly
fix this issue:
fit_repl <- graph_lme(y ~ -1, graph = graph, model = "WM",
parallel = TRUE, improve_hessian = TRUE)
We see that we did not receive any warning now, and the Std. errors were computed accordingly:
summary(fit_repl)
##
## Latent model - Whittle-Matern
##
## Call:
## graph_lme(formula = y ~ -1, graph = graph, model = "WM", parallel = TRUE,
## improve_hessian = TRUE)
##
## No fixed effects.
##
## Random effects:
## Estimate Std.error z-value
## alpha 1.280759 0.030016 42.669
## tau 0.054289 0.007189 7.552
## kappa 15.470222 0.901108 17.168
##
## Random effects (Matern parameterization):
## Estimate Std.error z-value
## nu 0.780759 0.030016 26.01
## sigma 1.323295 0.025305 52.29
## range 0.161550 0.004889 33.05
##
## Measurement error:
## Estimate Std.error z-value
## std. dev 0.301126 0.003117 96.62
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Log-Likelihood: -9741.752
## Number of function calls by 'optim' = 86
## Optimization method used in 'optim' = L-BFGS-B
##
## Time used to: fit the model = 43.45105 secs
## compute the Hessian = 10.00034 secs
## set up the parallelization = 2.68766 secs
Let us also take a glance of the fit:
glance(fit_repl)
## # A tibble: 1 × 9
## nobs sigma logLik AIC BIC deviance df.residual model alpha
## <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <chr> <dbl>
## 1 12000 0.301 -9742. 19492. 19521. 19484. 11996 Covariance-Based … 1.28
Let us compare the values of the parameters of the latent model with the true ones:
print(data.frame(sigma = c(sigma, fit_repl$matern_coeff$random_effects[2]),
range = c(range, fit_repl$matern_coeff$random_effects[3]),
nu = c(nu, fit_repl$matern_coeff$random_effects[1]),
row.names = c("Truth", "Estimates")))
## sigma range nu
## Truth 1.300000 0.15000 0.800000
## Estimates 1.323295 0.16155 0.780759
Let us do kriging. We will use the same prediction locations as in the previous example. Let us get prediction for replicate 10, then add the original observations on top of them:
Using the R-INLA implementation
We also have an R-INLA
implementation of the rational
SPDE approach for metric graphs.
We begin by defining the model by using the
rspde.metric_graph()
function. This function contains the
same arguments as the function rspde.matern()
. We refer the
reader to the R-INLA
implementation of the rational SPDE approach vignette for further
details.
We begin by clearing the previous observations and adding the observations (for the case without replicates) to the graph:
graph$clear_observations()
graph$add_observations(data = df_data, normalized = TRUE)
## Adding observations...
## list()
Let us create the model object:
library(INLA)
rspde_model <- rspde.metric_graph(graph)
By default, the order of the rational approximation is 2.
We can now create the auxiliary quantities that will be needed with
the graph_data_rspde()
function:
data_rspde <- graph_data_rspde(rspde_model, name = "field")
The remaining is standard: we create the formula object, the stack
object, and then fit the model by using the inla()
function. So, first we create the formula object:
f.s <- y ~ -1 + Intercept + x1 + x2 + f(field, model = rspde_model)
Now we create the inla.stack
object. To such an end,
observe that data_rspde
contains the dataset as the
data
component, the index as the index
component and the so-called A
matrix as the
basis
component. We will now create the stack using these
components:
stk.dat <- inla.stack(
data = data_rspde[["data"]]["y"], A = list(data_rspde[["basis"]],1), tag = "est",
effects =
list(c(
data_rspde[["index"]],
list(Intercept = 1)), list(x1 = data_rspde[["data"]]["x1"] ,
x2 = data_rspde[["data"]]["x2"])
)
)
Finally, we can fit the model:
rspde_fit <- inla(f.s, data = inla.stack.data(stk.dat),
control.inla = list(int.strategy = "eb"),
control.predictor = list(A = inla.stack.A(stk.dat), compute = TRUE)
)
We can use the same functions as the rspde
fitted models
in inla
. For instance, we can see the results in the
original scale by creating the result
object:
result_fit <- rspde.result(rspde_fit, "field", rspde_model)
summary(result_fit)
## mean sd 0.025quant 0.5quant 0.975quant mode
## std.dev 1.415360 0.1542600 1.1438500 1.403950 1.748880 1.377050
## range 0.150774 0.0329626 0.0980215 0.146573 0.226908 0.138074
## nu 0.871800 0.1277970 0.6227290 0.872179 1.121620 0.876478
Let us compare with the true values:
result_df <- data.frame(
parameter = c("std.dev", "range", "nu"),
true = c(sigma, range, nu),
mean = c(
result_fit$summary.std.dev$mean,
result_fit$summary.range$mean,
result_fit$summary.nu$mean
),
mode = c(
result_fit$summary.std.dev$mode,
result_fit$summary.range$mode,
result_fit$summary.nu$mode
)
)
print(result_df)
## parameter true mean mode
## 1 std.dev 1.30 1.4153576 1.3770544
## 2 range 0.15 0.1507739 0.1380736
## 3 nu 0.80 0.8718003 0.8764784
We can also plot the posterior marginal densities with the help of
the gg_df()
function:
posterior_df_fit <- gg_df(result_fit)
library(ggplot2)
ggplot(posterior_df_fit) + geom_line(aes(x = x, y = y)) +
facet_wrap(~parameter, scales = "free") + labs(y = "Density")
Kriging with the R-INLA
implementation
We will do kriging on the mesh locations:
pred_loc <- graph$mesh$VtE
Let us now add the observations for prediction:
graph$add_observations(data = data.frame(y=rep(NA,nrow(pred_loc)),
x1 = graph$mesh$VtE[,1],
x2 = graph$mesh$VtE[,2],
edge_number = pred_loc[,1],
distance_on_edge = pred_loc[,2]),
normalized = TRUE)
## Adding observations...
## list()
Let us now create a new model and, then, compute the auxiliary
components at the prediction locations. To this end, we set the argument
only_pred
to TRUE
, in which it will return the
data.frame
containing the NA
data.
rspde_model_prd <- rspde.metric_graph(graph)
data_rspde_prd <- graph_data_rspde(rspde_model_prd, only_pred = TRUE)
Let us build the prediction stack using the components of
data_rspde_prd
and gather it with the estimation stack.
ef.prd <-
list(c(data_rspde_prd[["index"]], list(Intercept = 1)),
list(x1 = data_rspde_prd[["data"]][["x1"]],
x2 = data_rspde_prd[["data"]][["x2"]]))
stk.prd <- inla.stack(
data = data.frame(y = data_rspde_prd[["data"]][["y"]]),
A = list(data_rspde_prd[["basis"]],1), tag = "prd",
effects = ef.prd
)
stk.all <- inla.stack(stk.dat, stk.prd)
Let us obtain the predictions:
rspde_fitprd <- inla(f.s,
data = inla.stack.data(stk.all),
control.predictor = list(
A = inla.stack.A(stk.all),
compute = TRUE, link = 1
),
control.compute = list(
return.marginals = FALSE,
return.marginals.predictor = FALSE
),
control.inla = list(int.strategy = "eb")
)
Let us now extract the indices of the predicted nodes and store the means:
id.prd <- inla.stack.index(stk.all, "prd")$data
m.prd <- rspde_fitprd$summary.fitted.values$mean[id.prd]
Finally, let us plot the predicted values. To this end we will use
the plot_function()
graph method.
graph$plot_function(m.prd, plotly = TRUE)
Using R-INLA
implementation to fit models with
replicates
Let us begin by cloning the graph and clearing the observations on the cloned graph:
graph_rep <- graph$clone()
graph_rep$clear_observations()
We will now add the data with replicates to the graph:
graph_rep$add_observations(data = data.frame(y=as.vector(Y.rep),
edge_number = rep(obs.loc[,1], n.rep),
distance_on_edge = rep(obs.loc[,2], n.rep),
repl = rep(1:n.rep, each = n.obs)),
group = "repl",
normalized = TRUE)
## Adding observations...
## list()
Let us create a new rspde
model object:
rspde_model_rep <- rspde.metric_graph(graph_rep)
To fit the model with replicates we need to create the auxiliary
quantities with the graph_data_rspde()
function, where we
set the repl
argument in the function
graph_data_spde
to .all
since we want to use
all replicates:
data_rspde_rep <- graph_data_rspde(rspde_model_rep, name = "field", repl = ".all")
Let us now create the corresponding inla.stack
object:
st.dat.rep <- inla.stack(
data = data_rspde_rep[["data"]],
A = data_rspde_rep[["basis"]],
effects = data_rspde_rep[["index"]]
)
Observe that we need the response variable y
to be a
vector. We can now create the formula
object, remembering
that since we gave the name argument field
, when creating
the index, we need to pass field.repl
to the
formula
:
f.rep <-
y ~ -1 + f(field,
model = rspde_model_rep,
replicate = field.repl
)
We can, finally, fit the model:
rspde_fit_rep <-
inla(f.rep,
data = inla.stack.data(st.dat.rep),
family = "gaussian",
control.predictor =
list(A = inla.stack.A(st.dat.rep))
)
We can obtain the estimates in the original scale with the
rspde.result()
function:
result_fit_rep <- rspde.result(rspde_fit_rep, "field", rspde_model_rep)
summary(result_fit_rep)
## mean sd 0.025quant 0.5quant 0.975quant mode
## std.dev 1.320840 0.02449400 1.273750 1.320500 1.369970 1.319560
## range 0.148385 0.00580176 0.137216 0.148313 0.160003 0.148228
## nu 0.847431 0.03720320 0.776637 0.846494 0.922634 0.843904
Let us compare with the true values of the parameters:
result_rep_df <- data.frame(
parameter = c("std.dev", "range", "nu"),
true = c(sigma, range, nu),
mean = c(
result_fit_rep$summary.std.dev$mean,
result_fit_rep$summary.range$mean,
result_fit_rep$summary.nu$mean
),
mode = c(
result_fit_rep$summary.std.dev$mode,
result_fit_rep$summary.range$mode,
result_fit_rep$summary.nu$mode
)
)
print(result_rep_df)
## parameter true mean mode
## 1 std.dev 1.30 1.3208387 1.3195620
## 2 range 0.15 0.1483851 0.1482285
## 3 nu 0.80 0.8474306 0.8439043
We can also plot the posterior marginal densities with the help of
the gg_df()
function:
Using inlabru
implementation
The inlabru
package allows us to fit models and do
kriging in a straighforward manner, without having to handle
A
matrices, indices nor inla.stack
objects.
Therefore, we suggest the reader to use this implementation when using
our implementation to fit real data.
Let us clear the graph, since it contains NA
observations we used for prediction, add the observations again, and
create a new rSPDE
model object:
graph$clear_observations()
graph$add_observations(data = df_data,
normalized = TRUE)
## Adding observations...
## list()
rspde_model <- rspde.metric_graph(graph)
Let us now load the inlabru
package and create the
component (which is inlabru
’s formula-like object). Since
we are using the data from the graph, inlabru
will also
obtain the locations from the graph, thus, there is no need to provide
the locations. However, we need a name for the locations for using
inlabru’s predict
method. Therefore, we can choose any name
for the location that is not a name being used in the graph’s data. In
our case we will use the name loc
:
Let us now build the auxiliary data to be used with the
graph_data_rspde()
function, where we pass the name of the
location variable in the above formula as the loc_name
argument, which in this case is "loc"
:
data_rspde_bru <- graph_data_rspde(rspde_model, loc_name = "loc")
Now, we can directly fit the model, by using the data
element of data_rspde_bru
:
rspde_bru_fit <-
bru(cmp,
data=data_rspde_bru[["data"]]
)
Let us now obtain the estimates of the parameters in the original
scale by using the rspde.result()
function:
result_bru_fit <- rspde.result(rspde_bru_fit, "field", rspde_model)
summary(result_bru_fit)
## mean sd 0.025quant 0.5quant 0.975quant mode
## std.dev 1.419390 0.1563540 1.1468900 1.406750 1.759890 1.376810
## range 0.150714 0.0330641 0.0973102 0.146699 0.226653 0.138791
## nu 0.869928 0.1302480 0.6201630 0.868495 1.128670 0.866517
Let us compare with the true values of the parameters:
result_bru_df <- data.frame(
parameter = c("std.dev", "range", "nu"),
true = c(sigma, range, nu),
mean = c(
result_bru_fit$summary.std.dev$mean,
result_bru_fit$summary.range$mean,
result_bru_fit$summary.nu$mean
),
mode = c(
result_bru_fit$summary.std.dev$mode,
result_bru_fit$summary.range$mode,
result_bru_fit$summary.nu$mode
)
)
print(result_bru_df)
## parameter true mean mode
## 1 std.dev 1.30 1.4193881 1.3768058
## 2 range 0.15 0.1507138 0.1387908
## 3 nu 0.80 0.8699281 0.8665170
We can also plot the posterior marginal densities with the help of
the gg_df()
function:
posterior_df_bru_fit <- gg_df(result_bru_fit)
ggplot(posterior_df_bru_fit) + geom_line(aes(x = x, y = y)) +
facet_wrap(~parameter, scales = "free") + labs(y = "Density")
Kriging with the inlabru
implementation
It is very easy to do kriging with the inlabru
implementation. We simply need to provide the prediction locations to
the predict()
method.
In this example we will use the mesh locations. To this end we will
use the get_mesh_locations()
method. We also set
bru=TRUE
and loc="loc"
to obtain a data list
suitable to be used with inlabru
.
data_prd_list <- graph$get_mesh_locations(bru = TRUE,
loc = "loc")
data_prd_list[["x1"]] <- data_prd_list$loc[,1]
data_prd_list[["x2"]] <- data_prd_list$loc[,2]
Now, we can simply provide these locations to the
predict
method along with the fitted object
rspde_bru_fit
:
y_pred <- predict(rspde_model, cmp, rspde_bru_fit, newdata=data_prd_list, ~Intercept + x1 + x2 + field)
Finally, let us plot the predicted values. To this end we will use
the plot()
method on y_pred
:
plot(y_pred)
We can also create the 3d plot, together with the true data:
p <- graph$plot(data = "y", plotly=TRUE)
plot(y_pred, plotly = TRUE, p = p)
Using inlabru to fit models with replicates
We can also use our inlabru
implementation to fit models
with replicates. We will consider the same data that was generated
above, where the number of replicates is 30.
For this implementation we will use the rspde_model_rep
object.
We can now create the component, passing the vector with the indices
of the replicates as the replicate
argument. To obtain the
auxiliary data, we will pass repl
argument we use the
function graph_data_rspde()
, where we set it to
.all
, since we want all replicates. Further, we also pass
the loc_name
argument.
data_rspde_rep <- graph_data_rspde(rspde_model_rep, repl = ".all", loc_name = "loc")
We can now define the bru
component formula, passing the
repl
as the replicate
argument:
cmp_rep <-
y ~ -1 + field(loc, model = rspde_model_rep,
replicate = repl)
Now, we are ready to fit the model:
rspde_bru_fit_rep <-
bru(cmp_rep,
data=data_rspde_rep[["data"]],
options=list(
family = "gaussian")
)
We can obtain the estimates in the original scale with the
rspde.result()
function:
result_bru_fit_rep <- rspde.result(rspde_bru_fit_rep, "field", rspde_model_rep)
summary(result_bru_fit_rep)
## mean sd 0.025quant 0.5quant 0.975quant mode
## std.dev 1.320840 0.02449400 1.273750 1.320500 1.369970 1.319560
## range 0.148385 0.00580178 0.137216 0.148313 0.160003 0.148228
## nu 0.847431 0.03720310 0.776637 0.846494 0.922633 0.843904
Let us compare with the true values of the parameters:
result_bru_rep_df <- data.frame(
parameter = c("std.dev", "range", "nu"),
true = c(sigma, range, nu),
mean = c(
result_bru_fit_rep$summary.std.dev$mean,
result_bru_fit_rep$summary.range$mean,
result_bru_fit_rep$summary.nu$mean
),
mode = c(
result_bru_fit_rep$summary.std.dev$mode,
result_bru_fit_rep$summary.range$mode,
result_bru_fit_rep$summary.nu$mode
)
)
print(result_bru_rep_df)
## parameter true mean mode
## 1 std.dev 1.30 1.3208387 1.3195620
## 2 range 0.15 0.1483851 0.1482285
## 3 nu 0.80 0.8474306 0.8439043
We can also plot the posterior marginal densities with the help of
the gg_df()
function:
posterior_df_bru_fit_rep <- gg_df(result_bru_fit_rep)
ggplot(posterior_df_bru_fit_rep) + geom_line(aes(x = x, y = y)) +
facet_wrap(~parameter, scales = "free") + labs(y = "Density")
Let us now do prediction for observations of replicate
10
. We start by building the data list with the prediction
locations:
data_prd_list_repl <- graph$get_mesh_locations(bru = TRUE,
loc = "loc")
data_prd_list_repl[["repl"]] <- rep(10, nrow(data_prd_list$loc))
Let us now obtain predictions for this replicate:
y_pred <- predict(rspde_model_rep, cmp_rep, rspde_bru_fit_rep,
newdata=data_prd_list_repl, ~field_eval(loc, replicate = repl))
We can now plot the predictions along with the observed values for
replicate 10
:
p <- plot(y_pred, plotly=TRUE)
graph_rep$plot(data = "y", group = 10, plotly=TRUE, p = p)
An example with a non-stationary model
Our goal now is to show how one can fit model with non-stationary (std. deviation) and non-stationary (a range parameter). One can also use the parameterization in terms of non-stationary SPDE parameters and .
We follow the same structure as INLA
. However,
INLA
only allows one to specify B.tau
and
B.kappa
matrices, and, in INLA
, if one wants
to parameterize in terms of range and standard deviation one needs to do
it manually. Here we provide the option to directly provide the matrices
B.sigma
and B.range
.
The usage of the matrices B.tau
and B.kappa
are identical to the corresponding ones in
inla.spde2.matern()
function. The matrices
B.sigma
and B.range
work in the same way, but
they parameterize the stardard deviation and range, respectively.
The columns of the B
matrices correspond to the same
parameter. The first column does not have any parameter to be estimated,
it is a constant column.
So, for instance, if one wants to share a parameter with both
sigma
and range
(or with both tau
and kappa
), one simply let the corresponding column to be
nonzero on both B.sigma
and B.range
(or on
B.tau
and B.kappa
).
Creating the graph and adding data
For this example we will consider the pems
data
contained in the MetricGraph
package. The data consists of
traffic speed observations on highways in the city of San Jose,
California. The variable y
contains the traffic speeds.
pems_graph <- metric_graph$new(edges = pems$edges)
pems_graph$add_observations(data = pems$data)
## list()
pems_graph$prune_vertices()
pems_graph$build_mesh(h=0.1)
The summary of this graph:
summary(pems_graph)
## A metric graph object with:
##
## Vertices:
## Total: 347
## Degree 1: 11; Degree 2: 16; Degree 3: 315; Degree 4: 5;
## With incompatible directions: 17
##
## Edges:
## Total: 504
## Lengths:
## Min: 0.01040218 ; Max: 7.677232 ; Total: 470.7559
## Weights:
## Columns: .weights
## That are circles: 0
##
## Graph units:
## Vertices unit: degree ; Lengths unit: km
##
## Longitude and Latitude coordinates: TRUE
## Which spatial package: sp
## CRS: +proj=longlat +datum=WGS84 +no_defs
##
## Some characteristics of the graph:
## Connected: TRUE
## Has loops: FALSE
## Has multiple edges: TRUE
## Is a tree: FALSE
## Distance consistent: FALSE
## Has Euclidean edges: FALSE
##
## Computed quantities inside the graph:
## Laplacian: FALSE ; Geodesic distances: TRUE
## Resistance distances: FALSE ; Finite element matrices: FALSE
##
## Mesh:
## Max h_e: 0.09998277 ; Min n_e: 0
##
## Data:
## Columns: y
## Groups: .group
##
## Tolerances:
## vertex-vertex: 0.001
## vertex-edge: 0.001
## edge-edge: 0
Observe that it is a non-Euclidean graph.
Let us create as non-stationary covariates, the position on the edge, which will capture if the traffic speed was taken close to the intersections. We will make this function symmetric around 0.5 by subtracting 0.5 for points larger than 0.5. That is, the covariate is zero close to intersections.
cov_pos <- 2 * ifelse(pems_graph$mesh$VtE[,2] > 0.5,
1-pems_graph$mesh$VtE[,2],
pems_graph$mesh$VtE[,2])
We will now build the non-stationary matrices to be used:
Let us also obtain the same covariate for the observations:
cov_obs <- pems$data[[".distance_on_edge"]]
cov_obs <- 2 * ifelse(cov_obs > 0.5,
1 - cov_obs,
cov_obs)
Let add this covariate to the data:
pems_graph$add_observations(data = pems_graph$mutate(cov_obs = cov_obs),
clear_obs = TRUE)
## Adding observations...
## The unit for edge lengths is km
## The current tolerance for removing distant observations is (in km): 3.83861599015656
## list()
Fitting the model with graph_lme
We are now in position to fit this model using the
graph_lme()
function. We will also add cov_obs
as a covariate for the model.
fit <- graph_lme(y ~ cov_obs, graph = pems_graph, model = list(type = "WhittleMatern",
B.sigma = B.sigma, B.range = B.range, fem = TRUE))
Let us now obtain a summary of the fitted model:
summary(fit)
##
## Latent model - Generalized Whittle-Matern
##
## Call:
## graph_lme(formula = y ~ cov_obs, graph = pems_graph, model = list(type = "WhittleMatern",
## B.sigma = B.sigma, B.range = B.range, fem = TRUE))
##
## Fixed effects:
## Estimate Std.error z-value Pr(>|z|)
## (Intercept) 47.173 2.779 16.97 < 2e-16 ***
## cov_obs 6.974 1.835 3.80 0.000145 ***
##
## Random effects:
## Estimate Std.error z-value
## alpha 2.0445 0.0297 68.838
## Theta 1 3.7147 0.7165 5.185
## Theta 2 2.8213 0.7764 3.634
## Theta 3 -1.7547 0.9409 -1.865
## Theta 4 -1.4858 1.0276 -1.446
##
## Measurement error:
## Estimate Std.error z-value
## std. dev 7.30787 0.05356 136.4
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Log-Likelihood: -1209.851
## Number of function calls by 'optim' = 501
## Optimization method used in 'optim' = Nelder-Mead
##
## Time used to: fit the model = 37.42936 secs
Let us plot the range parameter along the mesh, so we can see how it is varying:
est_range <- exp(B.range[,-1]%*%fit$coeff$random_effects[2:5])
pems_graph$plot_function(X = est_range, vertex_size = 0,
type = "mapview", mapview_caption = "Range")
## Warning: Found less unique colors (100) than unique zcol values (522)!
## Interpolating color vector to match number of zcol values.
## Warning: Found less unique colors (100) than unique zcol values (522)!
## Interpolating color vector to match number of zcol values.
Similarly, we have for sigma:
est_sigma <- exp(B.sigma[,-1]%*%fit$coeff$random_effects[2:5])
pems_graph$plot_function(X = est_sigma, vertex_size = 0,
type = "mapview", mapview_caption = "Sigma")
## Warning: Found less unique colors (100) than unique zcol values (520)!
## Interpolating color vector to match number of zcol values.
## Warning: Found less unique colors (100) than unique zcol values (520)!
## Interpolating color vector to match number of zcol values.
Our goal now is to plot the estimated marginal standard deviation of
this model. To this end, we start by creating the non-stationary Matérn
operator using the rSPDE
package:
rspde_object_ns <- rSPDE::spde.matern.operators(graph = pems_graph,
parameterization = "matern",
B.sigma = B.sigma,
B.range = B.range,
theta = fit$coeff$random_effects[2:5],
nu = fit$coeff$random_effects[1] - 0.5)
Now, we compute the estimated marginal standard deviation:
est_cov_matrix <- rspde_object_ns$covariance_mesh()
est_std_dev <- sqrt(Matrix::diag(est_cov_matrix))
We can now plot:
pems_graph$plot_function(X = est_std_dev, vertex_size = 0,
type = "mapview", mapview_caption = "Std. dev")
## Warning: Found less unique colors (100) than unique zcol values (4787)!
## Interpolating color vector to match number of zcol values.
## Warning: Found less unique colors (100) than unique zcol values (4787)!
## Interpolating color vector to match number of zcol values.
Fitting the inlabru rSPDE model
Let us then fit the same model using inlabru
now. We
start by defing the rSPDE
model with the
rspde.metric_graph()
function:
rspde_model_nonstat <- rspde.metric_graph(pems_graph,
B.sigma = B.sigma,
B.range = B.range,
parameterization = "matern")
Let us now create the data.frame()
and the vector with
the replicates indexes:
data_rspde_bru_ns <- graph_data_rspde(rspde_model_nonstat, loc_name = "loc")
Let us create the component and fit.
cmp_nonstat <-
y ~ -1 + Intercept(1) + cov_obs + field(loc,
model = rspde_model_nonstat
)
rspde_fit_nonstat <-
bru(cmp_nonstat,
data = data_rspde_bru_ns[["data"]],
family = "gaussian"
)
We can get the summary:
summary(rspde_fit_nonstat)
## inlabru version: 2.11.1
## INLA version: 24.10.13
## Components:
## Intercept: main = linear(1), group = exchangeable(1L), replicate = iid(1L)
## cov_obs: main = linear(cov_obs), group = exchangeable(1L), replicate = iid(1L)
## field: main = cgeneric(loc), group = exchangeable(1L), replicate = iid(1L)
## Likelihoods:
## Family: 'gaussian'
## Data class: 'metric_graph_data', 'list'
## Predictor: y ~ .
## Time used:
## Pre = 0.267, Running = 38.1, Post = 0.61, Total = 39
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## Intercept 49.737 2.910 43.911 49.753 55.465 49.751 0
## cov_obs 2.411 1.751 -1.028 2.411 5.846 2.411 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.020
## Theta1 for field 3.199 0.328 2.534 3.205
## Theta2 for field 2.038 0.274 1.487 2.042
## Theta3 for field 0.026 0.672 -1.244 0.010
## Theta4 for field 0.174 0.467 -0.711 0.163
## Theta5 for field 1.187 0.388 0.489 1.167
## 0.975quant mode
## Precision for the Gaussian observations 0.024 0.020
## Theta1 for field 3.826 3.232
## Theta2 for field 2.565 2.060
## Theta3 for field 1.398 -0.065
## Theta4 for field 1.128 0.113
## Theta5 for field 2.010 1.072
##
## Deviance Information Criterion (DIC) ...............: 2301.61
## Deviance Information Criterion (DIC, saturated) ....: 429.79
## Effective number of parameters .....................: 103.68
##
## Watanabe-Akaike information criterion (WAIC) ...: 2294.96
## Effective number of parameters .................: 79.52
##
## Marginal log-Likelihood: -1241.40
## 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)')
We can obtain outputs with respect to parameters in the original
scale by using the function rspde.result()
:
result_fit_nonstat <- rspde.result(rspde_fit_nonstat, "field", rspde_model_nonstat)
summary(result_fit_nonstat)
## mean sd 0.025quant 0.5quant 0.975quant mode
## Theta1.matern 3.1986600 0.328330 2.533660 3.20494000 3.82607 3.2323900
## Theta2.matern 2.0380500 0.273908 1.486530 2.04220000 2.56488 2.0600800
## Theta3.matern 0.0263043 0.671511 -1.244500 0.00958606 1.39762 -0.0647099
## Theta4.matern 0.1744900 0.467380 -0.711059 0.16317800 1.12801 0.1130010
## nu 1.5184700 0.134683 1.242420 1.52283000 1.76189 1.5275200
We can also plot the posterior densities. To this end we will use the
gg_df()
function, which creates ggplot2
user-friendly data frames: