Functional data analysis

Hilbertian random variables

In the Hilbert space viewpoint, one considers an ${\displaystyle H}$-valued random element ${\displaystyle X}$, where ${\displaystyle H}$ is a separable Hilbert space such as the space of square-integrable functions ${\displaystyle L^{2}[0,1]}$. Under the integrability condition that ${\displaystyle \mathbb {E} \|X\|_{L^{2}}^{2}=\mathbb {E} (\int _{0}^{1}|X(t)|^{2}dt)<\infty }$, one can define the mean of ${\displaystyle X}$ as the unique element ${\displaystyle \mu \in H}$ satisfying

${\displaystyle \mathbb {E} \langle X,h\rangle =\langle \mu ,h\rangle ,\qquad h\in H.}$

This formulation is the Pettis integral but the mean can also be defined as Bochner integral ${\displaystyle \mu =\mathbb {E} X}$. Under the integrability condition that ${\displaystyle \mathbb {E} \|X\|_{L^{2}}^{2}}$ is finite, the covariance operator of ${\displaystyle X}$ is a linear operator ${\displaystyle {\mathcal {C}}:H\to H}$ that is uniquely defined by the relation

${\displaystyle {\mathcal {C}}h=\mathbb {E} [\langle h,X-\mu \rangle (X-\mu )],\qquad h\in H,}$

or, in tensor form, ${\displaystyle {\mathcal {C}}=\mathbb {E} [(X-\mu )\otimes (X-\mu )]}$. The spectral theorem allows to decompose ${\displaystyle X}$ as the Karhunen-Loève decomposition

${\displaystyle X=\mu +\sum _{i=1}^{\infty }\langle X,\varphi _{i}\rangle \varphi _{i},}$

where ${\displaystyle \varphi _{i}}$ are eigenvectors of ${\displaystyle {\mathcal {C}}}$, corresponding to the nonnegative eigenvalues of ${\displaystyle {\mathcal {C}}}$, in a non-increasing order. Truncating this infinite series to a finite order underpins functional principal component analysis.

Stochastic processes

The Hilbertian point of view is mathematically convenient, but abstract; the above considerations do not necessarily even view ${\displaystyle X}$ as a function at all, since common choices of ${\displaystyle H}$ like ${\displaystyle L^{2}[0,1]}$ and Sobolev spaces consist of equivalence classes, not functions. The stochastic process perspective views ${\displaystyle X}$ as a collection of random variables

${\displaystyle \{X(t)\}_{t\in [0,1]}}$

indexed by the unit interval (or more generally interval ${\displaystyle {\mathcal {T}}}$). The mean and covariance functions are defined in a pointwise manner as

${\displaystyle \mu (t)=\mathbb {E} X(t),\qquad \Sigma (s,t)={\textrm {Cov}}(X(s),X(t)),\qquad s,t\in [0,1]}$

(if ${\displaystyle \mathbb {E} [X(t)^{2}]<\infty }$ for all ${\displaystyle t\in [0,1]}$).

Under the mean square continuity, ${\displaystyle \mu }$ and ${\displaystyle \Sigma }$ are continuous functions and then the covariance function ${\displaystyle \Sigma }$ defines a covariance operator ${\displaystyle {\mathcal {C}}:H\to H}$ given by

$${\displaystyle ({\mathcal {C}}f)(t)=\int _{0}^{1}\Sigma (s,t)f(s)\,\mathrm {d} s.}$$

1

The spectral theorem applies to ${\displaystyle {\mathcal {C}}}$, yielding eigenpairs ${\displaystyle (\lambda _{j},\varphi _{j})}$, so that in tensor product notation ${\displaystyle {\mathcal {C}}}$ writes

${\displaystyle {\mathcal {C}}=\sum _{j=1}^{\infty }\lambda _{j}\varphi _{j}\otimes \varphi _{j}.}$

Moreover, since ${\displaystyle {\mathcal {C}}f}$ is continuous for all ${\displaystyle f\in H}$, all the ${\displaystyle \varphi _{j}}$ are continuous. Mercer's theorem then states that

${\displaystyle \sup _{s,t\in [0,1]}\left|\Sigma (s,t)-\sum _{j=1}^{K}\lambda _{j}\varphi _{j}(s)\varphi _{j}(t)\right|\to 0,\qquad K\to \infty .}$

Finally, under the extra assumption that ${\displaystyle X}$ has continuous sample paths, namely that with probability one, the random function ${\displaystyle X:[0,1]\to \mathbb {R} }$ is continuous, the Karhunen-Loève expansion above holds for ${\displaystyle X}$ and the Hilbert space machinery can be subsequently applied. Continuity of sample paths can be shown using Kolmogorov continuity theorem.

1. Fully observed functions without noise at arbitrarily dense grid

Measurements ${\displaystyle Y_{it}=X_{i}(t)}$ available for all ${\displaystyle t\in {\mathcal {I}},\,i=1,\ldots ,n}$

Often unrealistic but mathematically convenient.

Real life example: Tecator spectral data.^[7]

2. Densely sampled functions with noisy measurements (dense design)

Measurements ${\displaystyle Y_{ij}=X_{i}(T_{ij})+\varepsilon _{ij}}$, where ${\displaystyle T_{ij}}$ are recorded on a regular grid,

${\displaystyle T_{i1},\ldots ,T_{iN_{i}}}$, and ${\displaystyle N_{i}\rightarrow \infty }$ applies to typical functional data.

Real life example: Berkeley Growth Study Data and Stock data

3. Sparsely sampled functions with noisy measurements (longitudinal data)

Measurements ${\displaystyle Y_{ij}=X_{i}(T_{ij})+\varepsilon _{ij}}$, where ${\displaystyle T_{ij}}$ are random times and their number ${\displaystyle N_{i}}$ per subject is random and finite.

Real life example: CD4 count data for AIDS patients.^[9]

Functional regression models with scalar response

Replacing the vector covariate ${\displaystyle X}$ and the coefficient vector ${\displaystyle \beta }$ in model (3) by a centered functional covariate ${\displaystyle X^{c}(t)=X(t)-\mu (t)}$ and coefficient function ${\displaystyle \beta =\beta (t)}$ for ${\displaystyle t\in [0,1]}$ and replacing the inner product in Euclidean space by that in Hilbert space ${\displaystyle L^{2}}$, one arrives at the functional linear model

$${\displaystyle Y=\beta _{0}+\langle X^{c},\beta \rangle +\varepsilon =\beta _{0}+\int _{0}^{1}X^{c}(t)\beta (t)\,dt+\varepsilon .}$$

4

The simple functional linear model (4) can be extended to multiple functional covariates, ${\displaystyle \{X_{j}\}_{j=1}^{p}}$, also including additional vector covariates ${\displaystyle Z=(Z_{1},\cdots ,Z_{q})}$, where ${\displaystyle Z_{1}=1}$, by

$${\displaystyle Y=\langle Z,\theta \rangle +\sum _{j=1}^{p}\int _{0}^{1}X_{j}^{c}(t)\beta _{j}(t)\,dt+\varepsilon ,}$$

5

where ${\displaystyle \theta \in \mathbb {R^{q}} }$ is regression coefficient for ${\displaystyle Z}$, the domain of ${\displaystyle X_{j}}$ is ${\displaystyle [0,1]}$, ${\displaystyle X_{j}^{c}}$ is the centered functional covariate given by ${\displaystyle X_{j}^{c}(t)=X_{j}(t)-\mu _{j}(t)}$, and ${\displaystyle \beta _{j}}$ is regression coefficient function for ${\displaystyle X_{j}^{c}}$, for ${\displaystyle j=1,\ldots ,p}$. Models (4) and (5) have been studied extensively.^[10][11][12]

Functional regression models with functional response

Consider a functional response ${\displaystyle Y(s)}$ on ${\displaystyle [0,1]}$ and multiple functional covariates ${\displaystyle X_{j}(t)}$, ${\displaystyle t\in [0,1]}$, ${\displaystyle j=1,\ldots ,p}$. Two major models have been considered in this setup.^[13][7] One of these two models, generally referred to as functional linear model (FLM), can be written as:

$${\displaystyle Y(s)=\alpha _{0}(s)+\sum _{j=1}^{p}\int _{0}^{1}\alpha _{j}(s,t)X_{j}^{c}(t)\,dt+\varepsilon (s),\ {\text{for}}\ s\in [0,1]}$$

6

where ${\displaystyle \alpha _{0}(s)}$ is the functional intercept, for ${\displaystyle j=1,\ldots ,p}$ , ${\displaystyle X_{j}^{c}(t)=X_{j}(t)-\mu _{j}(t)}$ is a centered functional covariate on ${\displaystyle [0,1]}$, ${\displaystyle \alpha _{j}(s,t)}$ is the corresponding functional slopes with same domain, respectively, and ${\displaystyle \varepsilon (s)}$ is usually a random process with mean zero and finite variance.^[13] In this case, at any given time ${\displaystyle s\in [0,1]}$, the value of ${\displaystyle Y}$, i.e., ${\displaystyle Y(s)}$, depends on the entire trajectories of ${\displaystyle \{X_{j}(t)\}_{j=1}^{p}}$. Model (6) has been studied extensively.^{[14][15][16][17][18]}

Concurrent regression models

This model is given by

$${\displaystyle Y(s)=\beta _{0}(s)+\sum _{j=1}^{p}\beta _{j}(s)X_{j}(s)+\varepsilon (s),\ {\text{for}}\ s\in [0,1],}$$

7

where ${\displaystyle X_{1},\ldots ,X_{p}}$ are functional covariates on ${\displaystyle [0,1]}$, ${\displaystyle \beta _{0},\beta _{1},\ldots ,\beta _{p}}$ are the coefficient functions defined on the same interval and ${\displaystyle \varepsilon (s)}$ is usually assumed to be a random process with mean zero and finite variance.^[13] This model assumes that the value of ${\displaystyle Y(s)}$ depends on the current value of ${\displaystyle \{X_{j}(s)\}_{j=1}^{p}}$ only and not the history ${\displaystyle \{X_{j}(t):t\leq s\}_{j=1}^{p}}$ or future value. Hence, it is a "concurrent regression model", which is also referred as "varying-coefficient" model. Further, various estimation methods have been proposed.^{[19][20][21][22][23][24]}

Functional polynomial regression models

Functional polynomial regression models may be viewed as a natural extension of the Functional Linear Models (FLMs) with scalar responses, analogous to extending linear regression model to polynomial regression model. For a scalar response ${\displaystyle Y}$ and a functional covariate ${\displaystyle X(\cdot )}$ with domain ${\displaystyle [0,1]}$ and the corresponding centered predictor processes ${\displaystyle X^{c}}$, the simplest and the most prominent member in the family of functional polynomial regression models is the quadratic functional regression^[25] given as follows,$${\displaystyle \mathbb {E} (Y|X)=\alpha +\int _{0}^{1}\beta (t)X^{c}(t)\,dt+\int _{0}^{1}\int _{0}^{1}\gamma (s,t)X^{c}(s)X^{c}(t)\,ds\,dt}$$where ${\displaystyle X^{c}(\cdot )=X(\cdot )-\mathbb {E} (X(\cdot ))}$ is the centered functional covariate, ${\displaystyle \alpha }$ is a scalar coefficient, ${\displaystyle \beta (\cdot )}$ and ${\displaystyle \gamma (\cdot ,\cdot )}$ are coefficient functions with domains ${\displaystyle [0,1]}$ and ${\displaystyle [0,1]\times [0,1]}$, respectively. In addition to the parameter function β that the above functional quadratic regression model shares with the FLM, it also features a parameter surface γ. By analogy to FLMs with scalar responses, estimation of functional polynomial models can be obtained through expanding both the centered covariate ${\displaystyle X^{c}}$ and the coefficient functions ${\displaystyle \beta }$ and ${\displaystyle \gamma }$ in an orthonormal basis.^[25][26]

Functional single and multiple index models

A functional multiple index model is given as below, with symbols having their usual meanings as formerly described,$${\displaystyle \mathbb {E} (Y|X)=g\left(\int _{0}^{1}X^{c}(t)\beta _{1}(t)\,dt,\ldots ,\int _{0}^{1}X^{c}(t)\beta _{p}(t)\,dt\right)}$$Here g represents an (unknown) general smooth function defined on a p-dimensional domain. The case ${\displaystyle p=1}$ yields a functional single index model while multiple index models correspond to the case ${\displaystyle p>1}$. However, for ${\displaystyle p>1}$, this model is problematic due to curse of dimensionality. With ${\displaystyle p>1}$ and relatively small sample sizes, the estimator given by this model often has large variance.^[27][28]

Functional additive models (FAMs)

For a given orthonormal basis ${\displaystyle \{\phi _{k}\}_{k=1}^{\infty }}$ on ${\displaystyle L^{2}[0,1]}$, we can expand ${\displaystyle X^{c}(t)=\sum _{k=1}^{\infty }x_{k}\phi _{k}(t)}$ on the domain ${\displaystyle [0,1]}$.

A functional linear model with scalar responses (see (3)) can thus be written as follows,$${\displaystyle \mathbb {E} (Y|X)=\mathbb {E} (Y)+\sum _{k=1}^{\infty }\beta _{k}x_{k}.}$$One form of FAMs is obtained by replacing the linear function of ${\displaystyle x_{k}}$ in the above expression ( i.e., ${\displaystyle \beta _{k}x_{k}}$) by a general smooth function ${\displaystyle f_{k}}$, analogous to the extension of multiple linear regression models to additive models and is expressed as,$${\displaystyle \mathbb {E} (Y|X)=\mathbb {E} (Y)+\sum _{k=1}^{\infty }f_{k}(x_{k}),}$$where ${\displaystyle f_{k}}$ satisfies ${\displaystyle \mathbb {E} (f_{k}(x_{k}))=0}$ for ${\displaystyle k\in \mathbb {N} }$.^[13][7] This constraint on the general smooth functions ${\displaystyle f_{k}}$ ensures identifiability in the sense that the estimates of these additive component functions do not interfere with that of the intercept term ${\displaystyle \mathbb {E} (Y)}$. Another form of FAM is the continuously additive model,^[29] expressed as,$${\displaystyle \mathbb {E} (Y|X)=\mathbb {E} (Y)+\int _{0}^{1}g(t,X(t))dt}$$for a bivariate smooth additive surface ${\displaystyle g:[0,1]\times \mathbb {R} \longrightarrow \mathbb {R} }$ which is required to satisfy ${\displaystyle \mathbb {E} [g(t,X(t))]=0}$ for all ${\displaystyle t\in [0,1]}$, in order to ensure identifiability.

Generalized functional linear model

An obvious and direct extension of FLMs with scalar responses (see (3)) is to add a link function leading to a generalized functional linear model (GFLM)^[30] in analogy to the generalized linear model (GLM). The three components of the GFLM are:

1. Linear predictor ${\displaystyle \eta =\beta _{0}+\int _{0}^{1}X^{c}(t)\beta (t)\,dt}$; [systematic component]
2. Variance function ${\displaystyle {\text{Var}}(Y|X)=V(\mu )}$, where ${\displaystyle \mu =\mathbb {E} (Y|X)}$ is the conditional mean; [random component]
3. Link function ${\displaystyle g}$ connecting the conditional mean ${\displaystyle \mu }$ and the linear predictor ${\displaystyle \eta }$ through ${\displaystyle \mu =g(\eta )}$. [systematic component]

Motivations

In addition to amplitude variation,^[50] time variation may also be assumed to present in functional data. Time variation occurs when the subject-specific timing of certain events of interest varies among subjects. One classical example is the Berkeley Growth Study Data,^[51] where the amplitude variation is the growth rate and the time variation explains the difference in children's biological age at which the pubertal and the pre-pubertal growth spurt occurred. In the presence of time variation, the cross-sectional mean function may not be an efficient estimate as peaks and troughs are located randomly and thus meaningful signals may be distorted or hidden.

Time warping, also known as curve registration,^[52] curve alignment or time synchronization, aims to identify and separate amplitude variation and time variation. If both time and amplitude variation are present, then the observed functional data ${\displaystyle Y_{i}}$ can be modeled as ${\displaystyle Y_{i}(t)=X_{i}[h_{i}^{-1}(t)],t\in [0,1]}$, where ${\displaystyle X_{i}{\overset {iid}{\sim }}X}$ is a latent amplitude function and ${\displaystyle h_{i}{\overset {iid}{\sim }}h}$ is a latent time warping function that corresponds to a cumulative distribution function. The time warping functions ${\displaystyle h}$ are assumed to be invertible and to satisfy ${\displaystyle \mathbb {E} (h^{-1}(t))=t}$.

The simplest case of a family of warping functions to specify phase variation is linear transformation, that is ${\displaystyle h(t)=\delta +\gamma t}$, which warps the time of an underlying template function by subjected-specific shift and scale. More general class of warping functions includes diffeomorphisms of the domain to itself, that is, loosely speaking, a class of invertible functions that maps the compact domain to itself such that both the function and its inverse are smooth. The set of linear transformation is contained in the set of diffeomorphisms.^[53] One challenge in time warping is identifiability of amplitude and phase variation. Specific assumptions are required to break this non-identifiability.

Methods

Earlier approaches include dynamic time warping (DTW) used for applications such as speech recognition.^[54] Another traditional method for time warping is landmark registration,^[55][56] which aligns special features such as peak locations to an average location. Other relevant warping methods include pairwise warping,^[57] registration using ${\displaystyle {\mathcal {L}}^{2}}$ distance^[53] and elastic warping.^[58]

Multidimensional domain of ${\displaystyle X(\cdot )}$

The domain of ${\displaystyle X(\cdot )}$ can be in ${\displaystyle R^{p}}$, for example the data could be a sample of random surfaces.^[59][60]

Multivariate stochastic process

The range set of the stochastic process may be extended from ${\displaystyle R}$ to ${\displaystyle R^{p}}$^[61][62][63] and further to nonlinear manifolds,^[64] Hilbert spaces^[65] and eventually to metric spaces.^[59]