# Robust Bayesian Functional Principal Component Analysis

Jiarui Zhang
Department of Statistics and Actuarial Science, Simon Fraser University
Jiguo Cao
Department of Statistics and Actuarial Science, Simon Fraser University
Liangliang Wang
Department of Statistics and Actuarial Science, Simon Fraser University

###### Abstract

We develop a robust Bayesian functional principal component analysis (FPCA) by
incorporating skew elliptical classes of distributions. The proposed method
effectively captures the primary source of variation among curves, even when
abnormal observations contaminate the data. We model the observations using
skew elliptical distributions by introducing skewness with transformation and
conditioning into the multivariate elliptical symmetric distribution. To
recast the covariance function, we employ an approximate spectral
decomposition. We discuss the selection of prior specifications and provide
detailed information on posterior inference, including the forms of the full
conditional distributions, choices of hyperparameters, and model selection
strategies. Furthermore, we extend our model to accommodate sparse functional
data with only a few observations per curve, thereby creating a more general
Bayesian framework for FPCA. To assess the performance of our proposed model,
we conduct simulation studies comparing it to well-known frequentist methods
and conventional Bayesian methods. The results demonstrate that our method
outperforms existing approaches in the presence of outliers and performs
competitively in outlier-free datasets. Furthermore, we illustrate the
effectiveness of our method by applying it to environmental and biological
data to identify outlying functional data. The implementation of our proposed
method and applications are available at https://github.com/SFU-Stat-
ML/RBFPCA.

Keywords: Functional data analysis, Robust estimation, Sparse functional data,
Multivariate skew elliptical distribution, Outlier detection

##  1 Introduction

The development of modern technology has resulted in the continuous recording
of data during a given period in many scientific fields, such as neuroscience,
biology, and environmental science. These data can be categorized as
functional data (Ramsay and Silverman, 2005b; Ferraty and Vieu, 2006; Horváth
and Kokoszka, 2012; Hsing and Eubank, 2015), which are usually observed over
time, space or any other continuous domain. For instance, Figure 1 displays
two examples of functional data. Figure 1a shows a collection of densely
observed trajectories from the Hawaii Ocean Oxygen dataset, wherein certain
outlying trajectories have been highlighted. Figure 1b presents one
illustration of sparse function data from the CD4 dataset, with only a few
observations per curve. Both datasets have a few curves demonstrating unusual
patterns in contrast to the rest.

(a)

(b)

Figure 1: (a) Hawaii Ocean Oxygen data contains 133 densely observed
trajectories for oxygen concentrations measuring at different depths below the
sea surface. The two blue trajectories exhibit atypical behaviors in
comparison with the remaining ones. (b) A subset of CD4 data with 10 sparsely
observed trajectories for CD4 cell counts are shown to demonstrate the sparse
nature. The two trajectories colored in blue display abnormal patterns.

As the number of functional data can grow rapidly, data compression is vital
for functional data. Functional principal component analysis (FPCA) plays a
crucial role in functional data analysis (FDA) as a dimension reduction
technique. FPCA can be viewed as an exploratory analysis tool to discover the
hidden structure of the data by capturing the optimal low-dimensional
representation and the major source of variations among curves. FPCA can also
be used to help reconstruct partially-observed functions or used as a
preprocessing step for regression analysis (Reiss et al., 2017) or clustering
tasks (Margaritella et al., 2021). The frequentist analysis of FPCA is a
mature field. Numerous works have been proposed to explore fully or densely
observed function data (Dauxois et al., 1982; Rice and Silverman, 1991;
Cardot, 2000; Hall and Hosseini-Nasab, 2006). Studies of the FPCA approach for
sparsely observed data, a more challenging situation, can also be seen in
James et al. (2000), Yao et al. (2005), and Paul and Peng (2009). Wang et al.
(2016) provided a detailed review of methods, open questions and applications
of FPCA.

In contrast, there has been limited investigation into the Bayesian viewpoint
of FPCA. The lack of Bayesian methodologies in FDA and FPCA may be due to the
challenges involved in specifying the full probability model and the limited
availability of inferential tools for implementation. Nonetheless, the
Bayesian framework is appealing from at least three perspectives. First, the
Bayesian methods are capable of quantifying uncertainty in a straightforward
manner, such as through credible intervals. Second, the Bayesian inferential
structure is flexible as it can incorporate subject matter expertise in the
definition of prior structures. Last but not least, the Bayesian framework
also provides a straightforward approach for conducting model selection by
marginal likelihood estimates.

One example of applying Bayesian methodologies in FPCA is presented by Behseta
et al. (2005), who utilize a hierarchical Gaussian process model to assess
variability among functions. Later, Van Der Linde (2008) proposed analyzing
the modes of variation of curves via variational inference. The prevalent
strategy has been to work within the hierarchical representation of the FPCA
model and specify the prior distributions for all model parameters
(Crainiceanu and Goldsmith, 2010; Margaritella et al., 2021). Alternatively,
Suarez and Ghosal (2017) proposed a Bayesian FPCA method using approximate
spectral decomposition and modelled the number of principal components with
truncated Poisson distributions. For partially observed functional data, Jiang
et al. (2020) considered a Bayesian model for sparse FPCA using a reduced rank
mixed-effects framework. To the best of our knowledge, no existing literature
has explored outlier detection within the context of a robust Bayesian FPCA
framework that is capable of handling both dense and sparse functional data.

Most previous works on the Bayesian FPCA model employed multivariate normal
distribution to model the discretized, noise-contaminated observations.
Despite the fact that the multivariate normal distribution has several
desirable properties for modelling data, it is common that observations are
not normally distributed in general. Various alternative distributions are
available to accommodate higher moments. For instance, the multivariate
Student-$t$ distribution works well for fat-tailed data but does not consider
asymmetry. The log-normal distribution has been used to model skewed data, but
its skewness is reflected as a function of the mean and variance. This is
somewhat inflexible because skewness is not a separate parameter. Branco and
Dey (2001) proposed a general class of skew elliptical distributions by
extending the previous work on a multivariate skew normal distribution as
described by Azzalini and Valle (1996). The general class developed by Branco
and Dey (2001) includes the multivariate normal, Student-$t$, and exponential
power, but with an additional parameter to control skewness. This family of
distributions was improved upon by Sahu et al. (2003) for capturing
multivariate asymmetry and adding more flexibility in adjusting the
correlation structure, which often leads to better-fitting models. Moreover,
the family of distributions proposed by Sahu et al. (2003) offers a convenient
implementation within a Bayesian framework.

Given the various advantages offered by Bayesian inference and the advancement
of the skew elliptical class of distributions, we propose a general method to
fit the observations with some multivariate skew elliptical distributions
within the Bayesian FPCA framework. Our proposed method has three main
advantages. First, our method produces more robust estimations of the
covariance function and the corresponding principal components by integrating
the class of skew elliptical distributions. In practice, many problems involve
a proper transformation for symmetry for skewed data. Such ad hoc
transformations can be avoided by using the proposed models. Rather than
transforming the data, our methods transform error distributions to handle
skewness. Second, adopting a Bayesian approach for the FPCA problem
facilitates the interpretation of uncertainty measures, such as credible
intervals. This aligns with the inherent view of making a probability
statement concerning the data after people have observed it. Finally, we
extend our method to fit functional data that are sparsely and irregularly
observed, making the Bayesian FPCA fully adapted to any type of functional
data.

The rest of the article is organized as follows: We introduce the setup of
Bayesian FPCA in Section 2.1. The proposed model is presented in Section 2.2,
along with computational details in Section 2.3, a description of outlier
detection method in Section 2.4, and an extension to sparse data in Section
2.5. Simulation studies are presented in Section 3, in which we compare the
proposed model to other Bayesian and frequentist FPCA approaches in the
presence of different data generation processes and noise levels. Three
examples including dense and sparse data are discussed in Section 4. A summary
and some directions for future work are given in Section 5.

##  2 Robust Bayesian FPCA

In this section, we will present the proposed robust Bayesian FPCA model and
describe the features that contribute to its improvement over the existing
FPCA methods.

###  2.1 Bayesian Functional Principal Component Analysis

Functional principal component analysis (FPCA) extends the idea of the
principal component analysis of multivariate data (see Jolliffe, 2002, for a
comprehensive introduction to classical multivariate PCA) to functional data.
FPCA is the most prevalent tool in FDA and serves as a standard first step
when processing functional data in most cases (Ramsay and Silverman, 2005a).
This is due to FPCA’s capability to convert inherently infinite-dimensional
functional data into a finite-dimensional vector of random scores. The general
goals of FPCA include finding patterns in data of high dimensions, capturing
the major sources of variation and reducing the dimensionality.

Denote $\\{X(t):t\in[0,T]\in\mathcal{R}\\}$ as a square integrable stochastic
process. Consider $n$ independent and identically distributed realizations,
$X_{1}(t),X_{2}(t),\ldots,X_{n}(t)$, of a mean-zero stochastic process with
covariance function $\text{Cov}_{X}(s,t)$ at a sequence of random points on
$\mathcal{T}=[0,T]$. The covariance function $\text{Cov}_{X}(s,t)$ specifies
the covariance between curve values $X_{i}(s)$ and $X_{i}(t)$ at times $s$ and
$t$, respectively. We assume the observations are noise-corrupted. That is,
the observed data $Y_{i}^{*}(t_{ij})$, $i=1,\ldots,n$, $j=1,\ldots,n_{i}$, is

$Y_{i}^{*}(t_{ij})=\mu(t)+X_{i}(t_{ij})+\epsilon_{ij},$ (1)

where $Y_{i}^{*}(t_{ij})$ is the discretized, noise-contaminated, observed
data for every curve $i=1,\ldots,n,$ at the time point $t_{ij}$; $\mu(t)$ is
the mean function; $X_{i}(t_{ij})$ is the true underlying detrended curve from
a Gaussian process with mean zero and covariance function, i.e.,
$X_{i}(t)\sim\text{GP}(\textbf{0},\text{Cov}_{X}(s,t))$ for
$s,t\in\mathcal{T}$, independently for $i=1,\ldots,n$; and $\epsilon_{ij}$’s
are the noise terms with $\mathbb{E}(\epsilon_{ij})=0$ and
$\text{Var}(\epsilon_{ij})=\sigma^{2}_{ij}$. We assume $\epsilon_{ij}$ are
independent across $i$ and $j$.

We first consider the case when the observed data $Y_{i}^{*}(t_{ij})$ are
dense and equally spaced, i.e., the number of measurements $n_{i}$ for each
curve is the same and the sequence $t_{i1},\ldots,t_{in_{i}}$ is equally
spaced for every curve $i=1,\ldots,n$. Let
$Y_{i}(t)=Y_{i}^{*}(t)-\hat{\mu}(t)$ be the detrended observations, where
$\hat{\mu}(t)=\frac{1}{n}\sum_{i=1}^{n}Y_{i}^{*}(t)$.

Mercer’s theorem states that, under mild assumptions, the spectral
decomposition of the covariance function can be approximated as

$\text{Cov}_{X}^{K}(s,t)=\sum^{K}_{k=1}\lambda_{k}\phi_{k}(s)\phi_{k}(t),\;\text{for}\>s,t\in\mathcal{T},$
(2)

where $\text{Cov}_{X}^{K}(s,t)$ converges uniformly to $\text{Cov}_{X}(s,t)$
as $K\rightarrow\infty$, $\\{\lambda_{1},\lambda_{2},\dots\\}$ are the non-
negative eigenvalues in non-increasing order of the covariance operator, and
$\\{\phi_{1},\phi_{2},\dots\\}$ are the corresponding orthogonal
eigenfunctions. We can simplify the notation for covariance function by
defining $Q(s,t)=\text{Cov}_{X}^{K}(s,t)$ throughout the rest of the article.
By the Karhunen-Loève expansion (Fukunaga and Koontz, 1970), each de-meaned
curve $X_{i}(t)$ can be approximated by a linear combination of the
eigenfunctions $\phi_{k}(t)$ and FPC scores $\xi_{ik}$:

$X_{i}(t)\approx\sum^{K}_{k=1}\xi_{ik}\phi_{k}(t).$ (3)

By truncating the infinity expansion at some level, $K\in\mathbb{N}$, the
resulting finitely truncated series provides a good approximation to
$X_{i}(t)$.

The eigenfunctions $\\{\phi_{1},\phi_{2},\dots\\}$ can be further represented
with the following basis expansion:

$\phi_{k}(t)=\sum^{\infty}_{p=1}c_{kp}h_{p}(t)\approx\sum^{P}_{p=1}c_{kp}h_{p}(t),$
(4)

where $h_{p}(t)$’s are some given basis functions and $c_{kp}$’s are the
corresponding basis coefficients. The expansion in equation (4) is truncated
at some level, $P\in\mathbb{N}$. Assigning priors on the truncating parameters
$P$ and $K$ ensure full support. Let $\boldsymbol{H}_{P}$ be the $T\times P$
matrix with columns consisting of the given basis functions evaluated at all
points. Let $\boldsymbol{C}_{KP}$ denote the coefficient matrix of size
$K\times P$ and $\boldsymbol{\Lambda}_{K}$ be a $K\times K$ diagonal matrix
whose diagonal is composed of the non-negative eigenvalues
$\\{\lambda_{1},\lambda_{2},\dots,\lambda_{K}\\}$ in non-increasing order.
Then the covariance function can be represented by the following relation:

$Q(s,t)=\boldsymbol{H}^{\prime}_{P}(s)\boldsymbol{C}^{\prime}_{KP}\boldsymbol{\Lambda}_{K}\boldsymbol{C}_{KP}\boldsymbol{H}_{P}(t),\;\text{for}\>s,t\in\mathcal{T}.$
(5)

###  2.2 Models and Prior Specification

In this section, we give details on the structure of the proposed robust
Bayesian FPCA (RB-FPCA) method and features that improve the previous Bayesian
FPCA model in Suarez and Ghosal (2017). The prior on the covariance function
$Q(s,t)$ in Equation 5 will be specified indirectly on
$\boldsymbol{\Omega}^{-1}=(\boldsymbol{C}^{\prime}_{KP}\boldsymbol{\Lambda}_{K}\boldsymbol{C}_{KP})^{-1}$,
and the prior is assumed to be a Wishart distribution. When the number of
principal components $(K)$ is less than the number of basis functions $(P)$
used for approximation, the Wishart matrix becomes singular (Uhlig, 1994).
Because the Wishart distribution gives zero density to singular matrices, to
allow the constraint of $K\leq P$, it is necessary to consider the inclusion
of singular Wishart matrices into the specification of the prior.

Let $\boldsymbol{\Psi}=\boldsymbol{U}\boldsymbol{L}\boldsymbol{U}^{\prime}$,
where $\boldsymbol{U}$ is a $P\times P$ orthogonal matrix and $\boldsymbol{L}$
is a $P\times P$ diagonal matrix with ordered eigenvalues along the diagonal.
Now define
$\boldsymbol{\Psi}_{K}=\boldsymbol{U}_{K}\boldsymbol{L}_{K}\boldsymbol{U}_{K}^{\prime}$,
where $\boldsymbol{U}_{K}$ is the first $K$ columns of $\boldsymbol{U}$ and
$\boldsymbol{L}_{K}$ is the first $K$ rows and $K$ columns of
$\boldsymbol{L}$. It follows that $\boldsymbol{\Psi}_{K}$ has rank $K$ even
though its dimension is $P\times P$. By defining the singular center matrix
$\boldsymbol{\Psi}_{K}$ of rank $K$, we allow for singular Wishart matrices in
specifying the prior (see Suarez and Ghosal (2017)). For any random matrix
$\boldsymbol{S}$, a Moore-Penrose inverse of $\boldsymbol{S}$ is defined as a
matrix $\boldsymbol{S}^{+}$. Suppose that $\boldsymbol{S}_{1}\sim
Wishart_{P}(\nu,\boldsymbol{\Psi}_{K}^{+})$ and $\boldsymbol{S}_{2}\sim
Wishart_{K}(\nu,\boldsymbol{L}_{k}^{-1})$, then

$\displaystyle\boldsymbol{S}_{1}$
$\displaystyle\xrightarrow{D}\boldsymbol{U}_{K}\boldsymbol{S}_{2}\boldsymbol{U}_{K}^{\prime},$
$\displaystyle\boldsymbol{S}_{1}^{+}$
$\displaystyle\xrightarrow{D}\boldsymbol{U}_{K}\boldsymbol{S}_{2}^{-1}\boldsymbol{U}_{K}^{\prime}.$

This property forms the desired singularity in the prior specification.

The model proposed in this work differs from that of Suarez and Ghosal (2017)
who model the discretized observations with multivariate normal distribution.
We propose to fit multivariate skew elliptical distributions to model the
observations, which produces more robust estimates in the Bayesian FPCA
framework. To construct the skew elliptical class of distributions, we first
present the general settings for multivariate symmetric elliptical
distributions, followed by an introduction to the general skew elliptical
class of distributions, where skewness is induced via transformation and
conditioning.

Let $\boldsymbol{\epsilon}$ and $\boldsymbol{z}$ be two $T$-dimensional random
vectors. Let $\boldsymbol{\mu}$ denote an $T$-dimensional vector and
$\boldsymbol{\Sigma}$ denote an $T\times T$ positive definite matrix. We will
use the following notation to represent the general class of elliptically
symmetric distributions:

$\boldsymbol{\eta}=\left(\begin{array}[]{c}\boldsymbol{\epsilon}\\\
\boldsymbol{z}\end{array}\right)\sim
El\left(\left(\begin{array}[]{c}\boldsymbol{\mu}\\\
\boldsymbol{0}\end{array}\right),\begin{bmatrix}\boldsymbol{\Sigma}&\boldsymbol{0}\\\
\boldsymbol{0}&\boldsymbol{I}\end{bmatrix};\>g^{(2T)}\right),$ (6)

where $\boldsymbol{0}$ and $\boldsymbol{I}$ are the null matrix and identity
matrix, respectively, and $g^{(2T)}$ is the generator of the probability
density function (PDF). Note that the generator $g^{(2T)}$ depends on the
dimensions of the random vector
$\boldsymbol{\eta}=(\boldsymbol{\epsilon}^{\prime},\boldsymbol{z}^{\prime})^{\prime}$,
and $2T$ is the dimension of $\boldsymbol{\eta}$. The density generator
$g^{(2T)}$ is a function from $\mathbb{R}^{+}$ to $\mathbb{R}^{+}$ which
satisfies

$\int_{o}^{\infty}w^{T-1}g^{(2T)}(w)\,dw=\pi^{-T}\Gamma\left(T\right),$ (7)

where $\Gamma(T)=(T-1)!$ is the gamma function. The choice of the density
generator function will determine the distribution of the random vector. This
flexible nature of the elliptical class of distributions allows for including
several widely recognized symmetric distributions as proper members, for
example, the multivariate normal distribution. Assume a random vector
$\boldsymbol{\eta}$ of dimension $2T$ follows a multivariate Normal
distribution with mean vector $\boldsymbol{\mu}$ and covariance matrix
$\boldsymbol{\Sigma}$, i.e., $\boldsymbol{\eta}\sim
MVN(\boldsymbol{\mu},\boldsymbol{\Sigma})$, the probability density function
can be formulated via the following generalized expression in terms of the
generator function $g^{(2T)}$:

$f(\boldsymbol{\eta}\>|\>\boldsymbol{\mu},\boldsymbol{\Sigma};g^{(2T)})=|\boldsymbol{\Sigma}|^{-1/2}g^{(2T)}((\boldsymbol{\eta}-\boldsymbol{\mu})^{\prime}\boldsymbol{\Sigma}^{-1}(\boldsymbol{\eta}-\boldsymbol{\mu}))$,

where the density generator function has the form of
$g^{(2T)}(u)=e^{-u/2}/(2\pi)^{T}$. The conventional form of the multivariate
normal distribution’s PDF can be retrieved with this density generator
function.

Next, we consider a general skew elliptical class of distributions by
implementing the transformation presented in Sahu et al. (2003):

$\boldsymbol{Y}=\boldsymbol{D}\boldsymbol{z}+\boldsymbol{\epsilon},$ (8)

where $\boldsymbol{Y}$ is the transformed variable which followed the skew
elliptical distribution, and $\boldsymbol{D}$ is a diagonal matrix with
diagonals $\boldsymbol{d}=(d_{1},\ldots,d_{T})^{\prime}$, which accommodates
skewness. Let $\text{vec}(\boldsymbol{D})$ denote a vector of diagonals of
$\boldsymbol{D}$. The skew elliptical class is developed by considering the
random variable $\boldsymbol{Y}\,|\,\boldsymbol{z}\\!>\\!0$, where
$\boldsymbol{z}=(z_{1},\ldots,z_{T})^{\prime}>0$ implies every element of
$\boldsymbol{z}$ is positive. The construction in equation (8) along with the
conditioning introduces skewness. Specifically, positive values of elements of
$\boldsymbol{d}$ result in positively skewed distributions, while negative
values lead to negatively skewed distributions. Note that the elliptically
symmetric distribution is retrieved if the diagonals of $\boldsymbol{D}$ are
zeros, i.e., $d_{1},\ldots,d_{T}=0$. With the transformation and conditioning,
the random variable $\boldsymbol{Y}$ follows the skew elliptical distribution
and we denote it with the notation $\boldsymbol{Y}\sim
SE(\boldsymbol{\mu},\boldsymbol{\Sigma},\boldsymbol{D};g^{(T)})$. We will
consider one example of the skew elliptical class of distributions throughout
this paper, namely the multivariate skew normal distribution.

We can now integrate the family of skew elliptical distributions within the
Bayesian FPCA framework and present the full RB-FPCA model. Let
$\boldsymbol{Y}_{i}=(Y_{i}(t))^{\prime}$ for $i=1,\ldots,n$. We assume
$\boldsymbol{Y}_{i}\sim
SE(\boldsymbol{\mu}_{i},\boldsymbol{\Sigma},\boldsymbol{D};g^{(n_{i})})$, and
further relax the assumption of homoscedasticity on noises. The hierarchical
representation of the proposed RB-FPCA model for $i=1,\ldots,n,$ is

$\displaystyle\boldsymbol{Y}_{i}\>|\>\cdot$ $\displaystyle\sim
MVN_{T}(\boldsymbol{H}_{P}\boldsymbol{U}_{K}\boldsymbol{\beta}_{i,K}+\boldsymbol{D}\boldsymbol{z}_{i},\>\boldsymbol{\Sigma}),$
$\displaystyle\boldsymbol{\beta}_{i,K}$ $\displaystyle\stackrel{{\scriptstyle
i.i.d}}{{\sim}}MVN_{K}(\boldsymbol{0},\>\boldsymbol{\Omega}),$
$\displaystyle\boldsymbol{\Omega}^{-1}$ $\displaystyle\sim
Wishart_{K}(\nu,\>\boldsymbol{L}_{k}^{-1}),$ $\displaystyle\boldsymbol{z}_{i}$
$\displaystyle\sim
MVN_{T}(\boldsymbol{0},\boldsymbol{I})I(\boldsymbol{z}_{i}>0),$
$\displaystyle\text{vec}(\boldsymbol{D})$ $\displaystyle\sim
MVN_{T}(\boldsymbol{0},\boldsymbol{\Gamma}),$
$\displaystyle\boldsymbol{\Sigma}^{-1}$ $\displaystyle\sim
Wishart_{T}(2r,2\boldsymbol{\kappa}),$

where $\nu,\boldsymbol{\Gamma},2r$ and $\boldsymbol{\kappa}$ are
hyperparameters. This proposed model takes singular Wishart matrices into
account in the prior specification because
$\boldsymbol{U}_{K}\boldsymbol{\Omega}^{-1}\boldsymbol{U}_{K}^{\prime}\sim
Wishart_{P}(\nu,\boldsymbol{\Psi}_{K}^{+})$. The main purpose of the inference
is to estimate the covariance function and the principal components from the
data. More specifically, the covariance function of interest is
$\boldsymbol{H}^{\prime}_{P}(s)\boldsymbol{U}_{K}^{\prime}\boldsymbol{\Omega}\boldsymbol{U}_{K}\boldsymbol{H}_{P}(t)$
and the principal components it induces are
$\boldsymbol{C}_{KP}\boldsymbol{H}^{\prime}_{P}(t)$.

###  2.3 Posterior Inference

Our goal is to sample the parameters from the posterior distribution by
deriving their full conditional distributions. All parameters’ conditional
distributions can be derived analytically in closed forms and are given as
follows. In each conditional distribution, we use $|\cdots$ to denote
conditioning on the data and all other parameters.

  * •

For $i=1,\ldots,n$, the full conditional distribution for
$\boldsymbol{\beta}_{i,K}$ is

$\begin{split}\boldsymbol{\beta}_{i,K}|\cdots\sim
MVN_{K}((\boldsymbol{U}_{K}^{\prime}\boldsymbol{H}_{P}^{\prime}&\boldsymbol{\Sigma}^{-1}\boldsymbol{H}_{P}\boldsymbol{U}_{K}+\boldsymbol{\Omega}^{-1})^{-1}(\boldsymbol{U}_{K}^{\prime}\boldsymbol{H}_{P}^{\prime}\boldsymbol{\Sigma}^{-1}\boldsymbol{Y}_{i}),\\\
&(\boldsymbol{U}_{K}^{\prime}\boldsymbol{H}_{P}^{\prime}\boldsymbol{\Sigma}^{-1}\boldsymbol{H}_{P}\boldsymbol{U}_{K}+\boldsymbol{\Omega}^{-1})^{-1}),\end{split}$
(9)

where $\boldsymbol{H}_{P}$ is picked to be the Legendre polynomial basis
functions and the choice of $\boldsymbol{U}_{K}$ is discussed in the
conditional distribution for $\boldsymbol{\Omega}^{-1}$.

  * •

For $\boldsymbol{\Omega}^{-1}$, the full conditional distribution is

$\boldsymbol{\Omega}^{-1}|\cdots\sim
Wishart_{K}(\nu+n+1,(\boldsymbol{L}_{k}+\sum^{n}_{i=1}\boldsymbol{\beta}_{i,K}\boldsymbol{\beta}_{i,K}^{\prime})^{-1}),$
(10)

where $\nu$ is the number of degrees of freedom. A sensible choice of
$\boldsymbol{\Psi}$ is given as
$\boldsymbol{\Psi}=(\boldsymbol{H}_{P}^{\prime}\boldsymbol{H}_{P})^{-1}\boldsymbol{H}_{P}^{\prime}\boldsymbol{\Omega}^{*}\boldsymbol{H}_{P}(\boldsymbol{H}_{P}^{\prime}\boldsymbol{H}_{P})^{-1}$,
where $\boldsymbol{\Omega}^{*}$ is the prior covariance function corresponding
to the time grid being used. $\boldsymbol{U}_{K}$ and $\boldsymbol{L}_{K}$ are
obtained by decomposing $\boldsymbol{\Psi}$ and subsetting the resulting
matrices $\boldsymbol{U}$ and $\boldsymbol{L}$ as described in Section 2.2.

  * •

For $\boldsymbol{z}_{i}$, the full conditional distribution is

$\boldsymbol{z}_{i}|\cdots\sim
MVN_{T}(\boldsymbol{A}_{i}^{-1}\boldsymbol{a}_{i},\boldsymbol{A}_{i}^{-1})I(\boldsymbol{z}_{i}>0),$
(11)

where
$\boldsymbol{A}_{i}=\boldsymbol{I}+\boldsymbol{D}\boldsymbol{\Sigma}^{-1}\boldsymbol{D}$
and
$\boldsymbol{a}_{i}=\boldsymbol{D}\boldsymbol{\Sigma}^{-1}(\boldsymbol{Y}_{i}-\boldsymbol{H}_{P}\boldsymbol{U}_{K}\boldsymbol{\beta}_{i,K})$.

  * •

For $\text{vec}(\boldsymbol{D})$, the full conditional distribution is

$\text{vec}(\boldsymbol{D})|\cdots\sim
MVN_{T}(\boldsymbol{B}^{-1}\boldsymbol{b},\boldsymbol{B}^{-1}),$ (12)

with
$\boldsymbol{B}=\boldsymbol{\Gamma}^{-1}+\sum^{n}_{i=1}\text{diag}(\boldsymbol{z}_{i})\boldsymbol{\Sigma}^{-1}\text{diag}(\boldsymbol{z}_{i})$,
$\boldsymbol{b}=\sum^{n}_{i=1}\text{diag}(\boldsymbol{z}_{i})\boldsymbol{\Sigma}^{-1}(\boldsymbol{Y}_{i}-\boldsymbol{H}_{P}\boldsymbol{U}_{K}\boldsymbol{\beta}_{i,K})$.

  * •

For $\boldsymbol{\Sigma}^{-1}$, the full conditional distribution is

$\begin{split}\boldsymbol{\Sigma}^{-1}|\cdots\sim&Wishart_{T}(2r+n,((2\boldsymbol{\kappa})^{-1}+\\\
&\sum^{n}_{i=1}(\boldsymbol{Y}_{i}-\boldsymbol{H}_{P}\boldsymbol{U}_{K}\boldsymbol{\beta}_{i,K}-\boldsymbol{D}\boldsymbol{z}_{i})(\boldsymbol{Y}_{i}-\boldsymbol{H}_{P}\boldsymbol{U}_{K}\boldsymbol{\beta}_{i,K}-\boldsymbol{D}\boldsymbol{z}_{i})^{\prime})^{-1}.\end{split}$
(13)

The conditional distributions for all parameters of interest are readily
accessible in closed forms, we therefore employ a Gibbs sampling algorithm to
sample from the posterior distribution. The algorithm samples new values for
each parameter iteratively based on the current values of the other unknown
parameters. Some advantages of the Gibbs sampling algorithm include its
simplicity and fast convergence in our implementation. For posterior
inference, the covariance function is estimated with the full posterior mean.
The principal components $\hat{\phi}_{k}(t)$ are estimated by decomposing the
estimated posterior mean of the covariance function. The corresponding FPC
scores are estimated with $\hat{\xi}_{ik}=$
$\int_{\mathcal{T}}\hat{\phi}_{k}(t)Y_{i}(t)\,dt$.

Next, we will discuss possible approaches to determine the number of
eigenfunctions $K$, which gives a reasonably good approximation to the
infinite-dimensional functional data. Model comparison through Bayes factors
(Jeffreys, 1935; Han and Carlin, 2001) is a standard way for selecting models
in the Bayesian framework. The calculation of the Bayes factor requires
estimating the marginal likelihood, which is the normalization factor of the
posterior density. There has been great difficulty in obtaining the marginal
likelihood estimates, which poses a challenge to the computation of the Bayes
factor. We employed the method by Chib (1995), and the marginal likelihood
estimates can be easily obtained since all conditional distributions in the
Gibbs sampling algorithm have closed-form expressions. For more efficient
computing, another straightforward approach is to select the number of
eigenfunctions $K$ which can explain a sufficiently large part of the total
variation. We predetermined the value of $K$ in all simulation studies and
data analysis, and results showed that the first four eigenfunctions are able
to explain at least 95% of the total variability in the data in all cases.

###  2.4 Outlier Detection

The estimated FPC scores can be used to detect any potential outlying
trajectories. We employ a robust outlier detection approach described in
Boente and Salibián-Barrera (2021). The first step is calculating the robust
estimates for multivariate location and scatter of the estimated FPC scores.
These estimates, named MM-estimates proposed by Yohai (1987), have a high
breakdown point while remaining efficient. The breakdown point is a measure of
robustness, defined as the maximum proportion of contamination or atypical
points the data may contain while the estimator remains informative about the
underlying parameter. Thus, estimators with higher breakdown-point are
desired. Then we compute the robust Mahalanobis distance, which is defined as

$d_{\text{Mahalanobis}}(\hat{\boldsymbol{\xi}}_{i\cdot})=\sqrt{(\hat{\boldsymbol{\xi}}_{i\cdot}-\hat{\boldsymbol{\mu}}_{MM})^{\prime}\hat{\boldsymbol{\Sigma}}_{MM}^{-1}(\hat{\boldsymbol{\xi}}_{i\cdot}-\hat{\boldsymbol{\mu}}_{MM})}$,

where
$\hat{\boldsymbol{\xi}}_{i\cdot}=(\hat{\xi}_{i1},\ldots,\hat{\xi}_{iK})^{\prime}$,
$\hat{\boldsymbol{\mu}}_{MM}$ and $\hat{\boldsymbol{\Sigma}}_{MM}$ denote the
MM-estimators for multivariate location and scatter, respectively. The
trajectory $i$ is flagged as an outlier if its FPC scores have a distance
larger than a threshold quantile of a $\chi^{2}_{K}$ distribution. Different
threshold values can be applied to detect mild and extreme outliers.

###  2.5 Extend to Sparse Longitudinal Data

In this section, we extend the robust Bayesian FPCA model in 2.2 to fit
longitudinal data that are sparsely and irregularly observed. We assume the
observations are contaminated by measurement errors and adopt the FDA model
given in equation 1. Now we consider the case when the number of measurements
$n_{i}$ made per subject is random due to the sparse and irregular designs. We
also assume that the random variables $n_{i}$ are i.i.d. The sparse data are
centered before fitting the RB-FPCA model, i.e.,
$Y_{i}(t_{ij})=Y_{i}^{*}(t_{ij})-\hat{\mu}(t)$ for $i=1,\ldots,n$,
$j=1,\ldots,n_{i}$. The mean function $\mu(t)$ can be estimated based on the
pooled data from all individuals.

Write
$\boldsymbol{Y}_{i\cdot}=(Y_{i}(t_{i1}),\ldots,Y_{i}(t_{in_{i}}))^{\prime}$.
The hierarchical representation of the proposed RB-FPCA model for
$i=1,\ldots,n$ with sparse longitudinal data is

$\displaystyle\boldsymbol{Y}_{i\cdot}\>|\>\cdot$ $\displaystyle\sim
MVN_{n_{i}}(\boldsymbol{H}_{P}^{(i)}\boldsymbol{U}_{K}^{(i)}\boldsymbol{\beta}_{i,K}+\boldsymbol{D}^{(i)}\boldsymbol{z}^{(i)},\>\boldsymbol{\Sigma}^{(i)}),$
$\displaystyle\boldsymbol{\beta}_{i,K}$ $\displaystyle\stackrel{{\scriptstyle
i.i.d}}{{\sim}}MVN_{K}(\boldsymbol{0},\>\boldsymbol{\Omega}),$
$\displaystyle\boldsymbol{\Omega}^{-1}$ $\displaystyle\sim
Wishart_{K}(\nu,\>\boldsymbol{L}_{k}^{-1}),$
$\displaystyle\boldsymbol{z}^{(i)}$ $\displaystyle\sim
MVN_{n_{i}}(\boldsymbol{0},\boldsymbol{I})I(\boldsymbol{z}^{(i)}>0),$
$\displaystyle\text{vec}(\boldsymbol{D}^{(i)})$ $\displaystyle\sim
MVN_{n_{i}}(\boldsymbol{0},\boldsymbol{\Gamma}^{(i)}),$
$\displaystyle(\boldsymbol{\Sigma}^{(i)})^{-1}$ $\displaystyle\sim
Wishart_{n_{i}}(2r,2\boldsymbol{\kappa}^{(i)}),$

where the superscript $(i)$ represents the parameters whose values and
dimensions depend on the curve index $i$. The full conditional distributions
in the Gibbs sampling for the sparse longitudinal data have similar forms
except for the hyperparameter $\boldsymbol{R}^{(i)}$ in the conditional
distribution of $\boldsymbol{\Sigma}^{-1}$, which is defined as a $n_{i}\times
n_{i}$ diagonal matrix where elements of the main diagonal are the squared
ranges of the corresponding components in the data. For the sparse
longitudinal data, the squared ranges of the components cannot be directly
calculated since data are observed at different time points. We estimate the
squared ranges with a search bandwidth denoted as $h_{R}$, and a default value
$h_{R}=0.05n$ is used in all simulations and examples. Specifically, a set of
$h_{R}$ observations that are closest to the target observation will be used
to estimate the squared ranges for the target individual at a specific time
point. Furthermore, to ensure the matrix $\boldsymbol{\Psi}$ described in the
full conditional distribution of $\boldsymbol{\Omega}^{-1}$ is valid, we
assume the number of principal components is less than or equal to the number
of basis functions, and the number of basis functions used for approximation
is less than or equal to the minimum number of observations per subject, i.e.,
$K\leq P\leq\min(1,\ldots,n_{i})$.

The posterior covariance function is
$\boldsymbol{H}^{\prime}_{P}(s)\boldsymbol{U}_{K}^{\prime}\boldsymbol{\Omega}\boldsymbol{U}_{K}\boldsymbol{H}_{P}(t)$,
where $\boldsymbol{H}_{P}$ and $\boldsymbol{U}_{K}$ are constructed using a
given number of support points in each direction of the covariance surface.
The default value of support points is set as 51 (same as in fdapace package
in R). The principal components are
$\boldsymbol{C}_{KP}\boldsymbol{H}^{\prime}_{P}(t)$. The traditional way to
estimate the FPC scores with numerical integration is used and performs well
when a sufficiently dense measurement grid for each subject is available.
However, when dealing with sparse longitudinal data, the numerical integration
approximated by sums could not produce reasonable approximations due to the
sparseness of the data. Thus, the alternative PACE method described by Yao et
al. (2005) is implemented here to find FPC scores for sparse longitudinal
data. Specifically, using the data from each subject, the prediction of the
FPC scores for the $i$th subject is given by the conditional expectation:

$\hat{\xi}_{ik}=\hat{\mathbb{E}}[\xi_{ik}|\boldsymbol{Y}_{i\cdot}]=\hat{\lambda}_{k}\hat{\phi}^{\prime}_{ik}(t)\hat{Q}^{-1}_{\boldsymbol{Y}_{i\cdot}}\boldsymbol{Y}_{i\cdot},$
(14)

where $\hat{\lambda}_{k}$ is the $k$th estimated eigenvalues of
$\hat{Q}(s,t)$, $\hat{\phi}_{ik}(t)$ is the estimates of the eigenfunctions
$\phi_{ik}(t)=(\phi_{k}(t_{i1}),\ldots,\phi_{k}(t_{in_{i}}))^{\prime}$, and
the $(j,l)$th element of $\hat{Q}_{\boldsymbol{Y}_{i\cdot}}$ is
$(\hat{Q}_{\boldsymbol{Y}_{i\cdot}})_{j,l}=\hat{Q}(t_{ij},t_{il})+\boldsymbol{\hat{\Sigma}}^{(i)}$.
This conditioning method works with sparse data with the existence of
measurement errors and gives the best prediction of the FPC scores under
Gaussian assumptions. The estimate $\hat{\xi}_{ik}$ in equation 14 is the best
linear prediction of $\xi_{ik}$ from the information in the $i$th subject,
regardless of whether the Gaussian assumption holds (as previously observed in
Yao et al. (2005)).

##  3 Simulation Studies

###  3.1 Simulation I: Dense functional data without outlisers

In the first simulation, we aim to compare the performance of the proposed RB-
FPCA method with different frequentist and Bayesian FPCA methods. We followed
a similar simulation setting as in Suarez and Ghosal (2017). Each set of
simulated data consists of 20 noisy observations over 50 time points. These
time points are evenly spaced in the interval $[-1,1]$. The underlying true
mean function is $\mu(t)=\text{sin}(2\pi t)$. Depending on the experimental
settings, the true underlying covariance function is of either
$\text{Cov}_{1}(s,t)=\exp\\{-3(t-s)^{2}\\}$ or
$\text{Cov}_{2}(s,t)=\text{min}\\{s+1,t+1\\}$. Next, independent noise is
added to the sampling process in the form of $N(0,0.3)$ random variables.

The four methods for comparison are listed below:

  1. i.

our proposed robust Bayesian FPCA (RB-FPCA) method

  2. ii.

Bayesian FPCA (BFPCA) method of Suarez and Ghosal (2017)

  3. iii.

Fast Covariance Estimation (FACE) method of Xiao et al. (2016), implemented in
refund package in R

  4. iv.

Principal Analysis by Conditional Estimation (PACE) method of Wang et al.
(2016), implemented in fdapace package in R

In RB-FPCA and BFPCA methods, we predetermined the numbers of basis functions
and eigenfunctions as $P=15$ and $K=5$. We established the values for the
hyperparameters as outlined in Sahu et al. (2003). We set $\nu=2K$ for
$\boldsymbol{\Omega}^{-1}$ in equation 10,
$\boldsymbol{\Gamma}=\text{diag}(10,\dots,10)$ for
$\text{vec}(\boldsymbol{D})$ in equation 12, $2r=T$ and
$\boldsymbol{\kappa}=100*\boldsymbol{R}^{-1}/(2r)$ for
$\boldsymbol{\Sigma}^{-1}$ in equation 13. $\boldsymbol{R}$ is a diagonal
matrix of dimensions $T\times T$, wherein the main diagonal elements are the
squared ranges of the corresponding components in the data. These chosen
hyperparameter values yield satisfactory results in the simulation studies and
are consistently employed in all examples unless explicitly stated otherwise.
We compared the methods in the following aspects with different metrics:
estimation of the covariance function with $L_{2}$ norms and estimation of the
principal components with mean squared errors. Each simulation represents a
choice between $\text{Cov}_{1}(s,t)$ and $\text{Cov}_{2}(s,t)$ as the true
covariance function and the prior covariance function, which results in four
experimental settings. For each setting, a total of 10000 iterations is run
and we discard the first 5000 as the burn-in and use the remainder 5000 for
inference. Each setting is repeated 30 times with different seeds and the
results are the average of the estimated posterior means from 30 runs.

| Truth | Prior | RB-FPCA | BFPCA | FACE | PACE | pct(RB-FPCA is best)
---|---|---|---|---|---|---|---
|  | Cov1 | 11.857 | 13.981 |  |  |
| Cov1 | Cov2 | 12.864 | 15.239 | 13.609 | 17.288 | 0.67
|  | Cov1 | 14.163 | 17.135 |  |  |
Covariance Function | Cov2 | Cov2 | 13.288 | 16.435 | 14.545 | 17.209 | 0.53
Table 1: Simulation I: Estimations of the covariance function are evaluated by the $L_{2}$ norm.  Figure 2: Simulation I: Boxplots of $L_{2}$ norms between the estimates and the true covariance function. | Truth | Prior | RB-FPCA | BFPCA | FACE | PACE | pct(RB-FPCA is the best)
---|---|---|---|---|---|---|---
|  | Cov1 | 0.041 | 0.044 |  |  |
| Cov1 | Cov2 | 0.048 | 0.046 | 0.043 | 0.039 | 0.50
|  | Cov1 | 0.017 | 0.019 |  |  |
PC1 | Cov2 | Cov2 | 0.035 | 0.035 | 0.037 | 0.046 | 0.30
|  | Cov1 | 0.028 | 0.029 |  |  |
| Cov1 | Cov2 | 0.040 | 0.041 | 0.045 | 0.048 | 0.63
|  | Cov1 | 0.036 | 0.046 |  |  |
PC2 | Cov2 | Cov2 | 0.029 | 0.028 | 0.033 | 0.041 | 0.20
|  | Cov1 | 0.026 | 0.027 |  |  |
| Cov1 | Cov2 | 0.051 | 0.055 | 0.031 | 0.040 | 0.50
|  | Cov1 | 0.024 | 0.021 |  |  |
PC3 | Cov2 | Cov2 | 0.042 | 0.047 | 0.032 | 0.038 | 0.40
Table 2: Simulation I: Estimations of the first 3 principal components are
compared with Mean Squared Errors (MSEs).

From Table 1 and Figure 2, covariance functions are estimated well by the
proposed RB-FPCA method under the condition that the true covariance function
is $\text{Cov}_{1}(s,t)$. The proposed RB-FPCA method has a competitive
performance against the two frequentist methods, given that the true
covariance function is $\text{Cov}_{2}(s,t)$. Generally, a prior covariance
function constructed based on the true covariance function can enhance
estimation. Some certain degrees of subject matter expertise in the structure
of the covariance function can be incorporated into the analysis to improve
the estimation further. The estimations of the principal components measured
by MSEs in Table 2 also indicate that the proposed RB-FPCA method can provide
comparable performance to the BFPCA and some well-known frequentist methods.

###  3.2 Simulation II: Dense functional data with outliers

The second simulation considers the presence of outliers. The goal is to
investigate the finite-sample performance and robustness of the proposed
model. Different percentages of atypical observations are considered. A total
of 100 observations over 50 time points are generated in each case.

For clean observations, the data are generated from the following model:

$Y_{i}(t)=\mu(t)+\sum_{k=1}^{K}\sqrt{\lambda_{k}}Z_{i,k}\phi_{k}(t),\;i=1,\ldots,n,$
(15)

with $Z_{i,k}\stackrel{{\scriptstyle i.i.d}}{{\sim}}\mathcal{N}(0,1)$ and
eigenvalues $\lambda_{1}\geq\lambda_{2}\geq\cdots\geq\lambda_{K}>0$. Time
points are evenly spaced in the interval $[0,1]$. The mean function is set to
$\mu(t)=10\sin(2\pi t)\exp(-3t)$. We set $K=4$ and choose the eigenfunctions
$\phi_{k}(t)$ as the first $K$ eigenfunctions of the Matérn Covariance
function (Rasmussen, 2003):

$\displaystyle\text{Cov}_{\text{Matérn}}(s,t)=\sigma^{2}\frac{2^{1-\nu}}{\Gamma(\nu)}\left(\frac{\sqrt{2\nu}|s-t|}{\rho}\right)^{\nu}K_{\nu}\left(\frac{\sqrt{2\nu}|s-t|}{\rho}\right),$

where $\Gamma(\cdot)$ is the Gamma function and $K_{\nu}$ is a modified Bessel
function of the second kind. We set the parameters as $\sigma^{2}=1$,
$\rho=3$, and $\nu=1/2$. The eigenvalues $\lambda_{1}=0.83$,
$\lambda_{2}=0.08$, $\lambda_{3}=0.03$, and $\lambda_{4}=0.015$ were selected
so that they have similar ratios to those of the first four eigenvalues of the
Matérn covariance function. In this model, the first principal direction
summarizes major sources of variation among the curves. At the same time, the
third and fourth eigenfunctions will tend to add some complexity to the
covariance function.

To add outliers to the clean observations, we introduced atypical observations
to change the order of the principal directions to affect the estimated
eigenvalues and eigenfunctions. We added outliers with a Bernoulli random
variable $B_{i}\sim Bernoulli(1,p)$, where $p$ corresponds to the percentage
of outliers. We contaminated the clean observations as follows:

  1. i.

For curve $i$, generate the Bernoulli random variable $B_{i}$. If $B_{i}$ = 0,
then proceed to (a). If $B_{i}$ = 1, then proceed to (b).

    1. (a)

Generate the curve $Y_{i}(t)$ as described in equation 15

    2. (b)

Introduce outliers by altering the order of the principal directions. Sample
the scores for the second and third principal directions with

$\displaystyle\left(\begin{array}[]{c}Z_{i,2}\\\ Z_{i,3}\end{array}\right)\sim
MVN_{2}\left(\left(\begin{array}[]{c}20\\\
25\end{array}\right),\begin{bmatrix}1/16&0\\\ 0&1/16\end{bmatrix}\right).$

  2. ii.

Repeat step i for all curves, i = 1,…,n.

Figure 3 shows one example of the 100 observations generated over 50 time
points with an outlier percentage of $p=0.1$. This figure illustrates how
contamination alters the pattern of the clean data. Black lines represent the
clean samples, and blue lines correspond to 10 outlying samples.

We examined the behaviour of the estimated covariance function and the
eigenfunctions for both clean and contaminated samples with
$p=0.05,0.10,0.15$. Each value of $p$ is repeated 30 times with different
seeds, and the final results are the average of the estimated posterior means
from 30 runs. The prior covariance function used in this simulation is
$\text{Cov}(s,t)=\exp\\{-3(t-s)^{2}\\}$. The numbers of basis functions and
eigenfunctions are set to $P=15$ and $K=4$.

Figure 3: Simulation II: Clean samples with contaminated samples from the data generation process described in 3.2. The outlier percentage is 10%. | p | RB-FPCA | BFPCA | FACE | PACE | pct(RB-FPCA is best)
---|---|---|---|---|---|---
| 0.05 | 43.956 | 44.057 | 44.082 | 44.074 | 0.30
| 0.10 | 43.991 | 48.131 | 44.141 | 44.135 | 0.93
Covariance Function | 0.15 | 44.035 | 64.065 | 44.271 | 44.269 | 1.00
Table 3: Simulation II: Estimations of the covariance function are evaluated
by the $L_{2}$ norms.

(a)

(b)

Figure 4: Simulation II: (a) Boxplots of $L_{2}$ norms (in log scale) between the estimates and the true covariance function. Comparison is between the proposed RB-FPCA method and the other Bayesian methods. (b) Boxplots of $L_{2}$ norms between the estimates and the true covariance function. Comparison is between the proposed RB-FPCA method and the two frequentist methods.  | p | PCs | RB-FPCA | BFPCA | FACE | PACE
---|---|---|---|---|---|---
|  | PC1 | 1.539 | 1.613 | 1.601 | 1.541
|  | PC2 | 1.463 | 1.464 | 1.590 | 1.554
| 0.05 | PC3 | 2.257 | 2.041 | 2.495 | 1.382
|  | PC1 | 1.572 | 1.574 | 1.593 | 1.549
|  | PC2 | 1.294 | 1.503 | 1.568 | 1.552
| 0.10 | PC3 | 2.083 | 2.234 | 2.489 | 1.507
|  | PC1 | 1.556 | 1.572 | 1.584 | 1.558
|  | PC2 | 1.244 | 1.497 | 1.569 | 1.565
PCs | 0.15 | PC3 | 2.027 | 2.346 | 2.487 | 1.568
Table 4: Simulation II: Estimations of the first 3 principal components are
compared with angles (in radians) between the truth and estimates.

Table 3 compares the estimations of the covariance function in terms of the
$L_{2}$ norms with varying outlier percentages, averaged across 30
replications. It can be observed that RB-FPCA shows significant improvement
over BFPCA, providing a smaller error, and a moderate improvement over the two
modern frequentist methods. Figures 4a and 4b display the comparison of the
boxplots of the $L_{2}$ norms. These boxplots confirm that RB-FPCA outperforms
the other three methods in terms of median and dispersion. We also observe
that the better performance of RB-FPCA is more pronounced when the outlier
percentage is high. Table 4 summarizes the PC estimation results for various
outlier percentage values. We can see that RB-FPCA outperforms BFPCA and FACE
in estimating all the first 3 PCs and has a competitive performance compared
with PACE.

###  3.3 Simulation III: Sparse functional data with outliers

The third simulation examines the performance of the Bayesian FPCA on the
sparse longitudinal data. Each time we first generated a dense dataset with
100 noisy observations over 50 time points with the true mean function
$\mu(t)=\text{sin}(2\pi t)$. The true underlying covariance function is of
either $\text{Cov}_{1}(s,t)=\exp\\{-3(t-s)^{2}\\}$ or
$\text{Cov}_{2}(s,t)=\text{min}\\{s+1,t+1\\}$ depending on the experimental
settings. Outlier percentage is chosen from $5\%,10\%\text{ or }15\%$ with
equal probability. For outliers, noises of $N(0,3)$ are added to the sampling
process. For clean samples, a small amount of noise of $N(0,0.3)$ is added. To
get sparsely observed samples, a discrete uniform distribution $U(5,10)$ was
used to select a random number of times $n_{i}$ each curve was observed. The
observed times $t_{ij}$ satisfy $t_{ij}\sim U(-1,1),\text{ i.i.d. for
}i=1,\ldots,100,j=1,\ldots,n_{i}$. This data generation process produces
sparse longitudinal data with some noisy samples.

We tested the performance of our proposed RB-FPCA method with different prior
covariance functions. When the true covariance function is
$\text{Cov}_{truth}(s,t)=\exp\\{-3(t-s)^{2}\\}$, we set prior covariance
functions with the following forms:

  * •

$\text{Cov}_{prior1}(s,t)=\exp\\{-3(t-s)^{2}\\}$ (same as truth)

  * •

$\text{Cov}_{prior2}(s,t)=\min\left(s+1,t+1\right)$

  * •

$\text{Cov}_{prior3}(s,t)=\exp\\{-(t-s)^{2}\\}$

  * •

$\text{Cov}_{prior4}(s,t)=\text{covariance estimation from PACE method}$

When the true covariance function is
$\text{Cov}_{truth}(s,t)=\min\left(s+1,t+1\right)$, we set prior covariance
functions with the follwing forms:

  * •

$\text{Cov}_{prior1}(s,t)=\exp\\{-3(t-s)^{2}\\}$

  * •

$\text{Cov}_{prior2}(s,t)=\min\left(s+1,t+1\right)$ (same as truth)

  * •

$\text{Cov}_{prior3}(s,t)=\left(s+1\right)*\left(t+1\right)$

  * •

$\text{Cov}_{prior4}(s,t)=\text{covariance estimation from PACE method}$

Figure 5 shows the contour plots of covariance functions. For each setting, we
selected four prior covariance functions with different characteristics that
may have an effect on the performance of the RB-FPCA method. Specifically,
except for using the true covariance function as the prior covariance
function, we also included one prior covariance function, which has a similar
surface structure as the true covariance function and one prior covariance
function with a different shape. In addition, we also involved the smoothed
covariance surface estimated from the PACE method.

(a)  (b)

(c)  (d)

Figure 5: Simulation III: Contour plots of prior covariance functions.

The numbers of basis functions and eigenfunctions are preset as $P=5$ and
$K=5$. For each setting, a total of 4000 iterations are run, and we discard
the first 2000 as the burn-in and use the remainder 2000 for inference. Each
setting is repeated 50 times with different seeds and noise levels, and the
results are the average of the estimated means from 50 runs. We compared the
performance of the proposed RB-FPCA method with two frequentist methods, which
are the PACE method of Wang et al. (2016), implemented in fdapace package in
R, and a robust FPCA (sparseFPCA) method of Boente and Salibián-Barrera
(2021), implemented in sparseFPCA package in R. The performance of each model
is compared in terms of estimations of the correlation surface evaluated by
$L_{2}$ norms and estimations of the first three principal components assessed
by the angle (in radians) between the truth and estimates. Results are
presented in Figure 6.

(a)

(b)

(c)

(d)

Figure 6: Simulation III: Comparison of RB-FPCA method with PACE and
sparseFPCA methods. (a) and (c) shows boxplots of $L_{2}$ norms between the
estimates and the true correlation function. (b) and (d) shows boxplots of the
angle (in radians) of the estimations for the first 3 principal components.
The top row shows the results given the true covariance function is
$\text{Cov}_{truth}(s,t)=\exp\\{-3(t-s)^{2}\\}$. The bottom row shows the
results when the true covariance function is
$\text{Cov}_{truth}(s,t)=\min\left(s+1,t+1\right)$. The black dots in each box
correspond to the values of the mean.

Figures 6a and 6b demonstrate the advantage of using the RB-FPCA method over
the other two frequentist methods by comparing the estimations of the
correlation functions and the PCs. The proposed RB-FPCA method produces better
measures of the correlation functions with smaller $L_{2}$ norms and admits
small variations. Figures 6c and 6d show different methods’ performances under
the other true covariance function. The proposed RB-FPCA method outperforms
one of the frequentist methods, namely sparseFPCA, in terms of the estimations
of the correlation functions. PACE performs well for estimating the
correlation function, whereas the RB-FPCA method still has a potential
advantage if considering the variations of the estimations. Note that even
though the sparseFPCA method produces estimations for PCs with the smallest
variations, it does not outperform PACE and RB-FPCA methods in estimating the
correlation functions and PCs. The robustness of the sparseFPCA method may
explain the small variations of the estimations of PCs. The overall picture is
that the proposed RB-FPCA method has the ability to deal with sparse
longitudinal data, perform competitively with some well-known frequentist
methods, and provide a Bayesian framework to enable incorporating the domain
knowledge into the analysis via flexible prior settings.

##  4 Data Applications

###  4.1 Hawaii Ocean Oxygen Data

We first illustrate our method on the Hawaii ocean oxygen data collected from
the Hawaii Ocean Time-series Data Organization & Graphical System. Scientists
participating in the Hawaii Ocean Time-series (HOT) program have been
monitoring and making continuous measurements of the water column’s
hydrography, chemistry and biology at a station near Oahu, Hawaii, since
October 1988. The primary objective of this research is to monitor and
interpret the variability of physical and biogeochemical processes at deep-
water hydrostations and deliver a comprehensive overview of the ocean at a
representative site in the North Pacific subtropical region. The deep-water
station is visited approximately once a month by cruises to obtain water
samples from desired ocean depths. Observations of ocean data over long
periods are extremely valuable for climate studies. Researchers require
repeated measurements of oceanographic data to investigate slow or irregular
changes in natural processes or phenomena and some rapid event-driven
variations. Oceanic data can be downloaded from the Hawaii Ocean Time-series
Data Organization & Graphical System website
(https://hahana.soest.hawaii.edu/hot/hot-dogs/cextraction.html).

In this study, we analyzed the oxygen concentrations in units of $\mu$mol/kg
measuring at different depths below the sea surface from January 1, 2008, to
December 31, 2021. The oxygen concentration is measured every 2 meters at a
depth of 0 to 200 meters under the ocean’s surface. Such data with a shorter
time frame have been analyzed by Shi and Cao (2022). We obtained 133
trajectories to study the functional relationship between the oxygen
concentrations and depth below the sea surface. Each trajectory has 100 data
points with no missing values. All trajectories have measurements at the same
depths.

We applied the RB-FPCA method for dense data to the Hawaii ocean oxygen data
with a predetermined set of numbers of basis functions and eigenfunctions,
i.e., $P=15$ and $K=5$. The data were first detrended, and the prior
covariance function was chosen as $\exp\\{-3(t-s)^{2}\\}$. A total of 5000
Gibbs sampling iterations was run, and the first 2500 iterations were
discarded as the burn-in. This study aims to find patterns in data of high
dimensions, capture the primary mode of variation, and flag potential
outliers.

Figure 7 shows the first four FPCs from the RB-FPCA model with variance
explained by each FPC. The 95% credible intervals are included with the dashed
lines. The credible intervals generated within the Bayesian approach offer a
direct interpretation, reflecting the probability of the true parameter being
contained within a specific range, for instance, the 95% credible interval
indicating a 95% probability of encompassing the true parameter. The credible
intervals provide valuable information for interpreting the significance of
each FPCs in the analysis. In general, the credible interval for the first FPC
is narrow and relatively narrow for the remaining FPCs except at the
boundaries. The first FPC is positive and significant over the whole range at
a level of 0.05. It represents the weighted average of the oxygen level of
each visit to the deep-water station over a depth from 0 to 200 meters. Most
of the variation is due to the amount of horizontal shift from the mean
function. The second FPC explains about 18% of the total variability. The
second FPC can be interpreted as a change in the oxygen level over two depth
intervals. The second FPC is positive when the depth is between 0 to 150
meters and negative when the depth is above 150 meters, indicating the
relative difference in oxygen levels between 0 to 150 meters and 150 to 200
meters. The credible intervals indicate significant variations over the depth
intervals [50, 134] and [152, 200]. FPCs at a higher level represent more
complex phenomena. The top 3 FPCs capture over 90 % of the total variation in
the data.

Figure 7: Hawaii Ocean Oxygen Data: The top four estimated FPCs from the full
posterior mean of the covariance function with the proportion of variance
explained by each FPC. The 95% credible intervals are represented by the
dashed lines. The shaded areas denote credible intervals that do not contain
zeros. The numeric labels along the vertical lines correspond to the depths
associated with the boundaries of each respective shaded area.

We implemented the proposed approach of identifying the outliers using the FPC
scores. Different threshold values result in different trajectories being
flagged as abnormal. We investigated four threshold values, and the results
are presented in Figure 8. The result follows the expectation that fewer
outliers are detected as the threshold value increases. The detected atypical
trajectories lie near the top or bottom. Compared with other curves, abnormal
trajectories appear to have different curvature, or maintain high values over
depth levels, or decline rapidly.

Figure 8: Hawaii Ocean Oxygen Data: Detected outliers and the oxygen
trajectories. The dashed curve is the mean function calculated with all
trajectories. The values on the top left corner is the threshold value used in
the quantile of the $\chi^{2}_{2}$ distribution to identify possible outliers.
For example, the top left plot shows outliers with FPC scores have distances
larger than the 99% quantile of a $\chi^{2}_{2}$ distribution.

###  4.2 Annual sea surface temperature data

The second example considers the annual sea surface temperature data. The sea
surface temperature is the temperature of the ocean’s surface at its top
millimeter. As a fundamental measure of global climate change, sea surface
temperature provides a glimpse into the overall trend in the climate system.
Understanding the behaviour of the sea surface temperature changes and
discovering potential anomalies are essential in studying global climate
change. Sea surface temperature anomalies are characteristic of El Niño and La
Niña climate cycles, which can affect weather patterns worldwide. Strong and
localized anomalies may identify ocean currents in sea surface temperature.
Anomalies in sea surface temperature over many years can be signs of regional
or global climate change, for example, global warming. In addition to their
scientific value, sea surface temperature anomalies are also useful for
practical purposes. For instance, an anomalous temperature (warm or cool) in
coastal areas can favor certain organisms in an ecosystem over others,
resulting in a thriving or declining population of bacteria, algae, or fish.
Such data have been analyzed by Hyndman and Shang (2010), Sun and Genton
(2011), Xie et al. (2017), and Dai et al. (2020). Annual sea surface
temperature data can be downloaded from the Climate Prediction Center website
(https://www.cpc.ncep.noaa.gov/data/indices/).

The dataset consists of monthly observations from Niño 1+2 region from January
1950 to December 2021. We obtained 72 functional observations observed on a
common grid of 12 time points. They correspond to the monthly sea surface
temperature over 72 years. The RB-FPCA method for dense data was applied to
find projections of maximum variance. The numbers of basis functions and
eigenfunctions are fixed with $P=10$ and $K=5$. The prior covariance function
that is estimated is $\exp\\{-3(t-s)^{2}\\}$. A total of 5000 MCMC iterations
were used for estimation after 2500 burn-in iterations. The estimations are
only calculated for the time points that have been sampled.

(a)

(b)

Figure 9: Annual sea surface temperature data: (a) FPC scores for the first
two FPCs. Possible outliers and normal observations are differentiated by
colors. (b) Plot of the robust Mahalanobis distances for observations.
Potential outliers are highlighted. These four suspected outliers correspond
to the same observations identified in both (a) and (b).

The first two FPCs explain 94% of the variance. The first FPC is almost
constant below the zero axis, representing that the major variation is the
degree of horizontal shift below the mean function. The second FPC crosses the
zero axis once near July, which can be interpreted as the relative change of
sea surface temperature between spring and fall months. To investigate the
existence of potential outliers, the plot of the FPC scores for the first two
FPCs is shown in Figure 9a. We have observed that four observations deviate
significantly from the majority of the observations, and have been identified
as outliers. The graphical representation of robust Mahalanobis distances in
Figure 9b provides further evidence of the outlying behaviors exhibited by
these four observations.

(a)

(b)

Figure 10: Annual sea surface temperature data: (a) Detected outliers and
annual temperature trajectories. The four outliers (blue) are annotated with
years, and the colors of the annotations are in line with the legends in panel
(b). The detected El Niño events are differentiated by square and triangle.
(b) Plot of the whole history of monthly sea surface temperature from January
1950 to December 2021. The locations of the four detected outliers are
colored.

From a practical aspect, detecting sea surface temperature anomalies is
essential to identify El Niño and La Niña events. El Niño events are widely
used to describe the warming of sea surface temperature that occurs every few
years. In contrast, La Niña events correspond to years with abnormally low sea
surface temperatures. We can determine the abnormal annual sea surface
temperature curves by FPC scores with robust Mahalanobis distances larger than
99.5% quantile of a $\chi^{2}_{2}$ distribution. Four trajectories (1982,
1983, 1997, and 1998) are identified as anomalies and shown in Figure 10. This
finding confirms the potential outliers found from the preliminary
investigation of the FPC scores, as depicted in Figure 9. The winters from
1982 to 1983 and 1997 to 1998 have the highest temperatures among all years,
and the temperatures remain at high values over the summers of 1983 and 1998.
According to National Climatic Data Center reports, these periods correspond
to two of the strongest El Niño events in history and are illustrated by
squares and triangles in 10a. Although our method cannot flag all significant
El Niño or La Niña years, we can detect some of the strongest activities. By
decreasing the threshold value in our outlier detection procedure, more
outlying curves are expected to appear.

###  4.3 Sparse CD4 Data

The last example considers the CD4 data, a sparse longitudinal dataset
collected within the Multicenter AIDS Cohort Study (MACS). Study participants
include 1809 HIV-infected men at the start of the study and 371 men who were
seronegative at entry and seroconverted during the study period (see Zeger and
Diggle (1994)). In our study, we used the data from the catdata package in R,
which includes 2376 measurements of the number of CD4 cells taken over time on
369 seroconverters. The variable of interest is the CD4 cell counts measured
in years since seroconversion, which can be used to assess disease
progression. Specifically, we are interested in the typical decay of CD4 cell
counts over time and the variability across subjects. The CD4 data are sparse,
owing to unequal numbers of repeated measurements for each subject and the
different timing of measurements for each subject. We obtained 241
trajectories after removing observations with less than three measurements,
and the number of measurements per subject ranges from 3 to 11, with a median
of 5. Figure 11a shows the data together with a smooth estimate of the mean
function from the PACE method. The overall trend in the CD4 cell counts is
decreasing.

We applied the model described in 2.5 to fit the CD4 data and compare the
estimates of the covariance function with the PACE method. We used a fixed
number of basis functions and eigenfunctions, i.e., $P=3$ and $K=3$, and ran a
total of 5000 Gibbs sampling iterations with 2500 burn-in iterations. We
utilized the PACE method’s estimates of the covariance function as the prior
covariance function because such an informative prior expresses one’s beliefs
about this quantity. The tuning parameters for the PACE method were set to
default values as described in Wang et al. (2016).

We identified some potential outliers using the estimated FPC scores from
conditional expectations. The top five most outlying trajectories are
identified and highlighted in Figure 11b. These outlying trajectories appear
to have either some rapid decreases or increases compared to the rest. They
also tend to have overall CD4 cell counts very high compared to the others
during the whole period. We estimated the covariance functions using all the
trajectories and only normal trajectories. The comparisons of the covariance
function estimates under RB-FPCA and PACE are displayed in Figure 12. The RB-
FPCA and PACE estimators have similar overall shapes, whereas RB-FPCA produces
a more smooth surface. When outliers are removed, RB-FPCA induces a similar
covariance surface compared with the estimation from complete data except near
the boundary, demonstrating the robustness of the RB-FPCA model.

(a)

(b)

Figure 11: Sparse CD4 Data: (a) Data contain 241 observed trajectories for CD4
cell counts. The dashed line is the smooth estimate of the mean function. (b)
Detected five most outlying trajectories in CD4 data.

(a) RB-FPCA (complete data)

(b) PACE (complete data)

(c) RB-PFCA (outliers removed)

(d) PACE (outliers removed)

Figure 12: Sparse CD4 Data: Comparison of smooth estimates of the covariance
functions under RB-FPCA and PACE. The covariance surfaces are estimated with
all trajectories in (a) and (b) and with the five most outlying trajectories
removed in (c) and (d).

##  5 Discussion

In this work, we developed a robust Bayesian functional principal components
analysis method that utilizes the class of skew elliptical distributions in
the modelling. The general and flexible class of skew elliptical distributions
can serve as an alternative to the symmetric distribution commonly assumed in
the previous works of Bayesian FPCA. The proposed method is able to handle the
sparse longitudinal data in which only a few observations per trajectory
(possibly sampled at irregular intervals) are available. Such data are
commonly seen and attracted interest in the area of functional data analysis.
We have shown in simulation studies and data applications that the proposed
method can effectively capture the major variation and provide useful
information in the presence of noises and outliers in the data.

Multiple future directions of research are presented as follows. The class of
skew elliptical distribution contains many standard families of distributions
such as the multivariate skew normal distribution, Student-$t$ and Pearson
type II distributions. We have only presented the case of multivariate skew
normal distribution in this work, and multivariate skew-t distribution is
another particular example which has the ability to deal with heavier-tailed
data. Extending the current framework to fit the multivariate skew-t
distribution is straightforward by modifying some equations presented in 2.2.
Specifically, we introduce $n$ i.i.d random variables $w_{i}$, and the
hierarchical representation of the model with multivariate skew-t distribution
for $i=1,\ldots,n$ is

$\displaystyle\boldsymbol{Y}_{i}\>|\>\cdot$ $\displaystyle\sim
MVN_{T}(\boldsymbol{H}_{P}\boldsymbol{U}_{K}\boldsymbol{\beta}_{i,K}+\boldsymbol{D}\boldsymbol{z}_{i},\>\boldsymbol{\Sigma}/w_{i}),$
$\displaystyle w_{i}$ $\displaystyle\sim Gamma(\nu_{w_{i}}/2,\nu_{w_{i}}/2),$
$\displaystyle\nu_{w_{i}}$ $\displaystyle\sim Gamma(1,0.1)I(\nu_{w_{i}}>2).$

The other distributional specifications remain the same. For the
specifications of the full conditional distributions, the only part that
depends on $w_{i}$ is for $\boldsymbol{z}_{i}$ in equation 11 where
$\boldsymbol{A}_{i}=\boldsymbol{I}+w_{i}\boldsymbol{D}\boldsymbol{\Sigma}^{-1}\boldsymbol{D}$
and
$\boldsymbol{a}_{i}=w_{i}\boldsymbol{D}\boldsymbol{\Sigma}^{-1}(\boldsymbol{Y}_{i}-\boldsymbol{H}_{P}\boldsymbol{U}_{K}\boldsymbol{\beta}_{i,K})$.
Implementing the multivariate skew-t distribution could improve the model
performance, which requires further exploration.

Furthermore, the standard FPCA operates in a static way, therefore, does not
provide an adequate dimension reduction when considering the functional time
series data. Hörmann et al. (2015) proposed a dynamic FPCA which takes into
account the serial dependence in the functional time series. A Bayesian
version of the robust FPCA for functional time series is worth investigating.
Regarding posterior inference, there exist popular alternative algorithms,
such as reversible jump Markov Chain Monte Carlo or sequential Monte Carlo,
which could provide some advantages in computation and model selection.
Further research into the computational aspects of the robust Bayesian FPCA
methods is necessary.

SUPPLEMENTARY MATERIAL

R code for RB-FPCA

R code implements the proposed model described in the article. All datasets
used as examples in the article are also included.

## References

  * Azzalini and Valle (1996) Azzalini, A. and A. D. Valle (1996).  The multivariate skew-normal distribution.  Biometrika 83(4), 715–726.
  * Behseta et al. (2005) Behseta, S., R. E. Kass, and G. L. Wallstrom (2005).  Hierarchical models for assessing variability among functions.  Biometrika 92(2), 419–434.
  * Boente and Salibián-Barrera (2021) Boente, G. and M. Salibián-Barrera (2021).  Robust functional principal components for sparse longitudinal data.  METRON 79(2), 159–188.
  * Branco and Dey (2001) Branco, M. D. and D. K. Dey (2001).  A general class of multivariate skew-elliptical distributions.  Journal of Multivariate Analysis 79(1), 99–113.
  * Cardot (2000) Cardot, H. (2000).  Nonparametric estimation of smoothed principal components analysis of sampled noisy functions.  Journal of Nonparametric Statistics 12(4), 503–538.
  * Chib (1995) Chib, S. (1995).  Marginal likelihood from the Gibbs output.  Journal of the American Statistical Association 90(432), 1313–1321.
  * Crainiceanu and Goldsmith (2010) Crainiceanu, C. M. and A. J. Goldsmith (2010).  Bayesian functional data analysis using WinBUGS.  Journal of Statistical Software 32, 1–33.
  * Dai et al. (2020) Dai, W., T. Mrkvička, Y. Sun, and M. G. Genton (2020).  Functional outlier detection and taxonomy by sequential transformations.  Computational Statistics & Data Analysis 149, 106960.
  * Dauxois et al. (1982) Dauxois, J., A. Pousse, and Y. Romain (1982).  Asymptotic theory for the principal component analysis of a vector random function: some applications to statistical inference.  Journal of Multivariate Analysis 12(1), 136–154.
  * Ferraty and Vieu (2006) Ferraty, F. and P. Vieu (2006).  Nonparametric Functional Data Analysis.  New York: Springer.
  * Fukunaga and Koontz (1970) Fukunaga, K. and W. L. Koontz (1970).  Representation of random processes using the finite Karhunen-Loeve expansion.  Information and Control 16(1), 85–101.
  * Hall and Hosseini-Nasab (2006) Hall, P. and M. Hosseini-Nasab (2006).  On properties of functional principal components analysis.  Journal of the Royal Statistical Society: Series B (Statistical Methodology) 68(1), 109–126.
  * Han and Carlin (2001) Han, C. and B. P. Carlin (2001).  Markov chain Monte Carlo methods for computing Bayes factors: A comparative review.  Journal of the American Statistical Association 96(455), 1122–1132.
  * Hörmann et al. (2015) Hörmann, S., Ł. Kidziński, and M. Hallin (2015).  Dynamic functional principal components.  Journal of the Royal Statistical Society: Series B (Statistical Methodology) 77(2), 319–348.
  * Horváth and Kokoszka (2012) Horváth, L. and P. Kokoszka (2012).  Inference for Functional Data with Applications.  New York: Springer.
  * Hsing and Eubank (2015) Hsing, T. and R. Eubank (2015).  Theoretical Foundations of Functional Data Analysis, with an Introduction to Linear Operators.  Chichester: John Wiley & Sons.
  * Hyndman and Shang (2010) Hyndman, R. J. and H. L. Shang (2010).  Rainbow plots, bagplots, and boxplots for functional data.  Journal of Computational and Graphical Statistics 19(1), 29–45.
  * James et al. (2000) James, G. M., T. J. Hastie, and C. A. Sugar (2000).  Principal component models for sparse functional data.  Biometrika 87(3), 587–602.
  * Jeffreys (1935) Jeffreys, H. (1935).  Some tests of significance, treated by the theory of probability.  In Mathematical Proceedings of the Cambridge Philosophical Society, Volume 31, pp. 203–222. Cambridge University Press.
  * Jiang et al. (2020) Jiang, L., Y. Zhong, C. Elrod, L. Natarajan, R. Knight, and W. K. Thompson (2020).  BayesTime: Bayesian functional principal components for sparse longitudinal data.  ArXiv Preprint ArXiv:2012.00579.
  * Jolliffe (2002) Jolliffe, I. T. (2002).  Principal component analysis for special types of data.  Springer.
  * Margaritella et al. (2021) Margaritella, N., V. Inácio, and R. King (2021).  Parameter clustering in Bayesian functional principal component analysis of neuroscientific data.  Statistics in Medicine 40(1), 167–184.
  * Paul and Peng (2009) Paul, D. and J. Peng (2009).  Consistency of restricted maximum likelihood estimators of principal components.  The Annals of Statistics 37(3), 1229–1271.
  * Ramsay and Silverman (2005a) Ramsay, J. and B. Silverman (2005a).  Principal components analysis for functional data.  pp. 147–172.
  * Ramsay and Silverman (2005b) Ramsay, J. O. and B. W. Silverman (2005b).  Functional Data Analysis (Second ed.).  New York: Springer.
  * Rasmussen (2003) Rasmussen, C. E. (2003).  Gaussian processes in machine learning.  In Summer School on Machine Learning, pp. 63–71. Springer.
  * Reiss et al. (2017) Reiss, P. T., J. Goldsmith, H. L. Shang, and R. T. Ogden (2017).  Methods for scalar-on-function regression.  International Statistical Review 85(2), 228–249.
  * Rice and Silverman (1991) Rice, J. A. and B. W. Silverman (1991).  Estimating the mean and covariance structure nonparametrically when the data are curves.  Journal of the Royal Statistical Society: Series B (Methodological) 53(1), 233–243.
  * Sahu et al. (2003) Sahu, S. K., D. K. Dey, and M. D. Branco (2003).  A new class of multivariate skew distributions with applications to Bayesian regression models.  Canadian Journal of Statistics 31(2), 129–150.
  * Shi and Cao (2022) Shi, H. and J. Cao (2022).  Robust functional principal component analysis based on a new regression framework.  Journal of Agricultural, Biological and Environmental Statistics, 1–21.
  * Suarez and Ghosal (2017) Suarez, A. J. and S. Ghosal (2017).  Bayesian estimation of principal components for functional data.  Bayesian Analysis 12(2), 311–333.
  * Sun and Genton (2011) Sun, Y. and M. G. Genton (2011).  Functional boxplots.  Journal of Computational and Graphical Statistics 20(2), 316–334.
  * Uhlig (1994) Uhlig, H. (1994).  On singular Wishart and singular multivariate Beta distributions.  The Annals of Statistics, 395–405.
  * Van Der Linde (2008) Van Der Linde, A. (2008).  Variational Bayesian functional PCA.  Computational Statistics & Data Analysis 53(2), 517–533.
  * Wang et al. (2016) Wang, J.-L., J.-M. Chiou, and H.-G. Müller (2016).  Functional data analysis.  Annual Review of Statistics and Its Application 3, 257–295.
  * Xiao et al. (2016) Xiao, L., V. Zipunnikov, D. Ruppert, and C. Crainiceanu (2016).  Fast covariance estimation for high-dimensional functional data.  Statistics and Computing 26(1), 409–421.
  * Xie et al. (2017) Xie, W., S. Kurtek, K. Bharath, and Y. Sun (2017).  A geometric approach to visualization of variability in functional data.  Journal of the American Statistical Association 112(519), 979–993.
  * Yao et al. (2005) Yao, F., H.-G. Müller, and J.-L. Wang (2005).  Functional data analysis for sparse longitudinal data.  Journal of the American Statistical Association 100(470), 577–590.
  * Yohai (1987) Yohai, V. J. (1987).  High breakdown-point and high efficiency robust estimates for regression.  The Annals of statistics, 642–656.
  * Zeger and Diggle (1994) Zeger, S. L. and P. J. Diggle (1994).  Semiparametric models for longitudinal data with application to CD4 cell numbers in HIV seroconverters.  Biometrics, 689–699.