# smooth.FEM: Spatial regression with differential regularization In fdaPDE: Statistical Analysis of Functional and Spatial Data, Based on Regression with PDE Regularization

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

• `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 `j`th 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), 681-703.

• Azzimonti, L., Sangalli, L. M., Secchi, P., Domanin, M., Nobile, F. (2015). Blood flow velocity field estimation via spatial regression with PDE penalization. Journal of the American Statistical Association, 110(511), 1057-1071.

## 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
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
 "C++ Code Execution"
Smoothing: isotropic and stationary case
 "C++ Code Execution"
[,1]
[1,] 1.887878
Smoothing: isotropic and stationary case
 "C++ Code Execution"
Smoothing: anysotropic and non-stationary case
 "C++ Code Execution"
Smoothing: anysotropic and non-stationary case
 "C++ Code Execution"
```

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