# Functional linear regression via canonical analysis

###### Abstract

We study regression models for the situation where both dependent and independent variables are square-integrable stochastic processes. Questions concerning the definition and existence of the corresponding functional linear regression models and some basic properties are explored for this situation. We derive a representation of the regression parameter function in terms of the canonical components of the processes involved. This representation establishes a connection between functional regression and functional canonical analysis and suggests alternative approaches for the implementation of functional linear regression analysis. A specific procedure for the estimation of the regression parameter function using canonical expansions is proposed and compared with an established functional principal component regression approach. As an example of an application, we present an analysis of mortality data for cohorts of medflies, obtained in experimental studies of aging and longevity.

0 \volume16 \issue3 2010 \firstpage705 \lastpage729 \doi10.3150/09-BEJ228 \newproclaimdefinitionDefinition[section] \newproclaimconditionCondition \newremarkexample[thm]Example

Regression via canonical analysis

1]\fnmsGuozhong \snmHe\thanksref1label=e1]ghe@dhs.ca.gov,
2]\fnmsHans-Georg \snmMüller\corref\thanksref2,e2label=e2,mark]mueller@wald.ucdavis.edu,
2]\fnmsJane-Ling \snmWang\thanksref2,e3label=e3,mark]wang@wald.ucdavis.edu

and 2]\fnmsWenjing \snmYang\thanksref2,e4label=e4,mark]wyang@wald.ucdavis.edu

canonical components \kwdcovariance operator \kwdfunctional data analysis \kwdfunctional linear model \kwdlongitudinal data \kwdparameter function \kwdstochastic process

## 1 Introduction

With the advancement of modern technology, data sets which contain repeated measurements obtained on a dense grid are becoming ubiquitous. Such data can be viewed as a sample of curves or functions and are referred to as functional data. We consider here the extension of the linear regression model to the case of functional data. In this extension, both predictors and responses are random functions rather than random vectors. It is well known (Ramsay and Dalzell (1991); Ramsay and Silverman (2005)) that the traditional linear regression model for multivariate data, defined as

(1) |

may be extended to the functional setting by postulating the model, for ,

(2) |

Writing all vectors as row vectors in the classical model (1), and are random vectors in , is a random vector in , and and are, respectively, and matrices containing the regression parameters. The vector has the usual interpretation of an error vector, with and , denoting the identity matrix. In the functional model (2), random vectors and in (1) are replaced by random functions defined on the intervals and . The extension of the classical linear model (1) to the functional linear model (2) is obtained by replacing the matrix operation on the right-hand side of (1) with an integral operator in (2). In the original approach of Ramsay and Dalzell (1991), a penalized least-squares approach using L-splines was adopted and applied to a study in temperature-precipitation patterns, based on data from Canadian weather stations.

The functional regression model (2) for the case of scalar responses has attracted much recent interest (Cardot and Sarda (2005); Müller and Stadtmüller (2005); Hall and Horowitz (2007)), while the case of functional responses has been much less thoroughly investigated (Ramsay and Dalzell (1991); Yao, Müller and Wang (2005b)). Discussions on various approaches and estimation procedures can be found in the insightful monograph of Ramsay and Silverman (2005). In this paper, we propose an alternative approach to predict from , by adopting a novel canonical representation of the regression parameter function . Several distinctive features of functional linear models emerge in the development of this canonical expansion approach.

It is well known that in the classical multivariate linear model, the regression slope parameter matrix is uniquely determined by , as long as the covariance matrix is invertible. In contrast, the corresponding parameter function , appearing in (2), is typically not identifiable. This identifiability issue is discussed in Section 2. It relates to the compactness of the covariance operator of the process which makes it non-invertible. In Section 2, we demonstrate how restriction to a subspace allows this problem to be circumvented. Under suitable restrictions, the components of model (2) are then well defined.

Utilizing the canonical decomposition in Theorem 3 below leads to an alternative approach to estimating the parameter function . The canonical decomposition links and through their functional canonical correlation structure. The corresponding canonical components form a bridge between canonical analysis and linear regression modeling. Canonical components provide a decomposition of the structure of the dependency between and and lead to a natural expansion of the regression parameter function , thus aiding in its interpretation. The canonical regression decomposition also suggests a new family of estimation procedures for functional regression analysis. We refer to this methodology as functional canonical regression analysis. Classical canonical correlation analysis (CCA) was introduced by Hotelling (1936) and was connected to function spaces by Hannan (1961). Substantial extensions and connections to reproducing kernel Hilbert spaces were recently developed in Eubank and Hsing (2008); for other recent developments see Cupidon et al. (2007).

Canonical correlation is known not to work particularly well for very high-dimensional multivariate data, as it involves an inverse problem. Leurgans, Moyeed and Silverman (1993) tackled the difficult problem of extending CCA to the case of infinite-dimensional functional data and discussed the precarious regularization issues which are faced; He, Müller and Wang (2003, 2004) further explored various aspects and proposed practically feasible regularization procedures for functional CCA. While CCA for functional data is worthwhile, but difficult to implement and interpret, the canonical approach to functional regression is here found to compare favorably with the well established principal-component-based regression approach in an example of an application (Section 5). This demonstrates a potentially important new role for canonical decompositions in functional regression analysis. The functional linear model (2) includes the varying coefficient linear model studied in Hoover et al. (1998) and Fan and Zhang (2000) as a special case, where ; here, is a delta function centered at and is the varying coefficient function. Other forms of functional regression models with vector-valued predictors and functional responses were considered by Faraway (1997), Shi, Weiss and Taylor (1996), Rice and Wu (2000), Chiou, Müller and Wang (2003) and Ritz and Streibig (2009).

The paper is organized as follows. Functional canonical analysis and functional linear models for -processes are introduced in Section 2. Sufficient conditions for the existence of functional normal equations are given in Proposition 2. The canonical regression decomposition and its properties are the theme of Section 3. In Section 4, we propose a novel estimation technique to obtain regression parameter function estimates based on functional canonical components. The regression parameter function is the basic model component of interest in functional linear models, in analogy to the parameter vector in classical linear models. The proposed estimation method, based on a canonical regression decomposition, is contrasted with an established functional regression method based on a principal component decomposition. These methods utilize a dimension reduction step to regularize the solution of the inverse problems posed by both functional regression and functional canonical analysis. As a selection criterion for tuning parameters, such as bandwidths or numbers of canonical components, we use minimization of prediction error via leave-one-curve-out cross-validation (Rice and Silverman (1991)). The proposed estimation procedures are applied to mortality data obtained for cohorts of medflies (Section 5). Our goal in this application is to predict a random trajectory of mortality for a female cohort of flies from the trajectory of mortality for a male cohort which was raised in the same cage. We find that the proposed functional canonical regression method gains an advantage over functional principal component regression in terms of prediction error.

## 2 Functional linear regression and the functional normal equation

In this section, we explore the formal setting as well as identifiability issues for functional linear regression models. Both response and predictor functions are considered to come from a sample of pairs of random curves. A basic assumption is that all random curves or functions are square-integrable stochastic processes. Consider a measure on a real index set and let be the class of real-valued functions such that . This is a Hilbert space with the inner product and we write if . The index set can be a set of time points, such as , a compact interval or even a rectangle formed by two intervals and , . We focus on index sets that are either compact real intervals or compact rectangles in and consider to be the Lebesgue measure on or . Extensions to other index sets and other measures are self-evident. An -process is a stochastic process , , with for all .

Let and .

Processes are subject to a functional linear model if

(3) |

where is the parameter function, is a random error process with for , and and are uncorrelated, in the sense that for all .

Without loss of generality, we assume from now on that all processes considered have zero mean functions, and for all , . We define the regression integral operator by

Equation (3) can then be rewritten as

(4) |

Denote the auto- and cross-covariance functions of and by

The autocovariance operator of is the integral operator , defined by

Replacing by , , we analogously define operators and , similarly . Then and are compact, self-adjoint and non-negative definite operators, and and are compact operators (Conway (1985)). We refer to He et al. (2003) for a discussion of various properties of these operators.

Another linear operator of interest is the integral operator ,

(5) |

The operator equation

(6) |

is a direct extension of the least-squares normal equation and may be referred to as the functional population normal equation.

###### Proposition \thedefinition

The following statements are equivalent for a function : {longlist}[(a)]

satisfies the linear model (4);

is a solution of the functional normal equation (6);

minimizes among all .

The proof is found Section 7. In the infinite-dimensional case, the operator is a Hilbert–Schmidt operator in the Hilbert space , according to Proposition 6 below. A problem we face is that it is known from functional analysis that a bounded inverse does not exist for such operators. A consequence is that the parameter function in (3), (4) is not identifiable without additional constraints. In a situation where the inverse of the covariance matrix does not exist in the multivariate case, a unique solution of the normal equation always exists within the column space of and this solution then minimizes on that space. Our idea to get around the non-invertibility issue in the functional infinite-dimensional case is to extend this approach for the non-invertible multivariate case to the functional case. Indeed, as is demonstrated in Theorem 2 below, under the additional Condition 2, the solution of (6) exists in the subspace defined by the range of . This unique solution indeed minimizes .

We will make use of the Karhunen–Loève decompositions (Ash and Gardner (1975)) for -processes and ,

(7) |

with random variables , , , and orthonormal families of -functions and . Here, , , and are the eigenvalues and eigenfunctions of the covariance operators and , respectively, with , . Note that is the Kronecker symbol with for , for .

We consider a subset of on which inverses of the operator can be defined. As a Hilbert–Schmidt operator, is compact and therefore not invertible on According to Conway (1985), page 50, the range of

is characterized by

(8) |

where Defining

we find that is a one-to-one mapping from the vector space onto the vector space Thus, restricting to a subdomain defined by the subspace we can define its inverse for as

(9) |

then satisfies the usual properties of an inverse, in the sense that for all and for all

The following Condition 2 for processes is of interest.

The -processes and with Karhunen–Loève decompositions (7) satisfy

If 2 is satisfied, then the solution to the non-invertibility problem as outlined above is viable in the functional case, as demonstrated by the following basic result on functional linear models.

[(Basic theorem for functional linear models)] A unique solution of the linear model (4) exists in if and only if and satisfy Condition 2. In this case, the unique solution is of the form

(10) |

As a consequence of Proposition 2, solutions of the functional linear model (4), solutions of the functional population normal equation (6) and minimizers of are all equivalent and allow the usual projection interpretation.

###### Proposition \thethm

It is well known that in a finite-dimensional situation, the linear model (6) always has a unique solution in the column space of , which may be obtained by using a generalized inverse of the matrix . However, in the infinite-dimensional case, such a solution does not always exist. The following example demonstrates that a pair of -processes does not necessarily satisfy Condition 2. In this case, the linear model (6) does not have a solution. {example} Assume processes and have Karhunen–Loève expansions (7), where the random variables , satisfy

(11) |

and let

(12) |

As shown in He et al. (2003), (11) and (12) can be satisfied by a pair of -processes with appropriate operators , and . Then

and, therefore, Condition 2 is not satisfied.

## 3 Canonical regression analysis

Canonical analysis is a time-honored tool for studying the dependency between the components of a pair of random vectors or stochastic processes; for multivariate stationary time series, its utility was established in the work of Brillinger (1985). In this section, we demonstrate that functional canonical decomposition provides a useful tool to represent functional linear models. The definition of functional canonical correlation for -processes is as follows.

The first canonical correlation and weight functions and for -processes and are defined as

(13) |

where and are subject to

(14) |

for . The th canonical correlation and weight functions , for processes and for are defined as

where and are subject to (14) for and

for We refer to and as the th canonical variates and to as the th canonical components.

It has been shown in He et al. (2003) that canonical correlations do not exist for all -processes, but that Condition 3 below is sufficient for the existence of canonical correlations and weight functions. We remark that Condition 3 implies Condition 2.

Let and be -processes, with Karhunen–Loève decompositions (7) satisfying

The proposed functional canonical regression analysis exploits features of functional principal components and of functional canonical analysis. In functional principal component analysis, one studies the structure of an -process via its decomposition into the eigenfunctions of its autocovariance operator, the Karhunen–Loève decomposition (Rice and Silverman (1991)). In functional canonical analysis, the relation between a pair of -processes is analyzed by decomposing the processes into their canonical components. The idea of canonical regression analysis is to expand the regression parameter function in terms of functional canonical components for predictor and response processes. The canonical regression decomposition (Theorem 3) below provides insights into the structure of the regression parameter functions and not only aids in the understanding of functional linear models, but also leads to promising estimation procedures for functional regression analysis. The details of these estimation procedures will be discussed in Section 4. We demonstrate in Section 5 that these estimates can lead to competitive prediction errors in a finite-sample situation.

We now state two key results. The first of these (Theorem 3) provides the canonical decomposition of the cross-covariance function of processes and . This result plays a central role in the solution of the population normal equation (6). This solution is referred to as canonical regression decomposition and it leads to an explicit representation of the underlying regression parameter function of the functional linear model (4). The decomposition is in terms of functional canonical correlations and canonical weight functions and . Given a predictor process , we obtain, as a consequence, an explicit representation for , where is as in (4). For the following main results, we refer to the definitions of , , , , in Definition 3. All proofs are found in Section 7.

[(Canonical decomposition of cross-covariance function)] Assume that -processes and satisfy Condition 3. The cross-covariance function then allows the following representation in terms of canonical correlations and weight functions and :

(15) |

[(Canonical regression decomposition)] Assume that the -processes and satisfy Condition 3. One then obtains, for the regression parameter function (10), the following explicit solution:

(16) |

To obtain the predicted value of the response process , we use the linear predictor

(17) |

This canonical regression decomposition leads to approximations of the regression parameter function and the predicted process via a finitely truncated version of the canonical expansions (16) and (17). The following result provides approximation errors incurred from finite truncation. Thus, we have a vehicle to achieve practically feasible estimation of and associated predictions (Section 4).

For , let be the finitely truncated version of the canonical regression decomposition (16) for and define . Then,

(18) |

with . Moreover,

and

(19) |

In finite-sample implementations, to be explored in the next two sections, truncation as in (18) is a practical necessity; this requires a choice of suitable truncation parameters.

## 4 Estimation procedures

### 4.1 Preliminaries

Estimating the regression parameter function and obtaining fitted processes from the linear model (2) based on a sample of curves is central to the implementation of functional linear models. In practice, data are observed at discrete time points and we temporarily assume, for simplicity, that the time points are the same for all observed predictor curves and are equidistantly spaced over the domain of the data. Analogous assumptions are made for the time points where the response curves are sampled. Thus, the original observations are , , where is an -dimensional vector sampled at time points , and is an -dimensional vector sampled at time points . We assume that and are both large. Without going into any analytical details, we compare the finite-sample behavior of two functional regression methods, one of which utilizes the canonical decomposition for regression and the other a well established direct principal component approach to implement functional linear regression.

The proposed practical version of functional regression analysis through functional canonical regression analysis (FCR) is discussed in Section 4.2. This method is compared with a more standard functional linear regression implementation that is based on principal components and referred to as functional principal regression (FPR), in Section 4.3. For the choice of the smoothing parameters for the various smoothing steps, we adopt leave-one-curve-out cross-validation (Rice and Silverman (1991)). Smoothing is implemented by local linear fitting for functions and surfaces (Fan and Gijbels (1996)), minimizing locally weighted least squares.

In a pre-processing step, all observed process data are centered by subtracting the cross-sectional means , and analogously for . If the data are not sampled on the same grid for different individuals, a smoothing step may be added before the cross-sectional average is obtained. As in the previous sections, we use in the following the notation to denote centered processes and trajectories.

When employing the Karhunen–Loève decomposition (7), we approximate observed centered processes by the fitted versions

(20) |

where and are the estimated first smoothed eigenfunctions for the random processes and , respectively, with the corresponding estimated eigenscores and for the th subject. We obtain these estimates as described in Yao et al. (2005a). Related estimation approaches, such as those of Rice and Silverman (1991) or Ramsay and Silverman (2005), could alternatively be used.

### 4.2 Functional canonical regression (FCR)

To obtain functional canonical correlations and the corresponding weight functions as needed for FCR, we adopt one of the methods proposed in He et al. (2004). In preliminary studies, we determined that the eigenbase method as described there yielded the best performance for regression applications, with the Fourier base method a close second. Adopting the eigenbase method, the implementation of FCR is as follows: {longlist}

Starting with the eigenscore estimates as in (20), estimated raw functional canonical correlations and -dimensional weight vectors , , are obtained by applying conventional numerical procedures of multivariate canonical analysis to the estimated eigenscore vectors and . This works empirically well for moderately sized values of , as typically obtained from automatic selectors.

Smooth weight function estimates are then obtained as

where

The estimated regression parameter function is obtained according to (16) by

where is an estimate of the covariance function of , obtained by two-dimensional smoothing of the empirical autocovariances of . This estimate is obtained as described in Yao et al. (2005a). Since the data are regularly sampled, the above integrals are easily obtained by the approximations , with defined analogously to in (22) below.

Fitted/predicted processes

(21) |

are obtained, where the integral is again evaluated numerically by

(22) |

Here, is chosen such that .

This procedure depends on two tuning parameters, a bandwidth for the smoothing steps (which are defined in detail, e.g., in Yao et al. (2005a)) and the number of canonical components that are included. These tuning parameters may be determined by leave-one-out cross-validation (Rice and Silverman (1991)) as follows. With , the th leave-one-out estimate for is

(23) |

where is the th canonical correlation, and and are the th weight function untransformed and transformed with the covariance operator, respectively, all obtained while leaving out the data for the th subject. Computation of these estimates follows steps (iii) and (iv) above, using tuning parameter , and omitting the th pair of observed curves . The average leave-one-out squared prediction error is then

(24) |

The cross-validation procedure then selects the tuning parameter that minimizes the approximate average prediction error,

where is obtained by replacing the integrals on the right-hand side of (24) by sums of the type (22).

### 4.3 Functional principal component regression (FPR)

Yao et al. (2005b) considered an implementation of functional linear regression whereby one uses functional principal component analysis for predictor and response functions separately, followed by simple linear regressions of the response principal component scores on the predictor scores. We adopt this approach as FPR.

Briefly, defining , this approach is based on representations

of the regression parameter function , where

(25) |

for all and .

For estimation, one first obtains a smooth estimate of the cross-covariance by smoothing sample cross-covariances, for example, by the method described in Yao et al. (2005b). This leads to estimates of by plugging in estimates for and for eigenfunctions (as described in Section 4.2), in combination with approximating the integrals in (25) by appropriate sums. One may then use these estimates in conjunction with estimates of eigenvalues to arrive at the estimate of the regression parameter function given by

For further details about numerical implementations, we refer to Yao et al. (2005b).

## 5 Application to medfly mortality data

In this section, we present an application to age-at-death data that were collected for cohorts of male and female medflies in a biodemographic study of survival and mortality patterns of cohorts of male and female Mediterranean fruit flies (Ceratitis capitata; for details, see Carey et al. (2002)). A point of interest in this study is the relation of mortality trajectories between male and female medflies which were raised in the same cage. One specifically desires to quantify the influence of male survival on female survival. This is of interest because female survival determines the number of eggs laid and thus reproductive success of these flies. We use a subsample of the data generated by this experiment, comprising 46 cages of medflies, to address these questions. Each cage contains both a male and a female cohort, consisting each of approximately 4000 male and 4000 female medflies. These flies were raised in the shared cage from the time of eclosion. For each cohort, the number of flies alive at the beginning of each day was recorded, simply by counting the dead flies on each day; we confined the analysis to the first 40 days. The observed processes and , are the estimated random hazard functions for male and female cohorts, respectively. All deaths are fully observed so that censoring is not an issue. In a pre-processing step, cohort-specific hazard functions were estimated nonparametrically from the lifetable data, implementing the transformation approach described in Müller et al. (1997a).

PE | |||
---|---|---|---|

FCR | 1.92 | 3 | 0.0100 |

FPR | 1.65 | 3 | 0.0121 |

A functional linear model was used to study the specific influence of male mortality on female mortality for flies that were raised in the same cage, with the hazard function of males as predictor process and that of females as response process. We applied both the proposed regression via canonical representation (FCR) and the more conventional functional regression based on principal components (FPR), implementing the estimation procedures described in the previous section. Tuning parameters were selected by cross-validation. Table 1 lists the average squared prediction error (PE) (24) obtained by the leave-one-out technique. For this application, the FCR procedure is seen to perform about 20% better than FPR in terms of PE.

The estimated regression parameter surface that is obtained for the FCR regression when choosing the cross-validated values for and , as given in Table 1, is shown in Figure 1. The shape of the regression surface indicates that female mortality at later ages is very clearly affected by male mortality throughout male lifespan, while female mortality at very early ages is not much influenced by male mortality. The effect of male mortality on female mortality is periodically elevated, as evidenced by the bumps visible in the surface. The particularly influential predictive periods are male mortality around days 10 and 20, which then has a particularly large influence on female mortality around days 15 and 25, that is, about five days later, and, again, around days 35 and 40, judging from the locations of the peaks in the surface of . In contrast, enhanced male mortality around day 30 leads to lessened female mortality throughout, while enhanced male mortality at age 40 is associated with higher older-age female mortality. These observations point to the existence of periodic waves of mortality, first affecting males and subsequently females. While some of the waves of increased male mortality tend to be associated with subsequently increased female mortality, others are associated with subsequently decreased female mortality.

These waves of mortality might be related to the so-called “vulnerable periods” that are characterized by locally heightened mortality (Müller et al. (1997b)). One such vulnerable period occurs around ages 10 and 20, and the analysis suggests that heightened male mortality during these phases is indicative of heightened female mortality. In contrast, heightened male mortality during a non-vulnerable period such as the time around 30 days seems to be associated with lower female mortality. A word of caution is in order as no inference methods are available to establish that the bumps observed in are real, so one cannot exclude the possibility that these bumps are enhanced by random fluctuations in the data.

Examples of observed, as well as predicted, female mortality trajectories for three randomly selected pairs of cohorts (male and female flies raised in the same cages) are displayed in Figure 2. The predicted female trajectories were constructed by applying both regression methods (FCR and FPR) with the leave-one-out technique. The prediction of an individual response trajectory from a predictor trajectory cannot, of course, be expected to be very close to the actually observed response trajectory, due to the extra random variation that is a large inherent component of response variability; this is analogous to the situation of predicting an individual response in the well-known simple linear regression case. Nevertheless, overall, FCR predictions are found to be closer to the target.

We note the presence of a “shoulder” at around day 20 for the three female mortality curves. This “shoulder” is related to the wave phenomenon visible in as discussed above and corresponds to a phase of elevated female mortality. The functional regression method based on FCR correctly predicts the shoulder effect and its overall shape in female mortality. At the rightmost points, for ages near 40 days, the variability of the mortality trajectories becomes large, posing extra difficulties for prediction in the right tail of the trajectories.

## 6 Additional results

Theorems 6 and 6 in this section provide a functional analog to the sums-of-squares decomposition of classical regression analysis. In addition, we provide two results characterizing the regression operators . We begin with two auxiliary results which are taken from He et al. (2003). The first of these characterizes the correlation operator between processes and .

###### Lemma \thethm

Assume that the -processes and satisfy Condition 3. The correlation operator can then be extended continuously to a Hilbert–Schmidt operator on to . Hence, is also a Hilbert–Schmidt operator with a countable number of non-zero eigenvalues and eigenfunctions , , Then: {longlist}

, and both and are -functions;

;

;

.

One of the main results in He et al. (2003) reveals that the -processes and can be expressed as sums of uncorrelated component functions and the correlation between the th components of the expansion is the th corresponding functional canonical correlation between the two processes.

###### Lemma \thethm ((Canonical decomposition))

Assume -processes and satisfy Condition 3. There then exists a decomposition: {longlist}

where

The index stands for canonical decomposition with components, and are as in Definition 3. Here, and share the same first canonical components, and and are uncorrelated, that is,

Let and Then

where . Here, and share the same canonical components, , and and are uncorrelated. Moreover, if forms a basis of the closure of the domain of and if forms a basis of the closure of the domain of .

Since the covariance operators of -processes are non-negative self-adjoint, they can be ordered as follows. The definitions of are in (17), (18) and Lemma 6(b), respectively.

For .

In multiple regression analysis, the ordering of the operators in Theorem 6 is related to the ordering of regression models in terms of a notion analogous to the regression sum of squares (SSR). The canonical regression decomposition provides information about the model in terms of its canonical components. Our next result describes the canonical correlations between observed and fitted processes. This provides an extension of the coefficient of multiple determination, an important quantity in classical multiple regression analysis, to the functional case; compare also Yao et al. (2005b).

Assume that -processes and satisfy Condition 3. The canonical correlations and weight functions for the pair of observed and fitted response processes are then and the corresponding -component (or -component) canonical decomposition for as defined in Lemma 6 for and denoted here by (or ), is equivalent to the process or given in Theorem 3, that is,

(26) |

We note that if is a scalar, then