R: Function for estimation with BRISC

BRISC_estimation {BRISC}

R Documentation

Function for estimation with BRISC

Description

The function BRISC_estimation fits univariate spatial regression models for large spatial data using Vecchia's approximate likelihood (Vecchia, 1988). BRISC_estimation uses the sparse Cholesky representation of Vecchia’s likelihood developed in Datta et al., 2016. The Maximum Likelihood Estimates (MLE) of the parameters are used later for calculating the confidence interval via the BRISC_bootstrap (BRISC, Saha & Datta, 2018).
We recommend using BRISC_estimation followed by BRISC_bootstrap to obtain the confidence intervals for the model parameters.

The optimization is performed with C library of limited-memory BFGS libLBFGS: a library of Limited-memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS),
http://www.chokkan.org/software/liblbfgs/ (Naoaki Okazaki). For user convenience the source codes of the package libLBFGS are provided in the package. The code for the coordinate ordering method, approximate Maximum Minimum Distance (Guinness, 2018) is available in
https://github.com/joeguinness/gp_reorder/tree/master/R and is adopted with minor modification. Some code blocks are borrowed from the R package: spNNGP: Spatial Regression Models for Large Datasets using Nearest Neighbor Gaussian Processes https://CRAN.R-project.org/package=spNNGP .

Usage

BRISC_estimation(coords, y, x = NULL, sigma.sq = 1,
                 tau.sq = 0.1, phi = 1,
                 nu = 1.5, n.neighbors = 15,
                 n_omp = 1, order = "Sum_coords",
                 cov.model = "exponential",
                 search.type = "tree",
                 stabilization = NULL,
                 pred.stabilization = 1e-8,
                 verbose = TRUE, eps = 2e-05,
                 nugget_status = 1, ordering = NULL,
                 neighbor = NULL, tol = 12)

Arguments

`coords`	an `n \times 2` matrix of the observation coordinates in `R^2` (e.g., easting and northing).
`y`	an `n` length vector of response at the observed coordinates.
`x`	an `n \times p` matrix of the covariates in the observation coordinates. Default value is `n \times 1` matrix of `1` to adjust for the mean(intercept).
`sigma.sq`	starting value of sigma square. Default value is 1.
`tau.sq`	starting value of tau square. Default value is 0.1.
`phi`	starting value of phi. Default value is 1.
`nu`	starting value of nu, only required for matern covariance model. Default value is 1.5.
`n.neighbors`	number of neighbors used in the NNGP. Default value is 15.
`n_omp`	number of threads to be used, value can be more than 1 if source code is compiled with OpenMP support. Default is 1.
`order`	keyword that specifies the ordering scheme to be used in ordering the observations. Supported keywords are: `"AMMD"` and `"Sum_coords"` for approximate Maximum Minimum Distance and sum of coordinate based ordering, respectively. Default value is `"Sum_coords"`. `n > 65` is required for `"AMMD"`. Ignored, if `"ordering"` in not `NULL`.
`cov.model`	keyword that specifies the covariance function to be used in modelling the spatial dependence structure among the observations. Supported keywords are: `"exponential"`, `"matern"`, `"spherical"`, and `"gaussian"` for exponential, Matern, spherical and Gaussian covariance function respectively. Default value is "exponential".
`search.type`	keyword that specifies type of nearest neighbor search algorithm to be used. Supported keywords are: `"brute"`, `"tree"` and `"cb"`. `"brute"` and `"tree"` provide the same result, though `"tree"` should be faster. `"cb"` implements fast code book search described in Ra and Kim (1993) modified for NNGP. If locations do not have identical coordinate values on the axis used for the nearest neighbor ordering (see `order` argument) then `"cb"` and `"brute"` should produce identical neighbor sets. However, if there are identical coordinate values on the axis used for nearest neighbor ordering, then `"cb"` and `"brute"` might produce different, but equally valid, neighbor sets, e.g., if data are on a grid. Default value is `"tree"`. Ignored, if `"neighbor"` in not `NULL`.
`stabilization`	when the spatial errors are generated from a very smooth covarince model (lower values of phi for spherical and Gaussian covariance and low phi and high nu for Matern covarinace), the estimation process may fail due to computational instability. If `stabilization = TRUE`, performs stabilization by adding a white noise to the reordered data with nugget `tau.sq =` `sigma.sq` * 1e-06. Estimation is performed on this new data with `nugget_status = 1` (see `nugget_status` argument below). Default value is `TRUE` for `cov.model = "expoenential"` and `FALSE` otherwise.
`pred.stabilization`	if not `NULL`, will truncate the estimated tau square to `pred.stabilization` * estimated sigma square. This provides additional stability in `BRISC_prediction`. Default value is `1e-8`.
`verbose`	if `TRUE`, model specifications along with information regarding OpenMP support and progress of the algorithm is printed to the screen. Otherwise, nothing is printed to the screen. Default value is `TRUE`.
`eps`	tolerance to be used in centred finite difference approximation of derivatives. Default value is 2e-05.
`nugget_status`	if `nugget_status = 0`, `tau.sq` is fixed to 0, if `nugget_status = 1` `tau.sq` is estimated. Default value is 1.
`ordering`	if `NULL`, the observed locations will be ordered following the scheme in `"order"`. If not `NULL`, the ordering step is bypassed and the input must denote the `n` length integer vector of ordering of the input coordinates that is to be used as the ordering of the coordinates for determination of the set of nearest neighbors. Output from `BRISC_order` can be used here.
`neighbor`	if `NULL`, neighbor set and corresponding information are created using the search type specified in `"search.type"`. If not `NULL`, the step of searching the neighbors is bypassed and the input must be an output from `BRISC_neighbor` with identical input in `"order"`, `"ordering"` and `"search.type"`.
`tol`	the input observation coordinates, response and the covariates are rounded to this many places after the decimal. The default value is 12.

Value

An object of class BRISC_Out, which is a list comprising:

`ord`	the vector of indices used to order data necessary for fitting the NNGP model.
`coords`	the matrix `coords[ord,]`.
`y`	If `stabilization = FALSE`, returns the vector `y[ord]`. If `stabilization = TRUE`, returns `y[ord]` + used white noise in stabilization process.
`X`	the matrix `x[ord,,drop=FALSE]`.
`n.neighbors`	the used value of `n.neighbors`.
`cov.model`	the used covariance model.
`eps`	value of used `eps` for approximate derivation. Default value is 2e-05.
`init`	initial values of the parameters of the covariance model; accounts for `stabilization`.
`Beta`	estimate of beta.
`Theta`	estimate of parameters of covarinace model.
`log_likelihood`	value of Vecchia’s approximate log likelihood with parameter estimates.
`estimation.time`	time (in seconds) required to perform the model fitting after ordering and preprocessing data in `R`, reported using `proc.time()`.
`BRISC_Object`	object required for bootstrap and prediction.

Author(s)

Arkajyoti Saha arkajyotisaha93@gmail.com,
Abhirup Datta abhidatta@jhu.edu

References

Saha, A., & Datta, A. (2018). BRISC: bootstrap for rapid inference on spatial covariances. Stat, e184, DOI: 10.1002/sta4.184.

Datta, A., S. Banerjee, A.O. Finley, and A.E. Gelfand. (2016) Hierarchical Nearest-Neighbor Gaussian process models for large geostatistical datasets. Journal of the American Statistical Association, 111:800-812.

Guinness, J. (2018) Permutation and Grouping Methods for Sharpening Gaussian Process Approximations, Technometrics, DOI: 10.1080/00401706.2018.1437476,
https://github.com/joeguinness/gp_reorder/tree/master/R .

Vecchia, A. V. (1988) Estimation and model identification for continuous spatial processes. Journal of the Royal Statistical Society. Series B (Methodological), 297-312.

Okazaki N. libLBFGS: a library of Limited-memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS),
http://www.chokkan.org/software/liblbfgs/ .

Andrew Finley, Abhirup Datta and Sudipto Banerjee (2020). spNNGP: Spatial Regression Models for Large Datasets using Nearest Neighbor Gaussian Processes. R package version 0.1.4. https://CRAN.R-project.org/package=spNNGP

Ra, S. W., & Kim, J. K. (1993). A fast mean-distance-ordered partial codebook search algorithm for image vector quantization. IEEE Transactions on Circuits and Systems II: Analog and Digital Signal Processing, 40(9), 576-579.

Examples


rmvn <- function(n, mu = 0, V = matrix(1)){
  p <- length(mu)
  if(any(is.na(match(dim(V),p))))
    stop("Dimension not right!")
  D <- chol(V)
  t(matrix(rnorm(n*p), ncol=p)%*%D + rep(mu,rep(n,p)))
}

set.seed(1)
n <- 1000
coords <- cbind(runif(n,0,1), runif(n,0,1))

beta <- c(1,5)
x <- cbind(rnorm(n), rnorm(n))

sigma.sq = 1
phi = 1
tau.sq = 0.1

B <- as.matrix(beta)
D <- as.matrix(dist(coords))
R <- exp(-phi*D)
w <- rmvn(1, rep(0,n), sigma.sq*R)

y <- rnorm(n, x%*%B + w, sqrt(tau.sq))

estimation_result <- BRISC_estimation(coords, y, x)
estimation_result$Theta ##Estimates of covariance model parameters.
estimation_result$Beta ##Estimates of Beta

[Package BRISC version 1.0.5 Index]