Spatial regression with differential regularization: stationary and isotropic case (Laplacian)

Share:

Description

This function implements a spatial regression model with differential regularization; isotropic and stationary case. In particular, the regularizing term involves the Laplacian of the spatial field. Space-varying covariates can be included in the model. The technique accurately handle data distributed over irregularly shaped domains. Moreover, various conditions can be imposed at the domain boundaries.

Usage

1
2
smooth.FEM.basis(locations = NULL, observations, FEMbasis, lambda, 
       covariates = NULL, BC = NULL, GCV = FALSE, CPP_CODE = TRUE)

Arguments

locations

A #observations-by-2 matrix where each row specifies the spatial coordinates x and y of the corresponding observations in the vector observations. This parameter can be NULL. In this case the spatial coordinates of the corresponding observations are assigned as specified in observations.

observations

A vector of length #observations with the observed data values over the domain. The locations of the observations can be specified with the locations argument. Otherwise if only the vector of observations is given, these are consider to be located in the corresponding node in the table nodes of the mesh. In this last case, an NA value in the observations vector indicates that there is no observation associated to the corresponding node.

FEMbasis

A FEMbasis object describing the Finite Element basis, as created by create.FEM.basis.

lambda

A scalar or vector of smoothing parameters.

covariates

A #observations-by-#covariates matrix where each row represents the covariates associated with the corresponding observed data value in observations.

BC

A list with two vectors: BC_indices, a vector with the indices in nodes of boundary nodes where a Dirichlet Boundary Condition should be applied; BC_values, a vector with the values that the spatial field must take at the nodes indicated in BC_indices.

GCV

Boolean. If TRUE the following quantities are computed: the trace of the smoothing matrix, the estimated error standard deviation, and the Generalized Cross Validation criterion, for each value of the smoothing parameter specified in lambda.

CPP_CODE

Boolean. If TRUE the computation relies on the C++ implementation of the algorithm. This usually ensures a much faster computation.

Value

A list with the following variables:

fit.FEM

A FEM object that represents the fitted spatial field.

PDEmisfit.FEM

A FEM object that represents the Laplacian of the estimated spatial field.

beta

If covariates is not NULL, a matrix with number of rows equal to the number of covariates and numer of columns equal to length of lambda. The jth column represents the vector of regression coefficients when the smoothing parameter is equal to lambda[j].

edf

If GCV is TRUE, a scalar or vector with the trace of the smoothing matrix for each value of the smoothing parameter specified in lambda.

stderr

If GCV is TRUE, a scalar or vector with the estimate of the standard deviation of the error for each value of the smoothing parameter specified in lambda.

GCV

If GCV is TRUE, a scalar or vector with the value of the GCV criterion for each value of the smoothing parameter specified in lambda.

References

Sangalli, L.M., Ramsay, J.O. & Ramsay, T.O., 2013. Spatial spline regression models. Journal of the Royal Statistical Society. Series B: Statistical Methodology, 75(4), pp. 681-703.

See Also

smooth.FEM.PDE.basis, smooth.FEM.PDE.sv.basis

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
library(fdaPDE)
## Load the Meuse data and a domain boundary for these data
data(MeuseData)
data(MeuseBorder)
## Create a triangular mesh for these data with the provided boundary and plot it
order=1
mesh <- create.MESH.2D(nodes = MeuseData[,c(2,3)], segments = MeuseBorder, order = order)
plot(mesh)
## Create the Finite Element basis 
FEMbasis = create.FEM.basis(mesh)
## Estimate zync field without using covariates, setting the smoothing parameter to 10^3.5
data = log(MeuseData[,"zinc"])
lambda = 10^3.5
ZincMeuse = smooth.FEM.basis(observations = data, 
                             FEMbasis = FEMbasis, lambda = lambda)
## Plot the estimated spatial field 
plot(ZincMeuse$fit.FEM)
# Now repeat the analysis using as covariates the square root of the log-distance 
# from river \code{sqrt(dist.log(m))} and the altitude \code{elev}
desmat = matrix(1,nrow=nrow(MeuseData),ncol=2)
desmat[,1] = sqrt(MeuseData[,"dist.log(m)"])
desmat[,2] = MeuseData[,"elev"]
ZincMeuseCovar = smooth.FEM.basis(observations = data, covariates = desmat, 
                                   FEMbasis = FEMbasis, lambda = lambda)
# Plot of the non parametric part (f) of the regression model y_i = beta_1 x_i1 + beta_2 x_i2 + f
plot(ZincMeuseCovar$fit.FEM)
# Print covariates' regression coefficients
print(ZincMeuseCovar$beta)