Optimal Penalized Function-on-Function Regression under a Reproducing Kernel Hilbert Space Framework

TL;DR

This paper introduces a penalized function-on-function regression model within a Reproducing Kernel Hilbert Space framework, providing a data-adaptive finite-dimensional solution with proven minimax convergence and demonstrated advantages through simulations and real data applications.

Contribution

It develops a novel penalized regression estimator for function-on-function data with a Representer Theorem, enabling efficient optimization and theoretical guarantees.

Findings

Estimator achieves minimax convergence rate in mean prediction.

Simulation studies show numerical advantages over existing methods.

Application to real data demonstrates practical utility.

Abstract

Many scientific studies collect data where the response and predictor variables are both functions of time, location, or some other covariate. Understanding the relationship between these functional variables is a common goal in these studies. Motivated from two real-life examples, we present in this paper a function-on-function regression model that can be used to analyze such kind of functional data. Our estimator of the 2D coefficient function is the optimizer of a form of penalized least squares where the penalty enforces a certain level of smoothness on the estimator. Our first result is the Representer Theorem which states that the exact optimizer of the penalized least squares actually resides in a data-adaptive finite dimensional subspace although the optimization problem is defined on a function space of infinite dimensions. This theorem then allows us an easy incorporation of the Gaussian quadrature into the optimization of the penalized least squares, which can be carried out through standard numerical procedures. We also show that our estimator achieves the minimax convergence rate in mean prediction under the framework of function-on-function regression. Extensive simulation studies demonstrate the numerical advantages of our method over the existing ones, where a sparse functional data extension is also introduced. The proposed method is then applied to our motivating examples of the benchmark Canadian weather data and a histone regulation study.

Figures

Tables

Table 1. Table 1: The mean, standard deviation and three quartiles of ISEs for the three approaches. The best result on each metric is in boldface.

Method	Median	Mean	Standard Deviation	1st Qu.	3rd Qu.
OPFFR	21.6400	40.2800	45.7631	13.8000	36.1700
FDA	25.9000	44.1600	56.9544	18.7400	40.6100
FPCA	30.7752	45.5065	45.7763	20.5031	52.1827

Table 2. Table 2: The mean, standard deviation and three quartiles of ISEs for the four approaches. The best result on each metric is in boldface.

Method	Median	Mean	Standard Deviation	1st Qu.	3rd Qu.
OPFFR	1.5700	7.7120	18.9180	0.5077	5.1900
OPFFR-S	1.4070	7.7150	18.6037	0.6972	5.5820
FDA	2.2060	7.9770	18.7004	0.5461	6.2750
FPCA	2.0170	8.4720	18.3978	0.9126	6.1790

Equations

$$Y_{i}(t)=\alpha(t)+\int_{I_{x}}\beta\big{(}t,s\big{)}X_{i}(s)ds+\epsilon_{i}(t),~{}~{}~{}~{}t\in I_{y},$$

$$Y_i(t)-\mathbb{E}{Y}_i(t)={\int }_{I_x}\beta (t,s)\{X_i(s)-\mathbb{E}{X}_i(s)\}ds+{ϵ}_i(t),t\in I_y,$$

$$Y_{i}(t)=\int_{I_{x}}\beta\big{(}t,s\big{)}X_{i}(s)ds+\epsilon_{i}(t),~{}~{}~{}~{}t\in I_{y}.$$

$${1\over n}\sum_{i=1}^{n}\int_{I_{y}}\Big{\{}Y_{i}(t)-\int_{I_{x}}\beta(t,s)X_{i}(s)ds\Big{\}}^{2}dt+\lambda J(\beta)$$

$$K_0((t_1,s_1),(t_2,s_2))$$

$$K_1((t_1,s_1),(t_2,s_2))$$

$$\beta (t,s)$$

$$\displaystyle=\Big{\{}d_{\beta_{y}}^{\top}\psi_{y}(t)+c_{\beta_{y}}^{\top}(K_{1y}Y)(t)\Big{\}}\Big{\{}d_{\beta_{x}}^{\top}\psi_{x}(s)+c_{\beta_{x}}^{\top}(K_{1x}X)(s)\Big{\}},$$

$$\langle f,g\rangle_{\mathcal{H}_{y}}=\Big{(}\int_{0}^{1}f\int_{0}^{1}g+\int_{0}^{1}f^{\prime}\int_{0}^{1}g^{\prime}\Big{)}+\int_{0}^{1}f^{\prime\prime}g^{\prime\prime}dt.$$

$$J(\beta)=\int_{0}^{1}\Big{[}\Big{\{}\int_{0}^{1}\frac{\partial^{2}}{\partial s^{2}}\beta(t,s)dt\Big{\}}^{2}+\Big{\{}\int_{0}^{1}\frac{\partial^{3}}{\partial t\partial s^{2}}\beta(t,s)dt\Big{\}}^{2}\Big{]}ds\\ +\int_{0}^{1}\Big{[}\Big{\{}\int_{0}^{1}\frac{\partial^{2}}{\partial t^{2}}\beta(t,s)ds\Big{\}}^{2}+\Big{\{}\int_{0}^{1}\frac{\partial^{3}}{\partial t^{2}\partial s}\beta(t,s)ds\Big{\}}^{2}\Big{]}dt+\int_{0}^{1}\int_{0}^{1}\Big{\{}\frac{\partial^{4}}{\partial t^{2}\partial s^{2}}\beta(t,s)\Big{\}}^{2}dtds$$

$$(Y_w-S_w\mathbf{d}-R_w\mathbf{c}{)}^T(Y_w-S_w\mathbf{d}-R_w\mathbf{c})+n\lambda {\mathbf{c}}^{T}Q\mathbf{c}$$

$$V(\lambda )=\frac{(nT{)}^{-1}{Y}_w^T(I-A(\lambda ){)}^2{Y}_w}{\{(nT{)}^{-1}tr(I-\alpha A(\lambda )){\}}^2}$$

$${\mathfrak{R}}_n({\hat{\beta }}_n)$$

$$=$$

$$\beta(t,s)=\Big{\langle}K\big{(}(t,s),\cdot\big{)},\beta\Big{\rangle}_{{\cal H}(K)},~{}~{}~{}\mbox{for any }\beta\in{\cal H}(K).$$

$$K(\beta)(\cdot)=\Big{\langle}K\big{(}(t,s),\cdot\big{)},\beta\Big{\rangle}_{L_{2}}=\int\int_{I}K((t,s),\cdot)\beta(t,s)dtds.$$

$$K((t_1,s_1),(t_2,s_2))=\sum _{k=1}^{\infty }{\kappa }_k{\zeta }_k(t_1,s_1){\zeta }_k(t_2,s_2),K({\zeta }_k)={\kappa }_k{\zeta }_k,k=1,2,\dots .$$

$$\Pi\big{(}(t_{1},s_{1}),(t_{2},s_{2})\big{)}=\int\int\int_{I_{x}\times I_{x}\times I_{y}}K^{1/2}\big{(}(t_{1},s_{1}),(z,u)\big{)}~{}C(u,v)~{}K^{1/2}\big{(}(t_{2},s_{2}),(z,v)\big{)}dudvdz.$$

$$\Pi\big{(}(t_{1},s_{1}),(t_{2},s_{2})\big{)}=\sum_{k=1}^{\infty}\rho_{k}\phi_{k}(t_{1},s_{1})\phi_{k}(t_{2},s_{2}),~{}~{}~{}\forall~{}(t_{1},s_{1}),(t_{2},s_{2})\in I_{y}\times I_{x}.$$

$$\int\mathbb{E}\Big{(}\int\beta(t,s)X(s)dt\Big{)}^{4}dt\leq c\int\Big{(}\mathbb{E}\Big{(}\int\beta(t,s)X(s)ds\Big{)}^{2}\Big{)}^{2}dt$$

$$\lim_{A\rightarrow\infty}\lim_{n\rightarrow\infty}\sup_{\beta\in{\cal H}(K)}\mathbb{P}\Big{\{}\mathfrak{R}_{n}\geq An^{-{2r\over 2r+1}}\Big{\}}=0,$$

$$\lim_{c\rightarrow 0}\lim_{n\rightarrow\infty}\inf_{\tilde{\eta}}\sup_{\beta\in{\cal H}(K)}\mathbb{P}\Big{(}{\mathfrak{R}}_{n}\geq cn^{-{2r\over 2r+1}}\Big{)}=1,$$

Full text

Keywords: Function-on-Function regression; Representer Theorem; Reproducing kernel Hilbert space; Penalized least squares; Minimax convergence rate.

1 Introduction

Functional data have attracted much attention in the past decades (Ramsay & Silverman, 2005). Most of the existing literature has only considered the regression models of a scalar response against one or more functional predictors, possibly with some scalar predictors as well. Some of them considered a reproducing kernel Hilbert space framework. For example, Yuan & Cai (2010) provided a thorough theoretical analysis of the penalized functional linear regression model with a scalar response. The paper laid the foundation for several theoretical developments including the representer theorem and minimax convergence rates for prediction and estimation for penalized functional linear regression models. In a follow-up, Cai & Yuan (2012) showed that the minimax rate of convergence for the excess prediction risk is determined by both the covariance kernel and the reproducing kernel. Then they designed a data-driven roughness regularization predictor that can achieve the optimal convergence rate adaptively without the knowledge of the covariance kernel. Du & Wang (2014) extended the work of Yuan & Cai (2010) to the setting of a generalized functional linear model, where the scalar response comes from an exponential family distribution.

In contrast to these functional linear regression models with a scalar response, the model with a functional response $Y(t)$ over a functional predictor $X(s)$ has only been scarcely investigated (Yao et al., 2005b; Ramsay & Silverman, 2005). Such data with functional responses and predictors are abundant in practice. We shall now present two motivating examples.

Example 1.1

**Canadian Weather Data

Daily temperature and precipitation at 35 different locations in Canada averaged over 1960 to 1994 were collected (Figure 1). The main interest is to use the daily temperature profile to predict the daily precipitation profile for a location in Canada.**

Example 1.2

**Histone Regulation Data

Extensive researches have been shown that histone variants, i.e. histones with structural changes compared to their primary sequence, play an important role in the regulation of chromatin metabolism and gene activity (Ausió, 2006). An ultra-high throughput time course experiment was conducted to study the regulation mechanism during heat stress in Arabidopsis thaliana. The genome-wide histone variant distribution was measured by ChIP sequencing (ChIP-seq) (Johnson et al., 2007) experiments. We computed histone levels over 350 base pairs (bp) on genomes from the ChIP-seq data, see left panel in Figure 2. The RNA sequencing (RNA-seq) (Wang et al., 2009) experiments measured the expression levels over seven time points within 24 hours, see right panel in Figure 2. Of primary interest is to study the regulation mechanism between gene expression levels over time domain and histone levels over spatial domain.**

Motivated by the examples, we now present the statistical model. Let \big{\{}(X(s),Y(t)):s\in I_{x},t\in I_{y}\big{\}} be two random processes defined respectively on $I_{x},I_{y}\subseteq\mathbb{R}$. Suppose $n$ independent copies of \big{(}X,Y\big{)} are observed: \big{(}X_{i}(s),Y_{i}(t)\big{)}, $i=1,\ldots,n$. The functional linear regression model of interest is

[TABLE]

where $\alpha(\cdot):I_{y}\rightarrow\mathbb{R}$ is the intercept function, $\beta(\cdot,~{}\cdot):I_{y}\times I_{x}\rightarrow\mathbb{R}$ is a bivariate coefficient function, and $\epsilon_{i}(t)$, independent of $X_{i}(s)$, are i.i.d. random error functions with $\mathbb{E}\epsilon_{i}(t)=0$ and $\mathbb{E}\|\epsilon_{i}(t)\|_{2}^{2}<\infty$. Here $\|\cdot\|_{2}$ denotes the $L_{2}$-norm. In Example 1.1, $Y_{i}(t)$ and $X_{i}(t)$ represent the daily precipitation and temperature at station $i$. In Example 1.2, the expression levels of gene $i$ over seven time points, $Y_{i}(t)$, from RNA-seq is used as the functional response. The histone levels of gene $i$ over 350 base pairs (bp), $X_{i}(s)$, from ChIP-seq is used as the functional predictor.

At a first look, model (1) might give the (wrong) impression of being an easy extension from the model with a scalar response, with the latter obtained from (1) by removing all the $t$ notation. However, the coefficient function in the scalar response case is univariate and thus can be easily estimated by most off-the-shelf smoothing methods. When extended to estimating a bivariate coefficient function $\beta(t,s)$ in (1), many of these smoothing methods may encounter major numerical and/or theoretical difficulties. This partly explains the much less abundance of research in this direction.

Some exceptions though are reviewed below. Cuevas et al. (2002) considered a fixed design case, a different setting from (1) with $Y_{i}(t)$ and $X_{i}(s)$ represented and analyzed as sequences. Nonetheless they provided many motivating applications in neuroscience, signal transmission, pharmacology, and chemometrics, where (1) can apply. The historical functional linear model in Malfait & Ramsay (2003) was among the first to study regression of a response functional variable over a predictor functional variable, or more precisely, the history of the predictor function. Ferraty et al. (2011) proposed a simple extension of the classical Nadaraya-Watson estimator to the functional case and derived its convergence rates. They provided no numerical results on the empirical performance of their kernel estimator. Benatia et al. (2015) extended ridge regression to the functional setting. However, their estimation relied on an empirical estimate of the covariance process of predictor functions. Theoretically sound as it is, this covariance process estimate is generally not reliable in practice. Consequently, their coefficient surface estimates suffered as shown in their simulation plots. Meyer et al. (2015) proposed a Bayesian function-on-function regression model for multi-level functional data, where the basis expansions of functional parameters were regularized by basis-space prior distributions and a random effect function was introduced to incorporate the with-subject correlation between functional observations.

A popular approach has been the functional principal component analysis (FPCA) as in Yao et al. (2005b) and Crambes & Mas (2013). The approach starts with a basis representation of $\beta(t,s)$ in terms of the eigenfunctions in the Karhunen-Loève expansions of $Y(t)$ and $X(s)$. Since this representation has infinitely many terms, it is truncated at certain point to obtain an estimable basis expansion of $\beta(t,s)$. Yao et al. (2005b) studied a general data setting where $Y(t)$ and $X(s)$ are only sparsely observed at some random points. They derived the consistency and proposed asymptotic point-wise confidence bands for predicting response trajectories. Crambes & Mas (2013) furthered the theoretical investigation of the FPCA approach by providing a minimax optimal rates in terms of the mean square prediction error. However, the FPCA approach has a couple of critical drawbacks. Firstly, $\beta(t,s)$ is a statistical quantity unrelated to $Y(t)$ or $X(s)$. Hence the leading eigenfunctions in the truncated Karhunen-Loève expansions of $Y(t)$ and $X(s)$ may not be an effective basis for representing $\beta(t,s)$. See, e.g., Cai & Yuan (2012) and Du & Wang (2014) for some scalar-response examples where the FPCA approach breaks down when the aforementioned situation happens. Secondly, the truncation point is integer-valued and thus only has a discrete control on the model complexity. This puts it at disadvantage against the roughness penalty regularization approach, which offers a continuous control via a positive and real-valued smoothing parameter (Ramsay & Silverman, 2005, Chapter 5).

In this paper, we consider a penalized function-on-function regression approach to estimating the bivariate coefficient function $\beta(t,s)$. There have been a few recent developments in the direction of penalized function-on-function regression. Lian (2015) studied the convergence rates of the function-on-function regression model under a reproducing kernel Hilbert space framework. Although his model resembled model (1), he developed everything with the variable $t$ fixed and did not enforce any regularization on the $t$ direction. Firstly, this lack of $t$-regularization can be problematic since this leaves the noisy errors on the $t$ direction completely uncontrolled and can result in an $\beta(s,t)$ estimate that is very rough on the $t$ direction. Secondly, this simplification of fixing $t$ essentially reduces the problem to a functional linear model with a scalar response and thus makes all the results in Yuan & Cai (2010) directly transferrable even without calling on any new proofs. The R package fda maintained by Ramsay et al. has implemented a version of penalized B-spline estimation of $\beta(t,s)$ with a fixed smoothing parameter. Ivanescu et al. (2015) considered a penalized function-on-function regression model where the coefficient functions were represented by expansions into some basis system such as tensor cubic B-splines. Quadratic penalties on the expansion coefficients were used to control the smoothness of the estimates. This work provided a nice multiple-predictor-function extension to the function-on-function regression model in the fda package. Scheipl & Greven (2016) studied the identifiability issue in these penalized function-on-function regression models. However, this penalized B-spline approach has several well-known drawbacks. First, it is difficult to show any theoretical optimality such as the minimax risk of mean prediction in Cai & Yuan (2012). So its theoretical soundness is hard to justify. Moreover, the B-spline expansion is only an approximate solution to the optimization of the penalized least squares score. Hence the penalized B-spline estimate is not numerically optimal from the beginning either. These drawbacks can have negative impacts on the numerical performance as we shall see from the simulation results in Section 4.

The penalized function-on-function regression method proposed in this paper obtains its estimator of $\beta(t,s)$ through the minimization of penalized least squares on a reproducing kernel Hilbert space that is naturally associated with the roughness penalty. Such a natural formulation through a reproducing kernel Hilbert space offers several advantages comparing with the existing penalized function-on-function regression methods. Firstly, it allows us to establish a Representer Theorem which states that, although the optimization of the penalized least squares is defined on an infinite dimensional function space, its solution actually resides in a data-adaptive finite dimensional subspace. This result guarantees an exact solution when the optimization is carried out on this finite dimensional subspace. This result itself is a nontrivial generalization of the Representer Theorems in the scenarios of nonparametric smooth regression model (Wahba, 1990) and the penalized functional regression model with a scalar response (Yuan & Cai, 2010). Based on the Representer Theorem, we propose an estimation algorithm which uses penalized least squares and Gaussian quadrature with the Gauss-Legendre rule to estimate the bivariate coefficient function. The smoothing parameter is selected by the generalized cross validation (GCV) method. Secondly, the reproducing kernel Hilbert space framework allows us to show that our estimator has the optimal rate of mean prediction since it achieves the minimax convergence rate in terms of the excess risk. This generalizes the results in Cai & Yuan (2012) and Du & Wang (2014) for functional linear regression with a scalar response to the functional response scenario. In the numerical study, we have also considered the problem with sparsely sampled data. Particularly, we introduce an extra pre-smoothing step before applying the proposed penalized functional regression model. The pre-smoothing step implements the principal-component-analysis-through-expectation (PACE) method in Yao et al. (2005a). Our extensive simulation studies demonstrate the numerical advantages of our method over the existing ones. In summary, our method has the following distinguishing features: (i) it makes no structural dependence assumptions of $\beta(t,s)$ over the predictor and response processes; (ii) the Representer Theorem guarantees an exact solution instead of an approximation to the optimization of the penalized score; (iii) benefited from the Representer Theorem, we develop a numerically reliable algorithm that has sound performance in simulations; (iv) we show theoretically the estimator achieves the optimal minimax convergence rate in mean prediction.

The rest of the paper is organized as follows. In Section 2, we first derive a Representer Theorem showing that the solution of the minimization of penalized least squares can be found in a finite-dimension subspace. In addition, an easily implementable estimation algorithm is considered in the Section 2. In Section 3, we prove that our method has the optimal rate of mean prediction. Numerical experiments are reported in Section 4, where we compare our method with the functional linear regressions in (Ramsay & Silverman, 2005; Yao et al., 2005b) in terms of prediction accuracy. Two real data examples, the Canadian weather data, and the histone regulation data are analyzed in Section 5. Discussion in Section 6 concludes the paper. Proofs of the theorems are collected in Supplementary Material.

2 Penalized Functional Linear Regression Method

We first introduce a simplification to model (1). Since model (1) implies that

[TABLE]

we may, for simplicity, only consider $X$ and $Y$ to be centered, i.e., $\mathbb{E}X=\mathbb{E}Y=0$. Thus, the functional linear regression model takes the form of

[TABLE]

2.1 The Representer Theorem

Assume that the unknown $\beta$ resides in a reproducing kernel Hilbert space ${\cal H}(K)$ with the reproducing kernel $K:I\times I\rightarrow\mathbb{R}$, where $I=I_{y}\times I_{x}$. The estimate $\hat{\beta}_{n}$ can be obtained by minimizing the following penalized least squares functional

[TABLE]

with respect to $\beta\in{\cal H}(K)$, where the sum of integrated squared errors represents the goodness-of-fit, $J$ is a roughness penalty on $\beta$, and $\lambda>0$ is the smoothing parameter balancing the trade-off. We now establish a Representer Theorem stating that $\hat{\beta}_{n}$ actually resides in a finite dimensional subspace of $\mathcal{H}(K)$. This result generalizes Theorem 1 in Yuan & Cai (2010) and facilitates the computation by reducing an infinite dimensional optimization problem to a finite dimensional one.

Note that the penalty functional $J$ is a squared semi-norm on ${\cal H}(K)$. Its null space ${\cal H}_{0}=\{\beta\in{\cal H}(K):J(\beta)=0\}$ is a finite-dimensional linear subspace of ${\cal H}(K)$. Denote by ${\cal H}_{1}$ its orthogonal complement in ${\cal H}(K)$ such that ${\cal H}(K)={\cal H}_{0}\oplus{\cal H}_{1}$. For any $\beta\in{\cal H}(K)$, there exists a unique decomposition $\beta=\beta_{0}+\beta_{1}$ where $\beta_{0}\in{\cal H}_{0}$ and $\beta_{1}\in{\cal H}_{1}$. Let $K_{0}(\cdot,\cdot)$ and $K_{1}(\cdot,\cdot)$ be the corresponding reproducing kernels of ${\cal H}_{0}$ and ${\cal H}_{1}$. Then $K_{0}$ and $K_{1}$ are both nonnegative definite operators on $L_{2}$, and $K=K_{0}+K_{1}$. In fact the penalty term $J(\beta)=\|\beta\|_{K_{1}}^{2}=\|\beta_{1}\|_{K_{1}}^{2}$. By the theory of reproducing kernel Hilbert spaces, $\mathcal{H}(K)$ has a tensor product decomposition $\mathcal{H}(K)=\mathcal{H}_{y}(K_{y})\otimes\mathcal{H}_{x}(K_{x})$. Here $\mathcal{H}_{y}(K_{y})$ is the reproducing kernel Hilbert space with a reproducing kernel $K_{y}:I_{y}\times I_{y}\rightarrow\mathbb{R}$, and $\mathcal{H}_{x}(K_{x})$ is the reproducing kernel Hilbert space with a reproducing kernel $K_{x}:I_{x}\times I_{x}\rightarrow\mathbb{R}$. For the reproducing kernels, we have $K(t,s)=K_{y}(t)K_{x}(s)$. Note that the functions in $\mathcal{H}_{y}(K_{y})$ and $\mathcal{H}_{x}(K_{x})$ are univariate and defined respectively on $I_{y}$ and $I_{x}$. Similar to the decomposition of $\mathcal{H}$ and $K$, we have the tensor sum decompositions of the marginal subspaces $\mathcal{H}_{y}(K_{y})=\mathcal{H}_{0y}\oplus\mathcal{H}_{1y}$ and $\mathcal{H}_{x}(K_{x})=\mathcal{H}_{0x}\oplus\mathcal{H}_{1x}$, and the orthogonal decompositions of the marginal reproducing kernels $K_{y}=K_{0y}+K_{1y}$ and $K_{x}=K_{0x}+K_{1x}$. Here $K_{*}$ is a reproducing kernel on $\mathcal{H}_{*}$ with $*$ running through the index set $\{0y,1y,0x,1x\}$.

Upon piecing the marginal decomposition parts back to the tensor product space, we obtain $\mathcal{H}_{0}=\mathcal{H}_{0y}\otimes\mathcal{H}_{0x}$ and $\mathcal{H}_{1}=(\mathcal{H}_{0y}\otimes\mathcal{H}_{1x})\oplus(\mathcal{H}_{1y}\otimes\mathcal{H}_{0x})\oplus(\mathcal{H}_{1y}\otimes\mathcal{H}_{1x})$. Correspondingly, the reproducing kernels satisfy that

[TABLE]

Let $N_{y}=\text{dim}(\mathcal{H}_{0y})$ and $N_{x}=\text{dim}(\mathcal{H}_{0x})$. Denote by $\{\psi_{k,y}:k=1,\ldots,N_{y}\}$ and $\{\psi_{l,x}:l=1,\ldots,N_{x}\}$ respectively the basis functions of $\mathcal{H}_{0y}$ and $\mathcal{H}_{0x}$. With some abuse of notation, define $(K_{1y}g)(\cdot)=\int_{I_{y}}K_{1y}(\cdot,t)g(t)dt$ and $(K_{1x}f)(\cdot)=\int_{I_{x}}K_{1x}(\cdot,s)f(s)ds$. Now we can state the Representer Theorem as follows with its proof collected in the Supplementary Material.

Theorem 2.1

Let $\hat{\beta}_{n}$ be the minimizer of (3) in $\mathcal{H}(K)$. Then $\hat{\beta}_{n}$ resides in the subspace of functions of the form

[TABLE]

where $d_{\beta_{y}}=(d_{1,\beta_{y}},\ldots,d_{N_{y},\beta_{y}})^{\top}$, $c_{\beta_{y}}=(c_{1,\beta_{y}},\ldots,c_{n,\beta_{y}})^{\top}$, $d_{\beta_{x}}=(d_{1,\beta_{x}},\ldots,d_{N_{x},\beta_{x}})^{\top}$ and $c_{\beta_{x}}=(c_{1,\beta_{x}},\ldots,c_{n,\beta_{x}})^{\top}$ are some coefficient vectors, and $\psi_{x},\psi_{y},K_{1y}Y$ and $K_{1x}X$ are vectors of functions.

For the purpose of illustration, we give a detailed example below.

Example 2.1

Consider the case of tensor product cubic splines with $I_{y}=I_{x}=[0,1]$. The marginal spaces $\mathcal{H}_{y}(K_{y})=\mathcal{H}_{x}(K_{x})=\{g:\int_{0}^{1}(g^{\prime\prime})^{2}<\infty\}$ with the inner product

[TABLE]

The marginal space $\mathcal{H}_{y}(K_{y})$ can be further decomposed into the tensor sum of $\mathcal{H}_{0y}=\{g:g^{\prime\prime}=0\}$ and $\mathcal{H}_{1y}=\{g:\int_{0}^{1}g=\int_{0}^{1}g^{\prime}=0,\int_{0}^{1}(g^{\prime\prime})^{2}<\infty\}$. The reproducing kernel $K_{y}$ is the orthogonal sum of $K_{0y}(t_{1},t_{2})=1+r_{1}(t_{1})r_{1}(t_{2})$ and $K_{1y}(t_{1},t_{2})=r_{2}(t_{1})r_{2}(t_{2})-r_{4}(|t_{1}-t_{2}|)$, where $r_{\nu}(t)=B_{\nu}(t)/\nu!$ is a scaled version of the Bernoulli polynomial $B_{\nu}$. The space $\mathcal{H}_{0y}$ has a dimension of $N_{y}=2$ and a set of basis functions $\{1,r_{1}(t)\}$.

The function space $\mathcal{H}(K)$ is defined as $\mathcal{H}(K)=\{\beta:J(\beta)<\infty\}$ with the reproducing kernel $K(t,s)=K_{y}(t)K_{x}(s)$ and the penalty functional

[TABLE]

We have $\mathcal{H}(K)=\mathcal{H}_{y}(K_{y})\otimes\mathcal{H}_{x}(K_{x})$ and $K=K_{y}K_{x}$; see, e.g., Chapter 2 of Gu (2013).

2.2 Estimation Algorithm

To introduce the computational algorithm, we first need some simplification of notation. Let $N=N_{y}N_{x}$ and $L=n(N_{y}+N_{x}+n)$. We rewrite the functions spanning the subspace in Theorem 2.1 as $\psi_{1}(t,s)=\psi_{1,y}(t)\psi_{1,x}(s)$, $\cdots$, $\psi_{N}(t,s)=\psi_{N_{y},y}(t)\psi_{N_{x},x}(s)$ and $\xi_{1}(t,s)=\psi_{1,y}(t)(K_{1x}X_{1})(s)$, $\cdots$, $\xi_{L}(t,s)=(K_{1y}Y_{n})(t)(K_{1x}X_{n})(s)$. Thus a function in this subspace has the form $\beta(t,s)=\mathbf{d}^{T}\psi(t,s)+\mathbf{c}^{T}\xi(t,s)$ for some coefficient vectors $\mathbf{d}$, $\mathbf{c}$ and vectors of functions $\psi(t,s)$, $\xi(t,s)$. To solve (3), we choose Gaussian quadrature with the Gauss-Legendre rule to calculate the integrals. Consider the Gaussian quadrature evaluation of an integral on $I_{y}$ with knots $\{t_{1},\cdots,t_{T}\}$ and weights $\{\alpha_{1},\cdots,\alpha_{T}\}$ such that $\int_{I_{y}}f(t)dt=\sum_{j=1}^{T}\alpha_{j}f(t_{j})$. Let $W$ be the diagonal matrix with $\alpha_{1},\cdots,\alpha_{T}$ repeating $n$ times on the diagonal. Then the estimation of $\beta$ in (3) reduces to the minimization of

[TABLE]

with respect to $\mathbf{d}$ and $\mathbf{c}$, where $Y_{w}=W^{1/2}Y$ with $Y=(Y_{1}(t_{1}),\ldots,Y_{1}(t_{T}),\ldots,Y_{n}(t_{1}),\ldots,Y_{n}(t_{T}))^{\top}$, $S_{w}=W^{1/2}S$ with $S$ being an $nT\times N$ matrix with the $((i-1)T+j,\nu)$th entry $\int_{I_{x}}\psi_{\nu}(t_{j},s)X_{i}(s)ds$, $R_{w}=W^{1/2}R$ with $R$ being an $nT\times L$ matrix with the $((i-1)T+j,k)$th entry $\int_{I_{x}}\xi_{k}(t_{j},s)X_{i}(s)ds$, and $Q$ is a $L\times L$ matrix with the $(i,j)$th entry $\langle\xi_{i},\xi_{j}\rangle_{\mathcal{H}_{1}}$. Let $Q_{x}=\left[\int_{0}^{1}\int_{0}^{1}X_{i}(u)K(u,v)X_{j}(v)dudv\right]_{i,j=1}^{n}$, $Q_{y}=\left[\int_{0}^{1}\int_{0}^{1}Y_{i}(u)K(u,v)Y_{j}(v)dudv\right]_{i,j=1}^{n}$, and $Q_{xy}=Q_{x}\otimes Q_{y}$, we have $Q=\text{diag}(Q_{x},Q_{x},Q_{y},Q_{y},Q_{xy})$.

We then utilize standard numerical linear algebra procudures such as the Cholesky decomposition with pivoting and forward and back substitutions, to calculate $\mathbf{c}$ and $\mathbf{d}$ in (5) (Gu, 2013, Section 3.5). To choose the smoothing parameter $\lambda$ in (5), a modified Generalized Cross-Validation (GCV) score (Wahba & Craven, 1979),

[TABLE]

is implemented, where $\alpha>1$ is a fudge factor curbing undersmoothing (Kim & Gu, 2004) and $A(\lambda)$ is the smoothing matrix bridging the prediction $\hat{Y}_{w}$ and the observation $Y_{w}$ as $\hat{Y}_{w}=A(\lambda)Y_{w}$, similar to the hat matrix in a general linear model.

3 Optimal Mean Prediction Risk

We are interested in the estimation of coefficient function $\beta$ and mean prediction, that is, to recover the functional $\eta_{\beta}(X,\cdot)=\int_{I_{x}}\beta(\cdot,s)X(s)ds$ based on the training sample $(X_{i},Y_{i})$, $i=1,\ldots,n$. Let $\hat{\beta}_{n}(t,s)$ be an estimate of $\beta(t,s)$. Suppose $(X_{n+1},Y_{n+1})$ is a new observation that has the same distribution as and is also independent of $(X_{i},Y_{i})$, $i=1,\ldots,n$. Then the prediction accuracy can be naturally measured by the excess risk

[TABLE]

where $\mathbb{E}^{*}$ represents the expectation taken over $(X_{n+1},Y_{n+1})$ only. We shall study the convergence rate of ${\mathfrak{R}}_{n}$ as the sample size $n$ increases.

This section collects two theorems whose combination indicates that our estimator achieves the optimal minimax convergence rate in mean prediction. We first establish the minimax lower bound for the convergence rate of the excess risk ${\mathfrak{R}}_{n}$. There is a one-to-one relationship between $K$ and ${\cal H}(K)$ which is a linear functional space endowed with an inner product $\langle\cdot,\cdot\rangle_{{\cal H}(K)}$ such that

[TABLE]

The kernel $K$ can also be treated as an integral operator such that

[TABLE]

It follows from the spectral theorem that there exist a set of orthonormal eigenfunctions $\{\zeta_{k}:k\geq 1\}$ and a sequence of eigenvalues $\kappa_{1}\geq\kappa_{2}\geq\cdots>0$ such that

[TABLE]

Denote K^{1/2}\big{(}(t_{1},s_{1}),(t_{2},s_{2}))\big{)}=\sum_{k=1}^{\infty}\kappa_{k}^{1/2}\zeta_{k}(t_{1},s_{1})\zeta_{k}(t_{2},s_{2}). Let C(t,s)=\mathrm{cov}\big{(}X(t),X(s)\big{)} be the covariance kernel of $X$. Define a new kernel $\Pi$ such that

[TABLE]

Let $\rho_{1}\geq\rho_{2}\geq\cdots>0$ be the eigenvalues of $\Pi$ and $\{\phi_{j}:j\geq 1\}$ be the corresponding eigenfunctions. Therefore,

[TABLE]

Theorem 3.1

Assume that for any $\beta\in L_{2}([0,1]^{2})$

[TABLE]

for a positive constant $c$. Suppose that the eigenvalues $\{\rho_{k}:k\geq 1\}$ of the kernel $\Pi$ in (7) satisfy $\rho_{k}\asymp k^{-2r}$ for some constant $0<r<\infty$. Then,

[TABLE]

when $\lambda$ is of order $n^{-2r/(2r+1)}$.

Theorem 3.1 indicates that the convergence rate is determined by the decay rate of the eigenvalues of this new operator $\Pi$, which is jointly determined by both reproducing kernel $K$ and the covariance kernel $C$ as well as the alignment between $K$ and $C$ in a complicated way. This result has not been reported in the literature before. A close and related result is from Yuan & Cai (2010) who studied an optimal prediction risk for functional linear models, where the optimal rate depends on the decay rate of the eigenvalues of $K^{1/2}CK^{1/2}$. It is interesting to see, on the other hand, whether the convergence rate of $\hat{\beta}_{n}$ in Theorem 3.1 is optimal. In the following, we derive a minimax lower bound for the risk.

Theorem 3.2

Let $r$ be as in Theorem 3.1. Then the excess prediction risk satisfies

[TABLE]

where the infimum is taken over all possible predictors $\tilde{\eta}$ based on $\{(X_{i},Y_{i}):i=1,\ldots,n\}$.

Theorem 3.2 shows that the minimax lower bound of the convergence rate for the prediction risk is $n^{-{2r/2r+1}}$, which is determined by $r$ and the decay rate of the eigenvalues of $\Pi$. We have shown that this rate is achieved by our penalized estimator, and therefore our estimator is rate-optimal.

4 Numerical Experiments

We compared the proposed optimal penalized function-on-function regression (OPFFR) method with existing function-on-function linear regression models under two different designs. In a dense design, each curve was densely sampled at regularly-spaced common time points. We compared the OPFFR with two existing models. In a sparse design, each curve was irregularly and sparsely sampled at possibly different time points. We extended the OPFFR to this design by adding an extra pre-smoothing step and compared it with the FPCA model. In the first model (Ramsay & Silverman, 2005) for comparison, the coefficient function is estimated by penalizing its B-spline basis function expansion. This approach does not have the optimal mean prediction property and partially implemented in the fda package of R (linmod function) for the case of a fixed smoothing parameter. We shall add a search on the grid $10^{(-2\,:\,0.4\,:\,2)}$ for smoothing parameter selection to their implementation and denote this augmented approach by FDA. The coefficient function is represented in terms of 10 basis functions each for the $t$ and $s$ directions. The second model for comparison was the functional principal component analysis (hence denoted by FPCA) approach proposed by Yao et al. (2005b). The coefficient function is represented in terms of the leading functional principal components. This is implemented in the MatLab package PACE (FPCreg function) maintained by the UC-Davis research group. The Akaike information criterion (AIC) and fraction of variance explained (FVE) criterion were used to select the number of principal components for predictor and response respectively. The cutoff value for FVE was $0.9$. The ‘regular’ parameter was set to 2 for the dense design and 0 for the sparse design. No binning was performed.

4.1 Simulation Study

4.1.1 Dense Design

We simulated data according to model (2) with three scenarios.

• •

Scenario 1: The predictor functions are $X_{i}(s)=\sum_{k=1}^{50}(-1)^{(k+1)}k^{-1}Z_{ik}\vartheta_{1}(s,k)$, where $Z_{ik}$ is from the uniform distribution $U(-\sqrt{3},\sqrt{3})$, and $\vartheta_{1}(s,k)=1$ if $k=1$ and $\sqrt{2}\cos((k-1)\pi s)$ otherwise. The coefficient function $\beta(t,s)=e^{-(t+s)}$ is the exponential function of $t$ and $s$.

• •

Scenario 2: The predictor functions $X_{i}(s)$ are the same as those in Scenario 1 and the coefficient function $\beta(t,s)=4\sum_{k=1}^{50}(-1)^{(k+1)}k^{-2}\vartheta_{1}(t,k)\vartheta_{1}(s,k)$.

• •

Scenario 3: The predictor functions $X_{i}(s)$ are generated as $X_{i}(s)=\sum_{k=1}^{3}(-1)^{(k+1)}k^{-1}Z_{ik}\vartheta_{2}(s,k)$, where $\vartheta_{2}(s,k)=1$ if $k=3$ and $\sqrt{2}\cos(k\pi s)$ otherwise. The coefficient function $\beta(t,s)=4\sum_{k=1}^{3}(-1)^{(k+1)}k^{-2}\vartheta_{2}(t,k)\vartheta_{2}(s,k)$.

For each simulation scenario, we generated $n=30$ samples, each with 20 time points on the interval $(0,1)$. The random errors $\epsilon(t)$ were from a normal distribution with a constant variance $\sigma^{2}$. The value of $\sigma$ was adjusted to deliver three levels of signal-to-noise ratio (SNR$=0.5$, 5, and 10) in each scenario. To assess the mean prediction accuracy, we generated an additional $n^{*}=30$ predictor curves $\tilde{X}$ and computed the mean integrated squared error $\text{MISE}=1/n^{*}\sum_{i=1}^{n^{*}}\int_{0}^{1}(\eta_{\hat{\beta}}(\tilde{X}_{i},t)-\eta_{\beta}(\tilde{X}_{i},t))^{2}dt$, where $\hat{\beta}$ was the estimator obtained from the training data. We had 100 runs for each combination of scenario and SNR.

We applied the OPFFR, FDA and FPCA methods to the simulated data sets. Figure 3 displayed the perspective plots of the true coefficient functions in the three scenarios as well as their respective estimates for a single run with SNR$=10$. In the first two scenarios, both OPFFR and FDA did a decent job in recovering the true coefficient function although the FDA estimates were slightly oversmoothed. In both scenarios the FPCA estimates clearly suffered since the true coefficient function could not be effectively represented by the eigen-functions of the predictor processes.

Figure 4 gave the summary reports of performances in terms of MISEs based on 100 runs. When the signal to noise ratio is low, the OPFFR and FDA approaches had comparable performances. But when the signal to noise ratio increases, OPFFR showed clear advantage against FDA. The FPCA method failed to deliver competitive performance against the other two methods in all the settings due to its restrictive requirement of the effective representation of the coefficient function.

4.1.2 Sparse Design

In this section, we compared the performance of the proposed OPFFR method and the FPCA method regarding prediction error on sparsely, irregularly, and noisily observed functional data. To extend our method to sparsely and noisily observed data, we first applied the principal-component-analysis-through-conditional-expectation (PACE) method in Yao et al. (2005a) to the sparse functional data. Then we obtained a dense version of functional data by computing the PACE-fitted response and predictor functions at 50 selected time points for each curve. We applied the OPFFR method to these densely generated data and called this sparse extension to the OPFFR by the OPFFR-S method. The original OPFFR method, FPCA and OPFFR-S methods were all applied to the simulated data for comparison.

We first generated $n=200$ samples for both response and predictor functions in Scenario 3, each with 50 time points on interval $(0,1)$. To obtain different sparsity levels, we then randomly chose 5, 10 and 15 time points from the 50 ones for each curve independently. Normally distributed random errors were added to functional response and predictor with the SNR set to 10 in generating each pair of noisy response and predictor. The mean integrated squared error (MISE) was calculated based on additional $n^{*}=50$ predictor curves without random noises.

Figure 5 displayed the perspective plots of the true coefficient functions in the sparse scenario as well as their respective estimates for a single run with 10 sampled time points per curve. The OPFFR-S method and FPCA performed well in estimating the coefficient function. The estimate recovered by the original OPFFR method was a little oversmoothed. In Figure 6, the performance in terms of MISEs based on 100 runs was compared. The OPFFR-S method always had the best prediction performances at all the three sparsity levels. When the sparsity level was high (5 time points per curve), the original OPFFR method had a worse prediction performance than the FPCA. However, its prediction performance quickly picked up as the data became denser. When the sparsity level was 15 time points per curve, it actually delivered a better prediction performance than the FPCA. Such an interesting phenomenon was referred to as the “phase transistion” (Cai & Yuan, 2011; Wang et al., 2016).

5 Real Data Examples

We analyzed two real example in this section. We showed that our method had the numerical advantage over other approaches in terms of prediction accuracy in the analysis of the Canadian weather and histone regulation data. The results in the Canadian weather data, a dense design case, and the histone regulation data, a sparse design case, echoed with our findings in the simulation study. The smoothing parameters used in FDA for Canadian weather data were taken from the example codes in Ramsay et al. (2009) and seven basis functions were used for the $t$ and $s$ directions respectively. In the histone regulation data we selected the smoothing parameter for FDA by a grid search on $10^{(-5\,:\,1\,:\,5)}$ and used six basis functions each for the $t$ and $s$ directions. For the FPCA method, the ‘regular’ parameter was set to $2$ for the Canadian weather data and [math] for the histone regulation data. The other parameters for FDA and FPCA approaches were the same as those used in the simulation study.

5.1 Canadian Weather Data

We first look at the Canadian weather data (Ramsay & Silverman, 2005), a benchmark data set in functional data analysis. The main goal is to predict the log daily precipitation profile based on the daily temperature profile for a geographic location in Canada. The daily temperature and precipitation data averaged over 1960 to 1994 were recorded at 35 locations in Canada. We compared OPFFR with FDA and FPCA in terms of prediction performance defined by integrated squared error (ISE) $\int_{0}^{365}(Y_{i}(t)-\eta_{\hat{\beta}_{-i}}(X_{i},t))^{2}dt$, where $i=1,\cdots,35$ and $\hat{\beta}_{-i}$ was estimated by the dataset without the $i$th observation. For the convenience of calculation, we computed $\lVert Y_{i}(t)-\eta_{\hat{\beta}_{-i}}(X_{i},t)\rVert_{2}^{2}$ at a grid of values $t$ as the surrogate of ISE. Since the findings through the coefficient function estimates were similar to those in Ramsay & Silverman (2005), we only focused on the comparison of prediction performance. The summary in Table  1 clearly showed the numerical advantage of the proposed OPFFR method over the FDA and FPCA methods.

5.2 Histone Regulation Data

Nucleosomes, the basic units of DNA packaging in eukaryotic cells, consist of eight histone protein cores including two copies of H2A, H2B, H3, and H4. Besides the role as DNA scaffold, histones provide a complex regulatory platform for regulating gene activity (Wollmann et al., 2012). Focused study of the interaction between histones and gene activity may reveal how the organisms respond to the environmental changes. There are multiple sequence variants of histone proteins, which have some amino acid changes compared to their primary sequence, coexist in the same nucleus. For instance, in both plants and animals, there exist three variants of H3, the H3.1, the H3.3, and the centromere-specific CENP-A (CENH3) (Deal & Henikoff, 2011). Each variant shows distinct regulatory mechanisms over gene expression.

In this paper, an ultra-high throughput time course study was conducted to explore the interaction mechanism between the gene activity and histone variant, H3.3, during heat stress in Arabidopsis thaliana. In this study, the 12-day-old Arabidopsis seedlings that had been grown at $22\,^{\circ}{\rm C}$ were subject to heat stress of $38\,^{\circ}{\rm C}$, and plants were harvested at 7 different time points within 24 hours for RNA sequencing (RNA-seq) (Wang et al., 2009) and ChIP sequencing (ChIP-seq) (Johnson et al., 2007) experiments. We were interested in the genes responding to the heat shock, therefore 160 genes in response to heat (GO:0006951) pathway (Ashburner et al., 2000) were chosen. We selected 55 genes with the fold change above 0.5 at at least two consecutive time points in RNA-seq data. In ChIP-seq experiments, we calculated the mean of normalized read counts by taking the average of normalized read counts over seven time points for the region of 350 base pairs (bp) in the downstream of transcription start sites (TSS) of selected 55 genes. The normalized read counts over 350 bp from ChIP-seq and the normalized fragments per kilobase of transcript per million mapped reads (FPKM) (Trapnell et al., 2010) over seven time points from RNA-seq were used to measure the histone levels and gene expression levels respectively.

We applied the OPFFR, FDA and FPCA methods to histone regulation data in example 1.2. Since the gene expression levels were sparsely observed, we also applied the OPFFR-S method to the data. The comparison of the four methods is shown in Table 2. In the table, the standard deviation of ISEs was the only measure that neither the OPFFR nor the OPFFR-S was the most optimal. This was caused by a few observations where all the methods failed to make a good prediction and the OPFFR methods happened to have larger ISEs. In terms of all the other measures, the proposed OPFFR and OPFFR-S methods clearly showed the advantage in prediction accuracy again. Since the results from the OPFFR and OPFFR-S were comparable to each other, we chose to present all the following results based on the OPFFR analysis.

Figure 7 is the plot of the fitted coefficient function generated from our OPFFR method. For region between 300 bp and 350 bp, there was a strong negative influence of H3.3 on genes activity from half hour to 8 hours. It indicted that the loss of H3.3 might have the biological influence on the up-regulation of heat-induced genes. This negative correlation phenomenon was also observed after 30 minutes on the region of 250 bp to 300 bp between H3.3 and gene activity. In addition, the region from 50 bp to 150 bp had a positive effect on genes activity over time domain from 0 hour to half hour and 4 hours to 8 hours. Therefore, we provided a numerical evidence that heat-shock-induced transcription of genes in response to heat stress might be regulated via the epigenetic changes of H3.3, especially on the downstream region of TSS. The sample plots in Figure 8 showed a nice match of the predicted gene expression curves with the observed values.

6 Conclusion

In this article, we have presented a new analysis tool for modeling the relationship of a functional response against a functional predictor. The proposed method is more flexible and generally delivers a better numerical performance than the FPCA approach since it does not have the restrictive structural dependence assumption on the coefficient function. When compared with the penalized B-splines method, the proposed method has the theoretical advantage of possessing the optimal rate for mean prediction as well as some numerical advantage as shown in the numerical studies. Moreover, the Representer Theorem guarantees an exact solution to the penalized least squares, a property that is not shared by the existing penalized function-on-function regression models. The application of our method to a histone regulation study provided numerical evidence that the changes in H3.3 might regulate some genes through transcription regulations. Although such a finding sheds light on the relationship between histone variant H3.3 and gene activity, the details of the regulation process are still unknown and merit further investigations. For instance, we may investigate how the H3.3 organizes the chromatins to up-regulate those active genes. Such investigations would call for more collaborations between statisticians and biologists.

When the regression model has a scalar response against one or more functional predictors, methods other than the roughness penalty approach are available to overcome the inefficient basis representation drawback in the FPCA method. For example, Delaigle et al. (2012) considered a partial least squares (PLS) based approach. Ferré & Yao (2003) and Yao et al. (2015) translated the idea of sufficient dimension reduction (SDR) into the setting of functional regression models. Intuitively, these methods might be more efficient in their selection of the principal component basis functions since they incorporate the response information into consideration. However, our experiments with a functional response version of the functional PLS (Preda & Saporta, 2005), not shown here due to space limit, did not look so promising. Therefore, further investigation in this direction is surely needed.