Type:	Package
Title:	Covariate-Based Covariance Functions for Nonstationary Spatial Modeling
Version:	0.1.4
Author:	Federico Blasi [aut, cre], Reinhard Furrer [ctb]
Maintainer:	Federico Blasi <federico.blasi@gmail.com>
Description:	Estimation, prediction, and simulation of nonstationary Gaussian process with modular covariate-based covariance functions. Sources of nonstationarity, such as spatial mean, variance, geometric anisotropy, smoothness, and nugget, can be considered based on spatial characteristics. An induced compact-supported nonstationary covariance function is provided, enabling fast and memory-efficient computations when handling densely sampled domains.
Encoding:	UTF-8
LazyData:	true
License:	GPL (≥ 3)
Depends:	R (≥ 3.5.0)
Imports:	Rcpp (≥ 1.0.10), spam (≥ 2.9.1), fields, optimParallel, methods, knitr
LinkingTo:	Rcpp, BH
BugReports:	https://github.com/blasif/cocons/issues
RoxygenNote:	7.3.2
VignetteBuilder:	knitr
NeedsCompilation:	yes
Packaged:	2024-12-12 12:49:31 UTC; blasi
Repository:	CRAN
Date/Publication:	2024-12-12 14:30:02 UTC

Covariate-based Covariance Functions for Nonstationary Gaussian Processes

Description

Provides routines and methods for estimating and predicting nonstationary Gaussian process models with modular covariate-based covariance functions. Several sources of nonstationarity can be modeled based on spatial information, including a spatial mean, marginal standard deviation, local geometric anisotropy, local nugget, and spatially varying smoothness. Each of these components is modeled separately. An induced compact-supported nonstationary covariance function is provided to speed up computations when handling densly sampled domains. Model parameters are estimated via maximum likelihood (and flavours of it, such as penalized and profile maximum likelihood). A variety of functions are also included to compute prediction metrics and to visualize, simulate, and summarize these types of models. Details of the models can be found in the vignette and in coco.

Disclaimer

This package is provided "as is" without warranty of any kind, either express or implied. Backwards compatibility will not be offered until later versions.

Author(s)

Federico Blasi [aut, cre], federico.blasi@gmail.com

Examples

## Not run: 
    vignette("cocons", package = "cocons")
    methods(class = "coco")
  
## End(Not run)

GetNeg2loglikelihood

Description

compute the negative 2 log likelihood based on theta

Usage

GetNeg2loglikelihood(theta, par.pos, locs, x_covariates, 
smooth.limits, z, n, lambda, safe = TRUE)

Arguments

Value

value

Author(s)

Federico Blasi

GetNeg2loglikelihoodProfile

Description

compute the negative 2 log likelihood based on theta

Usage

GetNeg2loglikelihoodProfile(theta, par.pos, locs, x_covariates, 
smooth.limits, z, n, x_betas,lambda, safe = TRUE)

Arguments

Value

value

Author(s)

Federico Blasi

GetNeg2loglikelihoodREML

Description

compute the negative 2 log REML likelihood based on theta

Usage

GetNeg2loglikelihoodREML(theta, par.pos, locs, x_covariates, x_betas,
smooth.limits, z, n, lambda, safe = TRUE)

Arguments

Value

value

Author(s)

Federico Blasi

GetNeg2loglikelihoodTaper

Description

compute the negative 2 log likelihood based on theta

Usage

GetNeg2loglikelihoodTaper(theta, par.pos, ref_taper, locs, 
x_covariates, smooth.limits, cholS, z, n, lambda, safe = TRUE)

Arguments

Value

value

Author(s)

Federico Blasi

GetNeg2loglikelihoodTaperProfile

Description

compute the negative 2 log likelihood based on theta

Usage

GetNeg2loglikelihoodTaperProfile(theta, par.pos, ref_taper, 
locs, x_covariates, smooth.limits, cholS, z, n, lambda, safe = TRUE)

Arguments

Value

(numeric)

Author(s)

Federico Blasi

Creates a coco S4 object

Description

Creates an S4 object of class coco, which is the centerpiece of the cocons package. The function provides a set of consistency checks for ensuring the suitability of the different objects involved.

Usage

coco(type, data, locs, z, model.list, info, output = list())

Arguments

Details

Two types of coco objects are available, each assuming a different type of covariance matrix for the Gaussian process. Type "dense" builds dense covariance matrices (non zero elements), while type "sparse" builds sparse covariance matrices by tapering the dense covariance matrix with a compact isotropic compact-supported correlation matrix [1]. Type "sparse" allows a set of efficient algorithms, thus making it more suitable for large sample sizes.

An important component of the coco S4 class is the model.list specification, involving individual formulas provided as a list, where each of them specifies a covariate-based parametric model for a specific source of nonstationarity. It involves "mean" for the spatial mean, the "std.dev" for the marginal standard deviation, "scale", "aniso" and "tilt", each of them shaping specific aspects of the local spatial geometrically anisotropy structure, "smooth" handling local smoothness, and "nugget" handling the local nugget effect. The models are defined as:

where \boldsymbol{\beta}, \boldsymbol{\alpha}, \boldsymbol{\theta}_1, \boldsymbol{\theta}_2, \boldsymbol{\theta}_3, \boldsymbol{\phi}, and \boldsymbol{\zeta} are the parameter vectors of each model, \nu_{l}, and \nu_{u} are the lower and upper bounds limiting the range of variation of the spatially-varying smoothness, and where \boldsymbol{X}_{\ell} relates to the design matrix defined by the specific models for each of the source of nonstationarity.

Lastly, arguments for the "info" list argument involve:

• "lambda": (numeric) a positive scalar specifying the regularization parameter. Larger values penalizes highly-smoothed long-tailed covariance functions.

• "smooth.limits": (numeric vector) specifying the range of variation for the spatially varying smoothness (e.g. c(0.5, 2.5)).

• "taper": (numeric) specifying the desired taper function from the spam package (only for "sparse" coco objects).

• "delta": (numeric) specifying the taper range/scale (only for "sparse" coco objects).

• "skip.scale": (integer vector) By default, all covariates are scaled. skip.scale allows to specify the index of those variables in data that should not be scaled during the optimization.

Value

(S4) An S4 object of class coco.

Author(s)

Federico Blasi

References

[1] Furrer, Reinhard, Marc G. Genton, and Douglas Nychka. "Covariance tapering for interpolation of large spatial datasets." Journal of Computational and Graphical Statistics 15.3 (2006): 502-523.

See Also

spam::cov.wend1()

Examples

## Not run: 
locs <- expand.grid(seq(0,1,length.out = 10),
seq(0,1,length.out = 10))

toydata <- data.frame('x' = locs[,1])

set.seed(1)
z <- rnorm(100)

model.list <- list('mean' = 0,
                   'std.dev' = formula( ~ 1),
                   'scale' = formula( ~ 1 + x),
                   'aniso' = 0,
                   'tilt' = 0,
                   'smooth' = 3/2,
                   'nugget' = -Inf)
                   
coco_object <- coco(type = 'dense',
                    data = toydata,
                    locs = as.matrix(locs),
                    z = z,
                    model.list = model.list)

coco_object

## End(Not run)

An S4 class to store information

Description

An S4 class to store information

Slots

type

(character) One of two available types "dense" or "sparse". See description.

data

(data.frame) A data.frame with covariates information, where colnames(data) matches model.list specification

locs

(numeric matrix) a matrix with locs matching data

z

(numeric matrix) A matrix of dimension n x p with response values

model.list

(list) A list specifying a model for each aspect of the spatial structure.

info

(list) a list with information about the coco object

output

(list) if building an already fitted coco object (not the standard approach), then requires an output from Optimparallel output, including as well boundaries, etc.

Author(s)

Federico Blasi

Optimizer for coco objects

Description

Estimation the spatial model parameters using the L-BFGS-B optimizer [1].

Usage

cocoOptim(coco.object, boundaries = list(), ncores = "auto", safe = TRUE,
optim.type, optim.control)

Arguments

Value

(S4) An optimized S4 object of class coco.

Author(s)

Federico Blasi

References

[1] Byrd, Richard H., et al. "A limited memory algorithm for bound constrained optimization." SIAM Journal on scientific computing 16.5 (1995): 1190-1208.

[2] Gerber, Florian, and Reinhard Furrer. "optimParallel: An R package providing a parallel version of the L-BFGS-B optimization method." R Journal 11.1 (2019): 352-358.

See Also

[optimParallel]

Examples

## Not run: 
model.list <- list('mean' = 0,
                   'std.dev' = formula( ~ 1 + cov_x + cov_y),
                   'scale' = formula( ~ 1 + cov_x + cov_y),
                   'aniso' = 0,
                   'tilt' = 0,
                   'smooth' = 3/2,
                   'nugget' = -Inf)
                   
coco_object <- coco(type = 'dense',
                    data = holes[[1]][1:100,],
                    locs = as.matrix(holes[[1]][1:100,1:2]),
                    z = holes[[1]][1:100,]$z,
                    model.list = model.list)
                    
optim_coco <- cocoOptim(coco_object,
boundaries = getBoundaries(coco_object,
lower.value = -3, 3))

plotOptimInfo(optim_coco)

plot(optim_coco)

plot(optim_coco, type = 'ellipse')

plot(optim_coco, type = 'correlations', index = c(2,3,5))

summary(optim_coco)
 
getEstims(optim_coco)


## End(Not run)

Prediction for coco objects

Description

Computes the conditional expectation and standard errors based on the conditional Gaussian distribution for nonstationary spatial models.

Usage

cocoPredict(coco.object, newdataset, newlocs, type = 'mean', ...)

Arguments

Value

A list containing:

• systematic: The systematic component of the conditional expectation.

• stochastic: The stochastic component of the conditional expectation.

• sd.pred: The standard errors, when type = 'pred' is specified.

Author(s)

Federico Blasi

Examples

## Not run: 

# Stationary model

model.list_stat <- list('mean' = 0,
'std.dev' = formula( ~ 1),
'scale' = formula( ~ 1),
'aniso' = 0,
'tilt' = 0,
'smooth' = 3/2,
'nugget' = -Inf)

 
model.list_ns <- list('mean' = 0,
'std.dev' = formula( ~ 1 + cov_x + cov_y),
'scale' = formula( ~ 1 + cov_x + cov_y),
'aniso' = 0,
'tilt' = 0,
'smooth' = 3/2,
'nugget' = -Inf)

coco_object <- coco(type = 'dense',
data = holes[[1]][1:100, ],
locs = as.matrix(holes[[1]][1:100, 1:2]),
z = holes[[1]][1:100, ]$z,
model.list = model.list_stat)

optim_coco_stat <- cocoOptim(coco_object,
boundaries = getBoundaries(coco_object,
lower.value = -3, 3))

coco_preds_stat <- cocoPredict(optim_coco_stat, newdataset = holes[[2]],
newlocs = as.matrix(holes[[2]][, 1:2]),
type = "pred")

# Update model
coco_object@model.list <- model.list_ns

optim_coco_ns <- cocoOptim(coco_object,
boundaries = getBoundaries(coco_object,
lower.value = -3, 3))

coco_preds_ns <- cocoPredict(optim_coco_ns, newdataset = holes[[2]],
newlocs = as.matrix(holes[[2]][, 1:2]),
type = "pred")

par(mfrow = c(1, 3))

fields::quilt.plot(main = "full data", holes[[1]][, 1:2], 
holes[[1]]$z, xlim = c(-1, 1), ylim = c(-1, 1))

fields::quilt.plot(main = "stationary se", holes[[2]][, 1:2], 
coco_preds_stat$sd.pred, xlim = c(-1, 1), ylim = c(-1, 1))
fields::quilt.plot(main = "nonstationary se", holes[[2]][, 1:2], 
coco_preds_ns$sd.pred, xlim = c(-1, 1), ylim = c(-1, 1))



## End(Not run)

Marginal and conditional simulation of nonstationary Gaussian processes

Description

draw realizations of stationary and nonstationary Gaussian processes with covariate-based covariance functions.

Usage

cocoSim(coco.object, pars, n, seed, standardize, 
type = 'classic', sim.type = NULL, cond.info = NULL)

Arguments

Details

The argument sim.type = 'cond' specifies a conditional simulation, requiring cond.info to be provided. cond.info is a list including newdataset, a data.frame containing covariates present in model.list at the simulation locations, and newlocs, a matrix specifying the locations corresponding to the simulation, with indexing that matches newdataset.

The argument type = 'classic' assumes a simplified parameterization for the covariance function, with log-parameterizations applied to the parameters std.dev, scale, and smooth.

Value

(matrix) a matrix dim(data)[1] x n.

Author(s)

Federico Blasi

See Also

coco

Examples

## Not run: 

model.list <- list('mean' = 0,
                   'std.dev' = formula( ~ 1 + cov_x + cov_y),
                   'scale' = formula( ~ 1 + cov_x + cov_y),
                   'aniso' = 0,
                   'tilt' = 0,
                   'smooth' = 0.5,
                   'nugget' = -Inf)
                   
coco_object <- coco(type = 'dense',
                    data = holes[[1]][1:1000,],
                    locs = as.matrix(holes[[1]][1:1000,1:2]),
                    z = holes[[1]][1:1000,]$z,
                    model.list = model.list)
                    
coco_sim <- cocoSim(coco.object = coco_object,
            pars = c(0,0.25,0.25,  # pars related to std.dev
            log(0.25),1,-1),       # pars related to scale
            n = 1, 
            standardize = TRUE) 

fields::quilt.plot(coco_object@locs,coco_sim)             

## End(Not run)

Dense covariance function (difference parameterization)

Description

Dense covariance function (difference parameterization)

Usage

cov_rns(theta, locs, x_covariates, smooth_limits)

Arguments

Value

dense covariance matrix

Dense covariance function (classic parameterization)

Description

Dense covariance function (classic parameterization)

Usage

cov_rns_classic(theta, locs, x_covariates)

Arguments

Value

dense covariance matrix with classic parameterization

Dense covariance function

Description

Dense covariance function

Usage

cov_rns_pred(
  theta,
  locs,
  locs_pred,
  x_covariates,
  x_covariates_pred,
  smooth_limits
)

Arguments

Value

dense covariance matrix

Sparse covariance function

Description

Sparse covariance function

Usage

cov_rns_taper(
  theta,
  locs,
  x_covariates,
  colindices,
  rowpointers,
  smooth_limits
)

Arguments

Value

sparse covariance matrix between locs and pred_locs

Sparse covariance function

Description

Sparse covariance function

Usage

cov_rns_taper_pred(
  theta,
  locs,
  locs_pred,
  x_covariates,
  x_covariates_pred,
  colindices,
  rowpointers,
  smooth_limits
)

Arguments

Value

sparse covariance matrix at locs

Retrieve AIC

Description

Retrieve the Akaike information criterion from a fitted coco object.

Usage

getAIC(coco.object)

Arguments

Value

(numeric) the associated AIC value

Author(s)

Federico Blasi

Retrieve BIC

Description

Retrieve BIC from a fitted coco object.

Usage

getBIC(coco.object)

Arguments

Value

(numeric) the associated BIC value

Author(s)

Federico Blasi

Simple build of boundaries

Description

provides a generic set of upper and lower bounds for the L-BFGS-B routine

Usage

getBoundaries(x, lower.value, upper.value)

Arguments

Value

(list) a list with boundaries and simple init values for the optim L-BFGS-B routine

Author(s)

Federico Blasi

Simple build of boundaries (v2)

Description

provides a generic set of upper and lower bounds for the L-BFGS-B routine

Usage

getBoundariesV2(coco.object, mean.limits, std.dev.limits, 
scale.limits, aniso.limits, tilt.limits, smooth.limits, nugget.limits)

Arguments

Value

(list) a list with boundaries for the optim L-BFGS-B routine

Author(s)

Federico Blasi

Simple build of boundaries (v3)

Description

provides a generic set of upper and lower bounds for the L-BFGS-B routine

Usage

getBoundariesV3(coco.object, mean.limits, global.lower, 
std.dev.max.effects, 
scale.max.effects, aniso.max.effects, tilt.max.effects, 
smooth.max.effects, nugget.max.effects)

Arguments

Value

(list) a list with boundaries for the optim L-BFGS-B routine

Author(s)

Federico Blasi

Compute approximate confidence intervals for a coco object

Description

Compute approximate confidence intervals for a (fitted) coco object.

Usage

getCIs(coco.object, inv.hess, alpha = 0.95)

Arguments

Value

(numeric matrix) a matrix with approximate confidence intervals for each parameter in the model.

Author(s)

Federico Blasi

Based on a set of predictions computes the Continuous Ranked Probability Score

Description

Retrieves the Continuous Ranked Probability Score (CRPS) [1].

Usage

getCRPS(z.pred, mean.pred, sd.pred)

Arguments

Value

(numeric vector) retrieves CRPS.

Author(s)

Federico Blasi

References

[1] Gneiting, Tilmann, and Adrian E. Raftery. "Strictly proper scoring rules, prediction, and estimation." Journal of the American statistical Association 102.477 (2007): 359-378.

Covariance matrix for "coco" class

Description

Compute the covariance matrix of coco.object.

Usage

getCovMatrix(coco.object, type = 'global', index = NULL)

Arguments

Value

(matrix or S4) a n x n covariance matrix (for 'dense' coco objects) or a S4 spam object (for 'sparse' coco objects).

Author(s)

Federico Blasi

Examples

## Not run: 
model.list <- list('mean' = 0,
                   'std.dev' = formula( ~ 1 + cov_x + cov_y),
                   'scale' = formula( ~ 1 + cov_x + cov_y),
                   'aniso' = 0,
                   'tilt' = 0,
                   'smooth' = 3/2,
                   'nugget' = -Inf)
                   
coco_object <- coco(type = 'dense',
                    data = holes[[1]][1:100,],
                    locs = as.matrix(holes[[1]][1:100,1:2]),
                    z = holes[[1]][1:100,]$z,
                    model.list = model.list)
                    
optim_coco <- cocoOptim(coco_object,
boundaries = getBoundaries(coco_object,
lower.value = -3, 3))

getCovMatrix(optim_coco)


## End(Not run)

Based on a specific taper scale (delta), retrieves the density of the covariance matrix.

Description

Based on a specific taper scale (delta), retrieves the density of the covariance matrix.

Usage

getDensityFromDelta(coco.object, delta)

Arguments

Value

(numeric vector) the associate density of the tapered covariance matrix.

Author(s)

Federico Blasi

Create an efficient design matrix based on a list of aspect models

Description

Creates a unique design matrix based on model specification for each of the different potentially spatially varying aspects.

Usage

getDesignMatrix(model.list, data)

Arguments

Value

(list) a list with two elements: a design matrix of dimension (n x p), and a par.pos object, indexing columns of the design matrix to each of the spatially-varying functions.

Author(s)

Federico Blasi

Retrieve estimates from a fitted coco object

Description

Retrieve estimates from a fitted coco object.

Usage

getEstims(coco.object)

Arguments

Value

(list) a list with the estimates parameters for the different aspects

Author(s)

Federico Blasi

getHessian

Description

numerically approximate the Hessian. Hessians of parameters based on "pml" are based on full likelihoods.

Usage

getHessian(coco.object, ncores = parallel::detectCores() - 1, 
eps = .Machine$double.eps^(1/4))

Arguments

Value

(numeric matrix) a symmetric matrix pxp of the approximated (observed) Hessian

Author(s)

Federico Blasi

Based on a set of predictions computes the Log-Score

Description

Computes the Log-Score [1].

Usage

getLogScore(z.pred, mean.pred, sd.pred)

Arguments

Value

(numeric vector) retrieves Log-Score.

Author(s)

Federico Blasi

References

[1] Gneiting, Tilmann, and Adrian E. Raftery. "Strictly proper scoring rules, prediction, and estimation." Journal of the American statistical Association 102.477 (2007): 359-378.

Retrieve the loglikelihood value

Description

Retrieve the loglikelihood value from a fitted coco object.

Usage

getLoglik(coco.object)

Arguments

Value

(numeric) wrap for value from a OptimParallel object

Author(s)

Federico Blasi

Retrieves the modified inverse of the hessian

Description

Based on the inverse of the Hessian (based on the difference parameterization for the std.dev and scale parameters), retrieves the modified inverse of the hessian (i.e. std.dev and scale).

Usage

getModHess(coco.object, inv.hess)

Arguments

Value

(numeric matrix) the modified inverse of the hessian matrix

Author(s)

Federico Blasi

Builds the necessary input for building covariance matrices

Description

Returns a list of parameter vectors for each of the aspects.

Usage

getModelLists(theta, par.pos, type = 'diff')

Arguments

Value

(list) a list of different spatial aspects and mean required for the cov.rns functions

Author(s)

Federico Blasi

Fast and simple standardization for the design matrix.

Description

Centers and scale the design matrix.

Usage

getScale(x, mean.vector = NULL, sd.vector = NULL)

Arguments

Value

(list) a list with a scaled design matrix of dimension n x (p+1), and a set of mean and sd vectors employed to scale the matrix

Author(s)

Federico Blasi

Evaluates the spatially-varying functions from a coco object at locs

Description

Evaluates the spatially-varying functions of the nonstationary spatial structure.

Usage

getSpatEffects(coco.object)

Arguments

Value

(list) a list with the different estimated surfaces.

Author(s)

Federico Blasi

Computes the spatial mean of a (fitted) coco object

Description

Computes the spatial mean of the (fitted) coco object.

Usage

getSpatMean(coco.object)

Arguments

Value

(numeric vector) a vector with the adjusted trend.

Author(s)

Federico Blasi

Holes Data Set

Description

The synthetic "holes" provides a set of training and test data.frame of a Gaussian process realization with a (inherently dense) nonstationary covariance function. Four holes are present in the training dataset, and the task is to predict them.

Usage

holes

Format

A list with training and test data.frame with rows and variables:

x

first spatial coordinate

y

second spatial coordinate

cox_x

first spatial characteristic

cov_y

second spatial characteristic

z

response variable

Source

Source of the data

Examples

data(holes)

Holes with trend + multiple realizations Data Set

Description

The synthetic "holes_bm" provides a set of training and test data.frame of a Gaussian process realization with a (inherently dense) nonstationary covariance function. Four holes are present in the training dataset, and the task is to predict them. This version provides ten independent realizations of the process, as well as considers a spatial mean effect.

Usage

holes_bm

Format

A list with training, training.z, test, and test.z data.frames with rows and variables:

x

first spatial coordinate

y

second spatial coordinate

cox_x

first spatial characteristic

cov_y

second spatial characteristic

cov_z

third spatial characteristic

z.i

i-th response variable

Source

Source of the data

Examples

data(holes_bm)

check whether an R object is a formula

Description

check whether an R object is a formula

Usage

is.formula(x)

Arguments

Value

TRUE/FALSE

Author(s)

Federico Blasi

Plot Method for coco objects

Description

This method plots objects of class coco.

Usage

## S4 method for signature 'coco,missing'
plot(x, y, type = NULL, index = NULL, factr = 0.1, ...)

Arguments

Value

Several plots are created.

Author(s)

Federico Blasi

Plot log info detailed

Description

plot output of optim

Usage

plotOptimInfo(coco.object, ...)

Arguments

Value

Outputs a sequence of plots detailing parameters during the optimization routine

Author(s)

Federico Blasi

See Also

cocoOptim()

Stripes Data Set

Description

The synthetic "stripes" provides a set of training and test data.frame of a Gaussian process realization with a (inherently sparse) nonstationary covariance function. Several stripes are present in the training dataset, and the task is to predict them.

Usage

stripes

Format

A list with training and test data.frame with rows and variables:

x

first spatial coordinate

y

second spatial coordinate

cox_x

first spatial characteristic

cov_y

second spatial characteristic

cov_xy

third spatial characteristic

z

response variable

Source

Source of the data

Examples

data(stripes)

Summary Method for Coco Class

Description

method summary for objects of class 'coco'.

Usage

## S4 method for signature 'coco'
summary(object, inv.hess = NULL)

Arguments

Value

summary the coco object

Author(s)

Federico Blasi