smooth.FEM: Spatial regression with differential regularization

Description Usage Arguments Value References Examples

View source: R/smoothing.R

Description

This function implements a spatial regression model with differential regularization. The regularizing term involves a Partial Differential Equation (PDE). In the simplest case the PDE involves only the Laplacian of the spatial field, that induces an isotropic smoothing. When prior information about the anisotropy or non-stationarity is available the PDE involves a general second order linear differential operator with possibly space-varying coefficients. The technique accurately handle data distributed over irregularly shaped domains. Moreover, various conditions can be imposed at the domain boundaries.

Usage

1
2
3
smooth.FEM(locations = NULL, observations, FEMbasis, lambda, 
     covariates = NULL, PDE_parameters=NULL, incidence_matrix = NULL, 
     BC = NULL, GCV = FALSE, GCVmethod = "Stochastic", nrealizations = 100)

Arguments

locations

A #observations-by-2 matrix in the 2D case and #observations-by-3 matrix in the 2.5D case, where each row specifies the spatial coordinates x and y (and z in 2.5D) of the corresponding observation in the vector observations. If the locations of the observations coincide with (or are a subset of) the nodes of the mesh in the FEMbasis, leave the parameter locations = NULL for a faster implementation.

observations

A vector of length #observations with the observed data values over the domain. If the locations argument is left NULL the vector of the observations have to be of length #nodes of the mesh in the FEMbasis. In this case, each observation is associated to the corresponding node in the mesh. If the observations are observed only on a subset of the mesh nodes, fill with NA the values of the vector observations in correspondence of unobserved data.

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 and each column is a different covariate.

PDE_parameters

A list specifying the parameters of the PDE in the regularizing term. Default is NULL, i.e. regularization is by means of the Laplacian (stationary, isotropic case). If the coefficients of the PDE are constant over the domain PDE_parameters must contain:

  • K, a 2-by-2 matrix of diffusion coefficients. This induces an anisotropic smoothing with a preferential direction that corresponds to the first eigenvector of the diffusion matrix K;

  • b, a vector of length 2 of advection coefficients. This induces a smoothing only in the direction specified by the vector b;

  • c, a scalar reaction coefficient. c induces a shrinkage of the surface to zero.

If the coefficients of the PDE are space-varying PDE_parameters must contain:

  • K, a function that for each spatial location in the spatial domain (indicated by the vector of the 2 spatial coordinates) returns a 2-by-2 matrix of diffusion coefficients. The function must support recycling for efficiency reasons, thus if the input parameter is a #point-by-2 matrix, the output should be an array with dimensions 2-by-2-by-#points.

  • b, a function that for each spatial location in the spatial domain returns a vector of length 2 of transport coefficients. The function must support recycling for efficiency reasons, thus if the input parameter is a #point-by-2 matrix, the output should be a matrix with dimensions 2-by-#points;

  • c, a function that for each spatial location in the spatial domain returns a scalar reaction coefficient. The function must support recycling for efficiency reasons, thus if the input parameter is a #point-by-2 matrix, the output should be a vector with length #points;

  • u, a function that for each spatial location in the spatial domain returns a scalar reaction coefficient. u induces a reaction effect. The function must support recycling for efficiency reasons, thus if the input parameter is a #point-by-2 matrix, the output should be a vector with length #points.

For 2.5D, only the Laplacian is available (PDE_parameters=NULL).

incidence_matrix

A #regions-by-#triangles matrix where the element (i,j) equals 1 if the j-th triangle is in the i-th region and 0 otherwise. This is needed only for areal data. In case of pointwise data, this parameter is set to NULL.

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.

GCVmethod

This parameter is considered only when GCV=TRUE. It can be either "Exact" or "Stochastic". If set to "Exact" the algoritm performs an exact (but possibly slow) computation of the GCV index. If set to "Stochastic" the GCV is approximated by a stochastic algorithm.

nrealizations

This parameter is considered only when GCV=TRUE and GCVmethod = "Stochastic". It is a positive integer that represents the number of uniform random variables used in stochastic GCV computation.

Value

A list with the following variables:

References

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
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
library(fdaPDE)

#### No prior information about anysotropy/non-stationarity (laplacian smoothing) ####
data(horseshoe2D)
boundary_nodes = horseshoe2D$boundary_nodes
boundary_segments = horseshoe2D$boundary_segments
locations = horseshoe2D$locations

mesh = create.mesh.2D(nodes = rbind(boundary_nodes, locations), segments = boundary_segments)
FEMbasis = create.FEM.basis(mesh)
lambda = 10^-1
# no covariate
data = fs.test(mesh$nodes[,1], mesh$nodes[,2]) + rnorm(nrow(mesh$nodes), sd = 0.5)

solution = smooth.FEM(observations = data, FEMbasis = FEMbasis, lambda = lambda)
plot(solution$fit.FEM)

# with covariates
covariate = covs.test(mesh$nodes[,1], mesh$nodes[,2])
data = fs.test(mesh$nodes[,1], mesh$nodes[,2]) +
             2*covariate + rnorm(nrow(mesh$nodes), sd = 0.5)

solution = smooth.FEM(observations = data, covariates = covariate, 
                      FEMbasis = FEMbasis, lambda = lambda)
# beta estimate:
solution$beta
# non-parametric estimate:
plot(solution$fit.FEM)

# Choose lambda with GCV:
lambda = 10^(-2:2)
solution = smooth.FEM(observations = data, 
                            covariates = covariate, 
                            FEMbasis = FEMbasis, 
                            lambda = lambda, 
                            GCV = TRUE)
bestLambda = lambda[which.min(solution$GCV)]


#### Smoothing with prior information about anysotropy/non-stationarity and boundary conditions ####
# See Azzimonti et al. for reference to the current exemple
data(quasicircle2D)
boundary_nodes = quasicircle2D$boundary_nodes
boundary_segments = quasicircle2D$boundary_segments
locations = quasicircle2D$locations
data = quasicircle2D$data

mesh = create.mesh.2D(nodes = rbind(boundary_nodes, locations), segments = boundary_segments)
FEMbasis = create.FEM.basis(mesh)
lambda = 10^-2

# Set the PDE parameters
R = 2.8
K1 = 0.1
K2 = 0.2
beta = 0.5
K_func<-function(points)
{
  output = array(0, c(2, 2, nrow(points)))
  for (i in 1:nrow(points))
    output[,,i] = 10*rbind(c(points[i,2]^2 + K1*points[i,1]^2 + 
                             K2*(R^2 - points[i,1]^2 - points[i,2]^2),
                             (K1-1)*points[i,1]*points[i,2]),
                           c((K1-1)*points[i,1]*points[i,2],
                             points[i,1]^2 + K1*points[i,2]^2 + 
                             K2*(R^2 - points[i,1]^2 - points[i,2]^2)))
  output
}

b_func<-function(points)
{
  output = array(0, c(2, nrow(points)))
  for (i in 1:nrow(points))
    output[,i] = 10*beta*c(points[i,1],points[i,2])
  output
}

c_func<-function(points)
{
  rep(c(0), nrow(points))
}

u_func<-function(points)
{
  rep(c(0), nrow(points))
}
PDE_parameters = list(K = K_func, b = b_func, c = c_func, u = u_func)

# Set the boundary conditions
BC = NULL
BC$BC_indices = which(mesh$nodesmarkers == 1) # b.c. on the complete boundary
BC$BC_values = rep(0,length(BC$BC_indices)) # homogeneus b.c.

# Since the data locations are a subset of the mesh nodes for a faster solution use:
dataNA = rep(NA, FEMbasis$nbasis)
dataNA[mesh$nodesmarkers == 0] = data

solution = smooth.FEM(observations = dataNA, 
                            FEMbasis = FEMbasis, 
                            lambda = lambda, 
                            PDE_parameters = PDE_parameters, 
                            BC = BC)
plot(solution$fit.FEM)
image(solution$fit.FEM)

#### Smoothing with areal data ####
# See Azzimonti et al. for reference to the current exemple
data(quasicircle2Dareal)
incidence_matrix = quasicircle2Dareal$incidence_matrix
data = quasicircle2Dareal$data
mesh = quasicircle2Dareal$mesh

FEMbasis = create.FEM.basis(mesh)
lambda = 10^-4

# Set the PDE parameters
R = 2.8
K1 = 0.1
K2 = 0.2
beta = 0.5
K_func<-function(points)
{
  output = array(0, c(2, 2, nrow(points)))
  for (i in 1:nrow(points))
    output[,,i] = 10*rbind(c(points[i,2]^2 + K1*points[i,1]^2 + 
                             K2*(R^2 - points[i,1]^2 - points[i,2]^2),
                             (K1-1)*points[i,1]*points[i,2]),
                           c((K1-1)*points[i,1]*points[i,2],
                             points[i,1]^2 + K1*points[i,2]^2 + 
                             K2*(R^2 - points[i,1]^2 - points[i,2]^2)))
  output
}

b_func<-function(points)
{
  output = array(0, c(2, nrow(points)))
  for (i in 1:nrow(points))
    output[,i] = 10*beta*c(points[i,1],points[i,2])
  output
}

c_func<-function(points)
{
  rep(c(0), nrow(points))
}

u_func<-function(points)
{
  rep(c(0), nrow(points))
}
PDE_parameters = list(K = K_func, b = b_func, c = c_func, u = u_func)

# Set the boundary conditions
BC = NULL
BC$BC_indices = which(mesh$nodesmarkers == 1) # b.c. on the complete boundary
BC$BC_values = rep(0,length(BC$BC_indices)) # homogeneus b.c.

solution = smooth.FEM(observations = data, 
                            incidence_matrix = incidence_matrix,
                            FEMbasis = FEMbasis, 
                            lambda = lambda, 
                            PDE_parameters = PDE_parameters, 
                            BC = BC)
plot(solution$fit.FEM)
image(solution$fit.FEM)

Example output

Loading required package: geometry
Loading required package: rgl
Loading required package: Matrix
Loading required package: plot3D
Loading required package: plot3Drgl
Warning messages:
1: In rgl.init(initValue, onlyNULL) : RGL: unable to open X11 display
2: 'rgl.init' failed, running with 'rgl.useNULL = TRUE'. 
3: no DISPLAY variable so Tk is not available 
Smoothing: isotropic and stationary case
[1] "C++ Code Execution"
Smoothing: isotropic and stationary case
[1] "C++ Code Execution"
         [,1]
[1,] 1.887878
Smoothing: isotropic and stationary case
[1] "C++ Code Execution"
Smoothing: anysotropic and non-stationary case
[1] "C++ Code Execution"
Smoothing: anysotropic and non-stationary case
[1] "C++ Code Execution"

fdaPDE documentation built on July 2, 2020, 2:22 a.m.