# Flexible joint models for multivariate longitudinal and time-to-event data
using multivariate functional principal components

Alexander Volkmann  Chair of Statistics, School of Business and Economics,
Humboldt-Universität zu Berlin, Germany  Nikolaus Umlauf  Department of
Statistics, Faculty of Economics and Statistics, Universität Innsbruck,
Austria  Sonja Greven  Chair of Statistics, School of Business and Economics,
Humboldt-Universität zu Berlin, Germany

###### Abstract

The joint modeling of multiple longitudinal biomarkers together with a time-
to-event outcome is a challenging modeling task of continued scientific
interest. In particular, the computational complexity of high dimensional
(generalized) mixed effects models often restricts the flexibility of shared
parameter joint models, even when the subject-specific marker trajectories
follow highly nonlinear courses. We propose a parsimonious multivariate
functional principal components representation of the shared random effects.
This allows better scalability, as the dimension of the random effects does
not directly increase with the number of markers, only with the chosen number
of principal component basis functions used in the approximation of the random
effects. The functional principal component representation additionally allows
to estimate highly flexible subject-specific random trajectories without
parametric assumptions. The modeled trajectories can thus be distinctly
different for each biomarker. We build on the framework of flexible Bayesian
additive joint models implemented in the R-package ‘bamlss’, which also
supports estimation of nonlinear covariate effects via Bayesian P-splines. The
flexible yet parsimonious functional principal components basis used in the
estimation of the joint model is first estimated in a preliminary step. We
validate our approach in a simulation study and illustrate its advantages by
analyzing a study on primary biliary cholangitis.

##  1 Introduction

More and more medical studies are gathering information on their participants’
longitudinal profiles of different biomarkers together with the time to an
event of interest such as rejection of a transplanted organ [10], onset of a
disease [13], or death [33]. It is then of interest to model the longitudinal
profiles together with the time-to-event data to gain a detailed understanding
of the association of the different processes and to improve prediction for
new subjects.[8, 54, 38] Shared parameter models are a popular approach to
jointly modeling the longitudinal and time-to-event processes.[35] They
typically combine mixed effects models for the longitudinal data with a hazard
model for the time-to-event data through shared random effects. For multiple
longitudinal outcomes, however, the generalized mixed models and thus the
shared parameter models can become quite computationally expensive due to a
large number of parameters to estimate. This can restrict the flexibility with
which the longitudinal outcomes are modeled, even though the respective
longitudinal trajectories often exhibit strong non-linearities in
practice.[35, 19, 22] Recent research has thus focused on lowering the
computational costs of these multivariate joint models (MJMs) by adapting two-
stage methods (fitting the longitudinal submodel and using its output to fit
the time-to-event submodel)[28, 29], improving numerical efficiency[36], or
using approximate expectation-maximization algorithms[31, 32] or approximate
Bayesian methods[44, 47].

Here, we propose a further means of reducing the computational complexity of
MJMs, which can be combined with the aforementioned approaches. By
understanding the longitudinal profiles as functions over the follow-up time,
we introduce multivariate functional principal component analysis (MFPCA) as a
popular tool for dimension reduction in functional data analysis to the MJM
framework. This allows us to represent the shared random effects using
multivariate functional principal components (MFPCs), which has several
advantages: An MFPCs basis representation allows to model highly flexible
functional random effects while being parsimonious at the same time. These
basis functions efficiently capture correlated nonlinear trends over the
different longitudinal outcomes, thus admitting a flexible modeling approach
with relatively few parameters, especially compared to a (penalized) spline
approach for each longitudinal outcome. As the MFPCA decomposes the entire
multivariate random process, the number of MFPC based random effects needed
for a good model fit does not directly depend on the number of longitudinal
outcomes but only on the complexity of the covariance structure. Furthermore,
the MFPCs form an orthogonal basis with uncorrelated random weights (due to
the Karhunen-Loève theorem), which simplifies the generally unstructured
random effects covariance matrix of e.g. random spline basis coefficients to a
diagonal matrix. Finally, the MFPCs can be ordered by the amount of variation
they explain, which provides a handy criterion for further dimension reduction
by dropping relatively unimportant modes of variation from the estimation and
thus reducing the number of parameters in the model.

There is a wide range of work that has combined functional data analysis with
joint models of (multivariate) longitudinal and time-to-event outcomes. Some
of this work has focused on incorporating functional data such as images as
baseline covariates or as a repeatedly measured longitudinal functional
outcome.[59, 23, 21, 24, 60, 18, 61] In this context, the functional domain is
different from the follow-up time and functional principal component analysis
is used to achieve a low dimensional representation of the functional data.
Here, we propose to view repeated scalar measurements as functional data so
that the functional domain coincides with the follow-up time. This has also
been considered by other authors, particularly in the context of two-stage
models and dynamic predictions.[16, 55, 56, 25, 17, 27, 7] These approaches
use univariate functional principal component analysis (UFPCA) or MFPCA to
estimate principal components and scores for the longitudinal submodel and
plug the scores into the survival submodel. While two-stage approaches are
known to yield biased parameter estimates[46, 40], the previously mentioned
works show that this estimation bias does not necessarily negatively affect
the predictive performance. Other approaches simultaneously estimate the
principal components and scores within the univariate[57] or multivariate
joint model.[26, 22] For multiple longitudinal outcomes, however, these
proposed MJMs feature strong assumptions. Li et al. (2021)[26] use UFPCA to
model a common latent trajectory shared by all longitudinal outcomes, which
allows to represent only highly correlated markers. Li et al. (2022)[22] and
Zou et al. (2023)[62] extend this approach by further including outcome
specific latent trajectories, but assume them to be independent and
identically distributed (i.i.d.) across outcomes and subjects. In both
approaches, the dependence between longitudinal outcomes is explained by only
one shared latent factor, which might not be adequate for multiple biomarkers
from potentially different biological pathways.

While our proposed estimation approach also involves two steps, it is
important to note that this differs greatly from the previously mentioned two-
stage procedures. In the preliminary step, we estimate MFPCs via MFPCA. These
are then used as an empirical, flexible, and parsimonious basis to model the
multiple longitudinal outcomes within an MJM, thus incorporating the time-to-
event information together with the longitudinal modelling and avoiding
estimation bias induced by the marginal two-stage approach. Volkmann et al.
(2023)[51] show that this two step estimation process using estimated MFPCs
works well and provides the flexibility to fit complex multivariate functional
data. We are thus extending their work to the scenario where the multivariate
functional data are subject to informative missingness.

Our implementation combines the flexibility and parsimoniousness of using an
MFPC basis representation of functional random effects with the framework of
flexible additive joint models proposed by Köhler et al. (2017)[19]. This
framework incorporates the shared parameter models into Bayesian additive
models for location, scale, and shape[48], and provides versatile, powerful,
and easy to use modeling options for covariate effects. The implementation
allows to include for example smooth, linear, spatial, time-varying, or random
effects terms of one or more covariates in the longitudinal and time-to-event
submodels, using the well developed and documented architecture of the R
package mgcv[52] for generalized additive models[53]. This supports more (and
arguably easier to use) options for the specification of nonlinear covariate
effects than other available R software[15, 12, 30], in particular the popular
JMbayes2 package[43]. We adopt a derivative based Markov chain Monte Carlo
(MCMC) algorithm to estimate the full MJM, thus avoiding the different
proposed approximations such as via two-stage approaches as in JMbayes2[28].

We illustrate our proposed approach by analyzing a well known and publicly
available primary biliary cholangitis (PBC) data set[40]. It is based on a
study conducted by the Mayo Clinic from January 1974 to May 1984 on the
survival of patients with PBC (previously called primary biliary cirrhosis),
which is a rare chronic liver disease causing the inflammation and destruction
of small bile ducts.[33] This data set is still widely used in recent
literature [1, 31, 32, 34, 44, 47], mainly due to two reasons: it is one of
the last studies to allow the analysis of the natural history of PBC patients
treated only with supportive care [9] and it features an interesting data
structure. Several biochemical measurements are available for 312 subjects at
different time points after entering the study, with missing observations and
dropout due to censoring, liver transplantation, or death. Furthermore, these
longitudinal measurements have previously been found to follow non-linear
trajectories and to be associated with the patients’ survival.[1, 20, 44]
Using an MFPCs based representation of the random effects can lower the
computational burden while increasing the flexibility in the modeling of the
longitudinal trajectories, which can lead to a more accurate estimation of the
association between the longitudinal and the event processes.

This paper is structured as follows. Section 2 formally introduces the
flexible MJM as well as the MFPC representation of the random effects, while
Section 3 presents the proposed estimation of the MFPC basis and the flexible
MJM. In Section 4 we examine our approach in different simulation scenarios
and illustrate its applicability using the PBC data in Section 5. Section 6
closes with a discussion and outlook.

##  2 Methodology

In the following, we present a flexible joint model for multivariate
longitudinal outcomes and time-to-event data. We follow and extend the
notation of Köhler et al. (2017)[19] to emphasize the high flexibility of all
model parts as structured additive predictors, including modeling nonlinear
covariate effects using Bayesian P-splines. We then introduce the idea of
representing the shared random effects of the joint model using multivariate
functional principal components, which gives a parsimonious representation of
a possibly very complex covariance structure.

###  2.1 Multivariate joint model

Let $T_{i}\in\mathcal{I}$ denote the potentially right-censored event time for
subject $i=1,...,n$ with $\mathcal{I}=[0,T_{max}]$ and $T_{max}$ the maximal
follow-up time. We model the hazard $h_{i}(t)$ of an event for subject $i$ at
time $t\in\mathcal{I}$ as a function of $K$ longitudinally observed markers
and further covariates as

$h_{i}(t)=\text{exp}\\{\eta_{i}(t)\\}=\text{exp}\\{\eta_{\lambda
i}(t)+\eta_{\gamma
i}+\sum_{k=1}^{K}\eta_{\alpha_{k}i}(t)\cdot\eta_{\mu_{k}i}(t)\\},$ (1)

where $\eta_{i}(t)$ is a structured additive predictor for subject $i$ at time
$t$ given a set of baseline and time-dependent covariates $\bm{x}_{i}$. Note
that we suppress the dependency on $t$ for the set of covariates for better
readability. The additive predictor $\eta_{i}(t)$ can be further separated
into structured additive predictors $\eta_{\lambda i}(t)$ for the baseline
hazard, time-dependent covariates, and time-dependent covariate effects,
$\eta_{\gamma i}$ for the baseline covariates, and the possibly time-dependent
association modeled by $\eta_{\alpha_{k}i}(t)$ of the hazard with the $K$
longitudinal outcomes $\eta_{\mu_{k}i}(t)$, $k=1,...,K$. We consider the
“true” underlying longitudinal outcomes $\eta_{\mu_{k}i}(t)$ as structured
additive functions over $\mathcal{I}$ and assume them to be square-integrable,
$\eta_{\mu_{k}i}(t)\in L^{2}(\mathcal{I})$. These functions, however, are only
observed at subject and outcome specific time-points $t_{ijk},j=1,...,n_{ik}$
with corresponding longitudinal observations $y_{ijk}$ modeled as

$y_{ijk}=\eta_{\mu_{k}i}(t_{ijk})+\epsilon_{ijk},$ (2)

where $\epsilon_{ijk}$ represents independent Gaussian noise

$\epsilon_{ijk}\overset{\text{ind}}{\sim}N(0,\exp\\{\eta_{\sigma_{k}i}(t_{ijk})\\}^{2}),$
(3)

whose variance depends on the structured additive predictor
$\eta_{\sigma_{k}i}(t_{ijk})$.

The association predictors $\eta_{\alpha_{k}i}(t)$ link the multivariate
longitudinal submodel of equations (2) and (3) to the time-to-event submodel
(1). The resulting MJM connects the current outcome levels to the risk of
experiencing the event at time $t$. Note that we estimate separate association
predictors $\eta_{\alpha_{k}i}(t)$ for the different longitudinal predictors
$\eta_{\mu_{k}i}(t)$, allowing us to identify longitudinal outcomes important
for the time-to-event process.

The additive predictors

$\eta_{li}(t)=\sum_{h=1}^{H_{l}}f_{lh}(\bm{x}_{i},t),l\in\mathcal{L}:=\\{\lambda,\gamma,\alpha_{1},...\alpha_{K},\mu_{1},...,\mu_{K},\sigma_{1},...,\sigma_{K}\\}$

comprise parametric or nonparametric functions of $\bm{x}_{i}$ as well as
random effects, all potentially varying over time $t$. In particular, the
effects $f_{h}(\bm{x}_{i},t)$ depend on subsets of covariates in $\bm{x}_{i}$;
they can specify linear or smooth covariate effects of single covariates in
$\bm{x}_{i}$ or interaction effects of two or more such covariates. As an
illustration, a scalar, time-constant covariate $z_{i}\in\bm{x}_{i}$ might be
modeled as having a linear effect $z_{i}\cdot\beta$, a linear effect varying
over time $z_{i}\cdot f(t)$, a smooth effect $f(z_{i})$, or a smooth effect
over time $f(z_{i},t)$. For a detailed overview of different possible effect
specifications including spatial or covariate interaction effects, see Table 2
of Köhler et al. (2017)[19]. We use suitable basis representations such as
spline basis matrices with corresponding penalties to model the structured
additive predictors as explained in detail in Section 3.2. In particular, we
further separate the longitudinal predictors

$\displaystyle\eta_{\mu_{k}i}(t)=\sum_{h=1}^{H_{\mu_{k}}-1}f_{kh}(\bm{x}_{i},t)+b_{i}^{(k)}(t)$
(4)

into outcome specific (functional, i.e. potentially time-varying) fixed
effects $f_{kh}(\bm{x}_{i},t)\in L^{2}(\mathcal{I}),h=1,...,H_{\mu_{k}}-1$,
and subject and outcome specific functional random effects $b_{i}^{(k)}(t)$,
which account for the correlation of longitudinal outcome values within a
subject over time and across different markers.[45, 51] To model this
correlation structure, we assume that the multivariate functional random
effects $\bm{b}_{i}(t)=(b_{i}^{(1)}(t),...,b_{i}^{(K)}(t))^{\top}\in
L_{K}^{2}(\mathcal{I})$ with
$L_{K}^{2}(\mathcal{I})=L^{2}(\mathcal{I})\times...\times L^{2}(\mathcal{I})$
are independent copies of a smooth zero-mean Gaussian random process with
multivariate covariance kernel $\mathcal{K}$. The associated covariance
operator $\Gamma:L^{2}_{K}(\mathcal{I})\rightarrow L^{2}_{K}(\mathcal{I})$
given by
$(\Gamma\bm{g})(t)=\langle\langle\bm{g},\mathcal{K}(t,\cdot)\rangle\rangle$ is
based on the (potentially weighted) scalar product

$\displaystyle\langle\langle\bm{g},\bm{h}\rangle\rangle:=\sum_{k=1}^{K}w_{k}\int_{\mathcal{I}}g^{(k)}(t)h^{(k)}(t)dt,\quad\bm{g},\bm{h}\in
L_{K}^{2}(\mathcal{I})$ (5)

for given positive weights $w_{k},k=1,...,K$, inducing a Hilbert space
structure[14]. Setting weights unequal to one allows to let variation from
different longitudinal outcomes contribute differently to the multivariate
kernel, which will be discussed in more detail in Section 3.1. The
corresponding covariance surfaces
$\mathcal{K}^{(d,e)}(s,t)=\text{Cov}\left[b_{i}^{(d)}(s),b_{i}^{(e)}(t)\right],d,e=1,..,K$
and $s,t\in\mathcal{I}$, are auto-covariances for $d=e$ and cross-covariances
for $d\neq e$.

The auto-covariances describe the correlation structure of the random effects
within a longitudinal outcome. The covariance implied by a subject specific
constant shift and linear time trend such as in a parametric random intercept
– random slope mixed model would be contained as a special case. We do not
parametrically restrict the auto-covariance structure, however, and it can
thus also represent more flexible nonparametric forms such as in Köhler et al.
(2017)[19]. The cross-covariances capture the correlation structure between
the different longitudinal outcomes. As these are hard to realistically
specify a priori, especially for nonparametric random effects, we allow the
cross-covariances to nonparametrically vary with the pair of time-points
considered. Furthermore, we assume the $\bm{b}_{i}(t)$ to be independent of
the Gaussian noise process (3), which allows to separate the error variance
from the diagonal of the smooth covariance kernel $\mathcal{K}$.[58]

###  2.2 Functional principal component representation of random effects

Köhler et al. (2017)[19] achieve a flexible (auto-)covariance structure for
one longitudinal marker by representing the functional random effects using
B-spline basis functions with random effect weights. For multiple longitudinal
outcomes, however, such an approach can become prohibitive as the number of
parameters quickly increases: Using $D$ basis functions for each
$b_{i}^{(k)}(t)$ results in $n\cdot K\cdot D$ random basis weights with a
$(K\cdot D)\times(K\cdot D)$ unstructured covariance matrix $\bm{\Sigma}$
containing $(K\cdot D+1)K\cdot D/2$ unknown parameters. Given that the desired
flexibility is in practice often achieved by penalizing a large number $D$ of
basis coefficients, a resulting MJM suffers from the considerable number of
parameters to estimate. On the other hand, reducing the number $D$ or assuming
a diagonal covariance matrix, or both, are simplifying but unrealistic
assumptions which could impair the estimation of the association with the
event hazard.

Instead, we propose to use a functional principal components based
representation of the multivariate functional random effects. We apply the
multivariate Karhunen-Loève (KL) theorem[14], which allows to represent
multivariate functional data in an infinite basis expansion using random
scores and the eigenfunctions of the corresponding covariance operator. For
the multivariate functional random effects we get

$\displaystyle\bm{b}_{i}(t)=\sum_{m=1}^{\infty}\rho_{im}\bm{\psi}_{m}(t)\approx\sum_{m=1}^{M}\rho_{im}\bm{\psi}_{m}(t)$
(6)

with random scores $\rho_{im}$ and the eigenfunctions
$\bm{\psi}_{m}(t)=(\psi_{m}^{(1)}(t),...,\psi_{m}^{(K)}(t))^{\top}\in
L_{K}^{2}(\mathcal{I})$ of the covariance kernel $\mathcal{K}$[14, 51], i.e.
due to Mercer’s theorem
$\mathcal{K}(s,t)=\sum_{m=1}^{\infty}\nu_{m}\bm{\psi}_{m}(s)\bm{\psi}_{m}(t)^{\top}$
with ordered eigenvalues $\nu_{1}\geq\nu_{2}\geq\dots\geq 0$. Note that the
eigenfunctions $\bm{\psi}_{m},m=1,...,\infty$ are orthonormal to each other
with respect to the scalar product (5). The scores $\rho_{im}\sim(0,\nu_{m})$
are uncorrelated across $m$ and, in particular, for the Gaussian process
$\bm{b}_{i}(t)$ are i.i.d. $\rho_{im}\sim N(0,\nu_{m})$. In practice, the
infinite expansion in (6) is truncated at a finite level using the leading $M$
eigenfunctions, where the truncation order $M$ controls the approximation
accuracy. This truncated multivariate KL representation has several benefits:
Interpretability, dimension reduction, and the simplification of the scalar
random effects covariance matrix to a diagonal matrix.

The leading eigenfunctions correspond to the main modes of variation in the
data due to decreasing corresponding eigenvalues and can be interpreted as
outcome specific characteristic trajectories corresponding to correlated
patterns that are simultaneously present across all longitudinal outcomes. The
direction and strength of expression of these characteristic trajectories for
each subject depend on their individual score $\rho_{im}$. The eigenvalues
$\nu_{m}$ corresponding to the eigenfunctions allow to quantify their
contributions to the overall variation in the data given that
$E(\langle\langle\bm{b}_{i},\bm{b}_{i}\rangle\rangle)=\sum_{k=1}^{K}w_{k}\int_{\mathcal{I}}\text{Var}(b_{i}^{(k)}(t))dt=\sum_{m=1}^{\infty}\nu_{m}$.
This can be used to determine the quality of the approximation in (6), e.g. by
chosing $M$ such that a pre-specified large percentage of variation is
retained. We estimate the eigenfunctions as MFPCs in a preliminary step
described in section 3.1, leaving only the random scores $\rho_{im}$ to be
estimated as random effects in the MJM. This considerably reduces the number
of parameters from $n\cdot K\cdot D$ to $n\cdot M$: The scores act as weights
for multivariate functions and are thus not outcome-specific. In addition, $M$
only depends on the complexity of the covariance kernel $\mathcal{K}$, which
is not necessarily connected to the number of longitudinal outcomes $K$, such
that typically in practice $M\ll K\cdot D$. Using the orthonormal MFPCs as
basis functions for the Gaussian random process also leads to i.i.d.
$\bm{\rho}_{i}=(\rho_{i1},...,\rho_{iM})^{\top}\sim
N_{M}(\bm{0},diag(\nu_{1},...,\nu_{M})),i=1,..,n$, where the mean of the
$M$-dimensional multivariate normal distribution is the vector of zeros
$\bm{0}$ and the covariance matrix is diagonal with eigenvalues $\nu_{m}$ as
elements. This represents a further simplification, reducing the covariance
parameters to $M\ll(K\cdot D+1)K\cdot D/2$ independent of $K$, and implies
computational advantages compared to an unstructured covariance matrix
$\bm{\Sigma}$ of a B-spline approach as described above.

The eigenfunction based representation of the random effects gives the best
(in terms of integrated squared error [39]) approximation to the longitudinal
outcomes $\bm{b}_{i}(t)$ for a given number of basis functions and allows a
parsimonious model formulation. While in practice we have to use estimates of
the basis functions $\bm{\psi}_{m}(t)$ in the MJM, we are less interested in
the MFPCs themselves, although they provide additional interpretations of main
directions of variation in marker trajectories. Our simulations show that even
imperfectly estimated versions still serve as a good empirical basis for the
longitudinal trajectories and give the MJM the flexibility needed for a good
model fit, yielding better results than more restrictive parametric or spline-
based approaches in settings with nonlinear trends over time.

##  3 Estimation

We use a derivative-based MCMC algorithm adapted from the one for univariate
joint models proposed by Köhler et al. (2017)[19] to estimate the posterior
distribution of the parameter vector $\bm{\theta}$ containing all model
parameters of interest in a Bayesian framework. Note that the MFPC basis is
estimated directly from the data and is thus not part of $\bm{\theta}$. This
might be interpreted as an empirical Bayes approach; however, we prefer to
view the MFPCs simply as a beneficial empirical basis for the random effects
(comparable to e.g. pre-defined B-spline basis functions). The following
sections outline the estimation of the MFPC basis and the Bayesian estimation
of $\bm{\theta}$.

###  3.1 Estimation of the MFPC basis

We estimate the MFPC basis following the MFPCA proposed by Happ and Greven
(2018)[14], which estimates the multivariate covariance kernel $\mathcal{K}$
from separate UFPCAs based on the different longitudinal outcomes of each
subject. This approach is very flexible as it allows to specify the UFPCA for
each longitudinal outcome individually and relies only on the covariance
between the estimated univariate scores of the subjects.

In a first step, we estimate a mean trajectory $\eta_{\mu_{k}}^{*}(t_{ijk})$
for each longitudinal outcome using the outcome specific fixed effects
$\sum_{h}^{H_{\mu_{k}}-1}f_{kh}(\bm{x}_{i},t)$ specified for the MJM under a
working independence assumption for the longitudinal observations ignoring
censoring. We then estimate the auto-covariance of each longitudinal outcome
$y_{ijk}(t)\in L^{2}(\mathcal{I})$ as a bivariate surface using fast symmetric
additive covariance smoothing based on the crossproducts of the centered
longitudinal observations
$y_{ijk}^{*}=y_{ijk}-\eta_{\mu_{k}}^{*}(t_{ijk})$.[5] A truncated
eigendecomposition of the covariance surface then yields estimates for the
first $M_{k}$ univariate eigenfunctions $\phi_{km_{k}}\in L^{2}(\mathcal{I})$
and eigenvalues corresponding to the variances of the scores $\xi_{km_{k}i}$
in the univariate KL expansion
$y_{ijk}(t)\approx\eta_{\mu_{k}}^{*}(t)+\sum_{m_{k}=1}^{M_{k}}\xi_{km_{k}i}\phi_{km_{k}}(t)$.
Given the estimated principal components $\hat{\phi}_{km_{k}}$, the univariate
scores $\xi_{km_{k}i}$ are predicted using a conditional expectation
approach.[58] Other approaches of estimating the UFPCA are also possible (cf.
Section 6).

For a Hilbert space induced by (5), Happ and Greven (2018)[14] show that the
multivariate eigenfunctions are given as linear combinations of their
univariate counterparts and derive the explicit formula for this relationship
as
$\psi^{(k)}_{m}(t)=(w_{k})^{-1/2}\sum_{m_{k}=1}^{M_{k}}c_{km_{k}}\phi_{km_{k}}(t)$,
with linear combination weights $c_{km_{k}}$ depending on (the eigenvectors
of) the covariance matrix of the univariate scores
$\xi_{km_{k}i},k=1,...,K$.[14] The MFPCA thus uses the estimated covariance
structure of the estimated univariate scores of all longitudinal outcomes to
estimate the MFPCs and consequently the multivariate covariance kernel
$\mathcal{K}$. Regarding the choice of weights in (5), setting
$w_{1}=...=w_{K}=1$ can be a reasonable starting point for the analysis in
many cases. When the range or variation of the longitudinal data differ
considerably between the longitudinal outcomes, we propose to use weights
inversely proportional to the integrated univariate variance as given by the
sum of univariate eigenvalues from the UFPCAs.[14] This approach is equivalent
to standardization of the different longitudinal markers to make them
comparable, as commonly also done in non-functional multivariate principal
component analysis, while avoiding rescaling of the data. The resulting MFPCs
are then often more balanced across the longitudinal outcomes, avoiding MFPCs
which are dominated by one longitudinal outcome essentially capturing
variation of that outcome alone. Refraining from the rescaling of longitudinal
outcomes has the additional advantage that the association structure of the
MJM via the current value of the longitudinal outcomes can be interpreted much
more intuitively.

As the MFPC basis is an important part of the subsequent MJM, its estimation
should be closely monitored in practice. It is important that the UFPCAs are
flexible enough to capture all the relevant modes of variation in the data. In
particular, the univariate truncation order $M_{k}$ should be high and a
sufficient number of basis functions for the estimation of the univariate
covariance surface should be specified, as variation lost by oversmoothing on
the univariate level is lost for all later stages of the analysis. In
addition, the penalization chosen by the fast symmetric additive covariance
smoothing directly influences the smoothness of the MFPCs and consequently
also the final fit of the MJM. Visually inspecting the estimated MFPC basis
and consulting a subject matter expert about their beliefs on systematic
variation in the data at this stage of the analysis can be helpful.

Note that the estimation of the MFPC basis uses the estimated covariance
between the estimated univariate scores. As univariate scores are estimated
with large uncertainty if only very few longitudinal observations are
available, we recommend a trimmed approach where we use all longitudinal
observations of subjects with at least one measurement after 10% of the time
interval $\mathcal{I}$ on all longitudinal outcomes. In our experience, this
approach avoids extremely outlying estimates of univariate scores in the UFPCA
for subjects with very early censoring, which via the estimated covariance
between the univariate scores can have a detrimental effect on the estimated
MFPC basis. Alternatively, outlying scores can be manually identified and
removed from the MFPCA, or robust alternatives to trimming such as weighting
would also be possible.

Given a reasonable estimate of the MFPC basis, the multivariate eigenvalues
$\nu_{m}$ of the MFPCA can be used to choose the truncation order $M$ of the
multivariate functional random effects representation (6). The maximal number
$M^{*}=\sum_{k=1}^{K}M_{k}$ of MFPCs corresponds to the total number of
univariate functional principal components (UFPCs) in the MFPCA. Reducing
$M<M^{*}$ allows further dimension reduction for the parameter vector
$\bm{\theta}$ with relatively little loss in the estimation quality for high
$M$ as shown in Section 4. A common criterion is to use a threshold for the
ratio of cumulative explained variance, e.g. $99\%$, and discard all trailing
MFPCs. Given the MFPC basis, choosing $M$ thus reduces to a model selection
problem and can also be solved using standard approaches such as comparing
information criteria for different model fits.

###  3.2 Basis representation of effects

The proposed MJM is very flexible with respect to the readily available types
of covariate effects for all structured additive predictors
$\eta_{li}(t)=\sum_{h=1}^{H_{l}}f_{lh}(\bm{x}_{i},t),l\in\mathcal{L}$. The
Bayesian additive mixed models framework allows a wide range of effects such
as linear, smooth, or random effects terms by using corresponding basis
function expansions for $f_{lh}(\bm{x}_{i},t)$ combined with associated
priors.[48] All of these flexible effects, including the functional principal
component based random effects, can be represented in matrix notation as

$\displaystyle\bm{f}_{lh}=\bm{X}_{lh}\bm{\beta}_{lh}$ (7)

with vector $\bm{f}_{lh}$ of stacked function evaluations over all subjects
for the survival submodel (vector length $n$) or over all subjects, markers
and time points for the longitudinal submodel (vector length
$N=\sum_{i}\sum_{k}n_{ik}$), design matrix $\bm{X}_{lh}$, and coefficient
vector $\bm{\beta}_{lh}$. This representation includes but is not limited to
non-linear covariate effects, which can be represented as such a matrix
product with the columns of $\bm{X}_{lh}$ containing e.g. evaluations of
spline basis functions. For the representation of further model terms, see
Umlauf et al. (2018)[48] and Köhler et al. (2017)[19]. Note that standard sum-
to-zero constraints are used to center the design matrix $\bm{X}_{lh}$ of non-
linear functions[45, 53] in order to ensure model identifiability in a model
also including an intercept.

In particular, the functional principal components based random effects
$b^{(k)}_{i}(t)$ for one longitudinal outcome in (4) can be expressed as

$\displaystyle\bm{b}^{(k)}=\bm{X}_{\mu_{k}H_{\mu_{k}}}\bm{\beta}_{\mu_{k}H_{\mu_{k}}}=blockdiag(\bm{\Psi}_{1}^{(k)},...,\bm{\Psi}_{n}^{(k)})\bm{\rho},$
(8)

where the $N_{k}=\sum_{i}n_{ik}$-vector $\bm{b}^{(k)}$ contains the
evaluations of the subject and outcome specific functional random effect at
the observed time points for outcome $k$ and all objects $i$ and
$\bm{\rho}=(\bm{\rho}_{1}^{\top},...,\bm{\rho}_{n}^{\top})^{\top}$ is the
vector containing $M$ random scores $\bm{\rho}_{i}$ per subject (corresponding
to the coefficient vector $\bm{\beta}_{\mu_{k}H_{\mu_{k}}}$). The design
matrix $\bm{X}_{\mu_{k}H_{\mu_{k}}}$ has a blockdiagonal structure with blocks
$\bm{\Psi}_{i}^{(k)}$ per $i$, which is an $n_{ik}\times M$ matrix containing
evaluations of the $k$-th component of the $M$ MFPCs at time points $t_{ijk}$,
$j=1,\dots,n_{ik}$. Note that $\bm{\rho}$ is not specific to outcome $k$ and
its dimension thus does not scale with $K$ due to our dimension reduction
approach. Consequently, stacking the longitudinal outcomes to
$\bm{b}=(\bm{b}^{(1)\top},...,\bm{b}^{(K)\top})^{\top}$ simply corresponds to
stacking the $K$ blockdiagonal design matrices to a full multivariate
longitudinal design matrix $\bm{X}_{\mu
H_{\mu}}=(\bm{X}_{\mu_{1}H_{\mu_{1}}}^{\top},...,\bm{X}_{\mu_{K}H_{\mu_{K}}}^{\top})^{\top}$
with a single coefficient vector $\bm{\rho}$. We point out that using a
standard (spline) random effects basis for $\bm{b}$ would not be as
parsimonious. While for each $k$ the design matrix
$\bm{X}_{\mu_{k}H_{\mu_{k}}}$ would have a similar blockdiagonal structure as
in (8) with blocks containing the evaluations of the random effect (spline)
basis, the random coefficient vector $\bm{\beta}_{\mu_{k}H_{\mu_{k}}}$ would
be different for each $k$. As a result, combining the longitudinal outcomes
would be equivalent to constructing a full multivariate longitudinal design
matrix $\bm{X}_{\mu
H_{\mu}}=blockdiag(\bm{X}_{\mu_{1}H_{\mu_{1}}},...,\bm{X}_{\mu_{K}H_{\mu_{K}}})$.

###  3.3 Priors

In order to achieve the penalization of flexible effects in a Bayesian
setting, appropriate priors for $\bm{\beta}_{lh}$ are specified. In general,
we assume multivariate normal priors

$\displaystyle
p(\bm{\beta}_{lh}\mid\tau^{2}_{lh})\propto(\tau_{lh}^{2})^{-\frac{\text{rank}(\bm{K}_{lh})}{2}}\exp\left(-\frac{1}{2\tau^{2}_{lh}}\bm{\beta}_{lh}^{\top}\bm{K}_{lh}\bm{\beta}_{lh}\right)$
(9)

with variance parameter $\tau_{lh}^{2}$, $\text{rank}(\bm{A})$ denoting the
rank of matrix $\bm{A}$, and $\bm{K}_{lh}$ the chosen precision matrix. The
variance parameter takes the role of an inverse smoothing parameter
controlling the trade-off between flexibility and smoothness of the flexible
effect, while the precision matrix determines the type of penalization and
contains e.g. difference penalties for Bayesian P-splines.[48] We assume
independent inverse Gamma hyperpriors for the variance parameters
$\tau_{lh}^{2}\sim IG(0.001,0.001)$, which results in inverse Gamma full
conditionals but other priors, such as half-Cauchy, are also possible.
Furthermore, priors for anisotropic smooths containing multiple variance
parameters can be specified as given in Köhler et al. (2017) [19]. For linear
or parametric time-constant terms, we approximate $\bm{K}_{lh}=\bm{0}$ by
using $\bm{\beta}_{lh}\sim N(\bm{0},1000^{2}\bm{I})$ with identity matrix
$\bm{I}$.

Following the truncated KL representation (6), we can assume independent
random scores over the subjects $i$ and the MFPCs $m$ so that we can separate
$\bm{\rho}=(\bm{\rho}_{(1)},...,\bm{\rho}_{(M)})^{\top}$ with
$\bm{\rho}_{(m)}=(\rho_{1m},...,\rho_{nm})^{\top}$ and use separate priors for
each $\bm{\rho}_{(m)}$. In particular, we assume the multivariate normal
priors (9) for each $m$, with precision matrix $\bm{K}_{(m)}=\bm{I}$ and
different variance parameters $\tau_{(m)}^{2}$. Note that the variances
$\tau_{(m)}^{2}$ correspond to the eigenvalues $\nu_{m}$ in the MFPCA. We
emphasize viewing the estimated MFPCs as an empirical basis by using the above
introduced weak Gamma priors for the $\tau_{(m)}^{2}$ instead of using the
estimated eigenvalues $\nu_{m}$ in the priors. With our specification, the
$\rho_{im}$ take the role of random effects for the basis functions
$\bm{\psi}_{m}$ \- defined simultaneously for all $k$ on the entire follow-up
time interval $\mathcal{I}$ \- as a flexible data-driven basis, replacing
parametric basis functions such as e.g. a constant and a linear function per
$k$ in a random intercept – random slope model.

###  3.4 Likelihood

We use the standard assumption of shared random effects models that given the
parameter vector $\bm{\theta}$, the longitudinal outcomes
$\bm{y}=(\bm{y}_{1}^{\top},...,\bm{y}_{K}^{\top})^{\top}$ and the survival
outcome $[\bm{T},\bm{\delta}]$ are conditionally independent. Here,
$\bm{y}_{k}=(\bm{y}_{1k}^{\top},...,\bm{y}_{nk}^{\top})^{\top},\bm{y}_{ik}=(y_{i1k},...,y_{in_{ik}k})^{\top}$,
$\bm{T}=(T_{1},...,T_{n})^{\top}$ the follow-up (event or censoring) times,
and the event indicators $\bm{\delta}=(\delta_{1},...,\delta_{n})^{\top}$ are
1 if the subject experiences the event and 0 for censoring. This allows to
factorize the likelihood of the MJM as the product of the Likelihood
$L^{surv}$ of survival submodel (1) and the Likelihood $L^{long}$ of the
longitudinal submodel (2)-(3). The log-likelihood of the MJM can then be
formulated as

$\displaystyle\ell(\bm{\theta}\mid\bm{T},\bm{\delta},\bm{y})=\ell^{surv}(\bm{\theta}\mid\bm{T},\bm{\delta})+\ell^{long}(\bm{\theta}\mid\bm{y})$
(10)

with the log-likelihoods $\ell^{surv}$ and $\ell^{long}$ of the survival and
longitudinal submodels, respectively.

For the survival submodel, the log-likelihood can be expressed as

$\displaystyle\ell^{surv}(\bm{\theta}\mid\bm{T},\bm{\delta})=\bm{\delta}^{\top}\bm{\eta}(\bm{T})-\bm{1}_{n}^{\top}\bm{\Lambda}(\bm{T}),$

where $\bm{\eta}(\bm{T})=(\eta_{1}(T_{1}),...,\eta_{n}(T_{n}))^{\top}$,
$\bm{1}_{u}$ is the $u$-vector of ones, and
$\bm{\Lambda}(\bm{T})=(\Lambda_{1}(T_{1}),...,\Lambda_{n}(T_{n}))^{\top}$ is
the vector of cumulative hazard rates
$\Lambda_{i}(T_{i})=\int_{0}^{T_{i}}\eta_{i}(s)ds$.

The log-likelihood of the longitudinal submodel is similar to a multivariate
functional additive mixed model [51] with the longitudinal outcomes stacked in
the outcome vector $\bm{y}$ and

$\displaystyle\ell^{long}(\bm{\theta}\mid\bm{y})$
$\displaystyle=-\frac{N}{2}\log(2\pi)-\bm{1}_{N}^{\top}\bm{\eta}_{\sigma}-\frac{1}{2}(\bm{y}-\bm{\eta}_{\mu})^{\top}\bm{R}^{-1}(\bm{y}-\bm{\eta}_{\mu})$

with similarly stacked vectors
$\bm{\eta}_{\sigma}=(\bm{\eta}_{\sigma_{1}}^{\top},...,\bm{\eta}_{\sigma_{K}}^{\top})^{\top}$,
$\bm{\eta}_{\sigma_{k}}=(\bm{\eta}_{\sigma_{k}1}^{\top},...,\bm{\eta}_{\sigma_{k}n}^{\top})^{\top}$,
$\bm{\eta}_{\sigma_{k}i}=(\eta_{\sigma_{k}i}(t_{i1k}),...,\eta_{\sigma_{k}i}(t_{in_{ik}k}))^{\top}$
and
$\bm{\eta}_{\mu}=(\bm{\eta}_{\mu_{1}}^{\top},...,\bm{\eta}_{\mu_{K}}^{\top})^{\top}$,
$\bm{\eta}_{\mu_{k}}=(\bm{\eta}_{\mu_{k}1}^{\top},...,\bm{\eta}_{\mu_{k}n}^{\top})^{\top}$,$\bm{\eta}_{\mu_{k}i}=(\eta_{\mu_{k}i}(t_{i1k}),...,\eta_{\mu_{k}i}(t_{in_{ik}k}))^{\top}$,
and diagonal matrix $\bm{R}=diag(\exp(\bm{\eta}_{\sigma})^{2})$.

###  3.5 Posterior mean estimation

The posterior of the full parameter vector $\bm{\theta}$ factorizes to

$\displaystyle p(\bm{\theta}\mid\mathbf{T},\bm{\delta},\mathbf{y})$
$\displaystyle\propto L^{surv}(\bm{\theta}\mid\mathbf{T},\bm{\delta})\cdot
L^{long}(\bm{\theta}\mid\bm{y})\cdot\prod_{l\in\mathcal{L}}\prod_{h=1}^{\tilde{H}_{l}}(p(\bm{\beta}_{lh}\mid\bm{\tau}_{lh}^{2})p(\bm{\tau}_{lh}^{2}))\cdot\prod_{m=1}^{M}(p(\bm{\rho}_{(m)}\mid\bm{\tau}_{(m)}^{2})p(\bm{\tau}_{(m)}^{2})),$

where $\tilde{H}_{l}=H_{l}-1$ for the longitudinal outcomes
$l\in\\{\mu_{1},...,\mu_{K}\\}$ and $\tilde{H}_{l}=H_{l}$ otherwise. We obtain
starting values for the derivative-based MCMC algorithm by approximating the
mode of the posterior using a blockwise Newton-Raphson procedure.[19] In
contrast to Köhler et al. (2017), however, we use the approximate mode only as
a starting value for mean estimation, thus prioritizing the fast convergence
of the Newton-Raphson algorithm. The step length is set to 1 except for the
survival predictors $l\in\\{\lambda,\gamma\\}$, where it is set to 0.1, and
the variance parameters are initialized by the bamlss algorithm so that the
empirical degrees of freedom are comparable across model terms.[49, 50]

Then, the MCMC algorithm uses the blockwise score vectors
$\bm{s}(\bm{\beta}_{lh})$ and Hessian $\bm{\mathcal{H}}(\bm{\beta}_{lh})$ of
the log-posterior given in Appendix A to approximate the full conditionals
$p(\bm{\beta}_{lh}\mid\cdot)$.[48] These approximations are based on a second-
order Taylor expansion of the log-posterior centered at the last state
$\bm{\beta}_{lh}^{[u]}$ for MCMC iteration $u$ and result in a proposal
density proportional to a multivariate normal distribution
$N(\bm{\mu}_{lh}^{[u]},\bm{\Sigma}_{lh}^{[u]})$ with mean
$\bm{\mu}_{lh}^{[u]}=\bm{\beta}_{lh}^{[u]}-\bm{\mathcal{H}}(\bm{\beta}_{lh}^{[u]})^{-1}\bm{s}(\bm{\beta}_{lh}^{[u]})$
and precision matrix
$(\bm{\Sigma}_{lh}^{[u]})^{-1}=-\bm{\mathcal{H}}(\bm{\beta}_{lh}^{[u]})$.
While the calculation of the derivatives increases computational cost, it has
the advantage that this MCMC algorithm approximates a Gibbs sampler, resulting
in high acceptance rates and good mixing. Gibbs sampling is used for variance
parameters $\tau_{lh}^{2}$ with conjugate inverse Gamma hyperpriors and slice
sampling otherwise.

###  3.6 Implementation

The implementation of the proposed approach builds on the flexibility provided
by the R package bamlss[49, 50]. We supply the package MJMbamlss on Github
(https://github.com/alexvolkmann/MJMbamlss) containing a new family-class for
the presented MJM as well as some convenience functions for estimating the
MFPC basis. The contained vignettes illustrate how to combine the two packages
for model estimation. Similar to the package JMbayes2[42] we find that
standardizing the survival design matrices within the estimation process
improves mixing and convergence of the MCMC algorithm as well as posterior
mode estimation, and is implemented by default. The integrals in the
cumulative hazard rates $\Lambda_{i}(T_{i})$ and in the corresponding
derivatives are of non-standard form and we use seven point Gaussian
quadrature to numerically approximate them.[37]

##  4 Simulations

We show the good performance of our proposed approach in two simulation
scenarios. Scenario I aims at comparing our approach to the JMbayes2[42]
package in a parametric setting that is favorable to JMbayes2, i.e. in a
setting with a simple random intercept – random slope covariance structure,
here for six longitudinal outcomes. Our approach, estimating the random
effects structure nonparametrically, should be able to also recover this truly
linear structure, even though a higher variance compared to correctly
specifying the parametric model in JMbayes2 is to be expected due to the
increased model flexibility. This example also illustrates the reduction of
parameters to estimate in our MFPC based representation, while keeping a
comparable estimation performance. Scenario II by contrast focuses on a
complex covariance structure, here between two nonlinear longitudinal markers,
and highlights the powerful and parsimonious modeling possibilities of MFPCs
for highly flexible longitudinal trajectories. In the following, the data
generating processes of the two simulation scenarios are described, the
general setup of the simulation is introduced, and the results are presented.
Full details of the data generation, such as the used parameter values, as
well as additional results are provided in Appendix B.

###  4.1 Scenario I

We generate a sequence of six different longitudinal trajectories for $n=150$
subjects at an equidistant grid of 101 time points $\mathcal{P}\in[0,1]$. The
longitudinal trajectories follow a linear predictor with a fixed linear time
trend and linear effect of a binary covariate $x_{gi}$, sampled from a
discrete uniform distribution, as well as a random intercept and random slope
over time. The fixed effects are chosen identical over the different
longitudinal outcomes, but the random effects differ considerably in
magnitude. We draw the random effects from a 12-dimensional multivariate
normal distribution with mean zero and given unstructured covariance matrix.
While the correlation between random intercept and slope within a longitudinal
outcome is the same for all outcomes, the scale of the variances and thus
covariances increases with $k$. At the same time, the correlation between
outcomes $k$ and $l$ is specified as a decreasing function of their
difference, i.e. “neighbouring” outcomes have a higher correlation, and equal
to zero for the furthest pairs (12 of the 78 covariance parameters are zero,
corresponding to independence of random effects). We specify the linear
predictors $\eta_{\gamma}$ as containing the binary covariate $x_{gi}$, the
baseline hazard $\eta_{\lambda}$ as following a Weibull distribution, and
$\eta_{\alpha_{1}},...,\eta_{\alpha_{6}}$ as constants with decreasing values
from $1.5$ to $-1.5$. The survival times are then generated by calculating the
hazard rate and deriving the corresponding survival probabilities.[3, 6]
Random censoring is introduced by drawing censoring times from a uniform
distribution on $[0,1.75]$. The subject specific follow-up time $T_{i}$ is the
minimum of the survival time, the censoring time and 1. We introduce subject
and outcome specific observation times by sampling $25\%$ of the remaining
subject specific subset of $\mathcal{P}$ with a maximum of $15$ longitudinal
observations per outcome. The longitudinal observations are then generated by
adding the Gaussian noise $\epsilon_{ijk}$ with variance constant across $k$
to the linear predictor evaluated at the observation times. In the generated
data sets, all subjects have baseline measurements of the longitudinal
outcomes at time point 0 so that the subjects have a minimum of six (one per
outcome) and up to $90$ (mean $63$) longitudinal observations across all
outcomes. The data sets show a mean event rate of $43\%$, with mean follow-up
time of $0.52$.

###  4.2 Scenario II

For simulation scenario II, we follow the same data generation structure to
obtain data sets with $n=300$ subjects and two longitudinal outcomes based on
different coefficient specifications for the additive predictors compared to
scenario I and a random censoring uniform distribution on $[0,3]$. Most
notably, the random intercept and slope components in the longitudinal
additive predictors are replaced by a more complex covariance structure. We
include multivariate functional random effects drawn from a multivariate
Gaussian process with a covariance kernel based on six multivariate
eigenfunctions. To do this, we follow the results for finite KL decompositions
of Happ and Greven (2018)[14], and define univariate non-orthogonal basis
functions and a covariance matrix for the coefficients of the corresponding
basis expansion to construct the multivariate covariance kernel. We specify
three basis functions per longitudinal outcome as splines with distinct
characteristics at the beginning, middle, or end of the interval $[0,1]$: For
outcome $1$ the basis functions correspond to level shifts, whereas for
outcome $2$ the basis functions describe short-term peaks (see Appendix B
Figure 3). The $6\times 6$ covariance matrix of the basis coefficients is
generated from an eigenvalue decomposition with specified eigenvalues and a
random orthonormal matrix. This gives a non-standard multivariate covariance
kernel with highly nonlinear eigenfunctions (see Appendix B Figure 4). To
generate multivariate functional random effects from this Gaussian process, we
draw i.i.d. samples $\bm{\rho}_{i}=(\rho_{i1},...,\rho_{i6})^{\top}\sim
N_{6}(\bm{0},diag(\nu_{1},...,\nu_{6})),i=1,...,300$ with corresponding
multivariate eigenvalues $\nu_{1},...,\nu_{6}$. The data sets in scenario II
have on average 24 longitudinal observations (minimum 2 and maximum 30) and a
mean event rate of $57\%$ with mean follow-up time 0.6.

###  4.3 Simulation Setup

For both scenarios we simulate 200 data sets to which we apply four different
model estimation approaches. First, we estimate a model using the true
underlying covariance kernel via the corresponding true MFPC basis in the
bamlss framework to have an oracle benchmark model (denoted as TRUE).
Comparing the benchmark to a model where the covariance structure is estimated
from the data as described in Section 3.1 (EST) allows us to quantify the loss
in estimation accuracy due to the MFPC basis estimation. We also visualize the
impact of further reducing the number of included MFPC based random effects by
lowering the truncation order $M$ to at least $99\%$ explained variance of the
Gaussian process in a third modeling approach (TRUNC). Finally, we compare our
proposed models to an appropriate model fit (depending on the scenario) in the
widely used JMbayes2 package which uses a corrected two-stage approach for
their Bayesian model estimation (JMB).[28] In particular, we incorporate the
knowledge of truly linear random effects and specify a random intercept –
random slope model with unstructured covariance matrix for the JMB model in
scenario I. In scenario II, the covariance structure is modeled by including a
subject and outcome specific random intercept as well as three cubic B-spline
basis functions resulting in an $8\times 8$ unstructured covariance matrix. We
acknowledge that this modeling approach might not be flexible enough for the
specified covariance structure of the data generating process. However,
increasing the number of basis functions to four already leads to numerical
problems in more than $75\%$ of the JMbayes2 simulation runs, whereas with the
proposed number only one model cannot be estimated in JMbayes2. Note that the
bamlss models do not show any numerical problems in our simulation. In both
frameworks, the baseline hazard is modeled using cubic P-splines with third
order difference penalty and 20 basis functions (before centering).

We run the MCMC algorithm for each model for 5500 iterations and use a burnin
of 500 samples and a thinnig of 5 to obtain a total of 1000 samples per model.
The model fits are evaluated using bias, root mean squared error (rMSE), and
frequentist coverage of the $95\%$ credible intervals of the additive
predictors. In particular, we define the average bias for the longitudinal
predictors $l\in\\{\mu_{k},\sigma_{k}|k=1,...,K\\}$ in simulation run
$u=1,...,200$ as
$Bias^{u}_{l}=\frac{1}{N_{k}^{u}}\sum_{i=1}^{n}\sum_{j=1}^{n^{u}_{ik}}(\hat{\eta}^{u}_{li}(t_{ijk}^{u})-\eta^{u}_{li}(t_{ijk}^{u}))$,
with estimate $\hat{\eta}^{u}_{li}$ and the superscript $u$ denoting
simulation run specific analogues. We also evaluate the model fit over time
using a time dependent bias defined as
$Bias^{u}_{l}(t)=\frac{1}{n_{risk}^{u}(t)}\sum_{i\in\mathcal{R}^{u}(t)}(\hat{\eta}^{u}_{li}(t)-\eta^{u}_{li}(t))$
for all $t$ in $\mathcal{P}$, where $n_{risk}^{u}(t)$ is the cardinality of
the risk set of subjects $\mathcal{R}^{u}(t)$ at time $t$. Note that for the
survival submodel, we evaluate the survival predictors $\eta_{\lambda}$ and
$\eta_{\gamma}$ jointly using their sum $\eta_{\lambda+\gamma}$, as the two
frameworks handle intercepts and sum-to-zero constraints for smooth effects
differently. For $l\in\\{\lambda+\gamma,\alpha_{1},...,\alpha_{K}\\}$ we
calculate the average bias
$Bias^{u}_{l}=\frac{1}{n}\sum_{i=1}^{n}(\eta_{li}^{u}(T_{i}^{u})-\hat{\eta}^{u}_{li}(T_{i}^{u}))$
at the subject specific follow-up times as well as the time dependent bias
$Bias^{u}_{l}(t)$ as above. We define the rMSE and the coverage analogously
for all predictors, for instance
$rMSE^{u}_{l}=(\frac{1}{n}\sum_{i=1}^{n}(\hat{\eta}^{u}_{li}(T_{i}^{u})-\eta^{u}_{li}(T_{i}^{u}))^{2})^{1/2}$
and
$Coverage_{l}^{u}=\frac{1}{n}\sum_{i=1}^{n}I(\hat{\eta}_{li}^{u,2.5}(T_{i}^{u})\leq\eta^{u}_{li}(T_{i}^{u})\leq\hat{\eta}_{li}^{u,97.5}(T_{i}^{u}))$
with indicator function $I(\cdot)$ and $\hat{\eta}^{u,\alpha}$ the $\alpha\%$
quantile of the MCMC samples for the posterior mean for the survival
predictors $l\in\\{\lambda+\gamma,\alpha_{1},...,\alpha_{K}\\}$. Naturally,
these definitions simplify correspondingly for the predictors
$l\in\\{\alpha_{1},...,\alpha_{K},\sigma_{1},...,\sigma_{K}\\}$, which are
constant over all subjects and time in our simulation. The evaluation criteria
of the 200 (or 199 for the JMbayes2 approach in scenario II) simulation runs
are then averaged for presentation.

Additionally, we evaluate the estimated MFPCs in both scenarios.
Eigenfunctions are generally only defined up to a sign change, which is why we
first use the norm $||\cdot||$ induced by (5) (with weights chosen as $1$
throughout the simulation) to determine whether a reflected version
$(-\hat{\bm{\psi}}_{m}(t))$ is closer to the true eigenfunction
$\bm{\psi}_{m}(t),m=1,...,M$. We then use the squared norm of the error,
$||\bm{\psi}_{m}-\hat{\bm{\psi}}_{m}||^{2}$, to evaluate the estimated MFPC.

###  4.4 Simulation results

Table 1 contains the results of simulation scenario I. Comparing the bamlss
model TRUE with true underlying eigenfunctions as basis functions to the
JMbayes2 model JMB allows to directly compare the two frameworks for MJMs,
contrasting the different choices of basis and correlation structure without
influence of the particular method used to estimate the MFPCs. More precisely,
while both approaches predict the two random effects for each subject in their
model fit, the MFPC based approach only estimates the random effect variances
whereas JMbayes2 also estimates the correlation between the random intercepts
and slopes on all longitudinal outcomes. Overall, we find that the two
approaches yield similar bias and rMSE values in the linear scenario I. In
more detail, bamlss performs slightly better for the survival predictors
$\eta_{\lambda}+\eta_{\gamma}$ and the longitudinal predictors
$\eta_{\mu_{1}},...\eta_{\mu_{6}}$, while the JMbayes2 corrected two-stage
approach results in better estimates of the association predictors
$\eta_{\alpha_{1}},...,\eta_{\alpha_{6}}$ compared to the fully Bayesian
approach in bamlss. Some of this difference in performance for the
$\eta_{\alpha_{k}}$ might also be caused by the relatively short MCMC chains
to lower the computational burden of the simulation. The frequentist coverage
of both frameworks is close to the nominal value of $95\%$ for all predictors
except the combined survival predictor containing the baseline hazard. Bender
et al. (2018)[2] point out that using P-splines with a single penalty/variance
parameter can be suboptimal when the baseline hazard changes quickly over the
follow-up period and suggest to use adaptive spline smooths instead, which is
only available (at increased computational cost) for bamlss and not for
JMbayes2, and thus not used in this comparison.

Table 1: Results of simulation scenario I. The bias and rMSE values for $\eta_{\mu_{1}},...,\eta_{\mu_{6}}$ are multiplied by factor $100$ for better readability. The values for $\eta_{\sigma_{1}},...,\eta_{\sigma_{6}}$ are averaged and reported as $\eta_{\sigma}$. | Bias | rMSE | Coverage
---|---|---|---
$\eta$ | TRUE | EST | TRUNC | JMB | TRUE | EST | TRUNC | JMB | TRUE | EST | TRUNC | JMB
$\lambda+\gamma$ | 0.811 | 0.812 | 0.812 | 1.006 | 0.907 | 0.911 | 0.909 | 1.124 | 0.438 | 0.433 | 0.436 | 0.465
$\alpha_{1}$ | 0.118 | 0.106 | 0.126 | 0.042 | 0.313 | 0.315 | 0.338 | 0.283 | 0.945 | 0.930 | 0.920 | 0.955
$\alpha_{2}$ | 0.014 | 0.014 | 0.012 | 0.003 | 0.256 | 0.261 | 0.267 | 0.242 | 0.950 | 0.955 | 0.965 | 0.970
$\alpha_{3}$ | 0.047 | 0.035 | 0.029 | 0.006 | 0.214 | 0.211 | 0.220 | 0.199 | 0.955 | 0.970 | 0.965 | 0.965
$\alpha_{4}$ | \- 0.053 | -0.046 | -0.052 | -0.046 | 0.208 | 0.213 | 0.216 | 0.198 | 0.945 | 0.945 | 0.935 | 0.955
$\alpha_{5}$ | -0.061 | -0.056 | -0.054 | -0.051 | 0.201 | 0.199 | 0.199 | 0.186 | 0.965 | 0.955 | 0.955 | 0.970
$\alpha_{6}$ | -0.168 | -0.157 | -0.158 | -0.106 | 0.259 | 0.251 | 0.256 | 0.227 | 0.905 | 0.895 | 0.905 | 0.960
$\mu_{1}$ | 0.003 | 0.004 | 0.003 | 0.004 | 2.510 | 4.785 | 9.786 | 2.516 | 0.949 | 0.843 | 0.673 | 0.948
$\mu_{2}$ | 0.001 | 0.000 | -0.000 | -0.004 | 2.516 | 4.872 | 9.005 | 2.522 | 0.948 | 0.842 | 0.715 | 0.946
$\mu_{3}$ | -0.004 | -0.003 | -0.003 | -0.010 | 2.507 | 5.079 | 8.250 | 2.511 | 0.950 | 0.839 | 0.755 | 0.947
$\mu_{4}$ | 0.004 | 0.004 | 0.004 | 0.002 | 2.523 | 5.363 | 6.944 | 2.528 | 0.949 | 0.834 | 0.796 | 0.947
$\mu_{5}$ | -0.002 | -0.002 | -0.002 | -0.003 | 2.518 | 5.463 | 6.377 | 2.522 | 0.950 | 0.831 | 0.807 | 0.949
$\mu_{6}$ | 0.011 | 0.012 | 0.011 | 0.011 | 2.527 | 5.914 | 6.240 | 2.531 | 0.950 | 0.827 | 0.817 | 0.948
$\sigma$ | 0.002 | 0.263 | 0.467 | 0.001 | 0.016 | 0.263 | 0.467 | 0.016 | 0.945 | 0.056 | 0.028 | 0.938

Estimating the MFPC basis functions from the data results in a clear increase
in the variance of the estimated longitudinal submodel, due to the higher
flexibility of our nonparametric approach, which gives an advantage to the
correctly specified parametric JMB and TRUE models (the latter using the true
linear MFPCs). The higher variability is reflected in higher rMSE values of
the EST models, while the bias values for the longitudinal predictors remain
stable. The reason is that flexibly estimating the eigenfunctions from the
data in scenario I does not necessarily result in linear estimated MFPCs (see
Appendix B Figure 6), which can lead to longitudinal model fits that are too
wiggly compared to the true parametric random intercept – random slope data
generating process assumed known in JMB and TRUE models. This mixing of the
correlation structure and the random noise also leads to an overestimation of
the variance predictors $\eta_{\sigma_{1}},...,\eta_{\sigma_{6}}$ (averaged in
Table 1) coupled with poor frequentist coverage in the longitudinal submodel.
This, however, does not negatively affect the estimation of the survival
submodel where the bias, rMSE, and coverage values are similar to the TRUE
models, including for the association predictors
$\eta_{\alpha_{1}},...,\eta_{\alpha_{6}}$. Note that the bias and rMSE values
for the longitudinal predictors in Table 1 have been multiplied by a factor
$100$, such that the longitudinal model fits using the estimated MFPCs still
are visually close to the true trajectories (see Appendix B Figure 5) and seem
close enough to correctly estimate the association predictors well.

Lowering the number of MFPC basis functions included in the model by
introducing truncation at a prespecified rate of explained variance of the
Gaussian process is expected to lower the computational burden, but to
deteriorate the longitudinal submodel fit. In particular, for a truncation at
$99\%$ explained variance, the TRUNC models offer a notable dimension
reduction, as most models only include ten MFPCs (minimum nine and maximum
eleven) instead of the full 12, thus cutting down the models by $2\times 150$
parameters. While the bias values of the longitudinal predictors are basically
unaffected by the truncation, we note that the rMSE values of the longitudinal
predictors are up to four times higher while the frequentist coverage drops to
about $70-80\%$. The variation of the Gaussian process now unaccounted for by
the excluded MFPCs is absorbed by the variance predictors
$\eta_{\sigma_{1}},...,\eta_{\sigma_{6}}$, which are overestimated throughout
all TRUNC models. The TRUNC models show, however, that the trailing MFPCs with
small variances mainly affect the fit of the longitudinal submodel (see
Appendix B Figure 5) while the association predictors are still estimated
reasonably well. Truncation thus allows an explicit choice and trade-off
between computational complexity and (longitudinal) model fit, with the
ordering of eigenvalues ensuring that components of little importance are
truncated first.

In simulation scenario II, the parametric covariance specification of scenario
I is replaced by a complex multivariate covariance structure. In this setting,
the JMB models based on the JMbayes2 framework no longer have the advantage of
“knowing” the correct parametric model. In fact, we find that the restricted
flexibility of the JMB models (as mentioned above) leads to a worse overall
model fit compared to TRUE models as measured by the bias, rMSE, and
frequentist coverage values for (almost) all predictors as given in Table 2.
The difference between the two models is especially pronounced for the
longitudinal predictors $\eta_{\mu_{1}},\eta_{\mu_{2}}$, and
$\eta_{\sigma_{1}},\eta_{\sigma_{2}}$, but also for the association predictors
$\eta_{\alpha_{1}},\eta_{\alpha_{2}}$. A negative bias for positive values of
the underlying association in simulation II indicates that the JMB models
underestimate the association parameters. Figure 1 illustrates for three
random subjects that the TRUE model using the underlying eigenfunctions as
basis functions fits closely to the true trajectories, while the JMB model
tends to oversmooth the data, which may explain its slight underestimation of
the $\alpha_{1},\alpha_{2}$ association parameters also known to occur in
regression settings with covariate measurement error[4]. Especially for
subjects with a short follow-up time (such as subject $56$), principal
component based random effects may be more capable of pooling the information
over all subjects and across markers to approximate the true trajectories,
while in the example of subject $56$ in Figure 1 the JMB model even fits an
incorrectly increasing time trend based on the few data points available.

Table 2: Results of simulation scenario II. The bias and rMSE values for $\eta_{\mu_{1}},\eta_{\mu_{2}}$ are multiplied by factor $100$ for better readability. | Bias | rMSE | Coverage
---|---|---|---
$\eta$ | TRUE | EST | TRUNC | JMB | TRUE | EST | TRUNC | JMB | TRUE | EST | TRUNC | JMB
$\lambda+\gamma$ | 1.087 | 1.129 | 1.135 | 1.473 | 1.185 | 1.223 | 1.228 | 1.545 | 0.154 | 0.147 | 0.148 | 0.103
$\alpha_{1}$ | 0.019 | 0.012 | 0.012 | -0.053 | 0.070 | 0.072 | 0.073 | 0.098 | 0.945 | 0.935 | 0.950 | 0.849
$\alpha_{2}$ | 0.010 | 0.005 | 0.002 | -0.043 | 0.070 | 0.074 | 0.079 | 0.104 | 0.990 | 0.960 | 0.955 | 0.859
$\mu_{1}$ | -0.004 | -0.004 | -0.006 | -0.002 | 2.894 | 11.271 | 12.808 | 13.013 | 0.944 | 0.795 | 0.770 | 0.810
$\mu_{2}$ | -0.027 | -0.028 | -0.027 | -0.036 | 2.850 | 8.917 | 9.752 | 20.246 | 0.934 | 0.783 | 0.776 | 0.791
$\sigma_{1}$ | 0.021 | 0.805 | 0.894 | 0.955 | 0.027 | 0.805 | 0.894 | 0.955 | 0.835 | 0.000 | 0.000 | 0.000
$\sigma_{2}$ | 0.019 | 0.625 | 0.662 | 1.364 | 0.023 | 0.625 | 0.662 | 1.364 | 0.790 | 0.000 | 0.000 | 0.000
Figure 1: Different model fits of the longitudinal submodel for three random
subjects of a simulated data set in scenario II. The solid black line
corresponds to the true underlying longitudinal trajectory (mu) and the grey
dots are the observed data points. The TRUNC model in this example includes
five MFPCs instead of the correct six.

Similar to scenario I, the estimation of the MFPC basis mainly increases the
variance of the estimates in the longitudinal submodel, with rMSE values still
lower than in the JMbayes2 framework. The corresponding coverage drops to
about $80\%$ for the longitudinal predictors and the variance predictors are
overestimated with a coverage of $0$, both similar to the JMB model. Figure 1
exemplifies that the longitudinal submodel fits are still sufficiently close
to the true trajectories for the estimation of the association predictors
$\eta_{\alpha_{1}},\eta_{\alpha_{2}}$ not to suffer from the estimation of the
MFPC basis. Truncating the basis to $99\%$ of explained variance again entails
a slight deterioration in the overall model fit as reported by the bias, rMSE,
and coverage. Note, however, that in $64\%$ of all simulation runs all six
MFPCs are included in the TRUNC models. Further decreasing the truncation
threshold to $95\%$ (results not shown) lowers the computational burden by
reducing the number of MFPC basis functions from six to on average four, thus
dropping $2\times 300$ parameters. The resulting decline in estimation
accuracy, however, is so considerable that the rMSE values are then higher
than for the oversmoothed JMB models for both the longitudinal and the
association predictors. We conclude that truncation needs to be carefully
considered and that a large percentage of variation should be retained unless
prohibited by computational considerations.

The time dependent evaluation in both simulation scenarios (results not shown)
reveals that the estimation of the baseline hazard with P-splines in both
frameworks leads to oversmoothing at the beginning of the time interval, where
the baseline hazard is quickly changing (cf. Bender et al., 2018[2]),
explaining results in the first rows of Tables 1 and 2. For the longitudinal
submodel, the drop-out process implies less information at the end of the time
interval, resulting in decreasing estimation accuracy. As the drop-out process
also affects the estimation of the MFPC basis functions and as this additional
uncertainty is not reflected in the EST and TRUNC models, this also results in
decreasing coverage over time.

The results from both simulation scenarios reinforce that the choice and for
MFPCs also estimation of the basis is an essential prerequisite for a
satisfactory model fit. We see in scenario I that if the covariance structure
is known to be parametric a priori, a parametric model will always perform
better than our more flexible nonparametric specification, which by estimating
the linear eigenfunctions introduces additional noise into the MFPCs. On the
other hand, parametric or less flexible spline specifications are not
sufficient to capture more complex trends and result in oversmoothing, such as
in scenario II, where our more flexible approach has a clear advantage. We
find that estimated MFPCs (see Appendix B Figure 7) in most simulation runs
recover the relevant characteristics of the eigenfunctions and thus time
trends well, with the leading MFPCs closer to the true eigenfunctions in
simulation scenario II (cf. Table 4 in the Appendix). This illustrates again
that the strength of our proposed MFPC based approach lies in the estimation
of complex covariance structures. For both simulation scenarios, we find that
lowering the truncation threshold can considerably decrease the estimation
accuracy of the longitudinal submodel and even negatively affect the survival
submodel. We therefore advocate to use a high rate of explained variance, such
as $99\%$, for determining the number of MFPCs in a model.

##  5 Application

We demonstrate our parsimonious but flexible modeling approach by analyzing a
PBC study, where we want to model the association between different
biochemical markers and the survival of the patients, and compare it to the
more parameter intensive modeling approach of using spline based random
effects to capture the (potentially) complex covariance structure in the data.
The code to reproduce the analysis is provided in the accompanying R package
MJMbamlss. The PBC data provided in the R package JMbayes2 include information
on $312$ subjects, in particular on the follow-up time (composite event of
liver transplantation or death), baseline covariates (sex, age, and treatment
group), and several longitudinally measured biochemical outcomes which are
suspected to be associated with disease progression. During the study period,
visits were scheduled but missing observations and drop out lead to irregular
measurements of the longitudinal outcomes with an average of $6.2$ visits per
subject. At the end of the study period, $169$ subjects had experienced the
event with median follow-up time of $4.5$ years (IQR: $2.2$, $6.5$) and $143$
subjects were event-free with median follow-up time of $8.7$ (IQR: $6.5$,
$10.5$). For our analysis, we include the continuous longitudinal outcomes
albumin in g/dl, serum bilirubin (serBilir) in mg/dl, serum cholesterol
(serChol) in mg/dl, and aspartate aminotransferase (SGOT) in U/liter, all log-
transformed to improve adherence to the normal assumption.[34] For a more
direct comparison to the JMbayes2 framework, eight subjects for whom no
measurements of serum cholesterol are available are excluded from the
analysis, giving a total of $304$ subjects. We use all available baseline
information on the subjects such as age at registration in years, sex (male or
female), and treatment group drug (placebo or D-penicillamine) for modeling
the hazard as well as the longitudinal outcomes. Note that we cap the follow-
up times $t$ at the maximum observation time of longitudinal measurements in
the data set ($14.1$ years, only affecting three event-free subjects) to avoid
extrapolation of the longitudinal trajectories.

We apply the proposed MFPC based approach (denoted as bamlss) to the PBC data
and compare the results to a similarly specified model fit of the JMbayes2
package (denoted as JMB). These models follow the specifications (1) – (4)
with

$\displaystyle\eta_{\gamma i}$ $\displaystyle=\beta_{\gamma 0}+\beta_{\gamma
1}\cdot\text{sex}_{i}+\beta_{\gamma
2}\cdot\text{drug}_{i}+f_{\gamma}(\text{age}_{i})\,,$
$\displaystyle\eta_{\alpha_{k}i}(t)$ $\displaystyle=\beta_{\alpha_{k}}\,,$
$\displaystyle\eta_{\mu_{k}i}(t)$
$\displaystyle=\beta_{k0}+\beta_{k1}\cdot\text{sex}_{i}+\beta_{k2}\cdot\text{drug}_{i}+f_{k1}(\text{age}_{i})+f_{k2}(t)+b_{i}^{(k)}(t)\,,\,\text{and}$
$\displaystyle\eta_{\sigma_{k}i}(t)$ $\displaystyle=\beta_{\sigma_{k}}$

for
$k\in\\{\textit{albumin},\textit{serBilir},\textit{serChol},\textit{SGOT}\\}$,
and subject index $i$. To include a smooth time trend for the longitudinal
outcomes and to allow for possibly nonlinear covariate effects of age, we use
P-splines with ten basis functions before centering to estimate the additive
functions $f_{\gamma}$, $f_{k1}$, and $f_{k2}$ with bamlss. For the JMbayes2
model we try to fully exploit the framework’s flexibility but for a natural
cubic B-spline basis, we reach the limit with three basis functions, as more
basis functions result in error messages. Similarly, we use the MFPC based
representation (6) to estimate the longitudinal random effects in our
approach, while the JMbayes2 model maximally allows random intercepts and a
natural cubic splines representation of the follow-up time with two basis
functions, giving three random effects per subject and longitudinal outcome.
The baseline hazard $\eta_{\lambda i}(t)$ is modeled in both frameworks using
cubic P-splines with second order difference penalty and 10 basis functions
(before centering).

For the estimation of the MFPC basis, we follow the outline given in section
3.1. In order to obtain stable estimates of univariate scores, the UFPCAs are
based on a subset of 145 subjects that have at least three longitudinal
measurements and one measurement after $10\%$ of the follow-up time interval.
We use seven marginal B-spline basis functions for the estimation of the
univariate covariance functions and truncate the UFPCA at $99\%$ explained
univariate variance, which results in three UFPCs per longitudinal outcome.
The resulting eigenvalues of the UFPCAs show a considerable imbalance in
univariate variation in the data with marker serBilir exhibiting an integrated
univariate variance more than 70 times larger than that of albumin. By using a
scalar product (5) with weights inversely proportional to the integrated
univariate variance for the MFPCA, we account for the different amounts of
variation equivalent to standardizing the longitudinal data without actually
transforming the outcome. The resulting MFPCs are then more balanced across
the longitudinal outcomes with the first two explaining $51.2\%$ and $22.8\%$
of the weighted variance, respectively, while using scalar product weights of
$1$ would give a first MFPC dominated by serBilir accounting for $79.5\%$ of
the (unweighted) multivariate variance. In a sensitivity analysis (see
Appendix C) we use the estimated MFPCs of the unweighted MFPCA and find that
the model fits better to the serBilir data but worse to the other markers
giving an overall worse model fit. Visual inspection of the estimated MFPCs
(see Appendix C Figure 8) shows that using seven B-spline basis functions in
the univariate specifications yields a good compromise between overly smooth
and too flexible (trailing) MFPCs. To decrease the computational burden of the
final model fit, we truncate the MFPC basis at $99\%$ explained variance,
which leads to nine included MFPCs for the weighted MFPCA models, i.e. roughly
a parsimonious two basis functions per marker.

Figure 2: MFPC based (bamlss) and spline based (JMB) model fits of the PBC
data for four subjects. The grey dots are the observed data points.
Longitudinal outcomes are log-transformed.

For each model, we run one MCMC chain with 12,000 iterations, a burnin of
2,000 and a thinning of 5, yielding 2,000 samples. Overall, we find that the
MFPC based bamlss model yields similar results to the spline based random
effects approach of the JMB model. In particular, Table 3 shows that we find
significant association of the time-to-event process with the current value of
each biochemical marker of our analysis, irrespective of the modeling approach
used. This finding aligns with other works in the literature[1, 44, 47], who
might have come to slightly different conclusions due to different model
specifications. Other estimated parameters such as the linear effects of sex
on the hazard or the longitudinal outcomes are also similar for the bamlss and
JMB models (see Appendix C Table 5). What is more, Figure 2 illustrates that
the longitudinal submodel fits of both models are generally very close. Given
that the estimated MFPC basis could provide much more flexibility than the
spline approach with only two basis functions, this suggests that the PBC data
might not require the added flexibility supplied by an MFPC basis
representation. Nevertheless, the bamlss model uses much fewer parameters for
the estimation of the random effects structure despite the larger flexibility.
While the spline based approach estimates $3\times 4\times 304$ random effects
with an unstructured covariance matrix of $78$ covariance parameters, the MFPC
based approach only needs $9\times 304$ random effects with a diagonal
covariance matrix of $9$ variance parameters, thus saving $3\times 304+69=981$
parameters total. We therefore show that using an MFPC basis not only reduces
the computational costs of fitting an MJM, but also constitutes an
alternative, highly flexible modelling approach to spline based random
effects, especially when numerical problems restrict their usage for nonlinear
covariance structures. Then the proposed approach also allows to check whether
the assumptions made in a spline model with few basis functions dictated by
computational possibilities are in fact warranted.

Table 3: Posterior mean and $95\%$ credible interval of the association predictors $\beta_{\alpha_{k}}$ for the two models also shown in Figure 2. Longitudinal outcomes are log-transformed. | bamlss | JMB
---|---|---
albumin | $-5.82$ |  $[-7.84$; $-3.89]$ | $-6.03$ |  $[-7.87$; $-4.24]$
serBilir | $1.48$ |  $[1.19$; $1.77]$ | $1.39$ |  $[1.11$; $1.68]$
serChol | $-0.76$ |  $[-1.36$; $-0.16]$ | $-0.68$ |  $[-1.26$; $-0.06]$
SGOT | $-0.74$ |  $[-1.41$; $-0.13]$ | $-0.76$ |  $[-1.4$; $-0.11]$

##  6 Discussion and Outlook

We have shown that shared parameter MJMs can greatly benefit from using an
MFPC based representation of the random effects structure. In particular, the
flexibility when modeling highly nonlinear subject specific trends in the
multivariate longitudinal outcomes is often restricted by the computational
complexity of the number of parameters to estimate, especially when a spline
approach is used. Using MFPCs as a basis for random effects allows a
parsimonious representation of increased flexibility, reducing the number of
parameters to estimate by decoupling it from the number of longitudinal
outcomes and simplifying the covariance structure of the random effects. Our
simulation study showcases that the MFPC based approach can achieve a better
model fit for nonlinear longitudinal outcome trajectories resulting also in a
more accurate estimation of the association structure with the event hazard.
Furthermore, the explicit ordering of the MFPCs in terms of explained variance
makes the trade-off between computational burden and goodness of fit explicit,
which allows its steering by the selection of the truncation level for the
number of MFPC based random effects. Note that while not focused on here,
dynamic prediction for subjects with incoming new longitudinal measurements is
directly available from our approach since the estimated MFPCs are used as
empirical basis functions also for all (new) subjects in the proposed MJM, and
the standard algorithms for dynamic prediction of the random scores can thus
be used.[41]

We also provide a working implementation for the framework of flexible MJMs in
the bamlss package, which supports a vast range of possible specifications for
covariate effects in a user friendly way. Note, however, that the idea of an
MFPC basis is more universal and can be combined with different approaches to
shared parameter MJMs such as approximations using adapted two-stage
models[28] in the JMbayes2 package or integrated nested Laplace (INL) [44].
For JMbayes2, both the simulation study and our application have shown that
the model flexibility is restricted by the number of spline basis functions
which can be included in the longitudinal submodel. This could be mitigated by
implementing our proposed MFPC approach also in JMbayes2. Additionally,
Rustand et al. (2020)[44] remark that simplifying the random effects
covariance matrix to be block-diagonal considerably reduces computation times.
While this corresponds to the strong assumption of independent random effects
across markers, an MFPC basis does in fact provide the further simplification
of a diagonal covariance matrix without the necessity of such an independence
assumption. On the other hand, both the adapted two-stage and the INL
approximations to the full MJM present computational advantages to the
derivative based MCMC sampling implemented in bamlss. In simulation scenario
I, for example, the fully Bayesian approach in bamlss estimates the MJM in
around $30$ minutes on a Linux system with $12$ cores at $3.5$ GHz and $32$ GB
memory while the highly optimized two-stage approach in JMbayes2 returns
results after roughly one minute with the same setup. This difference might in
future be diminished by increasing computational efficiency in the implemented
code via transfering more computations to other languages such as C or
exploiting sparse structures in the design matrices. A combination of an MFPC
approach and approximate approaches also holds promise.

It is important to note, however, that estimating the MFPC basis is an
integral part of our proposed model estimation approach. Besides model
selection criteria such as the Bayesian information criterion or the deviance
information criterion, we find that examination of the amount of variance
explained and visual inspection of the estimated UFPCs and MFPCs is helpful to
select the number of MFPCs and assure that the fitted longitudinal
trajectories in the MJM exhibit a reasonable amount of flexibility. Our
current approach of estimating the MFPC basis ignores the informative
missingness of the longitudinal process by treating the longitudinal outcomes
as sparsely observed functional data. While we are able to show that the
resulting estimated MFPCs serve as a sufficiently good empirical basis in the
MJM, we encourage further research on estimating UFPCs and MFPCs under
informative censoring. Several ad-hoc adjustments might improve the efficient
use of all available longitudinal information such as weighting univariate
outcomes depending on their length of follow-up in the estimation of the
univariate covariance structure instead of excluding subjects with short
follow-up time from the MFPCA entirely. Dong et al. (2023)[7] suggest that the
UFPC estimation under informative censoring might be improved by iteratively
estimating UFPCs and predicting the full longitudinal trajectories, as this
might stabilize estimation. Alternatively, separate univariate joint models
based on reduced rank models, which simultaneously estimate the time-to-event
process and the UFPCs as in Yao (2007)[57], could be employed, but would
considerably increase computational burden. Further research is needed to
evaluate and compare different approaches to estimate UFPCs under informative
missingness.

As a consequence of the proposed two step approach, all model based inference
is conditional on the estimated MFPCs. Since Happ and Greven [14] show that
the MFPC estimates asymptotically converge to the true eigenfunctions, the
inference can be expected to be asymptotically correct. Future work is needed
to account for the uncertainty in the estimation of the MFPC basis in finite
samples, for example using corrections along the lines of Goldsmith et al.
(2013)[11], which would improve the frequentist coverage properties of the
longitudinal submodel fits. Developing an estimation approach which jointly
estimates the scores and MFPCs within an MJM might further enhance modeling
and understanding of the longitudinal processes, also allowing for inference
regarding the estimated MFPCs.

Given that an MFPC based representation of random effects is a universal and
effective approach for MJMs, future work can follow several more research
avenues. From extending the presented joint model to incorporate functional
covariates, functional outcomes, or different parameterizations of the
association structure between the longitudinal and time-to-event submodels, to
developing efficient approximations to the full derivative based MCMC
algorithm in our bamlss implementation or advancing the MFPC basis approach to
non-normal longitudinal outcomes, MJMs remain an interesting and challenging
field of study.

##  7 Acknowledgements

We gratefully acknowledge funding by grant GR 3793/3-1 ‘Flexible regression
methods for curve and shape data’ from the German Research Foundation (DFG).

## References

  * Andrinopoulou and Rizopoulos [2016] E.-R. Andrinopoulou and D. Rizopoulos.  Bayesian shrinkage approach for a joint model of longitudinal and survival outcomes assuming different association structures.  _Statistics in Medicine_ , 35(26):4813–4823, 2016.
  * Bender et al. [2018] A. Bender, A. Groll, and F. Scheipl.  A generalized additive model approach to time-to-event analysis.  _Statistical Modelling_ , 18(3-4):299–321, 2018.
  * Bender et al. [2005] R. Bender, T. Augustin, and M. Blettner.  Generating survival times to simulate cox proportional hazards models.  _Statistics in Medicine_ , 24(11):1713–1723, 2005.
  * Carroll et al. [2006] R. J. Carroll, D. Ruppert, L. A. Stefanski, and C. M. Crainiceanu.  _Measurement error in nonlinear models: a modern perspective_.  Chapman and Hall/CRC, 2006.
  * Cederbaum et al. [2018] J. Cederbaum, F. Scheipl, and S. Greven.  Fast symmetric additive covariance smoothing.  _Computational Statistics & Data Analysis_, 120:25–41, 2018.
  * Crowther and Lambert [2013] M. J. Crowther and P. C. Lambert.  Simulating biologically plausible complex survival data.  _Statistics in Medicine_ , 32(23):4118–4134, 2013.
  * Dong et al. [2023] J. Dong, H. Shi, L. Wang, Y. Zhang, and J. Cao.  Jointly modelling multiple transplant outcomes by a competing risk model via functional principal component analysis.  _Journal of Applied Statistics_ , 50(1):43–59, 2023.
  * Faucett and Thomas [1996] C. L. Faucett and D. C. Thomas.  Simultaneously modelling censored survival data and repeatedly measured covariates: a gibbs sampling approach.  _Statistics in Medicine_ , 15(15):1663–1685, 1996.
  * Fleming and Harrington [2013] T. R. Fleming and D. P. Harrington.  _Counting processes and survival analysis_ , volume 625.  John Wiley & Sons, 2013.
  * Goldschmidt et al. [2018] I. Goldschmidt, A. Karch, R. Mikolajczyk, F. Mutschler, N. Junge, E. D. Pfister, T. Möhring, L. d’Antiga, P. McKiernan, D. Kelly, et al.  Immune monitoring after pediatric liver transplantation–the prospective chilsfree cohort study.  _BMC Gastroenterology_ , 18:1–9, 2018.
  * Goldsmith et al. [2013] J. Goldsmith, S. Greven, and C. Crainiceanu.  Corrected confidence bands for functional data using principal components.  _Biometrics_ , 69(1):41–51, 2013.
  * Goodrich et al. [2020] B. Goodrich, J. Gabry, I. Ali, and S. Brilleman.  rstanarm: Bayesian applied regression modeling via Stan.  ”https://mc-stan.org/rstanarm”, 2020.  R package version 2.21.1.
  * Hakola et al. [2019] L. Hakola, M. E. Miettinen, E. Syrjälä, M. Åkerlund, H.-M. Takkinen, T. E. Korhonen, S. Ahonen, J. Ilonen, J. Toppari, R. Veijola, et al.  Association of cereal, gluten, and dietary fiber intake with islet autoimmunity and type 1 diabetes.  _JAMA Pediatrics_ , 173(10):953–960, 2019.
  * Happ and Greven [2018] C. Happ and S. Greven.  Multivariate functional principal component analysis for data observed on different (dimensional) domains.  _Journal of the American Statistical Association_ , 113(522):649–659, 2018.
  * Hickey et al. [2018] G. L. Hickey, P. Philipson, A. Jorgensen, and R. Kolamunnage-Dona.  joinerml: a joint model and software package for time-to-event and multivariate longitudinal outcomes.  _BMC Medical Research Methodology_ , 18:1–14, 2018.
  * Holte et al. [2012] S. Holte, T. Randolph, J. Ding, J. Tien, R. McClelland, J. Baeten, and J. Overbaugh.  Efficient use of longitudinal cd4 counts and viral load measures in survival analysis.  _Statistics in Medicine_ , 31(19):2086–2097, 2012.
  * Hong et al. [2021] Y. Hong, L. Su, S. Song, and F. Yan.  Dynamic prediction of disease processes based on recurrent history and functional principal component analysis of longitudinal biomarkers: Application for ovarian epithelial cancer.  _Statistics in Medicine_ , 40(8):2006–2023, 2021.
  * Kang and Song [2023] K. Kang and X. Y. Song.  Joint modeling of longitudinal imaging and survival data.  _Journal of Computational and Graphical Statistics_ , 32(2):402–412, 2023.
  * Köhler et al. [2017] M. Köhler, N. Umlauf, A. Beyerlein, C. Winkler, A.-G. Ziegler, and S. Greven.  Flexible bayesian additive joint models with an application to type 1 diabetes research.  _Biometrical Journal_ , 59(6):1144–1165, 2017.
  * Köhler et al. [2018] M. Köhler, N. Umlauf, and S. Greven.  Nonlinear association structures in flexible bayesian additive joint models.  _Statistics in Medicine_ , 37(30):4771–4788, 2018.
  * Kong et al. [2018] D. Kong, J. G. Ibrahim, E. Lee, and H. Zhu.  Flcrm: Functional linear cox regression model.  _Biometrics_ , 74(1):109–117, 2018.
  * Li et al. [2022] C. Li, L. Xiao, and S. Luo.  Joint model for survival and multivariate sparse functional data with application to a study of alzheimer’s disease.  _Biometrics_ , 78(2):435–447, 2022.
  * Li and Luo [2017] K. Li and S. Luo.  Functional joint model for longitudinal and time-to-event data: an application to alzheimer’s disease.  _Statistics in Medicine_ , 36(22):3560–3572, 2017.
  * Li and Luo [2019a] K. Li and S. Luo.  Bayesian functional joint models for multivariate longitudinal and time-to-event data.  _Computational Statistics & Data Analysis_, 129:14–29, 2019a.
  * Li and Luo [2019b] K. Li and S. Luo.  Dynamic prediction of alzheimer’s disease progression using features of multiple longitudinal outcomes and time-to-event data.  _Statistics in Medicine_ , 38(24):4804–4818, 2019b.
  * Li et al. [2021] N. Li, Y. Liu, S. Li, R. M. Elashoff, and G. Li.  A flexible joint model for multiple longitudinal biomarkers and a time-to-event outcome: With applications to dynamic prediction using highly correlated biomarkers.  _Biometrical Journal_ , 63(8):1575–1586, 2021.
  * Lin et al. [2021] J. Lin, K. Li, and S. Luo.  Functional survival forests for multivariate longitudinal outcomes: Dynamic prediction of alzheimer’s disease progression.  _Statistical Methods in Medical Research_ , 30(1):99–111, 2021.
  * Mauff et al. [2020] K. Mauff, E. Steyerberg, I. Kardys, E. Boersma, and D. Rizopoulos.  Joint models with multiple longitudinal outcomes and a time-to-event outcome: a corrected two-stage approach.  _Statistics and Computing_ , 30:999–1014, 2020.
  * Mauff et al. [2021] K. Mauff, N. S. Erler, I. Kardys, and D. Rizopoulos.  Pairwise estimation of multivariate longitudinal outcomes in a bayesian setting with extensions to the joint model.  _Statistical Modelling_ , 21(1-2):115–136, 2021.
  * Murray [2023] J. Murray.  Package ’gmvjoint’: Joint models of survival and multivariate longitudinal data.  ”https://github.com/jamesmurray7/gmvjoint”, 2023.
  * Murray and Philipson [2022] J. Murray and P. Philipson.  A fast approximate em algorithm for joint models of survival and multivariate longitudinal data.  _Computational Statistics & Data Analysis_, 170:107438, 2022.
  * Murray and Philipson [2023] J. Murray and P. Philipson.  Fast estimation for generalised multivariate joint models using an approximate em algorithm.  _Computational Statistics & Data Analysis_, page 107819, 2023.
  * Murtaugh et al. [1994] P. A. Murtaugh, E. R. Dickson, G. M. Van Dam, M. Malinchoc, P. M. Grambsch, A. L. Langworthy, and C. H. Gips.  Primary biliary cirrhosis: prediction of short-term survival based on repeated patient visits.  _Hepatology_ , 20(1):126–134, 1994.
  * Pandolfi et al. [2023] S. Pandolfi, F. Bartolucci, and F. Pennoni.  A hidden markov model for continuous longitudinal data with missing responses and dropout.  _Biometrical Journal_ , page 2200016, 2023.
  * Papageorgiou et al. [2019] G. Papageorgiou, K. Mauff, A. Tomer, and D. Rizopoulos.  An overview of joint modeling of time-to-event and longitudinal outcomes.  _Annual Review of Statistics and Its Application_ , 6:223–240, 2019.
  * Philipson et al. [2020] P. Philipson, G. L. Hickey, M. J. Crowther, and R. Kolamunnage-Dona.  Faster monte carlo estimation of joint models for time-to-event and multivariate longitudinal data.  _Computational Statistics & Data Analysis_, 151:107010, 2020.
  * Press et al. [2007] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery.  _Numerical Recipes: The Art of Scientific Computing_.  Cambridge University Press, New York, 3rd edition, 2007.  ISBN 9780521884075.
  * Proust-Lima and Taylor [2009] C. Proust-Lima and J. M. Taylor.  Development and validation of a dynamic prognostic tool for prostate cancer recurrence using repeated measures of posttreatment psa: a joint modeling approach.  _Biostatistics_ , 10(3):535–549, 2009.
  * Ramsay and Silverman [2005] J. O. Ramsay and B. W. Silverman.  _Functional data analysis_.  Springer Science & Business Media, second edition, 2005.
  * Rizopoulos [2012] D. Rizopoulos.  _Joint models for longitudinal and time-to-event data: With applications in R_.  CRC press, 2012.
  * Rizopoulos [2016] D. Rizopoulos.  The r package JMbayes for fitting joint models for longitudinal and time-to-event data using mcmc.  _Journal of Statistical Software_ , 72(7):1–46, 2016.
  * Rizopoulos et al. [2022] D. Rizopoulos, G. Papageorgiou, and P. Miranda Afonso.  JMbayes2: Extended joint models for longitudinal and time-to-event data.  https://github.com/drizopoulos/JMbayes2, 2022.  https://drizopoulos.github.io/JMbayes2/.
  * Rizopoulos et al. [2023] D. Rizopoulos, G. Papageorgiou, and P. Miranda Afonso.  JMbayes2: Extended joint models for longitudinal and time-to-event data.  https://CRAN.R-project.org/package=JMbayes2, 2023.  URL https://CRAN.R-project.org/package=JMbayes2.  R package version 0.4-5.
  * Rustand et al. [2023] D. Rustand, J. van Niekerk, E. T. Krainski, H. Rue, and C. Proust-Lima.  Fast and flexible inference for joint models of multivariate longitudinal and survival data using integrated nested laplace approximations.  _Biostatistics_ , page kxad019, 2023.
  * Scheipl et al. [2015] F. Scheipl, A.-M. Staicu, and S. Greven.  Functional additive mixed models.  _Journal of Computational and Graphical Statistics_ , 24(2):477–501, 2015.
  * Tsiatis and Davidian [2004] A. A. Tsiatis and M. Davidian.  Joint modeling of longitudinal and time-to-event data: an overview.  _Statistica Sinica_ , pages 809–834, 2004.
  * Tu and Sun [2023] J. Tu and J. Sun.  Gaussian variational approximate inference for joint models of longitudinal biomarkers and a survival outcome.  _Statistics in Medicine_ , 42(3):316–330, 2023.
  * Umlauf et al. [2018] N. Umlauf, N. Klein, and A. Zeileis.  Bamlss: Bayesian additive models for location, scale, and shape (and beyond).  _Journal of Computational and Graphical Statistics_ , 27(3):612–627, 2018.
  * Umlauf et al. [2021] N. Umlauf, N. Klein, T. Simon, and A. Zeileis.  bamlss: A Lego toolbox for flexible Bayesian regression (and beyond).  _Journal of Statistical Software_ , 100(4):1–53, 2021.
  * Umlauf et al. [2023] N. Umlauf, N. Klein, A. Zeileis, and T. Simon.  bamlss: Bayesian additive models for location, scale, and shape (and beyond).  http://www.bamlss.org/, 2023.  URL http://www.bamlss.org/.  R package version 1.2-1.
  * Volkmann et al. [2023] A. Volkmann, A. Stöcker, F. Scheipl, and S. Greven.  Multivariate functional additive mixed models.  _Statistical Modelling_ , 23(4):303–326, 2023.
  * Wood [2021] S. Wood.  mgcv: Mixed GAM computation vehicle with automatic smoothness estimation.  https://CRAN.R-project.org/package=mgcv, 2021.  URL https://CRAN.R-project.org/package=mgcv.  R package version 1.8-34.
  * Wood [2017] S. N. Wood.  _Generalized additive models: An introduction with R_.  CRC press, 2017.
  * Wulfsohn and Tsiatis [1997] M. S. Wulfsohn and A. A. Tsiatis.  A joint model for survival and longitudinal data measured with error.  _Biometrics_ , pages 330–339, 1997.
  * Yan et al. [2017] F. Yan, X. Lin, and X. Huang.  Dynamic prediction of disease progression for leukemia patients by functional principal component analysis of longitudinal expression levels of an oncogene.  _The Annals of Applied Statistics_ , 11(3):1649–1670, 2017.
  * Yan et al. [2018] F. Yan, X. Lin, R. Li, and X. Huang.  Functional principal components analysis on moving time windows of longitudinal data: dynamic prediction of times to event.  _Journal of the Royal Statistical Society Series C: Applied Statistics_ , 67(4):961–978, 2018.
  * Yao [2007] F. Yao.  Functional principal component analysis for longitudinal and survival data.  _Statistica Sinica_ , 17(3):965–983, 2007.
  * Yao et al. [2005] F. Yao, H.-G. Müller, and J.-L. Wang.  Functional data analysis for sparse longitudinal data.  _Journal of the American Statistical Association_ , 100(470):577–590, 2005.
  * Ye et al. [2015] J. Ye, Y. Li, and Y. Guan.  Joint modeling of longitudinal drug using pattern and time to first relapse in cocaine dependence treatment data.  _The Annals of Applied Statistics_ , 9(3):1621–1642, 2015.
  * Zou et al. [2021] H. Zou, K. Li, D. Zeng, S. Luo, and A. D. N. Initiative.  Bayesian inference and dynamic prediction of multivariate joint model with functional data: An application to alzheimer’s disease.  _Statistics in Medicine_ , 40(30):6855–6872, 2021.
  * Zou et al. [2023a] H. Zou, L. Xiao, D. Zeng, S. Luo, and A. D. N. Initiative.  Multivariate functional mixed model with mri data: An application to alzheimer’s disease.  _Statistics in Medicine_ , 42(10):1492–1511, 2023a.
  * Zou et al. [2023b] H. Zou, D. Zeng, L. Xiao, and S. Luo.  Bayesian inference and dynamic prediction for multivariate longitudinal and survival data.  _The Annals of Applied Statistics_ , 17(3):2574–2595, 2023b.

##  Appendix A Derivatives

In the following, the blockwise score vectors $\bm{s}(\bm{\beta}_{lh})$ and
Hessian $\bm{\mathcal{H}}(\bm{\beta}_{lh})$ used in the posterior mode and
mean estimation are presented. Note that we show only the derivatives based on
the log-likelihood function as adding the derivatives for the respective log-
priors is straightforward.[19] To highlight the differences between the
marker-specific fixed and principal component based random effects, we denote
design matrices containing evaluations of the fixed effects as $\bm{X}_{\mu}$
and those containing evaluations of the MFPC basis as $\bm{\Psi}$. The $u$-th
row of a design matrix is given as $\bm{x}_{lu}^{\top}$ and
$\bm{\psi}_{u}^{\top}$, respectively. The additional index $(m)$ for a MFPC
basis matrix $\bm{\Psi}_{(m)}$ or a vector $\bm{\psi}_{(m)}$ denotes the
submatrix of $\bm{\Psi}$ or the subvector of $\bm{\psi}$ containing only
columns or elements corresponding to the $m$-th MFPC.

For the survival submodel, $\bm{X}_{\gamma}$ is the $n\times d_{\gamma}$
design matrix for the baseline covariates. The $n\times d_{l}$ design matrix
$\bm{X}_{l}(\bm{T}),l\in\\{\lambda,\alpha_{1},...,\alpha_{K},\mu_{1},...,\mu_{K}\\}$
contains evaluations of the corresponding $\tilde{H}_{l}$ model terms at the
event times $\bm{T}$, where here and in the following bold arguments of
matrices stand for row-wise argument evaluation. The principal component based
random effect design matrix
$\bm{\Psi}^{(k)}(\bm{T})=blockdiag(\bm{\Psi}_{1}^{(k)}(T_{1}),...,\bm{\Psi}_{n}^{(k)}(T_{n}))$
is blockdiagonal with blocks of size $1\times M$.

For the longitudinal submodel, the $N_{k}\times d_{l}$ design matrix
$\bm{X}_{l}(\bm{t}_{k}),l\in\\{\mu_{k},\sigma_{k}|k=1,...,K\\}$ and the
marker-specific $N_{k}\times Mn$ blockdiagonal matrix
$\bm{\Psi}^{(k)}(\bm{t}_{k})=blockdiag(\bm{\Psi}_{1}^{(k)}(\bm{t}_{1k}),...,\bm{\Psi}_{n}^{(k)}(\bm{t}_{nk}))$
contain evaluations of the corresponding $\tilde{H}_{l}$ model terms and the
MFPC based functional random effects at the $N_{k}=\sum_{i}n_{ik}$ repeated
measurement times
$\bm{t}_{k}=(\bm{t}_{1k}^{\top},...,\bm{t}_{nk}^{\top})^{\top}$ of individual
and outcome specific observation times
$\bm{t}_{ik}=(t_{ik1},...,t_{ikn_{ik}})^{\top}$, respectively. We also define
the fixed effects design matrix of the full multivariate longitudinal submodel
$\bm{X}_{\mu}(\bm{t})=blockdiag(\bm{X}_{\mu_{1}}(\bm{t}_{1}),...,\bm{X}_{\mu_{K}}(\bm{t}_{K}))$
comprising all longitudinal observations at respective time-points
$\bm{t}=(\bm{t}^{\top}_{1},...,\bm{t}_{K}^{\top})^{\top}$. The corresponding
$N\times Mn$ random effect matrix $\bm{\Psi}(\bm{t})$ stacks the marker-
specific $N_{k}\times Mn$ matrices. Analogously to $\bm{X}_{\mu}(\bm{t})$, the
design matrix for the variance of the Gaussian noise
$\bm{X}_{\sigma}(\bm{t})=blockdiag(\bm{X}_{\sigma_{1}}(\bm{t}_{1}),...,\bm{X}_{\sigma_{K}}(\bm{t}_{K}))$
is blockdiagonal with $N_{k}\times d_{\sigma_{k}}$ matrix
$\bm{X}_{\sigma_{k}}(\bm{t}_{k})$ and we define
$\bm{R}=blockdiag(\bm{R}_{1},...,\bm{R}_{K})$ with
$\bm{R}_{k}=diag(\exp(\bm{X}_{\sigma_{k}}(\bm{t}_{k})\bm{\beta}_{\sigma_{k}})^{2})$.
The corresponding parameter vectors
$\bm{\beta}_{\mu}=(\bm{\beta}_{\mu_{1}}^{\top},...,\bm{\beta}_{\mu_{K}}^{\top})^{\top}$
and
$\bm{\beta}_{\sigma}=(\bm{\beta}_{\sigma_{1}}^{\top},...,\bm{\beta}_{\sigma_{K}}^{\top})^{\top}$
are stacked over longitudinal outcomes, while the random score vector
$\bm{\rho}=(\bm{\rho}_{1}^{\top},...,\bm{\rho}_{n}^{\top})^{\top}$ is stacked
over subjects.

The log-likelihood of the MJM can be written as

$\displaystyle\ell(\bm{\theta}\mid\bm{T},\bm{\delta},\bm{y})=$
$\displaystyle-\frac{N}{2}\log(2\pi)-\bm{1}_{N}^{\top}\bm{X}_{\sigma}(\bm{t})\bm{\beta}_{\sigma}-\frac{1}{2}\left(\bm{y}-\bm{X}_{\mu}(\bm{t})\bm{\beta}_{\mu}-\bm{\Psi}\left(\bm{t})\bm{\rho}\right)^{\top}\bm{R}^{-1}(\bm{y}-\bm{X}_{\mu}(\bm{t})\bm{\beta}_{\mu}-\bm{\Psi}(\bm{t})\bm{\rho}\right)$
$\displaystyle+\bm{\delta}^{\top}\bigg{(}\bm{X}_{\lambda}(\bm{T})\bm{\beta}_{\lambda}+\bm{X}_{\gamma}\bm{\beta}_{\gamma}+\sum_{k=1}^{K}\big{(}\bm{X}_{\alpha_{k}}(\bm{T})\bm{\beta}_{\alpha_{k}}\big{)}\cdot\big{(}\bm{X}_{\mu_{k}}(\bm{T})\bm{\beta}_{\mu_{k}}+\bm{\Psi}^{(k)}(\bm{T})\bm{\rho}\big{)}\bigg{)}$
$\displaystyle-\sum_{i=1}^{n}\exp(\bm{x}_{\gamma
i}^{\top}\bm{\beta}_{\gamma})\int_{s=0}^{T_{i}}\exp\bigg{(}\bm{x}_{\lambda
i}(s)^{\top}\bm{\beta}_{\lambda}+\sum_{k=1}^{K}\big{(}\bm{x}_{\alpha_{k}i}(s)^{\top}\bm{\beta}_{\alpha_{k}}\big{)}\big{(}\bm{x}_{\mu_{k}i}(s)^{\top}\bm{\beta}_{\mu_{k}}+\bm{\psi}_{i}^{(k)}(s)^{\top}\bm{\rho}_{i}\big{)}\bigg{)}ds$

with $\cdot$ denoting elementwise multiplication of vectors. For notational
ease, write $w_{i}(s)=\exp\big{(}\bm{x}_{\lambda
i}(s)^{\top}\bm{\beta}_{\lambda}+\sum_{k=1}^{K}\big{(}\bm{x}_{\alpha_{k}i}(s)^{\top}\bm{\beta}_{\alpha_{k}}\big{)}\cdot\big{(}\bm{x}_{\mu_{k}i}(s)^{\top}\bm{\beta}_{\mu_{k}}+\bm{\psi}_{i}^{(k)}(s)^{\top}\bm{\rho}_{i}\big{)}\big{)}$.

Then the score vectors of the log-likelihood function are given as

$\displaystyle\bm{s}(\bm{\beta}_{\mu_{k}})=\frac{\partial\ell}{\partial\bm{\beta}_{\mu_{k}}}$
$\displaystyle=\bm{X}_{\mu_{k}}(\bm{t}_{k})^{\top}\bm{R}_{k}^{-1}\left(\bm{y}_{k}-\bm{X}_{\mu_{k}}(\bm{t}_{k})\bm{\beta}_{\mu_{k}}-\bm{\Psi}^{(k)}(\bm{t}_{k})\bm{\rho}\right)+\bm{X}_{\mu_{k}}(\bm{T})^{\top}diag(\bm{\delta})\bm{X}_{\alpha_{k}}(\bm{T})\bm{\beta}_{\alpha_{k}}$
$\displaystyle\quad-\sum_{i=1}^{n}\exp(\bm{x}_{\gamma
i}^{\top}\bm{\beta}_{\gamma})\int_{0}^{T_{i}}w_{i}(s)\bm{x}_{\alpha_{k}i}(s)^{\top}\bm{\beta}_{\alpha_{k}}\cdot\bm{x}_{\mu_{k}i}(s)ds$
$\displaystyle\bm{s}(\bm{\rho}_{(m)})=\frac{\partial\ell}{\partial\bm{\rho}_{(m)}}$
$\displaystyle=\bm{\Psi}_{(m)}(\bm{t})^{\top}\bm{R}^{-1}\left(\bm{y}-\bm{X}_{\mu}(\bm{t})\bm{\beta}_{\mu}-\bm{\Psi}(\bm{t})\bm{\rho}\right)+\sum_{k=1}^{K}\bm{\Psi}_{(m)}^{(k)}(\bm{T})diag(\bm{\delta})\bm{X}_{\alpha_{k}}(\bm{T})\bm{\beta}_{\alpha_{k}}$
$\displaystyle\quad-\sum_{i=1}^{n}\exp(\bm{x}_{\gamma
i}^{\top}\bm{\beta}_{\gamma
i})\int_{0}^{T_{i}}w_{i}(s)\bigg{(}\sum_{k=1}^{K}\big{(}\bm{x}_{\alpha_{k}i}(s)^{\top}\bm{\beta}_{\alpha_{k}}\big{)}\cdot\bm{\psi}_{(m)i}^{(k)}(s)\bigg{)}ds$
$\displaystyle\bm{s}(\bm{\beta}_{\alpha_{k}})=\frac{\partial\ell}{\partial\bm{\beta}_{\alpha_{k}}}$
$\displaystyle=\bm{X}_{\alpha_{k}}(\bm{T})^{\top}diag(\bm{\delta})\left(\bm{X}_{\mu_{k}}(\bm{T})\bm{\beta}_{\mu_{k}}+\bm{\Psi}^{(k)}(\bm{T})\bm{\rho}\right)$
$\displaystyle\quad-\sum_{i=1}^{n}\exp(\bm{x}_{\gamma
i}^{\top}\bm{\beta}_{\gamma
i})\int_{0}^{T_{i}}w_{i}(s)\big{(}\bm{x}_{\mu_{k}i}(s)^{\top}\bm{\beta}_{\mu_{k}}+\bm{\psi}_{i}^{(k)}(s)^{\top}\bm{\rho}_{i}\big{)}\cdot\bm{x}_{\alpha_{k}i}(s)ds$
$\displaystyle\bm{s}(\bm{\beta}_{\lambda})=\frac{\partial\ell}{\partial\bm{\beta}_{\lambda}}$
$\displaystyle=\bm{X}_{\lambda}(\bm{T})^{\top}\bm{\delta}-\sum_{i=1}^{n}\exp(\bm{x}_{\gamma_{i}}^{\top}\bm{\beta}_{\gamma})\int_{0}^{T_{i}}w_{i}(s)\bm{x}_{\lambda
i}(s)ds$
$\displaystyle\bm{s}(\bm{\beta}_{\gamma})=\frac{\partial\ell}{\partial\bm{\beta}_{\gamma}}$
$\displaystyle=\bm{X}_{\gamma}^{\top}\bm{\delta}-\sum_{i=1}^{n}\exp(\bm{x}_{\gamma_{i}}^{\top}\bm{\beta}_{\gamma})\bm{x}_{\gamma_{i}}\int_{0}^{T_{i}}w_{i}(s)ds$
$\displaystyle\bm{s}(\bm{\beta}_{\sigma})=\frac{\partial\ell}{\partial\bm{\beta}_{\sigma}}$
$\displaystyle=-\bm{X}_{\sigma}(\bm{t})^{\top}\bm{1}_{N}+\bigg{(}\bm{X}_{\sigma}(\bm{t})\odot\big{(}\bm{y}-\bm{X}_{\mu}(\bm{t})\bm{\beta}_{\mu}-\bm{\Psi}(\bm{t})\bm{\rho}\big{)}\bigg{)}^{\top}\bm{R}^{-1}\big{(}\bm{y}-\bm{X}_{\mu}(\bm{t})\bm{\beta}_{\mu}-\bm{\Psi}(\bm{t})\bm{\rho}\big{)}$

with row tensor product $\odot$ of a $p\times a$ matrix $A$ and a $p\times b$
matrix $B$ defined as the $p\times ab$ matrix $A\odot
B=(A\otimes\bm{1}_{b}^{\top})\cdot(\bm{1}_{a}^{\top}\otimes B)$ and $\otimes$
the Kronecker product.

The corresponding Hessians are

$\displaystyle\bm{\mathcal{H}}(\bm{\beta}_{\mu_{k}})$
$\displaystyle=\frac{\partial\ell}{\partial\bm{\beta}_{\mu_{k}}\partial\bm{\beta}_{\mu_{k}}^{\top}}=-\bm{X}_{\mu_{k}}(\bm{t}_{k})^{\top}\bm{R}_{k}^{-1}\bm{X}_{\mu_{k}}(\bm{t}_{k})-\sum_{i=1}^{n}\exp(\bm{x}_{\gamma
i}^{\top}\bm{\beta}_{\gamma
i})\int_{0}^{T_{i}}w_{i}(s)\big{(}\bm{x}_{\alpha_{k}i}(s)^{\top}\bm{\beta}_{\alpha_{k}}\big{)}^{2}\cdot\bm{x}_{\mu_{k}i}(s)\bm{x}_{\mu_{k}i}(s)^{\top}ds$
$\displaystyle\bm{\mathcal{H}}(\bm{\rho}_{(m)})$
$\displaystyle=\frac{\partial\ell}{\partial\bm{\rho}_{(m)}\partial\bm{\rho}_{(m)}^{\top}}=-\bm{\Psi}_{(m)}(\bm{t})^{\top}\bm{R}^{-1}\bm{\Psi}_{(m)}(\bm{t})$
$\displaystyle\qquad\qquad\qquad\quad-\sum_{i=1}^{n}\exp(\bm{x}_{\gamma
i}^{\top}\bm{\beta}_{\gamma
i})\int_{0}^{T_{i}}w_{i}(s)\bigg{(}\sum_{k=1}^{K}\big{(}\bm{x}_{\alpha_{k}i}(s)^{\top}\bm{\beta}_{\alpha_{k}}\big{)}\cdot\bm{\psi}_{(m)i}^{(k)}(s)\bigg{)}\bigg{(}\sum_{k=1}^{K}\big{(}\bm{x}_{\alpha_{k}i}(s)^{\top}\bm{\beta}_{\alpha_{k}}\big{)}\cdot\bm{\psi}_{(m)i}^{(k)}(s)\bigg{)}^{\top}ds$
$\displaystyle\bm{\mathcal{H}}(\bm{\beta}_{\alpha})$
$\displaystyle=\frac{\partial\ell}{\partial\bm{\beta}_{\alpha_{k}}\partial\bm{\beta}_{\alpha_{k}}^{\top}}=-\sum_{i=1}^{n}\exp(\bm{x}_{\gamma
i}^{\top}\bm{\beta}_{\gamma
i})\int_{0}^{T_{i}}w_{i}(s)\big{(}\bm{x}_{\mu_{k}i}(s)^{\top}\bm{\beta}_{\mu_{k}}+\bm{\psi}_{i}^{(k)}(s)^{\top}\bm{\rho}_{i}\big{)}^{2}\cdot\bm{x}_{\alpha_{k}i}(s)\bm{x}_{\alpha_{k}i}(s)^{\top}ds$
$\displaystyle\bm{\mathcal{H}}(\bm{\beta}_{\lambda})$
$\displaystyle=\frac{\partial\ell}{\partial\bm{\beta}_{\lambda}\partial\bm{\beta}_{\lambda}^{\top}}=-\sum_{i=1}^{n}\exp(\bm{x}_{\gamma_{i}}^{\top}\bm{\beta}_{\gamma})\int_{0}^{T_{i}}w_{i}(s)\bm{x}_{\lambda
i}(s)\bm{x}_{\lambda i}(s)^{\top}ds$
$\displaystyle\bm{\mathcal{H}}(\bm{\beta}_{\gamma})$
$\displaystyle=\frac{\partial\ell}{\partial\bm{\beta}_{\gamma}\partial\bm{\beta}_{\gamma}^{\top}}=-\sum_{i=1}^{n}\exp(\bm{x}_{\gamma_{i}}^{\top}\bm{\beta}_{\gamma})\bm{x}_{\gamma_{i}}\bm{x}_{\gamma_{i}}^{\top}\int_{0}^{T_{i}}w_{i}(s)ds$
$\displaystyle\bm{\mathcal{H}}(\bm{\beta}_{\sigma})$
$\displaystyle=\frac{\partial\ell}{\partial\bm{\beta}_{\sigma}\partial\bm{\beta}_{\sigma}^{\top}}=-2\bigg{(}\bm{X}_{\sigma}(\bm{t})\odot\big{(}\bm{y}-\bm{X}_{\mu}(\bm{t})\bm{\beta}_{\mu}-\bm{\Psi}(\bm{t})\bm{\rho}\big{)}\bigg{)}^{\top}\bm{R}^{-1}\bigg{(}\bm{X}_{\sigma}(\bm{t})\odot\big{(}\bm{y}-\bm{X}_{\mu}(\bm{t})\bm{\beta}_{\mu}-\bm{\Psi}(\bm{t})\bm{\rho}\big{)}\bigg{)}.$

##  Appendix B Simulation

We present details for the two simulation scenarios such as the used parameter
specifications and computational information. The code to reproduce the
simulation is provided in the accompanying R package MJMbamlss.

###  B.1 Description Scenario I

The data generating process of Scenario I is based on the following additive
predictors:

$\displaystyle\eta_{\lambda i}(t)$ $\displaystyle=1.37\cdot t^{0.37}$
$\displaystyle\eta_{\gamma i}$ $\displaystyle=-1.5+0.48\cdot x_{i}$
$\displaystyle\eta_{\mu_{k}i}(t)$ $\displaystyle=0.2\cdot t-0.25\cdot
x_{i}-0.05\cdot t\cdot x_{i}+b_{k1i}+b_{k2i}\cdot t,$
$\displaystyle\eta_{\alpha_{1}i}(t)$
$\displaystyle=1.5,\,\eta_{\alpha_{2}i}(t)=0.6,\,\eta_{\alpha_{3}i}(t)=0.3,\,\eta_{\alpha_{4}i}(t)=-0.3,\,\eta_{\alpha_{5}i}(t)=-0.6,\,\eta_{\alpha_{6}i}(t)=-1.5$
$\displaystyle\eta_{\sigma_{k}i}(t)$ $\displaystyle=\log(0.06),\quad
k=1,...,6$

where the random effects
$\bm{b}_{i}=(b_{11i},b_{12i},b_{21i},...,b_{61i},b_{62i})^{\top}$ are
independent draws from a 12-dimensional multivariate normal
$\mathcal{N}_{12}(\bm{0},\bm{\Sigma})$ with

$\displaystyle\bm{\Sigma}$
$\displaystyle=\begin{pmatrix}\bm{\Sigma}_{1}&\bm{\Sigma}_{2}^{\top}\\\
\bm{\Sigma}_{2}&\bm{\Sigma}_{3}\end{pmatrix}$

and

$\displaystyle\bm{\Sigma}_{1}=\begin{pmatrix}0.080&-0.070&0.030&0.030&0.022&0.022\\\
-0.070&0.900&0.030&0.030&0.022&0.022\\\
0.030&0.030&0.096&-0.084&0.030&0.030\\\
0.030&0.030&-0.084&1.080&0.030&0.030\\\
0.022&0.022&0.030&0.030&0.112&-0.098\\\
0.022&0.022&0.030&0.030&-0.098&1.260\end{pmatrix},\quad\bm{\Sigma}_{2}=\begin{pmatrix}0.015&0.015&0.022&0.022&0.030&0.030\\\
0.015&0.015&0.022&0.022&0.030&0.030\\\ 0.000&0.000&0.015&0.015&0.022&0.022\\\
0.000&0.000&0.015&0.015&0.022&0.022\\\ 0.000&0.000&0.000&0.000&0.015&0.015\\\
0.000&0.000&0.000&0.000&0.015&0.015\\\ \end{pmatrix},$
$\displaystyle\bm{\Sigma}_{3}=\begin{pmatrix}0.128&-0.112&0.030&0.030&0.022&0.022\\\
-0.112&1.440&0.030&0.030&0.022&0.022\\\
0.030&0.030&0.144&-0.126&0.030&0.030\\\
0.030&0.030&-0.126&1.620&0.030&0.030\\\
0.022&0.022&0.030&0.030&0.160&-0.140\\\
0.022&0.022&0.030&0.030&-0.140&1.800\\\ \end{pmatrix}.$

###  B.2 Description Scenario II

The data generating process of Scenario II is based on the following additive
predictors:

$\displaystyle\eta_{\lambda i}(t)$ $\displaystyle=1.65\cdot t^{0.65}$
$\displaystyle\eta_{\gamma i}$ $\displaystyle=-3+0.3\cdot x_{i}$
$\displaystyle\eta_{\mu_{k}i}(t)$ $\displaystyle=1\cdot t+0.3\cdot
x_{i}+0.3\cdot t\cdot x_{i}+\sum_{m=1}^{6}\rho_{im}\psi_{m}^{(k)}(t)$
$\displaystyle\eta_{\alpha_{k}i}(t)$ $\displaystyle=1.1$
$\displaystyle\eta_{\sigma_{k}i}(t)$ $\displaystyle=\log(0.06),\quad k=1,2.$

The covariance structure induced by $\sum_{m=1}^{6}\rho_{im}\psi_{m}^{(k)}(t)$
is based on a Gaussian process with covariance kernel $\mathcal{K}$ whose
eigenfunctions are constructed solving the eigenanalysis problem

$\displaystyle\bm{B}\bm{Q}\bm{c}=\nu\bm{c}$

with eigenvector $\bm{c}$, eigenvalue $\nu$, block diagonal matrix $\bm{B}$ of
scalar products of the univariate basis functions given in Figure 3, and the
covariance matrix $\bm{Q}$ of the weights of the univariate basis expansion.
The basis functions in Figure 3 reflect changes in the trajectories at the
beginning, middle, or end of the time interval. $\bm{Q}$ is generated as a
random covariance matrix with (rounded to three digits)

$\displaystyle\bm{Q}=\begin{pmatrix}3.124&-0.396&0.892&0.119&-0.668&0.005\\\
-0.396&1.657&0.162&-0.265&-0.495&-0.778\\\
0.892&0.162&1.980&-0.015&-0.906&0.491\\\
0.119&-0.265&-0.015&1.081&0.063&0.728\\\
-0.668&-0.495&-0.906&0.063&0.890&0.243\\\
0.005&-0.778&0.491&0.728&0.243&1.969\end{pmatrix}.$

The resulting eigenvectors $\hat{\bm{c}}$ and eigenvalues given as (rounded to
three digits) $\nu_{1}=1.376$, $\nu_{2}=0.531$, $\nu_{3}=0.149$,
$\nu_{4}=0.101$, $\nu_{5}=0.044$, and $\nu_{6}=0.019$ are then used to
calculate the corresponding multivariate eigenfunctions
$\bm{\psi}_{m},m=1,...,6$ (cf. Happ and Greven, 2018, Section 3.2[14]). The
resulting multivariate eigenfunctions are given in Figure 4 and the scores
$\rho_{im}$ are independent draws from $N(0,\nu_{m})$.

Figure 3: Non-orthogonal univariate basis functions used in simulation
scenario II to generate a multivariate covariance kernel based on a finite KL
decomposition. Figure 4: Multivariate eigenfunctions of the covariance kernel
used in simulation scenario II. Colors correspond to the same multivariate
functional principal components, $\psi_{m}^{(1)}$ and $\psi_{m}^{(2)}$.

###  B.3 Additional Results

Figure 5 shows simulated data and model fits for one simulation iteration in
scenario I. Figure 6 shows the estimated MFPCs in simulation scenario I, where
the true eigenfunctions are linear. Figure 7 shows the true multivariate
eigenfunctions and their estimated counterparts in simulation scenario II.
Table 4 presents the estimation accuracy of the MFPCs in both scenarios as
measured by the norm of the differences between true functions and estimates.

Figure 5: Exemplary model fits to the observed longitudinal data (grey dots) of three random subjects in one simulation iteration for simulation scenario I. The black lines (mu) correspond to the true trajectories. Figure 6: Estimated MFPC basis functions $\hat{\bm{\psi}}_{m}(t)$ for the first three eigenfunctions (red) in simulation scenario I. Estimates are flipped if a reflected version $-\hat{\bm{\psi}}_{m}(t)$ is closer to the true $\bm{\psi}_{m}(t)$. Figure 7: Estimated MFPC basis functions $\hat{\bm{\psi}}_{m}(t)$ for all six eigenfunctions (red) in simulation scenario II. Estimates are flipped if a reflected version $-\hat{\bm{\psi}}_{m}(t)$ is closer to the true $\bm{\psi}_{m}(t)$. Table 4: Mean multivariate norm of the differences between the eigenfunctions and the estimated MFPCs for both simulation scenarios. MFPC | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12
---|---|---|---|---|---|---|---|---|---|---|---|---
Scenario I | 0.57 | 1.12 | 1.07 | 1.12 | 0.96 | 0.52 | 0.65 | 0.90 | 1.02 | 1.04 | 1.00 | 0.67
Scenario II | 0.12 | 0.23 | 0.85 | 0.80 | 0.64 | 0.71 | – | – | – | – | – | –

##  Appendix C Application

Figure 8 shows the seven leading estimated MFPCs for different specifications
in the univariate covariance estimation and different weightings in the
multivariate scalar product. Table 5 presents the estimated parameters in the
PBC application for both the MFPC based (bamlss) and the spline based (JMB)
approaches.

Figure 8: Seven leading estimated MFPC basis functions depending on the number of B-spline basis functions used for the estimation of the univariate covariance functions: four (4B), seven (7B), or ten (10B) basis functions. Included are also the corresponding estimates from the unweighted MFPCA (7Bu). Table 5: Posterior mean and $95\%$ credible interval of the parameter estimates for the two modeling frameworks. | bamlss | JMB
---|---|---
$\beta_{\gamma 1}$ | $-0.20$ |  $[-0.71$; $0.31]$ | $-0.18$ |  $[-0.67$; $0.31]$
$\beta_{\gamma 2}$ | $-0.05$ |  $[-0.40$; $0.31]$ | $-0.02$ |  $[-0.39$; $0.35]$
$\beta_{albumin0}$ | $1.20$ |  $[1.17$; $1.24]$ | $1.36$ |  $[1.28$; $1.44]$
$\beta_{albumin1}$ | $-0.03$ |  $[-0.07$; $0.00]$ | $-0.04$ |  $[-0.08$; $0.00]$
$\beta_{albumin2}$ | $0.00$ |  $[-0.02$; $0.03]$ | $0.01$ |  $[-0.02$; $0.03]$
$\beta_{serBilir0}$ | $1.00$ |  $[0.72$; $1.26]$ | $0.39$ |  $[-0.26$; $1.01]$
$\beta_{serBilir1}$ | $-0.04$ |  $[-0.31$; $0.23]$ | $0.00$ |  $[-0.30$; $0.29]$
$\beta_{serBilir2}$ | $-0.02$ |  $[-0.22$; $0.17]$ | $-0.06$ |  $[-0.25$; $0.14]$
$\beta_{serChol0}$ | $5.61$ |  $[5.50$; $5.72]$ | $5.84$ |  $[5.58$; $6.10]$
$\beta_{serChol1}$ | $0.06$ |  $[-0.05$; $0.17]$ | $0.07$ |  $[-0.05$; $0.19]$
$\beta_{serChol2}$ | $0.05$ |  $[-0.02$; $0.12]$ | $0.04$ |  $[-0.04$; $0.11]$
$\beta_{SGOT0}$ | $4.75$ |  $[4.59$; $4.90]$ | $4.75$ |  $[4.46$; $5.04]$
$\beta_{SGOT1}$ | $0.02$ |  $[-0.13$; $0.18]$ | $-0.02$ |  $[-0.15$; $0.12]$
$\beta_{SGOT2}$ | $-0.06$ |  $[-0.15$; $0.05]$ | $-0.10$ |  $[-0.19$; $-0.01]$
$\beta_{\sigma albumin}$ | $-2.29$ |  $[-2.33$; $-2.26]$ | $-2.33$ |  $[-2.37$; $-2.30]$
$\beta_{\sigma serBilir}$ | $-1.19$ |  $[-1.23$; $-1.15]$ | $-1.21$ |  $[-1.25$; $-1.17]$
$\beta_{\sigma serChol}$ | $-1.71$ |  $[-1.76$; $-1.65]$ | $-1.73$ |  $[-1.79$; $-1.67]$
$\beta_{\sigma SGOT}$ | $-1.33$ |  $[-1.37$; $-1.29]$ | $-1.34$ |  $[-1.38$; $-1.30]$