# A Bayesian joint model of multiple longitudinal and categorical outcomes
with application to multiple myeloma using permutation-based variable
importance

Danilo Alvares1,∗, Jessica K. Barrett1, François Mercier2, Jochen Schulze2,
Sean Yiu3, Felipe Castro2, Spyros Roumpanis2, and Yajing Zhu2,∗∗

1MRC Biostatistics Unit, University of Cambridge, U.K.
2F. Hoffmann-La Roche Ltd, Basel, Switzerland
3Roche Products Ltd, Welwyn Garden City, U.K.

∗email<EMAIL_ADDRESS>
∗∗email<EMAIL_ADDRESS>
###### Abstract

Joint models have proven to be an effective approach for uncovering
potentially hidden connections between various types of outcomes, mainly
continuous, time-to-event, and binary. Typically, longitudinal continuous
outcomes are characterized by linear mixed-effects models, survival outcomes
are described by proportional hazards models, and the link between outcomes
are captured by shared random effects. Other modeling variations include
generalized linear mixed-effects models for longitudinal data and logistic
regression when a binary outcome is present, rather than time until an event
of interest. However, in a clinical research setting, one might be interested
in modeling the physician’s chosen treatment based on the patient’s medical
history in order to identify prognostic factors. In this situation, there are
often multiple treatment options, requiring the use of a multiclass
classification approach. Inspired by this context, we develop a Bayesian joint
model for longitudinal and categorical data. In particular, our motivation
comes from a multiple myeloma study, in which biomarkers display nonlinear
trajectories that are well captured through bi-exponential submodels, where
patient-level information is shared with the categorical submodel. We also
present a variable importance strategy for ranking prognostic factors. We
apply our proposal and a competing model to the multiple myeloma data, compare
the variable importance and inferential results for both models, and
illustrate patient-level interpretations using our joint model.
Keywords: Bayesian inference; Bi-exponential model; Categorical data; Model-
agnostic; M-protein dynamics.

##  1 Introduction

In recent years, there has been a growing interest in joint models that
simultaneously analyze longitudinal and time-to-event outcomes (Ibrahim et
al., 2010; Alsefri et al., 2020). A standard joint model comprises a single
longitudinal outcome characterized by a linear mixed-effects submodel, a
single survival outcome described by a proportional hazards submodel, and an
association structure that connects both submodels (Rizopoulos, 2012). In most
cases, this connection involves random effects and complex integrals
associated with the survival function, but current literature also presents
other alternatives, such as latent classes that stratify heterogeneous
subpopulations (Proust-Lima et al., 2014; Andrinopoulou et al., 2020),
marginal submodels connected by a copula (Ganjali and Baghfalaki, 2015; Cho et
al., 2024), and separate random effects without the need to perform numerical
integration (Alvares and Rubio, 2021).

Many authors have extended the standard joint model. For example, some studies
modeled multiple longitudinal outcomes and/or competing risks (Hickey et al.,
2016; Rué et al., 2017; Rustand et al., 2024), while others incorporated
spatial random effects dependence into joint modeling (Martins et al., 2016,
2017; Rappl et al., 2023). There were also proposals that considered
categorical longitudinal and survival outcomes (Choi et al., 2015; Chi et al.,
2024), two longitudinal submodels, one continuous and one ordinal (Delporte et
al., 2024), and joint models of longitudinal continuous and binary data (Wang
et al., 2000; Horrocks and van Den Heuvel, 2009; Lu et al., 2016; Zhou et al.,
2023). For the latter, to the best of our knowledge, there are currently no
available extensions that accommodate more than two categories. This extension
introduces at least two additional complexities: (i) more parameters to be
estimated, which requires more data, and (ii) imbalanced categories, leading
to poor predictive performance (Kolo, 2010).

To fill this gap in the literature, we developed a Bayesian joint model for
longitudinal and categorical outcomes. Our proposal is motivated by a
retrospective cohort study from a real-world database of patients who were
diagnosed with multiple myeloma in the US between January 2015 and February
2022. In this application, there are two nonlinear longitudinal biomarkers
that are modeled through bi-exponential submodels (Claret et al., 2009; Thai
et al., 2022). The bi-exponential model includes three components (baseline,
growth rate, and decay rate) that summarize the temporal dynamics of each
biomarker and potentially help to explain the physician’s chosen treatment.
Additionally, we considered three treatment options which can be described
using a (nominal) categorical submodel (Agresti, 2013). The longitudinal and
categorical submodels are connected by sharing the patient-level baseline,
growth rate, and decay rate parameters. The main research question here is
what are the key factors supporting the physician’s chosen treatment. This is
motivated by the clinical knowledge that there is a set of factors
contributing to treatment choice, but very rarely are there quantifications of
the relative contribution of each of them.

The adequacy of the bi-exponential submodels is verified through an analysis
of individual weighted residuals (Desmée et al., 2017), while performance of
the categorical submodel is evaluated using metrics for multiclass
classification (Grandini et al., 2020). Moreover, as a complement to the
posterior inference, we present a model-agnostic variable importance approach
with the aim of improving the transparency, trustability, and interpretability
of the results (Casalicchio et al., 2019).

The outline of the work is as follows. Section 2 presents the multiple myeloma
retrospective cohort study that serves as the foundation for the Bayesian
joint model of multiple longitudinal and categorical outcomes, introduced in
Section 3. Section 4 discusses a variable importance strategy for the proposed
joint model. In Section 5, we apply our joint model and a competing model to
the multiple myeloma data, compare classification metrics, variable
importance, and inferential results between both models, and illustrate
patient-level interpretations using our proposal. Finally, in Section 6,
concluding remarks are provided.

##  2 Multiple myeloma data

Multiple myeloma (MM) is a type of blood cancer that affects plasma cells.
These cancerous plasma cells multiply uncontrollably, eventually crowding out
healthy blood cells in the bone marrow, and secrete a monoclonal
protein/paraprotein, also known as M-protein (van de Donk et al., 2021). The
main MM symptoms are bone pain, weakness, fatigue, frequent infections, and
easily broken bones, as well as potentially causing kidney and immune system
problems (Kumar et al., 2017). MM treatments may include targeted therapy e.g.
proteasome inhibitors and immunomodulatory agents, chemotherapy, stem cell
transplant, and immune-based therapies including T-cell engaging antibodies
and CAR-T (Rajkumar and Kumar, 2020). While there is currently no cure for MM,
the advent of new treatments has led to an increase of progression free
survival and overall survival (Gulla and Anderson, 2020).

###  2.1 Data source and eligibility criteria

For illustrating our joint model proposal, we used a Flatiron Health
electronic health record (EHR)-derived de-identified database of patients who
were diagnosed with MM in the US between 1st January 2015 and 28th February
2022 (Birnbaum et al., 2020; Ma et al., 2023). During the study period, the
de-identified data came from over 280 US cancer clinics (approximately 800
sites of care). The majority of patients in the database originate from
community oncology settings; relative community/academic proportions may vary
depending on study cohort. Other eligibility criteria are: (i) at least 18
years of age at MM diagnosis, (ii) no longer than 60 days between initial
diagnosis and first activity (visit or line of therapy initiation), (iii) more
than three months on-treatment before the end of study follow-up, (iv) no
malignancies before MM diagnosis, and (v) patients treated with RVd (a three-
drug therapy composed of Lenalidomide, Bortezomib, and Dexamethasone) in the
first line of therapy, which was oncologist-defined and rule-based, and is one
of the most common regimens for patients with newly diagnosed MM (Punke et
al., 2017). The last criterion ensures the homogeneity of the selected
subpopulation with respect to the initial treatment received. Thus, 1579
patients met all eligibility criteria.

###  2.2 Outcomes and baseline covariates

Since all selected patients are initially treated with RVd, we are interested
in studying what are the key factors contributing to a physician’s chosen
treatment for the second line of therapy. Hence, the primary outcome in this
study is the treatment assigned to the patient after RVd, in which we
considered three options: (I) Carfilzomib-based therapy, (II) Pomalidomide-
based therapy, or (III) other therapy. Carfilzomib-based and Pomalidomide-
based therapies were chosen because they are the most frequently used
treatments. Other therapy add up to 121 regimens, including treatments that
combine Carfilzomib plus Pomalidomide (7%). Supplementary Appendix A shows the
15 most frequent regimens that make up the other therapy. We also have
secondary outcomes, which are longitudinal biomarkers –M-spike and serum free
light chains (FLC)– that measure the concentration of M-protein and the
involved light chain in serum, respectively (g/L for both). Our hypothesis is
that underlying characteristics of biomarker trajectories may support
treatment choice. Baseline covariates collected at initial diagnosis are also
available, such as sex, ethnicity, Eastern Cooperative Oncology Group (ECOG),
and International Staging System (ISS), age, albumin, beta-2-microglobulin
(B2M), creatinine, hemoglobin, lactate dehydrogenase (LDH), lymphocyte count,
neutrophil count, platelet count, immunoglobulin A (IgA), immunoglobulin G
(IgG), immunoglobulin M (IgM). Furthermore, we also considered the duration of
RVd treatment in the first line of therapy as a covariate. Descriptive
summaries for outcomes and baseline covariates can be found in Supplementary
Appendix A.

###  2.3 Treatment of continuous baseline covariates

Most continuous covariates contain missing data (10%-56%). In order to address
this issue, we apply a log transformation to reduce asymmetry, a
standardization (z-score) so that their scales are similar, and a mean
imputation to replace missing values. Supplementary Appendix A shows a
descriptive summary of these covariates before and after being treated.

##  3 Joint modeling framework

We propose a new Bayesian joint model of multiple longitudinal and categorical
outcomes for MM data presented in Section 2. We adopt the shared-parameter
specification (Wu and Carrol, 1988), where the longitudinal and categorical
processes depend on a hidden process defined by random effects. Below we
present the submodels and priors that make up the proposed joint model.

###  3.1 Longitudinal submodels

To describe the patient-level time evolution of M-spike ($k=1$) and FLC
($k=2$) biomarkers, we consider a bi-exponential model (Stein et al., 2008),
which has already been successfully applied in an MM framework (Alvares et
al., 2024) and is described as follows:

$y_{ki}(t)=B_{ki}\Big{[}\exp\left\\{G_{ki}t\right\\}+\exp\left\\{-D_{ki}t\right\\}-1\Big{]}+\epsilon_{ki}(t),$
(1)

where $y_{ki}(t)$ is the $k$th biomarker for patient $i=1,\ldots,n$ at time
$t$; $B_{ki}=\exp\left\\{\theta_{1k}+b_{1ki}\right\\}$,
$G_{ki}=\exp\left\\{\theta_{2k}+b_{2ki}\right\\}$, and
$D_{ki}=\exp\left\\{\theta_{3k}+b_{3ki}\right\\}$ are underlying
characteristics of the longitudinal process and represent baseline, growth
rate, and decay rate, respectively, where
${\bm{\theta}}_{k}=(\theta_{1k},\theta_{2k},\theta_{3k})^{\top}$ are
population parameters and ${\bm{b}}_{ki}=(b_{1ki},b_{2ki},b_{3ki})^{\top}$ are
random effects. Such random effects follow a zero-mean multivariate normal
distribution with unstructured variance-covariance matrix ${\bm{\Omega}}_{k}$.
The residual errors, $\epsilon_{k1}(t),\ldots,\epsilon_{kn}(t)$, are assumed
to be independently normally distributed with mean zero and variance
$\sigma_{k}^{2}$.

###  3.2 Categorical submodel

To model the treatment assigned to each patient in the second line of therapy,
we consider a (nominal) categorical model (Agresti, 2013). Specifically, let
$z_{i}\in\\{1,2,\ldots,J=3\\}$ be a treatment indicator for patient
$i=1,\ldots,n$. Then, assuming a log-odds specification, the categorical
submodel for $z_{i}$ can be written as

$\displaystyle z_{i}$ $\displaystyle\sim\text{Categorical}(J,{\bm{\phi}}_{i})$
(2) $\displaystyle\log\left(\frac{\phi_{ij}}{\phi_{iJ}}\right)$
$\displaystyle={\bm{X}}_{i}^{\top}{\bm{\beta}}_{j}+\sum_{k=1}^{2}\Big{(}\alpha_{1kj}B_{ki}^{\ast}+\alpha_{2kj}G_{ki}^{\ast}+\alpha_{3kj}D_{ki}^{\ast}\Big{)},$

where we contrast each $j=1,\ldots,J-1$ with category $J$. In our application,
$j=1$ and $j=2$ are Carfilzomib-based and Pomalidomide-based therapies,
respectively, and $J=3$ represent other therapies (reference). The probability
vector ${\bm{\phi}}_{i}=(\phi_{ij},\ldots,\phi_{iJ})^{\top}$ is defined on the
simplex $\mathcal{S}^{J-1}=\big{\\{}{\bm{\phi}}_{i}:\phi_{ij}\geq
0,\;\textstyle\sum_{j=1}^{J}\phi_{ij}=1\big{\\}}$; ${\bm{X}}_{i}$ is the
baseline covariate vector with coefficients ${\bm{\beta}}_{j}$;
$B_{ki}^{\ast}=\log(B_{ki})$, $G_{ki}^{\ast}=\log(G_{ki})$, and
$D_{ki}^{\ast}=\log(D_{ki})$ are the baseline, growth rate, and decay rate (in
log scale) of the $k$th biomarker for patient $i$, shared from the
longitudinal submodel (1), where $\alpha_{1kj}$, $\alpha_{2kj}$, and
$\alpha_{3kj}$ are the coefficients that measure the association of such
characteristics with the chosen treatment.

###  3.3 Prior elicitation

We consider default priors for parameters associated with biomarkers ($k=1,2$)
and category outcomes ($j=1,2$), but we acknowledge the possibility of using
alternative priors. In particular, we assume independent and weakly
informative marginal prior distributions. For the bi-exponential submodel (1):
population parameters, $\theta_{1k}$, $\theta_{2k}$, and $\theta_{3k}$, follow
Normal($0,10^{2}$) prior distributions; residual error variance,
$\sigma_{k}^{2}$, follows a half-Cauchy($0,5$) prior distribution; and random
effects variance-covariance matrix, ${\bm{\Omega}}_{k}$, follows an inverse-
Wishart($\mathbf{I}_{3},4$) prior distribution, where $\mathbf{I}_{3}$
represents a $3\times 3$ identity matrix. For the categorical submodel:
regression coefficients, ${\bm{\beta}}_{j}$, and association parameters,
$\alpha_{1kj}$, $\alpha_{2kj}$, and $\alpha_{3kj}$, follow Normal($0,10^{2}$)
prior distributions.

###  3.4 Posterior inference

Let
${\bm{\Psi}}_{1}=\big{(}\text{vec}({\bm{\Omega}}_{1}\big{)}^{\top},\sigma_{1}^{2}\big{)}^{\top}$,
${\bm{\Psi}}_{2}=\big{(}\text{vec}({\bm{\Omega}}_{2})^{\top},\sigma_{2}^{2}\big{)}^{\top}$,
and
${\bm{\Phi}}=\big{(}{\bm{\beta}}_{1}^{\top},{\bm{\beta}}_{2}^{\top},{\bm{\alpha}}_{\cdot
11}^{\top},{\bm{\alpha}}_{\cdot 21}^{\top},{\bm{\alpha}}_{\cdot
12}^{\top},{\bm{\alpha}}_{\cdot 22}^{\top}\big{)}^{\top}$ be parameter vectors
of the submodels associated with M-spike, FLC, and treatment outcomes,
respectively. Also, let ${\bm{\theta}}_{1}$, ${\bm{\theta}}_{2}$,
${\bm{b}}_{1}$, and ${\bm{b}}_{2}$ be the information shared between
submodels, where
${\bm{b}}_{k}=\big{(}{\bm{b}}_{k1}^{\top},\ldots,{\bm{b}}_{kn}^{\top}\big{)}^{\top}$
is the vector of all random effects of the $k$th biomarker. Denoting
${\mathcal{D}}=\\{{\bm{y}}_{1i},{\bm{y}}_{2i},z_{i},{\bm{X}}_{i};\,i=1,\ldots,n\\}$
as the observed data, under the shared-parameter assumption, the joint
posterior distribution of all parameters and random effects is proportionally
expressed as follows:

$\displaystyle\pi({\bm{\Psi}}_{1},{\bm{\Psi}}_{2},{\bm{\Phi}},{\bm{\theta}}_{1},{\bm{\theta}}_{2},{\bm{b}}_{1},{\bm{b}}_{2}\mid{\mathcal{D}})$
$\displaystyle\propto
f_{1}({\bm{y}}_{1}\mid{\bm{\Psi}}_{1},{\bm{\theta}}_{1},{\bm{b}}_{1})f_{2}({\bm{y}}_{2}\mid{\bm{\Psi}}_{2},{\bm{\theta}}_{2},{\bm{b}}_{2})f_{3}({\bm{z}}\mid{\bm{\Phi}},{\bm{\theta}}_{1},{\bm{\theta}}_{2},{\bm{b}}_{1},{\bm{b}}_{2})\times$
(3) $\displaystyle\quad\times
g_{1}({\bm{b}}_{1}\mid{\bm{\Psi}}_{1})g_{2}({\bm{b}}_{2}\mid{\bm{\Psi}}_{2})\pi({\bm{\Psi}}_{1},{\bm{\Psi}}_{2},{\bm{\Phi}},{\bm{\theta}}_{1},{\bm{\theta}}_{2}),$

where $f_{1}(\cdot)$ and $f_{2}(\cdot)$ are Normal density functions derived
from (1); $f_{3}(\cdot)$ is a probability mass function defined as
$\mathds{P}(Z_{i}=j)=\phi_{ij}$ for $j=1,2$, where $\phi_{ij}$ is calculated
as in (2); $g_{1}(\cdot)$ and $g_{2}(\cdot)$ are Normal density functions for
random effects ${\bm{b}}_{1}$ and ${\bm{b}}_{2}$, respectively; and
$\pi(\cdot)$ is the prior distribution specified as in Section 3.3.

##  4 Variable importance

Variable importance (VI) is a widely known tool in machine learning
(Casalicchio et al., 2019), but it has not been widely adopted in statistical
modeling. The concept of VI is to measure the importance of covariates by
their ability to improve a prediction model’s accuracy. VI and other
interpretability techniques can be divided into two groups: model-specific and
model-agnostic. Model-specific techniques are based on assumptions imposed by
a specific model, whereas model-agnostic ones do not rely on a predefined
structure and can be applied to various problems/models without needing
extensive customization or fine-tuning (Chen et al., 2023). In this work, we
focus on the model-agnostic approach due to its versatility and easy
adaptation to different contexts (Ribeiro et al., 2016). In particular, we use
the permutation-based VI strategy proposed by Fisher et al. (2019).

Taking a baseline performance metric (M${}_{\text{base}}$) for the fitted
model with all variables ($x_{1},x_{2},\ldots,x_{P}$), Fisher et al. (2019)
showed that the performance metric degrades when considering a random
permutation of variable $x_{p}$, while the rest of the variables remain
unchanged. Without loss of generality, assuming that a lower metric
performance value is better, the permutation-based importance score for
$x_{p}$ can be computed as
$\text{VI}_{p}=\text{M}_{p}-\text{M}_{\text{base}}$, with $\text{M}_{p}$ being
the performance metric evaluated on the data where only $x_{p}$ is randomly
permuted. This process is repeated for the $p$ variables and then their scores
are ranked, where the highest VI score indicates the most important variable.
Algorithm 1 summarizes this procedure.

Algorithm 1 Permutation-based variable importance algorithm.

1:Model fitted with all variables $x_{1},x_{2},\ldots,x_{P}$.

2:Compute the performance metric for the fitted model (M${}_{\text{base}}$).

3:for $p=1,\ldots,P$ do

4: Permute the values of $x_{p}$.

5: Compute the performance metric on the permuted data (Mp).

6: Compute the variable importance:
$\text{VI}_{p}=\text{M}_{p}-\text{M}_{\text{base}}$.

7:end for

8:Plot ranked $\text{VI}_{1},\text{VI}_{2},\ldots,\text{VI}_{P}$.

The choice of performance metric is an arbitrary decision and may be
conditioned on the type of model analyzed. In this work, we used the widely
applicable information criterion (WAIC) (Watanabe, 2010), but the VI ranking
was identical when applying the leave-one-out cross-validation (Vehtari et
al., 2017). Furthermore, due to randomness, calculating VI more than once
might greatly vary the results. This problem is mitigated by repeating
Algorithm 1 several times and averaging the VI score of each variable. Such a
strategy also allows for evaluating the variability of VI scores, for example,
by using their minimum and maximum values among the repetitions of the same
variable. Finally, it is worth mentioning that after each permutation the
model does not require refitting; we only recompute the performance metric.

##  5 MM data analysis

Our Bayesian joint model was implemented in Stan using the rstan package
version 2.26.23 (Stan Development Team, 2023) from the R language version
4.3.1 (R Core Team, 2023). All codes are publicly available at
www.github.com/daniloalvares/BJM-MBiExp-Categ. We set three Markov chains with
4000 posterior samples after 1000 warm-up iterations. Under this
configuration, posterior samples achieved convergence (Gelman-Rubin statistic,
R-hat $<$ 1.05) and efficiency (effective sample size, neff $>$ 100) (Vehtari
et al., 2021).

As a competing model, we propose the categorical model described in (2)
without considering the latent characteristics shared from the bi-exponential
submodels. This model helps to verify whether our joint model is really
needed.

We evaluate the goodness-of-fit of the joint and competing models, considering
metrics compatible with the nature of each outcome. For bi-exponential
submodels, individual weighted residuals (IWRES) across time by biomarker did
not suggest any model misspecification (Desmée et al., 2017), as the points
are evenly scattered around the horizontal zero-line (see Supplementary
Appendix C). For the categorical submodel and the competing model, class-
weighted classification metrics were applied, such as accuracy, precision,
recall, and F1-score (see formulas in Supplementary Appendix B). Table 1 shows
the percentage of each metric for both models and using a random classifier.

Table 1: Class-weighted classification metrics from the joint model, only
categorical model (no longitudinal information), and a random classifier
(average of 1000 runs).

| Accuracy & Recall∗ | Precision | F1-score
---|---|---|---
Joint model (1)-(2) | 72.13% | 75.19% | 61.23%
Only categorical model | 62.89% | 58.06% | 60.13%
Random classifier | 33.33% | 55.47% | 38.29%

  * ∗Mathematical equivalence in this particular setting (see Supplementary Appendix B).

The joint model proved to be better than the categorical model in all metrics,
with greater and smaller differences for accuracy and F1-score, respectively.
This corroborates our hypothesis that underlying characteristics of biomarker
trajectories help explain the physician’s chosen treatment. In addition, our
proposal was significantly superior to the random classifier in all metrics,
demonstrating that the joint model results are satisfactory. Still, it is
worth mentioning that the treatment decision is a medical choice that may take
into account factors not addressed in this study as well as qualitative
aspects of daily clinical practice that are not available. Finally, we also
compared both models in terms of WAIC and again the joint model (WAIC=2486)
was better than the categorical model (WAIC=2527).

Figure 1 shows the VI ranking for the joint model (1)-(2) based on 50 runs.
Note that the minimum and maximum values of VI scores for each variable
provide a margin of variation that is useful for evaluating potential
importance ties.

Figure 1: Variable importance (VI) ranking from the joint model (1)-(2). Bar
length, lower and upper limits of the horizontal lines represent average,
minimum and maximum VI scores, respectively, based on 50 runs.

We observed a clear higher importance for the top three, growth M-spike,
baseline FLC, and decay FLC, where all of them presented VI scores that are
considerably higher than the others, with minimal overlap between their
margins of variation. This means that, in terms of predictive robustness, the
growth rate of the M-spike trajectory, the initial FLC value and its decay
rate over time are the characteristics that most impact the predicted
treatment for the second line of therapy in a subpopulation initially treated
with RVd. Growth FLC, decay M-spike, and age occupy positions from 4 to 6 in
the ranking, where when analyzing their margins of variation there is no
obvious order between them, but their average scores are higher than the
subsequent ranked variables. The average score of the subsequent variables is
capable of ranking them, but their margins of variation present non-negligible
intersections. Supplementary Appendix D shows the VI ranking using the
categorical model (no longitudinal information). Except for the shared latent
characteristics, the competing model VI ranking produces results very similar
to Figure 1.

Table 2 shows the inferential results in terms of relative risk
(exponentialized coefficient) from the categorical submodel considering
Carfilzomib vs Other and Pomalidomide vs Other, as well as the VI ranking
described in Figure 1.

Table 2: Relative risk (RR) and its 95% credible interval (CI) from the categorical submodel parameters for Carfilzomib vs Other and Pomalidomide vs Other, and variable importance (VI) ranking. Statistically significant variables are shown in bold. Variable | Category | VI ranking | Carfilzomib vs Other | Pomalidomide vs Other
---|---|---|---|---
RR (95% CI) | RR (95% CI)
Sex | Female | 21 | 0.89 (0.68, 1.16) | 1.01 (0.75, 1.36)
Ethnicity | Non-Hisp. Black | 15 | 0.96 (0.67, 1.37) | 0.74 (0.48, 1.12)
Other | 1.18 (0.84, 1.62) | 1.04 (0.73, 1.49)
ECOG | 1 | 10 | 0.81 (0.57, 1.13) | 1.01 (0.68, 1.51)
2+ |  0.59 (0.37, 0.94) |  1.68 (1.08, 2.63)
ISS | Stage II | 8 | 1.23 (0.83, 1.82) | 0.75 (0.49, 1.13)
Stage III | 1.14 (0.68, 1.92) | 1.22 (0.70, 2.10)
Age | – | 6 |  0.76 (0.68, 0.86) |  1.25 (1.07, 1.47)
Albumin | – | 9 | 1.14 (0.96, 1.36) |  0.80 (0.70, 0.91)
B2M | – | 14 | 1.11 (0.88, 1.39) | 0.95 (0.73, 1.23)
Creatine | – | 13 | 1.15 (0.98, 1.34) | 0.95 (0.78, 1.14)
Hemoglobin | – | 18 | 0.92 (0.79, 1.07) | 1.00 (0.84, 1.18)
LDH | – | 12 |  1.33 (1.11, 1.60) | 1.05 (0.85, 1.30)
Lymphocyte | – | 20 | 0.96 (0.83, 1.11) | 0.91 (0.77, 1.07)
Neutrophil | – | 23 | 1.03 (0.87, 1.22) | 1.00 (0.83, 1.21)
Platelet | – | 22 | 0.98 (0.85, 1.15) | 0.99 (0.83, 1.19)
IgA | – | 17 | 1.00 (0.80, 1.25) | 0.95 (0.75, 1.20)
IgG | – | 16 | 1.12 (0.90, 1.39) | 0.97 (0.76, 1.23)
IgM | – | 19 | 0.94 (0.76, 1.14) | 0.86 (0.67, 1.07)
RVd duration | – | 7 |  0.76 (0.65, 0.89) |  1.19 (1.03, 1.37)
Baseline M-spike | – | 11 | 1.13 (0.89, 1.44) | 1.14 (0.87, 1.49)
Growth M-spike | – | 1 | 1.15 (0.80, 1.64) |  1.57 (1.03, 2.36)
Decay M-spike | – | 5 |  0.69 (0.53, 0.90) |  0.72 (0.54, 0.96)
Baseline FLC | – | 2 |  1.41 (1.10, 1.82) |  1.61 (1.24, 2.11)
Growth FLC | – | 4 |  1.67 (1.21, 2.30) |  1.66 (1.18, 2.35)
Decay FLC | – | 3 | 0.81 (0.64, 1.04) |  0.65 (0.51, 0.82)

In Table 2, we highlight ECOG 2+, age and RVd duration, which are significant
and present relative risks in opposite directions for Carfilzomib and
Pomalidomide compared to the reference treatment. For example, for RVd
duration, this means that the relative risk ratio for a one-year increase in
RVd treatment time in the first line of therapy is 1.32 (1/0.76) for being
treated with Other over Carfilzomib, while this ratio is 1.19 in favor of
Pomalidomide when contrasted with Other. Conversely, the significant M-spike
and FLC characteristics are consistent for both comparisons. Indeed,
increasing the M-spike decay rate favors the choice of Other over Carfilzomib
or Pomalidomide. On the other hand, starting RVd treatment with a high FLC
value or increasing the growth FLC rate throughout the first line of therapy
increases the likelihood of choosing Carfilzomib or Pomalidomide over Other
for the second line of therapy. Supplementary Appendix D shows that the
inferential results using the categorical model (no longitudinal information)
are consistent with those in Table 2, except for the M-spike and FLC
characteristics that are not present in the competing model.

Clinically, one may also be interested in identifying significant factors
between Carfilzomib and Pomalidomide. This analysis can be done by inspecting
non-intersecting 95% credible intervals per variable. For example, from Table
2, we can list ECOG 2+, age, albumin, and RVd duration as significant factors
that differentiate the choice between Carfilzomib and Pomalidomide. It is
worth noting that latent characteristics of M-spike and FLC trajectories are
not statistically relevant when contrasting Carfilzomib and Pomalidomide.

Estimated parameters for the bi-exponential M-spike and FLC submodels are
provided in Supplementary Appendix C. Population parameters show similar
shapes for the average trajectories of M-spike and FLC biomarkers, where
M-spike presents an estimated initial value lower than that of FLC, but its
growth and decay rates are higher. Residual error and random effects variances
indicate greater variability for FLC trajectories. In addition, significant
random effects covariances are consistent for both biomarkers, where the
random effects for baseline and growth rate are negatively correlated, while
baseline and decay rate are positively correlated.

In general, VI ranking and inference are consistent with each other (see Table
2). For example, both support that latent characteristics shared from the bi-
exponential submodels are important to explaining treatment choice. Still, it
is worth remembering that VI calculates its scores based on a performance
metric (in this work, WAIC), so the VI ranking prioritizes the magnitude of
the permutation impact on this metric. In other words, VI is not directly
based on statistical significance, so it is expected that both criteria
complement each other in the task of identifying relevant prognostic factors
for choosing treatment.

To illustrate individual-specific interpretations using our joint modeling, we
randomly selected two patients, named A and B, whose baseline covariates and
categorical outcomes are presented in Supplementary Appendix C. Figure 2(a)
shows the M-spike and FLC dynamics of patient A, the fit of these trajectories
using bi-exponential submodels, and the predicted probability of each
treatment option. Initial M-spike was estimated to be close to 25 g/L with a
slight decay over the first year followed by an increase until reaching 30
g/L, while FLC presented values below 20 g/L that were well estimated, but
with a beginning of trajectory with great uncertainty due to the absence of
FLC measurements for this patient over 40 days after starting RVd treatment.
Thereby, the predicted probabilities from the categorical submodel suggested
Pomalidomide (highest posterior mean) as a treatment for the second line of
therapy, which coincides with the physician’s chosen treatment.

(a) Patient A.

(b) Patient B.

Figure 2: Left side: observed values (asterisks), posterior mean trajectory
(solid line), and 95% credible interval (gray shadow) for each biomarker of
patients A and B. Right side: posterior mean (point) and 95% credible interval
(vertical line) for each treatment option in the second line of therapy.

Analogously, Figure 2(b) shows the estimated trajectories for M-spike and FLC
biomarkers of patient B, and the predicted probability of each treatment
option for the second line of therapy. Here, the range of longitudinal
observations is shorter, but we can notice an M-spike trajectory trend similar
to the previous case, but with lower values. On the other hand, the estimated
FLC trajectory of patient B is almost constant at 25 g/L, which is clearly
higher than that of patient A. This time, the categorical submodel estimated a
higher probability for Carfilzomib to be chosen in the second line of therapy,
which was also the physician’s selected treatment.

##  6 Discussion

Motivated by a clinical study on MM disease, we have proposed the first
Bayesian joint model that accommodates multiple nonlinear longitudinal
biomarkers and categorical outcomes. In addition, borrowing ideas from the
world of machine learning, we have revisited Fisher et al. (2019)’s proposal
to rank the importance of variables from the proposed model applied to MM
data.

Evaluation of model fit was performed by residual analysis and multiclass
classification metrics. IWRES showed good suitability of the bi-exponential
submodel to fit the longitudinal trajectories of M-spike and FLC biomarkers,
similar to the results obtained by Alvares et al. (2024). The classification
metrics also proved to be satisfactory for a three-class problem, mainly when
contrasted with the competing model and a random classifier (see Table 1). VI
ranking and inference were concise in pointing out the latent characteristics
of biomarker trajectories as relevant prognostic factors for choosing
treatment. However, it is also important to remember that VI is based on a
predictive performance metric, so some differences between VI ranking and
inferential significance are expected. Still, VI provides useful and
complementary information to the effect sizes, e.g. it proposes a clear
ranking of the contribution of each variable to the treatment decision. It is
worth mentioning that the treatment decision is a medical choice that may take
into account relevant clinical parameters not available in this study, such as
comorbidities, fitness, treatment toxicity, patient preference, among others.
Finally, we illustrated the use of our joint model for two random patients,
showing the longitudinal fit of their biomarkers, and their predicted
probabilities for each treatment option.

In conclusion, our study has filled a gap in the modeling literature, shed
light on topics related to the clinical practice of MM patients, and connected
a popular machine learning strategy to a joint model. As future work, we
envision several extensions of the proposed model. For example, in a high-
dimensional longitudinal biomarker framework, the robustness to identify
significant association parameters is reduced due to potentially strong
correlations among biomarkers, but this problem can be overcome using
regularization techniques (Andrinopoulou and Rizopoulos, 2016; Sun and Basu,
2024). In the presence of heterogeneous subgroups, we recommend adapting our
proposal to joint latent class models (Proust-Lima et al., 2014; Andrinopoulou
et al., 2020). In some scenarios, researchers may be interested in studying
the causal effect of treatment choice, requiring a causal inference approach
(Rizopoulos et al., 2024). Furthermore, the probability for each treatment can
be dynamically updated as new longitudinal patient information becomes
available (Rizopoulos, 2011; Andrinopoulou et al., 2017). Finally, machine
learning algorithms can be proposed for each submodel, such as ensembling
techniques for longitudinal data (Hu and Szymczak, 2023) and boosting methods
for multiclass classification (Tanha et al., 2020).

## Data availability statement

For eligible studies qualified researchers may request access to individual
patient-level clinical data through a data request platform. For up-to-date
details on Roche’s Global Policy on the Sharing of Clinical Information and
how to request access to related clinical study documents, see the website
(https://go.roche.com/data_sharing). Anonymised records for individual
patients across more than one data source external to Roche cannot, and should
not, be linked due to a potential increase in risk of patient re-
identification. The data that support the findings of this study were
originated by and are the property of Flatiron Health, Inc., which has
restrictions prohibiting the authors from making the data set publicly
available. Requests for data sharing by license or by permission for the
specific purpose of replicating results in this manuscript can be submitted to
<EMAIL_ADDRESS>The data are subject to a license
agreement with Flatiron Health to protect patient privacy and ensure
compliance with measures necessary to reduce the risk of re-identification.
For example, the data necessary to replicate the study include numerous
specific dates, including visit dates (i.e., laboratory or examination dates),
treatment start and stop dates, and month of death, as well as laboratory test
results. Other measures to maintain de-identification without contractual
agreements in place are not feasible due to the study question, methods used,
and data elements required.

## Acknowledgements

D.A. and J.K.B. were supported by the U.K. Medical Research Council grant
MC_UU_00002/5 and the collaboration grant jointly funded by Roche and the
University of Cambridge. We thank the following colleagues for helpful
discussions during the research and review stage of this work: Vallari Shah,
Mellissa Williamson, Sarwar Mozumder, Madlaina Breuleux, Pascal Chanu, and
Chris Harbron. For the purpose of open access, the author has applied a
Creative Commons Attribution (CC BY) license to any Author Accepted Manuscript
version arising from this submission.

## References

  * Agresti (2013) Agresti, A. (2013).  Categorical data analysis.  John Wiley & Sons, New Jersey, US, 3rd edition.
  * Alsefri et al. (2020) Alsefri, M., Sudell, M., García-Fiñana, M., and Kolamunnage-Dona, R. (2020).  Bayesian joint modelling of longitudinal and time to event data: A methodological review.  BMC Medical Research Methodology 20, 1–17.
  * Alvares et al. (2024) Alvares, D., Barrett, J. K., Mercier, F., Roumpanis, S., Yiu, S., Castro, F. Schulze, J., and Zhu, Y. (2024).  A bayesian joint model of multiple nonlinear longitudinal and competing risks outcomes for dynamic prediction in multiple myeloma: Joint estimation and corrected two-stage approaches.  arXiv:2405.20418 .
  * Alvares and Rubio (2021) Alvares, D. and Rubio, F. J. (2021).  A tractable bayesian joint model for longitudinal and survival data.  Statistics in Medicine 40, 4213–4229.
  * Andrinopoulou et al. (2020) Andrinopoulou, E. R., Nasserinejad, K., Szczesniak, R., and Rizopoulos, D. (2020).  Integrating latent classes in the Bayesian shared parameter joint model of longitudinal and survival outcomes.  Statistical Methods in Medical Research 29, 3294–3307.
  * Andrinopoulou and Rizopoulos (2016) Andrinopoulou, E. R. and Rizopoulos, D. (2016).  Bayesian shrinkage approach for a joint model of longitudinal and survival outcomes assuming different association structures.  Statistics in Medicine 35, 4813–4823.
  * Andrinopoulou et al. (2017) Andrinopoulou, E. R., Rizopoulos, D., Takkenberg, J. J. M., and Lesaffre, E. (2017).  Combined dynamic predictions using joint models of two longitudinal outcomes and competing risk data.  Statistical Methods in Medical Research 26, 1787–1801.
  * Birnbaum et al. (2020) Birnbaum, B., Nussbaum, N., Seidl-Rathkopf, K., Agrawal, M., Estevez, M., Estola, E., Haimson, J., He, L., Larson, P., and Richardson, P. (2020).  Model-assisted cohort selection with bias analysis for generating large-scale cohorts from the EHR for oncology research.  arXiv:2001.09765 .
  * Casalicchio et al. (2019) Casalicchio, G., Molnar, C., and Bischl, B. (2019).  Visualizing the feature importance for black box models.  In Machine learning and knowledge discovery in databases, pages 655–670, Cham, Switzerland. Springer International Publishing.
  * Chen et al. (2023) Chen, Z., Xiao, F., Guo, F., and Yan, J. (2023).  Interpretable machine learning for building energy management: A state-of-the-art review.  Advances in Applied Energy 9, 1–19.
  * Chi et al. (2024) Chi, M., Wang, X., Song, H., Peng, Y., and Tu, D. (2024).  Joint analysis of longitudinal ordinal categorical item response data and survival times with cure fraction.  Statistics in Biopharmaceutical Research pages 1–11.
  * Cho et al. (2024) Cho, S., Psioda, M. A., and Ibrahim, J. G. (2024).  Bayesian joint modeling of multivariate longitudinal and survival outcomes using Gaussian copulas.  Biostatistics pages 1–16.
  * Choi et al. (2015) Choi, J., Cai, J., Zeng, D., and Olshan, A. F. (2015).  Joint analysis of survival time and longitudinal categorical outcomes.  Statistics in Biosciences 7, 19–47.
  * Claret et al. (2009) Claret, L., Girard, P., Hoff, P. M., van Cutsem, E., Zuideveld, K. P., Jorga, K., Fagerberg, J., and Bruno, R. (2009).  Model-based prediction of phase III overall survival in colorectal cancer on the basis of phase II tumor dynamics.  Journal of Clinical Oncology 27, 4103–4108.
  * Delporte et al. (2024) Delporte, M., Molenberghs, G., Fieuws, S., and Verbeke, G. (2024).  A joint normal-ordinal (probit) model for ordinal and continuous longitudinal data.  Biostatistics pages 1–13.
  * Desmée et al. (2017) Desmée, S., Mentré, F., Veyrat-Follet, C., Sébastien, B., and Guedj, J. (2017).  Using the SAEM algorithm for mechanistic joint models characterizing the relationship between nonlinear PSA kinetics and survival in prostate cancer patients.  Biometrics 73, 305–312.
  * Donders et al. (2006) Donders, A. R. T., van der Heijden, G. J. M. G., Stijnen, T., and Moons, K. G. M. (2006).  Review: A gentle introduction to imputation of missing values.  Journal of Clinical Epidemiology 59, 1087–1091.
  * Fisher et al. (2019) Fisher, A., Rudin, C., and Dominici, F. (2019).  All models are wrong, but many are useful: Learning a variable’s importance by studying an entire class of prediction models simultaneously.  Journal of Machine Learning Research 20, 1–81.
  * Ganjali and Baghfalaki (2015) Ganjali, M. and Baghfalaki, T. (2015).  A copula approach to joint modeling of longitudinal measurements and survival times using Monte Carlo Expectation-Maximization with application to AIDS studies.  Journal of Biopharmaceutical Statistics 25, 1077–1099.
  * Grandini et al. (2020) Grandini, M., Bagli, E., and Visani, G. (2020).  Metrics for multi-class classification: An overview.  arXiv:2008.05756 .
  * Gulla and Anderson (2020) Gulla, A. and Anderson, K. C. (2020).  Multiple myeloma: The (r)evolution of current therapy and a glance into future.  Haematologica 105, 2358–2367.
  * Hickey et al. (2016) Hickey, G. L., Philipson, P., Jorgensen, A., and Kolamunnage-Dona, R. (2016).  Joint modelling of time-to-event and multivariate longitudinal outcomes: Recent developments and issues.  BMC Medical Research Methodology 16, 1–15.
  * Horrocks and van Den Heuvel (2009) Horrocks, J. and van Den Heuvel, M. J. (2009).  Prediction of pregnancy: A joint model for longitudinal and binary data.  Bayesian Analysis 4, 523–538.
  * Hu and Szymczak (2023) Hu, J. and Szymczak, S. (2023).  A review on longitudinal data analysis with random forest.  Briefings in Bioinformatics 24, 1–11.
  * Ibrahim et al. (2010) Ibrahim, J. G., Chu, H., and Chen, L. M. (2010).  Basic concepts and methods for joint models of longitudinal and survival data.  Journal of Clinical Oncology 28, 2796–2801.
  * Kolo (2010) Kolo, B. (2010).  Binary and multiclass classification.  Weatherford Press, Weatherford, OK, USA, 1st edition.
  * Kumar et al. (2017) Kumar, S. K., Rajkumar, V., Kyle, R. A., van Duin, M., Sonneveld, P., Mateos, M. V., Gay, F., and Anderson, K. C. (2017).  Multiple myeloma.  Nature Reviews Disease Primers 3, 1–20.
  * Lu et al. (2016) Lu, X., Huang, Y., and Zhou, R. (2016).  Joint analysis of nonlinear heterogeneous longitudinal data and binary outcome: An application to AIDS clinical studies.  Journal of Applied Statistics 43, 2713–2728.
  * Ma et al. (2023) Ma, X., Long, L., Moon, S., Adamson, B. J. S., and Baxi, S. S. (2023).  Comparison of population characteristics in real-world clinical oncology databases in the US: Flatiron Health, SEER, and NPCR.  medRxiv:10.1101/2020.03.16.20037143v3 .
  * Martins et al. (2016) Martins, R., Silva, G. L., and Andreozzi, V. (2016).  Bayesian joint modeling of longitudinal and spatial survival AIDS data.  Statistics in Medicine 35, 3368–3384.
  * Martins et al. (2017) Martins, R., Silva, G. L., and Andreozzi, V. (2017).  Joint analysis of longitudinal and survival AIDS data with a spatial fraction of long-term survivors: A Bayesian approach.  Biometrical Journal 59, 1166–1183.
  * Proust-Lima et al. (2014) Proust-Lima, C., Séne, M., Taylor, J. M. G., and Jacqmin-Gadda, H. (2014).  Joint latent class models for longitudinal and time-to-event data: A review.  Statistical Methods in Medical Research 23, 74–90.
  * Punke et al. (2017) Punke, A. P., Waddell, J. A., and Solimando, D. A. (2017).  Lenalidomide, Bortezomib, and Dexamethasone (RVD) regimen for multiple myeloma.  Hospital Pharmacy 52, 27–32.
  * R Core Team (2023) R Core Team (2023).  R: A language and environment for statistical computing.  R Foundation for Statistical Computing, https://www.R-project.org/.
  * Rajkumar and Kumar (2020) Rajkumar, S. V. and Kumar, S. (2020).  Multiple myeloma current treatment algorithms.  Blood Cancer Journal 10, 1–10.
  * Rappl et al. (2023) Rappl, A., Kneib, T., Lang, S., and Bergherr, E. (2023).  Spatial joint models through Bayesian structured piecewise additive joint modelling for longitudinal and time-to-event data.  Statistics and Computing 33, 1–16.
  * Ribeiro et al. (2016) Ribeiro, M. T., Singh, S., and Guestrin, C. (2016).  Model-agnostic interpretability of machine learning.  arXiv:1606.05386 .
  * Rizopoulos (2011) Rizopoulos, D. (2011).  Dynamic predictions and prospective accuracy in joint models for longitudinal and time-to-event data.  Biometrics 67, 819–829.
  * Rizopoulos (2012) Rizopoulos, D. (2012).  Joint models for longitudinal and time-to-event data: With applications in R.  Chapman & Hall/CRC, Boca Raton, FL, USA, 1st edition.
  * Rizopoulos et al. (2024) Rizopoulos, D., Taylor, J. M. G., Papageorgiou, G., and Morgan, T. M. (2024).  Using joint models for longitudinal and time-to-event data to investigate the causal effect of salvage therapy after prostatectomy.  Statistical Methods in Medical Research 33, 894–908.
  * Rué et al. (2017) Rué, M., Andrinopoulou, E. R., Alvares, D., Armero, C., Forte, A., and Blanch, L. (2017).  Bayesian joint modeling of bivariate longitudinal and competing risks data: An application to study patient-ventilator asynchronies in critical care patients.  Biometrical Journal 59, 1184–1203.
  * Rustand et al. (2024) Rustand, D., van Niekerk, J., Krainski, E. T., Rue, H., and Proust-Lima, C. (2024).  Fast and flexible inference for jointmodels of multivariate longitudinal and survival data using integrated nested Laplace approximations.  Biostatistics 25, 429–448.
  * Stan Development Team (2023) Stan Development Team (2023).  RStan: The R interface to Stan.  Stan, http://mc-stan.org/.
  * Stein et al. (2008) Stein, W. D., Figg, W. D., Dahut, W., Stein, A. D., Hoshen, M. B., Price, D., Bates, S. E., and Fojo, T. (2008).  Tumor growth rates derived from data for patients in a clinical trial correlate strongly with patient survival: A novel strategy for evaluation of clinical trial data.  The Oncologist 13, 1046–1054.
  * Sun and Basu (2024) Sun, J. and Basu, S. (2024).  Penalized joint models of high-dimensional longitudinal biomarkers and a survival outcome.  The Annals of Applied Statistics 18, 1490–1505.
  * Tanha et al. (2020) Tanha, J., Abdi, Y., Samadi, N., Razzaghi, N., and Asadpour, M. (2020).  Boosting methods for multi-class imbalanced data classification: An experimental review.  Journal of Big Data 7, 1–47.
  * Thai et al. (2022) Thai, H. T., Gaudel, N., Cerou, M., Ayral, G., Fau, J. B., Sebastien, B., van de Velde, H., Semiond, D., and Veyrat-Follet, C. (2022).  Joint modelling and simulation of M-protein dynamics and progression-free survival for alternative isatuximab dosing with pomalidomide/dexamethasone.  British Journal of Clinical Pharmacology 88, 2052–2064.
  * van de Donk et al. (2021) van de Donk, N. W. C. J., Pawlyn, C., and Yong, K. L. (2021).  Multiple myeloma.  The Lancet 397, 410–427.
  * Vehtari et al. (2017) Vehtari, A., Gelman, A., and Gabry, J. (2017).  Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC.  Statistics and Computing 27, 1413–1432.
  * Vehtari et al. (2021) Vehtari, A., Gelman, A., Simpson, D., Carpenter, B., and Bürkner, P. C. (2021).  Rank-normalization, folding, and localization: An improved $\hat{R}$ for assessing convergence of MCMC (with discussion).  Bayesian Analysis 16, 667–718.
  * Wang et al. (2000) Wang, C. Y., Wang, N., and Wang, S. (2000).  Regression analysis when covariates are regression parameters of a random effects model for observed longitudinal measurements.  Biometrics 56, 487–495.
  * Watanabe (2010) Watanabe, S. (2010).  Asymptotic equivalence of Bayes cross validation and widely applicable information criterion in singular learning theory.  Journal of Machine Learning Research 11, 3571–3594.
  * Wu and Carrol (1988) Wu, M. C. and Carrol, R. J. (1988).  Estimation and comparison of changes in the presence of informative right censoring by modeling the censoring problem.  Biometrics 44, 175–188.
  * Zhou et al. (2023) Zhou, G. C., Song, S., and Szczesniak, R. D. (2023).  Multilevel joint model of longitudinal continuous and binary outcomes for hierarchically structured data.  Statistics in Medicine 42, 2914–2927.

SUPPLEMENTARY MATERIAL

## Appendix A. Descriptive summaries for MM data

Table 1: The 15 most frequent regimens that make up the Other therapy. Regimen | Number of cases (%)
---|---
Dexamethasone + Lenalidomide | 123 (11)
Bortezomib + Dexamethasone + Lenalidomide | 123 (11)
Daratumumab + Dexamethasone + Lenalidomide | 75 (7)
Carfilzomib + Dexamethasone + Pomalidomide | 67 (6)
Lenalidomide | 56 (5)
Bortezomib + Dexamethasone | 54 (5)
Bortezomib + Daratumumab + Dexamethasone | 52 (5)
Dexamethasone + Ixazomib + Lenalidomide | 51 (5)
Bortezomib + Cyclophosphamide + Dexamethasone | 47 (4)
Dexamethasone | 39 (3)
Bortezomib | 28 (2)
Clinical Study Drug | 28 (2)
Daratumumab + Dexamethasone | 28 (2)
Bortezomib + Lenalidomide | 24 (2)
Bortezomib + Daratumumab + Dexamethasone + Lenalidomide | 24 (2)
Table 2: Summary of the distribution of the number of M-spike and FLC measurements per patient. Summary | M-spike | FLC
---|---|---
Min | 1 | 1
Quartile 1 | 2 | 3
Median | 4 | 5
Quartile 3 | 6 | 8
Max | 97 | 72
Table 3: Summary of categorical outcomes and covariates. Variable | Number of cases (%)
---|---
Treatment (outcome) |
Carfilzomib | 253 (16)
Pomalidomide | 194 (12)
Other | 1132 (72)
Sex |
Male | 861 (55)
Female | 718 (45)
Ethnicity |
Non-Hispanic white | 913 (58)
Non-Hispanic black | 246 (16)
Other | 291 (18)
Not reported | 129 (8)
ECOG (Eastern Cooperative Oncology Group) |
0 | 402 (25)
1 | 416 (26)
2+ | 199 (13)
Not reported | 562 (36)
ISS (International Staging System) |
Stage I | 395 (25)
Stage II | 377 (24)
Stage III | 326 (21)
Not reported | 481 (30)
Table 4: Summary of continuous variables in their original scales (Initial)
and after log transformation11, standardization (z-score), and imputation22
(Final). NA represents the number of missing observations before imputation.

Variable | Data | Mean | SD3 | Median | Min | Max | NA (%)
---|---|---|---|---|---|---|---
Age4 | Initial | 65.29 | 10.26 | 67.00 | 30.00 | 84.00 | –
(years) | Final | 0.00 | 1.00 | 0.23 | -4.47 | 1.55 | –
Albumin | Initial | 36.73 | 9.79 | 38.00 | 1.00 | 201.00 | 371 (23)
(serum, g/L) | Final | 0.00 | 0.87 | 0.04 | -7.43 | 3.80 | –
B2M (beta-2-microglobulin) | Initial | 7.67 | 71.22 | 3.85 | 0.20 | 2000.00 | 790 (50)
(serum, mg/L) | Final | 0.00 | 0.71 | 0.00 | -4.04 | 9.43 | –
Creatinine | Initial | 1.29 | 0.94 | 1.04 | 0.40 | 10.70 | 373 (24)
(serum, mg/dL) | Final | 0.00 | 0.87 | 0.00 | -2.10 | 5.03 | –
Hemoglobin | Initial | 10.83 | 2.20 | 10.80 | 4.30 | 17.70 | 157 (10)
(g/dL) | Final | 0.00 | 0.95 | 0.00 | -4.28 | 2.45 | –
LDH (lactate dehydrogenase) | Initial | 224.25 | 170.17 | 175.00 | 50.00 | 2377.00 | 891 (56)
(serum, U/L) | Final | 0.00 | 0.66 | 0.00 | -2.61 | 4.90 | –
Lymphocyte | Initial | 1.88 | 1.31 | 1.65 | 0.19 | 21.33 | 386 (24)
(count, $\times 10^{9}$/L) | Final | 0.00 | 0.87 | 0.00 | -3.73 | 5.21 | –
Neutrophil | Initial | 3.87 | 2.35 | 3.40 | 0.41 | 29.30 | 582 (37)
(count, $\times 10^{9}$/L) | Final | 0.00 | 0.79 | 0.00 | -3.57 | 4.01 | –
Platelet | Initial | 227.08 | 89.61 | 218.00 | 29.00 | 834.00 | 412 (26)
(count, $\times 10^{9}$/L) | Final | 0.00 | 0.86 | 0.00 | -4.88 | 3.40 | –
IgA (immunoglobulin A) | Initial | 8.09 | 16.46 | 0.54 | 0.00 | 136.00 | 717 (45)
(serum, g/L) | Final | 0.00 | 0.74 | 0.00 | -1.41 | 2.55 | –
IgG (immunoglobulin G) | Initial | 27.16 | 26.58 | 17.05 | 0.40 | 139.01 | 677 (43)
(serum, g/L) | Final | 0.00 | 0.76 | 0.00 | -2.88 | 1.88 | –
IgM (immunoglobulin M) | Initial | 0.60 | 4.55 | 0.22 | 0.00 | 97.43 | 857 (54)
(serum, g/L) | Final | 0.00 | 0.68 | 0.00 | -1.84 | 8.59 | –
RVd duration (years) | Initial | 1.11 | 0.97 | 0.79 | 0.08 | 5.82 | –

  * 1

$\log(x+0.1)$ to avoid numerical problems when $x=0$.

  * 2

Mean-value imputation strategy for missing values (Donders et al., 2006),
i.e., we replace them with zeros (after standardization, the mean of each
variable is zero).

  * 3

After standardization, the standard deviation is equal to 1, but the
imputation concentrates more values at zero and consequently reduces such
standard deviation.

  * 4

Patients with a BirthYear of [Data Cutoff Year - 85] or earlier may have an
adjusted BirthYear in Flatiron Health databases due to patient de-
identification requirements.

Figure 1: Correlation of continuous covariates after log transformation,
standardization, and imputation.

## Appendix B. Goodness-of-fit criteria

### Individual weighted residuals for longitudinal submodels:

$\text{IWRES}_{i}(t)=\frac{y_{i}(t)-\hat{y}_{i}(t)}{\hat{\sigma}}$

### Class-weighted classification metrics for categorical submodels:

|  | Observed |
---|---|---|---
|  | Treatment 1 | Treatment 2 | Treatment 3 | Total
| Treatment 1 | $\text{TP}_{1}$ | $\cdot$ | $\cdot$ | $n_{1(\text{pred})}$
Predicted | Treatment 2 | $\cdot$ | $\text{TP}_{2}$ | $\cdot$ | $n_{2(\text{pred})}$
| Treatment 3 | $\cdot$ | $\cdot$ | $\text{TP}_{3}$ | $n_{3(\text{pred})}$
| Total | $n_{1(\text{obs})}$ | $n_{2(\text{obs})}$ | $n_{3(\text{obs})}$ | $N$
| Weight | $w_{1}=n_{1(\text{obs})}/N$ | $w_{2}=n_{2(\text{obs})}/N$ | $w_{3}=n_{3(\text{obs})}/N$ |

  * •

Accuracy = $\displaystyle\frac{1}{N}\sum_{c=1}^{3}\text{TP}_{c}$

  * •

Precision =
$\displaystyle\sum_{c=1}^{3}w_{c}\text{Precision}_{c}=\sum_{c=1}^{3}w_{c}\left(\frac{\text{TP}_{c}}{n_{c(\text{pred})}}\right)$

  * •

Recall =
$\displaystyle\sum_{c=1}^{3}w_{c}\text{Recall}_{c}=\sum_{c=1}^{3}w_{c}\left(\frac{\text{TP}_{c}}{n_{c(\text{obs})}}\right)=\frac{1}{N}\sum_{c=1}^{3}\text{TP}_{c}$
= Accuracy

  * •

F1-score =
$\displaystyle\sum_{c=1}^{3}w_{c}\left(\frac{2}{\text{Precision}^{-1}_{c}+\text{Recall}^{-1}_{c}}\right)=\sum_{c=1}^{3}w_{c}\left(\frac{2}{\frac{n_{c(\text{pred})}}{\text{TP}_{c}}+\frac{n_{c(\text{obs})}}{\text{TP}_{c}}}\right)$

## Appendix C. Joint model applied to MM data

Figure 2: Individual weighted residuals (IWRES) from the bi-exponential submodel for each biomarker. Table 5: Posterior mean and 95% credible interval of the bi-exponential submodel parameters for each biomarker. Statistically significant variables are shown in bold, except for variance parameters. Interpretation | Parameter | M-spike | FLC
---|---|---|---
Baseline | $\exp\\{\theta_{1kl}\\}$ |  18.60 (17.62, 19.61) |  21.20 (19.40, 23.33)
Growth | $\exp\\{\theta_{2kl}\\}$ |  0.24 (0.22, 0.26) |  0.18 (0.16, 0.20)
Decay | $\exp\\{\theta_{3kl}\\}$ |  4.76 (4.42, 5.12) |  3.74 (3.33, 4.18)
Residual error variance | $\sigma^{2}_{kl}$ | 0.06 (0.05, 0.06) | 0.14 (0.13, 0.14)
| $\omega_{11kl}$ | 0.79 (0.71, 0.87) | 2.93 (2.70, 3.18)
Covariance | $\omega_{12kl}$ |  -0.22 (-0.30, -0.15) |  -1.37 (-1.55, -1.19)
matrix for | $\omega_{13kl}$ |  0.40 (0.31, 0.48) |  1.62 (1.41, 1.85)
random effects | $\omega_{22kl}$ | 0.70 (0.61, 0.79) | 1.60 (1.40, 1.81)
| $\omega_{23kl}$ | 0.05 (-0.03, 0.13) | -0.04 (-0.22, 0.14)
| $\omega_{33kl}$ | 0.87 (0.76, 1.00) | 2.35 (2.06, 2.66)
Table 6: Baseline covariates and categorical outcomes of the illustrative patients. Variable | Patient A | Patient B
---|---|---
Sex |  |
Ethnicity | Information not disclosed for the purpose of anonymization
Age |  |
ECOG | 1 | 1
ISS | Not reported | Stage III
Albumin (serum, g/L) | 33 | 40
B2B (serum, mg/L) | Not reported | Not reported
Creatinine (serum, mg/dL) | 0.7 | 1.21
Hemoglobin (g/dL) | 11.1 | 7
LDH (serum, U/L) | Not reported | 498
Lymphocyte (count, $\times 10^{9}$/L) | 0.6 | 0.9
Neutrophil (count, $\times 10^{9}$/L) | Not reported | 2.6
Platelet (count, $\times 10^{9}$/L) | 317 | 69
IgA (serum, g/L) | Not reported | Not reported
IgG (serum, g/L) | Not reported | Not reported
IgM (serum, g/L) | Not reported | Not reported
RVd duration (in years) | 2.092 | 0.326
Physician’s chosen treatment (outcome) | Pomalidomide | Carfilzomib

## Appendix D. Competing model applied to MM data

Figure 3: Variable importance (VI) ranking from the categorical model (no longitudinal information). Bar length, lower and upper limits of the horizontal lines represent average, minimum and maximum VI scores, respectively, based on 50 runs. Table 7: Relative risk (RR) and its 95% credible interval (CI) from the categorical model parameters (no longitudinal information) for Carfilzomib vs Other and Pomalidomide vs Other, and variable importance (VI) ranking. Statistically significant variables are shown in bold. Variable | Category | VI ranking | Carfilzomib vs Other | Pomalidomide vs Other
---|---|---|---|---
RR (95% CI) | RR (95% CI)
Sex | Female | 16 | 0.89 (0.69, 1.15) | 1.04 (0.79, 1.38)
Ethnicity | Non-Hisp. Black | 9 | 1.04 (0.73, 1.47) | 0.81 (0.53, 1.21)
Other | 1.22 (0.89, 1.67) | 1.10 (0.78, 1.54)
ECOG | 1 | 4 | 0.79 (0.57, 1.10) | 0.99 (0.68, 1.44)
2+ |  0.59 (0.37, 0.92) |  1.67 (1.10, 2.53)
ISS | Stage II | 3 | 1.26 (0.86, 1.84) | 0.78 (0.52, 1.17)
Stage III | 1.14 (0.69, 1.88) | 1.19 (0.70, 2.02)
Age | – | 1 |  0.77 (0.69, 0.86) |  1.28 (1.10, 1.49)
Albumin | – | 5 | 1.15 (0.98, 1.37) |  0.82 (0.72, 0.92)
B2M | – | 8 | 1.06 (0.85, 1.32) | 0.92 (0.71, 1.19)
Creatine | – | 7 | 1.16 (0.99, 1.35) | 0.99 (0.82, 1.18)
Hemoglobin | – | 12 | 0.93 (0.81, 1.08) | 1.02 (0.87, 1.19)
LDH | – | 6 |  1.33 (1.11, 1.58) | 1.06 (0.86, 1.29)
Lymphocyte | – | 13 | 0.95 (0.82, 1.09) | 0.89 (0.76, 1.04)
Neutrophil | – | 17 | 1.01 (0.86, 1.19) | 0.98 (0.82, 1.17)
Platelet | – | 15 | 0.97 (0.83, 1.13) | 0.99 (0.83, 1.17)
IgA | – | 11 | 0.98 (0.80, 1.21) | 0.92 (0.74, 1.14)
IgG | – | 10 | 0.95 (0.79, 1.15) | 0.95 (0.77, 1.18)
IgM | – | 14 | 0.94 (0.76, 1.14) | 0.89 (0.71, 1.10)
RVd duration | – | 2 |  0.76 (0.62, 0.84) | 1.08 (0.94, 1.23)