# Registration for Incomplete Non-Gaussian Functional Data

Alexander Bauer
Department of Statistics, Ludwig-Maximilians-Universität, Munich, Germany
and
Fabian Scheipl
Department of Statistics, Ludwig-Maximilians-Universität, Munich, Germany
and
Helmut Küchenhoff
Department of Statistics, Ludwig-Maximilians-Universität, Munich, Germany
and
Alice-Agnes Gabriel
Department of Earth and Environmental Sciences,
Ludwig-Maximilians-Universität, Munich, Germany  This work was supported by
the German Research Foundation (DFG) under Grant KU 1359/4-1; and the German
Federal Ministry of Education and Research (BMBF) under Grant No. 01IS18036A.

###### Abstract

Accounting for phase variability is a critical challenge in functional data
analysis. To separate it from amplitude variation, functional data are
registered, i.e., their observed domains are deformed elastically so that the
resulting functions are aligned with template functions. At present, most
available registration approaches are limited to datasets of complete and
densely measured curves with Gaussian noise. However, many real-world
functional data sets are not Gaussian and contain incomplete curves, in which
the underlying process is not recorded over its entire domain. In this work,
we extend and refine a framework for joint likelihood-based registration and
latent Gaussian process-based generalized functional principal component
analysis that is able to handle incomplete curves. Our approach is accompanied
by sophisticated open-source software, allowing for its application in diverse
non-Gaussian data settings and a public code repository to reproduce all
results. We register data from a seismological application comprising
spatially indexed, incomplete ground velocity time series with a highly
volatile Gamma structure. We describe, implement and evaluate the approach for
such incomplete non-Gaussian functional data and compare it to existing
routines.

Keywords: functional data analysis; phase variability; amplitude variability;
curve alignment; partially observed curves; seismology.

##  1 Introduction

Dealing with phase variability is crucial in functional data analysis. Many
different approaches exist for registering curves (see e.g. Marron et al.,
2015), but their application to diverse real world data settings remains
challenging. Most existing approaches target small to intermediate datasets of
completely observed curves with small amounts of Gaussian noise, evaluated on
dense, regular grids. Especially the registration of _incomplete_ curves,
i.e., curves whose underlying process is not observed from its natural
starting point all the way to its natural endpoint, has received only limited
attention so far (see e.g. Matuk et al., 2019; Bryner and Srivastava, 2021).
However, such data arise in many fields. Missing information about the initial
development of some processes, i.e. leading incompleteness, can be caused by
different starting conditions of subjects at the beginning of a study.
Trailing incompleteness towards the end of the underlying processes is present
in experiments and studies with a fixed endpoint that causes right-censoring,
or in panel studies with relevant dropout rates. If both types of
incompleteness are present, we use the term full incompleteness.
— Notation
In this functional data setting, we observe discretized incomplete curves
$Y_{i}(t_{i}^{*})$, $i=1,\ldots,N$ over individual chronological time domains
$T_{i}^{*}$, where the observed grids
$\boldsymbol{t}_{i}^{*}=[t^{*}_{ij}]_{j=1,\dots,n_{i}}$ may be irregular or
sparse, and chronological domains $T_{i}^{*}=[t_{\min,i}^{*},t_{\max,i}^{*}]$
are defined by the individual observation periods. W.l.o.g. we assume that all
observed curves are realizations of stochastic processes over a common
underlying internal time domain $T=[0,1]$ and $T_{i}^{*}\subseteq
T\;\forall\,i$. Registering the curves requires estimating inverse warping
functions $h^{-1}_{i}:T_{i}^{*}\mapsto T$ that account for the data’s phase
variation and map individual chronological times to the internal time of the
underlying process. The resulting registered curves
$Y_{i}(t)=Y_{i}(h^{-1}_{i}(t_{i}^{*}))$ only contain amplitude variation.

Figure 1: Simulated observed curves with trailing incompleteness (left
column), registered curves only comprising amplitude variation (center), and
inverse warping functions visualizing phase variation (right). Registration
was performed with the complete curve SRVF approach of Srivastava et al.
(2011) with function time_warping() from the R package fdasrvf (top row;
Tucker, 2020) and our incomplete curve approach (bottom row). Note the extreme
time dilation around the blue curve’s peak in the top row, which yields a
highly implausible registered curve. Figure 2: Typical seismic observations
recorded at different hypocentral distances. As a preprocessing step, the
leading zero measurements of the raw curves (left pane) are cut off and
absolute ground velocities are analyzed on the time since the first relevant
absolute ground velocity measurement $t^{*}_{0}$ (right).

Applying conventional “complete curve” registration methods to incomplete
curves often leads to nonsensical results (see Figure 1). This is caused by
the implicit, unwarranted assumption that the endpoints of the individual
observed chronological domains $T_{i}^{*}$ are identical to those of the
global internal time domain $T$.
— Data setting
Our approach is motivated by synthetic seismic data originating from large-
scale numerical _in silico_ experiments based on the 1994 magnitude $6.7$
earthquake in Northridge (California). The experiments were performed by the
Department of Earth and Environmental Sciences (LMU Munich, Germany) using the
software SeisSol (Pelties et al., 2014; Uphoff et al., 2017,
github.com/SeisSol/SeisSol) and are used to assess critical geophysical
parameters associated with high ground motion in the event of an earthquake.
The simulated ground motion curves contain highly relevant phase variation due
to different propagation velocities of the seismic waves and their varying
distance to the hypocenter of the earthquake. The complete data comprise 800
000 curves, each recorded over 30 seconds with a frequency of 2Hz. Seismic
activity has not subsided after 30 seconds in many cases, so the curves are
mostly incomplete. For more information on the data setting see Bauer et al.
(2017), Bauer et al. (2018) and Happ et al. (2019). In order to investigate
the structure of phase and amplitude variability of these highly variable and
spatially indexed data, we register them to spatially varying template
functions learnt from the data and represent the registered curves in a lower-
dimensional space. Note that these data do not have a simple structure with
additive Gaussian noise since absolute ground motion velocities are
nonnegative and higher values entail higher variability (see Figure 2).

We also apply our method to a version of the well-known Berkeley child growth
study data (Ramsay and Silverman, 2005). The data contain annual measurements
of the body heights of 39 boys and 54 girls from ages 1 to 18. The focus lies
on the first derivatives of the data, i.e., the speed of growth in different
stages of childhood and adolescence. We simulate full incompleteness in these
data by drawing an artificial initial age for every child in the first quarter
of the time domain as well as an individual drop-out year in the second half
of the domain. The observed curves with simulated incompleteness are
visualized in Figure 6.
— Study aim and contributions
Motivated by the challenges inherent in the seismic data, we require a curve
registration method that (i) is able to handle incomplete curves, (ii) is
applicable to non-Gaussian data and (iii) includes a lower-dimensional
representation of the registered curves. The latter is especially of interest
for further analyses of the estimated phase and amplitude variation structure.
To achieve these goals, we (i) derive and implement a novel penalized approach
for incomplete data registration and (ii, iii) derive and implement extensions
of the methods introduced by Wrobel et al. (2019) for exponential family
distributions beyond Binomial and Gaussian data. Furthermore, we implement
multiple computational improvements in the underlying software stack to
accelerate and stabilize the algorithm. Our implementation is available in the
registr package (Wrobel and Bauer, 2021) for the open-source software R (R
Core Team, 2020). All analyses in this paper can be reproduced based on the
code and data in our public GitHub repository (Bauer, 2021).

Before introducing our method in Section 3 we give an overview of prior work
in Section 2. Sections 4 and 5 comprise an extensive simulation study and the
applications. We end with a discussion in Section 6.

##  2 Related Work

###  2.1 Registration

Accounting for phase variation is often critical when analyzing functional
data. We refer to Marron et al. (2015) and the references therein for an
introduction to the general issue and a detailed overview of available
registration approaches. Many early approaches focused on the alignment of
curves towards given (salient) structures, i.e., landmark registration (e.g.
Kneip and Gasser, 1992). More recent proposals mostly perform alignment
towards template functions that are either based on domain knowledge or
estimated based on some measure of centrality, with such estimates often
iteratively refined over the course of the registration procedure. Warping
functions are commonly estimated purely nonparametrically (e.g. Chakraborty
and Panaretos, 2017; Tucker, 2020), as (piecewise) linear functions (e.g.
Sangalli et al., 2010; Vitelli, 2019; McDonnell et al., 2021) or in a basis
expansion. Common examples for the latter are the use of (penalized) B-spline
(Telesca and Inoue, 2008; Wrobel et al., 2019) or Fourier (Mattar et al.,
2009) bases or of warplets (Claeskens et al., 2010).
— SRVF-based registration approaches
One of the most popular approaches is the square-root velocity function (SRVF)
framework introduced by Srivastava et al. (2011), who showed that the warping-
invariant Fisher-Rao metric quantifies pure amplitude distances and is
equivalent to the simple $L_{2}$-metric in the SRVF space. Recently, Guo et
al. (2020) have extended this framework to jointly analyze amplitude, phase
and spatial variation.

Cheng et al. (2016), Kurtek (2017) and Lu et al. (2017) adapted the SRVF
approach to perform registration in a Bayesian setting. Bayesian approaches
were also introduced for data settings with stronger noise under informative
priors (Matuk et al., 2019; Tucker et al., 2021). These Bayesian approaches
can provide a full representation of the joint phase and amplitude
uncertainty, but are computationally very demanding. The method of Matuk et
al. (2019) handles sparse and fragmented Gaussian functional data where
measurements are only available over some parts of the observation domain, but
no software implementation was publicly available at the time of writing.

More recently, Nunez et al. (2021) and Chen and Srivastava (2021) introduced
the neural network-based frameworks SrvfNet and SrvfRegNet, respectively. Both
approaches build on the SRVF framework and enable a highly efficient
estimation and prediction of warping functions in large-scale data settings
for registering curves to their Karcher mean or fixed and pre-specified
template functions.

— Other registration approaches
All the above approaches are limited to continuous functional data, mostly
under the assumption of Gaussian errors. Some extensions to binary functional
data such as Wu and Srivastava (2014) and Panaretos and Zemel (2016) rely on a
pre-smoothing step to obtain a continuous representation of the curves. The
_congealing_ approach of Learned-Miller (2005) (adapted to functional data by
Mattar et al., 2009), which iteratively optimizes measures of alignment like
the integrated point-wise differential entropy via gradient descent, is
applicable in diverse data situations and computationally efficient. Wrobel et
al. (2019) utilizes a likelihood-based optimization approach for exponential
family data which is able to practically handle moderate-to-large scale
datasets (Wrobel et al., 2021).
— Joint registration and low-rank representations
For analyzing the main modes of amplitude and phase variation, it is common to
estimate a low-rank representation of the registered curves or, at least, the
template functions that serve as registration targets, and potentially also of
the estimated warping functions. Tucker et al. (2013), Hadjipantelis et al.
(2015), Lee and Jung (2016) and Happ et al. (2019) use (joint) functional
principal component analysis (FPCA) for finding compact representations of
both phase and amplitude modes of variation.

Tucker (2014) (utilizing the SRVF approach of Srivastava et al., 2011, for
registration), Wagner and Kneip (2019) and Kneip and Ramsay (2008) (optimizing
a least squares criterion) and Wrobel et al. (2019) (optimizing an exponential
family likelihood) all utilize iterative algorithms to successively refine
warping functions that lead to registered curves whose amplitude variation can
be represented in terms of a low-rank FPC basis.
— Incomplete curve registration
Comparatively few approaches have been developed so far for the registration
of incomplete curves. Some heuristic approaches are available in the field of
dynamic time warping (DTW) for time series analysis. Subsequence DTW offers a
framework to find a subsequence of a (fully observed) template curve to which
a partially observed curve can be matched (see Müller, 2015, 7.2). Tormene et
al. (2009) introduce an algorithm for “open-begin and open-end” DTW (OBE-DTW)
that is also able to handle full incompleteness.

More sophisticated registration approaches were introduced only recently.
Sangalli et al. (2010) and Vitelli (2019) make use of linear warping functions
with free starting points and endpoints. The observed curve domains are
constrained so that they dilate or extend the domain by a factor $0.9-1.1$ to
ensure reasonable warpings. As noted above, Matuk et al. (2019) allows to
analyze fragmented functional data, but their approach is feasible only for
small to intermediate sets of Gaussian data.

Bryner and Srivastava (2021) introduced an approach for elastic partial
matching to tackle trailing incompleteness. Before registering each curve to
its template using the complete curve SRVF approach of Srivastava et al.
(2011) they estimate the time scaling necessary to (partially) match the
observed domain of a specific curve to the domain of the template curve and
perform the registration only on the intersection of the curves’ domains. Both
steps are combined in a joint, gradient-based algorithm. At the time of
writing, no implementation of their method was available on request.
— Software implementations
Multiple packages for the statistical open-source software R (R Core Team,
2020) exist that implement registration approaches. Basic approaches outlined
in Ramsay and Silverman (2005) are implemented in package fda (Ramsay et al.,
2020). The OBE-DTW approach of Tormene et al. (2009) is available in dtw
(Giorgino, 2009). Package fdasrvf (Tucker, 2020) implements the SRVF
registration of Srivastava et al. (2011) and the iterative procedure of Tucker
(2014) for finding overall similar registered curves with a low-rank FPCA
representation. Code for Wagner and Kneip (2019) is available on GitHub
(Wagner, 2020). R package registr (Wrobel and Bauer, 2021) implements the
likelihood-based approach of Wrobel et al. (2019) and the incomplete curve
extensions presented in this work for various exponential family
distributions.

###  2.2 Generalized Functional Principal Component Analysis

Functional principal component analysis (FPCA) is a technique to analyze and
represent functional data in terms of their main modes of variation (Ramsay
and Silverman, 2005). To purely represent amplitude variation, FPCA is most
commonly applied to curves without phase variation, potentially after an
initial registration step. As noted above, the concept of FPCA was also
extended to separately (Tucker et al., 2013) or simultaneously (Happ et al.,
2019) analyze amplitude and phase variation. Multiple approaches exist for
FPCA on sparse or partially observed functional data, but mostly assume
Gaussianity (c.f. Stefanucci et al., 2018). Adaptations to the non-Gaussian
case for performing generalized FPCA (GFPCA) do exist, but have to be assessed
with care since marginal estimation of the overall mean can introduce bias
(Gertheiss et al., 2017).
— Probabilistic GFPCA
A popular method for multivariate non-Gaussian data is probabilistic FPCA
(Tipping and Bishop, 1999), based on likelihood optimization. For the case of
Gaussian functional data, James et al. (2000) and Rice and Wu (2001) use a
similar approach based on mixed-effects regression. Zhou et al. (2008) adapted
these methods in a paired-curve setting. Huang et al. (2014) extended the
ideas of James et al. (2000) and Zhou et al. (2008) to non-Gaussian functional
data and combined them with a clustering approach. Recently, Wrobel et al.
(2019) further adapted the mixed model-based approach by introducing a link
function and estimating the GFPCA based on computationally efficient
variational approximations. To the best of our knowledge, efficient
approximations are available only for Gaussian and binary data and are not
directly adaptable to further exponential family distributions.

Bayesian adaptations of probabilistic FPCA to non-Gaussian settings were
introduced by van der Linde (2009) for binary and count data, and Goldsmith et
al. (2015) with an extension to multilevel data. While these Bayesian
approaches can provide a full representation of the underlying uncertainty and
show good performance in sparse data situations (c.f. Gertheiss et al., 2017),
they are computationally demanding.
— Two-step GFPCA
A nonparametric approach to Gaussian FPCA was introduced by Yao et al. (2005)
and adapted by Hall et al. (2008) for the non-Gaussian case. Serban et al.
(2013) extended this method to further handle multilevel binary data with
potentially rare events. Li and Guan (2014) used a similar approach to model
point processes with a spatio-temporal correlation structure. Gertheiss et al.
(2017) showed that the marginal mean estimates proposed by Hall et al. (2008)
can introduce bias in non-Gaussian data settings and tackled this issue by
plugging the eigenfunction estimates into a generalized additive mixed model
to achieve an estimation of the mean structure conditional on FPC scores
represented as random effects.
— Some notes on consistency
Consistent estimators for the covariance operator are crucial for FPCA. The
consistency and convergence rates of estimators for covariance operators were
thoroughly studied for differently dense data settings (c.f. Wang et al.,
2016; Cao et al., 2016), their properties in the presence of stronger,
potentially non-Gaussian noise, however, remain an area of active research.
Standard techniques for covariance estimation quickly become computationally
infeasible in high-dimensional (Li et al., 2020) or irregular (Cederbaum et
al., 2018) data settings and algorithmic innovations are required. Recently,
Sarkar and Panaretos (2021) introduced promising neural network architectures
for the efficient, nonparametric approximation of (multidimensional)
covariance operators and their eigen-decomposition.

Several studies evaluated the consistency of covariance and FPCA estimators
for incomplete curve settings. When incompleteness originates from a missing
completely at random (MCAR) process and measurements are dense, established
estimators for the mean, the covariance and for eigenfunctions and eigenvalues
are consistent (Kraus, 2015). For the subject-specific functional principal
component (FPC) scores, Kraus (2015) introduced a consistent estimator. His
comparison to the PACE approach of Yao et al. (2005) indicates that the bias
of comparable conditional methods for estimating the FPC scores is likely to
be small.

Substantial bias can be caused by systematic missingness in the data. Liebl
and Rameseder (2019) review certain violations of the MCAR assumption and
motivate novel estimators for the mean and covariance structure for dense
incomplete curve settings. While the classical estimators for the mean and
covariance are consistent in regions of the domain with (virtually) no
missingness, they are prone to (severe) bias the stronger the violation from
the MCAR assumption and the fewer observations are available.

Estimation accuracy is also crucially affected by the coverage of the overall
domain, especially so for estimating the covariance operator. Only if a
sufficient number of observed curves overlap on the respective parts of their
observed domains can the corresponding regions of the covariance surface be
estimated reliably. For the setting of short observed domains (”functional
fragments”), Delaigle et al. (2020), Descary and Panaretos (2019) and Zhang
and Chen (2017) introduce conditions and approaches for consistently
estimating (parts of) the covariance surface.
— Software implementations
Some general FPCA methods are implemented in R packages fda (Ramsay et al.,
2020) and refund (Goldsmith et al., 2020). The multivariate FPCA approach of
Happ et al. (2019) is implemented in package MFPCA (Happ-Kurz, 2020); the PACE
algorithm of Yao et al. (2005) in fdapace (Carroll et al., 2020). The methods
outlined in Gertheiss et al. (2017) are available for the binary curves
setting in gfpca (Goldsmith, 2016), where the mixed regression in the two-step
approach is estimated with package gamm4 (Wood and Scheipl, 2020). Our
accompanying package registr (Wrobel and Bauer, 2021) implements the Gaussian
and binary curve GFPCA of Wrobel et al. (2019) as well as the two-step
approach of Gertheiss et al. (2017) for various exponential family
distributions.

##  3 Methods

As outlined, curves can have missing information at the beginning of their
domain (i.e., leading incompleteness), at the end of their domain (trailing
incompleteness), or both (full incompleteness). Our approach is able to handle
all three types of incompleteness for curves observed on potentially irregular
individual grids of evaluation points, without assuming Gaussianity of the
observed data.

We first introduce our registration approach and the approach for generalized
FPCA in full detail, and then present the main iterative algorithm to obtain a
solution where the registered curves are well represented by a low-rank GFPCA
basis. Potential identifiability issues and practical implications are
discussed at the end of this section. Computational details are given in
Appendix A1.

###  3.1 Registration for incomplete curves

We extend the likelihood-based framework for registering complete curves from
exponential family distributions of Wrobel et al. (2019). In the registration
step, the individual _chronological time domains_ $T_{i}^{*}$ are mapped onto
the registered _internal time domain_ $T$. This is achieved by estimating
inverse warping functions $h^{-1}_{i}$ that deform an unregistered curve
$Y_{i}(t_{i}^{*})$ toward a suitable template function $\mu_{i}(t)$ so that

$\displaystyle\mathsf{E}\left[Y_{i}\left(h^{-1}_{i}(t_{i}^{*})\right)|h^{-1}_{i}\right]$
$\displaystyle=\mu_{i}\left(t\right),$ (1) $\displaystyle\text{with}\ \ \ \
h^{-1}_{i}(t_{i}^{*})$
$\displaystyle=\boldsymbol{\Theta}_{h}(t^{*}_{i})\boldsymbol{\beta}_{i},$

with $Y_{i}(t)$ the registered curve. The inverse warping functions are
represented through a B-spline basis with design matrix
$\boldsymbol{\Theta}_{h}\in\mathbb{R}_{D_{i}\times K_{h}}$, $K_{h}$ basis
functions and a corresponding coefficient vector $\boldsymbol{\beta}_{i}$.

Given some distribution from the exponential family this yields the log-
likelihood for curve $i$:

$\ell\left(h^{-1}_{i}|y_{i},\mu_{i}\right)=\log\left(\prod_{j=1}^{D_{i}}f_{i,j}\left[y_{i}(t^{*}_{i,j})\right]\right),$
(2)

with $f_{i,j}(\cdot)$ the corresponding density with expected value
$\mu_{i}\left(h^{-1}_{i}(t^{*}_{i,j})\right)$ and observed vectors of function
evaluations $y_{i}(t^{*}_{ij})$, $j=1,\ldots,D_{i}$. We impose working
assumptions of mutual conditional independence across functions
$\left[Y_{i}\perp Y_{i^{\prime}}\right]|\mu_{i},\mu_{i^{\prime}}$ as well as
within functions $\left[Y_{i}(t_{ij})\perp Y_{i}(t_{ik})\right]|\mu_{i}$ .

— Constrained optimization
Warping functions must follow certain constraints so that they yield
reasonable transformations of the time domain. First, all warping functions
have to be strictly increasing to preserve the temporal order of a curve’s
measurements. Second, warping functions have to be domain-preserving with
regard to the maximal domain of the underlying process. We ensure both by
using a constrained optimization algorithm for the warping functions’ basis
coefficients (see Appendix A1).

— Circumventing the constraint of fixed time intervals
If all curves are observed over an identical time interval, domain
preservation requires that all warping functions map $t_{\min}^{*}$ and
$t_{\max}^{*}$ to themselves so that they begin and end on the diagonal line.
This assumption is made in most currently available registration procedures,
based on an implicit assumption that the process of interest (e.g., the growth
process of children) was observed from its very beginning up to its very end
for all subjects. For incomplete curves, however, forcing observed and
registered domain lengths to be identical is clearly unsuitable. We drop these
hard constraints on the warping functions’ basis coefficients and allow
warping functions to start and / or end at any point inside the overall time
domain. To avoid large deformations of the time domain that are not strongly
supported by the data, we penalize the total amount by which the registration
changes the duration of the (observed) time domain. In a setting with full
incompleteness, we use the following penalized log-likelihood for the
registration step:

$\displaystyle\ell_{\operatorname{pen}}\left(h^{-1}_{i}|y_{i},\mu_{i}\right)$
$\displaystyle=\ell\left(h^{-1}_{i}|y_{i},\mu_{i}\right)-\lambda\cdot
n_{i}\cdot\operatorname{pen}\left(h^{-1}_{i}\right),$ (3)
$\displaystyle\text{with}\ \ \ \ \ \
\operatorname{pen}\left(h^{-1}_{i}\right)$
$\displaystyle=\left(\left[h^{-1}_{i}(t_{\max,i}^{*})-h^{-1}_{i}(t_{\min,i}^{*})\right]-\left[t_{\max,i}^{*}-t_{\min,i}^{*}\right]\right)^{2}.$

For leading incompleteness with
$h^{-1}_{i}(t_{\max,i}^{*})=t_{\max,i}^{*}\;\forall\,i$, this simplifies to
$\operatorname{pen}\left(h^{-1}_{i}\right)=\left[h^{-1}_{i}(t_{\min,i}^{*})-t_{\min,i}^{*}\right]^{2}$
and for trailing incompleteness with
$h^{-1}_{i}(t_{\min,i}^{*})=t_{\min,i}^{*}\;\forall\,i$ to
$\operatorname{pen}\left(h^{-1}_{i}\right)=\left[h^{-1}_{i}(t_{\max,i}^{*})-t_{\max,i}^{*}\right]^{2}$.
In other words, the penalty for one-sided incompleteness represents the
squared distance of the respective endpoint of $h^{-1}_{i}$ to the diagonal.
The penalization parameter $\lambda$ controls how much overall time dilation
or compression the registration can perform and is scaled by the number of
measurements $n_{i}$ of curve $i$ to ensure that the impact of the
penalization relative to the likelihood is not affected by the number of
measurement points per function. Details on the choice of $\lambda$ are given
in Section 3.4.

###  3.2 Generalized Functional PCA for incomplete curves

We adapt the two-step approach of Gertheiss et al. (2017) to estimate a low-
rank GFPCA representation of the registered curves
$Y_{i}\left(t\right)=Y_{i}\left(h^{-1}_{i}(t_{i}^{*})\right)$. Following Hall
et al. (2008) and the groundwork of Yao et al. (2005), functional principal
components (FPCs) are estimated using a marginal, semiparametric method based
on assuming a latent Gaussian process $X_{i}(t)$ so that

$\displaystyle\mathsf{E}[Y_{i}(t)]$ $\displaystyle=\mu_{i}(t)=g[X_{i}(t)],$
(4) $\displaystyle X_{i}(t)$
$\displaystyle\approx\alpha(t)+\sum\limits_{k=1}^{K}c_{i,k}\cdot\psi_{k}(t),$

where each observed $Y_{i}(t)$ corresponds to the transformation of the latent
process $X_{i}(t)$ with some fixed response function $g(\cdot)$, and the
latent process itself can be decomposed into a smooth global mean $\alpha(t)$,
FPCs $\psi_{k}(t)$ with respective eigenvalues $\tau_{k}>0$, and FPC scores
$c_{i,k}\sim N(0,\tau_{k})$. With known $\psi_{k}(t)$, model (4) is a
generalized functional additive mixed model along the lines of Scheipl et al.
(2016a) with a smooth conditional latent mean function in a P-spline
representation, functional random effects in an FPC basis representation, and
functional random effect scores $c_{i,k}\sim N(0,\tau_{k})$. While the
derivation of the covariance approximation below assumes that $X_{i}(t)$ shows
only small variation around its mean, stronger variation only has “a modest
effect on the errors in individual predictions” (Hall et al., 2008, 4.2).

To derive the GFPCA solution, we first center the observed $Y_{i}(t)$ based on
their marginal mean $\mu_{Y}(t)=\mathsf{E}[Y_{i}(t)]$, estimated through a
simple smoother of all $(t_{ij},Y_{i}(t_{ij}))$-pairs via a generalized
additive model (GAM, Fahrmeir et al., 2013) with response function $g(\cdot)$
and the appropriate exponential family for the response. The covariance of the
latent process can then be approximated by

$\widehat{\mathsf{Cov}}\left[X_{i}(s),X_{i}(t)\right]\approx\frac{\hat{\sigma}_{Y}(s,t)}{g^{(1)}[\mu_{X}(s)]\cdot
g^{(1)}[\mu_{X}(t)]},$ (5)

with $\sigma_{Y}(s,t)=\mathsf{E}[Y_{c,i}(s)\cdot Y_{c,i}(t)]$ based on the
centered curves $Y_{c,i}(t)$, the marginal mean $\mu_{X}(t)$ estimated
accordingly to $\mu_{Y}(t)$ and $g^{(1)}(\cdot)$ the first derivative of the
response function. For given time points $s_{1}$ and $s_{2}$,
$\sigma_{Y}(s_{1},s_{2})$ is estimated as the mean of all pairwise products
$y_{c,i}(s_{1})\cdot y_{c,i}(s_{2})$. The estimated surface
$\hat{\sigma}_{Y}(s,t)$ is a smoothed version of $\sigma_{Y}(s,t)$ using a
bivariate tensor product P-spline basis (Fahrmeir et al., 2013). Since the
surface is expected to show some discontinuity along the diagonal this
smoothing step is performed under exclusion of the diagonal elements (c.f. Yao
et al., 2005). The FPCs $\psi_{k}(t)$ and their respective eigenvalues
$\tau_{k}$ are then estimated from the spectral decomposition of
$\widehat{\mathsf{Cov}}\left[X_{i}(t),X_{i}(s)\right]$. Our approach deviates
slightly from the method of Hall et al. (2008) it is based on: We mean-center
the data before taking their crossproducts instead of subtracting the
crossproduct of the estimated mean from the crossproducts of the data. In our
experience, this yields smoother estimates of the covariance surface which are
more amenable to a low-rank FPC representation.

###  3.3 Joint approach

We utilize the iterative algorithm of Wrobel et al. (2019) to combine the
outlined approaches for registration and GFPCA. Our aims are twofold: (i)
register all observed curves $Y_{i}(t_{i}^{*})$ to suitable template functions
and (ii) adequately represent the registered curves
$Y_{i}(t)=Y_{i}(h^{-1}_{i}(t_{i}^{*}))$ through a low-rank GFPCA basis. We
solve this problem by alternating the registration step (conditional on the
current GFPCA representations $\mu_{i}(t)$) and the GFPCA step (conditional on
the current estimates of the warping functions $h^{-1}_{i}$). The initial
registration step is performed with respect to a fixed common template
function ${\mu(t)}^{[0]}$ which has to be set by the user. Subsequent
iterations then use the low-rank GFPCA representations $\mu_{i}(t)$ as curve-
specific template functions. Full details on the iterative estimation are
given in Algorithm 1.

The number of FPCs in each iteration can be chosen based on the explained
proportion of variance. We adapt this criterion to account for peculiarities
of the covariance structure estimated with the two-step approach. Full details
are discussed at the end of Section 3.4.

Algorithm 1 Joint Registration & GFPCA

0: Observed curves $y_{i}(\boldsymbol{t}_{i}^{*})$; starting template
${\mu(t)}^{[0]}$; explained share of variance $\kappa_{\text{var}}$ of GFPCA
solution; convergence tolerance $\Delta_{h}$, iteration counter $q=0$.

1: Initial registration of observed curves $y_{i}(\boldsymbol{t}_{i}^{*})$ to
global initial template ${\mu(t)}^{[0]}$ to initialize
${\hat{h}^{-1}_{i}(t^{\star})}^{[0]}$;

2: while
$\sum_{i=1}^{N}\left(\sum_{j=1}^{D_{i}}\left[{\hat{h}^{-1}_{i}(t^{*}_{i,j})}^{[q]}-{\hat{h}^{-1}_{i}(t^{*}_{i,j})}^{[q-1]}\right]^{2}\right)>\Delta_{h}$
do

3: $q\to q+1$

4: Update GFPCA using registered curves
$y_{i}\left({\hat{h}^{-1}_{i}(\boldsymbol{t}_{i}^{*})}^{[q-1]}\right)$
(Section 3.2).

5: Re-estimate GFPCA representations ${\mu_{i}(t)}^{[q]}$ based on the first
$K^{[q]}$ FPCs that explain at least a share $\kappa_{\text{var}}$ of the
total variance (Model (4))

6: Update warping function estimates ${\hat{h}^{-1}_{i}(t^{\star})}^{[q]}$ by
re-registering observed curves $y_{i}(\boldsymbol{t}_{i}^{*})$ to
${\mu_{i}(t)}^{[q]}$.

7: end while

8: Final GFPCA estimation based on the registered curves
$y_{i}\left({\hat{h}^{-1}_{i}(\boldsymbol{t}_{i}^{*})}^{[q]}\right)$ to obtain
GFPCA representations $\mu_{i}(t)$ based on the first $K$ FPCs that explain at
least share $\kappa_{\text{var}}$ of the total variance.

###  3.4 Pitfalls and Practical Considerations

— Identifiability
A common issue in the separation of amplitude and phase variation is that
disentangling the two types of variation is an ill-posed problem in most
realistic settings. Structured variability in the curves can almost always be
attributed to either warpings of the time domain or superpositions of
principal components, or any combination of the two. Both Chakraborty and
Panaretos (2017) and Wagner and Kneip (2019) have shown that the general
registration problem has a unique solution only if the amplitude variation is
of rank 1, i.e. for FPC rank $K=1$. In practice, this non-identifiability can
be removed by introducing suitable inductive biases for estimates of the
warping and template functions through priors, penalties and/or limiting the
expressivity of model components, e.g. by choosing low-rank basis
representations. We assess the severity of this identifiability problem for
our method in a simulation study in Section 4. Note that the low-rank basis
representations of the warping functions we employ also seem to successfully
avoid the “pinching” problem (see e.g. Ramsay and Li, 1998, 4.2).

— Choice of the template function
As outlined above, the template function ${\mu(t)}^{[0]}$ for the initial
registration step in the joint estimation has to be set by the user. The
choice should be based on subject knowledge and can be crucial for obtaining
reasonable results (compare Appendix A7) and quick convergence in subsequent
iterations.

— Choice of the penalization parameter $\lambda$
Our registration approach controls the overall amount of compression or
dilation through the penalization parameter $\lambda$. The choice of $\lambda$
should be based on substantive knowledge so that estimated warping functions
represent realistic accelerations and/or decelerations of the observed
processes.

— Choice of the number of FPCs
The number of FPCs can be chosen based on the explained share of variance of
the low-dimensional FPC basis. In this regard, the two-step approach faces two
issues: First, since the spectral decomposition is applied to a smoothed
covariance surface and not the raw covariance of the data itself, the
“explained” share of variance is relative to this “structured” part of the
total observed variance. Second, based on our practical experience, spectral
decompositions of covariance surfaces often yield a large number of
subordinate FPCs which each explain only a very small amount of overall
variation, but jointly explain a relevant share. As we show in Appendix A6, it
can be argued that these subordinate FPCs often represent phase variation
rather than amplitude variation.

We deliberately avoid including such subordinate FPCs in the FPCA solution
since the FPCs in the joint approach should only represent the main amplitude
variation. Our goal is to find suitable template functions to register against
which don’t include modes of phase variation, i.e., the template functions do
not need to represent each individual registered curve with high fidelity.
Accordingly, we suggest a two-fold criterion for choosing the number of FPCs
based on our two-step approach: Choose as many FPCs as are needed to explain a
large portion (90%, by default) of the overall structured variation. However,
do not include such FPCs in the final solution that account for very little
variation (less than 2%, by default). In this way, we define a low-rank FPCA
representation for the template functions which might explain less than 90% of
the overall variation but does not include a multitude of subordinate modes of
(phase) variation.

##  4 Simulation Study

To assess the performance of our method and compare it to other established
approaches we perform a simulation study on both Gaussian and Gamma data,
motivated by our seismic application. We focus on the comparison of approaches
that jointly perform registration and FPCA and assess (i) their ability to
recover de-noised underlying curves, (ii) their performance in disentangling
and estimating the underlying amplitude and phase variation, and (iii) their
computational efficiency.

We compare our proposal (called “FGAMM” in the following) – combining two-step
FPCA with our (in)complete curve registration – with the earlier approach of
Wrobel et al. (2019) (“varEM”) – using an identical registration approach
combined with a variational EM-based FPCA – and the joint SRVF approach of
Tucker (2014) (Algorithm 4.1, “SRVF”) which combines the SRVF registration of
Srivastava et al. (2011) with the vertical fPCA introduced in Tucker et al.
(2013). The latter approach is only applied to complete curve settings since
the software implementation available at the time of writing (Tucker, 2020,
R-package fdasrvf) is not able to handle incomplete curves.

###  4.1 Simulation design

In each simulation setting, we first simulate $N=100$ complete curves on a
regular time grid on $[0,1]$ with length $D_{i}=50\ \forall i$ from model (4)
with FPC rank $K\in\\{1,3,4\\}$, with

  * •

mean function $\alpha(t)$ a Gaussian density function with $\mu=0.45$ and
$\sigma=0.2$,

  * •

eigenfunctions $\psi_{k}(t)$ as the $(k+1)$th orthonormal polynomial on
$[0,1]$,

  * •

mutually independent FPC scores $c_{i,k}\sim N(0,\tau_{k})$ and
$\boldsymbol{\tau}=1$ for $K=1$, $\boldsymbol{\tau}=(0.7,0.25,0.05)$ for
$K=3$, $\boldsymbol{\tau}=(0.4,0.3,0.2,0.1)$ for $K=4$,

  * •

Gaussian setting: $Y_{i}(t_{j})\sim
N\left(X_{i}(t_{j}),\sigma^{2}=0.03\right)$,

  * •

Gamma setting:
$Y_{i}(t_{j})\sim\Gamma\left(k=5,\theta=\frac{1}{5}\exp\left(X_{i}(t_{j})\right)\right)$,

with $X_{i}(t_{j})=\alpha(t)+\sum^{K}\psi_{k}(t)c_{i,k}$ the simulated
(latent) process.

Warping functions are simulated utilizing a B-spline basis using cubic splines
and three degrees of freedom. Their basis coefficients are drawn from a
uniform distribution over $[0,1]$ and cumulatively summed up to ensure
monotony. Three settings of (in)completeness are analyzed: Complete curves,
weak incompleteness and strong incompleteness. The latter two settings only
comprise trailing incompleteness. Weak incompleteness and strong
incompleteness are simulated by randomly drawing a cut-off time from a uniform
distribution over the last 40% and 70% of the time domain, respectively.

Regarding the correlation structure between the extents of (i) amplitude
variation, (ii) phase variation and (iii) incompleteness we analyze three
different settings. In the first, the three dimensions are mutually
uncorrelated. The second setting comprises a strong negative correlation
between amplitude and phase variation, shifting the peaks of curves with
larger amplitudes towards the beginning of the domain. The third setting
comprises a stronger positive correlation between amplitude and the amount of
incompleteness, resulting in stronger incompleteness for curves with lower
amplitudes.

Visualizations of the simulated data can be found in Appendix A3.1. For the
methods FGAMM and varEM, we use eight and four basis functions for the
estimation of the mean curve and the inverse warping functions, respectively.
In the Gaussian setting, the penalization parameter is set to $\lambda=0.025$.
In the Gamma setting, it is set to $1$ and $0.5$ for FGAMM (assuming a Gamma
distribution) and varEM (assuming a Gaussian distribution), respectively. The
observed overall mean curve is used as the initial template function. In the
FGAMM approach, the covariance surface is smoothed with ten marginal P-spline
basis functions. Since the implementation of the SRVF approach relies on the
curves being observed on a regular grid, the simulated curves are linearly
interpolated onto a regular grid. We perform 100 replications for each
simulation setting and method. The following results only cover the simulation
settings without correlation between phase, amplitude and incompleteness.
Unless noted otherwise, the results for the other simulation settings are
structurally similar (see Appendices A3.2 and A3.3).

— Adaptive estimation of the number of FPCs
The number of estimated FPCs was pre-set to the respective true simulated
amplitude rank. While our method includes adaptive, data-based estimation of
the number of FPCs (see Algorithm 1), we did not pursue this approach here
since this would jeopardize our ability to differentiate (i) its ability to
recover FPCs and their scores accurately and (ii) its ability to select a
suitable number of FPCs based on the data. Additional results based on the
more realistic use-case with adaptive estimation of the FPCs are given in
Appendices A3.4 and A3.5. The methods’ performances on the Gaussian simulation
settings with adaptively estimated FPC rank $K$ are structurally similar to
the ones with pre-specified rank. In the Gamma settings, while this is the
case for the estimated phase components, all methods struggle to recover the
correct number of FPCs and specifically the varEM approach performs worse in
terms of the estimation of amplitude variation.

###  4.2 Results

— Performance metrics
We base our method comparisons in Figures 3 and 4 on different performance
metrics, most based on the mean (integrated) squared error (MISE, MSE) for
functional and scalar estimates, respectively. Overall performance is
quantified using the difference between the simulated individual mean
structures (before adding random noise) and the respective representations
based on the final FPCA solution (measure $\operatorname{MISE}_{y}$). This
metric indicates how well the complete structured variation, i.e, the combined
phase and amplitude variation, of the observed data is recovered. The
performance regarding the separation and estimation of amplitude and phase
variation is quantified by (i) comparing the spans of the true and estimated
FPC bases with a measure introduced by Larsson and Villani (2001) and adapted
by Scheipl et al. (2016b) (amplitude variation, $\operatorname{LV}_{\psi}$)
and by (ii) comparing the true and estimated warping functions (phase
variation, $\operatorname{MISE}_{h}$). Following Scheipl et al. (2016b), the
measure $\operatorname{LV}_{\psi}$ quantifies the overlap of the spans of two
matrices $\boldsymbol{A}\in\mathbb{R}^{n\times p_{A}}$ and
$\boldsymbol{B}\in\mathbb{R}^{n\times p_{B}}$, $n>p_{A},p_{B}$:

$\operatorname{LV}_{\psi}(\boldsymbol{A},\boldsymbol{B})=\frac{1}{p_{A}}\cdot\operatorname{trace}\left(\boldsymbol{V}_{B}^{T}\boldsymbol{V}_{A}\boldsymbol{V}_{A}^{T}\boldsymbol{V}_{B}\right),$

with $\boldsymbol{V}_{Z}$, $Z\in\\{A,B\\}$, a matrix of the left singular
vectors of matrix $\boldsymbol{Z}$. We scale the measure by the dimension of
the true FPC basis $p_{A}$ to obtain a codomain of $[0,1]$ where value $1$
encodes perfect representation of the true amplitude variation space and $0$
represents completely orthogonal spans. In accordance with the other
performance measures we report $1-\operatorname{LV}_{\psi}$ so that smaller
values encode better performance. Note that $\operatorname{LV}_{\psi}$ cannot
be computed for the SRVF approach since that method is based on an FPCA of the
SRVF transforms of the original functions and does not yield orthonormal
eigenfunctions in the original function space. Finally, we compute the
estimation performance of the overall amount of time dilation or compression
by comparing the true and estimated domain lengths of the registered curves
($\operatorname{MSE}_{d}$).

— Results Gaussian settings
The results for the Gaussian settings are visualized in Figure 3. While
methods varEM and FGAMM do a good job in representing the underlying
structured variation ($\operatorname{MISE}_{y}$) and in estimating both
warping functions ($\operatorname{MISE}_{h}$) and original domain lengths
($\operatorname{MSE}_{d}$), amplitude variation ($\operatorname{LV}_{\psi}$)
is only estimated with higher accuracy for amplitude rank 1. Both FPCs and
warping functions are estimated more accurately if amplitude variation has
smaller rank.

Comparing the methods and focusing only on the complete curve settings (left
panels) for which it is applicable, the joint SRVF approach performs
consistently worse than the other two approaches. For
$\operatorname{MISE}_{y}$ and $\operatorname{MISE}_{h}$ the median performance
of FGAMM for amplitude rank 2–3 is better by 89% and 79% compared to SRVF,
respectively. The varEM approach performs slightly better than FGAMM for the
complete curve settings in terms of representing the overall variation, and
slightly worse in terms of recovering the space of amplitude variation.
Regarding the incomplete curve settings, the incomplete curve approaches
perform consistently best with respect to $\operatorname{MSE}_{d}$ and
$\operatorname{MISE}_{y}$. The estimated curve representations contain a much
higher share of the originally observed variation than is represented by
methods with assumed completeness. While the incomplete curve approaches
mostly perform better in terms of phase variation ($\operatorname{MISE}_{h}$
and $\operatorname{MSE}_{d}$), this is not consistently the case for the
estimation of amplitude variation ($\operatorname{LV}_{\psi}$). In summary,
among the evaluated incomplete curve methods, varEM performs somewhat better
than FGAMM, especially for representing the observed variation. We do not
observe a large drop in estimation performance between the settings with weak
and strong incompleteness.

Figure 3: Results for the simulation setting with Gaussian data and mutually
uncorrelated amplitude, phase and amount of incompleteness. All y scales are
$\log_{10}$ transformed.

— Results Gamma settings
For the Gamma settings, the varEM and SRVF approachesfall back on a
misspecified Gaussian or ”least squares” approach since neither are
implemented for Gamma data. FGAMM utilizes the appropriate Gamma likelihood
for both registration and GFPCA steps. The results are displayed in Figure 4.
All in all, for the setting without correlation between amplitude, phase and
incompleteness, the performance with regard to $\operatorname{MISE}_{h}$ and
$\operatorname{MSE}_{d}$ is similar to the Gaussian case. All methods show
consistently worse performance than in the Gaussian setting in terms of
$\operatorname{MISE}_{y}$ and $\operatorname{LV}_{\psi}$, also for small
amplitude ranks.

On complete data, the SRVF approach again performs worst in terms overall
representation and warping function estimation, with the FGAMM median
performance for amplitude rank 2–3 being better by 83% and 82%, respectively.
Regarding the estimation of the overall representation and the warping
functions, the (incomplete) FGAMM approach assuming the Gamma structure
performs consistently better than varEM. However, FGAMM performs worse in
recovering the FPC space, especially for the largest amplitude rank. For the
estimation of the amplitude structure and the domain lengths, varEM leads to
consistently better results. Comparing these results to the settings with weak
and strong incompleteness, the latter only show a structurally worse
estimation performance for $\operatorname{MSE}_{d}$.

Figure 4: Results for the simulation setting with Gamma data and uncorrelated
amplitude, phase and amount of incompleteness. All y scales are $\log_{10}$
transformed.

###  4.3 Runtime analysis

We evaluate the efficiency of the approaches on one simulation setting with a
Gaussian structure, amplitude rank 2–3 and complete curves. Only FGAMM is
additionally applied to the respective setting with Gamma data. The median
runtimes of each method, based on 20 runs, are visualized in Figure 5.

For the comparison, methods FGAMM and varEM (“varEM 2.1”) are based on
function register_fpca in version 2.1.5 of the registr package, which uses
methods from packages gamm4 (Wood and Scheipl, 2020) (v0.2.7) and lme4 (Bates
et al., 2015) (v1.1.26). These methods are compared to the old version of
registr (v1.0.0, based on gamm4 v0.2.6 and lme4 v1.1.23), which does not
contain the algorithmic improvements outlined in Appendix A1. We also compare
our methods to function align_fPCA of package fdasrvf (Tucker, 2020) (v1.9.4)
which implements the joint SRVF approach. All methods except version 1.0 of
package registr (“varEM 1.0”) were run in parallel mode using ten cores.

— Main findings
As can be seen in Figure 5, the optimized algorithm in varEM 2.1 is clearly
the most efficient method. For the setting with 50 measurements per curve and
3 000 curves (“$D_{i}=50,N=3000$”) varEM 2.1 (runtime 23 seconds) is on
average 86% faster than varEM 1.0 (159 seconds). The estimation of FGAMM is
computationally much more expensive. For the setting “$D_{i}=50,N=3000$” it
takes about 14 min, i.e., 37 times longer than varEM 2.1. Also, the runtime of
FGAMM scales quadratically in both the number of curves and the number of
measurements per curve. The efficiency of the SRVF approach lies between the
other methods for smaller samples. However, it becomes computationally
demanding for densely observed datasets with higher numbers of measurements
per curve.

Figure 5: Median runtimes for one setting of the simulation study with
amplitude rank 2-3 and no incompleteness, based on 20 runs for each parameter
combination. For the analysis, we vary the number of curves $N$ and the number
of measurements per curve $D_{i}$. Dashed curves are runtimes for FGAMM
(Gamma).

##  5 Application

###  5.1 Berkeley growth study

We compare FGAMM and varEM results on the well-known Berkeley growth data with
simulated strong full incompleteness as outlined in Section 1 and visualized
in Appendix A4.1. That is, we randomly remove both leading and trailing
segments of the curves, with starting points and endpoints drawn at random in
the first quarter and the last half of the time domain, respectively. Both
methods are then applied with and without the assumption of completely
observed curves, using a Gaussian likelihood and the same hyperparameters as
used for the simulation study. The number of FPCs to be used in each iteration
of the joint registration and FPCA algorithm was estimated adaptively, based
on the criterion outlined at the end of Section 3.

While the FGAMM approach chose 5 (assuming completeness) and 4 (incomplete)
FPCs, the varEM method chose 7 and 6 FPCs, respectively. In the comparison in
Figure 6, we focus on the first two FPCs estimated by each method. Results in
full detail are given in Appendix A4.2. Appendix A4.3 shows how the results of
the incomplete curve FGAMM approach changes when different values for the
penalization parameter $\lambda$ are used.

While the first two FPCs estimated by the methods with assumed incompleteness
show some differences, they represent similar main modes of amplitude
variation. The first FPC mainly represents variation at the very beginning of
the domain along with the information that the peak in growth in adolescent
age appears earlier on if the growth rate in the very first year was stronger.
The second FPC represents the information that if the initial growth rate was
higher, the peak in adolescent age is more attenuated.

These first two FPCs as estimated by the incomplete curve approaches differ
from the first FPCs estimated with assumed completeness. This is mainly due to
the fact that the “completeness-assumed” approaches are not able to adequately
align the structures observed in the last third of the domain (c.f. top row of
Figure 6). In this data setting, the incomplete curve approaches are clearly
better able to recover the underlying phase variation in the curves. In terms
of computation time, both FGAMM variants and the incomplete varEM take about a
minute, while varEM assuming completeness takes almost 2 minutes.

Figure 6: Observed curves with simulated incompleteness (top left pane),
registered curves (top row) and the first two estimated FPCs based on the
different approaches. The FPCs $\psi_{k}(t)$ are visualized by the overall
mean curve (solid line) plus (blue line) and minus (red line)
$2\cdot\sqrt{\hat{\tau}_{k}}\cdot\psi_{k}(t)$, with $\sqrt{\hat{\tau}_{k}}$
the standard deviation of the estimated scores for the $k$’th FPC.

###  5.2 Seismic ground motion propagation

We analyze a subset of seismological interest of the seismic data outlined in
Section 1, comprising 2 484 curves from various earthquakes simulated with
different physical parameters. Sepcifically, we use data from simulated quakes
characterized by (i) a sedimentary subsurface structure amplifying ground
shaking, (ii) a geologically well-oriented direction of the tectonic
background stress between 27∘ and 35∘, and (iii) the friction parameter of
critical linear slip weakening distance between 1.1m and 1.5m. We also
restrict the data to those seismograms most relevant for seismic hazard
assessment, which (i) lie in forward directivity direction (between cardinal
directions 280∘ and 342∘) to focus on wave propagations to the northwest in
the direction of the main rupture pulse, and (ii) with a hypocentral distance
shorter than 35km. Previous analyses show that the ground velocity curves we
study are primarily shaped by the hypocentral distance of the measurement
station and the _dynamic coefficient of friction_ which resembles the
frictional resistance of the geological fault during earthquake propagation
(Bauer et al., 2017). For our analysis, we focus on how these two parameters
and the topography of the evaluated region are associated with phase and
amplitude variation.

As shown in Figure 2, all curves are pre-processed by cutting off any leading
zero measurements below $0.01$, leading to the observed time domain
$t_{0}^{*}$ which – being the time since the first relevant absolute ground
velocity measurement – begins with the arrival time of seismic P-waves and
comprises trailing incompleteness only towards the end of the domain after
$t_{0}^{*}=23.5$ seconds. Since this induces a MAR structure, where short
observed domain lengths are caused by higher hypocentral distances (causing
later P-wave arrival times and smaller amplitudes), the results towards the
end of the domain must be interpreted with great care.

We apply the FGAMM approach assuming a Gamma structure and trailing
incompleteness, and using a similar parametrization as in the simulation
study. The mean curve of all observed curves was used as the template function
for the initial registration step. We use a penalization parameter of
$\lambda=0.004$ to discourage extreme distortions of the time domain.
Estimation of the joint approach took ten joint iterations and a runtime of
3:31h using a parallelized call for the registration steps with 5 cores. Two
FPCs were chosen based on the selection criterion outlined in Section 3 when
aiming to explain 95% of amplitude variation.

The two estimated FPCs along with the observed, registered and represented
curves are visualized in Figure 7. The full estimated warping functions are
shown in Appendix A5.1. The first FPC as the main mode of amplitude variation
represents the overall magnitude of the ground velocities, shaped by two
salient peaks that resemble the shaking caused by surface wave phase arrivals.
The second FPC represents a subsequent mode of variation and mainly shapes how
pronounced the initial peak is.

Figure 7: Lasagna plots of observed and registered curves and of the curves as
represented by the final GFPCA solution based on the FGAMM approach, on log10
scale (left pane). Curves are sorted by their maximum observed value. The FPCs
$\psi_{k}(t)$ are visualized by the overall mean curve (black line) plus (blue
line) and minus (red line) $2\cdot\sqrt{\hat{\tau}_{k}}\cdot\psi_{k}(t)$, with
$\sqrt{\hat{\tau}_{k}}$ the standard deviation of the estimated scores for the
$k$’th FPC (right pane).

The associations of phase and amplitude variation with the hypocentral
distance and the dynamic coefficient of friction are visualized in Figure 8.
Amplitude variation shows a very pronounced association structure with both
parameters. Focusing on the first FPC, ground velocities are overall stronger
the closer the measurement was taken to the hypocenter and the smaller the
dynamic friction. The strongest ground motion is observed at hypocentral
distances between 20 and 25km, caused by the nonlinear interaction of rupture
propagation and the radiated seismic wavefield with topography and the
subsurface structure. We find that in this region source effects (rupture
directivity) and seismic wave path effects (surface waves) unleash the most
energy. The second FPC’s scores show a somewhat similar association structure
but are more strongly shaped by the hypocentral distance. The highest scores
were estimated at around 25km of distance, representing the most pronounced
initial peak structure, especially in simulations with low dynamic friction
values.

Figure 8: Estimated phase and amplitude variation conditional on the
hypocentral distance of the virtual seismometer and the dynamic coefficient of
friction $\mu_{d}$ of the simulation. Phase variation and amplitude variation
are shown by displaying the mean of the overall time distortion after 5 and 20
seconds (left pane, with positive and negative values representing time
dilation and compression, respectively) and of the curves’ mean scores for the
FPCs shown in Figure 7, respectively.

Phase variation is also strongly associated with both evaluated parameters.
The estimated time distortions at time $t$ given by $\hat{h}^{-1}_{i}(t)-t$
for $t\in\\{5,20\\}$ show somewhat similar patterns to the association
structures of the FPC scores. This corroborates the structure displayed in
Figure 2 which indicates a strong coupling between amplitude and phase
variation since the initial peak is generally observed later for smaller
overall observed ground velocities (an effect known in seismology as
geometrical spreading). Stronger time distortion of the initial five seconds
was mostly estimated for medium-to-large friction values, which are
accompanied by smaller ground velocities. While these initial five seconds for
curves observed between 20 and 25km of hypocentral distance (shaped by a more
pronounced structure of the initial peak around $t_{0}=5$, see Figure 7) were
mainly compressed, for curves closer to and farther away from the hypocenter
(shaped by a less salient initial peak around $t_{0}=2$) they were mainly
dilated. The time distortion of the initial 20 seconds shows a very similar
structure to the scores for the first FPC. For curves with higher ground
velocities, these initial 20 seconds are mainly compressed, and mainly dilated
for lower ground velocities. Finally, the estimated inverse warping functions
more often tend to more extreme distortions for higher hypocentral distance
and higher dynamic friction values (see Appendix A5, Figure 20). This is due
to curves under these conditions showing very small ground velocities with a
less salient structure that is hard to align to the estimated template
functions.

As expected due to the dominant effects of source directivity and surface
waves in the evaluated region around the hypocenter, no structural association
of amplitude or phase variation with the local topography was detected (see
Appendix A5.2). The obtained results are geophysically plausible and in line
with previous analyses of the seismic experiments (Bauer et al., 2017).

##  6 Discussion

Incomplete data are very common in longitudinal settings but remain under-
discussed in many fields of functional data analysis. Our likelihood-based
approach for joint registration and generalized FPCA allows for analyzing
curves with leading, trailing or full incompleteness in the presence of
substantial phase variation and is able to handle non-Gaussian data. All
methods are implemented in the open-source R package registr.

Our simulation study results indicate that accounting for incompleteness
improves the performance in different data settings. While the FGAMM approach
shows some bias in the estimation of the underlying FPC structure in the Gamma
settings, its substantially better estimation of the warping functions leads
to improved overall performance in terms of the representation of the joint
phase and amplitude variation structure of the individual curves. Stronger
incompleteness does not seem to structurally harm the overall performance.
Applications to incomplete Berkeley growth curves and a seismic data setting
showcase the practical utility of our new approach.

— Comparison to SRVF-based approaches
In contrast to methods based on the SRVF framework, we do not utilize the
warping-invariant Fisher-Rao metric. Instead, our flexible penalized
likelihood-based approach allows for representing more complex structures of
variation in diverse non-Gaussian data situations and is backed by robust
optimization algorithms. While SRVF approaches rely on the availability of
functional derivatives evaluated on a common, regular grid and may struggle in
the presence of stronger (non-Gaussian) noise, this is generally not the case
for our method. We utilize a low-dimensional B-spline basis for the inverse
warping functions. In our applications, this seemed sufficient to avoid the
pinching problem. Extreme time distortions were only estimated for few seismic
curve outliers without a pronounced shape.

— Covariance estimation
One central topic for future research on GFPCA is a thorough evaluation of the
consistency and robustness of different covariance estimators. This comprises
questions like at what point in the estimation procedure smoothing and
centering (of the raw curves or the final covariance surface) should be
performed to obtain the best estimator. Covariance estimators should be
evaluated for common practical data settings entailing relevant non-Gaussian
noise in combination with small numbers of curves and measurements per curve
and different levels of their respective density over the domain.

— Computational efficiency
A practical constraint for the application of the evaluated methods remains
their computational efficiency in large-scale data settings. In this regard, a
promising strain of research are recently proposed neural network based
frameworks like Nunez et al. (2021) and Chen and Srivastava (2021) for
registration and Sarkar and Panaretos (2021) for covariance estimation.

## References

  * Bates et al. (2015) Bates, D., Mächler, M., Bolker, B., and Walker, S. (2015), “Fitting Linear Mixed-Effects Models Using lme4,” Journal of Statistical Software, 67, 1–48.
  * Bates et al. (2014) Bates, D., Mullen, K. M., Nash, J. C., and Varadhan, R. (2014), minqa: Derivative-free optimization algorithms by quadratic approximation, R package version 1.2.4.
  * Bauer (2021) Bauer, A. (2021), “bauer-alex/IncFunDatRegistration_supp: Full code and data supplement,” .
  * Bauer et al. (2017) Bauer, A., Scheipl, F., Küchenhoff, H., and Gabriel, A.-A. (2017), “Modeling spatio-temporal earthquake dynamics using generalized functional additive regressions,” in International Workshop on Statistical Modelling 2017.
  * Bauer et al. (2018) Bauer, A., Scheipl, F., Küchenhoff, H.— (2018), “An introduction to semiparametric function-on-scalar regression,” Statistical Modelling, 18, 346–364.
  * Bryner and Srivastava (2021) Bryner, D., and Srivastava, A. (2021), “Shape Analysis of Functional Data with Elastic Partial Matching,” arXiv preprint arXiv:2105.08604.
  * Cao et al. (2016) Cao, G., Wang, L., Li, Y., and Yang, L. (2016), “Oracle-efficient confidence envelopes for covariance functions in dense functional data,” Statistica Sinica, 26, 359–383.
  * Carroll et al. (2020) Carroll, C., Gajardo, A., Chen, Y., Dai, X., Fan, J., Hadjipantelis, P. Z., Han, K., Ji, H., Mueller, H.-G., and Wang, J.-L. (2020), fdapace: Functional Data Analysis and Empirical Dynamics, R package version 0.5.5.
  * Cederbaum et al. (2018) Cederbaum, J., Scheipl, F., and Greven, S. (2018), “Fast symmetric additive covariance smoothing,” Computational Statistics & Data Analysis, 120, 25–41.
  * Chakraborty and Panaretos (2017) Chakraborty, A., and Panaretos, V. M. (2017), “Functional Registration and Local Variations: Identifiability, Rank, and Tuning,” arXiv preprint arXiv:1702.03556.
  * Chen and Srivastava (2021) Chen, C., and Srivastava, A. (2021), “SrvfRegNet: Elastic Function Registration Using Deep Neural Networks,” in Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 4462–4471.
  * Cheng et al. (2016) Cheng, W., Dryden, I. L., Huang, X. et al. (2016), “Bayesian registration of functions and curves,” Bayesian Analysis, 11, 447–475.
  * Claeskens et al. (2010) Claeskens, G., Silverman, B. W., and Slaets, L. (2010), “A multiresolution approach to time warping achieved by a Bayesian prior–posterior transfer fitting strategy,” Journal of the Royal Statistical Society: Series B (Statistical Methodology), 72, 673–694.
  * Delaigle et al. (2020) Delaigle, A., Hall, P., Huang, W., and Kneip, A. (2020), “Estimating the covariance of fragmented and other related types of functional data,” Journal of the American Statistical Association, 1–19.
  * Descary and Panaretos (2019) Descary, M.-H., and Panaretos, V. M. (2019), “Recovering covariance from functional fragments,” Biometrika, 106, 145–160.
  * Fahrmeir et al. (2013) Fahrmeir, L., Kneib, T., Lang, S., and Marx, B. (2013), Regression: models, methods and applications, Springer Science & Business Media.
  * Gertheiss et al. (2017) Gertheiss, J., Goldsmith, J., and Staicu, A.-M. (2017), “A note on modeling sparse exponential-family functional response curves,” Computational statistics & data analysis, 105, 46–52.
  * Giorgino (2009) Giorgino, T. (2009), “Computing and Visualizing Dynamic Time Warping Alignments in R: The dtw Package,” Journal of Statistical Software, 31, 1–24.
  * Goldsmith (2016) Goldsmith, J. (2016), gfpca: Generalized Functional Principal Components Analysis, R package version 1.1.0. URL https://github.com/jeff-goldsmith/gfpca.
  * Goldsmith et al. (2020) Goldsmith, J., Scheipl, F., Huang, L., Wrobel, J., Di, C., Gellar, J., Harezlak, J., McLean, M. W., Swihart, B., Xiao, L., Crainiceanu, C., and Reiss, P. T. (2020), refund: Regression with Functional Data, R package version 0.1-23.
  * Goldsmith et al. (2015) Goldsmith, J., Zipunnikov, V., and Schrack, J. (2015), “Generalized multilevel function-on-scalar regression and principal component analysis,” Biometrics, 71, 344–353.
  * Guo et al. (2020) Guo, X., Kurtek, S., and Bharath, K. (2020), “Variograms for spatial functional data with phase variation,” arXiv preprint arXiv:2010.09578.
  * Hadjipantelis et al. (2015) Hadjipantelis, P. Z., Aston, J. A., Müller, H.-G., and Evans, J. P. (2015), “Unifying amplitude and phase analysis: A compositional data approach to functional multivariate mixed-effects modeling of Mandarin Chinese,” Journal of the American Statistical Association, 110, 545–559.
  * Hall et al. (2008) Hall, P., Müller, H.-G., and Yao, F. (2008), “Modelling sparse generalized longitudinal observations with latent Gaussian processes,” Journal of the Royal Statistical Society: Series B (Statistical Methodology), 70, 703–723.
  * Happ et al. (2019) Happ, C., Scheipl, F., Gabriel, A.-A., and Greven, S. (2019), “A general framework for multivariate functional principal component analysis of amplitude and phase variation,” Stat, 8, e220.
  * Happ-Kurz (2020) Happ-Kurz, C. (2020), MFPCA: Multivariate Functional Principal Component Analysis for Data Observed on Different Dimensional Domains, R package version 1.3-6.
  * Huang et al. (2014) Huang, H., Li, Y., and Guan, Y. (2014), “Joint modeling and clustering paired generalized longitudinal trajectories with application to cocaine abuse treatment data,” Journal of the American Statistical Association, 109, 1412–1424.
  * James et al. (2000) James, G. M., Hastie, T. J., and Sugar, C. A. (2000), “Principal component models for sparse functional data,” Biometrika, 87, 587–602.
  * Johnson (2020) Johnson, S. G. (2020), The NLopt nonlinear-optimization package, uRL https://github.com/stevengj/nlopt.
  * Kneip and Gasser (1992) Kneip, A., and Gasser, T. (1992), “Statistical tools to analyze data representing a sample of curves,” The Annals of Statistics, 1266–1305.
  * Kneip and Ramsay (2008) Kneip, A., and Ramsay, J. O. (2008), “Combining registration and fitting for functional models,” Journal of the American Statistical Association, 103, 1155–1165.
  * Kraus (2015) Kraus, D. (2015), “Components and completion of partially observed functional data,” Journal of the Royal Statistical Society: Series B: Statistical Methodology, 777–801.
  * Kurtek (2017) Kurtek, S. (2017), “A geometric approach to pairwise Bayesian alignment of functional data using importance sampling,” Electronic Journal of Statistics, 11, 502–531.
  * Larsson and Villani (2001) Larsson, R., and Villani, M. (2001), “A distance measure between cointegration spaces,” Economics Letters, 70, 21–27.
  * Learned-Miller (2005) Learned-Miller, E. G. (2005), “Data driven image models through continuous joint alignment,” IEEE Transactions on Pattern Analysis and Machine Intelligence, 28, 236–250.
  * Lee and Jung (2016) Lee, S., and Jung, S. (2016), “Combined analysis of amplitude and phase variations in functional data,” arXiv preprint arXiv:1603.01775.
  * Li et al. (2020) Li, C., Xiao, L., and Luo, S. (2020), “Fast covariance estimation for multivariate sparse functional data,” Stat, 9, e245.
  * Li and Guan (2014) Li, Y., and Guan, Y. (2014), “Functional principal component analysis of spatiotemporal point processes with applications in disease surveillance,” Journal of the American Statistical Association, 109, 1205–1215.
  * Liebl and Rameseder (2019) Liebl, D., and Rameseder, S. (2019), “Partially observed functional data: The case of systematically missing parts,” Computational Statistics & Data Analysis, 131, 104–115.
  * Lu et al. (2017) Lu, Y., Herbei, R., and Kurtek, S. (2017), “Bayesian registration of functions with a Gaussian process prior,” Journal of Computational and Graphical Statistics, 26, 894–904.
  * Marron et al. (2015) Marron, J. S., Ramsay, J. O., Sangalli, L. M., and Srivastava, A. (2015), “Functional data analysis of amplitude and phase variation,” Statistical Science, 468–484.
  * Mattar et al. (2009) Mattar, M. A., Ross, M. G., and Learned-Miller, E. G. (2009), “Nonparametric curve alignment,” in 2009 IEEE International Conference on Acoustics, Speech and Signal Processing, IEEE, pp. 3457–3460.
  * Matuk et al. (2019) Matuk, J., Bharath, K., Chkrebtii, O., and Kurtek, S. (2019), “Bayesian Framework for Simultaneous Registration and Estimation of Noisy, Sparse and Fragmented Functional Data,” arXiv preprint arXiv:1912.05125.
  * McDonnell et al. (2021) McDonnell, E. I., Zipunnikov, V., Schrack, J. A., Goldsmith, J., and Wrobel, J. (2021), “Registration of 24-hour accelerometric rest-activity profiles and its application to human chronotypes,” Biological Rhythm Research, 1–21.
  * Müller (2015) Müller, M. (2015), Fundamentals of music processing: Audio, analysis, algorithms, applications, Springer.
  * Nunez et al. (2021) Nunez, E., Lizarraga, A., and Joshi, S. H. (2021), “SrvfNet: A Generative Network for Unsupervised Multiple Diffeomorphic Functional Alignment,” in Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 4481–4489.
  * Panaretos and Zemel (2016) Panaretos, V. M., and Zemel, Y. (2016), “Amplitude and phase variation of point processes,” The Annals of Statistics, 44, 771–812.
  * Pelties et al. (2014) Pelties, C., Gabriel, A.-A., and Ampuero, J.-P. (2014), “Verification of an ADER-DG method for complex dynamic rupture problems,” Geoscientific Model Development, 847–866.
  * Powell (2009) Powell, M. J. (2009), “The BOBYQA algorithm for bound constrained optimization without derivatives,” Cambridge NA Report NA2009/06, University of Cambridge, Cambridge, 26–46.
  * R Core Team (2020) R Core Team (2020), R: A Language and Environment for Statistical Computing, R Foundation for Statistical Computing, Vienna, Austria, URL https://www.R-project.org/.
  * Ramsay et al. (2020) Ramsay, J. O., Graves, S., and Hooker, G. (2020), fda: Functional Data Analysis, R package version 5.1.9.
  * Ramsay and Li (1998) Ramsay, J. O., and Li, X. (1998), “Curve registration,” Journal of the Royal Statistical Society: Series B (Statistical Methodology), 60, 351–363.
  * Ramsay and Silverman (2005) Ramsay, J. O., and Silverman, B. W. (2005), Functional Data Analysis, Springer.
  * Rice and Wu (2001) Rice, J. A., and Wu, C. O. (2001), “Nonparametric mixed effects models for unequally sampled noisy curves,” Biometrics, 57, 253–259.
  * Sangalli et al. (2010) Sangalli, L. M., Secchi, P., Vantini, S., and Vitelli, V. (2010), “K-mean alignment for curve clustering,” Computational Statistics & Data Analysis, 54, 1219–1233.
  * Sarkar and Panaretos (2021) Sarkar, S., and Panaretos, V. M. (2021), “CovNet: Covariance Networks for Functional Data on Multidimensional Domains,” arXiv preprint arXiv:2104.05021.
  * Scheipl et al. (2016a) Scheipl, F., Gertheiss, J., Greven, S. et al. (2016a), “Generalized functional additive mixed models,” Electronic Journal of Statistics, 10, 1455–1492.
  * Scheipl et al. (2016b) Scheipl, F., Greven, S. et al. (2016b), “Identifiability in penalized function-on-function regression models,” Electronic Journal of Statistics, 10, 495–526.
  * Serban et al. (2013) Serban, N., Staicu, A.-M., and Carroll, R. J. (2013), “Multilevel cross-dependent binary longitudinal data,” Biometrics, 69, 903–913.
  * Srivastava et al. (2011) Srivastava, A., Wu, W., Kurtek, S., Klassen, E., and Marron, J. S. (2011), “Registration of functional data using Fisher-Rao metric,” arXiv preprint arXiv:1103.3817.
  * Stefanucci et al. (2018) Stefanucci, M., Sangalli, L. M., and Brutti, P. (2018), “PCA-based discrimination of partially observed functional data, with an application to AneuRisk65 data set,” Statistica Neerlandica, 72, 246–264.
  * Telesca and Inoue (2008) Telesca, D., and Inoue, L. Y. T. (2008), “Bayesian hierarchical curve registration,” Journal of the American Statistical Association, 103, 328–339.
  * Tipping and Bishop (1999) Tipping, M. E., and Bishop, C. M. (1999), “Probabilistic principal component analysis,” Journal of the Royal Statistical Society: Series B (Statistical Methodology), 61, 611–622.
  * Tormene et al. (2009) Tormene, P., Giorgino, T., Quaglini, S., and Stefanelli, M. (2009), “Matching incomplete time series with dynamic time warping: an algorithm and an application to post-stroke rehabilitation,” Artificial intelligence in medicine, 45, 11–34.
  * Tucker (2014) Tucker, J. D. (2014), “Functional component analysis and regression using elastic methods,” Ph.D. thesis, The Florida State University.
  * Tucker (2020) — (2020), fdasrvf: Elastic Functional Data Analysis, R package version 1.9.4.
  * Tucker et al. (2021) Tucker, J. D., Shand, L., and Chowdhary, K. (2021), “Multimodal Bayesian registration of noisy functions using Hamiltonian Monte Carlo,” Computational Statistics & Data Analysis, 107298.
  * Tucker et al. (2013) Tucker, J. D., Wu, W., and Srivastava, A. (2013), “Generative models for functional data using phase and amplitude separation,” Computational Statistics & Data Analysis, 61, 50–66.
  * Uphoff et al. (2017) Uphoff, C., Rettenberger, S., Bader, M., Madden, E. H., Ulrich, T., Wollherr, S., and Gabriel, A.-A. (2017), “Extreme scale multi-physics simulations of the tsunamigenic 2004 sumatra megathrust earthquake,” in Proceedings of the international conference for high performance computing, networking, storage and analysis, pp. 1–16.
  * van der Linde (2009) van der Linde, A. (2009), “A Bayesian latent variable approach to functional principal components analysis with binary and count data,” AStA Advances in Statistical Analysis, 93, 307–333.
  * Vitelli (2019) Vitelli, V. (2019), “A novel framework for joint sparse clustering and alignment of functional data,” arXiv preprint arXiv:1912.00687.
  * Wagner (2020) Wagner, H. (2020), Nonparametric Registration to Low Dimensional Function Spaces, URL https://github.com/heikowagner/Nonparametric-Registration-to-Low-Dimensional-Function-Spaces.
  * Wagner and Kneip (2019) Wagner, H., and Kneip, A. (2019), “Nonparametric registration to low-dimensional function spaces,” Computational Statistics & Data Analysis, 138, 49–63.
  * Wang et al. (2016) Wang, J.-L., Chiou, J.-M., and Müller, H.-G. (2016), “Functional data analysis,” Annual Review of Statistics and Its Application, 3, 257–295.
  * Wood and Scheipl (2020) Wood, S., and Scheipl, F. (2020), gamm4: Generalized Additive Mixed Models using ’mgcv’ and ’lme4’, R package version 0.2-6. URL https://CRAN.R-project.org/package=gamm4.
  * Wood (2017) Wood, S. N. (2017), Generalized additive models: an introduction with R, CRC press.
  * Wood et al. (2017) Wood, S. N., Li, Z., Shaddick, G., and Augustin, N. H. (2017), “Generalized additive models for gigadata: modeling the UK black smoke network daily data,” Journal of the American Statistical Association, 112, 1199–1210.
  * Wrobel (2018) Wrobel, J. (2018), “register: Registration for Exponential Family Functional Data,” Journal of Open Source Software, 3, 557.
  * Wrobel and Bauer (2021) Wrobel, J., and Bauer, A. (2021), “registr 2.0: Incomplete Curve Registration for Exponential Family Functional Data,” Journal of Open Source Software, 6, 2964.
  * Wrobel et al. (2021) Wrobel, J., Muschelli, J., and Leroux, A. (2021), “Diurnal Physical Activity Patterns across Ages in a Large UK Based Cohort: The UK Biobank Study,” Sensors, 21, 1545.
  * Wrobel et al. (2019) Wrobel, J., Zipunnikov, V., Schrack, J., and Goldsmith, J. (2019), “Registration for exponential family functional data,” Biometrics, 75, 48–57.
  * Wu and Srivastava (2014) Wu, W., and Srivastava, A. (2014), “Analysis of spike train data: Alignment and comparisons using the extended Fisher-Rao metric,” Electronic Journal of Statistics, 8, 1776–1785.
  * Yao et al. (2005) Yao, F., Müller, H.-G., and Wang, J.-L. (2005), “Functional data analysis for sparse longitudinal data,” Journal of the American statistical association, 100, 577–590.
  * Ypma and Johnson (2020) Ypma, J., and Johnson, S. G. (2020), nloptr, R package version 1.2.2.2.
  * Zhang and Chen (2017) Zhang, A., and Chen, K. (2017), “Nonparametric covariance estimation for mixed longitudinal studies, with applications in midlife women’s health,” arXiv preprint arXiv:1711.00101.
  * Zhou et al. (2008) Zhou, L., Huang, J. Z., and Carroll, R. J. (2008), “Joint modelling of paired sparse functional data using principal components,” Biometrika, 95, 601–619.

Appendix

##  A1 Computational Details

We implemented our approach in the package registr (Wrobel and Bauer, 2021)
for the statistical open-source software R (R Core Team, 2020). The registr
package allows for the estimation of the joint registration and GFPCA approach
both for complete and incomplete curves. All three types of incompleteness
(leading, trailing and full incompleteness) and irregular grids are supported.
Additional to the methods outlined in this work the package comprises the
methods of Wrobel et al. (2019). Several exponential family distributions are
available. In the following, we give details on some computational aspects of
our method.

###  A1.1 Registration

The registration codebase builds on the implementation outlined in Wrobel et
al. (2019) and Wrobel (2018). We extended the methods by allowing the observed
curves to be incomplete. Since the estimation of warping functions in the
registration step is performed separately for each curve, we added the option
of a parallelized call over the individual curves.

Constrained optimization for the spline coefficients representing the warpings
is performed with function constrOptim() by inducing linear inequality
constraints of the form

$\boldsymbol{u}_{i}\cdot\boldsymbol{\beta}_{i}-\boldsymbol{c}_{i}\geq 0,$

with parameter vector $\boldsymbol{\beta}_{i}$ and constraints given by matrix
$\boldsymbol{u}_{i}$ and vector $\boldsymbol{c}_{i}$. Further details on the
constraint matrices are given in Appendix A2. Alternative optimization
algorithms from the NLopt library (Johnson, 2020) and made available by
package nloptr (Ypma and Johnson, 2020) were evaluated as well, but did not
improve the overall results or lead to a more efficient estimation.

###  A1.2 Generalized Functional Principal Component Analysis

As outlined in Section 3.2, our adaptation of the two-step GFPCA approach of
Gertheiss et al. (2017) is based on an additive regression model with random
intercept terms for the individual FPCs. We build on robust and highly
efficient software for these kinds of models, available in packages gamm4
(Wood and Scheipl, 2020) and lme4 (Bates et al., 2015). The algorithms of lme4
are highly efficient for estimating models with random intercept terms with
several thousand individual categories. The estimation of the marginal mean of
the process $X_{i}(t)$ in the case of very large data with $>100\,000$ rows is
performed with the discretization-based estimation algorithm of function
mgcv::bam (Wood et al., 2017) rather than the estimation algorithm of
mgcv::gam (Wood, 2017).

Our implementation of the two-step GFPCA approach of Gertheiss et al. (2017)
is based on their accompanying package gfpca (Goldsmith, 2016). Additionally
we made several adjustments to their codebase to improve overall efficiency:
First, while functions lmer() and glmer() from the lme4 package default to the
optimization routine implemented in function bobyqa (package minqa, Bates et
al., 2014), we make use of the more efficient optimizer NLOPT_LN_BOBYQA from
the NLopt library (Johnson, 2020) as described in Powell (2009).

Second, we tackle one major issue in the building of the covariance structure.
In principle, the covariance matrix comprises the pairwise covariances between
all unique observed time points per functional datum $y_{i}(t)$. In real data
situations with highly irregular grids, the number of unique combinations of
time points can explode in size even for settings with a relatively low number
of curves. We utilize a binning strategy to handle this problem. Before
building the covariance matrix, we round the vector of observed time points to
$k$ significant digits. E.g., $k=3$ then leads to at most $1000^{2}$ unique
combinations and a covariance matrix with maximal size $1000\times 1000$.
Similar to the estimation of the marginal mean of $X_{i}(t)$, the smoothing of
the covariance surface is performed with the discretization-based estimation
algorithm mgcv::bam rather than mgcv::gam if the crossproduct matrix comprises
$>100\,000$ elements.

Third, we updated the codebase of gamm4 to make the initial construction of
the random effect model matrices much more efficient by fully exploiting their
sparse structure. Our patched version is currently available on GitHub
(https://github.com/r-gam/gamm4) and will in future be integrated into the
main codebase of the gamm4 package.

###  A1.3 Joint approach

To make the overall algorithm more efficient, we introduce two major changes.
First, all intermediate iterations regarding the GFPCA step, apart from the
very first and the very last one, are performed with less accuracy. Above all
else, we use larger tolerance values and a simple Laplace approximation to the
GLMM likelihood (i.e., nAGQ = 0 and nAGQinitStep = FALSE in function
gamm4::gamm4) for these iterations. Secondly, we use the solution of the
previous GFPCA step as starting values for the subsequent GFPCA step.

##  A2 Constraint Matrices for constrOptim()

As outlined in Appendix A1, we estimate the warping functions using function
constrOptim(). In the estimation step for one warping function, the parameter
vector is constrained s.t. the resulting warping function is monotone and does
not exceed the overall time domain $[t_{min},t_{max}]$.

In the following the constraint matrices are listed for the different settings
of (in)completeness and assuming a parameter vector of length $p$:

$\boldsymbol{\beta}_{i}=\left(\begin{array}[]{c}\beta_{i1}\\\ \beta_{i2}\\\
\vdots\\\ \beta_{ip}\end{array}\right)\in\mathbb{R}_{p\times 1}$

###  A2.1 Complete curve setting

When all curves were observed completely – i.e. the underlying processes of
interest were all observed from the beginning until the end – warping
functions can typically be assumed to start and end on the diagonal, since
each process is completely observed in its observation interval
$[t^{*}_{min,i},t^{*}_{max,i}]\subset[t_{min},t_{max}]$.

Assuming that both the starting point and the endpoint lie on the diagonal, we
set $\beta_{i1}=t^{*}_{min,i}$ and $\beta_{ip}=t^{*}_{max,i}$ and only perform
the estimation for

$\left(\begin{array}[]{c}\beta_{i2}\\\ \beta_{i3}\\\ \vdots\\\
\beta_{i(p-1)}\end{array}\right)\in\mathbb{R}_{(p-2)\times 1}$

This results in the following constraint matrices, that allow a mapping from
the observed domain $[t^{*}_{min,i},t^{*}_{max,i}]$ to the domain itself
$[t^{*}_{min,i},t^{*}_{max,i}]\subset[t_{min},t_{max}]$:

$\displaystyle\boldsymbol{u}_{i}$
$\displaystyle=\left(\begin{array}[]{cccccccc}1&0&0&0&\ldots&0&0&0\\\
-1&1&0&0&\ldots&0&0&0\\\ 0&-1&1&0&\ldots&0&0&0\\\
\vdots&\ddots&\ddots&\ddots&\ddots&\ddots&\ddots&\ddots\\\
0&0&0&0&\ldots&0&-1&1\\\
0&0&0&0&\ldots&0&0&-1\end{array}\right)\in\mathbb{R}_{(p-1)\times(p-2)}$
$\displaystyle\boldsymbol{c}_{i}$
$\displaystyle=\left(\begin{array}[]{c}t^{*}_{min,i}\\\ 0\\\ 0\\\ \vdots\\\
0\\\ -1\cdot t^{*}_{max,i}\end{array}\right)\in\mathbb{R}_{(p-1)\times 1}$

###  A2.2 Leading incompleteness only

In the case of leading incompleteness – i.e. the underlying processes of
interest were all observed until their very end but not necessarily starting
from their beginning – warping functions can typically be assumed to end on
the diagonal, s.t. one assumes $\beta_{ip}=t^{*}_{max,i}$ to let the warping
functions end at the last observed time point $t^{*}_{max,i}$. The estimation
is then performed for the remaining parameter vector

$\left(\begin{array}[]{c}\beta_{i1}\\\ \beta_{i3}\\\ \vdots\\\
\beta_{i(p-1)}\end{array}\right)\in\mathbb{R}_{(p-1)\times 1}$

This results in the following constraint matrices, that allow a mapping from
the observed domain $[t^{*}_{min,i},t^{*}_{max,i}]$ to the domain
$[t_{min},t^{*}_{max,i}]\subset[t_{min},t_{max}]$:

$\displaystyle\boldsymbol{u}_{i}$
$\displaystyle=\left(\begin{array}[]{cccccccc}1&0&0&0&\ldots&0&0&0\\\
-1&1&0&0&\ldots&0&0&0\\\ 0&-1&1&0&\ldots&0&0&0\\\
\vdots&\ddots&\ddots&\ddots&\ddots&\ddots&\ddots&\ddots\\\
0&0&0&0&\ldots&0&-1&1\\\
0&0&0&0&\ldots&0&0&-1\end{array}\right)\in\mathbb{R}_{p\times(p-1)}$
$\displaystyle\boldsymbol{c}_{i}$
$\displaystyle=\left(\begin{array}[]{c}t_{min}\\\ 0\\\ 0\\\ \vdots\\\ 0\\\
-1\cdot t^{*}_{max,i}\end{array}\right)\in\mathbb{R}_{p\times 1}$

###  A2.3 Trailing incompleteness only

In the case of trailing incompleteness – i.e. the underlying processes of
interest were all observed from the beginning but not necessarily until their
very end – warping functions can typically be assumed to start on the
diagonal, s.t. one assumes $\beta_{i1}=t^{*}_{min,i}$ to let the warping
functions start at the first observed time point $t^{*}_{min,i}$. The
estimation is then performed for the remaining parameter vector

$\left(\begin{array}[]{c}\beta_{i2}\\\ \beta_{i3}\\\ \vdots\\\
\beta_{ip}\end{array}\right)\in\mathbb{R}_{(p-1)\times 1}$

This results in the following constraint matrices, that allow a mapping from
the observed domain $[t^{*}_{min,i},t^{*}_{max,i}]$ to the domain
$[t^{*}_{min,i},t_{max}]\subset[t_{min},t_{max}]$:

$\displaystyle\boldsymbol{u}_{i}$ identical to the version for leading
incompleteness $\displaystyle\boldsymbol{c}_{i}$
$\displaystyle=\left(\begin{array}[]{c}t^{*}_{min,i}\\\ 0\\\ 0\\\ \vdots\\\
0\\\ -1\cdot t_{max}\end{array}\right)\in\mathbb{R}_{p\times 1}$

###  A2.4 Leading and trailing incompleteness

In the case of both leading and trailing incompleteness – i.e. the underlying
processes of interest were neither necessarily observed from their very
beginnings nor to their very ends – warping functions can typically only be
assumed to map the observed domains $[t^{*}_{min,i},t^{*}_{max,i}]$ to the
overall domain $[t_{min},t_{max}]$.

This results in the following constraint matrices:

$\displaystyle\boldsymbol{u}_{i}$
$\displaystyle=\left(\begin{array}[]{cccccccc}1&0&0&0&\ldots&0&0&0\\\
-1&1&0&0&\ldots&0&0&0\\\ 0&-1&1&0&\ldots&0&0&0\\\
\vdots&\ddots&\ddots&\ddots&\ddots&\ddots&\ddots&\ddots\\\
0&0&0&0&\ldots&0&-1&1\\\
0&0&0&0&\ldots&0&0&-1\end{array}\right)\in\mathbb{R}_{(p+1)\times p}$
$\displaystyle\boldsymbol{c}_{i}$
$\displaystyle=\left(\begin{array}[]{c}t_{min}\\\ 0\\\ 0\\\ \vdots\\\ 0\\\
-1\cdot t_{max}\end{array}\right)\in\mathbb{R}_{(p+1)\times 1}$

##  A3 Simulation study

###  A3.1 Simulation setting

This subsection contains figures for all relevant components of the curves
simulated in the simulation study.

####  A3.1.1 Distribution of the data

Figure 9: Structure of simulated Gaussian and Gamma data (top and bottom row,
respectively), before adding amplitude variation

####  A3.1.2 Rank of amplitude variation

Figure 10: Simulated eigenfunctions / functional principal components (FPCs),
visualized by adding and subtracting them from a some mean curve (black line)

Figure 11: Simulated curves with Gaussian structure, including amplitude
variation and random warping.

####  A3.1.3 Strength of incompleteness

Figure 12: Simulated curves with Gaussian structure and different strengths of
incompleteness.

####  A3.1.4 Correlation structure

Figure 13: Simulated curves with Gaussian structure and correlated amplitude
and phase variation. Figure 14: Simulated curves with Gaussian structure and
correlated amplitude variation and amount of incompleteness.

###  A3.2 Simulation results – Gaussian with correlation structure

Figure 15: Results for the simulation setting with Gaussian data and a
correlation between amplitude and phase.  Figure 16: Results for the
simulation setting with Gaussian data and a correlation between amplitude and
the amount of incompleteness.

###  A3.3 Simulation results – Gamma with correlation structure

Figure 17: Results for the simulation setting with Gamma data and a
correlation between amplitude and phase.  Figure 18: Results for the
simulation setting with Gamma data and a correlation between amplitude and the
amount of incompleteness.

###  A3.4 Simulation results – Gaussian with adaptive FPC estimation

In contrast to the previous settings of the simulation study, the number of
Functional Principal Components (FPCs) in the following settings is not fixed
to the simulated rank of amplitude variation. Instead, in each iterative FPCA
step (i) the varEM approach uses as many FPCs as are needed to explain 90% of
the overall amplitude variation, and (ii) the FGAMM approach uses as many FPCs
as are needed to explain 90% of the overall amplitude variation, while
dropping such FPCs that explain $<2\%$ of the variation (see criterion
outlined at the end of Section 3).

For the varEM approach, the explained share of variance and accordingly the
number of FPCs in each iteration is estimated before the main iteration’s
estimation step by once running the FPCA with 20 FPCs and correspondingly 20
B-spline basis functions to represent the FPC basis. Doing so, we approximate
the overall variance in the varEM approach with the variance represented by
this FPC basis with 20 FPCs. In contrast to the simulation results in the main
part of our paper, we accordingly use 20 instead of eight basis functions for
the estimation of the FPC basis in the varEM approach.

Note that the third and fourth FPC in the simulation settings with amplitude
rank 2–3 and 3–4 only explain 5% and 10% of the overall amplitude variation,
respectively (see Section 4.1). Accordingly, it is not unreasonable if fewer
than 3 and 4 FPCs are chosen based on the $\geq 90\%$ criterion, respectively.

Figure 19: Results for the simulation setting with Gaussian data and mutually
uncorrelated amplitude, phase and amount of incompleteness, where the number
of FPCs was adaptively estimated. All y scales are $\log_{10}$ transformed.
Figure 20: Results for the simulation setting with Gaussian data and a
correlation between amplitude and phase, where the number of FPCs was
adaptively estimated.  Figure 21: Results for the simulation setting with
Gaussian data and a correlation between amplitude and the amount of
incompleteness, where the number of FPCs was adaptively estimated.

###  A3.5 Simulation results – Gamma with adaptive FPC estimation

Note our remarks at the beginning of Appendix A3.4. The only difference to the
Gaussian setting is that the FGAMM approach assumes a Gamma distribution
instead of a Gaussian structure.

Figure 22: Results for the simulation setting with Gamma data and mutually
uncorrelated amplitude, phase and amount of incompleteness, where the number
of FPCs was adaptively estimated.  Figure 23: Results for the simulation
setting with Gamma data and a correlation between amplitude and phase, where
the number of FPCs was adaptively estimated.  Figure 24: Results for the
simulation setting with Gamma data and a correlation between amplitude and the
amount of incompleteness, where the number of FPCs was adaptively estimated.

##  A4 Berkeley application

###  A4.1 Curves with simulated incompleteness

Figure 25: Lasagna plot of observed curves (left pane) and curves with
simulated incompleteness (right) for the first derivative of the Berkeley
child growth data.

###  A4.2 Detailed results

Figure 26: Observed curves (top left pane), registered curves (top row),
estimated inverse warping functions (second row) and the first four estimated
FPCs based on the different approaches. The FPCs $\psi_{k}(t)$ are visualized
by displaying the overall mean curve (solid line) plus (dashed line, $+$) and
minus (dotted line, $-$) $x\cdot\psi_{k}(t)$, with $x$ twice the standard
deviation of the individual FPC’s scores.

###  A4.3 Varying the penalization parameter $\lambda$

The results based on different $\lambda$ values are shown in Figure 27. The
example is based on the Berkeley data discussed in Section 5.

While the overall domain dilation of the warping functions is not penalized
with $\lambda=0$, this is the case the higher the penalization parameter
$\lambda$ is chosen. With value $\lambda=1$ the penalization is strong enough
to cause all registered domain lengths to be (quasi) identical to the observed
domain lengths.

Figure 27: Results for varying values of the penalization parameter $\lambda$
after joint registration and Gaussian FPCA with the FGAMM approach. The
graphic shows spaghetti plots of the registered curves (first row), estimated
warping functions (second row) and the difference between the observed domain
lengths and the registered domain lengths (bottom row).

##  A5 Seismic application

###  A5.1 Estimated inverse warping functions

Figure 28: Estimated inverse warping functions displayed against the
hypocentral distance of the seismometers (x-axis) and the dynamic coefficient
of friction $\mu_{d}$ (y-axis). In each panel, a solid blue curve marks the
mean curve based on all respective warping functions.

###  A5.2 Estimated phase and amplitude variation over space

Figure 29: Estimated phase and amplitude variation visualized over the
evaluated region. Phase variation and amplitude variation are shown by
displaying the mean of the overall time distortion after 5 and 20 seconds
(left pane, with positive and negative values representing time dilation and
compression, respectively) and of the curves’ mean scores for the FPCs,
respectively. The right plot shows the topography of the region, with the
epicenter marked as red dot. The grey dashed lines mark the distance to the
epicenter in 5km steps.

##  A6 GFPCA: Structure of subordinate FPCs

This section evaluates one Gaussian data setting from the simulation study to
showcase the issue of subordinate functional principal components (FPCs) which

  1. 1.

individually explain a very small share of the overall amplitude variation,
but jointly explain a relevant share, and

  2. 2.

often tend to represent phase variation rather than amplitude variation.

For this evaluation, 100 curves are simulated similarly to the simulation
setting with complete curves, Gaussian noise, amplitude rank 2–3 and no
correlation between amplitude variation and phase variation. The only
differences to the simulation study are the following:

  * •

the curves are not randomly warped,

  * •

a regular time grid with length 100 is used.

The simulated curves are visualized in Figure 30. We estimate a solution with
20 FPCs with the two-step approach. These first 20 FPCs and their explained
shares of variance are visualized in Figure 31.

Figure 30: Lasagna plot of the simulated curves.  Figure 31: Visualization of
the first 20 FPCs including their percentage of explained variance (PVE).

##  A7 Choosing the initial template function

To check how much the results of the joint registration and GFPCA approach
vary based on different template functions for the initial registration step,
we run the application on the Berkeley data (from Section 5.1) with four
different template functions:

  * •

Template 1: Overall mean curve (similar to the application in the main paper)

  * •

Template 2: Curve where the main peak in the second half of the domain is not
very salient and occurrs quite early on

  * •

Template 3: Curve where the main peak occurs a bit later on and is a bit more
salient

  * •

Template 4: Curve where the main peak is even more salient

The template functions and the results of the application of the FGAMM
approach to the data can be found in the following Figure.

Figure 32: Results of the FGAMM approach based on the different initial
template functions (one column per template function). The rows contain the
observed curves with the template function in blue (first row), the registered
curves (second row), the estimated inverse warping functions (third row) and
the first two estimated FPCs (last two rows). The FPCs $\psi_{k}(t)$ are
visualized by the overall mean curve (solid line) plus (blue line) and minus
(red line) $2\cdot\sqrt{\hat{\tau}_{k}}\cdot\psi_{k}(t)$, with
$\sqrt{\hat{\tau}_{k}}$ the standard deviation of the estimated scores for the
$k$’th FPC.