SpatialBSS (version 0.16-0)

robsbss: Robust Spatial Blind Source Separation

Description

robsbss is a robust variant of sbss. It estimates the unmixing matrix assuming a spatial blind source separation model by jointly diagonalizing the Hettmansperger-Randles scatter matrix and one/many generalized local sign covariance matrices. These local generalized sign covariance matrices are determined by spatial kernel functions and radial functions. Three types of such kernel functions and three types of radial functions are supported.

Usage

robsbss(x, ...)
# S3 method for default
robsbss(x, coords, kernel_type = c('ring', 'ball', 'gauss'), 
     kernel_parameters, lcov = c('norm', 'winsor', 'qwinsor'), 
     ordered = TRUE, kernel_list = NULL, ...)
# S3 method for SpatialPointsDataFrame
robsbss(x, ...)
# S3 method for sf
robsbss(x, ...)

Value

robsbss returns a list of class 'sbss' with the following entries:

s

object of class(x) containing the estimated source random field.

coords

coordinates of the observations. Is NULL if x was a matrix and the argument kernel_list was not NULL at the robsbss call.

w

estimated unmixing matrix.

weights

numeric vector of length(n) giving the weights for each observation for the robust local covariance estimation.

w_inv

inverse of the estimated unmixing matrix.

pevals

(pseudo-)eigenvalues for each latent field entry.

d

matrix of stacked (jointly) diagonalized local covariance matrices with dimension c(length(kernel_parameters)*p,p) for 'ball' and 'gauss' kernel or c( (length(kernel_parameters) / 2)*p,p) for 'ring' kernel.

diags

matrix of dimension c(length(kernel_parameters),p) where the rows contain the diagonal of the diagonalized local autocovariance matrices.

x_mu

robustly estimated columnmeans of x.

cov_inv_sqrt

square root of the inverse sample covariance matrix of x.

Arguments

x

either a numeric matrix of dimension c(n, p) where the p columns correspond to the entries of the random field and the n rows are the observations, an object of class SpatialPointsDataFrame or an object of class sf.

coords

a numeric matrix of dimension c(n,2) where each row represents the coordinates of a point in the spatial domain. Only needed if x is a matrix and the argument kernel_list is NULL.

kernel_type

a string indicating which kernel function to use. Either 'ring' (default), 'ball' or 'gauss'.

kernel_parameters

a numeric vector that gives the parameters for the kernel function. At least length of one for 'ball' and 'gauss' or two for 'ring' kernel, see details.

lcov

a string indicating which radial function or type of robust local covariance matrix to use. Either 'norm' (default), 'winsor' or 'qwinsor'. See also
local_gss_covariance_matrix for details.

ordered

logical. If TRUE the entries of the latent field are ordered by the sum of squared (pseudo-)eigenvalues of the diagonalized local covariance matrix/matrices. Default is TRUE.

kernel_list

a list of spatial kernel matrices with dimension c(n,n), see details. Usually computed by the function spatial_kernel_matrix.

...

further arguments for the fast real joint diagonalization algorithm that jointly diagonalizes the local covariance matrices. See details and frjd.

Details

robsbss is a robust variant of sbss which uses Hettmansperger-Randles (HR) location and scatter estimates (Hettmansperger & Randles, 2002) for whitening (see white_data for details) and jointly diagonalizes HR scatter matrix and generalized local sign matrices to estimate the unmixing matrix. The generalized local sign matrices are determined by radial functions \(w(l_i)\), where \(l_i = ||x(s_i)-T(x)||\) and \(T(x)\) is HR location estimator, and kernel functions \(f(d_{i,j})\), where \(d_{i,j}=||s_i - s_j||\). Generalized local sign covariance (gLSCM) matrix is then calculated as

$$gLSCM(f,w) = 1/(n F^{1/2}_{f,n}) \sum_{i,j} f(d_{i,j}) w(l_i)w(l_j)(x(s_i)-T(x)) (x(s_j)-T(x))'$$ with $$ F_{f,n} = 1 / n \sum_{i,j} f^2(d_{i,j}).$$ Three radial functions (Raymaekers & Rousseeuw, 2019) \(w(l_i)\) are implemented, the parameter lcov defines which is used:

• 'norm': $$ w(l_i) = 1/l_i $$

• 'winsor': $$ w(l_i) = Q/l_i $$

• 'qwinsor': $$ w(l_i) = Q^2/l_i^2.$$

The cutoff \(Q\) is defined as \(Q = l_{(h)}\), where \(l_{(h)}\) is \(h\)th order statistic of \(\{l_1, ..., l_n\}\) and \(h = (n + p + 1)/2\). In addition, three kernel functions \(f(d)\) are implemented, the parameter kernel_type defines which is used:

• 'ring': parameters are inner radius \(r_{in}\) and outer radius \(r_{out}\), with \(r_{in} < r_{out}\), and \(r_{in}, r_{out} \ge 0\): $$f(d;r_{in}, r_{out}) = I(r_{in} < d \le r_{out})$$

• 'ball': parameter is the radius \(r\), with \(r \ge 0\): $$f(d;r) = I(d \le r)$$

• 'gauss': Gaussian function where 95% of the mass is inside the parameter \(r\), with \(r \ge 0\): $$f(d;r) = exp(-0.5 (\Phi^{-1}(0.95) d/r)^2).$$

The argument kernel_type determines the used kernel function as presented above, the argument kernel_parameters gives the corresponding parameters for the kernel function. Specifically, if kernel_type equals 'ball' or 'gauss' then kernel_parameters is a numeric vector where each entry corresponds to one parameter. Hence, length(kernel_parameters) local covariance matrices are used. Whereas, if kernel_type equals 'ring', then kernel_parameters must be a numeric vector of even length where subsequently the inner and outer radii must be given (informally: c(r_in1, r_out1, r_in2, r_out2, ...)). In that case length(kernel_parameters) / 2 local covariance matrices are used.

robsbss calls spatial_kernel_matrix internally to compute a list of c(n,n) kernel matrices based on the parameters given, where each entry of those matrices corresponds to \(f(d_{i,j})\). Alternatively, such a list of kernel matrices can be given directly to the function robsbss via the kernel_list argument. This is useful when robsbss is called numerous times with the same coordinates/kernel functions as the computation of the kernel matrices is then done only once prior the actual robsbss calls. For details see also spatial_kernel_matrix.

If more than one generalized local sign covariance matrix is used robsbss jointly diagonalizes these matrices with the function frjd. ... provides arguments for frjd, useful arguments might be:

• eps: tolerance for convergence.

• maxiter: maximum number of iterations.

References

Hettmansperger, T. P., & Randles, R. H. (2002). A practical affine equivariant multivariate median. Biometrika, 89 , 851-860. tools:::Rd_expr_doi("10.1093/biomet/89.4.851").

Raymaekers, J., & Rousseeuw, P. (2019). A generalized spatial sign covariance matrix. Journal of Multivariate Analysis, 171 , 94-111. tools:::Rd_expr_doi("10.1016/j.jmva.2018.11.010").

Sipila, M., Muehlmann, C. Nordhausen, K. & Taskinen, S. (2024). Robust second order stationary spatial blind source separation using generalized sign matrices. Spatial Statistics, 59, 100803. tools:::Rd_expr_doi("10.1016/j.spasta.2023.100803").

See Also

spatial_kernel_matrix, local_gss_covariance_matrix, sp, sf, frjd

Examples

# simulate coordinates
coords <- runif(1000 * 2) * 20
dim(coords) <- c(1000, 2)
coords_df <- as.data.frame(coords)
names(coords_df) <- c("x", "y")
# simulate random field
if (!requireNamespace('gstat', quietly = TRUE)) {
  message('Please install the package gstat to run the example code.')
} else {
  library(gstat)
  model_1 <- gstat(formula = z ~ 1, locations = ~ x + y, dummy = TRUE, beta = 0, 
                   model = vgm(psill = 0.025, range = 1, model = 'Exp'), nmax = 20)
  model_2 <- gstat(formula = z ~ 1, locations = ~ x + y, dummy = TRUE, beta = 0, 
                   model = vgm(psill = 0.025, range = 1, kappa = 2, model = 'Mat'), 
                   nmax = 20)
  model_3 <- gstat(formula = z ~ 1, locations = ~ x + y, dummy = TRUE, beta = 0, 
                   model = vgm(psill = 0.025, range = 1, model = 'Gau'), nmax = 20)
  field_1 <- predict(model_1, newdata = coords_df, nsim = 1)$sim1
  field_2 <- predict(model_2, newdata = coords_df, nsim = 1)$sim1
  field_3 <- predict(model_3, newdata = coords_df, nsim = 1)$sim1
  field <- cbind(field_1, field_2, field_3)
  # Generate 5 % local outliers to data
  field_cont <- gen_loc_outl(field, coords, radius = 2, 
                             swap_order = "regular")[,1:3]
  X <- as.matrix(field_cont)
  
  # apply sbss with three ring kernels
  kernel_parameters <- c(0, 1, 1, 2, 2, 3)
  robsbss_result <- 
    robsbss(X, coords, kernel_type = 'ring', kernel_parameters = kernel_parameters)
  
  # print object
  print(robsbss_result)
  
  # plot latent field
  plot(robsbss_result, colorkey = TRUE, as.table = TRUE, cex = 1)
  
  # predict latent fields on grid
  predict(robsbss_result, colorkey = TRUE, as.table = TRUE, cex = 1)
  
  # unmixing matrix
  w_unmix <- coef(robsbss_result)
}