Efficient Estimation of Partially Linear Models for Spatial Data over Complex Domains

arXiv:1605.08737v2 [] 3 Jun 2016

Li Wanga , Guannan Wangb , Min-Jun Laic and Lei Gaoa
a Iowa

State University, b College of William & Mary and c University of Georgia

Abstract: In this paper, we study the estimation of partially linear models for spatial data distributed over complex domains. We use bivariate splines over triangulations to represent the
nonparametric component on an irregular two-dimensional domain. The proposed method is formulated as a constrained minimization problem which does not require constructing finite elements
or locally supported basis functions. Thus, it allows an easier implementation of piecewise polynomial representations of various degrees and various smoothness over an arbitrary triangulation.
Moreover, the constrained minimization problem is converted into an unconstrained minimization
via a QR decomposition of the smoothness constraints, which leads to a penalized least squares
method to estimate the model. The estimators of the parameters are proved to be asymptotically
normal under some regularity conditions. The estimator of the bivariate function is consistent, and
its rate of convergence is also established. The proposed method enables us to construct confidence
intervals and permits inference for the parameters. The performance of the estimators is evaluated
by two simulation examples and by a real data analysis.
Key words and phrases: Bivariate splines, Penalty, Semiparametric regression, Spatial data, Triangulation.

1. Introduction
In many geospatial studies, spatially distributed covariate information is available. For
example, geographic information systems may contain measurements obtained from satellite
images at some locations. These spatially explicit data can be useful in the construction
and estimation of regression models, but sometimes they are distributed over irregular twodimensional (2-D) domains that may have complex boundaries or holes inside. It is well
known that many conventional smoothing tools suffer from the problem of “leakage” across
the complex domains, which refers to the poor estimation over difficult regions by smoothing inappropriately across boundary features, such as peninsulas; see excellent discussions in
Ramsay (2002) and Wood, Bravington and Hedley (2008). In this paper, we propose to use
bivariate splines (smooth piecewise polynomial functions over a triangulation of the domain
of interest) to model spatially explicit datasets which enable us to overcome the “leakage”
problem and find the estimation more accurately.
Address for correspondence: Li Wang, Department of Statistics and the Statistical Laboratory, Iowa State
University, Ames, IA, USA. Email: lilywang@iastate.edu

2

Li Wang, Guannan Wang, Min-Jun Lai and Lei Gao

We focus here on the partially linear model (PLM), popularized by Hardle, Liang and
Gao (2000), for data randomly distributed over 2-D domains. To be more specific, let Xi =
(Xi1 , Xi2 )T be the location of i-th point, i = 1, . . . , n, which ranges over a bounded domain
Ω ⊆ R2 of arbitrary shape, for example, a domain with polygon boundary. Let Yi be the
response variable and Zi = (Zi1 , . . . , Zip )T be the predictors at location Xi . Suppose that
{(Zi , Xi , Yi )}ni=1 satisfies the following model
Y i = ZT
i β + g (Xi ) + i ,

i = 1, . . . , n,

(1.1)

where β = (β1 , . . . , βp )T are unknown parameters, g(·) is some unknown but smooth bivariate
function, and i ’s are i.i.d random noises with E (i ) = 0 and Var (i ) = σ 2 . Each i is
independent of Xi and Zi . In many situations, our main interest is in estimating and making
inference for the regression parameters β, which provides measures of the effect of the covariate
Z after adjusting for the location effect of X.
If g(·) is a univariate function, model (1.1) becomes a typical PLM. In the past three
decades, flexible and parsimonious PLMs have been extensively studied and widely used in
many statistical applications, from biostatistics to econometrics, from engineering to social science; see Chen, Liang and Wang (2011), Huang, Zhang and Zhou (2007), Liang and Li (2009),
Liu, Wang and Liang (2011), Ma, Song and Wang (2013), Ma and Yang (2011), Wang, et al
(2011), Wang, et al (2014), Zhang, Cheng and Liu (2011) for some recent works on PLMs.
When g(·) is a bivariate function, there are two popular estimation tools: bivariate P-splines
Marx and Eilers (2005) and thin plate splines Wood (2003). Later, Xiao, Li and Ruppert (2013)
proposed an efficient sandwich smoother, which has a tensor product structure that simplifies
an asymptotic analysis and can be fast computed. Their method has been applied to quantifying the lifetime circadian rhythm of physical activity Xiao, et al (2015). The application to
spatial data analysis over complex domains, however, has been hampered due to the scarcity
of bivariate smoothing tools that are not only computationally effficient but also theoretically
reliable to solve the problem of “leakage” across the domain. Traditional smoothing methods
in practical data analysis, such as kernel smoothing, wavelet-based smoothing, tensor product
splines and thin plate splines, usually perform poorly for those data, since they do not take
into account the shape of the domain and also smooth across concave boundary regions.
There are several challenges when going from rectangular domains to irregular domains
with complex boundaries or holes. Some efforts have recently been devoted to studying the
smoothing over irregular domains, and significant progress has been made. Most of them are
based on the roughness penalization approach Green and Silverman (1994). To deal with
irregular domains, Wang and Ranalli (2007) applied low-rank thin-plate splines defined as
functions of the geodesic distance instead of the Euclidean distance. Eilers (2006) utilized the

Efficient Estimation of Partially Linear Models for Spatial Data over Complex Domains

3

Schwarz-Christoffel transform to convert the complex domains to regular domains. To solve the
“leakage” problem, in a pioneering paper, Ramsay (2002) suggested a penalized least squares
approach with a Laplacian penalty and transformed the problem to that of solving a system of
partial differential equations (PDEs). Wood, Bravington and Hedley (2008) provided an elegant
solution and developed the soap film smoothing estimator for smoothing over difficult regions
that can be represented by a low-rank basis and one or two quadratic penalties. Recently,
Sangalli, Ramsay and Ramsay (2013) extended the method in Ramsay (2002) to the PLMs,
which allows for spatially distributed covariate information to be in the models. The data
smoothing problem in Sangalli, Ramsay and Ramsay (2013) is solved using finite element
method (FEM), a method mainly developed and used to solve PDEs. Although their method
is practically useful, the theoretical properties of the estimation are not investigated.
In this paper, we tackle the estimation problem differently from Sangalli, Ramsay and
Ramsay (2013). Instead of using FEM, we approximate the nonparametric function g(·) by
using the bivariate splines over triangulations in Lai and Schumaker (2007). An important
feature of this approach is that it uses splines for applications without constructing locally
supported splines or finite elements and without computing the dimension. This method has
been shown to be more efficient and flexible than the conventional FEM in data fitting problems
and solving PDEs; see Awanou, Lai and Wenston (2005), Ettinger, Guillas and Lai (2015),
Guillas and Lai (2010), Lai and Wang (2013) and Liu, Guillas and Lai (2015). For example,
the users can choose spline functions of flexible degrees and various smoothness across any
given domain. Another advantage is that the linear systems arising in this approach are more
easy to assemble than those from the finite elements or locally supported spline basis functions.
The linear systems are sparser than that from any macro-FEM method. In addition, due to
the great scalability, the assembling can be computed in parallel.
To the best of the authors’ knowledge, statistical aspects of smoothing for PLMs by using
bivariate splines have not been discussed in the literature so far. This paper presents the
first attempt at investigating the asymptotic properties of the PLMs for data distributed on
a non-rectangular complex region. We study the asymptotic properties of the least squares
estimators of β and g(·) by using bivariate splines over triangulations with a penalty term.
We show that our estimator of β is root-n consistent and asymptotically normal, although the
convergence rate of the estimator of the nonparametric component g(·) is slower than root-n.
A standard error formula for the estimated coefficients is provided and tested to be accurate
enough for practical purposes. Hence, the proposed method enables us to construct confidence
intervals for the regression parameters. We also obtain the convergence rate for the functional
estimator of g(·). We show, by using numerical studies, that our method is competitive with
existing methods such as the soap lm smoother and the thin-plate regression spline.

4

Li Wang, Guannan Wang, Min-Jun Lai and Lei Gao

The rest of the paper is organized as follows. In Section 2, we give a brief review of the
triangulations and propose our estimation method based on penalized bivariate splines. We
also discuss the details on how to choose the penalty parameters. Section 3 is devoted to the
asymptotic analysis of the proposed estimators. Section 4 provides a detailed simulation to
compare several methods in two different scenarios and explores the estimation and prediction
accuracy. Section 5 studies a real dataset on house values over California. Some concluding
remarks are given in Section 6. Technical details are provided in Appendix A.
2. Triangulations and Penalized Spline Estimators
Our estimation method is based on penalized bivariate splines over triangulations. The
idea is to approximate the function g(·) by bivariate splines that are piecewise polynomial
functions over a 2D triangulated domain. We use this approximation to construct least squares
estimators of the linear and nonlinear components of the model with a penalization term. In
the following of this section, we describe the background of triangulations, B-form bivariate
splines and introduce the penalized spline estimators.
2.1. Triangulations
Triangulation is an effective strategy to handle data distribution on irregular regions with
complex boundaries and/or interior holes. Recently, it has attracted substantial recent attention in many applied areas, such as geospatial studies, numerical solutions of PDEs, image
enhancements, and computer aided geometric design. See, for example, the recent comprehensive book by Lai and Schumaker (2007) and the article by Lai (2008).
We use τ to denote a triangle which is a convex hull of three points not located in one line.
A collection 4 = {τ1 , ..., τN } of N triangles is called a triangulation of Ω = ∪N
i=1 τi provided
that if a pair of triangles in 4 intersect, then their intersection is either a common vertex or
a common edge. In general, any kind of polygon shapes can be used for the partition of Ω. In
this paper we consider triangulations of Ω because any polygonal domain of arbitrary shape
can be partitioned into finitely many triangles. In the following, we assume that all Xi s are
inside triangles of 4. That is, they are not on edges or vertices of triangles in 4. Otherwise,
we can simply count them twice or multiple times if any observation is located on an edge or
at a vertex of 4. Given a triangle τ ∈ 4, we let |τ | be its longest edge length, and denote the
size of 4 by |4| = max{|τ |, τ ∈ 4}, i.e., the length of the longest edge of 4. Furthermore, let
ρτ be the radius of the largest circle inscribed in τ . We measure the quality of a triangulation
4 by δ4 = maxτ ∈4 |τ |/ρτ < ∞, which is equivalent to the smallest angle of 4. The study in
Lai and Schumaker (2007) shows that the approximation of spline spaces over 4 is dependent
on δ4 , i.e., the larger the δ4 is, the worse the spline approximation is. In the rest of the paper,
we restrict our attention to the triangulations satisfying δ4 < δ for a positive constant δ.

5

Efficient Estimation of Partially Linear Models for Spatial Data over Complex Domains

2.2. B-form bivariate splines
In this section we give a brief introduction to the bivariate splines. More in-depth description can be found in Lai and Schumaker (2007), Lai (2008), as well as Zhou and Pan (2014)
and the details of the implementation is provided in Awanou, Lai and Wenston (2005). Let
τ = hv1 , v2 , v3 i be a non-degenerate (i.e. with non-zero area) triangle with vertices v1 , v2 ,
and v3 . Then for any point v ∈ R, there is a unique representation in the form
v = b1 v1 + b2 v2 + b3 v3
with b1 + b2 + b3 = 1, where b1 , b2 and b3 are called the barycentric coordinates of the point
v relative to the triangle τ . The Bernstein polynomials of degree d relative to triangle τ is
defined as
τ,d
Bijk
(v) =

d i j k
b b b .
i!j!k! 1 2 3

Then for any τ ∈ 4, we can write the polynomial piece of spline s restricted on τ ∈ 4 as
X

s|τ =

τ,d
τ
,
Bijk
γijk

i+j+k=d
τ , i + j + k = d} are called B-coefficients of s.
where the coefficients γ τ = {γijk

For a nonnegative integer r, let Cr (Ω) be the collection of all r-th continuously differentiable
functions over Ω. Given a triangulation 4, let Srd (4) = {s ∈ Cr (Ω) : s|τ ∈ Pd (τ ), τ ∈ 4} be
a spline space of degree d and smoothness r over triangulation 4, where Pd is the space of all
polynomials of degree less than or equal to d. For notation simplicity, let S = Sr3r+2 (4) for a
fixed smoothness r ≥ 1, and we know that such a spline space has the optimal approximation
order (rate of convergence) for noise-free datasets; see Lai and Schumaker (1998) and Lai and
Schumaker (2007).
For notation simplicity, let {Bξ }ξ∈K be the set of degree-d bivariate Bernstein basis polynomials for S, where K stands for an index set of K Bernstein basis polynomials. Then for any
function s ∈ S, we can represent it by using the following basis expansion:
s(x) =

X

Bξ (x)γξ = B(x)T γ,

(2.1)

ξ∈K

where γ T = (γξ , ξ ∈ K) is the spline coefficient vector. To meet the smoothness requirement of
the splines, we need to impose some linear constraints on the spline coefficients γ in (2.1). We
require that γ satisfies Hγ = 0 with H being the matrix for all smoothness conditions across
shared edges of triangles, which depends on r and the structure of the triangulation.
2.3. Penalized Spline Estimators

6

Li Wang, Guannan Wang, Min-Jun Lai and Lei Gao

To define the penalized spline method, for any direction xj , j = 1, 2, let Dxq j f (x) denote
the q-th order derivative in the direction xj at the point x = (x1 , x2 ). Let
X Z X υ 
(Dxi 1 Dxj 2 f )2 dx1 dx2
Eυ (f ) =
i
τ
τ ∈4

(2.2)

i+j=υ

be the energy functional for a fixed integer υ ≥ 1. Although all partial derivatives up to the
chosen order υ can be included in the penalty of (2.2), for simplicity, in the remaining part
of the paper, we use υ = 2, and one can study the similar problem for general υ ≥ 2. When
υ = 2,
Z
E2 (f ) =
Ω


(Dx21 f )2 + 2(Dx1 Dx2 f )2 + (Dx22 f )2 dx1 dx2 ,

(2.3)

which is similar to the thin-plate spline penalty (Green and Silverman, 1994) except the latter
is integrated over the entire plane R2 . Sangalli, Ramsay and Ramsay (2013) used a different
roughness penalty from (2.3), specifically, they use the integral of the square of the Laplacian
R
of f , that is, λ Ω (Dx21 f + Dx22 f )2 dx1 dx2 . Both forms of penalties are invariant with respect
to Euclidean transformations of spatial co-ordinates, thus, the bivariate smoothing does not
depend on the choice of the coordinate system.
Given λ > 0 and {(Zi , Xi , Yi )}ni=1 , we consider the following minimization problem:
min
s∈S

n
X


Yi − ZT
i β − s (Xi )

2

+ λEυ (s).

(2.4)

i=1

Let Y = (Y1 , . . . , Yn )T be the vector of n observations of the response variable. Denote
by Xn×2 = {(Xi1 , Xi2 )}ni=1 the location design matrix and Zn×p = {(Zi1 , . . . , Zip )}ni=1 the
collection of all covariates. Denote by B the n × K evaluation matrix of Bernstein basis
polynomials whose i-th row is given by BT
i = {Bξ (Xi ), ξ ∈ K}. Then the minimization
problem in (2.4) reduces to

min L(β, γ) = min kY − Zβ − Bγk2 + λγ T Pγ
β,γ

β,γ

subject to Hγ = 0,

(2.5)

where P is the block diagonal penalty matrix satisfying that γ T Pγ = Eυ (Bγ).
To solve the constrained minimization problem (2.5), we first remove the constraint via
QR decomposition of the transpose of the constrain matrix H. Specifically, we write
!
R
1
HT = QR = (Q1 Q2 )
,
R2

(2.6)

where Q is an orthogonal matrix and R is an upper triangle matrix, the submatrix Q1 is the
first r columns of Q, where r is the rank of matrix H, and R2 is a matrix of zeros. It is easy
to see the following result.

7

Efficient Estimation of Partially Linear Models for Spatial Data over Complex Domains

Lemma 1. Let Q1 , Q2 be submatrices as in (2.6). Let γ = Q2 θ for a vector θ of appropriate
size. Then Hγ = 0. On the other hand, if Hγ = 0, then there exists a vector θ such that
γ = Q2 θ.
T
Proof. By (2.6), we have HT = Q1 R1 since R2 = 0. That is, H = RT
1 Q1 . Thus,
T
Hγ =HQ2 θ = RT
1 Q1 Q2 θ = 0

since QT
1 Q2 = 0. On the other hand, if
T
0 = Hγ = RT
1 Q1 γ,

we have QT
1 γ = 0 since R1 is invertible. Thus, γ is in the perpendicular subspace of the space
spanned by the columns of Q1 . That is, γ is in the space spanned by the columns of Q2 . Thus,
there exists a vector θ such that γ = Q2 θ. These complete the proof. 
The problem (2.5), is now converted to a conventional penalized regression problem without
any constraints:

min kY − Zβ − BQ2 θk2 + λ(Q2 θ)T P(Q2 θ) .
β,θ

For a fixed penalty parameter λ, we have
b
β
b
θ

!

(
=

ZT Z

!

ZT BQ2

+λ

T
T T
QT
2 B Z Q2 B BQ2

!)−1

0
QT
2 PQ2

ZT Y
T
QT
2B Y

!
.

Letting
V=

V11 V12

!
=

V21 V22

ZT Z

!

ZT BQ2

T
T
T
QT
2 B Z Q2 (B B + λP)Q2

,

(2.7)

we have
b
β
b
θ

!
= V−1

ZT Y

!

T
QT
2B Y

.

Next let us write
V

−1

≡U=

U11 U12
U21 U22

!
=

U11

−1
−U11 V12 V22

−1
−U22 V21 V11

U22

!
,

(2.8)

where


−1
T
−1 T T
U−1
= V11 − V12 V22
V21 = ZT I − BQ2 {QT
Z, (2.9)
2 (B B + λP)Q2 } Q2 B
11



−1
T
U−1
= V22 − V21 V11
V12 = QT
I − Z(ZT Z)−1 ZT B + λP Q2 .
(2.10)
2 B
22

8

Li Wang, Guannan Wang, Min-Jun Lai and Lei Gao

Then the minimizers of (2.7) can be given precisely as follows:


b = U11 ZT I − BQ2 V−1 QT BT Y=U11 ZT I − BQ2 {QT (BT B + λP)Q2 }−1 QT BT Y,
β
2
2
2
22


−1 T
T T
T T
T
−1 T
b
θ = U22 Q2 B I − ZV11 Z Y=U22 Q2 B I − Z(Z Z) Z Y.
Therefore, one obtains the estimators for γ and g(·), respectively:

b = Q2 U22 QT BT I − Z(ZT Z)−1 ZT Y,
b = Q2 θ
γ
2
b=
gb(x) = B(x)T γ

X

Bξ (x)b
γξ .

(2.11)

ξ∈K

The fitted values at the n data points are
b + Bb
b = Zβ
Y
γ = S(λ)Y,
where the smoothing or hat matrix is


−1 T
−1 T T
T
Z .
Q2 B + BQ2 U22 QT
I − ZV11
S(λ) = ZU11 ZT I − BQ2 V22
2B
In nonparametric regression, the trace tr(S(λ)) of smoothing matrix S(λ) is often called the
degrees of freedom of the model fit. It has the rough interpretation as the equivalent number
of parameters and can be thought as a generalization of the definition in linear regression.
Finally, we can estimate the variance of the error term, σ 2 by
σ̂ 2 =

b 2
kY − Yk
.
n − tr(S(λ))

(2.12)

2.4. Choosing the Triangulation
Triangulation has been extensively investigated in the past few decades, and various packages have been developed. For example, one can use the “Delaunay” algorithm to find a
triangulation; see MATLAB program delaunay.m or MATHEMATICA function DelaunayTriangulation. “Triangle” (Shewchuk, 1996) is also widely used in many applications, and one can
download it for free from http://www.cs.cmu.edu/ quake/triangle.html. It is a C++ program for
two-dimensional mesh generation and construction of Delaunay triangulations. “DistMesh” is
another method to generate unstructured triangular and tetrahedral meshes; see the DistMesh
generator on http://persson.berkeley.edu/distmesh/. A detailed description of the program is
provided by Persson and Strang (2004). We used our own triangulation code in simulation
studies and real data analysis below.
As is usual with the one-dimensional (1-D) penalized least squares (PLS) splines, the
number of knots is not important given that it is above some minimum depending upon the

Efficient Estimation of Partially Linear Models for Spatial Data over Complex Domains

9

degree of the smoothness; see Li and Ruppert (2008). For bivariate PLS splines, Lai and Wang
(2013) also observed that the number of triangles N is not very critical, provided N is larger
than some threshold. In fact, one of the main advantages of using PLS splines over discrete
least squares (DLS) splines is the flexibility of choosing knots in the 1-D setting and choosing
triangles in the 2-D setting. For DLS splines, one has to have large enough sample according to
the requirement of the degree of splines on each subinterval in the 1-D case or each triangle in
the 2-D case to guarantee that a solution can be found. However, there is no such requirement
for PLS splines. When the smoothness r ≥ 1, the only requirement for bivariate PLS splines is
that there is at least one triangle containing three points which are not in one line (Lai, 2008).
Also, PLS splines perform similar to DLS splines as long as the penalty parameter λ is very
small. So in summary, the proposed bivariate PLS splines are very flexible and convenient for
data fitting, even for smoothing sparse and unevenly sampled data.
In practice, we recommend that the user first constructs a polygon domain Ω containing all
the design points of the data and makes a simple triangulation 40 of Ω by hand or computer,
then refines 40 several times to have a triangulation of desired size.
2.5. Penalty Parameter Selection
Selecting a suitable value of smoothing parameter λ is critical to good model fitting. A
large value of λ enforces a smoother fitted function with potentially larger fitting errors, while
a small value yields a rougher fitted function and potentially smaller fitting errors. Since the
in-sample fitting errors can not gauge the prediction property of the fitted function, one should
target a criterion function that mimics the out-of-sample performance of the fitted model. The
generalized cross-validation (GCV) (Craven and Wahba, 1979; Wahba, 1990) is such a criterion
and is widely used for choosing the penalty parameter. We choose the smoothing parameter λ
by minimizing the following generalized cross-validation (GCV) criterion

GCV(λ) =

nkY − S(λ)Yk2
,
{n − tr(S(λ))}2

over a grid of values of λ. We use the 10-point grid where the values of log10 (λ) are equally
spaced between −6 and 7.
3. IAsymptotic Results
This section studies the asymptotic properties for the proposed estimators. To discuss
these properties, we first introduce some notation. For any function f over the closure of
domain Ω, denote kf k∞ = supx∈Ω |f (x)| the supremum norm of function f and |f |υ,∞ =
maxi+j=υ kDxi 1 Dxj 2 f (x)k∞ the maximum norms of all the υth order derivatives of f over Ω.

10

Li Wang, Guannan Wang, Min-Jun Lai and Lei Gao

Let
W `,∞ (Ω) = {f on Ω : |f |k,∞ < ∞, 0 ≤ k ≤ `}

(3.1)

be the standard Sobolev space. For any j = 1, . . . , p, let zj be the coordinate mapping that
maps z to its j-th component so that zj (Zi ) = Zij , and let
hj = argminh∈L2 kzj − hk2L2 = argminh∈L2 E{(Zij − h(Xi ))2 }

(3.2)

be the orthogonal projection of zj onto L2 .
Before we state the results, we make the following assumptions:
(A1) The random variables Zij are bounded, uniformly in i = 1, . . . , n, j = 1, . . . , p.


T 

T
T
(A2) The eigenvalues of E
Xi are bounded away from 0.
1 Zi
1 Zi


(A3) The noise  satisfies that limη→∞ E 2 I( > η) = 0.
Assumptions (A1)–(A3) are typical in semiparametric smoothing literature, see for instance
Huang, Zhang and Zhou (2007) and Wang, et al (2011). The purpose of Assumption (A2) is
to ensure that the vector (1, ZT
i ) is not multicolinear. We next introduce some assumptions
on the properties of the true bivariate function in model (1.1) and the bivariate spline space
introduced in Section 2.2.
(C1) The bivariate functions hj (·), j = 1, . . . , p, and the true function in model (1.1) g(·) ∈
W `+1,∞ (Ω) in (3.1) for an integer ` ≥ 2.
(C2) For every s ∈ S and every τ ∈ 4, there exists a positive constant F1 , independent of s
and τ , such that
1/2


X

F1 ksk∞,τ ≤ 

s (Xi )2 

,

for all τ ∈ 4.

(3.3)

Xi ∈τ, i=1,··· ,n

(C3) Let F2 be the largest among the numbers of observations in triangles τ ∈ 4 in the sense
that
1/2


X


s (Xi )2 

≤ F2 ksk∞,τ ,

for all τ ∈ 4,

(3.4)

Xi ∈τ, i=1,··· ,n

where ksk∞,τ denotes the supremum norm of s over triangle τ . The constants F1 and F2
defined in (3.3) and (3.4) satisfy F2 /F1 = O(1).

Efficient Estimation of Partially Linear Models for Spatial Data over Complex Domains

11

(C4) The number N of the triangles and the sample size n satisfy that N = Cnγ for some
constant C > 0 and 1/(` + 1) ≤ γ ≤ 1/3.
(C5) The penalized parameter λ satisfies λ = o(n1/2 N −1 ).
Condition (C1) describes the requirement for the true bivariate function as usually used
in the literature of nonparametric or semiparametric estimation. Condition (C2) ensures the
existence of a least squares spline. Although one can get a decent penalized least squares spline
fitting with F1 = 0 for some triangles, we need (C2) to study the convergence of bivariate
penalized least squares splines. Condition (C3) suggests that we should not put too many
observations in one triangle. Condition (C4) requires that the number of triangles is above
some minimum depending upon the degree of the spline, which is similar to the requirement
of Li and Ruppert (2008) in the univariate case. It also ensures the asymptotic equivalence
of the theoretical and empirical inner products/norms defined at the beginning of Section .
Condition (C5) is required to reduce the bias of the bivariate spline approximation through
“under smoothing” and “choosing smaller λ”.
To avoid confusion, in the following let β 0 and g0 be the true parameter value and function
b is root-n and β
b is
in model (1.1). The following theorem states that the rate convergence of β
asymptotically normal.
b is asympTheorem 1. Suppose Assumptions (A1)-(A3), (C1)-(C5) hold, then the estimator β
totically normal, that is,
b − β ) → N (0, I),
(nΣ)1/2 (β
0
where I is a p × p identity matrix,
e i )(Zi − Z
e i )T }
Σ = σ −2 E{(Zi − Z

(3.5)

e i = {h1 (Xi ), . . . , hp (Xi )}T , for hj (·) defined in (3.2), j = 1, . . . , p. In addition, Σ can
with Z
be consistently estimated by
Σn =

n
1 X
b i )(Zi − Z
b i )T = 1 (Z − Z)
b T (Z − Z).
b
(Zi − Z
2
nb
σ
nb
σ2

(3.6)

i=1

b i is the i-th column of Z
b T = ZT BQ2 V−1 QT BT and σ
where Z
b2 is given by (2.12).
2
22
The results in Theorem 1 enable us to construct confidence intervals for the parameters.
The next theorem provides the global convergence of the nonparametric estimator gb(·).
Theorem 2. Suppose Assumptions (A1)-(A3), (C1)-(C5) hold, then the bivariate penalized
estimator gb(·) in ( 2.11) is consistent with the true function g0 , and satisfies that




λ
λ
F2
1
`+1
√
kb
g − g0 kL2 = OP
|g
|
+
1
+
|4|
|g
|
+
.
0 2,∞
0 `,∞
n|4|
n |4|3
n |4|5 F1

12

Li Wang, Guannan Wang, Min-Jun Lai and Lei Gao

The proofs of the above two theorems are given in Appendix. We notice that the rate
of convergence given in Theorem 2 is the same as those for nonparametric spline regression
without including the covariate information obtained in Lai and Wang (2013).
4. Simulation
In this section, we carry out two numerical studies to assess the performance of our proposed method. We compare the performance of bivariate penalized splines over triangulations
(BPST) with filtered kriging (KRIG), thin plate splines (TPS), soap film smoothing (SOAP)
in Wood, Bravington and Hedley (2008), linear finite elements (LFE) and quadratic finite
elements (QFE) in Sangalli, Ramsay and Ramsay (2013).
Example 1. In this example, we consider a modified horseshoe domain with the surface test
function g(·) used by Wood, Bravington and Hedley (2008) and Sangalli, Ramsay and Ramsay
(2013). In particular, for 201 × 501 grid points over the domain, we simulate data as follows,
the response variable Y is generated from the following PLM:
Y = β1 Z1 + β2 Z2 + g(X1 , X2 ) + .
The random error, , is generated from an N (0, σ2 ) distribution with σ = 0.5. In addition,
we set the parameters as β1 = −1, β2 = 1. For the design of the explanatory variables, Z1
and Z2 , two scenarios are considered based on the relationship between the location variables
(X1 , X2 ) and covariates (Z1 , Z2 ). Under both scenarios, Z1 ∼ unif orm[−1, 1]. On the other
hand, the variable Z2 = cos[π(ρ(X12 + X22 ) + (1 − ρ)U )] where U ∼ unif orm[−1, 1] and is
independent from (X1 , X2 ) as well as Z1 . We consider both independent design: ρ = 0.0 and
dependent design: ρ = 0.7 in this example. Under both scenarios, 100 Monte Carlo replicates
are generated. For each replication, we randomly sample n = 200 locations uniformly from the
grid points inside the horseshoes domain.
Figure 1 (a) and (b) show the surface and the contour map of the true function g(·),
respectively. Figure 1 (c) demonstrates the sampled location points of replicate 1, and Figure1
(d) and (e) illustrates two different triangulations used in the BPST method. In the first
triangulation (41 ), we use 91 triangles and 74 vertices, while in the second one (42 ) we use
109 triangles and 95 vertices.
[Figure 1 about here]
To implement the TPS and KRIG methods, we use the R package fields under the standard
implementation setting of Furrer, Nychka and Sainand (2011). For KRIG, we try different
covariance structures, and we choose the Matérn covariance with smoothness parameter ν = 1,
which gives the best prediction. For the SOAP method, we implement it by using R package

Efficient Estimation of Partially Linear Models for Spatial Data over Complex Domains

13

mgcv (Wood, Bravington and Hedley, 2008) with 32 interior knots. In addition, a rank 39
(40-knot) cyclic penalized cubic regression spline is used as the boundary curve. For all the
methods requiring a smoothing or roughness parameter, GCV is used to choose the values of
the parameter.
To see the accuracy of the estimators, we compute the root mean squared error (RMSE)
for each of the components based on 100 Monte carlo samples. Table 1 shows the RMSEs of
the estimate of the parameters β1 , β2 , σε as well as the average mean squared prediction error
(MSPE) of the nonlinear function g(·) over all the grid points in the horseshoes shaped domain.
From Table 1, one sees that methods SOAP and BPST give very comparable estimates of the
parameters β1 and β2 , respectively, while BPST method produces the best prediction of the
nonlinear function g(·), regardless of the choice of triangulation.
[Table 1 about here]
Figure 2 shows the estimated functions via different methods for the first replicate. For
the test function on Figure 2, the BPST estimate looks visually better than the other four
estimates. In addition, one sees there is a “leakage effect” in KRIG and TPS estimates and
the SOAP, LFE, QFE and BPST methods improve the model fitting of KRIG and TPS. The
poor performance of KRIG and TPS is because they do not take the complex boundary into
any account and smooth across the gap inappropriately. In addition, from Figure 2, one also
sees that the BPST estimators based on 41 and 42 are very similar, which agrees our findings
for penalized splines that the number of triangles is not very critical for the fitting as long as
it is large enough.
[Figure 2 about here]
Next we test the accuracy of the standard error formula in (3.6) for β̂1 and β̂2 , and the
results are listed in Table 2. The standard deviations of the estimated parameters are computed based on 100 replications, which can be regarded as the true standard errors (column
labeled “SEmc ”) and compared with the mean and median of the 100 estimated standard errors
calculated using (3.6) (columns labeled “SEmean ” and “SEmedian ”, respectively). The column
labeled “SEmad ” is the interquartile range of the 100 estimated standard errors divided by
1.349, which is a robust estimate of the standard deviation. From Table 2 one observes that
the averages or medians of the standard errors calculated using the formula are very close
to the true standard deviations, which confirms the accuracy of the proposed standard error
formula.
[Table 2 about here]

14

Li Wang, Guannan Wang, Min-Jun Lai and Lei Gao

Finally, for computing time, it takes about 60 seconds to compute BPST in one run of
simulation with even 2000 observations over 319 triangles under the computing environment
of x64 PC with Intel Dual Core i5. Overall, the proposed algorithm is fast to compute.
Example 2. In this example, we consider a rectangular domain, [0, 1]2 , where there is no
irregular shape or complex boundaries problem. In this case, classical methods for spatial data
analysis, such as KRIG and TPS, will not encounter any difficulty. We obtain the true signal
and noisy observation for each coordinate pair lying on a 101 × 101-grid over [0, 1]2 using the
following model:
Y = ZT β + g(X1 , X2 ) + ξ(X1 , X2 ) + ,
where β = (−1, 1)T and g(x1 , x2 ) = 5 sin(2π(x21 + x22 )). The random error, , is generated from
an N (0, σ2 ) distribution with σ = 0.5, and the process ξ(·) is generated from a stationary
gaussian random field with the Matérn(1,1) covariance structure. The components of Z and
 are standard normal, and Z, ξ and  are independent of each other. Similar to Example
1, we simulate Z1 ∼ unif orm[−1, 1], and Z2 = cos[π(ρ(X12 + X22 ) + (1 − ρ)U )], where ρ =
0.0 or 0.7, U ∼ unif orm[−1, 1] and is independent from (X1 , X2 ) and Z1 . Next we take 100
Monte Carlo random samples of size n = 200 from the 101 × 101 points.
Figure 4 (a) and (b) display the true sinusoidal surface and the contour map, respectively.
4 (c) represents the Matérn structure. We use the triangulation in Figure 4 (e) and (f), and
there are 8 triangles and 9 vertices as well as 32 triangles and 25 vertices, respectively. In
addition, the points in Figure 4 (d) demonstrate the sampled location points of replicate 1.
[Figure 4 about here]
Similar to the study in Example 4.1, we also compare the proposed BPST estimator with
estimators from the KRIG, TPS, SOAP, LFE and QFE methods, which are implemented in
the same way as Example 4.1. Specifically, for KRIG, we choose the true Matérn covariance
with smoothness parameter ν = 1. To see the accuracy of the estimators, we compute the
RMSEs of the coefficient estimators and the estimator of σ . To see the overall prediction
accuracy, we make prediction over the 101 × 101 grid points on the domain for each replication
using different methods, and compare the predicted values with the true observations of Y at
these grid points, and we report the average mean squared prediction errors (MSPE) over all
replications.
All the results are summarized in Table 3. As expected, KRIG works pretty well since the
true covariance structure is used in the KRIG fitting. When ρ = 0, TPS, KRIG and BPST
all perform very well. When ρ = 0.7, KRIG and BPST are the best among all the estimators,
and both of them outperform the rest of the estimators. In both scenarios, the differences
between BPST and KRIG are almost unnoticeable. One also notices that, compared with the

Efficient Estimation of Partially Linear Models for Spatial Data over Complex Domains

15

FEM estimators (LEM and QEM), our BPST estimator shows much better performance in
terms of both estimation and prediction, because BPST provides a more flexible and easier
construction of splines with piecewise polynomials of various degrees and smoothness than the
FEM method. Finally as pointed out in Wood, Bravington and Hedley (2008), the standard
(linear) FEM method requires a very fine triangulation in order to reach certain approximation
power, however, BPST doesn’t need such a strict fineness requirement as it uses piecewise
polynomials of degree ≥ 5 yielding an higher order approximation power. In fact, we have
used 64 times more triangles in the FEM than that for the BPST in our simulation. That is,
the BPST is computationally more efficient than the FEM to approximate smooth functions.
[Table 3 about here]
Table 4 lists the accuracy results of the standard error formula in (3.6) for β̂1 and β̂2 using
BPST. From Table 4, one sees that the estimated standard errors based on sample size n = 200
are very accurate.
[Table 4 about here]

5. Application to California House Value Data
In this section we apply the proposed method to analyze the California house value data
from the StatLib repository. The data appeared in Pace and Barry (1997). The spatial data
consists of information of all the block groups in California defined by centroid of census
enumeration areas. In the data set, a block group on average includes 1425.5 individuals living
in a geographically compact area and there are 20, 640 blocks in the final data.
In this paper, we study how different features and factors influence the real estate property
prices. The response variable is the median house value (Value). The investigated factors
include the median income (MedInc), median house age (Age), the average number of bedrooms
(AveBedrms), housing density as reflected by the number of households (Hhd) in each block,
and the average occupancy in each household (AveOccup). It is obvious that the location
of a house is very crucial for making an accurate prediction, so we also include the latitude
(Latitude) and longitude (Longitude) of the block. We model the median house value using
the following PLM:
log(Value) = β0 + β1 MedInc + β2 log(Age) + β3 log(AveBedrms)
+β4 log(AveOccup) + β5 log(Hhd) + g(Latitude, Longitude).

(5.1)

To fit model (5.1), we use six different methods: KRIG, TPS, SOAP, LFE, QFE, and BPST.

16

Li Wang, Guannan Wang, Min-Jun Lai and Lei Gao

For comparison, we also consider the purely linear model without using the spatial information:
log(Value) = β0 + β1 MedInc + β2 log(Age) + β3 log(AveBedrms)
+β4 log(AveOccup) + β5 log(Hhd),

(5.2)

and fit it using the ordinary linear least squares (OLS) method.
In the following, we report our estimation results for the linear and nonlinear components
in model (5.1) by using the OLS, KRIG, TPS, SOAP, LFE, and QFE methods. Table 6
summarizes the coefficient estimation results of all procedures, where the BPST confidence
intervals are constructed based on the standard error (s.e.) calculated using (3.6). Note that
the coefficient on “log(Age)” is positive 0.165 (s.e. = 0.006) in the fitted linear regression
model (5.2) via OLS, which is obviously against the common sense in household real estate
market. Using the PLM (5.1) with the spatial coordinates, we obtain negative coefficients
for “log(Age)” regardless which spatial based method is employed. The coefficient for the
variable “log(AveBedrms)” is −0.037 with a standard error 0.018 using OLS. In contrast, all
the semiparametric methods unanimously suggest that the coefficient for “log(AveBedrms)” is
positive in the PLM in (5.1) after the location effect is controlled. In addition, the coefficient on
“log(Hhd)” is positive 0.088 and very significant (s.e. = 0.004) in the OLS, but it is statistically
insignificant when we apply the PLM (5.1) to the data. In summary, compared with the
OLS method, our results from PLM are more consistent with intuitions in the real estate
market study. The above counter-intuitive phenomenon in OLS is perhaps due to the model
misspecification of (5.2) that completely ignores the location factor of the property.
[Table 6 about here]
Based on the median house values, we classify the houses in the dataset into six different
groups: (1) less than 50K, (2) 50K–100K, (3) 100K–150K, (4)150K–200K, (5) 200K–300K, and
(6) greater than 300K, and these groups are plotted in Figure 5. We plot the estimated median
house values using different methods in Figure 6, where different colors are used to represent
the value of houses as in Figure 5. All the plots in Figures 5 and 6 show that expensive houses
are clustered around the major cities and inland house values are lower than coastal house
values. The OLS significantly underestimates the coastal enclaves of expensive houses and
overestimates the house values in the central valley; see Figure 6 (a). In contrast, the PLM
methods in Figure 6 (b)-(g) provide much more accurate estimates.
[Figures 5, 6 about here]
Figure 7 further demonstrates the differences between the estimated house values and
the observed house values using seven different methods. As seen in Figure 7, some methods,

Efficient Estimation of Partially Linear Models for Spatial Data over Complex Domains

17

especially the OLS, have difficulties in estimating the house values in major metropolitan areas
and the central valley. The proposed BPST method dramatically reduces the estimation errors
in those areas.
[Figure 7 about here]
To evaluate different methods, we also estimate the out-of-sample prediction errors of each
method using the 5-fold cross validation. We randomly split all the observations into five
roughly equal-sized parts. For each k = 1, . . . , 5, we leave out part k, fit the model to the
other four parts (combined), and then obtain predictions for the left-out kth part. Table 5
summarizes the mean squared prediction errors of the logarithm of the median house value
based on different methods. As expected, incorporating the spatial information dramatically
reduces the prediction errors, and the proposed BPST method provides the most accurate
prediction on the median housing values among all the methods.
[Table 5 about here]
6. Concluding Remarks
In this paper, we have considered PLMs for analyzing spatial data. We introduce a framework of bivariate penalized splines over triangulations in the semiparametric estimation. Our
work differs from existing works in four major aspects.
First, the proposed estimator solves the problem of “leakage” across the complex domains
where many conventional smoothing tools suffer from. The numerical results of the simulation
show our method is very effective on both regular and irregular domains.
Secondly, we provide new statistical theories for estimating the PLM for data distributed
over complex spatial domains. It is shown that our estimates of both parametric part and
non-parametric part of the model enjoy excellent asymptotic properties. In particular, we
have shown that our estimates of the coefficients in the parametric part are asymptotically
normal and derived the convergence rate of the nonparametric component under regularity
conditions. We have also provided a standard error formula for the estimated parameters and
our simulation studies show that the standard errors are estimated with good accuracy. The
theoretical results provide measures of the effect of covariates after adjusting for the location
effect. In addition, they give valuable insights into the accuracy of our estimate of the PLM
and permit joint inference for the parameters.
Thirdly, comparing with the conventional FEM, our method provides a more flexible and
easier construction of splines with piecewise polynomials of various degrees and various smoothness.

18

Li Wang, Guannan Wang, Min-Jun Lai and Lei Gao

Finally, the proposed method greatly enhances the application of PLMs to spatial data
analysis. We don’t require the data to be evenly distributed or on regular-spaced grid. When
we have regions of sparse data, PLS splines provides a more convenient tool for data fitting
than DLS splines since the roughness penalty helps regularize the estimation. Our estimation
is computationally fast and efficient since it can formulate a penalized regression problem using
a QR decomposition.
This paper leaves open the problem of how to choose a triangulation for a given data set.
An optimal triangulation 4 enables us to achieve the best approximation of the given data
using a spline space S(4) of a fixed degree and fixed smoothness. We leave the problem for
our future exploration.

Acknowledgment
The first author’s research was supported in part by National Science Foundation grants
DMS-11-06816 and DMS-15-42332, and the third author’s research was supported in part by
National Science Foundation grant DMS-15-21537.
Appendices
A.1. Preliminaries
In the following, we use c, C, c1 , c2 , C1 , C2 , etc. as generic constants, which may be
different even in the same line. For functions f1 and f2 on Ω × Rp , we define the empirical
P
inner product and norm as hf1 , f2 in = n1 ni=1 f1 (Xi , Zi )f2 (Xi , Zi ) and kf1 k2n = hf1 , f1 in . If
f1 and f2 are L2 -integrable, we define the theoretical inner product and theoretical L2 norm
as hf1 , f2 iL2 = E {f1 (Xi , Zi )f2 (Xi , Zi )} and kf1 kL2 = hf1 , f1 iL2 . Furthermore, let k·kEυ be the
norm introduced by the inner product h·, ·iEυ , where
1/2 
1/2

Z  X  
  X υ 

υ
hg1 , g2 iEυ =
(Dxi 1 Dxj 2 g1 )2
(Dxi 1 Dxj 2 g2 )2
dx1 dx2
 

i
i
Ω
i+j=υ

i+j=υ

for g1 and g2 on Ω.
Lemma A.1. [Lai and Schumaker (2007)] Let {Bξ }ξ∈K be the Bernstein polynomial basis for
spline space S with smoothness r, where K stands for an index set. Then there exist positive
constants c, C depending on the smoothness r and the shape parameter δ such that
X
X
X
c|4|2
|γξ |2 ≤ k
γξ Bξ k2L2 ≤ C|4|2
|γξ |2
ξ∈K

ξ∈K

ξ∈K

for all γξ , ξ ∈ K.
With the above stability condition, Lai and Wang (2013) established the following uniform
rate at which the empirical inner product approximates the theoretical inner product.

19

Efficient Estimation of Partially Linear Models for Spatial Data over Complex Domains

Lemma A.2. [Lemma 2 of the Supplement of Lai and Wang (2013)] Let g1 =
P
g2 = ζ∈K e
cζ Bζ be any spline functions in S. Under Condition (C4), we have
sup
g1 ,g2 ∈S

P

ξ∈K cξ Bξ ,

n
o
hg1 , g2 in − hg1 , g2 iL2
= OP (N log n)1/2 /n1/2 .
kg1 kL2 kg2 kL2

Lemma A.3. [Corollary of Theorem 6 in Lai (2008)] Assume g(·) is in Sobolev space W `+1,∞ (Ω).
For bi-integer (α1 , α2 ) with 0 ≤ α1 + α2 ≤ υ, there exists a unique spline fit g ∗ (·) ∈ S such that
kDxα11 Dxα22 (g − g ∗ ) k∞ ≤ C

F2
|4|`+1−α1 −α2 |g|`+1,∞ ,
F1

where C is an absolute constant depending on r and δ.
For any smooth bivariate function g(·) and λ > 0, define
sλ,g = argmins∈S

n
X

{g(Xi ) − s(Xi )}2 + λEυ (s)

(A.1)

i=1

the penalized least squares splines of g(·). Then s0,g is the nonpenalized spline estimator of
g(·).
Lemma A.4. Assume g(·) is in Sobolev space W `+1,∞ (Ω). For bi-integer (α1 , α2 ) with 0 ≤
α1 +α2 ≤ υ, there exists an absolute constant C depending on r and δ, such that with probability
approaching 1,
kDxα11 Dxα22 (g − s0,g ) k∞ ≤ C

F2
|4|`+1−α1 −α2 |g|`+1,∞ .
F1

Proof. Note that
kDxα11 Dxα22 (g − s0,g ) k∞ ≤ kDxα11 Dxα22 (g − g ∗ ) k∞ + kDxα11 Dxα22 (g ∗ − s0,g ) k∞
≤ kDxα11 Dxα22 (g − g ∗ ) k∞ + kDxα11 Dxα22 (s0,g∗ −g ) k∞ .
The desired result follows from Lemma A.3, the projection bound result in Golitschek and
Schumaker (2002), and the differentiation properties of bivariate splines over triangulatins
given in Lai and Schumaker (2007). 
Lemma A.5. Suppose g(·) is in the Sobolev space W `+1,∞ (Ω), and let sλ,g be its penalized
spline estimator defined in (A.1). Under Conditions (C2) and (C3),

kg − sλ,g kn = OP

F2
λ
|4|`+1 |g|`+1,∞ +
F1
n |4|2



F2
|g|υ,∞ +
|4|`+1−υ |g|`+1,∞
.
F1

20

Li Wang, Guannan Wang, Min-Jun Lai and Lei Gao

Proof. Note that sλ,g is characterized by the orthogonality relations
n hg − sλ,g , uin = λ hsλ,g , uiEυ ,

for all u ∈ S,

(A.2)

while s0,g is characterized by
hg − s0,g , uin = 0,

for all u ∈ S.

(A.3)

By (A.2) and (A.3), n hs0,g − sλ,g , uin = λ hsλ,g , uiEυ , for all u ∈ S. Replacing u by s0,g − sλ,g
yields that
n ks0,g − sλ,g k2n = λ hsλ,g , s0,g − sλ,g iEυ .

(A.4)

Thus, by Cauchy-Schwarz inequality,
n ks0,g −

sλ,g k2n


≤ λ ksλ,g kEυ ks0,g − sλ,g kEυ ≤ λ ksλ,g kEυ sup
f ∈S


kf kEυ
, kf kn 6= 0 ks0,g − sλ,g kn .
kf kn

Similarly, using (A.4), we have
n
o
n ks0,g − sλ,g k2n = λ hsλ,g , s0,g iEυ − hsλ,g , sλ,g iEυ ≥ 0.
Therefore, by Cauchy-Schwarz inequality,
ksλ,g k2Eυ ≤ hsλ,g , s0,g iEυ ≤ ksλ,g kEυ ks0,g kEυ ,
which implies that ksλ,g kEυ ≤ ks0,g kEυ . Therefore,
ks0,g − sλ,g kn ≤ n

−1


λ ks0,g kEυ sup
f ∈S

By Lemma A.4, with probability approaching 1,
(
X
ks0,g kEυ ≤ C1 AΩ |g|υ,∞ +


kf kEυ
, kf kn 6= 0 .
kf kn

)
Dxα11 Dxα22

(g − s0,g )

∞

α1 +α2 =υ


≤ C 2 AΩ


F2
`+1−υ
|g|υ,∞ +
|4|
|g|`+1,∞ ,
F1

(A.5)

where AΩ denotes the area of Ω. By Markov’s inequality, for any f ∈ S, kf kEυ ≤ C |4|−2 kf kL2 .

Lemma (A.2) implies that supf ∈S { kf kn / kf kL2 } ≥ 1−OP (N log n)1/2 /n1/2 . Thus, we have

sup
f ∈S


h
n
oi−1/2


kf kEυ
, kf kn 6= 0 ≤ C |4|−2 1 − OP (N log n)1/2 /n1/2
= OP |4|−2 . (A.6)
kf kn

21

Efficient Estimation of Partially Linear Models for Spatial Data over Complex Domains

Therefore,

ks0,g − sλ,g kn = OP

λ
n |4|2



F2
`+1−υ
|g|υ,∞ +
|4|
|g|`+1,∞
,
F1

and
kg − sλ,g kn ≤ kg − s0,g kn + ks0,g − sλ,g kn .
By Lemma A.4,

kg − s0,g kn ≤ kg − s0,g k∞ = OP


F2
`+1
|4|
|g|`+1,∞ .
F1

Thus, the desired result is established. 
Lemma A.6. Under Assumptions (A1), (A2), (C4), (C5), there exist constants 0 < cU <
CU < ∞, such that with probability approaching 1 as n → ∞, cU Ip×p ≤ nU11 ≤ CU Ip×p , where
U11 is given in (2.8).
Proof. Denote by
" n
#
X

1
1
λ
Γλ =
BT B + λP =
Bξ (Xi ) Bζ (Xi ) + hBξ , Bζ iEυ
n
n
n
i=1

ξ,ζ∈K

a symmetric positive definite matrix. Then for V22 defined in (2.7), we can rewrite it as
V22 = nQT
2 Γλ Q2 . Let αmin (λ) and αmax (λ) be the smallest and largest eigenvalues of Γλ .
As shown in the proof of Theorem 2 in the Supplement of Lai and Wang (2013), there exist
positive constants 0 < c3 < C3 such that under Conditions (C4) and (C5), with probability
approaching 1, we have

c3 |4| ≤ αmin (λ) ≤ αmax (λ) ≤ C3 |4|2 +
2

λ
n|4|2


.

Therefore, we have

c4 |4|2 +

λ
n|4|2

−1

−1
−1
−2
2
kak2 ≤ naT V22
a = aT (QT
2 Γλ Q2 ) a ≤ C4 |4| kak .

Thus, by Assumption (A2), we have with probability approaching 1

c5 |4|2 +

λ
n|4|2

−1

−1
−1 T T
|4|2 kak2 ≤ aT V12 V22
V21 a = aT ZT BQ2 V22
Q2 B Za ≤ C5 kak2 .

(A.7)
According to (2.8) and (2.9), we have
−1
−1
(nU11 )−1 = n−1 (V11 − V12 V22
V21 ) = n−1 (ZT Z − V12 V22
V21 ).

22

Li Wang, Guannan Wang, Min-Jun Lai and Lei Gao

2

The desired result follows from Assumptions (A1), (A2) and (A.7).
A.2. Proof of Theorem 1
T
T
Let µi = ZT
i β 0 + g0 (Xi ), µ = (µ1 , . . . , µn ), and let  = (1 , . . . , n ). Define


e = U11 ZT I − BQ2 V−1 QT BT µ,
β
µ
2
22

e = U11 ZT I − BQ2 V−1 QT BT .
β

2
22

(A.8)
(A.9)

b − β = (β
e −β )+β
e .
Then β
0
µ
0


e − β k = oP n−1/2 for β
e in
Lemma A.7. Under Assumptions (A1), (A2), (C1)-(C5), kβ
µ
0
µ
(A.8).
Proof. Let g0 = (g0 (X1 ), . . . , g0 (Xn ))T . It is clear that

e − β = U11 ZT I − BQ2 V−1 QT BT g0
β
µ
0
2
22


T
−1 T T
= U11 ZT g0 − BQ2 {QT
2 (B B + λP)Q2 } Q2 B g0
= nU11 A,
where A = (A1 , . . . , Ap )T , with


T
T
−1 T T
Aj = n−1 ZT
j g0 − BQ2 {Q2 (B B + λP)Q2 } Q2 B g0
for ZT
j = (Z1j , ..., Znj ). Next we derive the order of Aj , 1 ≤ j ≤ p, as follows. For any gj ∈ S,
by (A.2) we have
Aj = hzj , g0 − sλ,g0 in = hzj − gj , g0 − sλ,g0 in +

λ
hsλ,g0 , gj iEυ .
n

For any j = 1, . . . , p, let hj (·) be the function h(·) that minimizes E{Zij − h(Xi )}2 as defined
in (3.2). According to Lemma A.3, there exists a function e
hj ∈ S satisfy
F2
ke
hj − hj k∞ ≤ C
|4|`+1 |hj |`+1,∞ ,
F1

(A.10)

then
Aj = hzj − hj , g0 − sλ,g0 in + hhj − e
hj , g0 − sλ,g0 in +

λ
hsλ,g0 , e
hj iEυ = Aj,1 + Aj,2 + Aj,3 .
n

Since hj satisfies hzj − hj , ψiL2 (Ω) = 0 for any ψ ∈ L2 (Ω), E (Aj,1 ) = 0. According to Proposition 1 in Lai and Wang (2013),



F2
λ
F2
`+1
`−1
kg0 − sλ,g0 k∞ = OP
|4|
|g0 |`+1,∞ +
|g0 |2,∞ +
|4|
|g0 |`+1,∞
.
F1
F1
n |4|3

Efficient Estimation of Partially Linear Models for Spatial Data over Complex Domains

23

Next,
Var (Aj,1 ) =

n
kg0 − sλ,g0 k2∞
1 X
2
kzj − hj k2L2 ,
E
[{Z
−
h
(X
)}
(g
−
s
)]
≤
ij
j
i
0
λ,g0
n2
n
i=1

which together with E (Aj,1 ) = 0 implies that



F2
λ
F2
`+1
`−1
|Aj,1 | = OP
|4|
|g
|
+
|g
|
+
|4|
|g
|
.
0 `+1,∞
0 2,∞
0 `+1,∞
F1
n1/2 F1
n3/2 |4|3
Cauchy-Schwartz inequality, Lemma A.5 and (A.10) imply that
|Aj,2 | ≤ khj − e
hj kn kg0 − sλ,g0 kn


F2
`+1
|4|
|hj |`+1,∞
= OP
F1



λ
F2
F2
`+1
`−1
|4|
|g0 |`+1,∞ +
|4|
|g0 |`+1,∞
.
|g0 |2,∞ +
× OP
F1
F1
n |4|2
Finally, by (A.5), we have
|Aj,3 | ≤
≤

λ
λ
ksλ,g0 kEυ ke
hj kEυ ≤ ks0,g0 kEυ ke
hj kEυ
n 
n


λ
F2
F2
`−1
`−1
C1 |g0 |2,∞ +
|4|
|g0 |`+1,∞
|hj |2,∞ +
|4|
|hj |`+1,∞ .
n
F1
F1

Combining all the above results yields that




λ
F2
F2
`−1
|4|`+1 |g0 |`+1,∞ +
|g
|
+
|4|
|g
|
|Aj | = OP n−1/2
0 2,∞
0 `+1,∞
F1
F1
n |4|3
for j = 1, . . . , p. By Assumptions (C3)-(C5), |Aj | = oP (n−1/2 ), for j = 1, . . . , p. In addition,

e − β k = oP n−1/2 .
we have nU11 = OP (1) according to Lemma A.6. Therefore, kβ
2
µ
0
Lemma A.8. Under Assumptions (A1)-(A3) and (C1)-(C5), as n → ∞,
h


i−1/2
e −→ N (0, Ip×p ) ,
e |{(Xi , Zi ) , i = 1, . . . , n}
Var β
β



e is given in (A.9).
where β

Proof. Note that

e = U11 ZT I − BQ2 V−1 QT BT .
β

2
22
e =
For any b ∈ Rp with kbk = 1, we can write bT β


n
P

αi i , where

i=1
−1 T
αi2 = n−2 bT (nU11 ) ZT
i − V12 V22 Q2 Bi




−1
Zi − B T
i Q2 V22 V21 (nU11 )b,

24

Li Wang, Guannan Wang, Min-Jun Lai and Lei Gao

and conditioning on {(Xi , Zi ) , i = 1, . . . , n}, αi i ’s are independent. By Lemma A.6, we have
n
o
2
−1 T
max αi2 ≤ Cn−2 max kZi k2 + V12 V22
Q2 Bi
,
1≤i≤n

1≤i≤n

where for any a ∈ Rp ,
−1 T
−1 T
−1
−2 T T
aT V12 V22
Q2 Bi a = n−1 aT V12 (QT
2 Γλ Q2 ) Q2 Bi a ≤ Cn |4| a Z BBi a,

P
P
and the j-th component of n−1 ZT BBi is n1 ni0 =1 Zi0 j ξ∈K Bξ (Xi0 ) Bξ (Xi ). Using Assumpn P
o2
P
= O(1), for large n,
tions (A1) and (A2), we have E n1 ni0 =1 Zi0 j ξ∈K Bξ (Xi0 ) Bξ (Xi )
thus with probability approaching 1,

max

1≤i≤n

n
1 XX
Zi0 j Bξ (Xi0 ) Bξ (Xi ) = OP (1),
n 0
i =1 ξ∈K

−1 T
Q2 Bi k2 = OP (|4|−2 ).
max kV12 V22

1≤i≤n


Therefore, max1≤i≤n αi2 = OP n−2 |4|−2 . Next, with probability approaching 1,
n
X

i
h
e |{(Xi , Zi ) , i = 1, . . . , n}
αi2 = Var bT β


i=1



−1 T T
−1 T T
Q2 B ZU11 bσ 2
Q2 B
I − BQ2 V22
= bT U11 ZT I − BQ2 V22
(
)
n
X
−1 T
−1
T
b i )(Zi − Z
b i)
= n b (nU11 ) n
(Zi − Z
(nU11 ) bσ 2 ,

(A.11)

i=1

b i is the i-th column of ZT BQ2 V−1 QT BT . Using Lemma A.6 again, we have Pn α2 ≥
where Z
2
22
i=1 i

P
cn−1 . So max αi2 / ni=1 αi2 = OP n−1 |4|−2 = oP (1) from Assumption (C4). By Linderberg1≤i≤n

P
Pn
2 −1/2 −→ N (0, 1). Then the desired result follows.
Feller CLT, we have ni=1 αi i /
i=1 αi
2
For any j = 1, . . . , p and λ > 0, define
sλ,zj = argmins∈S

n
X
{zj (Xi ) − s(Xi )}2 + λEυ (s),

(A.12)

i=1

where zj is the coordinate mapping that maps z to its j-th component.
Lemma A.9. Under Assumptions (A1), (A2), (C2), (C3), for sλ,zj defined in (A.12), s0,zj − sλ,zj
OP

(λn−1 |4|−5 ),

j = 1, . . . , p.

Proof. Note that
n zj − sλ,zj , u

n

= λ sλ,zj , u

Eυ

,

zj − s0,zj , u

n

=0

for all u ∈ S,

n

=

25

Efficient Estimation of Partially Linear Models for Spatial Data over Complex Domains

Inserting u = s0,zj − sλ,zj in the above yields that
n s0,zj − sλ,zj

2
n

= λ sλ,zj , s0,zj − sλ,zj
2
Eυ

By Cauchy-Schwarz inequality, sλ,zj

Eυ

= λ(hsλ,zj , s0,zj iEυ − hsλ,zj , sλ,zj iEυ ).

≤ sλ,zj , s0,zj

sλ,zj

Eυ

≤ s0,zj

≤ sλ,zj

Eυ
Eυ

Eυ

s0,zj

Eυ

, which implies

.

(A.13)

By (A.6), we have for large n
2
n

n s0,zj − sλ,zj
thus, s0,zj − sλ,zj

n

≤ s0,zj

Eυ

≤ λ sλ,zj

Eυ

s0,zj − sλ,zj

n

× OP (|4|−2 ),

× OP (λn−1 |4|−2 ). Markov’s inequality implies that
s0,zj

Eυ

≤ C|4|−2 maxξ∈K n−1

≤

C1
s0,zj
|4|2

∞

.

(A.14)

Pn

(Xi ) Zij with probability approaching one.
P
According to Assumptions (A1) and (A2), maxξ∈K n−1 ni=1 Bξ (Xi ) Zij = OP (|4|). The
Note that s0,zj

∞

i=1 Bξ

2

desired results follows.

Lemma A.10. Under Assumptions (A1)-(A3) and 
(C1)-(C5), for the covariance
 matrix Σ
∗
∗
e
defined in (3.5), we have cΣ Ip ≤ Σ ≤ CΣ Ip , and Var β  |{(Xi , Zi ) , i = 1, . . . , n} = n−1 Σ +
oP (1).
Proof. According to (A.11),




(

e |{(Xi , Zi ) , i = 1, . . . , n} = n−1 (nU11 ) n−1
Var β


n
X

)
b i )(Zi − Z
b i )T
(Zi − Z

(nU11 ) σ 2 .

i=1

By the definition of U−1
11 in (2.9), we have
n

−1

(nU11 )



1X
b i )T = hzj , zj 0 − sλ,z 0 in
=
Zi (Zi − Z
.
j
n
1≤j,j 0 ≤p
i=1

As in the proof of Lemma A.7, let e
hj ∈ S and hj satisfy (A.10). Then we have
hzj , zj 0 − sλ,zj 0 in = hzj − e
hj , zj 0 − sλ,zj 0 in +

λ
hsλ,zj 0 , e
hj iEυ .
n

(A.15)

By (A.5), (A.13) and (A.14), we have
hsλ,zj 0 , e
hj 0 iEυ ≤ ksλ,zj 0 kEυ ke
hj 0 kEυ ≤ ks0,zj 0 kEυ ke
hj 0 kEυ

C
CC ∗
F2
e
0
≤
ks
k
k
h
k
≤
h0j 2,∞ +
|4|`+1−υ h0j
0,zj 0 ∞
j Eυ
2
3
|4|
|4|
F1


`+1,∞

.

26

Li Wang, Guannan Wang, Min-Jun Lai and Lei Gao

Note that
hzj − e
hj , zj 0 − sλ,zj 0 in = hzj − hj , zj 0 − hj 0 in + hhj − e
hj , hj 0 − e
hj 0 in + hzj − hj , hj 0 − e
hj 0 in
+ hhj − e
hj , zj 0 − hj 0 in + hzj − hj , e
hj 0 − sλ,zj 0 in + hhj − e
hj , e
hj 0 − sλ,zj 0 in .

(A.16)

According to (A.10), the second term on the right side of (A.16) satisfies that
hj k∞ khj 0 − e
hj 0 k∞ = oP (1).
hhj − e
hj , hj 0 − e
hj 0 i∞ ≤ khj − e
By Lemma A.2 and (A.10), the third term on the right side of (A.16) satisfies that
hzj − hj , hj 0 − e
hj 0 in ≤ {kzj − hj kL2 (1 + oP (1))} khj 0 − e
hj 0 k∞ = oP (1).
Similarly, we have
hhj − e
hj , zj 0 − hj 0 in = oP (1).
From the triangle inequality, we have
hj − hj kn + khj − s0,zj kn + ks0,zj − sλ,zj kn .
ke
hj − sλ,zj kn ≤ ke
According to (A.10) and Lemma A.9, we have
ke
hj − sλ,zj kn ≤ khj − s0,zj kn + oP (1).
Define h∗j,n = argminh∈S kzj − hkL2 , then, from the triangle inequality, we have
khj − s0,zj kn ≤ khj − h∗j,n kn + kh∗j,n − s0,zj kn
Note that khj − h∗j,n kL2 = oP (1). Lemma A.2 implies that khj − h∗j,n kn = oP (1). Next note
that ks0,zj − h∗j,n k2L2 = kzj − s0,zj k2L2 − kzj − h∗j,n k2L2 and kzj − s0,zj kn ≤ kzj − h∗j,n kn . Using
Lemma A.2 again, we have
ks0,zj − h∗j,n k2L2 = oP (kzj − h∗j,n k2L2 ) + oP (kzj − s0,zj k2L2 ).
Since there exists a constant C such that kzj − h∗j,n kL2 ≤ C, so we have
kzj − s0,zj kL2 ≤ kzj − h∗j,n kL2 + kh∗j,n − s0,zj kL2 ≤ C + kh∗j,n − s0,zj kL2 .
Therefore, kh∗j,n − s0,zj kL2 = oP (1). Lemma A.2 implies that kh∗j,n − s0,zj kn = oP (1). As a
consequence,
ks0,zj − hj kn = oP (1).

(A.17)

For the fifth item, by Lemma A.2 and (A.17), we have

hzj − hj , e
hj 0 − sλ,zj 0 in ≤ {kzj − hj kL2 (1 + oP (1))} khj − s0,zj kn + oP (1) = oP (1).

27

Efficient Estimation of Partially Linear Models for Spatial Data over Complex Domains

Similarly, for the sixth item, we have

hhj − e
hj , e
hj 0 − sλ,zj 0 in ≤ khj − e
hj kn khj − s0,zj kn + oP (1) = oP (1).

(A.18)

Combining the above results from (A.15) to (A.18) gives that
hzj , zj 0 − sλ,zj 0 in = hzj − hj , zj 0 − h∗j 0 in + oP (1).
Therefore,
n

(nU11 )

−1

1X
e i )(Zi − Z
e i )T + oP (1) = E[(Zi − Z
e i )(Zi − Z
e i )T ] + oP (1).
(Zi − Z
=
n
i=1

Hence,


e |{(Xi , Zi ) , i = 1, . . . , n} = n−1 Σ−1 + oP (1) .
Var β




A.3. Proof of Theorem 2
Let HZ = I − Z(ZT Z)−1 ZT , then
b = U22 QT BT HZ Y = U22 QT BT HZ g0 + U22 QT BT HZ  = θ
eµ + θ
e .
θ
2
2
2
According to Lemma A.3, ks0,g0 − g0 k∞ ≤ C FF21 |4|`+1 |g0 |`+1,∞ . Denote by γ 0 = Q2 θ 0 the
b − θ0 = θ
eµ − θ 0 + θ
e .
spline coefficients of s0,g0 . Then we have the following decomposition: θ
Note that
eµ − θ 0 = U22 QT BT HZ g0 − θ 0 = U22 QT BT HZ (g0 − BQ2 θ 0 ) − λU22 QT PQ2 θ 0 .
θ
2
2
2
According to (2.10), for any a

T T
T
aT U−1
22 a = a Q2 B HZ B + λP Q2 a.
Since HZ is idempotent, so its eigenvalues πj is either 0 or 1. Without loss of generality we
can arrange the eigenvalues in decreasing order so that πj = 1, j = 1, . . . , m and πj = 1,
j = m + 1, . . . , n. Therefore, we have
m

aT (nU22 )−1 a =

λ T T
1X
T
T
πj aT QT
2 B ej ej BQ2 a + a Q2 PQ2 a,
n
n
j=1

where ej be the indicator vector which is a zero vector except for an entry of one at position
j. Using Markov’s inequality, we have


X
λ 
λ C1
Eυ
aξ Bξ  ≤
C2 kak2 .
n
n |4|2
ξ∈K

28

Li Wang, Guannan Wang, Min-Jun Lai and Lei Gao

Thus, by Conditions (C4) and (C5), naT U22 a ≤ C|4|−2 . Next
T
1/2
kU22 QT
|4|−1 n−1 kBT HZ (g0 − BQ2 θ 0 )k
2 B HZ (g0 − BQ2 θ 0 )k ≤ C

1/2


X
F2
`
1/2
−1 −1 
T
2
|4| |g0 |`+1,∞ ,
≤ C |4| n
= OP
{Bξ HZ (g0 − BQ2 θ 0 )}
F1
ξ∈K

and
λkU22 QT
2 PQ2 θ 0 k

Cλ
Cλ
≤
ks0,g0 kEυ ≤
4
n|4|
n|4|4



F2
`−1
|g0 |2,∞ +
|4| |g0 |`+1,∞ .
F1

Thus,





λ
λ
F2
`
|g
|
+
1
+
|4|
|g
|
0 2,∞
0 `+1,∞ .
n |4|4
n |4|5 F1
e = Pn αi i and α2 = αT U22 QT BT HZ BQ2 U22 α.
For any α with kαk = 1, we write αT θ
2
i
i=1
e − θ 0 k = OP
kθ

Following the same arguments as those in Lemma A.8, we have max1≤i≤n αi2 = OP (n−2 |4|−4 ).
e k ≤ |4|−1 |αT θ
e | = |4|−1 | Pn αi i | = OP (|4|−2 n−1/2 ). Therefore,
Thus, kθ
i=1


b − θ 0 k = OP
kθ




λ
1
λ
F2
`
.
|4| |g0 |`+1,∞ + √
|g0 |2,∞ + 1 +
n|4|2
n |4|4
n |4|5 F1

b we have
Observing that gb(x) = B(x)b
γ = B(x)Q2 θ,
g − s0,g0 kL2 + C
kb
g − g0 kL2 ≤ kb

F2
|4|`+1 |g0 |`+1,∞ .
F1

According to Lemma A.1, we have.


F2
`+1
|4| |g0 |`+1,∞
kb
g − g0 kL2 ≤ C |4|kb
γ − γ 0k +
F1




λ
λ
F2
1
`+1
= OP
|g0 |2,∞, + 1 +
|4|
|g0 |`+1,∞ + √
.
n|4|
n |4|3
n |4|5 F1
The proof is completed.

Efficient Estimation of Partially Linear Models for Spatial Data over Complex Domains

29

References
Awanou, G. and Lai, M. J. and Wenston, P. (2005). The multivariate spline method for
scattered data fitting and numerical solutions of partial differential equations. Wavelets
and splines: Athens 2005 24–74.
Chen, R. and Liang, H. and Wang, J. (2011). On determination of linear components in
additive models. Journal of Nonparametric Statistics 23 367–383.
Craven, P. and Wahba, G. (1979). Smoothing noisy data with spline functions. Numerische
Mathematik 31 377–403.
Eilers, P. (2006). P-spline smoothing on difficult domains. [online] Available at
http://www.statistik.lmu.de/sfb386/workshop/smcs2006/slides/eilers.pdf
Ettinger, B. and Guillas, S. and Lai, M. J. (2015). Bivariate splines for functional regression
models with application to ozone forecasting. Environmetrics Accepted.
Furrer, R. and Nychka, D. and Sainand, S. (2011). Package ‘fields’. R package version 6.6.1..
[online] Available at http://cran.r-project.org/web/packages/fields/fields.pdf.
Golitschek, M. and Schumaker, L. L. (2002). Bounds on projections onto bivariate polynomial
spline spaces with stable local bases. Constructive Approximation 18 241–254.
Green, P. J. and Silverman, B. W. (1994). Nonparametric regression and generalized linear
models. Chapman and Hall, London.
Guillas, S. and Lai, M. J. (2010). Bivariate splines for spatial functional regression models.
Journal of Nonparametric Statistics 22 477–497.
Härdle and W., Liang, H. and Gao, J. T. (2000). Partially linear models. Heidelberg: Springer
Physica-Verlag.
Huang, J. Z. and Zhang, L. and Zhou, L. (2007). Efficient estimation in marginal partially
linear models for longitudinal/clustered data using splines. Scandinavian Journal of
Statistics 34 451–477.
Lai, M. J. (2008). Multivariate splines for data fitting and approximation. Conference Proceedings of the 12th Approximation Theory 210–228.
Lai, M. J. and Schumaker, L. L. (2007). Spline functions on triangulations. Cambridge
University Press.

30

Li Wang, Guannan Wang, Min-Jun Lai and Lei Gao

Lai, M. J. and Schumaker, L. L. (1998). Approximation power of bivariate splines. Advances
in Computational Mathematics 29 595-623.
Lai, M. J. and Wang, L. (2013). Bivariate penalized splines for regression. Statistica Sinica
23 1399-1417.
Li, Y. and Ruppert, D. (2008). On the asymptotics of penalized splines. Biometrika 95
415–436.
Liang, H. and Li, R. (2009). Variable selection for partially linear models with measurement
errors. Journal of the American Statistical Association 104 234–248.
Liu, X. and Guillas, S. and Lai, M. J. (2015). Efficient spatial modeling using the SPDE
approach with bivariate splines. Journal of Computational and Graphical Statistics Accepted.
Liu, X. and Wang, L. and Liang, H. (2011). Estimation and variable selection for semiparametric additive partial linear models. Statistica Sinica 21 1225–1248.
Ma, S. and Song, Q. and Wang, L. (2013). Simultaneous variable selection and estimation in
semiparametric modeling of longitudinal/clustered data. Bernoulli 19 252–274.
Ma, S. and Yang, L. (2011). Spline-backfitted kernel smoothing of partially linear additive
model. Journal of Statistical Planning and Inference 141 204–219.
Marx, B. and Eilers, P. (2005) Multidimensional penalized signal regression. Technometrics
47 13–22.
Pace, R. K. and Barry, R. (1997) Sparse spatial autoregressions. Statistics & Probability
Letters 33 291–297.
Persson, P. O. and Strang, G. (2002). A simple mesh generator in MATLAB. SIAM Review
64 329–345.
Ramsay, T. (2002). Spline smoothing over difficult regions. Journal of the Royal Statistical
Society, Series B 64 307–319.
Sangalli, L. and Ramsay, J. and Ramsay, T. (2013). Spatial spline regression models. Journal
of the Royal Statistical Society, Series B 75 681–703.
Shewchuk, J. R. (1996). Applied Computational Geometry: Towards Geometric Engineering.
Springer-Verlag.

Efficient Estimation of Partially Linear Models for Spatial Data over Complex Domains

31

Wahba, G. (1990). Spline models for observational data. SIAM Publications, Philadelphia.
136, 2506-2534.
Wang, L. and Liu, X. and Liang, H. and Carroll, R. (2011). Estimation and variable selection
for generalized additive partial linear models. Annals of Statistics 39 1827–1851
Wang, H. and Ranalli, M. G. (2007). Low-rank smoothing splines on complicated domains.
Biometrics 63 209–217.
Wang, L. and Xue, L. and Qu, A. and Liang, H. (2014). Estimation and model selection in
generalized additive partial linear models for correlated data with diverging number of
covariates. Annals of Statistics 42 592–624.
Wood, S. N. (2003). Thin plate regression splines. Journal of the Royal Statistical Society,
Series B 65 95–114.
Wood, S. N. and Bravington, M. V. and Hedley, S. L. (2008). Soap film smoothing. Journal
of the Royal Statistical Society, Series B 70 931–955.
Xiao, L. and Huang, L. and Schrack, J. A. and Ferrucci, L. and Zipunnikov, V. and Crainiceanu,
C. (2015). Quantifying the life-time circadian rhythm of physical activity: a covariatedependent functional approach. Biostatistics 16 352–367.
Xiao, L. and Li, Y. and Ruppert, D. (2013). Fast bivariate P-splines: the sandwich smoother.
Journal of the Royal Statistical Society, Series B 75 577–599.
Zhang, H. and Cheng, G. and Liu, Y. (2011). Linear or nonlinear? Automatic structure
discovery for partially linear models. Journal of American Statistical Association 106
1099–1112.
Zhou, L. and Pan, H. (2014). Smoothing noisy data for irregular regions using penalized
bivariate splines on triangulations. Computational Statistics 29 263–281.

Table 1: Root mean squared errors of the estimates in Example 1.
ρ

Method

β1

β2

σε

g(·)

0.0

KRIG
TPS
SOAP
LFE
QFE
BPST (41 )
BPST (42 )

0.0850
0.0801
0.0706
0.0712
0.0720
0.0705
0.0702

0.0591
0.0541
0.0511
0.0515
0.0520
0.0510
0.0509

0.0497
0.0415
0.0284
0.0305
0.0341
0.0285
0.0283

0.4240
0.3268
0.2164
0.1713
0.1800
0.1686
0.1669

0.7

KRIG
TPS
SOAP
LFE
QFE
BPST (41 )
BPST (42 )

0.0867
0.0788
0.0705
0.0746
0.0770
0.0706
0.0702

0.0589
0.0565
0.0545
0.0537
0.0533
0.0556
0.0549

0.0454
0.0377
0.0241
0.0288
0.0341
0.0237
0.0236

0.4207
0.3224
0.2016
0.1673
0.1765
0.1648
0.1624

Table 2: Standard error estimates of the BPST (42 ) coefficients in Example 1.
ρ

Parameter

SEmc

SEmean

SEmedian

SEmad

0.0

β1
β2

0.0702
0.0501

0.0660
0.0539

0.0662
0.0536

0.0032
0.0027

0.7

β1
β2

0.0723
0.0551

0.0660
0.0540

0.0661
0.0540

0.0037
0.0028

Table 3: Root mean squared errors of the estimates in Example 2.
ρ

0.0

0.7

Method

β1

β2

σε

MSPE

KRIG
TPS
SOAP
LFE
QFE
BPST (41 )
BPST (42 )

0.0901
0.0835
0.1106
0.1130
0.0940
0.0839
0.0799

0.0754
0.0690
0.0920
0.0930
0.0821
0.0699
0.0703

0.3369
0.0656
0.3223
0.4999
0.1223
0.0633
0.0425

0.7700
0.7648
0.9140
0.9649
0.8802
0.7647
0.7652

KRIG
TPS
SOAP
LFE
QFE
BPST (41 )
BPST (42 )

0.0983
0.0906
0.1269
0.1075
0.0988
0.0893
0.0891

0.0733
0.0690
0.0882
0.1067
0.0766
0.0692
0.0681

0.3353
0.0612
0.3404
0.4999
0.1316
0.0556
0.0450

0.7552
0.8758
0.9116
0.9496
0.8694
0.7574
0.7553

Table 4: Standard error estimates of the BPST (42 ) coefficients in Example 2.
ρ

Parameter

SEmc

SEmean

SEmedian

SEmad

0.0

β1
β2

0.0803
0.0701

0.0828
0.0673

0.0821
0.0673

0.0059
0.0054

0.7

β1
β2

0.0893
0.0679

0.0893
0.0681

0.0843
0.0683

0.0046
0.0054

Table 5: 5-fold cross validation prediction errors for California housing data.
Method

OLS

KRIG

TPS

SOAP

LFE

QFE

BPST

Prediction Error
(in log scale)

0.157

0.086

0.083

0.080

0.067

0.059

0.054

Variable
MedInc
(95% CI)
log(Age)
(95% CI)
log(AveBedrms)
(95% CI)
log(AveOccup)
(95% CI)
log(Hhd)
(95% CI)

OLS
0.202 (s.e. = 0.001)
(0.200, 0.204)
0.165 (s.e. = 0.005)
(0.155, 0.175)
−0.037 (s.e. = 0.018)
(−0.073, −0.001)
−0.465 (s.e. = 0.022)
(−0.509, −0.421)
0.088 (s.e. = 0.004)
(0.080, 0.096)

0.094
NA
−0.039
NA
0.083
NA
−0.081
NA
−0.003
NA

0.148
NA
−0.022
NA
0.133
NA
−0.276
NA
0.024
NA

Table 6: Coefficients.
KRIG
TPS
SOAP
0.085
NA
−0.037
NA
0.068
NA
−0.063
NA
−0.004
NA

† Estimates and confidence intervals (CI) are reported in the California house value study.

0.139
NA
−0.038
NA
0.112
NA
−0.249
NA
0.011
NA

LFE

0.122
NA
−0.040
NA
0.100
NA
−0.169
NA
0.006
NA

QFE

0.114 (s.e. = 0.001)
(0.112, 0.116)
−0.035 (s.e. = 0.004)
(−0.043, −0.027)
0.095 (s.e. = 0.016)
(0.062, 0.128)
−0.133 (s.e. = 0.018)
(−0.169, −0.097)
0.003 (s.e. = 0.003)
(−0.002, 0.008)

BPST

‡ Note that the standard errors of β̂ are not available for KRIG, TPS, SOAP, LFE and QFE methods, so the corresponding CIs are not available,
and we leave them as “NA”.

1

4

0.8

5
4

3

0.6
2

3

0.4
2

1

0.2

1
0

0

0

-1

-0.2

-1

-2

-0.4

-3

-2

-4

-0.6
-3

-5

-0.8

1
0.8
0.6

3.5

0.4

3

0.2

2.5
0

-1
-1

0.5

1

1.5

2

2.5

3

3.5

4

1

-0.4

0.5

-0.6

0

-0.8

-0.5
-1

(a)

(b)

-1

1

1

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0

0

-0.2

-0.2

-0.4

-0.4

-0.6

-0.6

-0.8

-0.8

-0.5

0

1.5

-0.2

-1
-1

-4

-0.5

2

0

0.5

1

1.5

2

2.5

3

3.5

(c)

-1
-1

-0.5

0

0.5

1

1.5

2

2.5

3

3.5

(d)

1
0.8
0.6
0.4
0.2
0
-0.2
-0.4
-0.6
-0.8
-1
-1

-0.5

0

0.5

1

1.5

2

2.5

3

3.5

(e)
Figure 1: Example 1: (a) true function of g(·); (b) contour map of true function g(·); (c)
sampled location points of replicate 1; (d) first triangulation (41 ) and (e) second triangulation
(42 ) over the domain.

1

4

0.8

3

0.6

1

3

0.8
2

0.6
2

0.4

1

0.4
1

0.2

0.2
0
0

0

0
-1

-0.2

-1

-0.4

-0.2
-0.4

-2

-2

-0.6

-0.6
-3

-0.8
-1
-1

-3

-0.8
-4

-0.5

0

0.5

1

1.5

2

2.5

3

3.5

4

-1
-1

-0.5

0

0.5

1

(a)

1.5

2

2.5

3

3.5

4

(b)

1

4

0.8

3

1

3

0.8
2

0.6

0.6
2

1

0.4

0.4
1

0.2

0.2
0

0

0

-0.2

-1

0
-1

-0.2
-0.4

-0.4

-2

-2

-0.6

-0.6
-3

-0.8
-1
-1

-0.5

0

0.5

1

1.5

2

2.5

3

3.5

-1
-1

4

-3

-0.8

-0.5

0

0.5

1

(c)

1.5

2

2.5

3

3.5

4

(d)

1

4

0.8

3

0.6

1

4

0.8

3

0.6
2

2

0.4

0.4
1

1

0.2

0.2
0

0
-0.2

-1

0

0
-0.2

-1

-0.4

-0.4

-2

-2

-0.6

-0.6
-3

-0.8
-1
-1

-4

-0.5

0

0.5

1

1.5

2

2.5

3

3.5

4

(e)

-3

-0.8
-1
-1

-4

-0.5

0

0.5

1

1.5

2

2.5

3

3.5

4

(f)

1

4

0.8

3

0.6
2

0.4
1

0.2
0

0
-0.2

-1

-0.4
-2

-0.6
-3

-0.8
-1
-1

-4

-0.5

0

0.5

1

1.5

2

2.5

3

3.5

4

(g)
Figure 2: Contour maps for the estimators (ρ = 0.0) in Example 1: (a) KRIG; (b) TPS; (c)
SOAP; (d) LFE; (e) QFE; (f) BPST (41 ) and (g) BPST (42 ).

1

4

0.8

3

0.6

1

4

0.8

3

0.6
2

0.4

2

0.4
1

0.2

1

0.2
0

0

0

0

-0.2

-1

-0.2

-0.4

-2

-0.4

-1

-2

-0.6

-0.6
-3

-3

-0.8

-0.8
-4

-1
-1

-0.5

0

0.5

1

1.5

2

2.5

3

3.5

-4

-1
-1

4

-0.5

0

0.5

1

(a)

1.5

2

2.5

3

3.5

4

(b)

1
0.8

4

1

3

0.8

0.6

3

2

0.6
2
1

0.4

0.4
1

0.2

0.2

0

0

0

0

-1

-1

-0.2
-0.4
-0.6
-0.8
-1
-1

0

0.5

1

1.5

2

2.5

3

3.5

-2

-2

-0.4

-3

-0.6

-4

-0.5

-0.2

-4

-1
-1

4

-3

-0.8

-0.5

0

0.5

1

(c)

1.5

2

2.5

3

3.5

4

(d)

1

3

1

4

0.8

0.8

3

2

0.6

0.6

2

1

0.4
0.2

0

0.4
1

0.2
0

0

0

-1

-0.2

-0.2

-1

-2

-0.4

-0.4
-2

-0.6

-0.6

-3
-3

-0.8

-0.8
-4

-1
-1

-0.5

0

0.5

1

1.5

2

2.5

3

3.5

-1
-1

4

(e)

-0.5

0

0.5

1

1.5

2

2.5

3

3.5

4

(f)

1

4

0.8

3

0.6
2

0.4
1

0.2
0

0

-0.2
-1

-0.4
-2

-0.6
-3

-0.8
-1
-1

-0.5

0

0.5

1

1.5

2

2.5

3

3.5

4

(g)
Figure 3: Contour maps for the estimators (ρ = 0.7) in Example 1: (a) KRIG; (b) TPS; (c)
SOAP; (d) LFE; (e) QFE; (f) BPST (41 ) and (g) BPST (42 ).

1

5

0.9

4

0.8

3

0.7

2

0.6

1

0.5

0

0.4

-1

0.3

-2

0.2

-3

0.1

-4

5

0

-5
0

1
0.8

0.2
0.6

0.4
0.4

0.6
0.2

0.8

0

1

0

0

0.1

0.2

0.3

0.5

0.6

0.7

0.8

0.9

1

(b)

0.0

0.0

0.2

0.2

0.4

0.4

Y Coord

0.6

0.6

0.8

0.8

1.0

1.0

(a)

0.4

0.0

0.2

0.4

0.6

0.8

0.0

1.0

0.2

0.4

0.6

0.8

1.0

X Coord

(c)

(d)

1.2

1.2

1

1

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0

0

-0.2
-0.2

0

0.2

0.4

0.6

(e)

0.8

1

1.2

-0.2
-0.2

0

0.2

0.4

0.6

0.8

1

1.2

(f)

Figure 4: Example 2: (a) true function of g(·); (b) contour map of g(·); (c) gaussian random
field ξ(·) and (d) sampled location points of replicate 1; (e) first triangulation (41 ); and (f)
second triangulation (42 ) over the domain.

<50K
50K-100K
100K-150K
150K-200K
200K-300K
>300K
boundary

(a)

<50K
50K-100K
100K-150K
150K-200K
200K-300K
>300K
boundary

(b)

Figure 5: The median house values: (a) data location points (the colors of dots indicate different
values of the houses); (b) domain triangulation.

(a)

(e)

<50K
50K-100K
100K-150K
150K-200K
200K-300K
>300K
boundary

<50K
50K-100K
100K-150K
150K-200K
200K-300K
>300K
boundary

(b)

(f)

<50K
50K-100K
100K-150K
150K-200K
200K-300K
>300K
boundary

<50K
50K-100K
100K-150K
150K-200K
200K-300K
>300K
boundary

(c)

(g)

<50K
50K-100K
100K-150K
150K-200K
200K-300K
>300K
boundary

<50K
50K-100K
100K-150K
150K-200K
200K-300K
>300K
boundary

(d)

<50K
50K-100K
100K-150K
150K-200K
200K-300K
>300K
boundary

Figure 6: Scatter plots for the estimated house values using methods: (a) OLS; (b) KRIG; (c) TPS; (d) SOAP; (e) LFE; (f)
QFE and (g) BPST.

<-300K
-300K~-200K
-200K~-100K
100K~200K
200K~300K
>300K
boundary

(f)

(b)
<-300K
-300K~-200K
-200K~-100K
100K~200K
200K~300K
>300K
boundary

<-300K
-300K~-200K
-200K~-100K
100K~200K
200K~300K
>300K
boundary

(g)

(c)
<-300K
-300K~-200K
-200K~-100K
100K~200K
200K~300K
>300K
boundary

<-300K
-300K~-200K
-200K~-100K
100K~200K
200K~300K
>300K
boundary

(d)

<-300K
-300K~-200K
-200K~-100K
100K~200K
200K~300K
>300K
boundary

Figure 7: Scatter plots for the differences between the estimated house values and the observed house values using methods:
(a) OLS; (b) KRIG; (c) TPS; (d) SOAP; (e) LFE; (f) QFE and (g) BPST.

(e)

(a)

<-300K
-300K~-200K
-200K~-100K
100K~200K
200K~300K
>300K
boundary

