stcos: Space-Time Change of Support

Description

An R Package for Space-Time Change of Support (STCOS) modeling.

Details

Supports building and running STCOS and related models. A guide on package use is given in <arXiv:1904.12092>.

References

Jonathan R. Bradley, Christopher K. Wikle, and Scott H. Holan (2015). Spatio-temporal change of support with application to American Community Survey multi-year period estimates. STAT 4 pp.255-270. doi:10.1002/sta4.94.

Andrew M. Raim, Scott H. Holan, Jonathan R. Bradley, and Christopher K. Wikle (2017). A model selection study for spatio-temporal change of support. In JSM Proceedings, Government Statistics Section. Alexandria, VA: American Statistical Association, pp.1524-1540.

Andrew M. Raim, Scott H. Holan, Jonathan R. Bradley, and Christopher K. Wikle (2020+). Spatio-Temporal Change of Support Modeling with R. https://arxiv.org/abs/1904.12092.

Best Approximation to Covariance Structure

Description

Compute the best positive approximant for use in the STCOS model, under several prespecified covariance structures.

Usage

cov_approx_randwalk(Delta, S)

cov_approx_blockdiag(Delta, S)

Arguments

Details

Let \bm{\Sigma} be an N \times N symmetric and positive-definite covariance matrix and \bm{S} be an N \times r matrix with rank r. The objective is to compute a matrix \bm{K} which minimizes the Frobenius norm

\Vert \bm{\Sigma} - \bm{S} \bm{C} \bm{S}^\top {\Vert}_\textrm{F},

over symmetric positive-definite matrices \bm{C}. The solution is given by

\bm{K} = (\bm{S}^\top \bm{S})^{-1} \bm{S}^\top \bm{\Sigma} \bm{S} (\bm{S}^\top \bm{S})^{-1}.

In the STCOS model, \bm{S} represents the design matrix from a basis function computed from a fine-level support having n areas, using T time steps. Therefore N = n T represents the dimension of covariance for the fine-level support.

We provide functions to handle some possible structures for target covariance matrices of the form

\bm{\Sigma} = \left( \begin{array}{ccc} \bm{\Gamma}(1,1) & \cdots & \bm{\Gamma}(1,T) \\ \vdots & \ddots & \vdots \\ \bm{\Gamma}(T,1) & \cdots & \bm{\Gamma}(T,T) \end{array} \right),

where each \bm{\Gamma}(s,t) is an n \times n matrix.

• cov_approx_randwalk assumes \bm{\Sigma} is based on the autocovariance function of a random walk

  \bm{Y}_{t+1} = \bm{Y}_{t} + \bm{\epsilon}_t, \quad \bm{\epsilon}_t \sim \textrm{N}(\bm{0}, \bm{\Delta}).

  so that

  \bm{\Gamma}(s,t) = \min(s,t) \bm{\Delta}.

• cov_approx_blockdiag assumes \bm{\Sigma} is based on

  \bm{Y}_{t+1} = \bm{Y}_{t} + \bm{\epsilon}_t, \quad \bm{\epsilon}_t \sim \textrm{N}(\bm{0}, \bm{\Delta}).

  which are independent across t, so that

  \bm{\Gamma}(s,t) = I(s = t) \bm{\Delta},

The block structure is used to reduce the computational burden, as N may be large.

Deviance Information Criterion

Description

Generic function to calculate Deviance Information Criterion (DIC) for a given model object.

Usage

DIC(object, ...)

Arguments

Value

A numeric value of the DIC.

Shapes and ACS estimates for Boone County, MO.

Description

An sf object with ACS estimates for:

• Boone County, Missouri

• Table B19013

• Block group level geography

• Years 2013 - 2017

Usage

acs5_2013

acs5_2014

acs5_2015

acs5_2016

acs5_2017

Format

sf objects.

An object of class sf (inherits from data.frame) with 87 rows and 9 columns.

An object of class sf (inherits from data.frame) with 87 rows and 9 columns.

An object of class sf (inherits from data.frame) with 85 rows and 9 columns.

An object of class sf (inherits from data.frame) with 87 rows and 9 columns.

An object of class sf (inherits from data.frame) with 87 rows and 9 columns.

Details

This dataset was assembled using shapefiles from the tigris package and ACS estimates from the American FactFinder website on 2/28/2019. See data-prep-aff.R in the Columbia example code for details. American FactFinder has since been deprecated, and similar data are available at http://data.census.gov.

Sparse adjacency matrix between two sets of areas.

Description

A convenience function to convert output from sf::st_touches to a sparse matrix as defined in the Matrix package.

Usage

adjacency_matrix(dom)

Arguments

Details

Returns a matrix A whose (i,j)th entry contains a 1 if areal units dom[i,] and dom[j,] are adjacent; 0 otherwise.

Value

An adjacency matrix

Examples

data("acs_sf")
dom = acs5_2013[1:4,]
A = adjacency_matrix(dom)

Areal Space-Time Bisquare Basis

Description

Space-Time bisquare basis on areal data.

Usage

areal_spacetime_bisquare(dom, period, knots, w_s, w_t, control = NULL)

Arguments

Details

Notes about arguments:

• knots may be provided as either an sf or sfc object, or as a matrix of points.

• If an sf or sfc object is provided for knots, r three-dimensional POINT entries are expected in st_geometry(knots). Otherwise, knots will be interpreted as an r \times 3 numeric matrix.

• If knots is an sf or sfc object, it is checked to ensure the coordinate system matches dom.

For each area A in the given domain, and time period \bm{v} = (v_1, \ldots, v_m) compute the basis functions

\psi_j^{(m)}(A, \bm{v}) = \frac{1}{m} \sum_{k=1}^m \frac{1}{|A|} \int_A \psi_j(\bm{u},v_k) d\bm{u},

for j = 1, \ldots, r. Here, \varphi_j{(\bm{u},v)} represent spacetime_bisquare basis functions defined at the point level using \bm{c}_j, g_j, w_s, and w_t.

The basis requires an integration which may be computed using one of two methods. The mc method uses a Monte Carlo approximation

\psi_j^{(m)}(A, \bm{v}) \approx \frac{1}{m} \sum_{k=1}^m \frac{1}{Q} \sum_{q=1}^Q \psi_j(\bm{u}_q, v_k),

based on a random sample of locations \bm{u}_1, \ldots, \bm{u}_Q from a uniform distribution on area A. The rect method uses a simple quadrature approximation

\psi_j^{(m)}(A, \bm{v}) \approx \frac{1}{m} \sum_{k=1}^m \frac{1}{|A|} \sum_{a=1}^{n_x} \sum_{b=1}^{n_y} \psi_j(\bm{u}_{ab}, v_k) I(\bm{u}_{ab} \in A) \Delta_x \Delta_y.

Here, the bounding box st_bbox(A) is divided evenly into a grid of n_x \times n_y rectangles, each of size \Delta_x \times \Delta_y. Each \bm{u}_{ab} = (u_a, u_b) is a point from the (a,b)th rectangle, for a = 1, \ldots, n_x and b = 1, \ldots, n_y.

Due to the treatment of A_i and \bm{c}_j as objects in a Euclidean space, this basis is more suitable for coordinates from a map projection than coordinates based on a globe representation.

The control argument is a list which may provide any of the following:

• method specifies computation method: mc or rect. Default is mc.

• mc_reps is number of repetitions to use for mc. Default is 1000.

• nx is number of x-axis points to use for rect method. Default is 50.

• ny is number of y-axis points to use for rect method. Default is 50.

• report_period is an integer; print a message with progress each time this many areas are processed. Default is Inf so that message is suppressed.

• verbose is a logical; if TRUE print descriptive messages about the computation. Default is FALSE.

• mc_sampling_factor is a positive number; an oversampling factor used to compute blocksize in the rdomain function. I.e., blocksize = ceiling(mc_sampling_factor * mc_reps). Default is 1.2.

Value

A sparse n \times r matrix whose ith row is \bm{s}_i^\top = \Big( \psi_1^{(m)}(A_i), \ldots, \psi_r^{(m)}(A_i) \Big).

See Also

Other bisquare: areal_spatial_bisquare(), spacetime_bisquare(), spatial_bisquare()

Examples

set.seed(1234)

# Create knot points
seq_x = seq(0, 1, length.out = 3)
seq_y = seq(0, 1, length.out = 3)
seq_t = seq(0, 1, length.out = 3)
knots = expand.grid(x = seq_x, y = seq_y, t = seq_t)
knots_sf = st_as_sf(knots, coords = c("x","y","t"), crs = NA, dim = "XYM", agr = "constant")

# Create a simple domain (of rectangles) to evaluate
shape1 = matrix(c(0.0,0.0, 0.5,0.0, 0.5,0.5, 0.0,0.5, 0.0,0.0), ncol=2, byrow=TRUE)
shape2 = shape1 + cbind(rep(0.5,5), rep(0.0,5))
shape3 = shape1 + cbind(rep(0.0,5), rep(0.5,5))
shape4 = shape1 + cbind(rep(0.5,5), rep(0.5,5))
sfc = st_sfc(
   st_polygon(list(shape1)),
   st_polygon(list(shape2)),
   st_polygon(list(shape3)),
   st_polygon(list(shape4))
)
dom = st_sf(data.frame(geoid = 1:length(sfc), geom = sfc))

rad = 0.5
period = c(0.4, 0.7)
areal_spacetime_bisquare(dom, period, knots, w = rad, w_t = 1)
areal_spacetime_bisquare(dom, period, knots_sf, w_s = rad, w_t = 1)

# Plot the (spatial) knots and the (spatial) domain at which we evaluated
# the basis.
plot(knots[,1], knots[,2], pch = 4, cex = 1.5, col = "red")
plot(dom[,1], col = NA, add = TRUE)

# Draw a circle representing the basis' radius around one of the knot points
tseq = seq(0, 2*pi, length=100) 
coords = cbind(rad * cos(tseq) + seq_x[2], rad * sin(tseq) + seq_y[2])
lines(coords, col = "red")

Areal Spatial Bisquare Basis

Description

Spatial bisquare basis on areal data.

Usage

areal_spatial_bisquare(dom, knots, w, control = NULL)

Arguments

Details

Notes about arguments:

• knots may be provided as either an sf or sfc object, or as a matrix of points.

• If an sf or sfc object is provided for knots, r two-dimensional POINT entries are expected in st_geometry(knots). Otherwise, knots will be interpreted as an r \times 2 numeric matrix.

• If knots is an sf or sfc object, it is checked to ensure the coordinate system matches dom.

For each area A in the given domain, compute an the basis functions

\bar{\varphi}_j(A) = \frac{1}{|A|} \int_A \varphi_j(\bm{u}) d\bm{u}

for j = 1, \ldots, r. Here, \varphi_j(\bm{u}) represent spatial_bisquare basis functions defined at the point level using \bm{c}_j and w.

The basis requires an integration which may be computed using one of two methods. The mc method uses

\bar{\varphi}_j(A) \approx \frac{1}{Q} \sum_{q=1}^Q \varphi_j(\bm{u}_q),

based on a random sample of locations \bm{u}_1, \ldots, \bm{u}_Q from a uniform distribution on area A. The rect method uses a simple quadrature approximation

\bar{\varphi}_j(A) \approx \frac{1}{|A|} \sum_{a=1}^{n_x} \sum_{b=1}^{n_y} \varphi_j(\bm{u}_{ab}) I(\bm{u}_{ab} \in A) \Delta_x \Delta_y.

Here, the bounding box st_bbox(A) is divided evenly into a grid of n_x \times n_y rectangles, each of size \Delta_x \times \Delta_y. Each \bm{u}_{ab} = (u_a, u_b) is a point from the (a,b)th rectangle, for a = 1, \ldots, n_x and b = 1, \ldots, n_y.

Due to the treatment of A_i and \bm{c}_j as objects in a Euclidean space, this basis is more suitable for coordinates from a map projection than coordinates based on a globe representation.

The control argument is a list which may provide any of the following:

• method specifies computation method: mc or rect. Default is mc.

• mc_reps is number of repetitions to use for mc. Default is 1000.

• nx is number of x-axis points to use for rect method. Default is 50.

• ny is number of y-axis oints to use for rect method. Default is 50.

• report_period is an integer; print a message with progress each time this many areas are processed. Default is Inf so that message is suppressed.

• verbose is a logical; if TRUE print descriptive messages about the computation. Default is FALSE.

• mc_sampling_factor is a positive number; an oversampling factor used to compute blocksize in the rdomain function. I.e., blocksize = ceiling(mc_sampling_factor * mc_reps). Default is 1.2.

Value

A sparse n \times r matrix whose ith row is \bm{s}_i^\top = \Big( \bar{\varphi}_1(A_i), \ldots, \bar{\varphi}_r(A_i) \Big).

See Also

Other bisquare: areal_spacetime_bisquare(), spacetime_bisquare(), spatial_bisquare()

Examples

set.seed(1234)

# Create knot points
seq_x = seq(0, 1, length.out = 3)
seq_y = seq(0, 1, length.out = 3)
knots = expand.grid(x = seq_x, y = seq_y)
knots_sf = st_as_sf(knots, coords = c("x","y"), crs = NA, agr = "constant")

# Create a simple domain (of rectangles) to evaluate
shape1 = matrix(c(0.0,0.0, 0.5,0.0, 0.5,0.5, 0.0,0.5, 0.0,0.0), ncol=2, byrow=TRUE)
shape2 = shape1 + cbind(rep(0.5,5), rep(0.0,5))
shape3 = shape1 + cbind(rep(0.0,5), rep(0.5,5))
shape4 = shape1 + cbind(rep(0.5,5), rep(0.5,5))
sfc = st_sfc(
   st_polygon(list(shape1)),
   st_polygon(list(shape2)),
   st_polygon(list(shape3)),
   st_polygon(list(shape4))
)
dom = st_sf(data.frame(geoid = 1:length(sfc), geom = sfc))

rad = 0.5
areal_spatial_bisquare(dom, knots, rad)
areal_spatial_bisquare(dom, knots_sf, rad)

# Plot the knots and the points at which we evaluated the basis
plot(knots[,1], knots[,2], pch = 4, cex = 1.5, col = "red")
plot(dom[,1], col = NA, add = TRUE)

# Draw a circle representing the basis' radius around one of the knot points
tseq = seq(0, 2*pi, length=100) 
coords = cbind(rad * cos(tseq) + seq_x[2], rad * sin(tseq) + seq_y[2])
lines(coords, col = "red")

Compute the autocovariance matrix for a VAR(1) process.

Description

Compute the autocovariance matrix for a VAR(1) process.

Usage

autocov_VAR1(A, Sigma, lag_max)

Arguments

Details

Computes the autocovariance matrix \bm{\Gamma}(h) of the m-dimensional VAR(1) process

\bm{Y}_t = \bm{A} \bm{Y}_{t-1} + \bm{\epsilon}_t, \quad \bm{\epsilon}_t \sim \textrm{N}(\bm{0}, \bm{\Sigma})

For the required computation of \bm{\Gamma}(0), this function solves the m^2 \times m^2 system

\textrm{vec}(\bm{\Gamma}(0)) = [\bm{I} - \bm{A} \otimes \bm{A}]^{-1} \textrm{vec}(\bm{\Sigma}).

without directly computing m^2 \times m^2 matrices.

Value

An array Gamma of dimension c(m, m, lag_max + 1), where the slice Gamma[,,h] represents the autocovariance at lag h = 0, 1, ..., lag_max.

Examples

U = matrix(NA, 3, 3)
U[,1] = c(1, 1, 1) / sqrt(3)
U[,2] = c(1, 0, -1) / sqrt(2)
U[,3] = c(0, 1, -1) / sqrt(2)
B = U %*% diag(c(0.5, 0.2, 0.1)) %*% t(U)
A = (B + t(B)) / 2
Sigma = diag(x = 2, nrow = 3)
autocov_VAR1(A, Sigma, lag_max = 5)

CAR Precision Matrix

Description

A convenience function to compute the CAR precision matrix based on a given adjacency matrix.

Usage

car_precision(A, tau = 1, scale = FALSE)

Arguments

Details

Suppose \bm{A} is an n \times n adjacency matrix and \bm{D} = \textrm{Diag}(\bm{A} \bm{1}) = \textrm{Diag}(a_{1+}, \ldots, a_{n+}). If scale is FALSE, return the CAR precision matrix

\bm{Q} = \bm{D} - \tau \bm{A}.

If scale is TRUE, return a scaled version

\tilde{\bm{Q}} = \bm{D}^{-1} \bm{Q}.

An error is thrown if scale = TRUE and any of \{ a_{1+}, \ldots, a_{n+} \} are equal to 0. Taking \tau = 1 corresponds to the Intrinsic CAR precision matrix.

Typically in a modeling context, the precision matrix will be multiplied by a scaling parameter; e.g., a CAR model for random effects \bm{\phi} could be

f(\bm{\phi} \mid \alpha) \propto \alpha^{-q} \exp\left\{ -\frac{1}{2 \alpha^2} \bm{\phi}^\top \bm{Q} \bm{\phi} \right\}.

where q = \textrm{rank}(Q).

Value

CAR precision matrix.

Examples

data("acs_sf")
dom = acs5_2013[1:4,]
A = adjacency_matrix(dom)
Q = car_precision(A)

City of Columbia neighborhoods.

Description

An sf object containing the geometry of four neighborhoods in the City of Columbia, Boone County, Missouri. Based on shapefiles provided by the Office of Information Technology / GIS, City of Columbia, Missouri.

Usage

columbia_neighbs

Format

An sf object with 4 features (neighborhoods).

Gibbs Sampler for STCOS Model

Description

Gibbs Sampler for STCOS Model

Usage

gibbs_stcos(
  z,
  v,
  H,
  S,
  Kinv,
  R,
  report_period = R + 1,
  burn = 0,
  thin = 1,
  init = NULL,
  fixed = NULL,
  hyper = NULL
)

## S3 method for class 'stcos_gibbs'
logLik(object, ...)

## S3 method for class 'stcos_gibbs'
DIC(object, ...)

## S3 method for class 'stcos_gibbs'
print(x, ...)

## S3 method for class 'stcos_gibbs'
fitted(object, H, S, ...)

## S3 method for class 'stcos_gibbs'
predict(object, H, S, ...)

Arguments

Details

Fits the model

\bm{Z} = \bm{H} \bm{\mu}_B + \bm{S} \bm{\eta} + \bm{\xi} + \bm{\varepsilon}, \quad \bm{\varepsilon} \sim \textrm{N}(0, \bm{V}),

\bm{\eta} \sim \textrm{N}(\bm{0}, \sigma_K^2 \bm{K}), \quad \bm{\xi} \sim \textrm{N}(0, \sigma_{\xi}^2 \bm{I}),

\bm{\mu}_B \sim \textrm{N}(\bm{0}, \sigma_\mu^2 \bm{I}), \quad \sigma_\mu^2 \sim \textrm{IG}(a_\mu, b_\mu),

\sigma_K^2 \sim \textrm{IG}(a_K, b_K), \quad \sigma_\xi^2 \sim \textrm{IG}(a_\xi, b_\xi),

using a Gibbs sampler with closed-form draws.

Helper functions produce the following outputs:

• logLik computes the log-likelihood for each saved draw.

• DIC computes the Deviance information criterion for each saved draw.

• print displays a summary of the draws.

• fitted computes the mean E(Y_i) for each observation i = 1, \ldots, n, for each saved draw.

• predict draws Y_i for each observation i = 1, \ldots, n, using the parameter values for each saved Gibbs sampler draw.

Value

gibbs_stcos returns an stcos object which contains draws from the sampler. Helper functions take this object as an input and produce various outputs (see details).

Examples

## Not run: 
demo = prepare_stcos_demo()
out = gibbs_stcos(demo$z, demo$v, demo$H, demo$S, solve(demo$K),
    R = 100, burn = 0, thin = 1)
print(out)
logLik(out)
DIC(out)
fitted(out, demo$H, demo$S)
predict(out, demo$H, demo$S)

## End(Not run)

licols

Description

Extract a linearly independent set of columns of a matrix.

Usage

licols(X, tol = 1e-10, quiet = FALSE)

Arguments

Details

An R version of a Matlab licols function given in this MathWorks forum post.

Value

Xsub contains the extracted columns of X and idx contains the indices (into X) of those columns. The elapsed time is stored in elapsed.sec.

Examples

x = 0:19 %% 3 + 1
Z = model.matrix(~ as.factor(x) - 1)
X = cbind(1, Z)
licols(X)

MLE for STCOS Model

Description

MLE for STCOS Model

Usage

mle_stcos(
  z,
  v,
  H,
  S,
  K,
  init = NULL,
  optim_control = list(),
  optim_method = "L-BFGS-B"
)

Arguments

Details

Maximize the likelihood of the STCOS model

f(\bm{z} \mid \bm{\mu}_B, \sigma_K^2, \sigma_\xi^2) = \textrm{N}(\bm{z} \mid \bm{H} \bm{\mu}_B, \bm{\Delta} ), \quad \bm{\Delta} = \sigma_\xi^2 \bm{I} + \bm{V} + \sigma_K^2 \bm{S} \bm{K} \bm{S}^\top,

by numerical maximization of the profile likelihood

\ell(\sigma_K^2, \sigma_\xi^2) = -\frac{N}{2} \log(2 \pi) -\frac{1}{2} \log |\bm{\Delta}| -\frac{1}{2} (\bm{z} - \bm{H} \hat{\bm{\mu}}_B)^\top \bm{\Delta}^{-1} (\bm{z} - \bm{H} \hat{\bm{\mu}}_B)

using \hat{\bm{\mu}}_B = (\bm{H}^\top \bm{\Delta}^{-1} \bm{H})^{-1} \bm{H}^\top \bm{\Delta}^{-1} \bm{z}.

Value

A list containing maximum likelihood estimates.

Examples

## Not run: 
demo = prepare_stcos_demo()
mle_out = mle_stcos(demo$z, demo$v, demo$S, demo$H, demo$K)
sig2K_hat = mle_out$sig2K_hat
sig2xi_hat = mle_out$sig2xi_hat
mu_hat = mle_out$mu_hat

## End(Not run)

Matrix of overlaps between two sets of areas.

Description

A convenience function to convert output from sf::st_intersection to a sparse matrix as defined in the Matrix package.

Usage

overlap_matrix(dom1, dom2, proportion = TRUE)

Arguments

Details

Returns a matrix H whose (i,j)th entry represent the area of the overlap between areal units dom1[i,] and dom2[j,].

Value

An matrix of overlaps.

Examples

data("acs_sf")
dom1 = acs5_2013[1:10,]
dom2 = acs5_2016[1:10,]
H1 = overlap_matrix(dom1, dom2)
H2 = overlap_matrix(dom1, dom2, proportion = FALSE)

Prepare Demo Data for STCOS Model

Description

Create demo data based on ACS example, making a few simple model choices. The purpose of this function is to facilitate examples in other functions. Uses functions in the package to create model terms from shapefiles.

Usage

prepare_stcos_demo(num_knots_sp = 200, basis_mc_reps = 200, eigval_prop = 0.65)

Arguments

Value

A list containing the following:

• z direct estimates.

• v direct variance estimates.

• H overlap matrix.

• S design matrix of basis expansion.

• K covariance matrix of the random effect.

Examples

## Not run: 
out = prepare_stcos_demo()

## End(Not run)

Draw uniformly distributed points from a set of areas

Description

An alternative to sf::st_sample which draws uniformly distributed points using a simple accept-reject method.

Usage

rdomain(n, dom, blocksize = n, itmax = Inf)

Arguments

Details

Draws a sample of blocksize points uniformly from a bounding box on dom, and accepts only the points which belong to dom. This yields a uniform sample on dom. The process is repeated until n accepted draws are obtained, or until it has been attempted itmax times. If itmax iterations are reached without accepting n draws, an error is thrown.

This seems to be an order of magnitude faster than the current implementation of st_sample, although the latter can accomplish the same objective and is more general. The improved performance is worthwhile when used in the areal basis functions, which sample repeatedly from the domain.

Performance will degrade when areal units have small area relative to their bounding box, as many candidate points may need to be discarded. For example, this will occur if dom contains a set of small scattered islands in an ocean. In this case, it would be more efficient to sample from each island at a time.

Value

An sf object with 2-dimensional points.

Examples

dom = acs5_2013[c(1,5,8,12),]
pts = rdomain(10000, dom)

Space-Time Bisquare Basis

Description

Space-time bisquare basis on point data.

Usage

spacetime_bisquare(dom, knots, w_s, w_t)

Arguments

Details

Notes about arguments:

• Both dom and knots may be provided as either sf or sfc objects, or as matrices of points.

• If an sf or sfc object is provided for dom, n three-dimensional POINT entries are expected in st_geometry(dom). Otherwise, dom will be interpreted as an n \times 3 numeric matrix.

• If an sf or sfc object is provided for knots, r three-dimensional POINT entries are expected in st_geometry(knots). Otherwise, knots will be interpreted as an r \times 3 numeric matrix.

• If both dom and knots_s are given as sf or sfc objects, they will be checked to ensure a common coordinate system.

For each (\bm{u}_i,v_i), compute the basis functions

\psi_j(\bm{u},v) = \left[ 2 - \frac{\Vert \bm{u} - \bm{c}_j \Vert^2}{w_s^2}- \frac{|v - g_j|^2}{w_t^2} \right]^2 \cdot I(\Vert \bm{u} - \bm{c}_j \Vert \leq w_s) \cdot I(|v - g_j| \leq w_t)

for j = 1, \ldots, r.

Due to the treatment of \bm{u}_i and \bm{c}_j as points in a Euclidean space, this basis is more suitable for coordinates from a map projection than coordinates based on a globe representation.

Value

A sparse n \times r matrix whose ith row is

\bm{s}_i^\top = \Big( \psi_1(\bm{u}_i,v_i), \ldots, \psi_r(\bm{u}_i,v_i) \Big).

See Also

Other bisquare: areal_spacetime_bisquare(), areal_spatial_bisquare(), spatial_bisquare()

Examples

set.seed(1234)

# Create knot points
seq_x = seq(0, 1, length.out = 3)
seq_y = seq(0, 1, length.out = 3)
seq_t = seq(0, 1, length.out = 3)
knots = expand.grid(x = seq_x, y = seq_y, t = seq_t)
knots_sf = st_as_sf(knots, coords = c("x","y","t"), crs = NA, dim = "XYM", agr = "constant")

# Points to evaluate
x = runif(50)
y = runif(50)
t = sample(1:3, size = 50, replace = TRUE)
pts = data.frame(x = x, y = y, t = t)
dom = st_as_sf(pts, coords = c("x","y","t"), crs = NA, dim = "XYM", agr = "constant")

rad = 0.5
spacetime_bisquare(cbind(x,y,t), knots, w_s = rad, w_t = 1)
spacetime_bisquare(dom, knots_sf, w_s = rad, w_t = 1)

# Plot the (spatial) knots and the points at which we evaluated the basis
plot(knots[,1], knots[,2], pch = 4, cex = 1.5, col = "red")
text(x, y, labels = t, cex = 0.75)

# Draw a circle representing the basis' radius around one of the knot points
tseq = seq(0, 2*pi, length=100) 
coords = cbind(rad * cos(tseq) + seq_x[2], rad * sin(tseq) + seq_y[2])
lines(coords, col = "red")

Spatial Bisquare Basis

Description

Spatial bisquare basis on point data.

Usage

spatial_bisquare(dom, knots, w)

Arguments

Details

Notes about arguments:

• Both dom and knots may be provided as either sf or sfc objects, or as matrices of points.

• If an sf or sfc object is provided for dom, n two-dimensional POINT entries are expected in st_geometry(dom). Otherwise, dom will be interpreted as an n \times 2 numeric matrix.

• If an sf or sfc object is provided for knots, r two-dimensional POINT entries are expected in st_geometry(knots). Otherwise, knots will be interpreted as an r \times 2 numeric matrix.

• If both dom and knots are given as sf or sfc objects, they will be checked to ensure a common coordinate system.

For each \bm{u}_i, compute the basis functions

\varphi_j(\bm{u}) = \left[ 1 - \frac{\Vert\bm{u} - \bm{c}_j \Vert^2}{w^2} \right]^2 \cdot I(\Vert \bm{u} - \bm{c}_j \Vert \leq w)

for j = 1, \ldots, r.

Due to the treatment of \bm{u}_i and \bm{c}_j as points in a Euclidean space, this basis is more suitable for coordinates from a map projection than coordinates based on a globe representation.

Value

A sparse n \times r matrix whose ith row is \bm{s}_i^\top = \Big( \varphi_1(\bm{u}_i), \ldots, \varphi_r(\bm{u}_i) \Big).

See Also

Other bisquare: areal_spacetime_bisquare(), areal_spatial_bisquare(), spacetime_bisquare()

Examples

set.seed(1234)

# Create knot points
seq_x = seq(0, 1, length.out = 3)
seq_y = seq(0, 1, length.out = 3)
knots = expand.grid(x = seq_x, y = seq_y)
knots_sf = st_as_sf(knots, coords = c("x","y"), crs = NA, agr = "constant")

# Points to evaluate
x = runif(50)
y = runif(50)
pts = data.frame(x = x, y = y)
dom = st_as_sf(pts, coords = c("x","y"), crs = NA, agr = "constant")

rad = 0.5
spatial_bisquare(cbind(x,y), knots, rad)
spatial_bisquare(dom, knots, rad)

# Plot the knots and the points at which we evaluated the basis
plot(knots[,1], knots[,2], pch = 4, cex = 1.5, col = "red")
points(x, y, cex = 0.5)

# Draw a circle representing the basis' radius around one of the knot points
tseq = seq(0, 2*pi, length=100) 
coords = cbind(rad * cos(tseq) + seq_x[2], rad * sin(tseq) + seq_y[2])
lines(coords, col = "red")