Smoothing Spline ANOVA Models and their Applications in Complex and Massive Datasets

Abstract

Complex and massive datasets can be easily accessed using the newly developed data acquisition technology. In spite of the fact that the smoothing spline ANOVA models have proven to be useful in a variety of fields, these datasets impose the challenges on the applications of the models. In this chapter, we present a selected review of the smoothing spline ANOVA models and highlight some challenges and opportunities in massive datasets. We review two approaches to significantly reduce the computational costs of fitting the model. One real case study is used to illustrate the performance of the reviewed methods.

1. Introduction

Among the nonparametric models, smoothing splines have been widely used in many real applications. There has been a rich body of literature in smoothing splines such as the additive smoothing spline [1, 2, 3, 4, 5, 6], the interaction smoothing spline [7, 8, 9, 10], and smoothing spline ANOVA (SSANOVA) models [11, 12, 13, 14].

In this chapter, we focus on studying the SSANOVA models. Suppose that the data $\left(y_i,x_i\right)$and $i=1,2,\dots ,n$are independent and identically distributed (i.i.d.) copies of $\left(Y,X\right)$, where $Y\in \mathcal{Y}\subset \mathbb{R}$is the response variable and $X\in \mathcal{X}\subset {\mathbb{R}}^d$is the covariate variable. We consider the regression model:

$$y_i=\eta \left(x_i\right)+{\mathrm{ϵ}}_i,i=1,2,\dots ,n,$$E1

where $y_i$is the response, $\eta$is the nonparametric function varying in an infinite-dimensional space, $x_i={\left(x_{i⟨1⟩},\dots ,x_{i⟨d⟩}\right)}^T$is on the domain $\mathcal{X}\subset {\mathbb{R}}^d$, and ${\mathrm{ϵ}}_i{\displaystyle \stackrel{\mathrm{i}.\mathrm{i}.\mathrm{d}.}{\sim }}N\left(0,{\sigma }^2\right)$. More general cases, in which the conditional distribution of $Y$given $x$, denoted as $Y\mid x$, which follows different distributions instead of the Gaussian distribution, will be discussed later. The nonparametric function $\eta$in (1) can be decomposed into

$$\eta \left(x\right)={\eta }_c+\sum _{j=1}^d{\eta }_j\left(x_{⟨j⟩}\right)+\sum _{k<j}{\eta }_{\mathit{kj}}\left(x_{⟨k⟩},x_{⟨j⟩}\right)+\dots$$

through the functional ANOVA, where ${\eta }_c$is a constant function, ${\eta }_j$is the main effect function of $x_{⟨j⟩}$, ${\eta }_{\mathit{kj}}$is the interaction effect of $x_{⟨k⟩}$and $x_{⟨j⟩}$, and so on.

In the model (1), $\eta$can be estimated by minimizing the following penalized likelihood functional:

$$L\left(\eta \right)+\mathit{\lambda J}\left(\eta \right),$$E2

where $L\left(\eta \right)$is a log likelihood measuring the goodness of fit of $\eta$, $J\left(\eta \right)$is a quadratic functional on $\eta$to quantify its smoothness, and $\lambda$is the smoothing parameter balancing trade-offs between the goodness of fit and the smoothness of $\eta$[11, 12, 13]. The computational complexity of estimating $\eta$by minimizing (2) is of the order $O\left(n^3\right)$for the sample of size $n$. Therefore, the high computational costs render SSANOVA models impractical for massive datasets. In this chapter, we review two methods to lower the computational costs. One approach is through the adaptive basis selection algorithm [14]. By carefully sampling a smaller set of basic functions conditional on the response variables, the adaptive sampling reduces the computational costs to $O\left({\mathit{nn}}^{\ast 2}\right)$, where $n^{\ast }\ll n$is the number of the sampled basis functions. The computational costs can also be reduced by the rounding algorithm [15]. This algorithm can significantly decrease the sample size to $\mu$by rounding the data with a given precision, where $\mu \ll n$. After rounding, the computational costs can be dramatically reduced to $O\left({\mu }^3\right)$.

The rest of the chapter is organized as follows. Section 2 provides a detailed introduction to SSANOVA models and the model estimation. The details of adaptive basis selection algorithm and rounding algorithm are reviewed in Section 3. In Appendix, we demonstrate the numerical implementations using the R software.

2. Smoothing spline ANOVA models

In this section, we first review smoothing spline models and the reproducing kernel Hilbert space. Second, we present how to decompose a nonparametric function on tensor product domains, which lays the theoretical foundation for SSANOVA models. In the end, we show the estimation of SSANOVA models and illustrate the model with a real data example.

2.1. Introduction of smoothing spline models

In the model (1), $\eta$is located in an infinite-dimensional space. One way to estimate it is to add some constraints and estimate $\eta$in a finite-dimensional space. With the smoothness constraint, we estimate $\eta$by minimizing the penalized likelihood functional (2), and the minimizer of (2) is called a smoothing spline.Example 1.

Cubic smoothing splines

Suppose that$Y\mid x$follows a normal distribution, that is,$Y\mid x_i\sim N\left(\eta \left(x_i\right),{\sigma }^2\right)$. Then, the penalized likelihood functional(2) can be reduced as the penalized least squares:

$$\frac{1}{n}\sum _{i=1}^n{\left(y_i-\eta \left(x_i\right)\right)}^2+\lambda {\int }_{\mathcal{X}}{\left(\ddot{\eta }\left(x\right)\right)}^2\mathit{dx},$$E3

where$\ddot{\eta }=d^2\eta /{\mathit{dx}}^2$. The minimizer of(3) is called a cubic smoothing spline[16, 17, 18]. In(3), the first term quantifies the fidelity to the data, and the second term controls the roughness of the function.

Example 2.

Exponential family smoothing splines

Suppose that$Y\mid x$follows an exponential family distribution:

$$Y\mid x_i\sim \text{exp}\left\{\left(\mathit{y\eta }\left(x_i\right)-b\left(\eta \left(x_i\right)\right)\right)/a\left(\varphi \right)+c\left(y,\varphi \right)\right\},$$

where$a>0$, $b$, and$c$are known functions and$\varphi$is either known or a nuisance parameter. Then,$\eta$can be estimated by minimizing the following penalized likelihood functional[19, 20]:

$$-\frac{1}{n}\sum _{i=1}^n\left\{y_i\eta \left(x_i\right)-b\left(\eta \left(x_i\right)\right)\right\}+\mathit{\lambda J}\left(\eta \right).$$E4

Note that the cubic smoothing spline inExample 1 is a special case of exponential family smoothing splines when$Y\mid x$follows the Gaussian distribution.

The smoothing parameter $\lambda$is sensitive to the estimation of $\eta$(see Figure 1). Therefore, it is crucial to implement some proper smoothing parameter selection methods to estimate $\lambda$. One of the most popular methods is the generalized cross validation (GCV) [21, 22]. More details will be discussed in Section 2.6.2.

Figure 1.

The data are generated by the model ϵ.$y=5+e^{3x_1}+{10}^6{x}_2^{11}{\left(1-x_2\right)}^6+{10}^4{x}_2^2{\left(1-x_2\right)}^{10}+5\text{cos}\left(2\pi \left(x_1-x_2\right)\right)+\mathrm{ϵ}$, where$\mathrm{ϵ}\sim N\left(0,1\right)$. Panels (a) and (b) show the data and the true function, respectively. The estimated functions depending on different smoothing parameters are shown in panels (c), (d), and (e). We set$\lambda \to 0$in panel (c) and$\lambda \to \infty$in panel (e). The proper$\lambda$selected by generalized cross validation (GCV) is used in panel (d).

2.2. Reproducing kernel Hilbert space

We assume that readers are familiar with Hilbert space, which is a complete vector space with an inner product well defined that allows length and angle to be measured [23]. In a general Hilbert space, the continuity of a functional, which is required in minimizing (2) on $\mathcal{H}=\left\{J\left(\eta \right)<\infty \right\}$, is not always satisfied. To circumvent the problem, we optimize (2) in a special Hilbert space named reproducing kernel Hilbert space [24].

For each $g\in \mathcal{H}$, there exists a corresponding continuous linear functional $L_g$such that $L_g\left(f\right)=⟨g,f⟩$, where $f\in \mathcal{H}$and $⟨\cdot ,\cdot ⟩$defines the inner product in $\mathcal{H}$. Conversely, an element $g_L\in \mathcal{H}$can also be found such that $⟨g_L,f⟩=L\left(f\right)$for any continuous linear functional $L$of $\mathcal{H}$[23]. This is known as the Riesz representation theorem.Theorem 2.1.

Riesz representation

Let$\mathcal{H}$be a Hilbert space. For any functional$L$of$\mathcal{H}$, there uniquely exists an element$g_L\in \mathcal{H}$such that

$$L\left(\cdot \right)=⟨g_L,\cdot ⟩,$$

where$g_L$is called the representer of$L$. The uniqueness is in the sense that$g_1$and$g_2$are considered as the same representer for any$g_1$and$g_2$satisfying$\Vert g_1-g_2\Vert =0$, where$\parallel \cdot \parallel =⟨\cdot ,\cdot ⟩$defines the norm in$\mathcal{H}$.

For a better construction of estimator minimizing (2), one needs the continuity of evaluation functional $\left[x\right]f=f\left(x\right)$. Roughly speaking, this means that if two functions $f$and $g$are close in norm, that is, $\Vert f-g\Vert$is small, then $f$and $g$are also pointwise close, that is, $\mid f\left(x\right)-g\left(x\right)\mid$is small for all $x$.Definition 1.

Reproducing kernel Hilbert space

Consider a Hilbert space$\mathcal{H}$consisting of functions on domain$\mathcal{X}$. For every element$x\in \mathcal{X}$, define an evaluation functional$\left[x\right]$such that$\left[x\right]f=f\left(x\right)$. If all the evaluation functional$\left[x\right]$s are continuous,$\forall x\in \mathcal{X}$, then$\mathcal{H}$is called a reproducing kernel Hilbert space.

By Theorem 2.1, for every evaluation functional $\left[x\right]$, there exists a corresponding function $R_x\in \mathcal{H}$on $\mathcal{X}$as the representer of $\left[x\right]$, such that $⟨R_x,f⟩=\left[x\right]f=f\left(x\right)$and $\forall f\in \mathcal{H}$. By the definition of evaluation functional, it follows

$$R_x\left(y\right)=⟨R_x,R_y⟩=R_y\left(x\right).$$E5

The bivariate function $R\left(x,y\right)=⟨R_x,R_y⟩$is called the reproducing kernel of $\mathcal{H}$, which is unique if it exists. The essential meaning of the name “reproducing kernel” comes from its reproducing property

$$⟨R\left(x,\cdot \right),f⟩=⟨R_x\left(\cdot \right),f⟩=f\left(x\right)$$

for any $f\in \mathcal{H}$. In general, a reproducing kernel Hilbert space defines a reproducing kernel function that is both symmetric and positive definite. In addition, Moore-Aronszajn theorem states that every symmetric, positive definite kernel defines a unique reproducing kernel Hilbert space [25], and hence one can construct a reproducing kernel Hilbert space simply by specifying its reproducing kernel.

We now introduce the concept of tensor sum decomposition. Suppose that $\mathcal{H}$is a Hilbert space and $\mathcal{G}$is a linear subspace of $\mathcal{H}$. The linear subspace ${\mathcal{G}}^c=\left\{f\in \mathcal{H},⟨f,g⟩=0,\forall g\in \mathcal{G}\right\}$is called the orthogonal complement of $\mathcal{G}$. It is easy to verify that for any $f\in \mathcal{H}$, there exists a unique decomposition $f=f_{\mathcal{G}}+f_{{\mathcal{G}}^c}$, where $f_{\mathcal{G}}\in \mathcal{G}$and $f_{{\mathcal{G}}^c}\in {\mathcal{G}}^c$. This decomposition is called a tensor sum decomposition, denoted by $\mathcal{H}=\mathcal{G}\oplus {\mathcal{G}}^c$. Suppose that ${\mathcal{H}}_1$and ${\mathcal{H}}_2$are two Hilbert spaces with inner products ${⟨\cdot ,\cdot ⟩}_1$and ${⟨\cdot ,\cdot ⟩}_2$. If the only common element of these two spaces is $\mathbf{0}$, one can also define a tensor sum Hilbert space $\mathcal{H}={\mathcal{H}}_1\oplus {\mathcal{H}}_2$. For any $f,g\in \mathcal{H}$, one has unique decompositions $f=f_1+f_2$and $g=g_1+g_2$, where $f_1,g_1\in {\mathcal{H}}_1$and $f_2,g_2\in {\mathcal{H}}_2$. Moreover, the inner product defined on $\mathcal{H}$would be $⟨f,g⟩={⟨f_1,g_1⟩}_1+{⟨f_2,g_2⟩}_2$. The following theorem provides the rules in the tensor sum decomposition of a reproducing kernel Hilbert space.Theorem 2.2

Suppose that$R_1$and$R_2$are the reproducing kernel Hilbert spaces${\mathcal{H}}_1$and${\mathcal{H}}_2$, respectively. If${\mathcal{H}}_1\cap {\mathcal{H}}_2=\left\{\mathbf{0}\right\}$, then$\mathcal{H}={\mathcal{H}}_1\oplus {\mathcal{H}}_2$has a reproducing kernel$R=R_1+R_2$.

Conversely, if the reproducing kernel$R$of$\mathcal{H}$can be decomposed into$R=R_1+R_2$, where both$R_1$and$R_2$are positive definite, and they are orthogonal to each other, that is, $⟨R_1\left(x_1,\cdot \right),R_2\left(x_2,\cdot \right)⟩=0$for$\forall x_1,x_2\in \mathcal{X}$, then the spaces${\mathcal{H}}_1$and${\mathcal{H}}_2$corresponding to the kernels$R_1$and$R_2$form a tensor sum decomposition$\mathcal{H}={\mathcal{H}}_1\oplus {\mathcal{H}}_2$.

2.3. Representer theorem

In (2), the smoothness penalty term $J\left(\eta \right)=J\left(\eta ,\eta \right)$is nonnegative definite, that is, $J\left(\eta ,\eta \right)\ge 0$, and hence it is a squared semi-norm on the reproducing kernel Hilbert space $\mathcal{H}=\left\{\eta :J\left(\eta \right)<\infty \right\}$. Denote ${\mathcal{N}}_J=\left\{\eta \in \mathcal{H}:J\left(\eta \right)=0\right\}$as the null space of $J\left(\eta \right)$and ${\mathcal{H}}_J$as its orthogonal complement. By the tensor sum decomposition $\mathcal{H}={\mathcal{N}}_J\oplus {\mathcal{H}}_J$, one may decompose the $\eta$into two parts: one in the null space ${\mathcal{N}}_J$that has no contribution on the smoothness penalty and the other in ${\mathcal{H}}_J$“reproduced” by the reproducing kernel $R\left(\cdot ,\cdot \right)$[12].Theorem 2.3.

There exist coefficient vectors$\mathbf{d}={\left(d_1,\dots ,d_M\right)}^T\in {\mathbb{R}}^M$and$\mathbf{c}={\left(c_1,\dots ,c_n\right)}^T\in {\mathbb{R}}^n$such that

$$\eta \left(x\right)=\sum _{k=1}^M{d}_k{\xi }_k\left(x\right)+\sum _{i=1}^n{c}_{i}R\left(x_i,x\right),$$E6

where$\left\{{\xi }_k,k=1,..,M\right\}$is the basis of null space${\mathcal{N}}_J$and$R\left(\cdot ,\cdot \right)$is the reproducing kernel of${\mathcal{H}}_J$.

This theorem indicates that although the minimization problem is in an infinite-dimensional space, the minimizer of (2) lies in a data-adaptive finite-dimensional space.

2.4. Function decomposition

The decomposition of a multivariate function is similar to the classical ANOVA. In this section, we present the functional ANOVA which lays the foundation for SSANOVA models.

2.4.1. One-way ANOVA decomposition

We consider a classical one-way ANOVA model $y_{\mathit{ij}}={\mu }_i+{\mathrm{ϵ}}_{\mathit{ij}}$, where $y_{\mathit{ij}}$is the observed data, ${\mu }_i$is the treatment mean for $i=1,\dots ,K$and $j=1,\dots ,J$, and ${\mathrm{ϵ}}_{\mathit{ij}}$s are the random errors. The treatment mean ${\mu }_i$can be further decomposed as ${\mu }_i=\mu +{\alpha }_i$, where $\mu$is the overall mean and ${\alpha }_i$is the treatment effect with the constraint ${\sum }_{i=1}^K{\alpha }_i=0$.

Similar to the classical ANOVA decomposition, a univariate function $f\left(x\right)$can be decomposed as

$$f=\mathit{Af}+\left(I-A\right)f=f_c+f_x,$$E7

where $A$is an averaging operator that averages the effect of $x$and $I$is an identity operator. The operator $A$averages a function $f$to a constant function $f_c$satisfying $A\left(I-A\right)=0$. For example, one can take $\mathit{Af}={\int }_0^{1}f\left(x\right)\mathit{dx}$in ${\mathcal{L}}_1\left[0,1\right]=\left\{f:{\int }_0^1|f\left(x\right)|\mathit{dx}<\infty \right\}$. In (7), $f_c=\mathit{Af}$is the mean function, and $f_x=\left(I-A\right)f$is the treatment effect.

2.4.2. Multiway ANOVA decomposition

On a $d$-dimensional product domain $\mathcal{X}={\prod }_{j=1}^d{\mathcal{X}}_j\in {\mathbb{R}}^d$, a multivariate function $f\left(x_{⟨1⟩},\dots ,x_{⟨d⟩}\right)$can be decomposed similarly to the one-way ANOVA decomposition. Let $A_j$, $j=1,\dots ,d$, be the average operator on ${\mathcal{X}}_j$, and then $A_{j}f$is a constant function on ${\mathcal{X}}_j$. One can define the ANOVA decomposition on $\mathcal{X}$as

$$\begin{array}{ll}f& =\left\{\prod _{j=1}^d\left(I-A_j+A_j\right)\right\}f\\ & =\sum _{\mathcal{S}}\left\{\prod _{j\in \mathcal{S}}\left(I-A_j\right)\prod _{j\notin \mathcal{S}}{A}_j\right\}f=\sum _{\mathcal{S}}{f}_{\mathcal{S}},\end{array}$$E8

where $\mathcal{S}\subseteq \left\{1,\dots ,d\right\}$. The term $f_c={\prod }_{j=1}^d{A}_{j}f$is the constant function, $f_j=\left(I-A_j\right){\prod }_{\alpha \ne j}{A}_{\alpha }f$is the main effect term of $x_{⟨j⟩}$, the term $f_{\mathit{\mu \nu }}=\left(I-A_{\mu }\right)\left(I-A_{\nu }\right){\prod }_{\alpha \ne \mu ,\nu }{A}_{\alpha }f$is the interaction of $x_{⟨\mu ⟩}$and $x_{⟨\nu ⟩}$, and so on.

2.5. Some examples of model conduction

Smoothing splines on${\mathcal{C}}^{\left(m\right)}\left[0,1\right]$. If we define

$$J\left(\eta \right)={\int }_0^1{\left({\eta }^{\left(m\right)}\right)}^2\mathit{dx}$$

in the space ${\mathcal{C}}^{\left(m\right)}\left[0,1\right]=\left\{f:f^{\left(m\right)}\in {\mathcal{L}}_2\left[0,1\right]\right\}$, where $f^{\left(m\right)}$denotes the $m$th differentiation of $f$, ${\mathcal{L}}_2=\left\{f:{\int }_0^1{\left(f\left(x\right)\right)}^2\mathit{dx}<\infty \right\}$, then the minimizer of (2) is called a polynomial smoothing spline.

Here, we use an inner product

$$⟨f,g⟩=\sum _{\nu =0}^{m-1}\left({\int }_0^1{f}^{\left(\nu \right)}\left(x\right)\mathit{dx}\right)\left({\int }_0^1{g}^{\left(\nu \right)}\left(x\right)\mathit{dx}\right)+{\int }_0^1{f}^{\left(m\right)}\left(x\right)g^{\left(m\right)}\left(x\right)\mathit{dx}.$$E9

One can easily check that (9) is a well-defined inner product in ${\mathcal{C}}^{\left(m\right)}\left[0,1\right]$with ${\mathcal{H}}_0=\left\{f:f^{\left(m\right)}=0\right\}$equipped with the inner product ${\sum }_{\nu =0}^{m-1}\left({\int }_0^1{f}^{\left(\nu \right)}\left(x\right)\mathit{dx}\right)\left({\int }_0^1{g}^{\left(\nu \right)}\left(x\right)\mathit{dx}\right)$[21]. To construct the reproducing kernel, define

$$k_{\nu }\left(x\right)=-\left(\sum _{\mu =-\infty }^{-1}+\sum _{\mu =1}^{\infty }\right)\frac{\text{exp}\left(2\pi \mathbf{i}\mathit{\mu x}\right)}{{\left(2\pi \mathbf{i}\mu \right)}^{\nu }},$$E10

where $\mathbf{i}=\sqrt{-1}$. One can verify that ${\int }_0^1{k}_{\nu }^{\left(\mu \right)}\left(x\right)\mathit{dx}={\delta }_{\mathit{\nu \mu }}$and $\nu ,\mu =0,1,\dots ,m-1$, where ${\delta }_{\mathit{\nu \mu }}$is the Kronecker delta [26]. Indeed, $\left\{k_0,\dots ,k_{m-1}\right\}$forms an orthonormal basis of ${\mathcal{H}}_0$. Then, the reproducing kernel in ${\mathcal{H}}_0$is

$$R_0\left(x,y\right)=\sum _{\nu =0}^{m-1}{k}_{\nu }\left(x\right)k_{\nu }\left(y\right).$$

For space

$${\mathcal{H}}_1=\left\{f:{\int }_0^1{f}^{\left(\nu \right)}\left(x\right)\mathit{dx}=0,\nu =0,1,\dots ,m-1,f^{\left(m\right)}\in {\mathcal{L}}_2\left[0,1\right]\right\},$$

one can check that the reproducing kernel in ${\mathcal{H}}_1$is

$$R_1\left(x,y\right)=k_m\left(x\right)k_m\left(y\right)+{\left(-1\right)}^{m-1}{k}_{2m}\left(x-y\right),$$

(See details in [11]; Section~2.3).

SSANOVA models on product domains: A natural way to construct reproducing kernel Hilbert space on product domain ${\prod }_{j=1}^d{\mathcal{X}}_j$is taking the tensor product of spaces constructed on the marginal domains ${\mathcal{X}}_j$s. According to the Moore-Aronszajn theorem, every nonnegative definite function $R$corresponds to a reproducing kernel Hilbert space with $R$as its reproducing kernel. Therefore, the construction of the tensor product reproducing kernel Hilbert space is induced by constructing its reproducing kernel.Theorem 2.4.

Suppose that if$R_{⟨1⟩}\left(x_{⟨1⟩},y_{⟨1⟩}\right)$is nonnegative definite on${\mathcal{X}}_1$and$R_{⟨2⟩}\left(x_{⟨2⟩},y_{⟨2⟩}\right)$is nonnegative definite on${\mathcal{X}}_2$, then$R\left(x,y\right)=R_{⟨1⟩}\left(x_{⟨1⟩},y_{⟨1⟩}\right)R_{⟨2⟩}\left(x_{⟨2⟩},y_{⟨2⟩}\right)$is nonnegative definite on$\mathcal{X}={\mathcal{X}}_1\times {\mathcal{X}}_2$.

Theorem 2.4 implies that a reproducing kernel $R$on tensor product reproducing kernel Hilbert space can be derived from the reproducing kernels on marginal domains. Indeed, let ${\mathcal{H}}_{⟨j⟩}$be the space on ${\mathcal{X}}_j$with reproducing kernel $R_{⟨j⟩}$, where $j=1,2$. Then, $R=R_{⟨1⟩}{R}_{⟨2⟩}$is nonnegative definite on ${\mathcal{X}}_1\times {\mathcal{X}}_2$. The space $\mathcal{H}$corresponding to $R\left(\cdot ,\cdot \right)$is called the tensor product space of ${\mathcal{H}}_{⟨1⟩}$and ${\mathcal{H}}_{⟨2⟩}$, denoted by $\mathcal{H}={\mathcal{H}}_{⟨1⟩}\otimes {\mathcal{H}}_{⟨2⟩}$.

One can decompose each marginal space ${\mathcal{H}}_{⟨j⟩}$into ${\mathcal{H}}_{⟨j⟩}={\mathcal{H}}_{⟨j⟩0}\oplus {\mathcal{H}}_{⟨j⟩1}$, where ${\mathcal{H}}_{⟨j⟩0}$denotes the averaging space and ${\mathcal{H}}_{⟨j⟩1}$denotes the orthogonal complement. Then, by the discussion in Section 2.4, the one-way ANOVA decomposition on each marginal space can be generalized to a multidimensional space $\mathcal{H}={\otimes }_{j=1}^d{\mathcal{H}}_{⟨j⟩}$as

$$\begin{array}{ll}\mathcal{H}& ={\otimes }_{j=1}^d\left({\mathcal{H}}_{⟨j⟩0}\oplus {\mathcal{H}}_{⟨j⟩1}\right)\\ & ={\oplus }_{\mathcal{S}}\left\{\left({\otimes }_{j\in \mathcal{S}}{\mathcal{H}}_{⟨j⟩1}\right)\otimes \left({\otimes }_{j\notin \mathcal{S}}{\mathcal{H}}_{⟨j⟩0}\right)\right\}\\ & ={\oplus }_{\mathcal{S}}{\mathcal{H}}_{\mathcal{S}},\end{array}$$E11

where $\mathcal{S}$denotes all the subsets of $\left\{1,\dots ,d\right\}$. The component $f_{\mathcal{S}}$in (8) is in the space ${\mathcal{H}}_{\mathcal{S}}$. Based on the decomposition, the minimizer of (2) is called a tensor product smoothing spline. One can construct a tensor product smoothing spline following Theorem 2.3, in which the reproducing kernel term may be calculated in the same way as the tensor product (11).

In the following, we will give some examples of tensor product smoothing splines on product domains.

2.5.1. Smoothing splines on $\left\{\mathit{1},\dots ,K\right\}\times \left[\mathit{0},\mathit{1}\right]$

We construct the reproducing kernel Hilbert space by using

$$R_{⟨1⟩0}=1/K\text{and}{R}_{⟨1⟩1}=I_{\left[x_{⟨1⟩}=y_{⟨1⟩}\right]}$$

on $\left\{1,\dots ,K\right\}$. On $\left[0,1\right]$, assume that if $m=2$, then we have

$$R_{⟨2⟩0}=1+k_1\left(x_{⟨2⟩}\right)k_1\left(y_{⟨2⟩}\right)$$

and

$$R_{⟨2⟩1}=k_2\left(x_{⟨2⟩}\right)k_2\left(y_{⟨2⟩}\right)-k_4\left(x_{⟨2⟩}-y_{⟨2⟩}\right).$$

In this case, the space $\mathcal{H}$can be further decomposed as

$$\mathcal{H}=\left({\mathcal{H}}_{⟨1⟩0}\oplus {\mathcal{H}}_{⟨1⟩1}\right)\otimes \left({\mathcal{H}}_{⟨2⟩00}\oplus {\mathcal{H}}_{⟨2⟩01}\oplus {\mathcal{H}}_{⟨2⟩1}\right).$$E12

The reproducing kernels of tensor product cubic spline on $\left\{1,\dots ,K\right\}\times \left[0,1\right]$are listed in Table 1.

Subspace	Reproducing kernel
${\mathcal{H}}_{⟨1⟩0}\otimes {\mathcal{H}}_{⟨2⟩00}$	$1/K$
${\mathcal{H}}_{⟨1⟩0}\otimes {\mathcal{H}}_{⟨2⟩01}$	$k_1\left(x_{⟨2⟩}\right)k_1\left(y_{⟨2⟩}\right)/K$
${\mathcal{H}}_{⟨1⟩0}\otimes {\mathcal{H}}_{⟨2⟩1}$	$\left[k_2\left(x_{⟨2⟩}\right)k_2\left(y_{⟨2⟩}\right)-k_4\left(x_{⟨2⟩}-y_{⟨2⟩}\right)\right]/K$
${\mathcal{H}}_{⟨1⟩1}\otimes {\mathcal{H}}_{⟨2⟩00}$	$I_{\left[x_{⟨1⟩}=y_{⟨1⟩}\right]}-1/K$
${\mathcal{H}}_{⟨1⟩1}\otimes {\mathcal{H}}_{⟨2⟩01}$	$\left[I_{\left[x_{⟨1⟩}=y_{⟨1⟩}\right]}-1/K\right]k_1\left(x_{⟨2⟩}\right)k_1\left(y_{⟨2⟩}\right)$
${\mathcal{H}}_{⟨1⟩1}\otimes {\mathcal{H}}_{⟨2⟩1}$	$\left[I_{\left[x_{⟨1⟩}=y_{⟨1⟩}\right]}-1/K\right]\left[k_2\left(x_{⟨2⟩}\right)k_2\left(y_{⟨2⟩}\right)-k_4\left(x_{⟨2⟩}-y_{⟨2⟩}\right)\right]$

Table 1.

Reproducing kernels of (12) on $\left\{1,\dots ,K\right\}\times \left[0,1\right]$when $m=2$.

On other product domains, for example, ${\left[0,1\right]}^2$, the tensor product reproducing kernel Hilbert space can be decomposed in a similar way. More examples are available in ([11], Section~2.4).

2.5.1.1. General form

In general, a tensor product reproducing kernel Hilbert space can be specified as $\mathcal{H}={\oplus }_j{\mathcal{H}}_j$, where $j\in B$is a genetic index. Suppose that ${\mathcal{H}}_j$is equipped with a reproducing kernel $R_j$and an inner product ${⟨f,g⟩}_j$. Denote $f_j$as the projection of $f$onto ${\mathcal{H}}_j$. Then, an inner product in $\mathcal{H}$can be defined as

$$J\left(f,g\right)=\sum _j{\theta }_j^{-1}{⟨f_j,g_j⟩}_j,$$E13

where ${\theta }_j\ge 0$are the tuning parameters. If a penalty $J$in (2) has the form (13), the SSANOVA models can be defined on the space $\mathcal{H}={\oplus }_j{\mathcal{H}}_j$with the reproducing kernel:

$$R=\sum _j{\theta }_j{R}_j.$$E14

2.6. Estimation

In this section, we show the procedure of estimating the minimizer $\widehat{\eta }$of (2) under the Gaussian assumption and selecting the smoothing parameters.

2.6.1. Penalized least squares

We consider the same model shown in (1), and then the $\eta$can be estimated by minimizing the penalized least squares:

$$\frac{1}{n}\sum _{i=1}^n{\left(y_i-\eta \left(x_i\right)\right)}^2+\mathit{\lambda J}\left(\eta \right).$$E15

Let $S$denote the $n\times M$matrix with the $\left(i,j\right)$th entry ${\xi }_j\left(x_i\right)$as in (6) and $R$denote the $n\times n$matrix with the $\left(i,j\right)$th entry $R\left(x_i,x_j\right)$with the form (14). Then, based on Theorem 2.3, $\mathit{\eta }$can be expressed as

$$\mathit{\eta }=S\mathbf{d}+R\mathbf{c},$$

where $\mathit{\eta }={\left(\eta \left(x_1\right),\dots ,\eta \left(x_n\right)\right)}^T$, $\mathbf{d}={\left(d_1,\dots ,d_M\right)}^T$, and $\mathbf{c}={\left(c_1,\dots ,c_n\right)}^T$. The least squares term in (15) becomes

$$\frac{1}{n}\sum _{i=1}^n{\left(y_i-\eta \left(x_i\right)\right)}^2=\frac{1}{n}{\left(\mathbf{y}-S\mathbf{d}-R\mathbf{c}\right)}^T\left(\mathbf{y}-S\mathbf{d}-R\mathbf{c}\right),$$

where $\mathbf{y}={\left(y_1,\dots ,y_n\right)}^T$.

By the reproducing property (5), the roughness penalty term can be expressed as

$$J\left(\eta \right)=\sum _{i=1}^n\sum _{j=1}^n{c}_{i}R\left(x_i,x_j\right)c_j={\mathbf{c}}^{T}R\mathbf{c}.$$

Therefore, the penalized least squares criterion (15) becomes

$$\frac{1}{n}{\left(\mathbf{y}-S\mathbf{d}-R\mathbf{c}\right)}^T\left(\mathbf{y}-S\mathbf{d}-R\mathbf{c}\right)+\lambda {\mathbf{c}}^{T}R\mathbf{c}.$$E16

The penalized least squares (16) is a quadratic form of both $\mathbf{d}$and $\mathbf{c}$. By differentiating (16), one can obtain the linear system:

$$\left(\begin{array}{ll}{S}^{T}S& S^{T}R\\ R^{T}S& R^{T}R+\mathit{n\lambda R}\end{array}\right)\left(\begin{array}{l}\mathbf{d}\\ \mathbf{c}\end{array}\right)=\left(\begin{array}{l}{S}^T\mathbf{y}\\ R^T\mathbf{y}\end{array}\right).$$E17

Note that (17) only works for penalized least squares (15), and hence a normal assumption is needed in this case.

2.6.2. Selection of smoothing parameters

In SSANOVA models, properly selecting smoothing parameters is important to estimate $\eta$[9, 27, 28], as shown in Figure 1. Here, we introduce the generalized cross validation (GCV) method for the smoothing parameter selection.

For the multivariate predictors, the penalty term in (15) has the form

$$\mathit{\lambda J}\left(f\right)=\lambda \sum _{j=1}^S{\theta }_j^{-1}{⟨f_j,f_j⟩}_j,$$

where $S$is the number of smoothing parameters, which is related to the functional ANOVA decomposition, and ${\theta }_j$'s are the extra smoothing parameters. To avoid overparameterization, we treat $\lambda ={\left(\lambda /{\theta }_1,\cdots ,\lambda /{\theta }_S\right)}^T$as the effective smoothing parameters.

A GCV score is defined as

$$V\left(\lambda \right)=\frac{n^{-1}{\mathbf{y}}^T{\left(I-A\left(\lambda \right)\right)}^2\mathbf{y}}{{\left[n^{-1}\mathrm{tr}\left(I-A\left(\lambda \right)\right)\right]}^2},$$

where $A\left(\lambda \right)$is a symmetric matrix similar to the hat matrix in linear regression. We can select a proper $\lambda$by minimizing the GCV score [21].

2.7. Case study: Twitter data

Tweets in the contiguous United States were collected over five weekdays in January 2014. The dataset contains information of time, GPS location, and tweet counts (see Figure 2). To illustrate the application of SSANOVA models, we study the time and spatial patterns in this data.

Figure 2.

Heatmaps of tweet counts in the contiguous United States. (a) Tweet counts at 2:00 a.m. (b) Tweet counts at 6:00 p.m.

The bivariate function $\eta \left(x_{⟨1⟩},x_{⟨2⟩}\right)$is a function of time and location, where $x_{⟨1⟩}$denotes the time and $x_{⟨2⟩}$represents the longitude and latitude coordinates. We use the thin-plate spline for the spatial variable and cubic spline for the time variable. As a rotation-free method, the thin-plate spline is popular for modeling spatial data [29, 30, 31]. For a better interpretation, we decompose the function $\eta$as

$$\eta \left(x_{⟨1⟩},x_{⟨2⟩}\right)={\eta }_c+{\eta }_1\left(x_{⟨1⟩}\right)+{\eta }_2\left(x_{⟨2⟩}\right)+{\eta }_{12}\left(x_{⟨1⟩},x_{⟨2⟩}\right),$$

where ${\eta }_c$is a constant function; ${\eta }_1$and ${\eta }_2$are the main effects of time and location, respectively; and ${\eta }_{12}$is the spatial-time interaction effect.

The main effects of time and location are shown in Figure 3. Obviously, in panel (a), the number of tweets has the periodic effect, where it attains the maximum value at 8:00 p.m. and the minimum value at 5:00 a.m. The main effect of time shows the variations of Twitter usages in the United States. In addition, we can infer how the tweet counts vary across different locations based on panel (b) in Figure 3. There tend to be more tweets in the east than those in the west regions and more tweets in the coastal zone than those in the inland. We use the scaled dot product

$$\pi ={\left({\widehat{\mathit{\eta }}}_{12}\right)}^T\widehat{\mathbf{y}}/\parallel \widehat{\mathbf{y}}{\parallel }^2$$

to quantify the percentage decomposition of the sum of squares of $\widehat{\mathbf{y}}$[11], where $\widehat{\mathbf{y}}={\left({\widehat{y}}_1,\dots ,{\widehat{y}}_n\right)}^T$is the predicted values of $\text{log}\left(\#\text{Tweets}\right)$, and ${\widehat{\mathit{\eta }}}_{12}={\left({\eta }_{12}\left(x_1\right),\dots ,{\eta }_{12}\left(x_n\right)\right)}^T$is the estimated interaction effect term, where ${\eta }_{12}\left(x\right)={\eta }_{12}\left(x_{⟨1⟩},x_{⟨2⟩}\right)$. In our fitted model, $\pi =3\times {10}^{-16}$, which is so small that the interaction term is negligible. This indicates that there is no significant difference for the Twitter usages across time in the contiguous United States.

Figure 3.

(a) The main effect function of time (hours). (b) The main effect function of location.

3. Efficient approximation algorithm in massive datasets

In this section, we consider SSANOVA models under the big data settings. The computational cost of solving (17) is of the order $O\left(n^3\right)$and thus gives rise to a challenge on the application of SSANOVA models when the volume of data grows. To reduce the computational load, an obvious way is to select a subset of basis functions randomly. However, it is hard to keep the data features by uniform sampling. In the following section, we present an adaptive basis selection method and show its advantages over uniform sampling [14]. Instead of selecting basis functions, another approach to reduce the computational cost is shrinking the original sample size by rounding algorithm [15].

3.1. Adaptive basis selection

A natural way to select the basis functions is through uniform sampling. Suppose that we randomly select a subset $\stackrel{⌣}{x}=\left\{{\stackrel{⌣}{x}}_1,\dots ,{\stackrel{⌣}{x}}_{\stackrel{⌣}{n}}\right\}$from $\left\{x_1,\dots ,x_n\right\}$, where $\stackrel{⌣}{n}$is the subsample size. Thus, the kernel matrix would be $R\left({\stackrel{⌣}{x}}_i,x\right)$, $i=1,\dots ,\stackrel{⌣}{n}$. Then, one minimizes (17) in the effective model space:

$${\mathcal{H}}_E=\mathcal{N}\oplus \text{span}\left\{R\left({\stackrel{⌣}{x}}_i,x\right),i=1,2,\dots ,\stackrel{⌣}{n}\right\}.$$

The computational cost will be reduced significantly to $O\left(n{\stackrel{⌣}{n}}^2\right)$if $\stackrel{⌣}{n}$is much smaller than $n$. Furthermore, it can be proven that the minimizer of (2), $\stackrel{⌣}{\eta }$, by uniform sampling basis selection, has the same asymptotic convergence rate as the full basis minimizer $\widehat{\eta }$.

Although the uniform basis selection reduces the computational cost and the corresponding $\stackrel{⌣}{\eta }$achieves the optimal asymptotic convergence rate, it may fail to retain the data features occasionally. For example, when the data are not evenly distributed, it is hard for uniform sampling to capture the feature where there are only a few data points. In [14], an adaptive basis selection method is proposed. The main idea is to sample more basis functions where the response functions change largely and fewer basis functions on those flat regions. More details of adaptive basis selection method are shown in the following procedure:

Step 1Divide the range of responses ${\left\{y_i\right\}}_{i=1}^n$into $K$disjoint intervals, $S_1,\dots ,S_K$. Denote $\mid S_k\mid$as the number of observations in $S_k$.

Step 2For each $S_k$, draw a random sample of size $n_k$from this collection. Let $x^{\ast \left(k\right)}=\left(x_1^{\ast \left(k\right)},\dots ,x_{n_k}^{\ast \left(k\right)}\right)$be the predictor values.

Step 3Combine $x^{\ast \left(1\right)},\dots ,x^{\ast \left(K\right)}$together to form a set of sampled predictor values $\left\{x_1^{\ast },\dots ,x_{n\ast }^{\ast }\right\}$, where $n^{\ast }={\sum }_{k=1}^K{n}_k$.

Step 4Define

$${\mathcal{H}}_E={\mathcal{H}}_0\oplus \text{span}\left\{R\left(x_i^{\ast },\cdot \right),i=1,2,\dots ,n^{\ast }\right\}$$

as the effective model space.

By adaptive basis selection, the minimizer of (2) keeps the same form as that in Theorem 2.3:

$${\eta }_A\left(x\right)=\sum _{k=1}^M{d}_k{\xi }_k\left(x\right)+\sum _{i=1}^{n^{\ast }}{c}_{i}R\left(x_i^{\ast },x\right).$$

Let $R_{\ast }$be an $n\times n^{\ast }$matrix, and its $\left(i,j\right)$th entry is $R\left(x_i,x_j^{\ast }\right)$. Let $R_{\ast \ast }$be an $n^{\ast }\times n^{\ast }$matrix, and its $\left(i,j\right)$th entry is $R\left(x_i^{\ast },x_j^{\ast }\right)$. Then, the estimator ${\mathit{\eta }}_{\mathit{A}}$satisfies

$${\mathit{\eta }}_{\mathit{A}}=S{\mathbf{d}}_{\mathbf{A}}+R_{\ast }{\mathbf{c}}_{\mathbf{A}},$$

where ${\mathit{\eta }}_{\mathit{A}}={\left({\eta }_A\left(x_1\right),\cdots ,{\eta }_A\left(x_{n^{\ast }}\right)\right)}^T$, ${\mathbf{d}}_{\mathbf{A}}={\left(d_1,\dots ,d_M\right)}^T$, and ${\mathbf{c}}_{\mathbf{A}}={\left(c_1,\dots ,c_{n^{\ast }}\right)}^T$. Similar to (17), the linear system of equations in this case is

$$\left(\begin{array}{ll}{S}^{T}S& S^T{R}_{\ast }\\ R_{\ast }^{T}S& R_{\ast }^T{R}_{\ast }+\mathit{n\lambda }{R}_{\ast \ast }\end{array}\right)\left(\begin{array}{l}{\mathbf{d}}_{\mathbf{A}}\\ {\mathbf{c}}_{\mathbf{A}}\end{array}\right)=\left(\begin{array}{l}{S}^T\mathbf{y}\\ R_{\ast }^T\mathbf{y}\end{array}\right).$$E18

The computational complexity of solving (18) is of the order $O\left({\mathit{nn}}^{\ast 2}\right)$, so the method decreases the computational cost significantly. It can also be shown that the adaptive sampling basis selection smoothing spline estimator ${\eta }_A$has the same convergence property as the full basis method. More details about the consistency theory can be found in [14]. Moreover, adaptive sampling basis selection method for exponential family smoothing spline models was developed in [32].

3.2. Rounding algorithm

Other than sampling a smaller set of basis functions to save the computational resources, for example, the adaptive basis selection method presented previously, [15] proposed a new rounding algorithm to fit SSANOVA models in the context of big data.

Rounding algorithm: The details of rounding algorithm can be shown in the following procedure:

Step 1Assume that all predictors are continuous.

Step 2Convert all predictors to the interval $\left[0,1\right]$.

Step 3Round the raw data by using the transformation:

$$z_{i⟨j⟩}=\mathit{RD}\left(x_{i⟨j⟩}/r_{⟨j⟩}\right)r_{⟨j⟩},\text{for}i\in \left\{1,\cdots ,n\right\},j\in \left\{1,\cdots ,d\right\},$$

where the rounding parameter $r_{⟨j⟩}\in \left(0,1\right]$and rounding function $\mathit{RD}\left(\cdot \right)$transform input data to the nearest integer.

Step 4After replacing $x_{i⟨j⟩}$with $z_{i⟨j⟩}$, we redefine $S$and $R$in (16) and then estimate $\eta$by minimizing the penalized least squares (16).Remark 1

In Step 3, if $r_{⟨j⟩}$is the rounding parameter for $j$th predictor and its value is 0.03, then each $z_{i⟨j⟩}$is formed by rounding the corresponding $x_{i⟨j⟩}$to its nearest 0.03.

Remark 2

It is evident that the value of rounding parameter can influence the precision of approximation. The smaller the rounding parameter, the better the model estimation and the higher the computational cost.

Computational benefits: We now briefly explain why the implementation of rounding algorithm can reduce the computational loads. For example, if the rounding parameter $r=0.01$, it is obvious that $u\le 101$, where $u$denotes the number of uniquely observed values. In conclusion, using user-tunable rounding algorithm can dramatically reduce the computational burden of fitting SSANOVA models from the order of $O\left(n^3\right)$to $O\left(u^3\right)$, where $u\ll n$.

Case study: To illustrate the benefit of the rounding algorithm, we apply the algorithm to the electroencephalography (EEG) dataset. Note that EEG is a monitoring method to record the electrical activity of the brain. It can be used to diagnose sleep disorders, epilepsy, encephalopathies, and brain death.

The dataset [33] contains 44 controls and 76 alcoholics. Each subject was repeatedly measured 10 times by using visual stimulus at a frequency of 256 Hz. This brings about $n=10$replications $\times 120$subjects $\times 256$time points $=\mathrm{307,200}$observations. There are two predictors, time and group (control vs. alcoholic). We apply the cubic spline to the time effect and the nominal spline to the group effect.

After applying the model to the unrounded data, rounded data with rounding parameter $r=0.01$and $r=0.05$for time covariate, we can obtain a summary table about GCV, AIC [34], BIC [35], and running time in Table 2.

	GCV	AIC	BIC	CPU time (seconds)
Unrounded data	85.9574	2,240,019	2,240,562	15.65
Rounded data with $r=0.01$	86.6667	2,242,544	2,242,833	1.22
Rounded data with $r=0.05$	86.7654	2,242,893	2,243,089	1.13

Table 2.

Fit statistics and running time for SSANOVA models.

Based on Table 2, we can easily see that there are no significant difference among the GCV scores and AIC/BIC. In addition using rounding algorithm reduces $92\%$CPU time compared to using unrounded dataset.

4. Conclusion

Smoothing spline ANOVA (SSANOVA) models are widely used in applications [11, 20, 36, 37]. In this chapter, we introduced the general framework of the SSANOVA models in Section 2. In Section 3, we discussed the models under the big data settings. When the volume of data grows, fitting the models is computing-intensive [11]. The adaptive basis selection algorithm [14] and rounding algorithm [15] we presented can significantly reduce the computational cost.

Acknowledgments

This work is partially supported by the NIH grants R01 GM122080 and R01 GM113242; NSF grants DMS-1222718, DMS-1438957, and DMS-1228288; and NSFC grant 71331005.

Conflict of interest

The authors whose names are listed immediately below certify that they have NO affiliations with or involvement in any organization or entity with any financial interest (such as honoraria; educational grants; participation in speakers’ bureaus; membership, employment, consultancies, stock ownership, or other equity interests; and expert testimony or patent-licensing arrangements) or nonfinancial interest (such as personal or professional relationships, affiliations, knowledge, or beliefs) in the subject matter or materials discussed in this manuscript.

In this appendix, we use two examples to illustrate how to implement smoothing spline ANOVA (SSANOVA) models in R. The gss package in R, which can be downloaded on the CRAN https://cran.r-project.org/, is utilized.

We now load the gss package:

library(gss)

Example I: Apply the smoothing spline to a simulated dataset.

Suppose that the predictor $x$follows a uniform distribution on $\left[0,1\right]$, and the response $y$is generated based on $y=5+2\text{cos}\left(3\mathit{\pi x}\right)+\mathrm{ϵ}$, where $\mathrm{ϵ}\sim N\left(0,1\right)$.

x<-runif(100);y<-5+2*cos(3*pi*x)+rnorm(x)

Then, fit cubic smoothing spline model:

cubic.fit<-ssanova(y^˜x)

To evaluate the predicted values, one uses:

new<-data.frame(x=seq(min(x),max(x),len=50))

est<-predict(cubic.fit,new,se=TRUE)

These.fitparameter indicates if one can get the pointwise standard errors for the predicted values. The predicted values and Bayesian confidence interval, shown in Figure 4, are generated by:

Figure 4.

The solid red line represents the fitted values. The green lines represent the 95% Bayesian confidence interval. The raw data are shown as the circles.

plot(x,y,col=1)

lines(new$x,est$fit,col=2)

lines(new$x,est$fit+1.96*est$se,col=3)

lines(new$x,est$fit-1.96*est$se,col=3)

Example II:Apply the SSANOVA model to a real dataset.

In this example, we illustrate how to implement the SSANOVA model using the gss package. The data is from an experiment in which a single-cylinder engine is run with ethanol to see how the nox concentrationnoxin the exhaust depends on the compression ratiocompand the equivalence ratioequi. The fitted model contains two predictors (compandequi) and one interaction term.

data(nox)

nox.fit <- ssanova(log10(nox)^˜comp*equi,data=nox)

The predicted values are shown in Figure 5.

Figure 5.

The x-axis, y-axis, and z-axis represent the compression ratio, the equivalence ratio, and the predicted values, respectively.

x=seq(min(nox$comp),max(nox$comp),len=50)

y=seq(min(nox$equi),max(nox$equi),len=50)

temp <- function(x, y){

new=data.frame(comp=x,equi=y)

return(predict(nox.fit,new,se=FALSE))

}

z=outer(x, y, temp)

persp(x,y,z,theta = 30).