# Multivariate Functional Additive Mixed Models

Alexander Volkmann  Chair of Statistics, School of Business and Economics,
Humboldt-Universität zu Berlin, Germany  Almond Stöcker  Chair of Statistics,
School of Business and Economics, Humboldt-Universität zu Berlin, Germany
Fabian Scheipl  Department of Statistics, Ludwig-Maximilians-Universität
München, Germany  Sonja Greven  Chair of Statistics, School of Business and
Economics, Humboldt-Universität zu Berlin, Germany

###### Abstract

Multivariate functional data can be intrinsically multivariate like movement
trajectories in 2D or complementary like precipitation, temperature, and wind
speeds over time at a given weather station. We propose a multivariate
functional additive mixed model (multiFAMM) and show its application to both
data situations using examples from sports science (movement trajectories of
snooker players) and phonetic science (acoustic signals and articulation of
consonants). The approach includes linear and nonlinear covariate effects and
models the dependency structure between the dimensions of the responses using
multivariate functional principal component analysis. Multivariate functional
random intercepts capture both the auto-correlation within a given function
and cross-correlations between the multivariate functional dimensions. They
also allow us to model between-function correlations as induced by e.g.
repeated measurements or crossed study designs. Modeling the dependency
structure between the dimensions can generate additional insight into the
properties of the multivariate functional process, improves the estimation of
random effects, and yields corrected confidence bands for covariate effects.
Extensive simulation studies indicate that a multivariate modeling approach is
more parsimonious than fitting independent univariate models to the data while
maintaining or improving model fit.

Keywords: functional additive mixed model; multivariate functional principal
components; multivariate functional data; snooker trajectories; speech
production

##  1 Introduction

With the technological advances seen in recent years, functional data sets are
increasingly multivariate. They can be multivariate with respect to the domain
of a function, its codomain, or both. Here, we focus on multivariate functions
with a one-dimensional domain
$\bm{f}=(f^{(1)},...,f^{(D)})\colon\mathcal{I}\subset\mathbb{R}\to\mathbb{R}^{D}$
with square-integrable components $f^{(d)}\in L^{2}(\mathcal{I}),d=1,...,D$.
For this type of data, we can distinguish two subclasses: One has
interpretable separate dimensions and can be seen as several complementary
modes of a common phenomenon (“multimodal” data, cf. Uludağ and Roebroeck;
2014) as in the analysis of acoustic signals and articulation processes in
speech production in one of our data examples. The codomain then simply is the
Cartesian product
$\mathcal{S}=\mathcal{S}^{(1)}\times...\times\mathcal{S}^{(D)}$ of
interpretable univariate codomains $\mathcal{S}^{(d)}\subset\mathbb{R}$. The
other subclass is more “intrinsically” multivariate insofar as univariate
analyses would not yield meaningful results. Consider for example two-
dimensional movement trajectories as in one of our motivating applications,
where the function measures Cartesian coordinates over time: for fixed
trajectories, rotation or translation of the essentially arbitrary coordinate
system would change the results of univariate analyses. For intrinsically
multivariate functional data a multivariate approach is the natural and
preferred mode of analysis, yielding interpretable results on the observation
level. Even for multimodal functional data, a joint analysis may generate
additional insight by incorporating the covariance structure between the
dimensions. This motivates the development of statistical methods for
multivariate functional data. We here propose multivariate functional additive
mixed models to model potentially sparsely observed functions with flexible
covariate effects and crossed or nested study designs.

Multivariate functional data have been the interest in different statistical
fields such as clustering (Jacques and Preda; 2014; Park and Ahn; 2017),
functional principal component analysis (Chiou et al.; 2014; Happ and Greven;
2018; Backenroth et al.; 2018; Li et al.; 2020), and registration (Carroll et
al.; 2021; Steyer et al.; 2021). There is also ample literature on
multivariate functional data regression such as graphical models (Zhu et al.;
2016), reduced rank regression (Liu et al.; 2020), or varying coefficient
models (Zhu et al.; 2012; Li et al.; 2017). Yet, so far, there are only few
approaches that can handle multilevel regression when the functional response
is multivariate. In particular, Goldsmith and Kitago (2016) propose a
hierarchical Bayesian multivariate functional regression model that can
include subject level and residual random effect functions to account for
correlation between and within functions. They work with bivariate functional
data observed on a regular and dense grid and assume _a priori_ independence
between the different dimensions of the subject-specific random effects. Thus,
they model the correlation between the dimensions only in the residual
function. As our approach explicitely models the dependencies between
dimensions for multiple functional random effects and also handles data
observed on sparse and irregular grids on more than two dimensions, the model
proposed by Goldsmith and Kitago (2016) can be seen as a special case of our
more general model class.

Alternatively, Zhu et al. (2017) use a two stage transformation with basis
functions for the multivariate functional mixed model. This allows the
estimation of scalar regression models for the resulting basis coefficients
that are argued to be approximately independent. The proposed model is part of
the so called functional mixed model (FMM) framework (Morris; 2017). While
FMMs use basis transformations of functional responses (observed on equal
grids) at the start of the analysis, we propose a multivariate model in the
functional additive mixed model (FAMM) framework, which uses basis
representations of all (effect) functions in the model (Scheipl et al.; 2015).
The differences between these two functional regression frameworks have been
extensively discussed before (Greven and Scheipl; 2017; Morris; 2017).

The main advantages of our multivariate regression model, also compared to
Goldsmith and Kitago (2016) and Zhu et al. (2017), are that it is readily
available for sparse and irregular functional data and that it allows to
include multiple nested or crossed random processes, both of which are
required in our data examples. Another important contribution is that our
approach directly models the multivariate covariance structure of all random
effects included in the model using multivariate functional principal
components (FPCs) and thus implicitly models the covariances between the
dimensions. This makes the model representation more parsimonious, avoids
assumptions difficult to verify, and allows further interpretation of the
random effect processes, such as their relative importance and their
dominating modes. As part of the FAMM framework, our model provides a vast
toolkit of modeling options for covariate and random effects, of estimation
and inference (Wood; 2017). The proposed multivariate functional additive
mixed model (multiFAMM) extends the FAMM framework combining ideas from
multilevel modeling (Cederbaum et al.; 2016) and multivariate functional data
(Happ and Greven; 2018) to account for sparse and irregular functional data
and different study designs.

We illustrate the multiFAMM on two motivating examples. The first
(intrinsically multivariate) data stem from a study on the effect of a
training programme for snooker players with a nested study design (shots
within sessions within players) (Enghofer; 2014). The movement trajectories of
a player’s elbow, hand, and shoulder during a snooker shot are recorded on
camera, yielding six-dimensional multivariate functional data (see Figure 1).
In the second data example, we analyse multimodal data from a speech
production study with a crossed study design (speakers crossed with words)
(Pouplier and Hoole; 2016) on so-called “assimilation” of consonants. The two
measured modes (acoustic and articulatory, see Figure 3) are expected to be
closely related but joint analyses have not yet incorporated the functional
nature of the data.

These two examples motivate the development of a regression model for sparse
and irregularly sampled multivariate functional data that can incorporate
crossed or nested functional random effects as required by the study design in
addition to flexible covariate effects. The proposed approach is implemented
in R (R Core Team; 2020) in package multifamm(Volkmann; 2021). The paper is
structured as follows: Section 2 specifies the multiFAMM and section 3 its
estimation process. Section 4 presents the application of the multiFAMM to the
data examples and section 5 shows the estimation performance of our proposed
approach in simulations. Section 6 closes with a discussion and outlook.

##  2 Multivariate Functional Additive Mixed Model

###  2.1 General Model

Let $\bm{y}^{*}_{i}(t)=(y_{i}^{*(1)}(t),...,y_{i}^{*(D)}(t))^{\top}$ be the
multivariate functional response of unit $i=1,...,N$ over $t\in\mathcal{I}$,
consisting of dimensions $d=1,...,D$. Without loss of generality, we assume a
common one-dimensional interval domain $\mathcal{I}=[0,1]$ for all dimensions,
and square-integrable $y_{i}^{*(d)}\in L^{2}(\mathcal{I})$. Define
$L^{2}_{D}(\mathcal{I}):=L^{2}(\mathcal{I})\times...\times L^{2}(\mathcal{I})$
so that $\bm{y}^{*}_{i}\in L^{2}_{D}(\mathcal{I})$. The underlying smooth
function $\bm{y}_{i}^{*}$, however, is only evaluated at (potentially sparse
or dimension-specific) points
$\bm{y}_{it}^{*}=(y_{it}^{*(1)},...,y_{it}^{*(D)})^{\top}$ and the evaluation
is subject to white noise, i.e.
$\bm{y}_{it}=\bm{y}_{it}^{*}+\bm{\epsilon}_{it}$. The residual term
$\bm{\epsilon}_{it}$ reflects additional uncorrelated white noise measurement
error, following a $D$-dimensional multivariate normal distribution
$\mathcal{N}_{D}$ with zero-mean and diagonal covariance matrix
$\tilde{\bm{\Sigma}}=\operatorname{diag}(\sigma_{1}^{2},...,\sigma_{D}^{2})$
with dimension specific variances $\sigma_{d}^{2}$. We construct a
multivariate functional mixed model as

$\displaystyle\begin{split}\bm{y}_{it}&=\bm{y}_{i}^{*}(t)+\bm{\epsilon}_{it}=\bm{\mu}(\bm{x}_{i},t)+\bm{U}(t)\bm{z}_{i}+\bm{\epsilon}_{it}\\\
&=\bm{\mu}(\bm{x}_{i},t)+\sum_{j=1}^{q}\bm{U_{j}}(t)\bm{z}_{ij}+\bm{E}_{i}(t)+\bm{\epsilon}_{it},\quad
t\in\mathcal{I},\end{split}$ (1)

where

$\displaystyle\bm{U}_{j}(t)$
$\displaystyle=(\bm{U}_{j1}(t),...,\bm{U}_{jV_{U_{j}}}(t));j=1,...,q,$
$\displaystyle\bm{U}_{jv}(t)$
$\displaystyle\stackrel{{\scriptstyle\text{ind.c.}}}{{\sim}}MGP\left(\bm{0},K_{U_{j}}\right);v=1,...,V_{U_{j}};\forall
j=1,...,q,$ $\displaystyle\bm{E}_{i}(t)$
$\displaystyle\stackrel{{\scriptstyle\text{ind.c.}}}{{\sim}}MGP\left(\bm{0},K_{E}\right);i=1,...,N,\text{
and }$ $\displaystyle\bm{\epsilon}_{it}$
$\displaystyle\stackrel{{\scriptstyle\text{i.i.d.}}}{{\sim}}\mathcal{N}_{D}\left(\bm{0},\tilde{\bm{\Sigma}}=\operatorname{diag}(\sigma_{1}^{2},...,\sigma_{D}^{2})\right);i=1,...,N;\quad
t\in\mathcal{I}.$

We assume an additive predictor
$\bm{\mu}(\bm{x}_{i},\cdot)=\sum_{l=1}^{p}\bm{f}_{l}(\bm{x}_{i},\cdot)$ of
fixed effects, which consists of partial predictors
$\bm{f}_{l}(\bm{x}_{i},\cdot)=(f_{l}^{(1)}(\bm{x}_{i},\cdot),...,f_{l}^{(D)}(\bm{x}_{i},\cdot))^{\top}\in
L^{2}_{D}(\mathcal{I}),\ l=1,...,p$, that are multivariate functions depending
on a subset of the vector of scalar covariates $\bm{x}_{i}$. This allows to
include linear or smooth covariate effects as well as interaction effects
between multiple covariates as in the univariate FAMM (Scheipl et al.; 2015).
Partial predictors may also depend on dimension specific subsets of
covariates.

For random effects $\bm{U}$, we focus on model scenarios with $q$ independent
multivariate functional random intercepts for crossed and/or nested designs.
For group level $v=1,\dots,V_{U_{j}}$ within grouping layer $j=1,\dots,q$,
these take the value $\bm{U}_{jv}\in L^{2}_{D}(\mathcal{I})$. For each layer,
the $\bm{U}_{j1},...,\bm{U}_{jV_{U_{j}}}$ present independent copies of a
multivariate smooth zero-mean Gaussian random process. Analogously to scalar
linear mixed models, the $\bm{U}_{jv}$ model correlations between different
response functions $\bm{y}_{i}^{*}$ within the same group as well as variation
across groups. By arranging them in a $(D\times V_{U_{j}})$ matrix
$\bm{U}_{j}(t)$ per $t$, the $j$th random intercept can be expressed in the
common mixed model notation in (1) using appropriate group indicators
$\bm{z}_{ij}=(z_{ij1},\dots,z_{ijV_{U_{j}}})^{\top}$ for the respective
design.

Although technically a curve-specific functional random intercept, we
distinguish the smooth residuals $\bm{E}_{i}\in L^{2}_{D}(\mathcal{I})$ in the
notation, as they model correlation within rather than between response
functions. We write $\bm{E}_{v}\in L^{2}_{D}(\mathcal{I}),v=1,...,V_{E}$ with
$V_{E}=N$. The $\bm{E}_{i}$ capture smooth deviations from the group-specific
mean $\bm{\mu}(\bm{x}_{i},\cdot)+\sum_{j=1}^{q}\bm{U}_{j}(\cdot)\bm{z}_{ij}$.

For a more compact representation, we can arrange all $\bm{U}_{j}(t)$ and
$\bm{E}_{i}(t)$ together in a $(D\times(\sum_{j=1}^{q}V_{U_{j}}+N))$ matrix
$\bm{U}(t)$ per $t$, and the group indicators for all layers in a
corresponding vector
$\bm{z}_{i}=(\bm{z}_{i1}^{\top},\dots,\bm{z}_{iq}^{\top},\bm{e}_{i}^{\top})^{\top}$
with $\bm{e}_{i}$ the $i$-th unit vector. The resulting model term
$\bm{U}(t)\bm{z}_{i}$ then comprises all smooth random functions, accounting
for all correlation between/within response functions $\bm{y}_{i}^{*}$ given
the covariates $\bm{x}_{i}$ as required by the respective experimental design.

$\bm{E}_{i}$ and $\bm{U}_{jv}$ are independent copies (ind. c.) of random
processes having multivariate $D\times D$ covariance kernels
$K_{E},K_{U_{j}}$, with univariate covariance surfaces
$K_{E}^{(d,e)}(t,t^{\prime})=\text{Cov}\left[E_{i}^{(d)}(t),E_{i}^{(e)}(t^{\prime})\right]$
and
$K_{U_{j}}^{(d,e)}(t,t^{\prime})=\text{Cov}\left[U_{jv}^{(d)}(t),U_{jv}^{(e)}(t^{\prime})\right]$
reflecting the covariance between the process dimensions $d$ and $e$ at $t$
and $t^{\prime}$. We call these auto-covariances for $d=e$ and cross-
covariance otherwise. The multivariate Gaussian processes are uniquely defined
by their multivariate mean function, here the null function $\bf{0}$, and the
multivariate covariance kernels $K_{g}$ and we write
$MGP\left(\bm{0},K_{g}\right),g\in\\{U_{1},\dots,U_{q},E\\}$. Note that
vectorizing the matrix $\bm{U}(t)$ allows to formulate the joint distribution
assumption $\operatorname{vec}(\bm{U}(t))\sim MGP\left(\bm{0},K_{U}\right)$
with $K_{U}(t,t^{\prime})$ having a block-diagonal structure repeating each
$K_{U_{j}}(t,t^{\prime})$ for $V_{U_{j}}$ times and $K_{E}(t,t^{\prime})$ for
$N$ times.

We assume that the different sources of variation
$\bm{U}_{j}(t),j=1,...,q,\bm{E}_{i}(t)$, and $\bm{\epsilon}_{it}$ are mutually
uncorrelated random processes to assure model identification. Assuming
smoothness of the covariance kernel $K_{E}$ further guarantees that the
residual process $\bm{E}_{i}(t)$ can be separated from the white noise
$\bm{\epsilon}_{it}$, removing the error variance from the diagonal of the
smooth covariance kernel (e.g., Yao et al.; 2005).

###  2.2 FPC Representation of the Random Effects

Model (1) specifies a univariate functional linear mixed model (FLMM) as given
in Cederbaum et al. (2016) for each dimension $d$. The main difference lies in
the multivariate random processes that introduce dependencies between the
dimensions. In order to avoid restrictive assumptions about the structure of
these multivariate covariance operators, which would typically be very
difficult to elicit _a priori_ or verify _ex post_ , we estimate them directly
from the data. The main difficulty then becomes computationally efficient
estimation, which is already costly in the univariate case. Especially for
higher dimensional multivariate functional data, accounting for the cross-
covariances can become a complex task, which we tackle with multivariate
functional principal component analysis (MFPCA).

Given the covariance operators (see Section 3), we represent the multivariate
random effects in model (1) using truncated multivariate Karhunen-Loève (KL)
expansions

$\displaystyle\bm{U}_{jv}(t)$
$\displaystyle\approx\sum_{m=1}^{M_{U_{j}}}\rho_{U_{j}vm}\bm{\psi}_{U_{j}m}(t),\;j=1,...,q;\,v=1,...,V_{U_{j}},$
(2) $\displaystyle\bm{E}_{v}(t)$
$\displaystyle\approx\sum_{m=1}^{M_{E}}\rho_{Evm}\bm{\psi}_{Em}(t),\;v=1,...,N,$

where the orthonormal multivariate eigenfunctions
$\bm{\psi}_{gm}=(\psi_{gm}^{(1)},...,\psi_{gm}^{(D)})^{\top}\in
L^{2}_{D}(\mathcal{I})$, $m=1,...,M_{g}$, $g\in\\{U_{1},...,U_{q},E\\}$ of the
corresponding covariance operators with truncation order $M_{g}$ are used as
basis functions and the random scores $\rho_{gvm}\sim N(0,\nu_{gm})$ are
independent and identically distributed (i.i.d.) with ordered eigenvalues
$\nu_{gm}$ of the corresponding covariance operator. Note that the assumption
of Gaussianity for the random processes can be relaxed. For non-Gaussian
random processes, the KL expansion still gives uncorrelated (but non-normal)
scores and estimation based on a penalized least squares (PLS) criterion (see
Section 3.2) remains reasonable.

Using KL expansions gives a parsimonious representation of the multivariate
random processes that is an optimal approximation with respect to the
integrated squared error (cf. Ramsay and Silverman; 2005), as well as
interpretable basis functions capturing the most prominent modes of variation
of the respective process. The distinct feature of this approach is that the
multivariate FPCs directly account for the dependency structure of each random
process across the dimensions. If, by contrast, e.g. splines were used in the
basis representation of the random effects, it would be necessary to
explicitly model the cross-covariances of each random process in the model
(cf. Li et al.; 2020). Multivariate eigenfunctions, however, are designed to
incorporate the dependency structure between dimensions and allow the
assumption of independent (univariate) basis coefficients $\rho_{gvm}$ via the
KL theorem (see e.g. Happ and Greven; 2018). This leads to a parsimonious
multivariate basis for the random effects, where a typically small vector of
scalar scores $\rho_{gvm}$ common to all dimensions represents nearly the
entire information about these $D$-dimensional processes.

##  3 Estimation

We use a two-step approach to estimate the multiFAMM and the respective
multivariate covariance operators. In a first step (section 3.1), the
D-dimensional eigenfunctions $\bm{\psi}_{gm}(t)$ with their corresponding
eigenvalues $\nu_{gm}$ are estimated from their univariate counterparts
following Cederbaum et al. (2018) and Happ and Greven (2018). These estimates
are then plugged into (2) and we represent the multiFAMM as part of the
general FAMM framework (section 3.2) by suitable re-arrangement. We can view
the estimated $\bm{\psi}_{gm}(t)$ simply as an empirically derived basis that
parsimoniously represents the patterns in the observed data. While their
estimation adds uncertainty, we are not interested in inferential statements
for the variance modes and our simulations (see Section 5) suggest that the
estimated eigenfunctions are reasonable approximations that work well as a
basis.

###  3.1 Step 1: Estimation of the Eigenfunction Basis

#### Step 1 (i): Univariate Mean Estimation

In a first step, we obtain preliminary estimates of the dimension specific
means $\mu^{(d)}(\bm{x}_{i},t)=\sum_{l=1}^{p}f_{l}^{(d)}(\bm{x}_{il},t)$ using
univariate FAMMs. We model

$\displaystyle y_{it}^{(d)}$
$\displaystyle=\mu^{(d)}(\bm{x}_{i},t)+\epsilon^{(d)}_{it};\;d=1,\dots,D$ (3)

independently for all $d$ with i.i.d. Gaussian random variables
$\epsilon^{(d)}_{it}$. The estimation of $\mu^{(d)}(\bm{x}_{i},t)$ proceeds
analogously to the estimation of the multiFAMM described in section 3.2. It is
based on the evaluation points of the $y_{i}^{*(d)}(t)$, whose locations on
the interval $\mathcal{I}$ can vary across dimensions. Model (3) thus
accommodates sparse and irregular multivariate functional data and implies a
working independence assumption across scalar observations within and across
functions.

#### Step 1 (ii): Univariate Covariance Estimation

This preliminary mean function is used to centre the data
$\tilde{y}_{it}^{(d)}=y_{it}^{(d)}-\hat{\mu}^{(d)}(\bm{x}_{i},t)$ in order to
obtain noisy evaluations of the detrended functions
$\tilde{y}_{i}^{*(d)}(t)=y_{i}^{*(d)}(t)-{\mu}^{(d)}(\bm{x}_{i},t)$ for
covariance estimation. Cederbaum et al. (2016) already find that for this
purpose, the working independence assumption within functions across
evaluation points in (3) gives reasonable results. The expectation of the
crossproducts of the centred functions then coincides with the auto-
covariance, i.e.
$\mathbb{E}\left(\tilde{y}_{it}^{(d)}\tilde{y}_{i^{\prime}t^{\prime}}^{(d)}\right)\approx\text{Cov}\left[y_{it}^{(d)},y_{i^{\prime}t^{\prime}}^{(d)}\right]$.
For the independent random components specified in model (1), this overall
covariance decomposes additively into contributions from each random process
as

$\displaystyle\mathbb{E}\left(\tilde{y}_{it}^{(d)}\tilde{y}_{i^{\prime}t^{\prime}}^{(d)}\right)\approx\sum_{j=1}^{q}K_{U_{j}}^{(d,d)}(t,t^{\prime})\delta_{v_{j}v_{j}^{\prime}}+\big{(}K_{E}^{(d,d)}(t,t^{\prime})+\sigma^{2}_{d}\delta_{tt^{\prime}}\big{)}\delta_{ii^{\prime}},$
(4)

using indicators $\delta_{xx^{\prime}}$ that equal one for $x=x^{\prime}$ and
zero otherwise. The indicator $\delta_{v_{j}v_{j}^{\prime}}$ thus identifies
if the curves in the crossproduct belong to the same group $v_{j}$ of the
$j$th layer. Using $t$, $t^{\prime}$, and the indicators
$\delta_{v_{j}v_{j}^{\prime}},\delta_{tt^{\prime}},\delta_{ii^{\prime}}$ as
covariates and the crossproducts of the centred data as responses, we can
estimate the auto-covariances $K_{U_{1}}^{(d,d)},...,K_{U_{q}}^{(d,d)},$ and
$K_{E}^{(d,d)}$ of the random processes using symmetric additive covariance
smoothing (Cederbaum et al.; 2018). This extends the univariate approach
proposed by Cederbaum et al. (2016). In particular, we also allow a nested
random effects structure as required for the snooker training application in
section 4.1 by specifying the indicator of the nested effect as the product of
subject and session indicators. Note that estimating (4) also yields estimates
of the dimension specific error variances $\sigma_{d}^{2}$ as a byproduct.

#### Step 1 (iii): Univariate Eigenfunction Estimation

Based on the covariance kernel estimates, we apply separate univariate
functional principal component analyses (FPCAs) for each random process by
conducting an eigendecomposition of the respective linear integral operator.
Practically, each estimated process- and dimension-specific auto-covariance is
re-evaluated on a dense grid so that a univariate FPCA can be conducted.
Alternatively, Reiss and Xu (2020) provide an explicit spline representation
of the estimated eigenfunctions. Eigenfunctions with non-positive eigenvalues
are removed to ensure positive definiteness, and further regularization by
truncation based on the proportion of variance explained is possible (see e.g.
Di et al.; 2009; Peng and Paul; 2009; Cederbaum et al.; 2016). However, we
suggest to keep all univariate FPCs with positive eigenvalues for the
computation of the MFPCA in order to preserve all important modes of variation
and cross-correlation in the data.

#### Step 1 (iv): Multivariate Eigenfunction Estimation

The estimated univariate eigenfunctions and scores are then used to conduct an
MFPCA for each of the $g$ multivariate random processes separately. The MFPCA
exploits correlations between univariate FPC scores across dimensions to
reduce the number of basis functions needed to sufficiently represent the
random processes. We base the MFPCA on the following definition of a
(weighted) scalar product

$\displaystyle\langle\langle\bm{f},\bm{g}\rangle\rangle:=\sum_{d=1}^{D}w_{d}\int_{\mathcal{I}}f^{(d)}(t)g^{(d)}(t)dt,\quad\bm{f},\bm{g}\in
L^{2}_{D}(\mathcal{I}),$ (5)

for the response space with positive weights $w_{d},d=1,...,D$ and the induced
norm denoted by $|||\cdot|||$. The corresponding covariance operators
$\Gamma_{g}:L^{2}_{D}(\mathcal{I})\rightarrow L^{2}_{D}(\mathcal{I})$ of the
multivariate random processes $\bm{U}_{jv}$ and $\bm{E}_{v}$ are then given by
$(\Gamma_{g}\bm{f})(t)=\langle\langle\bm{f},K_{g}(t,\cdot)\rangle\rangle$,
$g\in\\{U_{1},...,U_{q},E\\}$. The standard choice of weights in our
applications is $w_{1}=...=w_{D}=1$ (unweighted scalar product) but other
choices are possible. Consider for example a scenario where dimensions are
observed with different amounts of measurement error. If variation in
dimensions with a large proportion of measurement error is to be downweighted,
we propose to use $w_{d}=\frac{1}{\hat{\sigma}^{2}_{d}}$ with the dimension
specific measurement error variance estimates $\hat{\sigma}^{2}_{d}$ obtained
from (4).

Happ and Greven (2018) show that estimates of the multivariate eigenvalues
$\nu_{gm}$ of $\Gamma_{g}$ can be obtained from an eigenanalysis of a
covariance matrix of the univariate random scores. The corresponding
multivariate eigenfunctions $\bm{\psi}_{gm}$ can be obtained as linear
combinations of the univariate eigenfunctions with the weights given by the
resulting eigenvectors. The estimates $\hat{\bm{\psi}}_{gm}$ are then
substituted for the basis functions of the truncated multivariate KL
expansions of the random effects $\bm{U}_{jv}$ and $\bm{E}_{v}$ in (2). Note
that for each random process $g$, the maximum number of FPCs is given by the
total number of univariate eigenfunctions included in the estimation process
of the MFPCA of $g$. To achieve further regularization and analogously to
Cederbaum et al. (2016), we propose to choose truncation orders $M_{g}$ for
each KL expansion of the multivariate random processes using a prespecified
proportion of explained variation.

#### Step 1 (v): Multivariate Truncation Order

We offer two different approaches for the choice of truncation orders $M_{g}$
based on different variance decompositions (derivation in Appendix A):

$\displaystyle\mathbb{E}\left(|||\bm{y}_{i}-\bm{\mu}(\bm{x}_{i})|||^{2}\right)=\sum_{d=1}^{D}w_{d}\int_{\mathcal{I}}\mathrm{Var}\big{(}y_{i}^{(d)}(t)\big{)}dt=\sum_{g}\sum_{m=1}^{\infty}\nu_{gm}+\sum_{d=1}^{D}w_{d}\sigma_{d}^{2}|\mathcal{I}|,$
(6)
$\displaystyle\text{and}\quad\int_{\mathcal{I}}\mathrm{Var}\big{(}y_{i}^{(d)}(t)\big{)}dt=\sum_{g}\sum_{m=1}^{\infty}\nu_{gm}||\psi_{gm}^{(d)}||^{2}+\sigma_{d}^{2}|\mathcal{I}|$
(7)

with $|\mathcal{I}|$ the length of the interval $\mathcal{I}$ (here equal to
one) and $||\cdot||$ the $L^{2}$ norm. Multivariate variance decomposition (6)
uses the (weighted) sum of total variation in the data across dimensions. We
select the FPCs with highest associated eigenvalues $\nu_{gm}$ over all random
processes $g$ until their sum reaches a prespecified proportion (e.g. 0.95) of
the total variation, thus approximating the infinite sums in (6) with $M_{g}$
summands. For the approach based on the univariate variance (7), we require
$M_{g}$ to be the smallest truncation order for which at least a prespecified
proportion of variance is explained on every dimension $d$. This second choice
of $M_{g}$ might be preferable in situations where the variation is
considerably different (in amount or structure) across dimensions, whereas the
first approach gives a more parsimonious representation of the random effects.
Note that both approaches can lead to a simplification of the multiFAMM if
$M_{g}=0$ is chosen for some $g$. The simulation results of section 5 suggest
that increasing the number of FPCs improves model accuracy which is why
sensitivity analyses with regard to the truncation order are recommended.

###  3.2 Step 2: Estimation of the multiFAMM

In the following, we discuss estimating the multiFAMM given the estimated
multivariate FPCs. We base the proposed model on the general FAMM framework of
Scheipl et al. (2015), which models functional responses using basis
representations. To make the extension of the FAMM framework to multivariate
functional data more apparent, the multivariate response vectors and the
respective model matrices are stacked over dimensions, so that every block has
the structure of a univariate FAMM over all observations $i$. This gives
concatenated basis functions with discontinuities between the dimensions. The
fixed effects are modelled analogously to the univariate case by interacting
all covariate effects with a dimension indicator. The random effects are based
on the parsimonious, concatenated multivariate FPC basis.

#### Matrix Representation

For notational simplicity we assume that the functions are evaluated on a
fixed grid of time points $\bm{t}=(\bm{t}^{(1)},...,\bm{t}^{(D)})^{\top}$ with
$\bm{t}^{(d)}=(\bm{t}_{1}^{(d)},...,\bm{t}_{N}^{(d)})$ and identical
$\bm{t}_{i}^{(d)}\equiv(t_{1},...,t_{T})$ over all $N$ individuals and $D$
dimensions. However, our framework allows for sparse functional data using
different grids per dimension and per observed function as in the two
applications (Section 4). Correspondingly,
$\bm{y}=(\bm{y}^{(1)\top},...,\bm{y}^{(D)\top})^{\top}$ is the $DNT$-vector of
stacked evaluation points with
$\bm{y}^{(d)}=(\bm{y}_{1}^{(d)\top},...,\bm{y}_{N}^{(d)\top})^{\top}$ and
$\bm{y}_{i}^{(d)}=(y_{i1}^{(d)},...,y_{iT}^{(d)})^{\top}$. Model (1) on this
grid can be written as

$\displaystyle\bm{y}=\bm{\Phi\theta}+\bm{\Psi\rho}+\bm{\epsilon}$ (8)

with $\bm{\Phi},\bm{\Psi}$ the model matrices for the fixed and random
effects, respectively, $\bm{\theta},\bm{\rho}$ the vectors of coefficients and
random effect scores to be estimated, and
$\bm{\epsilon}=(\bm{\epsilon}^{(1)\top},...,\bm{\epsilon}^{(D)\top})^{\top}$,
$\bm{\epsilon}^{(d)}=(\epsilon_{11}^{(d)},...,\epsilon_{1T}^{(d)},...,\epsilon_{NT}^{(d)})^{\top}$
the vector of residuals. We have $\bm{\epsilon}\sim N(\bm{0},\bm{\Sigma})$
with
$\bm{\Sigma}=\operatorname{diag}(\sigma_{1}^{2},...,\sigma_{D}^{2})\otimes\bm{I}_{NT}$,
the Kronecker product denoted by $\otimes$, and the $(NT\times NT)$ identity
matrix $\bm{I}_{NT}$.

We estimate $\bm{\theta}$ and $\bm{\rho}$ by minimizing the PLS criterion

$\displaystyle(\bm{y}-\bm{\Phi\theta}-\bm{\Psi\rho})\bm{\Sigma}^{-1}(\bm{y}-\bm{\Phi\theta}-\bm{\Psi\rho})^{\top}+\sum_{l=1}^{p}\bm{\theta}_{l}^{\top}\bm{P}_{l}(\bm{\lambda}_{xl},\bm{\lambda}_{tl})\bm{\theta}_{l}+\sum_{g}\lambda_{g}\bm{\rho}_{g}^{\top}\bm{P}_{g}\bm{\rho}_{g}$
(9)

using appropriate penalty matrices
$\bm{P}_{l}(\bm{\lambda}_{xl},\bm{\lambda}_{tl})$ and $\bm{P}_{g}$ for the
fixed effects and random effects, respectively, and smoothing parameters
$\bm{\lambda}_{xl}=\left(\lambda_{xl}^{(1)},...,\lambda_{xl}^{(D)}\right),\bm{\lambda}_{tl}=\left(\lambda_{tl}^{(1)},...,\lambda_{tl}^{(D)}\right)$,
and $\lambda_{g}$. The model and penalty matrices as well as the parameter
vectors of (8) and (9) are discussed in detail below.

#### Modeling of Fixed Effects

The block diagonal matrix
$\bm{\Phi}=\operatorname{diag}\left(\bm{\Phi}^{(1)},...,\bm{\Phi}^{(D)}\right)$
models the fixed effects separately on each dimension as in a FAMM (Scheipl et
al.; 2015). The $(DNT\times b)$ matrix $\bm{\Phi}$ consists of the design
matrices
$\bm{\Phi}^{(d)}=(\bm{\Phi}_{1}^{(d)}\,|\,...\,|\,\bm{\Phi}_{p}^{(d)})$ that
are constructed for the partial predictors
$f_{l}^{(d)}(\bm{x},\bm{t}^{(d)}),l=1,...,p$, which correspond to the
$NT$-vectors of evaluations of the effect functions $f_{l}^{(d)}$. The vectors
of scalar covariates $\bm{x}_{i}$ are repeated $T$ times to form the matrix of
covariate information
$\bm{x}=(\bm{x}_{1},...,\bm{x}_{1},...,\bm{x}_{N})^{\top}$. We use the basis
representations

$\displaystyle
f_{l}^{(d)}(\bm{x},\bm{t}^{(d)})\approx\bm{\Phi}_{l}^{(d)}\bm{\theta}_{l}^{(d)}=(\bm{\Phi}_{xl}^{(d)}\odot\bm{\Phi}_{tl}^{(d)})\bm{\theta}_{l}^{(d)},$

where $\bm{A}\odot\bm{B}$ denotes the row tensor product
$(\bm{A}\otimes\bm{1}_{b}^{\top})\cdot(\bm{1}_{a}^{\top}\otimes\bm{B})$ of the
$(h\times a)$ matrix $\bm{A}$ and the $(h\times b)$ matrix $\bm{B}$ with
element-wise multiplication $\cdot$ and $\bm{1}_{c}$ the $c$-vector of ones.
This modeling approach combines the $(NT\times b_{xl}^{(d)})$ basis matrix
$\bm{\Phi}_{xl}^{(d)}$ with the $(NT\times b_{tl}^{(d)})$ basis matrix
$\bm{\Phi}_{tl}^{(d)}$. These matrices contain the evaluations of suitable
marginal bases in $\bm{x}$ and $\bm{t}^{(d)}$, respectively. For a linear
effect, for example, the basis matrix $\bm{\Phi}_{xl}^{(d)}$ is specified as
the familiar linear model design matrix $\bm{x}$ for the linear effect
$f_{l}^{(d)}(\bm{x},\bm{t}^{(d)})=\bm{x}\bm{\beta}_{l}^{(d)}(\bm{t}^{(d)})$
with coefficient function $\bm{\beta}_{l}^{(d)}(\bm{t}^{(d)})$. For a
nonlinear effect
$f_{l}^{(d)}(\bm{x},\bm{t}^{(d)})=g_{l}^{(d)}(\bm{x},\bm{t}^{(d)})$, the basis
matrix $\bm{\Phi}_{xl}^{(d)}$ contains an (e.g. B-spline) basis representation
analogously to a scalar additive model. For the functional intercept,
$\bm{\Phi}_{xl}^{(d)}$ is a vector of ones, and we generally use a spline
basis for $\bm{\Phi}_{tl}^{(d)}$. For a complete list of possible effect
specifications with examples, we refer to Scheipl et al. (2015). The tensor
product basis is weighted by the $b_{xl}^{(d)}b_{tl}^{(d)}$ unknown basis
coefficients in $\bm{\theta}_{l}^{(d)}$. Stacking the vectors
$\bm{\theta}_{l}^{(d)}$ gives
$\bm{\theta}^{(d)}=(\bm{\theta}_{1}^{(d)\top},...,\bm{\theta}_{p}^{(d)\top})^{\top}$
and finally the $b$-vector
$\bm{\theta}=(\bm{\theta}^{(1)\top},...,\bm{\theta}^{(D)\top})^{\top}$ with
$b=\sum_{d}\sum_{l}b_{xl}^{(d)}b_{tl}^{(d)}$.

Choosing the number of basis functions is a well known challenge in the
estimation of nonlinear or functional effects. We introduce regularization by
a corresponding quadratic penalty term in (9). Let $\bm{\theta}_{l}$ contain
the coefficients corresponding to the partial predictor $l$ and order it by
dimensions. The penalty $\bm{P}_{l}(\bm{\lambda}_{xl},\bm{\lambda}_{tl})$ is
then constructed from the penalty on the marginal basis for the covariate
effect, $\bm{P}_{xl}^{(d)}$, and the penalty on the marginal basis over the
functional index, $\bm{P}_{tl}^{(d)}$. Specifically,
$\bm{P}_{l}(\bm{\lambda}_{xl},\bm{\lambda}_{tl})$ is a block-diagonal matrix
with blocks for each $d$ corresponding to the Kronecker sums of the marginal
penalty matrices
$\lambda_{xl}^{(d)}\bm{P}_{xl}^{(d)}\otimes\bm{I}_{b_{tl}^{(d)}}+\lambda_{tl}^{(d)}\bm{I}_{b_{xl}^{(d)}}\otimes\bm{P}_{tl}^{(d)}$
(Wood; 2017). A standard choice for these marginal penalty matrices given a
B-splines basis representation are second or third order difference penalties,
thus approximately penalizing squared second or third derivatives of the
respective functions (Eilers and Marx; 1996). For unpenalized effects such as
a linear effect of a scalar covariate, the corresponding $\bm{P}_{xl}^{(d)}$
is simply a matrix of zeroes.

#### Modeling of Random Effects

We represent the $DNT$-vectors
$\bm{U}_{j}(\bm{t})=(\bm{U}_{j}(\bm{t}^{(1)})^{\top},...,\bm{U}_{j}(\bm{t}^{(D)})^{\top})^{\top}$,
$\bm{E}(\bm{t})=(\bm{E}(\bm{t}^{(1)})^{\top},...,\bm{E}(\bm{t}^{(D)})^{\top})^{\top}$
with $\bm{U}_{j}(\bm{t}^{(d)})$, $\bm{E}(\bm{t}^{(d)})$ containing the
evaluations of the univariate random effects for the corresponding groups and
time points using the basis approximations

$\displaystyle\bm{U}_{j}(\bm{t})\approx\bm{\Psi}_{U_{j}}\bm{\rho}_{U_{j}}=(\bm{\delta}_{U_{j}}\odot\tilde{\bm{\Psi}}_{U_{j}})\bm{\rho}_{U_{j}},\quad\bm{E}(\bm{t})\approx\bm{\Psi}_{E}\bm{\rho}_{E}=(\bm{\delta}_{E}\odot\tilde{\bm{\Psi}}_{E})\bm{\rho}_{E}.$

The $v$th column in the $(DNT\times V_{g}),g\in\\{U_{1},...,U_{q},E\\}$
indicator matrix $\bm{\delta}_{g}$ indicates whether a given row is from the
$v$th group of the corresponding grouping layer. Thus, the rows of the
indicator matrix $\bm{\delta}_{g}$ contain repetitions of the group indicators
$\bm{z}_{ij}^{\top}$ and $\bm{e}_{i}^{\top}$ in model (1). For the smooth
residual, $\bm{\delta}_{E}$ simplifies to
$\bm{1}_{D}\otimes(\bm{I}_{N}\otimes\bm{1}_{T})$. The $(DNT\times M_{g})$
matrix
$\tilde{\bm{\Psi}}_{g}=(\tilde{\bm{\Psi}}_{g}^{(1)\top}|...|\tilde{\bm{\Psi}}_{g}^{(D)\top})^{\top}$
comprises the evaluations of the $M_{g}$ multivariate eigenfunctions
$\psi_{gm}^{(d)}(t)$ on dimensions $d=1,...,D$ for the $NT$ time points
contained in the $(NT\times M_{g})$ matrix $\tilde{\bm{\Psi}}_{g}^{(d)}$. The
$M_{g}V_{g}$ vector
$\bm{\rho}_{g}=(\bm{\rho}_{g1}^{\top},...,\bm{\rho}_{gV_{g}}^{\top})^{\top}$
with $\bm{\rho}_{gv}=(\rho_{gv1},...,\rho_{gvM_{g}})^{\top}$ stacks all the
unknown random scores for the functional random effect $g$. The
$(DNT\times\sum_{g}M_{g}V_{g})$ model matrix
$\bm{\Psi}=(\bm{\Psi}_{U_{1}}|...|\bm{\Psi}_{U_{q}}|\bm{\Psi}_{E})$ then
combines all random effect design matrices. Stacking the vectors of random
scores in a $\sum_{g}M_{g}V_{g}$ vector
$\bm{\rho}=(\bm{\rho}_{U_{1}}^{\top},...,\bm{\rho}_{U_{q}}^{\top},\bm{\rho}_{E}^{\top})^{\top}$
lets us represent all functional random intercepts in the model via
$\bm{\Psi}\bm{\rho}$.

For a given functional random effect, the penalty takes the form
$\bm{\rho}_{g}^{\top}\bm{P}_{g}\bm{\rho}_{g}=\bm{\rho}_{g}^{\top}(\bm{I}_{V_{g}}\otimes\tilde{\bm{P}}_{g})\bm{\rho}_{g}$,
where $\bm{I}_{V_{g}}$ corresponds to the assumed independence between the
$V_{g}$ different groups. The diagonal matrix
$\tilde{\bm{P}}_{g}=\operatorname{diag}(\nu_{g1},...,\nu_{gM_{g}})^{-1}$
contains the (estimated) eigenvalues $\nu_{gm}$ of the associated multivariate
FPCs. This quadratic penalty is mathematically equivalent to a normal
distribution assumption on the scores $\bm{\rho}_{gv}$ with mean zero and
covariance matrix $\tilde{\bm{P}}_{g}^{-1}$, as implied by the KL theorem for
Gaussian random processes. Note that the smoothing parameter $\lambda_{g}$
allows for additional scaling of the covariance of the corresponding random
process.

#### Estimation

We estimate the unknown smoothing parameters in
$\bm{\lambda}_{xl},\bm{\lambda}_{tl}$, and $\lambda_{g}$ using fast restricted
maximum likelihood (REML)-estimation (Wood; 2017). The standard
identifiability constraints of FAMMs are used (Scheipl et al.; 2015). In
particular, in addition to the constraints for the fixed effects, the
multivariate random intercepts are subject to a sum-to-zero constraint over
all evaluation points as given by e.g. Goldsmith et al. (2016).

We propose a weighted regression approach to handle the heteroscedasticity
assumption contained in $\bm{\Sigma}$. We weigh each observation
proportionally to the inverse of the estimated univariate measurement error
variances $\hat{\sigma}_{d}^{2}$ from the estimation of the univariate
covariances (4). Alternatively, updated measurement error variances can be
obtained from fitting separate univariate FAMMs on the dimensions using the
univariate components of the multivariate FPCs basis. In practice, we found
that the less computationally intensive former option gives reasonable
results.

As our proposed model is part of the FAMM framework, inference for the
multiFAMM is readily available based on inference for scalar additive mixed
models (Wood; 2017). Note, however, that all inferential statements do not
incorporate uncertainty due to the estimated multivariate eigenfunction bases,
nor in the chosen smoothing parameters. The estimation process readily
provides i.a. standard errors for the construction of point-wise univariate
confidence bands (CBs).

###  3.3 Implementation

We provide an implementation of the estimation of the proposed multiFAMM in
the multifamm R-package (Volkmann; 2021). It is possible to include up to two
functional random intercepts in $\bm{U}(t)$, which can have a nested or
crossed structure, in addition to the curve-specific random intercept
$\bm{E}_{i}(t)$. While including e.g. functional covariates is conceptually
straightforward (see Scheipl et al.; 2015), our implementation is restricted
to scalar covariates and interactions thereof. We provide different
alternatives for specifying the multivariate scalar product, the multivariate
cut-off criterion, and the covariance matrix of the white noise error term.
Note that the estimated univariate error variances have been proposed as
weights for two separate and independent modeling decisions: as weights in the
scalar product of the MFPCA and as regression weights under heteroscedasticity
across dimensions.

##  4 Applications

We illustrate the proposed multiFAMM for two different data applications
corresponding to intrinsically multivariate and multimodal fuctional data. The
presentation focuses on the first application with a detailed description of
the multimodal data application in Appendix C. We provide the data and the
code to produce all analyses in the github repository multifammPaper.

###  4.1 Snooker Training Data

#### Data Set and Preprocessing

In a study by Enghofer (2014), 25 recreational snooker players split into two
groups, one of which had instructions to follow a self-administered training
schedule over the next six weeks consisting of exercises aimed at improving
snooker specific muscular coordination. The second was a control group. Before
and after the training period, both groups were recorded on high speed digital
camera under similar conditions to investigate the effects of the training on
their snooker shot of maximal force. In each of the two recording sessions,
six successful shots per participant were videotaped. The recordings were then
used to manually locate points of interest (a participant’s shoulder, elbow,
and hand) and track them on a two-dimensional grid over the course of the
video. This yields a six dimensional functional observation per snooker shot
$\bm{y}^{*}=(y^{*(\text{elbow.x})},...,y^{*(\text{shoulder.y})}):\mathcal{I}=[0,1]\to\mathbb{R}^{6}$,
i.e. a two-dimensional movement trajectory for each point of interest (see
Figure 1).

Figure 1: Left: Screenshot of software for tracking (blue lines) the points of
interest (circles). Right: Two-dimensional trajectories of the snooker
training data set (grey curves, right). For both groups of skilled and
unskilled participants, three randomly selected observations are highlighted.
Every colour corresponds to one multivariate observation, i.e. one observation
consists of three trajectories: elbow (top), shoulder (right), hand (bottom).
The start of the exemplary trajectories are marked with a black asterisk with
the hand trajectory centred at the origin.

In their starting position (hand centred at the origin), the snooker players
are positioned centrally in front of the snooker table aiming at the cue ball.
From their starting position, the players draw back the cue, then accelerate
it forwards and hit the cue ball shortly after their hands enter the positive
range of the horizontal $x$ axis. After the impulse onto the cue ball, the
hand movement continues until it is stopped at a player’s chest. Enghofer
(2014) identify two underlying techniques that a player can apply: dynamic and
fixed elbow. With a dynamic elbow, the cue can be moved in an almost straight
line (piston stroke) whereas additionally fixing the elbow results in a
pendular motion (pendulum stroke). In both cases, the shoulder serves as a
fixed point and should be positioned close to the snooker table.

We adjust the data for differences in body height and relative speed (Steyer
et al.; 2021) and apply a coarsening method to reduce the number of redundant
data points, thereby lowering computational demands of the analysis. Appendix
B provides a detailed description of the data preprocessing. As some
recordings and evaluations of bivariate trajectories are missing, the final
data set contains 295 functional observations with a total of 56,910
evaluation points. These multivariate functional data are irregular and
sparse, with a median of 30 evaluation points per functional observation
(minimum 8, maximum 80) for each of the six dimensions.

#### Model Specification

We estimate the following model

$\displaystyle\bm{y}_{ijht}=\bm{\mu}(\bm{x}_{ij},t)+\bm{B}_{i}(t)+\bm{C}_{ij}(t)+\bm{E}_{ijh}(t)+\bm{\epsilon}_{ijht},$
(10)

with $i=1,...,25$ the index for the snooker player, $j=1,2$ the index for the
session, $h=1,...,H_{ij}$ the index for the typically six snooker shot
repetitions in a session, and $t\in[0,1]$ relative time. Correspondingly,
$\bm{B}_{i}(t)$ is a subject-specific random intercept, $\bm{C}_{ij}(t)$ is a
nested subject-and-session-specific random intercept, and $\bm{E}_{ijh}(t)$ is
the shot-specific random intercept (smooth residual). The nested random effect
$\bm{C}_{ij}(t)$ is supposed to capture the variation within players between
sessions (e.g. differences due to players having a good or bad day). Different
positioning of participants with respect to the recording equipment or the
snooker table as well as shot to shot variation are captured by the smooth
residual $\bm{E}_{ijh}(t)$. The white noise measurement error
$\bm{\epsilon}_{ijht}$ is assumed to follow a zero-mean multivariate normal
distribution with covariance matrix $\sigma^{2}\bm{I}_{6}$, as all six
dimensions are measured with the same set-up. The additive predictor is
defined as

$\displaystyle\bm{\mu}(\bm{x}_{ij},t)$
$\displaystyle=\bm{f}_{0}(t)+\texttt{skill}_{i}\cdot\bm{f}_{1}(t)+\texttt{group}_{i}\cdot\bm{f}_{2}(t)+\texttt{session}_{j}\cdot\bm{f}_{3}(t)$
$\displaystyle\quad+\texttt{group}_{i}\cdot\texttt{session}_{j}\cdot\bm{f}_{4}(t).$

The dummy covariates $\texttt{skill}_{i}$ and $\texttt{group}_{i}$ indicate
whether player $i$ is an advanced snooker player and belongs to the treatment
group (i.e. receives the training programme), respectively. Note that the
snooker players self-select into training and control group to improve
compliance with the training programme, which is why we include a group effect
in the model. The dummy covariate $\texttt{session}_{j}$ indicates whether the
shot $j$ is recorded after the training period. The effect function
$\bm{f}_{4}(t)$ can thus be interpreted as the treatment effect of the
training programme.

Cubic P-splines with first order difference penalty, penalizing deviations
from constant functions over time, with 8 basis functions are used for all
effect functions in the preliminary mean estimation as well as in the final
multiFAMM. For the estimation of the auto-covariances of the random processes,
we use cubic P-splines with first order difference penalty on five marginal
basis functions. We use an unweighted scalar product (5) for the MFPCA to give
equal weight to all spatial dimensions, as we can assume that the measurement
error mechanism is similar across dimensions. Additionally, we find that hand,
elbow, and shoulder contribute roughly the same amount of variation to the
data, cf. Table 1 in Appendix B.3, where we also discuss potential weighting
schemes for the MFPCA. The multivariate truncation order is chosen such that
95% of the (unweighted) sum of variation (6) is explained.

#### Results

The MFPCA gives sets of five (for $\bm{C}$ and $\bm{E}$) and six (for
$\bm{B}$) multivariate FPCs that explain 95% of the total variation. The
estimated eigenvalues allow to quantify their relative importance.
Approximately 41% of the total variation (conditional on covariates) can be
attributed to the nested subject-and-session-specific random intercept
$\bm{C}_{ij}(t)$, 33% to the subject-specific random intercept
$\bm{B}_{i}(t)$, 14% to the shot-specific $\bm{E}_{ijh}(t)$, and 7% to white
noise. This suggests that day to day variation within a snooker player is
larger than the variation between snooker players. Note that these proportions
are based on estimation step 1 (see Section 3.1).

Figure 2: Left: Dominant mode ($\bm{\psi}_{C1}$) of the subject-and-session-
specific random effect, explaining $27.7\%$ of total variation and shown as
mean trajectory (black solid) plus (red $+$) or minus (blue $-$)
$2\sqrt{\nu_{C1}}$ times the first FPC. An asterisk marks the start of a
trajectory. Right: Estimated covariate effect functions for skill. The central
plot shows the effect of the coefficient function (red) on the two-dimensional
trajectories for the reference group (blue). The marginal plots show the
estimated univariate effect functions (black) with pointwise 95% CBs (dotted)
and the baseline (blue dashed).

The left plot of Figure 2 displays the first FPC for $\bm{C}$, which explains
about $28\%$ of total variation. A suitable multiple of the FPCs is added
($+$) to and subtracted ($-$) from the overall mean function (black solid
line, all covariate values set to $0.5$). We find that the dominant mode of
the random subject-and-session specific effect influences the relative
positioning of a player’s elbow, shoulder, and hand, thus suggesting a strong
dependence between the dimensions. Enghofer (2014) argue from a theoretical
viewpoint that the ideal starting position should place elbow and hand in a
line perpendicular to the plane of the snooker table (corresponding to the x
axis). The most prominent mode of variation captures deviations from this
ideal starting position found in the overall mean. The next most important FPC
$\bm{\psi}_{B1}$ of the subject-specific random effect, which explains about
$15\%$ of total variation, represents a subject’s tendency towards the piston
or pendulum stroke (see Appendix Figure 8). This additional insight into the
underlying structure of the variance components might be helpful for e.g.
developing personalized training programmes.

The central plot on the right of Figure 2 compares the estimated mean movement
trajectory for advanced snooker players (red) to that in the reference group
(blue). It suggests that more experienced players tend towards the dynamic
elbow technique, generating a hand trajectory resembling a straight line
(piston stroke). Uncertainties in the trajectory could be represented by
pointwise ellipses, but inference is more straightforward to obtain from the
univariate effect functions. The marginal plots display the estimated
univariate effects with pointwise 95% confidence intervals. Even though we
find only little statistical evidence for increased movement of the elbow
(horizontal-left and vertical-top marginal panels), the hand and shoulder
movements (horizontal centre and right, vertical centre and bottom) strongly
suggest that the skill level indeed influences the mean movement trajectory of
a snooker player. Further results indicate that the mean hand trajectories
might slightly differ between treatment and control group at baseline as well
as between sessions ($\bm{f}_{2}(t)$ and $\bm{f}_{3}(t)$, see Appendix Figure
12). The estimated treatment effect $\bm{f}_{4}(t)$ (Appendix Figure 11),
however, suggests that the training programme did not change the participants’
mean movement trajectories substantially. Appendix B.3 contains a detailed
discussion of all model terms as well as some model diagnostics and
sensitivity analyses.

###  4.2 Consonant Assimilation Data

#### Data Set and Model Specification

Pouplier and Hoole (2016) study the assimilation of the German /s/ and /sh/
sounds such as the final consonant sounds in “Kürbis” (English example:
“haggis”) and “Gemisch” (English example: “dish”), respectively. The research
question is how these sounds assimilate in fluent speech when combined across
words such as in “Kürbis-Schale” or “Gemisch-Salbe”, denoted as /s#sh/ and
/sh#s/ with # the word boundary. The 9 native German speakers in the study
repeated a set of 16 selected word combinations five times. Two different
types of functional data, i.e. acoustic (ACO) and electropalatographic (EPG)
data, were recorded for each repetition to capture the acoustic (produced
sound) and articulatory (tongue movements) aspects of assimilation over
(relative) time $t$ within the consonant combination.

Each functional index varies roughly between $+1$ and $-1$ and measures how
similar the articulatory or acoustic pattern is to its reference patterns for
the first ($+1$) and second ($-1$) consonant at every observed time point
(Cederbaum et al.; 2016). Without assimilation, the data are thus expected to
shift from positive to negative values in a sinus-like form (see Figure 3).
The data set contains 707 bivariate functional observations with differently
spaced grids of evaluation points per curve and dimension, with the number of
evaluation points ranging from 22 to 59 with a median of 35. Note that the
consonant assimilation data are unaligned as registration of the time domain
would mask transition speeds between the consonants, which are an interesting
part of assimilation.

Figure 3: Index curves of the consonant assimilation data set for both ACO and
EPG data as a function of standardized time $t$ (grey curves). For every
consonant order, three randomly selected observations have been highlighted.
Every colour corresponds to one multivariate observation, i.e. one observation
consists of two index curves.

For comparability, we follow the model specification of Cederbaum et al.
(2016), who analyse only the ACO dimension and ignore the second mode EPG. Our
specified multivariate model is similar to (10) with $i=1,...,9$ the speaker
index, $j=1,...,16$ the word combination index, $h=1,...,H_{ij}$ the
repetition index and $t\in[0,1]$ relative time. Note that the nested effect
$\bm{C}_{ij}(t)$ is replaced by the crossed random effect $\bm{C}_{j}(t)$
specific to the word combinations. The additive predictor
$\bm{\mu}(\bm{x}_{j},t)$ now contains eight partial effects: the functional
intercept plus main and interaction effects of scalar covariates describing
characteristics of the word combination such as the order of the consonants
/s/ and /sh/. The white noise measurement error $\bm{\epsilon}_{ijht}$ is
assumed to follow a zero-mean bivariate normal distribution with diagonal
covariance matrix
$\operatorname{diag}(\sigma^{2}_{\text{ACO}},\sigma^{2}_{\text{EPG}})$. The
basis and penalty specifications follow the univariate analysis in Cederbaum
et al. (2016). Given different sampling mechanisms, we also compare the
multiFAMM based on weighted and unweighted scalar products for the MFPCA.

#### Results

The multivariate analysis supports the findings of Cederbaum et al. (2016)
that assimilation is asymmetric (different mean patterns for /s#sh/ and
/sh#s/). Overall, the estimated fixed effects are similar across dimensions as
well as comparable to the univariate analysis. Hence, the multivariate
analysis indicates that previous results for the acoustics are consistently
found also for the articulation. Compared to univariate analyses, our approach
reduces the number of FPC basis functions and thus the number of parameters in
the analysis. The multiFAMM can improve the model fit and can provide smaller
CBs for the ACO dimension compared to the univariate model in Cederbaum et al.
(2016) due to the strong cross-correlation between the dimensions. We find
similar modes of variation for the multivariate and the univariate analysis as
well as across dimensions. In particular, the word combination-specific random
effect $\bm{C}_{j}(t)$ is dropped from the model as much of the between-word
variation is already explained by the included fixed effects. The definition
of the scalar product has little effect on the estimated fixed effects but
changes the interpretation of the FPCs. Appendix C contains a more in depth
description of this application.

##  5 Simulations

###  5.1 Simulation Set-Up

We conduct an extensive simulation study to investigate the performance of the
multiFAMM depending on different model specifications and data settings (over
20 scenarios total), and to compare it to univariate regression models as
proposed by Cederbaum et al. (2016), estimated on each dimension
independently. Given the broad scope of analysed model scenarios, we refer the
interested reader to Appendix D for a detailed report and restrict the
presentation here to the main results.

We mimic our two presented data examples (Section 4) and simulate new data
based on the respective multiFAMM-fit. Each scenario consists of model fits to
500 generated data sets, where we randomly draw the number and location of the
evaluation points, the random scores, and the measurement errors according to
different data settings. The accuracy of the estimated model components is
measured by the root relative mean squared error (rrMSE) based on the
unweighted multivariate norm but otherwise as defined by Cederbaum et al.
(2016), see Appendix D.1. The rrMSE takes on (unbounded) positive values with
smaller values indicating a better fit.

###  5.2 Simulation Results

Figure 4 compares the rrMSE values over selected modeling scenarios based on
the consonant assimilation data. We generate a benchmark scenario (dark blue),
which imitates the original data without misspecification of any model
component. In particular, the number of FPCs is fixed to avoid truncation
effects. Comparing this scenario to the other blue scenarios illustrates the
importance of the number of FPCs in the accuracy of the estimation. Choosing
the truncation order via the proportion of univariate variance explained (Cut-
Off Uni) as in (7) gives models with roughly the same number of FPCs (mean
$\bm{B}:2.8,\bm{E}:5$) as is used for the data generation
($\bm{B}:3,\bm{E}:5$). The cut-off criterion based on the multivariate
variance (Cut-Off Mul) given by (6) results in more parsimonious models (mean
$\bm{B}:2.15,\bm{E}:4$) and thus considerably higher rrMSE values. The
increased variation in the rrMSE values can also be attributed to variability
in the truncation orders (cf. Appendix Figure 23), leading to a mixture
distribution. Comparing the benchmark scenario to more sparsely observed
functional data (ceteris paribus) suggests a lower estimation accuracy for the
Sparse Data scenario, especially for the curve-specific random effect
$\bm{E}_{ijh}(t)$ and resultingly the fitted curves $\bm{y}_{ijh}(t)$, but
pooling the information across functions helps the estimation of
$\bm{\mu}(\bm{x}_{ij},t)$ and $\bm{B}_{i}(t)$. In particular, the estimation
of the mean $\bm{\mu}(\bm{x}_{ij},t)$ is quite robust against the increased
uncertainty of these three scenarios. Only when the random scores are not
centred and decorrelated as in the benchmark scenario do we find an increase
in rrMSE values for the mean (Uncentred Scores). This corresponds to a
departure from the modeling assumptions likely to occur in practice when only
few levels of a random effect are available (here for the subject-specific
$\bm{B}_{i}(t)$). The model then has difficulties to correctly separate the
intercept in $\bm{\mu}(\bm{x}_{ij},t)$ and the random effects $\bm{B}_{i}(t)$.
The empirical (non-zero) mean of the $\bm{B}_{i}(t)$ is then absorbed by the
intercept in $\bm{\mu}(\bm{x}_{ij},t)$, resulting in higher rrMSE values for
both of these model terms. However, this shift does not affect the overall fit
to the data $\bm{y}_{ijh}(t)$ nor the estimation of the other fixed effects
(cf. Appendix Figure 31). Note that the rrMSE values of the Sparse Data and
Uncentred Scores scenarios are based on slightly different normalizing
constants (i.e. different true data) and cannot be directly compared except
for the mean.

Figure 4: rrMSE values of the fitted curves $\bm{y}_{ijh}(t)$, the mean
$\bm{\mu}(\bm{x}_{ij},t)$, and the random effects $\bm{B}_{i}(t)$ and
$\bm{E}_{ijh}(t)$ for different modeling scenarios. Scenarios coloured with
blue correspond to different model specifications in the same data setting.

Our simulation study thus suggests that basing the truncation orders on the
proportion of explained variation on each dimension (7) gives parsimonious and
well fitting models. If interest lies mainly in the estimation of fixed
effects, the alternative cut-off criterion based on the total variation in the
data (6) allows even more parsimonious models at the cost of a less accurate
estimation of the random effects and overall model fit. Furthermore, the
results presented in Appendix D show that the mean estimation is relatively
stable over different model scenarios including misspecification of the
measurement error variance structure or of the multivariate scalar product, as
well as in scenarios with strong heteroscedasticity across dimensions. In our
benchmark scenario, the CBs cover the true effect $89-94\%$ of the time but
coverage can further decrease with additional uncertainty e.g. about the
number of FPCs. Overall, the covariance structure such as the leading FPCs can
be recovered well, also for a nested random effect such as in the snooker
training application. The comparison to the univariate modeling approach
suggests that the multiFAMM can improve the mean estimation but is especially
beneficial for the prediction of the random effects while reducing the number
of parameters to estimate. In some cases like strong heteroscedasticity,
including weights in the multivariate scalar product might further improve the
modeling.

##  6 Discussion

The proposed multivariate functional regression model is an additive mixed
model, which allows to model flexible covariate effects for sparse or
irregular multivariate functional data. It uses FPC based functional random
effects to model complex correlations within and between functions and
dimensions. An important contribution of our approach is estimating the
parsimonious multivariate FPC basis from the data. This allows us to account
not only for auto-covariances, but also for non-trivial cross-covariances over
dimensions, which are difficult to adequately model using alternative
approaches such as parametric covariance functions like the Matèrn family or
penalized splines, which imply a parsimonious covariance only within but not
necessarily between functions. As a FAMM-type regression model, a wide range
of covariate effect types is available, also providing pointwise CBs. Our
applications show that the multiFAMMs can give valuable insight into the
multivariate correlation structure of the functions in addition to the mean
structure.

An apparent benefit of multivariate modeling is that it allows to answer
research questions simultaneously relating to different dimensions. In
addition, using multivariate FPCs reduces the number of parameters compared to
fitting comparable univariate models while improving the random effects
estimation by incorporating the cross-covariance in the multivariate analysis.
The added computational costs are small: For our multimodal application, the
multivariate approach prolongs the computation time by only five percent (104
vs. 109 minutes on a 64-bit Linux platform).

We find that the average point-wise coverage of the point-wise CBs can in some
cases lie considerably below the nominal value. There are two main reasons for
this: One, the CBs presented here do not incorporate the uncertainty of the
eigenfunction estimation nor of the smoothing parameter selection. Two,
coverage issues can arise in (scalar) mixed models, if effect functions are
estimated as constant when in truth they are not (e.g. Wood; 2017; Greven and
Scheipl; 2016). To resolve these issues, further research on the level of
scalar mixed models might be needed. A large body of research covering CB
estimation for functional data (e.g. Goldsmith et al.; 2013; Choi and
Reimherr; 2018; Liebl and Reimherr; 2019) suggests that the construction of
CBs is an interesting and complex problem, also outside of the FAMM framework.

It would be interesting to extend the multiFAMM to more general scenarios of
multivariate functional data such as observations consisting of functions with
different dimensional domains, e.g. functions over time and images as in Happ
and Greven (2018). This would require adapting the estimation of the
univariate auto-covariances for spatial arguments $t,t^{\prime}$. Exploiting
properties of dense functional data, such as the block structure of design
matrices for functions observed on a grid, could help to reduce computational
cost in this case. Future research could further generalize the covariance
structure of the multiFAMM by allowing for additional covariate effects. In
our snooker training application, for example, a treatment effect of the
snooker training might show itself in the form of reduced intra-player
variance (cf. Backenroth et al.; 2018). Ideas from distributional regression
could be incorporated to jointly model the mean trajectories and covariance
structure conditional on covariates.

## Acknowledgements

We thank Timon Enghofer, Phil Hoole, and Marianne Pouplier for providing
access to their data and for fruitful discussions. We also thank Lisa Steyer
for contributing the data registration of the snooker training data and two
anonymous reviewers for their helpful suggestions. Sonja Greven, Almond
Stöcker, and Alexander Volkmann were funded by grant GR 3793/3-1 from the
German research foundation (DFG). Fabian Scheipl was funded by the German
Federal Ministry of Education and Research (BMBF) under Grant No. 01IS18036A.

## References

  * (1)
  * Backenroth et al. (2018) Backenroth, D., Goldsmith, J., Harran, M. D., Cortes, J. C., Krakauer, J. W. and Kitago, T. (2018).  Modeling motor learning using heteroscedastic functional principal components analysis, Journal of the American Statistical Association 113(523): 1003–1015.
  * Carroll et al. (2021) Carroll, C., Müller, H.-G. and Kneip, A. (2021).  Cross-component registration for multivariate functional data, with application to growth curves, Biometrics 77(3): 839–851.
  * Cederbaum et al. (2016) Cederbaum, J., Pouplier, M., Hoole, P. and Greven, S. (2016).  Functional linear mixed models for irregularly or sparsely sampled data, Statistical Modelling 16(1): 67–88.
  * Cederbaum et al. (2018) Cederbaum, J., Scheipl, F. and Greven, S. (2018).  Fast symmetric additive covariance smoothing, Computational Statistics & Data Analysis 120: 25–41.
  * Chiou et al. (2014) Chiou, J.-M., Chen, Y.-T. and Yang, Y.-F. (2014).  Multivariate functional principal component analysis: A normalization approach, Statistica Sinica 24(4): 1571–1596.
  * Choi and Reimherr (2018) Choi, H. and Reimherr, M. (2018).  A geometric approach to confidence regions and bands for functional parameters, Journal of the Royal Statistical Society: Series B (Statistical Methodology) 80(1): 239–260.
  * Di et al. (2009) Di, C.-Z., Crainiceanu, C. M., Caffo, B. S. and Punjabi, N. M. (2009).  Multilevel functional principal component analysis, The Annals of Applied Statistics 3(1): 458.
  * Eilers and Marx (1996) Eilers, P. H. and Marx, B. D. (1996).  Flexible smoothing with B-splines and penalties, Statistical Science 11(2): 89–102.
  * Enghofer (2014) Enghofer, T. (2014).  Überblick über die Sportart Snooker, Entwicklung eines Muskeltrainings und Untersuchung dessen Einflusses auf die Stoßtechnik.  Unpublished Zulassungsarbeit, Technische Universität München.
  * Goldsmith et al. (2013) Goldsmith, J., Greven, S. and Crainiceanu, C. (2013).  Corrected confidence bands for functional data using principal components, Biometrics 69(1): 41–51.
  * Goldsmith and Kitago (2016) Goldsmith, J. and Kitago, T. (2016).  Assessing systematic effects of stroke on motor control by using hierarchical function-on-scalar regression, Journal of the Royal Statistical Society: Series C (Applied Statistics) 65(2): 215–236.
  * Goldsmith et al. (2016) Goldsmith, J., Scheipl, F., Huang, L., Wrobel, J., Gellar, J., Harezlak, J., McLean, M. W., Swihart, B., Xiao, L., Crainiceanu, C. and Reiss, P. T. (2016).  refund: Regression with Functional Data.  R package version 0.1-16.
https://CRAN.R-project.org/package=refund

  * Greven and Scheipl (2016) Greven, S. and Scheipl, F. (2016).  Comment, Journal of the American Statistical Association 111(516): 1568–1573.
  * Greven and Scheipl (2017) Greven, S. and Scheipl, F. (2017).  A general framework for functional regression modelling, Statistical Modelling 17(1-2): 1–35, 100–115.
  * Happ and Greven (2018) Happ, C. and Greven, S. (2018).  Multivariate functional principal component analysis for data observed on different (dimensional) domains, Journal of the American Statistical Association 113(522): 649–659.
  * Jacques and Preda (2014) Jacques, J. and Preda, C. (2014).  Model-based clustering for multivariate functional data, Computational Statistics & Data Analysis 71: 92–106.
  * Li et al. (2020) Li, C., Xiao, L. and Luo, S. (2020).  Fast covariance estimation for multivariate sparse functional data, Stat 9(1): e245.
  * Li et al. (2017) Li, J., Huang, C., Hongtu, Z. and Alzheimer’s Disease Neuroimaging Initiative (2017).  A functional varying-coefficient single-index model for functional response data, Journal of the American Statistical Association 112(519): 1169–1181.
  * Liebl and Reimherr (2019) Liebl, D. and Reimherr, M. (2019).  Fast and fair simultaneous confidence bands for functional parameters, arXiv preprint arXiv:1910.00131 .
  * Liu et al. (2020) Liu, Y., Yan, B., Merikangas, K. and Shou, H. (2020).  Graph-fused multivariate regression via total variation regularization, arXiv preprint arXiv:2001.04968 .
  * Morris (2017) Morris, J. S. (2017).  Comparison and contrast of two general functional regression modelling frameworks, Statistical Modelling 17(1-2): 59–85.
  * Park and Ahn (2017) Park, J. and Ahn, J. (2017).  Clustering multivariate functional data with phase variation, Biometrics 73(1): 324–333.
  * Peng and Paul (2009) Peng, J. and Paul, D. (2009).  A geometric approach to maximum likelihood estimation of the functional principal components from sparse longitudinal data, Journal of Computational and Graphical Statistics 18(4): 995–1015.
  * Pouplier and Hoole (2016) Pouplier, M. and Hoole, P. (2016).  Articulatory and acoustic characteristics of german fricative clusters, Phonetica 73(1): 52–78.
  * R Core Team (2020) R Core Team (2020).  R: A Language and Environment for Statistical Computing, R Foundation for Statistical Computing, Vienna, Austria.
https://www.R-project.org

  * Ramsay and Silverman (2005) Ramsay, J. O. and Silverman, B. W. (2005).  Functional data analysis, second edn, Springer Science & Business Media.
  * Reiss and Xu (2020) Reiss, P. T. and Xu, M. (2020).  Tensor product splines and functional principal components, Journal of Statistical Planning and Inference 208: 1–12.
  * Scheipl et al. (2015) Scheipl, F., Staicu, A.-M. and Greven, S. (2015).  Functional additive mixed models, Journal of Computational and Graphical Statistics 24(2): 477–501.
  * Steyer et al. (2021) Steyer, L., Stöcker, A. and Greven, S. (2021).  Elastic analysis of irregularly or sparsely sampled curves, arXiv preprint arXiv:2104.11039 .
  * Uludağ and Roebroeck (2014) Uludağ, K. and Roebroeck, A. (2014).  General overview on the merits of multimodal neuroimaging data fusion, Neuroimage 102: 3–10.
  * Volkmann (2021) Volkmann, A. (2021).  multifamm: Multivariate Functional Additive Mixed Models.  R package version 0.1.1.
https://cran.r-project.org/web/packages/multifamm

  * Wood (2017) Wood, S. N. (2017).  Generalized Additive Models: An Introduction with R, 2 edn, Chapman and Hall/CRC Press.
  * 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(470): 577–590.
  * Zhu et al. (2012) Zhu, H., Li, R. and Kong, L. (2012).  Multivariate varying coefficient model for functional responses, Annals of Statistics 40(5): 2634.
  * Zhu et al. (2017) Zhu, H., Morris, J. S., Wei, F. and Cox, D. D. (2017).  Multivariate functional response regression, with application to fluorescence spectroscopy in a cervical pre-cancer study, Computational Statistics & Data Analysis 111: 88–101.
  * Zhu et al. (2016) Zhu, H., Strawn, N. and Dunson, D. B. (2016).  Bayesian graphical models for multivariate functional data, The Journal of Machine Learning Research 17(1): 7157–7183.

##  Appendix A Derivation of Variance Decompositions

This is the derivation of the variance decompositions (6) and (7) from the
model equations (1) and (2). Following the argumentation of Cederbaum et al.
(2016), the variation of the response can be decomposed per dimension using
iterated expectations as

$\displaystyle\int_{\mathcal{I}}\mathrm{Var}\big{(}y_{i}^{(d)}(t)\big{)}dt=$
$\displaystyle\sum_{j=1}^{q}\sum_{v=1}^{V_{U_{j}}}z_{ijv}^{2}\int_{\mathcal{I}}\mathrm{Var}\big{(}U_{jv}^{(d)}(t)\big{)}dt+\int_{\mathcal{I}}\mathrm{Var}\big{(}E_{i}^{(d)}(t)\big{)}dt+\int_{\mathcal{I}}\mathrm{Var}\big{(}\epsilon_{it}^{(d)}\big{)}dt$
$\displaystyle=$
$\displaystyle\sum_{j=1}^{q}\sum_{m=1}^{\infty}\nu_{U_{j}m}\underbrace{\int_{\mathcal{I}}\psi_{U_{j}m}^{(d)}(t)\psi_{U_{j}m}^{(d)}(t)dt}_{\textstyle=||\psi_{U_{j}m}^{(d)}||^{2}}$
$\displaystyle+\sum_{m=1}^{\infty}\nu_{Em}\underbrace{\int_{\mathcal{I}}\psi_{Em}^{(d)}(t)\psi_{Em}^{(d)}(t)dt}_{\textstyle=||\psi_{Em}^{(d)}||^{2}}+\int_{\mathcal{I}}\sigma_{d}^{2}dt$
$\displaystyle=$
$\displaystyle\sum_{j=1}^{q}\sum_{m=1}^{\infty}\nu_{U_{j}m}||\psi_{U_{j}m}^{(d)}||^{2}+\sum_{m=1}^{\infty}\nu_{Em}||\psi_{Em}^{(d)}||^{2}+\sigma_{d}^{2}|\mathcal{I}|$

assuming in the second equality that each observation is in exactly one group
in the $j$th grouping layer, i.e. exactly one indicator
$z_{ijv}=z_{ijv}^{2},v=1,\dots,V_{U_{j}}$ is one for each $i$ and $j$, and
also using in the second equality that the scores are uncorrelated across
eigenfunctions, in particular
$\text{Cov}(\rho_{U_{j}vm},\rho_{U_{j}vm^{\prime}})=0$ if $m\neq m^{\prime}$.

Similarly, the total variance as measured by the sum of the univariate
variances can be decomposed as

$\displaystyle\mathbb{E}\left(|||\bm{y}_{i}-\bm{\mu}(\bm{x}_{i})|||^{2}\right)=$
$\displaystyle\sum_{d=1}^{D}w_{d}\int_{\mathcal{I}}\mathrm{Var}\big{(}y_{i}^{(d)}(t)\big{)}dt$
$\displaystyle=$
$\displaystyle\sum_{j=1}^{q}\sum_{d=1}^{D}\sum_{v=1}^{V_{U_{j}}}z_{ijv}^{2}w_{d}\int_{\mathcal{I}}\mathrm{Var}\big{(}U_{jv}^{(d)}(t)\big{)}dt+\sum_{d=1}^{D}w_{d}\int_{\mathcal{I}}\mathrm{Var}\big{(}E_{i}^{(d)}(t)\big{)}dt$
$\displaystyle+\sum_{d=1}^{D}w_{d}\int_{\mathcal{I}}\mathrm{Var}\big{(}\epsilon_{it}^{(d)}\big{)}dt$
$\displaystyle=$
$\displaystyle\sum_{j=1}^{q}\sum_{m=1}^{\infty}\nu_{U_{j}m}\underbrace{\sum_{d=1}^{D}w_{d}\int_{\mathcal{I}}\psi_{U_{j}m}^{(d)}(t)\psi_{U_{j}m}^{(d)}(t)dt}_{\textstyle=|||\bm{\psi}_{U_{j}m}|||^{2}=1}$
$\displaystyle+\sum_{m=1}^{\infty}\nu_{Em}\underbrace{\sum_{d=1}^{D}w_{d}\int_{\mathcal{I}}\psi_{Em}^{(d)}(t)\psi_{Em}^{(d)}(t)dt}_{\textstyle=|||\bm{\psi}_{Em}|||^{2}=1}+\sum_{d=1}^{D}w_{d}\int_{\mathcal{I}}\sigma_{d}^{2}dt$
$\displaystyle=$
$\displaystyle\sum_{j=1}^{q}\sum_{m=1}^{\infty}\nu_{U_{j}m}+\sum_{m=1}^{\infty}\nu_{Em}+\sum_{d=1}^{D}w_{d}\sigma_{d}^{2}|\mathcal{I}|.$

##  Appendix B Additional Results for the Snooker Data

###  B.1 Data Description

In this section, we provide additional information of the study by Enghofer
(2014). The participants of the snooker training study volunteered if they
wanted to take part in the training group. The training schedule for this
treatment group recommended to perform a training two to three times a week
which consisted of three exercises that are supposed to improve snooker
specific muscular coordination. Note that these were autonomous trainings so
that there is no reliable information if the participants followed the
recommendations.

In order to analyse the snooker shot of maximal force, the participants were
asked to place themselves in the typical playing position. Yet no exact
measure of the distance between the tip of their cue and the cue ball is
available and thus the time of the impulse onto the cue ball cannot be exactly
determined from the data. The participants were videotaped and the open source
software Kinovea (https://www.kinovea.org) was used to manually locate points
of interest (a participant’s shoulder, elbow, and hand) and track them on a
two-dimensional grid over the course of the video. Figure 5 shows the
univariate functions that make up the two-dimensional trajectories of the
participants’ elbow, hand, and shoulder. The univariate functions for the
x-axis of the hand show that for most observed curves, the time of impact
might lie around $0.75$ of the relative time of the snooker shot. That is when
the hand moves into the positive range on the $x$ dimension. For their shot,
the snooker players can either try to fix their elbow or move the elbow
dynamically. In order to move the hand in a straight line (piston stroke), the
elbow has to move dynamically down and up when the hand is drawn back and
accelerated towards the cue ball. The latter part of the movement trajectory
(after the cue hits the cue ball) does not impact the snooker shot itself.
Note that a long final downwards movement of an elbow trajectory indicates
that the hand is not stopped at the chest but is rather pulled through towards
the shoulder. This can give insight into a player’s stance at the snooker
table, e.g. if the upper body is close to the table or if it is more erect.

Figure 5: Univariate functions of the snooker trajectories as functions of the
standardized time $t$. The value of the function corresponds to the position
of the elbow, hand, and shoulder on the $x$ and $y$ axes at the associated
time point.

Note that the univariate functions of Figure 5 show different characteristics
over the different dimensions, and a first-order difference penalty penalizing
deviations from a constant function seems to be a sensible default choice.

###  B.2 Data Preprocessing

In order to account for differences in body height between the participants,
we first rescale the observed coordinate locations by the median distance
between hand and elbow. As we are interested in the movement trajectory
irrespective of phase variation, i.e. independent of the exact timing of
different parts of the stroke, we apply an elastic registration approach to
the data which aligns the two dimensional trajectory of the hand to its
respective mean trajectory across all players and shots (Steyer et al.; 2021).
We then reparameterize the time for the elbow and shoulder trajectories
according to the results of the hand alignment. Other registration approaches
are possible and we provide both the original and the reparameterized time of
the data. Due to the high frame rate of the high speed camera and a
comparatively rough grid of the tracking software, the resulting multivariate
functional data are dense with redundant information. The following subsection
describes the coarsening method in detail. The coarsening reduces the data set
from roughly 400,000 to only 56,910 scalar observations. Note that 5
functional observations are lost due to technical problems during the
recording. Even bivariate trajectories within one multivariate functional
observation can have different numbers of scalar observations due to
inconsistencies in the tracking programme.

#### Coarsening Method for the Snooker Data

Functional response variables as given by the snooker trajectories typically
exhibit very high autocorrelation. Thus, up to a certain point, removing
single curve measurements will inflict only very limited information loss
while reducing computational cost in the model fit at least linearly. We
therefore coarsen the snooker trajectories as a preprocessing step, as they
were originally sampled “overly” densely in that sense. In particular, there
are also evaluation points that do not carry additional information as they
measure the same location as the evaluation points right before and after
(e.g. because the hand does not move noticeably over these three or more time
points). The coarsening should a) retain original data points via subsampling
to avoid pre-smoothing losses and b) be optimal in the sense that as much
information as possible is preserved given the subsample size.

Finding an optimal subsample, with respect to a general criterion and for all
possible subsamples, is however typically very computationally demanding in
itself. We therefore propose a fast greedy coarsening algorithm recursively
discarding single points from the sample. The ultimate aim of the algorithm is
to find a subsample polygon (Fig. 6) representing the whole data as well as
possible.

Figure 6: An example subsample polygon of length four through some curve
measurements. The polygon canonically connects points in the subsample
according to the order of the indices $t_{i}$ of the curve measurements
$y(t_{i})$ in the subsample. Figure 7: One iteration of the proposed fast
greedy coarsening algorithm: Given the current subsample polygon (top), the
loss $\Delta^{(k)}(t)$ inflicted when approximating $y(t)$ by the line segment
between its neighbors is computed for each index $t_{2},\dots,t_{5}$ (middle).
The data point with the least approximation loss is discarded (bottom). Note
that $\Delta^{(k+1)}(t)$ has to be re-computed only for the two adjacent
points. All others can be passed forward.

One step of the coarsening algorithm described in the following is illustrated
in Figure 7. Let $\mathcal{T}_{0}=\\{t_{1},\dots,t_{m}\\}\subset\mathcal{I}$
be the set of evaluation points $t_{1}<\dots<t_{m}$ of a curve
$y:\mathcal{I}\rightarrow\mathbb{R}^{2}$. We select a sequence of subsamples
$\mathcal{T}_{0}\supset\mathcal{T}_{1}\supset\mathcal{T}_{2}\supset\dots$ by
recursively discarding the point $t$ from $\mathcal{T}_{k}$ where $y(t)$ is
closest to the line segment $[y(t_{-}),y(t_{+})]$ between its adjacent points
at $t_{-}=\max\\{s\in\mathcal{T}_{k}:s<t\\}$ and
$t_{+}=\min\\{s\in\mathcal{T}_{k}:s>t\\}$. Note that this is only well-defined
for
$t\in\mathring{\mathcal{T}}_{k}=\mathcal{T}_{k}\setminus\\{\min\mathcal{T}_{k},\max\mathcal{T}_{k}\\}$,
in the “interior” of $\mathcal{T}_{k}$, and the first and last point $t_{1}$
and $t_{m}$ are kept throughout the algorithm.

In more detail, we proceed in the $k$th step of the algorithm as follows

  1. 1.

$\forall t\in\mathring{\mathcal{T}}_{k}$ compute the step-to-step quadratic
error

$\displaystyle\Delta^{(k)}(t)=\min_{p\in[y(t_{-}),y(t_{+})]}(y(t),p)^{2}.$

For computing the quadratic distance $\Delta^{(k)}(t)$, we have to distinguish
the case where the projection of $y(t)$ on the line through $y(t_{-})$ and
$y(t_{+})$ lies between the two points and the case where it lies on either
side of them. To do so, define $y^{\prime}(t)=y(t)-y(t_{-})$ and
$y^{\prime}_{+}=y^{\prime}(t_{+})$, as well as the unit vector
$u=\frac{y^{\prime}_{+}}{\|y^{\prime}_{+}\|}$ pointing from $y(t_{-})$ to
$y(t_{+})$. Then we can compute

$\displaystyle\Delta^{(k)}(t)=\begin{cases}\|y^{\prime}(t)\|^{2}&\text{ if
}\langle y^{\prime}(t),u\rangle\leq 0\\\ \|y^{\prime}(t)-\langle
y^{\prime}(t),u\rangle u\|^{2}&\text{ if }\|y^{\prime}_{+}\|\geq\langle
y^{\prime}(t),u\rangle>0\\\ \|y^{\prime}(t)-y^{\prime}_{+}\|^{2}&\text{ if
}\langle y^{\prime}(t),u\rangle>\|y^{\prime}_{+}\|\end{cases}.$

  2. 2.

Choose
$t^{(k)}\in\operatorname{argmin}_{t\in\mathring{\mathcal{T}}_{k}}\Delta^{(k)}(t)$
and set $\mathcal{T}_{k+1}=\mathcal{T}_{k}\setminus\\{t^{(k)}\\}$.

This procedure may be repeated until $\mathcal{T}_{k+1}$ is of the desired
sample size. However, we consider an error based stopping criterion more
convenient in most cases. We suggest to specify a threshold for the cumulative
step-to-step error $S_{k}=\sum_{\kappa=1}^{k}\Delta^{(\kappa)}(t^{(\kappa)})$
or a relative version $R_{k}=S_{k}/\bar{\Delta}$ of it, relative to
$\bar{\Delta}=\frac{1}{m-2}\sum_{\mathring{\mathcal{T}}}\Delta^{(\infty)}(t)$
the mean quadratic distance to the line segment $[y(t_{1}),y(t_{m})]$ between
the two only points left at the ultimate subsample $\mathcal{T}_{\infty}$.

In our experience, a threshold $R^{*}=10^{-4}$ for $R_{k}$ yields a close to
lossless subsample from a visual point of view. For the model, $R^{*}=0.003$
was chosen to limit computational demands. While the multiFAMM would allow for
dimension specific evaluations of the curves, the warping algorithm applied as
part of the preprocessing does not. Thus, we decided to apply the coarsening
algorithm to the hand trajectories, which we consider the most informative
component with the best signal to noise ratio. Unselected observed time points
for each hand trajectory are then also dropped for the corresponding
trajectories of the shoulder and elbow. As a result, the coarsened data retain
the characteristic of the snooker trajectory data such that for a given time
point, we observe the location of (almost) all points of interest (i.e. the
grid of evaluation points is identical over the dimensions for a given
observation).

###  B.3 Additional Results of the Analysis

#### Analysis of the Variance Components

The multiFAMM estimates a total of 61 univariate FPCs (20 each for $B$ and
$C$, and 21 for $E$), where each process is represented by three to five
univariate FPCs on each dimension. Table 1 shows that the three points of
interest (hand, elbow, and shoulder) all show a similar amount of variation.
The source of variation, however, differs in interpretation. While the
variation in the hand trajectories stems from differing movements, the
variation in the shoulder mainly reflects different positioning. Applying
higher weights to dimensions associated with the hand, for example, can shift
the focus of the MFPCA to favor movement based variation of the hand
trajectories in the analysis. From Table 1 it is also apparent that movements
or positions along the horizontal $x$-axis contribute more to the variation in
the data than movements or positions along the vertical $y$-axis. This also
suggests that variation in the $x$ direction is the driving factor of the
MFPCA. By assigning higher weights in the scalar product to dimensions
associated with the $y$-axis, the analyst can “distort” the natural grid of
observations to balance out the different variation in the axes. Note that we
use an unweighted scalar product in the analysis. Also keep in mind that the
variation in Table 1 is calculated based on estimation step 1.

Based on the cut-off criterion (6), 16 multivariate FPCs are chosen to explain
$95\%$ of total variation. For a similar amount of variance explained, 23
univariate FPCs would be needed, which would give a more complex model that
contains redundancies and ignores the multivariate nature of the data. Figures
8 and 9 show i.a. the two most prominant modes of variation in the snooker
training data. As described in the main part, the FPC $\bm{\psi}_{C1}$ seems
to represent variation in the starting position (explaining about $27\%$ of
total variation. The second leading FPC $\bm{\psi}_{B1}$ (explaining about
$15\%$ of total variation) shows subject-specific variation related to the two
snooker techniques identified by Enghofer (2014): The red trajectory ($+$)
shows that the elbow moves first down and then up to draw the hand back and
accelerate it towards the cue ball. This then vertically contracts the hand
trajectory compared to fixing the elbow during the acceleration phase (blue
$-$), which allows the hand to swing more freely (pendulum stroke). We also
find that players with a personal tendency towards the pendulum stroke seem to
not stop their hands at their chest. Note that for $\bm{\psi}_{B1}$ in Figure
8, the red trajectory is overlapping and thus masks the down-up-and-down
movement of the elbow.

Table 1: Total variation of the centred responses per univariate dimension and overall. The variation is calculated as the sum of non-negative univariate eigenvalues and the dimension-specific measurement error variance as given from Step $1$ in the estimation of the multiFAMM as described in Section 3.1 | elbow.x | elbow.y | hand.x | hand.y | shoulder.x | shoulder.y | Total
---|---|---|---|---|---|---|---
Variation | 0.012 | 0.004 | 0.015 | 0.001 | 0.014 | 0.008 | 0.055
Proportion | 22% | 7% | 28% | 2% | 26% | 14% | 100%

Figure 8 shows the multivariate FPCs of the subject-specific functional random
intercept included in the multiFAMM. $\bm{\psi}_{B2}(t)$ (explaining about 9%
of the total variation) seems to capture variation in the positioning of the
shoulder, size differences of the upper arm, and the final part of the
movement trajectories after the cue ball is hit. Figure 9 shows the
multivariate FPCs of the subject-and-session-specific random intercept. The
second FPC (explaining about $6\%$ of the total variation) shifts the relative
positioning of the elbow and shoulder. Figure 10 shows the multivariate FPCs
of the curve-specific random intercept. The leading FPC (explaining about
$7\%$ of total variation) captures variation in the starting position of elbow
and shoulder, and in the final hand movement.

Figure 8: FPCs for the subject-specific functional random effect
$\bm{B}_{i}(t)$ with the respective proportion of variance explained. The
black solid line represents the mean trajectory to which a suitable multiple
($2\sqrt{\nu_{B\cdot}}$) of the FPC is added (red $+$) and subtracted (blue
$-$). The start of the trajectories are marked with an asterisk. Figure 9:
FPCs for the subject-and-session-specific functional random effect
$\bm{C}_{ij}(t)$ with the respective proportion of variance explained. The
black solid line represents the mean trajectory to which a suitable multiple
($2\sqrt{\nu_{C\cdot}}$) of the FPC is added (red $+$) and subtracted (blue
$-$). The start of the trajectories are marked with an asterisk. Figure 10:
FPCs for the curve-specific functional random effect $\bm{E}_{ijh}(t)$ with
the respective proportion of variance explained. The black solid line
represents the mean trajectory to which a suitable multiple
($2\sqrt{\nu_{E\cdot}}$) of the FPC is added (red $+$) and subtracted (blue
$-$). The start of the trajectories are marked with an asterisk.

#### Analysis of the Estimated Effect Functions

Figure 11: Estimated covariate effect functions for skill (left) and treatment
effect (right). The central plot shows the effect of the coefficient function
in parenthesis in the title (red) on the two dimensional trajectories for the
reference (without paranthesis, blue). The start of the trajectories are
marked with an asterisk. The marginal plots show the estimated univariate
effect functions (black) with pointwise 95% confidence intervals (dotted) and
the reference (blue dashed).

The left panel of Figure 11 shows the estimated covariate effect functions for
the covariate skill. In addition to the effect described in the main part, we
point out that for skilled players, the shoulder is positioned closer to the
table as well as to the body all other things being equal. The movement
trajectories in the right panel of Figure 11 are composed of the estimated
effect functions of intercept, treatment group, and session (blue
trajectories), so that the interaction effect can be separated and interpreted
(red trajectories). We do not find strong evidence for differences in the
displayed mean trajectories, nor in the univariate effects (marginal plots).
This suggests that the training programme did not considerably change the
participants’ mean movement trajectories.

Figure 12 shows the estimated intercept (left) and effect functions of the
covariates for treatment group (center) and session (right). The intercept
(scalar plus functional) gives the mean movement trajectories (dark blue) in
the reference group, i.e. an unskilled snooker player in the control group
with a shot in the first session. In the middle panel, the red trajectories
show the estimated effect of the treatment group added to the intercept. The
marginal plots around the movement trajectories show the estimated univariate
effects and their pointwise 95% confidence intervals. We find only minor
differences between movement trajectories in treatment and control group. The
hand trajectories of the treatment group seem to be slightly higher and
further along the $x$-axis than the control group, all other things being
equal. The right panel compares the mean trajectory for the reference group of
players in the first (blue) to the second (red) session. The univariate
estimated covariate effects on the $y$-axis seem to indicate a slight shift in
the vertical direction of the trajectories. Keep in mind that the trajectories
of the hand have been centred to the origin before the analysis.

Figure 12: Estimated functional fixed effects of intercept (left), covariate
treatment group (center), and covariate session (right). The main plots show
the effect of the coefficient function in parenthesis (red) on the two
dimensional trajectories for the intercept (blue). The start of the
trajectories are marked with an asterisk. The marginal plots show the
estimated univariate effect functions (black) with pointwise 95% confidence
intervals (dotted) and the reference (grey/blue dashed).

#### Model Diagnostics and Sensitivity Analysis

We use the univariate root relative mean squared error (urrMSE) defined in
Section 5 as a criterion for model fit between the different dimensions. Table
2 shows that the model fits comparatively well on all dimensions except for
the $x$-axis of the elbow and the $y$-axis of the hand. In order to get an
impression of the quality of the model fit, Figure 13 shows selected fitted
trajectories in solid red together with the observed trajectories in dashed
grey. We choose to present the quantiles from the contribution of each
observation to the multivariate root relative mean squared error (mrrMSE),
i.e. quantiles from $|||\bm{\zeta}_{i}-\hat{\bm{\zeta}}_{i}|||^{2}$ of the
multivariate functional trajectory $\bm{\zeta}_{i},i=1,...,N$ and the
corresponding estimate $\hat{\bm{\zeta}}_{i}$ (see Section 5).

The presented fitted trajectories suggest that the estimates might not always
overlap with the observed trajectories, which suggests residual structure in
the model residuals. The top panels of Figure 14 plot the scalar residuals
over time for each of the dimensions. Especially on the $y$-axis, we can
discern patterns in the residuals even though this type of graphic is prone to
overplotting. In order to investigate residual autocorrelation, we apply the
following ad-hoc method, which allows us to approximate the well-known concept
of an autocorrelation function in time series analysis. We first use fast
symmetric additive covariance smoothing (Cederbaum et al.; 2018) to estimate a
smoothed correlation matrix for the residuals. Then we calculate the means of
the off-diagonals, which corresponds to an approximated mean autocorrelation
for a given time lag. Figure 15 shows the results based on the residuals for
each dimension separately. Overall, we find that the model residuals (red)
show clear signs of autocorrelation, in particular up to a time lag of $0.25$.
However, compared to the autocorrelation of the original data (solid grey
line), we see a considerable reduction, which can be attributed to the
functional random effects in the model.

Table 3 also underlines the importance of the random effects in modeling the
snooker trajectories. We use the predictor variance as an indicator for
quantifying the effect of the fixed and random effects on the model fit.
Separately on each dimension, we calculate the variance of the partial
predictors on the data and give the respective proportions in Table 3. We find
that with exception of the reference mean ($\bm{f}_{0}(t)$), the partial
predictors of the fixed effects contribute generally little to the overall
predictor variance (highest proportions for $\bm{f}_{1}(t)$ with around $5\%$
on most dimensions). Even the reference mean contributes little to the
predictor variance on the shoulder and the $x$-axis of the elbow.
Consequently, we find that most variance in the partial predictors can be
assigned to the random effects. This suggests that the snooker training data
might be dominated by random processes with little explanatory power of the
available covariates.

Table 2: urrMSE values for the model fit on the snooker training data. | elbow.x | elbow.y | hand.x | hand.y | shoulder.x | shoulder.y
---|---|---|---|---|---|---
Main Model | 0.172 | 0.028 | 0.088 | 0.370 | 0.025 | 0.039
Sensitivity Analysis | 0.121 | 0.019 | 0.084 | 0.233 | 0.019 | 0.025
Figure 13: Fitted (red) and observed (blue) snooker training trajectories for selected observations (grey dashed). The qantiles given in the subtitles correspond to the contribution to the mrrMSE. Table 3: Proportion of predictor variance for each dimension. | $f_{0}^{(d)}(t)$ | $f_{1}^{(d)}(t)$ | $f_{2}^{(d)}(t)$ | $f_{3}^{(d)}(t)$ | $f_{4}^{(d)}(t)$ | $B_{i}^{(d)}(t)$ | $C_{ij}^{(d)}(t)$ | $E_{ijh}^{(d)}(t)$
---|---|---|---|---|---|---|---|---
elbow.x | 0.018 | 0.020 | 0.053 | 0.003 | 0.001 | 0.210 | 0.120 | 0.574
elbow.y | 0.381 | 0.048 | 0.006 | 0.002 | 0.000 | 0.144 | 0.261 | 0.157
hand.x | 0.862 | 0.022 | 0.014 | 0.000 | 0.000 | 0.022 | 0.050 | 0.030
hand.y | 0.755 | 0.037 | 0.012 | 0.004 | 0.000 | 0.103 | 0.076 | 0.012
shoulder.x | 0.007 | 0.047 | 0.013 | 0.001 | 0.019 | 0.184 | 0.144 | 0.585
shoulder.y | 0.009 | 0.068 | 0.007 | 0.018 | 0.001 | 0.105 | 0.512 | 0.280

Figure 14: Model residuals per dimension over time. The top panels correspond
to the model presented in the main analysis, the lower panels to the model
from the sensitivity analysis. Five observations have been removed from the
lower panels so that the scale of the plots is identical. Figure 15:
Approximated mean autocorrelation function over the time lag for the observed
data (grey), the fit of the model presented in the main analysis (red), and
the model from the sensitivity analysis (blue). The dashed line marks zero
autocorrelation.

We conduct a sensitivity analysis, where we increase the prespecified
proportion of variance to $0.99$. Additionally, we assume heteroscedasticity
in the model given the residual plot in Figure 14. The model then estimates
different measurement error variances with a predominantly large variance on
dimension $hand.x$. The model almost doubles the number of FPCs and includes
30 multivariate FPCs (B: $9$, C: $10$, E:$11$). The estimation of fixed
effects hardly change (results not shown) but we find a better fit to the data
(compare Table 2 and Figure 13). The scalar residuals in the lower panels of
Figure 14 seem to exhibit less structure for the model of the sensitivity
analysis but we can still discern patterns in the residuals. Considering the
residual autocorrelation, Figure 15 overall shows a reduction of
autocorrelation for small time lags compared to the main model. For larger
time lags we find a somewhat surprising increase of the approximated
autocorrelation function. Including a large number of comparatively trivial
FPCs might thus overall improve model fit while introducing new dependencies
in the residuals. Given that our main interest lies in the leading FPCs and
the fixed effects, we conclude that the sensitivity analysis yields similar
results with added computational complexity.

##  Appendix C Detailed Analysis of the Consonant Assimilation Data

###  C.1 Data and Model Description

#### Data Description

In this section we present more information on the analysis of the consonant
assimilation data. The ACO data were recorded via microphone at 32,768 Hz. In
addition, all speakers wore a custom-fit palate sensor with 62 electrodes that
measured the area where the tongue touched the palate during the articulation
in high temporal resolution (EPG data). These primary data were summarized and
transformed into a functional index over the time of articulation for the two
dimensions ACO and EPG separately. Each functional index measures how similar
the articulatory or acoustic pattern is to its reference patterns for the
first and second consonant at every observed time point. These similarity
indices vary roughly between $+1$ and $-1$, where the value $+1$ indicates
patterns close to the reference for the final target consonant (i.e. the
consonant ending the first word) and $-1$ for the initial target consonant
(i.e. the consonant beginning the second word). Due to the specifics of the
index calculation the index values can lie slightly outside the interval
$[-1,1]$. The curves on the ACO dimension are generally smoother than on the
EPG dimension, which reflects the index calculation: The ACO signal is much
richer in information because 256 continuously valued points are considered in
the calculation of the similarity index for every time point while only 62
binary points enter the index calculation for the EPG signal. Details on the
data generation and functional index computation can be found in Pouplier and
Hoole (2016) and Cederbaum et al. (2016).

#### Model Specification

For this application, we follow the model specification of Cederbaum et al.
(2016), who analyse only the ACO dimension of the data and ignore the second
mode of the consonant assimilation data. Our specified multivariate model is

$\displaystyle\bm{y}_{ijh}(t)=\bm{\mu}(\bm{x}_{ij},t)+\bm{B}_{i}(t)+\bm{C}_{j}(t)+\bm{E}_{ijh}(t)+\bm{\epsilon}_{ijht},$
(11)

with $i=1,...,9$ the speaker index, $j=1,...,16$ the word combination index,
and $h=1,...,H_{ij}$ the repetition index. Again, $\bm{B}_{i}(t)$ and
$\bm{E}_{ijh}(t)$ are the person-specific and curve-specific random
intercepts. $\bm{C}_{j}(t)$ is the word combination-specific random intercept,
which gives a crossed random effects structure. Note that here, the curve-
specific random component $\bm{E}_{ijh}(t)$ also captures speaker and word
interactions. Taking the different smoothness of the dimensions into account,
the white noise measurement error $\bm{\epsilon}_{ijht}$ is assumed to follow
a zero-mean bivariate normal distribution with diagonal covariance matrix
$\operatorname{diag}(\sigma^{2}_{\text{ACO}},\sigma^{2}_{\text{EPG}})$.

The additive predictor of the multiFAMM is specified as

$\displaystyle\bm{\mu}(\bm{x}_{ij},t)$
$\displaystyle=\bm{f}_{0}(t)+\texttt{order}_{j}\cdot\bm{f}_{1}(t)+\texttt{stress1}_{j}\cdot\bm{f}_{2}(t)+\texttt{stress2}_{j}\cdot\bm{f}_{3}(t)$
$\displaystyle\quad+\texttt{vowel}_{j}\cdot\bm{f}_{4}(t)+\texttt{order}_{j}\cdot\texttt{stress1}_{j}\cdot\bm{f}_{5}(t)$
$\displaystyle\quad+\texttt{order}_{j}\cdot\texttt{stress2}_{j}\cdot\bm{f}_{6}(t)+\texttt{order}_{j}\cdot\texttt{vowel}_{j}\cdot\bm{f}_{7}(t),$

with dummy covariates
$\texttt{order}_{j},\texttt{stress1}_{j},\texttt{stress2}_{j}$, and
$\texttt{vowel}_{j}$ indicating whether in the word combination /sh/ is
followed by /s/ (instead of the reference /s#sh/), the final target syllable
is not stressed, the initial target syllable is not stressed, and the vowel
context (i. e. the adjacent vowel for each consonant sound of interest) is
/a#i/ (instead of the reference /i#a/ as in “Gemisch Salbe”), respectively.
The functional intercept $\bm{f}_{0}(t)$ and the effect function
$\bm{f}_{1}(t)$ and their deviation from a sinus-like form or zero,
respectively, are especially interesting for studying assimilation.

For comparability, we follow the univariate analysis in Cederbaum et al.
(2016) by specifying cubic P-splines with third order difference penalty with
eight and five (marginal) basis functions for the covariate effect functions
and auto-covariance estimation, respectively. The MFPCA is based on an
unweighted scalar product. We choose the multivariate MFPCA truncation orders
so that at least 95% of the total variation in the data is explained as
presented in (6). To handle the heteroscedasticity, we use the weighted
regression approach. Given the different sampling mechanisms on the
dimensions, Appendix C.3 also contains an alternative analysis based on a
weighted scalar product for the MFPCA.

###  C.2 Results of the Model Presented in the Main Part

#### Analysis of the Variance Components

The univariate decomposition of the auto-covariances yields five univariate
FPCs for the random components $\bm{B}$ and $\bm{E}$ on both dimensions and
one and two FPCs for $\bm{C}$ on ACO and EPG, respectively. The cut-off
criterion based on the sum of total variation selects a total of eight FPCs.
The estimated multiFAMM then contains five multivariate FPCs for the smooth
residual (representing $64\%$ of total variation) and three for the subject-
specific random effects ($20\%$ of total variation) with $12\%$ of total
variation due to measurement error. With the chosen truncation order
$M_{C}=0$, the crossed random component $\bm{C}$ is effectively dropped from
the model as a lot of variation in the data is already explained by the
included fixed effects, all of which capture characteristics of the word
combinations ($8$ fixed effects for $16$ word combinations). Table 4 shows the
contribution of each random process to the total variation in the consonant
assimilation data according to the fitted multiFAMM. The first row gives the
eigenvalues of the random components $\bm{B}$ and $\bm{E}$ as well as the
estimated univariate error variances and the total amount of variation in the
data as computed in (6). The second and third row show the univariate $L^{2}$
norm of the multivariate eigenfunctions. The proportion of explained variance
$\pi$ is given in the final three rows. The proportions displayed in the
fourth and fifth row are computed according to (7) and the last row is
computed according to (6). Note that the fitted model uses the latter
multivariate cut-off criterion to determine the number of FPCs. With different
sampling mechanisms for different dimensions and thus different measurement
errors, it is advisable to check whether the proportion of explained variance
on each dimension (7) is adequate. For the consonant assimilation data, the
selected number of FPCs explains about $97\%$ of variation on EPG but only
$94\%$ on ACO, which we still deem acceptable. If a greater disparity is
revealed in an application, we recommend to use the proportion of explained
univariate variation as a cut-off criterion (in our case this would lead to
including a sixth FPC for component $E$).

Table 4: Variance components included in the multiFAMM. First row: Estimates
of eigenvalues, univariate error variances, total variation. Second (third)
row: Univariate norms of estimated FPCs on dimension ACO (EPG). Fourth (fifth)
row: Proportion of univariate variation explained on dimension ACO (EPG) by
eigenfunctions and error variance, total univariate variation explained. Sixth
row: Proportion of multivariate variation explained by eigenfunctions and
error variances, total multivariate variation explained.

| $B_{1}$ | $B_{2}$ | $B_{3}$ | $E_{1}$ | $E_{2}$ | $E_{3}$ | $E_{4}$ | $E_{5}$ | $\sigma_{ACO}^{2}$ | $\sigma_{EPG}^{2}$ | Total
---|---|---|---|---|---|---|---|---|---|---|---
Variation | $0.018$ | $0.009$ | $0.004$ | $0.060$ | $0.017$ | $0.012$ | $0.007$ | $0.003$ | $0.004$ | $0.014$ | $0.153$
$||\psi^{(ACO)}||^{2}$ | $0.169$ | $0.585$ | $0.642$ | $0.153$ | $0.217$ | $0.849$ | $0.178$ | $0.713$ | $--$ | $--$ | $--$
$||\psi^{(EPG)}||^{2}$ | $0.831$ | $0.415$ | $0.358$ | $0.847$ | $0.783$ | $0.151$ | $0.822$ | $0.287$ | $--$ | $--$ | $--$
$\pi^{(ACO)}$ | $0.068$ | $0.126$ | $0.052$ | $0.210$ | $0.082$ | $0.226$ | $0.026$ | $0.056$ | $0.094$ | $--$ | $0.940$
$\pi^{(EPG)}$ | $0.134$ | $0.036$ | $0.012$ | $0.467$ | $0.119$ | $0.016$ | $0.049$ | $0.009$ | $--$ | $0.127$ | $0.969$
$\pi$ | $0.115$ | $0.062$ | $0.023$ | $0.393$ | $0.108$ | $0.077$ | $0.042$ | $0.023$ | $0.027$ | $0.091$ | $0.961$

We present the estimated multivariate FPCs for component $\bm{B}$ in Figure
16. The estimated leading multivariate FPC $\bm{\psi}_{B1}(t)$ of the subject-
specific random effect $\bm{B}_{i}(t)$ (see left panels) shows a personal
tendency to pronounce the final or initial target syllable distinctly, which
explains roughly $12\%$ of the total variation. It changes the relative time
spent pronouncing each syllable and pulls the entire curve closer to the
reference sound of the preferred consonant. The second FPC (accounting for
$6\%$ variation), however, shows the individual tendency for assimilation. It
flattens or amplifies the sinus-like shape of the overall mean (black solid
line) thus rendering the speaker more or less prone to distinguish between the
two sounds. Note that in the univariate analysis for ACO presented in
Cederbaum et al. (2016), this mode of variation is identified as the leading
FPC of the person-specific random effect. The three leading FPCs
$\bm{\psi}_{E1}(t)$, $\bm{\psi}_{B1}(t)$, and $\bm{\psi}_{E2}(t)$ all show
greater univariate norms on the EPG dimension (Table 4) suggesting that this
dimension is the driving source of variation in the model – but this dimension
was ignored in the analysis of Cederbaum et al. (2016). Note that only
$\bm{\psi}_{B3}(t)$ impacts the dimensions differently in that more
assimilation on one dimension equals less assimilation on the other.

Figure 16: FPCs for the subject-specific functional random effect
$\bm{B}_{i}(t)$ with the respective proportion of variance explained. The
black solid line represents the mean trajectory to which a suitable multiple
($2\sqrt{\nu_{B\cdot}}$) of the FPC is added (red $+$) and subtracted (blue
$-$).

For the curve-specific random effect $\bm{E}_{ijh}(t)$ (Figure 17), the
leading FPC impacts primarily the final target consonant pulling it towards or
pushing it away from its reference sound. $\bm{\psi}_{E2}(t)$ has a similar
shape to $\bm{\psi}_{B2}(t)$. Note that the third leading FPC of $E$ affects
the mean function in opposite directions on to the two dimensions shifting the
curve up on one dimension and down on the other.

Figure 17: FPCs for the curve-specific functional random effect
$\bm{E}_{ijh}(t)$ with the respective proportion of variance explained. The
black solid line represents the mean trajectory to which a suitable multiple
($2\sqrt{\nu_{E\cdot}}$) of the FPC is added (red $+$) and subtracted (blue
$-$).

Figure 18 shows the estimated auto-and cross-covariance surface of components
$\bm{B}$ and $\bm{E}$ . It is evident that the cross-covariances are about as
large in magnitude as the auto-covariances for the dimension ACO. In
univariate analyses, however, the cross-covariance is completely ignored. Note
that when no covariates are included in the model, the estimated multiFAMM
contains one FPC for component $C$ (results not shown).

Figure 18: Auto- and cross-covariance surface of the subject-specific
functional random effect $\bm{B}_{i}(t)$ (left) and curve-specific random
effect $\bm{E}_{ijh}(t)$ (right). The bottom-left and top-right panels show
the auto-covariance and the top-left and bottom-right panels show the cross-
covariance function as a function of normalized time.

#### Analysis of the Estimated Effect Functions

Figure 19 presents the estimated covariate effects (black, top plot for
consonant order /s#sh/, bottom for /sh#s/). Overall, we find similar shapes
for the estimated effects on both dimensions. In the reference group, the
functional intercept $\bm{f}_{0}(t)$ shows signs of assimilation for consonant
/s/. The effect $\bm{f}_{1}(t)$ of covariate order, however, pushes the final
target syllable (in this case /sh/) towards its reference, all other things
being equal. The positive effect at the end of the observed interval then
pushes the initial target syllable /s/ towards the center, thus indicating a
more assimilated pronunciation. Shortly before, we find a negative effect on
the dimension EPG which might indicate that for a brief section, the
articulatory pattern of /s/ is indeed close to its reference but this does not
necessarily translate to the produced sound. Thus, similar to Cederbaum et al.
(2016) (red), we find that the final target syllable /sh/ is pulled towards
the reference while the initial target syllable /s/ seems to be less affected.
Given the similar shape on the dimension EPG, this supports their finding that
assimilation is asymmetric. Since the estimated effects are similar across
dimensions and similar to the univariate results, we refer to Cederbaum et al.
(2016) for interpretation of the other fixed effects.

Figure 19: Estimated covariate effects (black) and comparison to univariate
model (red). The corresponding $95\%$ confidence intervals are given by the
dashed line. Upper (from left to right): reference mean and covariate effects
of stress1, stress2, vowel (/s#sh/). Lower (from left to right): covariate
effect of order and interactions of order with stress1, stress2, vowel
(/sh#s/). Areas shaded in red indicate where the standard errors of the
univariate analysis are smaller.

#### Comparison to Univariate Analysis

Compared to separate univariate FAMMs, the multivariate model incorporates the
dependency between the dimensions, thus reducing the number of FPC bases in
the analysis (for a similar amount of variance explained, ten univariate FPCs
would be needed). The shaded areas in Figure 19 indicate where the standard
errors of the univariate analysis presented in Cederbaum et al. (2016) are
smaller than the corresponding standard errors of the multiFAMM. We find that
overall, the multivariate analysis seems to give smaller standard errors. In
order to compare the model fit between the univariate and the multivariate
analysis, we can use the urrMSE defined in Section 5. We then compare the
fitted values with the observed values for the consonant assimilation data.
The multivariate analysis yields urrMSE values of $0.970$ (ACO) and $0.966$
(EPG), whereas independent univariate FAMMs give values of $0.978$ (ACO) and
$0.961$ (EPG), respectively. Consequently, a univariate model analogously
specified to the ACO model in Cederbaum et al. (2016) gives a slightly better
model fit on the EPG data as measured by the urrMSE. However, on the ACO
dimension, the multivariate analysis is slightly preferred.

The added computational complexity of multivariate analyses is also negligible
in our case: Fitting a univariate FAMM as proposed by Cederbaum et al. (2016)
takes around 52 minutes on a 64-bit Linux platform and requires a considerable
amount of RAM memory (32 GB is sufficient). The multivariate analysis
maintains the requirements for internal memory, while the duration to fit the
multivariate model (109 minutes on the aforementioned platform) only slightly
increases compared to sequentially fitting two univariate models.

###  C.3 Results of the Model Using a Weighted Scalar Product

This section gives a short description of the considered multiFAMM when a
weighted scalar product is used. We only present the results of estimating the
eigenfunction basis as the effect on the estimated covariate effects is
negligible. Table 5 shows the contribution of each random process to the total
weighted variation and is structured similarly to Table 4. The number of
eigenfunctions is chosen by weighted sum of total variation (6) which results
in the same number of eigenfunctions for each random component as with the
unweighted approach. However, the total share of weighted variation explained
is slightly different (about $22\%$ and $65\%$ explained by the components
$\bm{B}$ and $\bm{E}$, respectively). Consequently, the proportion of
(weighted) variation assigned to the measurement error is smaller when
applying the weighted scalar product. We also find that the univariate norms
of the eigenfunctions are now more evenly distributed between the dimensions
ACO and EPG, with most of the norms bigger on the ACO dimension. We thus see
that weighting the scalar product shifts emphasis to the dimension ACO which
has a lower measurement error. With this emphasis on ACO, the proportion of
univariate variation explained on EPG falls now short of $95\%$ while ACO has
a proportion of univariate variation explained of about $98\%$. For this
model, too, we have specified the cut-off using the weighted sum of total
variation (6) and the model explains $96\%$ of variation in the data.

Table 5: Variance components included in the multiFAMM using a weighted scalar
product. First row: Estimates of eigenvalues, univariate error variances,
total variation. Second (third) row: Univariate norms of estimated FPCs on
dimension ACO (EPG). Fourth (fifth) row: Proportion of univariate variation
explained on dimension ACO (EPG) by eigenfunctions and error variance, total
univariate variation explained. Sixth row: Proportion of multivariate
variation explained by eigenfunctions and error variances, total multivariate
variation explained.

| $B_{1}$ | $B_{2}$ | $B_{3}$ | $E_{1}$ | $E_{2}$ | $E_{3}$ | $E_{4}$ | $E_{5}$ | $\sigma_{ACO}^{2}$ | $\sigma_{EPG}^{2}$ | Total
---|---|---|---|---|---|---|---|---|---|---|---
Variation | $1.937$ | $1.625$ | $0.493$ | $6.501$ | $2.168$ | $1.823$ | $0.694$ | $0.507$ | $0.004$ | $0.014$ | $18.476$
$||\psi^{(ACO)}||^{2}$ | $0.630$ | $0.743$ | $0.442$ | $0.574$ | $0.741$ | $0.404$ | $0.564$ | $0.455$ | $--$ | $--$ | $--$
$||\psi^{(EPG)}||^{2}$ | $0.370$ | $0.257$ | $0.558$ | $0.426$ | $0.259$ | $0.596$ | $0.436$ | $0.545$ | $--$ | $--$ | $--$
$\pi^{(ACO)}$ | $0.115$ | $0.114$ | $0.021$ | $0.352$ | $0.152$ | $0.070$ | $0.037$ | $0.022$ | $0.094$ | $--$ | $0.975$
$\pi^{(EPG)}$ | $0.091$ | $0.053$ | $0.035$ | $0.352$ | $0.071$ | $0.138$ | $0.038$ | $0.035$ | $--$ | $0.127$ | $0.941$
$\pi$ | $0.105$ | $0.088$ | $0.027$ | $0.352$ | $0.117$ | $0.099$ | $0.038$ | $0.027$ | $0.054$ | $0.054$ | $0.960$

Figure 20 shows the estimated multivariate FPCs of random component $\bm{B}$.
We find that the leading FPC again depicts an individual tendency to spend
more time pronouncing the final or the initial target consonant. For
$\bm{\psi}_{B2}(t)$, almost the same mode of variation is captured on the
dimension ACO as with the unweighted scalar product. On the dimension EPG,
however, we find the assimilation primarily in the final target syllable with
the initial target syllable rather unaffected (more or less time spent in the
pronunciation). The third FPC for the subject-specific random effect is
comparable to the scenario of the unweighted scalar product.

Figure 20: FPCs for the subject-specific functional random effect
$\bm{B}_{i}(t)$ with the respective proportion of weighted variance explained.
The black solid line represents the mean trajectory to which a suitable
multiple ($2\sqrt{\nu_{B\cdot}}$) of the FPC is added (red $+$) and subtracted
(blue $-$).

Figure 21 shows the estimated multivariate FPCs of random component $\bm{E}$.
The leading FPC is somewhat unchanged by the weighting of the scalar product.
While in the unweighted case the effect of $\bm{\psi}_{E2}$ was distributed
equally across the two consonants, we find a stronger effect on the initial
target consonant for dimension ACO and on the final target consonant for
dimension EPG in the weighted case. This seems to be compensated by the third
FPC, where in the weighted scenario the emphasis lies on differences in the
final target consonant on ACO and in the initial target consonant on EPG. The
remaining FPCs are comparable for both the weighted and unweighted scalar
product.

Figure 21: FPCs for the curve-specific functional random effect
$\bm{E}_{ijh}(t)$ with the respective proportion of weighted variance
explained. The black solid line represents the mean trajectory to which a
suitable multiple ($2\sqrt{\nu_{E\cdot}}$) of the FPC is added (red $+$) and
subtracted (blue $-$).

We conclude that the model based on the weighted scalar product gives similar
results. However, if some dimension were considerably more noisy than the
rest, weighting the scalar product might be advisable (see results of the
simulation).

##  Appendix D Detailed Analysis of the Simulation Study

We conduct an extensive simulation study in order to answer the following
questions: How does the performance of the proposed multiFAMM depend on
different model specifications? How does the model perform in different data
settings? And how does the multiFAMM compare to a univariate modeling
approach, where univariate regression models as proposed by Cederbaum et al.
(2016) are estimated on each dimension independently? Additionally, we
evaluate the estimation of the covariance and the fixed effects. We simulate
data based on the presented model fits (10) of the snooker training data (main
part) and (11) of the consonant assimilation data (Appendix C). Six different
settings of the data generating process are analysed, where for each data
setting, we additionally compare up to seven different model specifications.
One modeling scenario (data setting plus model specification) consists of
independent model fits based on 500 generated data sets. Table 6 provides an
overview of the analysed data settings and model specifications, which are
described in the following section. Note that we use data setting 1 and model
specification A, respectively, to generate benchmark estimates useful for
comparison to other settings and specifications.

Table 6: Summary of the different data settings and model specifications analysed in the simulation study. Each analysed combination of data setting and model specification comprises 500 model fits. Data Settings | Model Specification
---|---
1 | Consonant assimilation data | A | True model
| 2 |  | Strong heteroscedasticity | B | Truncation via total variation (TV)
| 3 |  | Sparse data | C | Truncation via univariate variation (UV)
| 4 |  | Uncentred scores | D | TV with alternate scalar product
| 5 |  | Weighted scalar product | E | UV with alternate scalar product
6 | Snooker training data | F | Scedastic misspecification
|  |  |  | U | Univariate models

###  D.1 Simulation Description

#### Data Settings

Table 7 provides detailed description of the data generating process for all
data settings. For each data setting, 500 data sets have been simulated. For
data setting $1$, we generate new observations based on the model fit (11)
(Appendix C) of the consonant assimilation data by randomly drawing the
evaluation points, the random scores, and the measurement errors. Note that we
center and decorrelate the random scores so that the empirical means and
covariances coincide with the theoretical counterparts. As in the original
data, each simulated data set contains the bivariate functional observations
of $i=1,..,9$ individuals, each repeating the $j=1,...,16$ different word
combinations five times.

Setting 2 (strong heteroscedasticity) is analogously to setting 1 but with
larger difference between the measurement error variances of the two
dimensions. Setting 2 mimics an application setting with (multimodal)
multivariate data, where some dimensions are much noisier than others, in
order to evaluate whether this imbalance interferes with the variance
decomposition of the multiFAMM. Setting 3 (sparse data) focuses on the
estimation quality for sparse functional data. Compared to setting 1, the
number of evaluation points is reduced to three to ten measurements per
dimension. In setting 4 (uncentred scores), the random scores are not
decorrelated or centred (otherwise identical to setting 1). Especially for
covariance components with few grouping levels (particularly $\bm{B}$), this
can result in a considerable departure from the modeling assumptions, which is
likely to occur in such settings. This setting explores the sensitivity of the
multiFAMM to violations of its assumptions. Setting 5 (weighted scalar
product) is identical to setting 1 but all model components are simulated
based on the estimated model using a weighted scalar product for the MFPCA
(see Appendix C.3). This data setting helps to understand the impact of
weights in the scalar product. For setting 6 (snooker training data), we
generate new trajectory data according to the model fit (10) of the snooker
training data. This allows us to evaluate a higher dimensional multiFAMM as
well as the estimation quality of nested random effects.

Table 7: Description of data settings.

Data Setting |  Description
---|---
1 |  Data based on consonant assimilation data (bivariate functions on $dim1,dim2$, each on $[0,1]$). $9$ individuals, each repeating $16$ crossed grouping levels $5$ times. Number of observed time points is an independent random draw for each dimension from a uniform discrete distribution of natural numbers in $[20,50]$. Observed time point is an independent random draw from uniform distribution on $[0,1]$. Mean consists of functional intercept and $7$ covariate effect functions as given in Appendix C.2. Covariates are given for each observation based on the crossed grouping level (word combination). Functional random effects $\bm{B},\bm{E}$ are based on estimated eigenfunctions of model in Appendix C.2 ($3$ and $5$ FPCs). Corresponding random scores are independent draws from $N(0,\hat{\nu}_{\cdot\cdot})$, then demeaned and decorrelated. Measurement errors are independent draws from $N(0,\hat{\sigma}^{2}_{\cdot})$ for each observed time point on the respective dimension.
2 |  Measurement errors are independent draws from $N(0,\hat{\sigma}^{2}_{1})$ for each observed time point on $dim1$ and independent draws from $N(0,16\cdot\hat{\sigma}^{2}_{1})$ on $dim2$. Rest as in 1.
3 |  Number of observed time points is a random draw for each dimension from a uniform discrete distribution of natural numbers in $[3,10]$. Rest as in 1.
4 |  Random scores are independent draws from $N(0,\hat{\nu}_{\cdot\cdot})$ (no centering or decorrelation). Rest as in 1.
5 |  All estimated components (eigenfunctions, eigenvalues, covariate effect functions, measurement error variances) stem from the model based on a weighted scalar product for the MFPCA as presented in Appendix C.3. Rest as in 1.
6 |  Data based on snooker training data (six-dimensional functions on $dim1,...,dim6$, each on $[0,1]$). $25$ individuals, each repeating $2$ nested grouping levels $6$ times. Number of observed points is an independent random draw from a uniform discrete distribution of natural numbers in $[10,50]$ per multivariate curve with observed time points identical over the dimensions. Observed time point is an independent random draw from uniform distribution on $[0,1]$. Mean consists of functional intercept and $4$ covariate effect functions as given in Appendix B. Covariates skill, group are independent random draws from Bernoulli distribution with probability $0.5$ per subject and session is given for each observation based on the nested grouping level. Functional random effects $\bm{B},\bm{C},\bm{E}$ are based on estimated eigenfunctions and eigenvalues of model in Appendix B ($6,5$ and $5$ FPCs). Measurement errors are independent draws from $N(0,\hat{\sigma}^{2})$ for each observed time point on each dimension.

#### Model Specifications

Table 8 provides a detailed description of the different modeling scenarios
used in the simulation study. We denote the most accurate approach of modeling
as specification A (true model). This standard scenario mirrors the data
generation so that there is no model misspecification. Most notably, we fix
the number of FPCs to the number used for generating the data in order to