Type:	Package
Title:	Spatial and Spatio-Temporal Models using 'INLA'
Version:	0.1.12
Description:	Prepare objects to implement models over spatial and spacetime domains with the 'INLA' package (https://www.r-inla.org). These objects contain data to for the 'cgeneric' interface in 'INLA', enabling fast parallel computations. We implemented the spatial barrier model, see Bakka et. al. (2019) <doi:10.1016/j.spasta.2019.01.002>, and some of the spatio-temporal models proposed in Lindgren et. al. (2023) https://www.idescat.cat/sort/sort481/48.1.1.Lindgren-etal.pdf. Details are provided in the available vignettes and from the URL bellow.
URL:	https://github.com/eliaskrainski/INLAspacetime
Additional_repositories:	https://inla.r-inla-download.org/R/testing
BugReports:	https://github.com/eliaskrainski/INLAspacetime/issues
License:	GPL-2 | GPL-3 [expanded from: GPL (≥ 2)]
Encoding:	UTF-8
RoxygenNote:	7.3.2
NeedsCompilation:	yes
Depends:	R (≥ 4.3), Matrix, fmesher
Imports:	graphics, grDevices, methods, stats, utils, sp, sf, terra
Suggests:	INLA (≥ 24.02.09), inlabru (≥ 2.10.1), knitr, ggplot2, rmarkdown, parallel, data.table, rnaturalearth
VignetteBuilder:	knitr
BuildVignettes:	true
Packaged:	2025-03-14 10:16:47 UTC; eliask
Author:	Elias Teixeira Krainski [cre, aut, cph], Finn Lindgren [aut], Haavard Rue [aut]
Maintainer:	Elias Teixeira Krainski <eliaskrainski@gmail.com>
Repository:	CRAN
Date/Publication:	2025-03-14 12:10:02 UTC

Function to define the boundary Earch polygon in longlat projection for a given resolution.

Description

Function to define the boundary Earch polygon in longlat projection for a given resolution.

Usage

Earth_poly(resol = 300, crs = "+proj=moll +units=km")

Arguments

Value

a 'st_sfc' object with the Earth polygon.

Internal util functions for polygon properties.

Description

This computes the area of a triangle given its three coordinates.

Usage

Heron(x, y)

Area(x, y)

s2trArea(tr, R = 1)

flatArea(tr)

Stiffness(tr)

Arguments

Details

Function used internally to compute the area of a triangle.

Value

the area of a 2d triangle

the area of a 2d polygon

the area of a triangle in S2

the area of a triangle

the stiffness matrix for a triangle

Warning

Internal functions, not exported.

Spatial and Spatio-Temporal Models using INLA

Description

This package main purpose is to provide user friendly functions to fit temporal, spatial and space-time models using the INLA software available at www.r-inla.org as well the inlabru package available

Usage

INLAspacetime()

Value

opens the Vignettes directory on a browser

The 2nd order temporal matrices with boundary correction

Description

The 2nd order temporal matrices with boundary correction

Usage

Jmatrices(tmesh)

Arguments

Details

Temporal GMRF representation with stationary boundary conditions as in Appendix E of the paper.

Value

return a list of temporal finite element method matrices for the supplied mesh.

Illustrative code to compute the covariance of the second order autoregression (AR2) model.

Description

Computes the auto-covariance for given coefficients.

Usage

ar2cov(a1, a2, k = 30, useC = FALSE)

Arguments

Value

the autocorrelation as a vector or matrix, whenever a1 or a2 are scalar or vector.

Details

Let the second order auto-regression model defined as ⁠x_t + a_1 x_{t-1} + a_2 x_{t-2} = w_t⁠ where w_t ~ N(0, 1).

See Also

ar2precision

Examples

ar2cov(c(-1.7, -1.8), 0.963, k = 5)
plot(ar2cov(-1.7, 0.963), type = "o")

Precision matrix for an AR2 model.

Description

Creates a precision matrix as a sparse matrix object considering the specification stated in Details.

Usage

ar2precision(n, a)

Arguments

Value

the precision matrix as a sparse matrix object with edge correction.

Details

Let the second order auto-regression model be defined as

a_0 x_t + a_1 x_{t-1} + a_2 x_{t-2} = w_t, w_t ~ N(0, 1).

The n times n symmetric precision matrix Q for x_1, x_2, ..., x_n has the following non-zero elements:

Q_{1,1} = Q_{n,n} = a_0^2

Q_{2,2} = Q_{n-1,n-1} = a_0^2 + a_1^2

Q_{1,2} = Q_{2,1} = Q_{n-1,n} = Q_{n,n-1} = a_0 a_1

Q_{t,t} = q_0 = a_0^2 + a_1^2 + a_2^2, t = 3, 4, ..., n-2

Q_{t,t-1} = Q_{t-1,t} = q_1 = a_1(a_0 + a_2), t = 3, 4, ..., n-1

Q_{t,t-2} = Q_{t-2,t} = q_2 = a_2 a_0, t = 3, 4, ..., n

See Also

ar2cov

Examples

ar2precision(7, c(1, -1.5, 0.9))

Define a spacetime model object for the f() call.

Description

Define a spacetime model object for the f() call.

Usage

barrierModel.define(
  mesh,
  barrier.triangles,
  prior.range,
  prior.sigma,
  range.fraction = 0.1,
  constr = FALSE,
  debug = FALSE,
  useINLAprecomp = TRUE,
  libpath = NULL
)

Arguments

Details

See the paper.

Value

objects to be used in the f() formula term in INLA.

Mapper object for automatic inlabru interface

Description

Return an inlabru bru_mapper object that can be used for computing model matrices for the space-time model components. The bru_get_mapper() function is called by the inlabru methods to automatically obtain the needed mapper object (from inlabru ⁠2.7.0.9001⁠; before that, use mapper = bru_get_mapper(model) explicitly).

Usage

bru_get_mapper.stModel_cgeneric(model, ...)

Arguments

Value

A bru_mapper object of class bru_mapper_multi with sub-mappers space and time based on the model smesh and tmesh or mesh objects.

See Also

inlabru::bru_get_mapper()

Computes the Whittle-Matern correlation function.

Description

This computes the correlation function as derived in Matern model, see Matern (1960) eq. (2.4.7). For nu=1, see Whittle (1954) eq. (68). For the limiting case of nu=0, see Besag (1981) eq. (14-15).

Usage

cWhittleMatern(x, range, nu, kappa = sqrt(8 * nu)/range)

Arguments

Value

the correlation.

Details

Whittle, P. (1954) On Stationary Processes in the Plane. Biometrika, Vol. 41, No. 3/4, pp. 434-449. http://www.jstor.org/stable/2332724

Matern, B. (1960) Spatial Variation: Stochastic models and their application to some problems in forest surveys and other sampling investigations. PhD Thesis.

Besag, J. (1981) On a System of Two-Dimensional Recurrence Equations. JRSS-B, Vol. 43 No. 3, pp. 302-309. https://www.jstor.org/stable/2984940

Examples

plot(function(x) cWhittleMatern(x, 1, 5),
  bty = "n", las = 1,
  xlab = "Distance", ylab = "Correlation"
)
plot(function(x) cWhittleMatern(x, 1, 1), add = TRUE, lty = 2)
plot(function(x) cWhittleMatern(x, 1, 0.5), add = TRUE, lty = 3)
abline(h = 0.139, lty = 3, col = gray(0.5, 0.5))

Define the stationary SPDE cgeneric model for INLA.

Description

Define the stationary SPDE cgeneric model for INLA.

Usage

cgeneric_sspde(
  mesh,
  alpha,
  control.priors,
  constr = FALSE,
  debug = FALSE,
  useINLAprecomp = TRUE,
  libpath = NULL
)

Arguments

Value

objects to be used in the f() formula term in INLA.

Note

This is the stationary case of INLA::inla.spde2.pcmatern() with slight change on the marginal variance when the domain is the sphere, following Eq. (23) in Lindgren et. al. (2024).

References

Geir-Arne Fuglstad, Daniel Simpson, Finn Lindgren & Håvard Rue (2019). Constructing Priors that Penalize the Complexity of Gaussian Random Fields. Journal of the American Statistical Association, V. 114, Issue 525.

Finn Lindgren, Haakon Bakka, David Bolin, Elias Krainski and Håvard Rue (2024). A diffusion-based spatio-temporal extension of Gaussian Matérn fields. SORT 48 (1), 3-66 <doi: 10.57645/20.8080.02.13>

check package version and load

Description

check package version and load

Usage

check_package_version_and_load(pkg, minimum_version, quietly = FALSE)

Arguments

Note

this is the same code as in the inlabru package

Download files from the NOAA's GHCN daily data

Description

Download files from the NOAA's GHCN daily data

Usage

downloadUtilFiles(data.dir, year = 2022, force = FALSE)

Arguments

Value

a named character vector with the local file names: daily.data, stations.all, elevation.

See Also

ghcndSelect()

Select data from the daily dataset

Description

Select data from the daily dataset

Usage

ghcndSelect(
  gzfile,
  variable = c("TMIN", "TAVG", "TMAX"),
  station = NULL,
  qflag = "",
  verbose = TRUE,
  astype = as.integer
)

Arguments

Value

if more than one variable, it returns an array whose dimentions are days, stations, variables. If one variable, then it returns a matrix whose dimentions are days, stations.

Details

The default selects TMIN, TAVG and TMAX and return it as integer because the original data is also integer with units in 10 Celcius degrees.

Warning

It can take time to execute if, for example, the data.table package is not available.

References

Menne, M., Durre, I., Vose, R., Gleason, B. and Houston, T. (2012) An overview of the global historical climatology network-daily database. Journal of Atmospheric and Oceanic Technology, 897–910.

Extracts the dual of a mesh object.

Description

Extracts the dual of a mesh object.

Usage

mesh.dual(
  mesh,
  returnclass = c("list", "sf", "sv", "SP"),
  mc.cores = getOption("mc.cores", 2L)
)

Arguments

Value

one of the three in 'returnclass'

Illustrative code for building a mesh in 2d domain.

Description

Creates a mesh object. This is just a test code. For efficient, reliable and general code use the fmesher package.

Usage

mesh2d(loc, domain, max.edge, offset, SP = TRUE)

Arguments

Value

a mesh object.

Warning

This is just for illustration purpose and one should consider the efficient function available a the INLA package.

Illustrative code for Finite Element matrices of a mesh in 2d domain.

Description

Illustrative code for Finite Element matrices of a mesh in 2d domain.

Illustrative code for Finite Element matrices when some triangles are in a barrier domain.

Usage

mesh2fem(mesh, order = 2, barrier.triangles = NULL)

mesh2fem.barrier(mesh, barrier.triangles = NULL)

Arguments

Value

a list object containing the FE matrices.

a list object containing the FE matrices for the barrier problem.

Illustrative code to build the projector matrix for SPDE models.

Description

Creates a projector matrix object.

Usage

mesh2projector(
  mesh,
  loc = NULL,
  lattice = NULL,
  xlim = NULL,
  ylim = NULL,
  dims = c(100, 100)
)

Arguments

Value

the projector matrix as a list with sparse matrix object at x$proj$A..

Warning

This is just for illustration purpose and one should consider the efficient functions available in the INLA and inlabru packages, e.g. inlabru::fm_evaluator.

Detect outliers in a time series considering the raw data and a smoothed version of it.

Description

Detect outliers in a time series considering the raw data and a smoothed version of it.

Usage

outDetect(x, weights = NULL, ff = c(7, 7))

Arguments

Value

logical vector indicating if the data is an outlier with attributes as detailed bellow.

• attr(, 'm') is the mean of x.

• attr(, 's') is the standard devation of x.

• attr(, 'ss') is the standard deviation for the smoothed data y_t that is defined as

y_t = \sum_{k=j}^h w_j * (x_{t-j}+x_{t+j})/2

Both s and ss are used to define outliers if

|x_t-m|/s>ff_1 or |x_t-y_t|/ss>ff_2

• attr(, 'xs') the smoothed time series y_t

Functions to help converting from/to user/internal parametrization. The internal parameters are 'gamma_s, 'gamma_t', 'gamma_E' The user parameters are 'r_s', 'r_t', 'sigma'

Description

Functions to help converting from/to user/internal parametrization. The internal parameters are 'gamma_s, 'gamma_t', 'gamma_E' The user parameters are 'r_s', 'r_t', 'sigma'

Convert from user parameters to SPDE parameters

Convert from SPDE parameters to user parameters

Usage

lgsConstant(lg.s, alpha, smanifold)

params2gammas(
  lparams,
  alpha.t,
  alpha.s,
  alpha.e,
  smanifold = "R2",
  verbose = FALSE
)

gammas2params(lgammas, alpha.t, alpha.s, alpha.e, smanifold = "R2")

Arguments

Details

See equation (23) in the paper.

See equations (19), (20) and (21) in the paper.

See equations (19), (20) and (21) in the paper.

Value

the part of sigma from the spatial constant and ⁠\gamma_s⁠.

log(gamma.s, gamma.t, gamma.e)

log(spatial range, temporal range, sigma)

Examples

params2gammas(log(c(1, 1, 1)), 1, 2, 1, "R2")
gammas2params(log(c(0, 0, 0)), 1, 2, 1, "R2")

Illustrative code to build the precision matrix for SPDE kind models.

Description

Creates a precision matrix as a sparse matrix object. For general code look at the functions in the INLA package.

Usage

spde2precision(kappa, fem, alpha)

Arguments

Value

the precision matrix as a sparse matrix object.

Warning

This is just for illustration purpose and one should consider the efficient function available a the INLA package.

Define a spacetime model object for the f() call.

Description

Define a spacetime model object for the f() call.

Usage

stModel.define(
  smesh,
  tmesh,
  model,
  control.priors,
  constr = FALSE,
  debug = FALSE,
  useINLAprecomp = TRUE,
  libpath = NULL
)

Arguments

Details

The matrices of the kronecker product in Theorem 4.1 of Lindgren et. al. (2024) are computed with the stModel.matrices and the parameters are as Eq (19-21), but parametrized in log scale.

Value

objects to be used in the f() formula term in INLA.

References

Finn Lindgren, Haakon Bakka, David Bolin, Elias Krainski and Håvard Rue (2024). A diffusion-based spatio-temporal extension of Gaussian Matérn fields. SORT 48 (1), 3-66 <doi: 10.57645/20.8080.02.13>

Define the spacetime model matrices.

Description

This function computes all the matrices needed to build the precision matrix for spatio-temporal model, as in Lindgren et. al. (2023)

Usage

stModel.matrices(smesh, tmesh, model, constr = FALSE)

Arguments

Details

See the paper for details.

Value

a list containing needed objects for model definition.

1. 'manifold' to spedify the a string with the model identification

2. a length three vector with the constants c1, c2 and c3

3. the vector d

4. the matrix T

5. the model matrices M_1, ..., M_m

Spacetime precision matrix.

Description

To build the the precision matrix for a spacetime model given the temporal and the spatial meshes.

Usage

stModel.precision(smesh, tmesh, model, theta, verbose = FALSE)

Arguments

Value

a (sparse) precision matrix, as in Lindgren et. al. (2023)

To retrieve goodness of fit statistics for a specific model class.

Description

Extracts dic, waic and log-cpo from an output returned by the inla function from the INLA package or by the bru function from the inlabru package, and computes log-po, mse, mae, crps and scrps for a given input. A summary is applied considering the user imputed function, which by default is the mean.

Usage

stats.inla(m, i = NULL, y, fsummarize = mean)

Arguments

Value

A named numeric vector with the extracted statistics.

Details

It assumes Gaussian posterior predictive distributions! Considering the defaults, for n observations, y_i, i = 1, 2, ..., n, we have

. dic

\sum_i d_i/n

where d_i is the dic computed for observation i.

. waic

\sum_i w_i/n

where w_i is the waic computed for observation i.

. lcpo

-\sum_i \log(p_i)/n

where p_i is the cpo computed for observation i.

For the log-po, crps, and scrps scores it assumes a Gaussian predictive distribution for each observation y_i which the following definitions: z_i = (y_i-\mu_i)/\sigma_i, \mu_i is the posterior mean for the linear predictor, \sigma_i = \sqrt{v_i + 1/\tau_y}, \tau_y is the observation posterior mean, v_i is the posterior variance of the linear predictor for y_i.

Then we consider \phi() the density of a standard Gaussian variable and \psi() the corresponding Cumulative Probability Distribution.

. lpo

-\sum_i \log(\phi(z_i))/n

. crps

\sum_i r_i/n

where

r_i=\sigma_i/\sqrt{\pi} - 2\sigma_i\phi(z_i) + (y_i-\mu_i)(1-2\psi(z_i))

. scrps

\sum_i s_i/n

where

s_i=-\log(2\sigma_i/\sqrt{\pi})/2 -\sqrt{\pi}(\phi(z_i)-\sigma_iz_i/2+z_i\psi(z_i))

Warning

All the scores are negatively oriented which means that smaller scores are better.

References

Held, L. and Schrödle, B. and Rue, H. (2009). Posterior and Cross-validatory Predictive Checks: A Comparison of MCMC and INLA. Statistical Modelling and Regression Structures pp 91–110. https://link.springer.com/chapter/10.1007/978-3-7908-2413-1_6.

Bolin, D. and Wallin, J. (2022) Local scale invariance and robustness of proper scoring rules. Statistical Science. doi:10.1214/22-STS864.

To check unusual low/high variance segments

Description

To check unusual low/high variance segments

Usage

stdSubs(x, nsub = 12, fs = 15)

Arguments

Value

logical indicating if any of the st are fs times lower/higher the average of st, where is returned as an attribute:

• attr(, 'st') numeric vector with the standard deviation at each segment of the data.

To visualize time series over space.

Description

To visualize time series over space.

Usage

stlines(
  stdata,
  spatial,
  group = NULL,
  nmax.group = NULL,
  xscale = 1,
  yscale = 1,
  colour = NULL,
  ...
)

stpoints(
  stdata,
  spatial,
  group = NULL,
  nmax.group = NULL,
  xscale = 1,
  yscale = 1,
  colour = NULL,
  ...
)

Arguments

Details

Scaling the times series is needed before drawing it over the map. The area of the bounding box for the spatial object divided by the number of locations is the standard scaling factor. This is further multiplied by the user given xcale and yscale.

Value

add lines to an existing plot

Functions

• stlines(): each time series over the map centered at the location.

• stpoints(): each time series over the map centered at the location.

Warning

if there are too many geographical locations, it will not look good

Prepare a matrix or a list of matrices for use in some 'cgeneric' code.

Description

Define a graph of the union of the supplied matrices and return the row ordered diagonal plus upper triangle after padding with zeroes each one so that all the returned matrices have the same pattern.

Usage

upperPadding(M, relative = FALSE)

Arguments

Details

If relative=FALSE, each columns of 'xx' is the elements of the corresponding matrix after being padded to fill the pattern of the union graph. If relative=TRUE, each element of 'xx' would be a list with a relative index, 'r', for each non-zero elements of each matrix is returned relative to the union graph, the non-lower elements, 'x', of the corresponding matrix, and a vector, 'o', with the number of non-zero elements for each line of each resulting matrix.

Value

If a unique matrix is given, return the upper triangle considering the 'T' representation in the Matrix package. If a list of matrices is given, return a list of two elements: 'graph' and 'xx'. The 'graph' is the union of the graph from each matrix. If relative=FALSE, 'xx' is a matrix with number of column equals the the number of matrices imputed. If relative=TRUE, it is a list of length equal the number of matrices imputed. See details.

Examples

A <- sparseMatrix(
  i = c(1, 1, 2, 3, 3, 5),
  j = c(2, 5, 3, 4, 5, 5),
  x = -(0:5), symmetric = TRUE
)
A
upperPadding(A)
B <- Diagonal(nrow(A), -colSums(A))
list(a = A, a = B)
upperPadding(list(a = A, b = B))
upperPadding(list(a = A, b = B), relative = TRUE)

Helper functions to retrieve the world map, a world polygon, and create grid centers.

Description

Retrieve the map of the countries

Usage

worldMap(
  crs = "+proj=moll +units=km",
  scale = "medium",
  returnclass = c("sf", "sv")
)

Arguments

References

The land and ocean maps are obtained with the 'rnaturalearth' package.

Define a regular grid in 'Mollweide' projection, with units in kilometers.

Description

Define a regular grid in 'Mollweide' projection, with units in kilometers.

Usage

world_grid(size = 50, domain)

Arguments

Value

a 'sf' points object with the centers of a grid set within Earth (and the supplied domain)