Movatterモバイル変換

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 guideon 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 CommunitySurvey 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 ofsupport. 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 STCOSmodel, under several prespecified covariance structures.

Usage

cov_approx_randwalk(Delta, S)cov_approx_blockdiag(Delta, S)

Arguments

Details

Let\bm{\Sigma} be anN \times N symmetric and positive-definitecovariance matrix and\bm{S} be anN \times r matrix withrankr. The objective is to compute a matrix\bm{K} which minimizesthe Frobenius norm

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

over symmetric positive-definite matrices\bm{C}. Thesolution 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 basisfunction computed from a fine-level support havingn areas, usingT time steps. ThereforeN = n T represents the dimension ofcovariance for the fine-level support.

We provide functions to handle some possible structures for targetcovariance 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 ann \times n matrix.

• cov_approx_randwalk assumes\bm{\Sigma} is based on theautocovariance 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 acrosst, so that

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

The block structure is used to reduce the computational burden, asNmay be large.

Deviance Information Criterion

Description

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

Usage

DIC(object, ...)

Arguments

Value

A numeric value of the DIC.

Shapes and ACS estimates for Boone County, MO.

Description

Ansf object with ACS estimates for:

• Boone County, Missouri

• Table B19013

• Block group level geography

• Years 2013 - 2017

Usage

acs5_2013acs5_2014acs5_2015acs5_2016acs5_2017

Format

sf objects.

An object of classsf (inherits fromdata.frame) with 87 rows and 9 columns.

An object of classsf (inherits fromdata.frame) with 87 rows and 9 columns.

An object of classsf (inherits fromdata.frame) with 85 rows and 9 columns.

An object of classsf (inherits fromdata.frame) with 87 rows and 9 columns.

An object of classsf (inherits fromdata.frame) with 87 rows and 9 columns.

Details

This dataset was assembled using shapefiles from thetigris packageand ACS estimates from the American FactFinder website on 2/28/2019.Seedata-prep-aff.R in the Columbia example code for details.American FactFinder has since been deprecated, and similar data areavailable athttp://data.census.gov.

Sparse adjacency matrix between two sets of areas.

Description

A convenience function to convert output fromsf::st_touchesto a sparse matrix as defined in theMatrix package.

Usage

adjacency_matrix(dom)

Arguments

Details

Returns a matrixA whose (i,j)th entry contains a 1 ifareal unitsdom[i,] anddom[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 ansf orsfc object, or as amatrix of points.

• If ansf orsfc object is provided forknots,rthree-dimensionalPOINT entries are expected inst_geometry(knots).Otherwise,knots will be interpreted as anr \times 3 numeric matrix.

• Ifknots is ansf orsfc object, it is checkedto ensure the coordinate system matchesdom.

For each areaA in the given domain, and time period\bm{v} = (v_1, \ldots, v_m) compute the basisfunctions

\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},

forj = 1, \ldots, r. Here,\varphi_j{(\bm{u},v)}representspacetime_bisquare basis functions defined at the pointlevel using\bm{c}_j,g_j,w_s, andw_t.

The basis requires an integration which may be computed using oneof two methods. Themc 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 froma uniform distribution on areaA. Therect method usesa 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 boxst_bbox(A) is divided evenly into a grid ofn_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)threctangle, fora = 1, \ldots, n_x andb = 1, \ldots, n_y.

Due to the treatment ofA_i and\bm{c}_j as objects in aEuclidean space, this basis is more suitable for coordinates from a mapprojection than coordinates based on a globe representation.

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

• method specifies computation method:mc orrect.Default ismc.

• mc_reps is number of repetitions to use formc.Default is 1000.

• nx is number of x-axis points to use forrectmethod. Default is 50.

• ny is number of y-axis points to use forrectmethod. Default is 50.

• report_period is an integer; print a message with progress eachtime this many areas are processed. Default isInf so that messageis suppressed.

• verbose is a logical; ifTRUE print descriptivemessages about the computation. Default isFALSE.

• mc_sampling_factor is a positive number; an oversampling factorused to computeblocksize in therdomain function. I.e.,blocksize = ceiling(mc_sampling_factor * mc_reps). Defaultis 1.2.

Value

A sparsen \times r matrix whoseith rowis\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 pointsseq_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 evaluateshape1 = 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.5period = 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 pointstseq = 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 ansf orsfc object, or as amatrix of points.

• If ansf orsfc object is provided forknots,rtwo-dimensionalPOINT entries are expected inst_geometry(knots).Otherwise,knots will be interpreted as anr \times 2 numeric matrix.

• Ifknots is ansf orsfc object, it is checkedto ensure the coordinate system matchesdom.

For each areaA 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}

forj = 1, \ldots, r. Here,\varphi_j(\bm{u}) representspatial_bisquare basis functions defined at the point levelusing\bm{c}_j andw.

The basis requires an integration which may be computed using oneof two methods. Themc 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 froma uniform distribution on areaA. Therect method usesa 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 boxst_bbox(A) is divided evenly into a grid ofn_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)threctangle, fora = 1, \ldots, n_x andb = 1, \ldots, n_y.

Due to the treatment ofA_i and\bm{c}_j as objects in aEuclidean space, this basis is more suitable for coordinates from a mapprojection than coordinates based on a globe representation.

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

• method specifies computation method:mc orrect.Default ismc.

• mc_reps is number of repetitions to use formc.Default is 1000.

• nx is number of x-axis points to use forrectmethod. Default is 50.

• ny is number of y-axis oints to use forrectmethod. Default is 50.

• report_period is an integer; print a message with progress eachtime this many areas are processed. Default isInf so that messageis suppressed.

• verbose is a logical; ifTRUE print descriptivemessages about the computation. Default isFALSE.

• mc_sampling_factor is a positive number; an oversampling factorused to computeblocksize in therdomain function. I.e.,blocksize = ceiling(mc_sampling_factor * mc_reps). Defaultis 1.2.

Value

A sparsen \times r matrix whoseith rowis\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 pointsseq_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 evaluateshape1 = 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.5areal_spatial_bisquare(dom, knots, rad)areal_spatial_bisquare(dom, knots_sf, rad)# Plot the knots and the points at which we evaluated the basisplot(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 pointstseq = 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 them-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 functionsolves them^2 \times m^2 system

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

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

Value

An arrayGamma of dimensionc(m, m, lag_max + 1),where the sliceGamma[,,h] represents the autocovariance at lagh = 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)) / 2Sigma = diag(x = 2, nrow = 3)autocov_VAR1(A, Sigma, lag_max = 5)

CAR Precision Matrix

Description

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

Usage

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

Arguments

Details

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

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

Ifscale isTRUE, return a scaled version

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

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

Typically in a modeling context, the precision matrix will bemultiplied by a scaling parameter; e.g., a CAR model forrandom 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\}.

whereq = \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

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

Usage

columbia_neighbs

Format

Ansf 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 meanE(Y_i) for each observationi = 1, \ldots, n, for each saved draw.

• predict drawsY_i for each observationi = 1, \ldots, n, using the parameter values for each savedGibbs sampler draw.

Value

gibbs_stcos returns anstcos object which containsdraws from the sampler. Helper functions take this object as an inputand 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 Matlablicols function given inthis MathWorks forum post.

Value

Xsub contains the extracted columns ofX andidxcontains the indices (into X) of those columns. The elapsed time is stored inelapsed.sec.

Examples

x = 0:19 %% 3 + 1Z = 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_hatsig2xi_hat = mle_out$sig2xi_hatmu_hat = mle_out$mu_hat## End(Not run)

Matrix of overlaps between two sets of areas.

Description

A convenience function to convert output fromsf::st_intersectionto a sparse matrix as defined in theMatrix package.

Usage

overlap_matrix(dom1, dom2, proportion = TRUE)

Arguments

Details

Returns a matrixH whose (i,j)th entry represent the area of the overlapbetween areal unitsdom1[i,] anddom2[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 tosf::st_sample which draws uniformly distributedpoints using a simple accept-reject method.

Usage

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

Arguments

Details

Draws a sample ofblocksize points uniformly from a bounding box ondom, and accepts only the points which belong todom. Thisyields a uniform sample ondom. The process is repeated untilnaccepted draws are obtained, or until it has been attempteditmaxtimes. Ifitmax iterations are reached without acceptingndraws, an error is thrown.

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

Performance will degrade when areal units have small area relative to theirbounding box, as many candidate points may need to be discarded. Forexample, this will occur ifdom contains a set of small scatteredislands in an ocean. In this case, it would be more efficient to samplefrom each island at a time.

Value

Ansf 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:

• Bothdom andknots may be provided as eithersf orsfc objects, or as matrices of points.

• If ansf orsfc object is provided fordom,nthree-dimensionalPOINT entries are expected inst_geometry(dom).Otherwise,dom will be interpreted as ann \times 3 numeric matrix.

• If ansf orsfc object is provided forknots,rthree-dimensionalPOINT entries are expected inst_geometry(knots).Otherwise,knots will be interpreted as anr \times 3 numeric matrix.

• If bothdom andknots_s are given assf orsfc 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 \cdotI(\Vert \bm{u} - \bm{c}_j \Vert \leq w_s) \cdotI(|v - g_j| \leq w_t)

forj = 1, \ldots, r.

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

Value

A sparsen \times r matrix whoseith rowis

\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 pointsseq_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 evaluatex = 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.5spacetime_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 basisplot(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 pointstseq = 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:

• Bothdom andknots may be provided as eithersf orsfc objects, or as matrices of points.

• If ansf orsfc object is provided fordom,ntwo-dimensionalPOINT entries are expected inst_geometry(dom).Otherwise,dom will be interpreted as ann \times 2 numeric matrix.

• If ansf orsfc object is provided forknots,rtwo-dimensionalPOINT entries are expected inst_geometry(knots).Otherwise,knots will be interpreted as anr \times 2 numeric matrix.

• If bothdom andknots are given assf orsfc 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 \cdotI(\Vert \bm{u} - \bm{c}_j \Vert \leq w)

forj = 1, \ldots, r.

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

Value

A sparsen \times r matrix whoseith rowis\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 pointsseq_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 evaluatex = 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.5spatial_bisquare(cbind(x,y), knots, rad)spatial_bisquare(dom, knots, rad)# Plot the knots and the points at which we evaluated the basisplot(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 pointstseq = seq(0, 2*pi, length=100) coords = cbind(rad * cos(tseq) + seq_x[2], rad * sin(tseq) + seq_y[2])lines(coords, col = "red")