# Seemingly Unrelated Multi-State processes: a Bayesian semiparametric
approach

Andrea Cremaschi  Singapore Institute for Clinical Sciences, A*STAR, Singapore
Raffele Argiento  Yong Loo Lin School of Medicine, National University of
Singapore, Singapore  Singapore Institute for Clinical Sciences, A*STAR,
Singapore  Division of Science, Yale-NUS College, Singapore  Department of
Statistical Science, University College London, UK  Maria De Iorio  Yong Loo
Lin School of Medicine, National University of Singapore, Singapore  Singapore
Institute for Clinical Sciences, A*STAR, Singapore  Division of Science, Yale-
NUS College, Singapore  Department of Statistical Science, University College
London, UK  Cai Shirong  Singapore Institute for Clinical Sciences, A*STAR,
Singapore  Yap Seng Chong  Yong Loo Lin School of Medicine, National
University of Singapore, Singapore  Singapore Institute for Clinical Sciences,
A*STAR, Singapore  Michael J. Meaney  Singapore Institute for Clinical
Sciences, A*STAR, Singapore  Yong Loo Lin School of Medicine, National
University of Singapore, Singapore  Department of Psychiatry, McGill
University, Montreal, Canada  Michelle Z. L. Kee  Singapore Institute for
Clinical Sciences, A*STAR, Singapore

Abstract: Many applications in medical statistics as well as in other fields
can be described by transitions between multiple states (e.g. from health to
disease) experienced by individuals over time. In this context, multi-state
models are a popular statistical technique, in particular when the exact
transition times are not observed. The key quantities of interest are the
transition rates, capturing the instantaneous risk of moving from one state to
another. The main contribution of this work is to propose a joint
semiparametric model for several possibly related multi-state processes
(Seemingly Unrelated Multi-State, SUMS, processes), assuming a Markov
structure for the transitions over time. The dependence between different
processes is captured by specifying a joint random effect distribution on the
transition rates of each process. We assume a flexible random effect
distribution, which allows for clustering of the individuals, overdispersion
and outliers. Moreover, we employ a graph structure to describe the dependence
among processes, exploiting tools from the Gaussian Graphical model
literature. It is also possible to include covariate effects. We use our
approach to model disease progression in mental health. Posterior inference is
performed through a specially devised MCMC algorithm.

Keywords: Multi-State Models, Normalized Point Processes, Graphical Models,
Mixture Models, Markov Chain Monte Carlo

##  1 Introduction

Biomedical data are characterized by a high number of different variables, in
many cases mostly categorical and recorded on a (nowadays often large) set of
subjects. This is mainly due to the practice in clinical settings to record
the absence/presence of symptoms and/or to use ordinal scales to represent
disease markers. Typically, we only observe clinical variables at fixed time
points (usually corresponding to follow up or hospital visits), and as such
these variables are interval-censored (i.e., panel data). The objective of
clinical studies is often to model a patient’s disease progression, as
captured by the evolution over time of one or more responses of interest, e.g.
representing the disease status, and associated clinical markers. A popular
approach to disease progression is to use multi-state models describing the
transition of individuals among multiple states in continuous time (see, for
instance, Cook, 1999; Jackson et al., 2011; van den Hout et al., 2015; De
Iorio et al., 2018). In this framework, it is straightforward to include time-
homogeneous covariates and time varying continuous covariates (leading to a
Markov regression model).

In this work, we propose a joint modelling approach for several categorical
variables evolving simultaneously through time. More in details, our approach
is based on a combination of ideas from multi-state models, seemingly
unrelated regression (Zellner, 1963; Wang, 2010), Gaussian Graphical models
and Product Partition Models with Covariates (PPMx) (Müller et al., 2011). In
a Bayesian framework, we define a joint model for several multi-state
processes, which represent the evolution of, for instance, clinical markers of
interest as in the disease progression application of Section 3. The main idea
is that the different markers provide complementary information on the
underlying health status and, as such, they are regarded as stochastic
processes defined on a finite state-space, evolving in continuous time
according to dependent Markov processes. We link the different Markov
processes through the specification of a flexible prior distribution on the
instantaneous transition rates, specifically a mixture distribution with
random number of components (Argiento and De Iorio, 2019). In this way, we
build a robust modelling strategy, which leads to covariate-driven clustering
of the subjects and enables the inclusion of different types of
covariates/responses in a natural and efficient way (Barcella et al., 2017).
Each multi-state process is then, conditionally on the vector of instantaneous
transition rates, independent from the other processes, resembling the
seemingly unrelated regressions (SUR) setting of Zellner (1963). Furthermore,
we allow the dependence structure between the transition rates to be encoded
into a random graph, which is also object of posterior inference, as it is
done in sparse SUR models (SSUR, Wang, 2010). Thus, the nature of the
dependence is learnt from the data, spanning from independence to full inter-
dependence. As such, we refer to our model as Seemingly Unrelated Multi-State
(SUMS) processes. Briefly, the proposed model allows for: (i) multiple
responses; (ii) processes with more than two states; (iii) patient- and
process-specific times of observation; (iv) inclusion of mixed-type
covariates; (v) covariate-driven clustering of the subjects; (vi) missing
initial state information.

One of the main advantages of our modelling strategy is that the relationship
between different multi-state processes is encoded into a graph structure.
Indeed, if there is an edge linking two processes, it means that they are
conditionally dependent, while the lack of an edge implies conditional
independence. This gives insight into the co-regulatory mechanisms of the
different processes. This is relevant in many application as often it is of
interest also to identify important factors affecting disease progression, for
better prognosis and therapeutic choices. Moreover, the model allows for the
inclusion of time-homogeneous covariates (of any type) and time-varying
continuous covariates in a regression component, for which standard variable
selection techniques (e.g. shrinkage, spike and slab priors) can be employed.

The manuscript is organised as follows: Section 2 introduces the SUMS model,
by presenting how its key components - the joint multi-state model, the
mixture prior with unknown number of components and the graphical structure
describing the dependence among processes - interrelate, as well as the
specifically designed MCMC algorithm. Section 3 presents an application of the
model to the analysis of mental health indicators obtained from the GUSTO
cohort study. Section 4 concludes the work. In Supplemetary Material we
include a detailed description of the algorithm and of the GUSTO dataset, a
sensitivity analysis and a simulation study, as well as further results from
the analysis of the GUSTO data.

##  2 SUMS: Seemingly Unrelated Multi-State Processes

###  2.1 Modelling of multi-state processes

Multi-state models can be used to describe how an individual moves between a
set of states in continuous time. In this work, we focus on multi-state
processes for panel data, where the states of several processes are observed
only at certain time points, and their exact transition times are not known.
For each $h=1,\dots,p$, let $\\{Y^{(h)}(t),t\in\mathbb{R}^{+}\\}$ be a
continuous time Markov process, where $Y^{(h)}(t)$ represents the state of the
$h-$th process over time, with state-space
$\mathcal{S}^{(h)}=\\{1,\dots,d^{(h)}\\}$ of dimension $d^{(h)}$, i.e.
$Y^{(h)}(t)\in\mathcal{S}^{(h)}$; the elements of $\mathcal{S}^{(h)}$
represent the states that the $h$-th process can visit between transitions.
The exact times of transition of the processes $Y^{(h)}(t)$ are not known, but
in applications the processes are observed on a discrete set of time points,
$\bm{t}_{i}^{(h)}=\left(t_{i1}^{(h)},\ldots,t_{in^{(h)}_{i}}^{(h)}\right)$,
where $n^{(h)}_{i}$ denotes the number of observed time points for the $i$-th
individual and $h$-th process. Notice that the times of observation and their
number are both process- and subject-specific. We indicate with
$Y^{(h)}_{ij}=Y^{(h)}(t_{ij}^{(h)})$ the value of the $h$-th process
$Y^{(h)}(t)$ at the $j$-th observed time $t_{ij}^{(h)}$ for the $i$-th
subject. Hence, for each subject $i=1,\dots,N$, we observe the random vector
$\bm{Y}^{(h)}_{i}=\left(Y^{(h)}_{i1},\dots,Y^{(h)}_{in^{(h)}_{i}}\right)$,
whose joint distribution is the finite-dimensional law of the process
$Y^{(h)}(t)$ at the times of observation, for $h=1,\dots,p$. The aim of this
work is to jointly model the processes $Y^{(h)}(t)$, capturing their time
evolution and possible dependencies. For each process, we assume that the
Markov property holds, i.e. conditionally on current and past events, future
transitions only depend on the current state. The probability law of the
$h$-th process $Y^{(h)}(t)$ is assigned via the matrix of instantaneous
transition rates $\bm{Q}^{(h)}(t)=[\lambda^{(h)}(r,s;t)]_{r,s}$, which is also
time-dependent, and whose elements are the instantaneous transition rates
$\lambda^{(h)}(r,s;t)$ with $r,s\in\mathcal{S}^{(h)}$. In what follows, the
vector
$\bm{\lambda}^{(h)}(t)=\\{\lambda^{(h)}(r,s;t):r,s\in\mathcal{S}^{(h)},r\neq
s\\}$, of dimension $D_{p}=\sum_{h=1}^{p}d^{(h)}(d^{(h)}-1)$, indicates the
off-diagonal transition rates of the matrix $\bm{Q}^{(h)}(t)$ at time $t>0$,
concatenated by row from top to bottom. For simplicity, we indicate a
transition of the $h$-th process between different states of
$\mathcal{S}^{(h)}$ with the notation $r\rightarrow s$. The instantaneous
transition rates can be made covariate-dependent by specifying a Cox
proportional hazard model. This allows the inclusion of both time-homogeneous
covariates as well as time-varying continuous covariates. Alternatively, a
semi-proportional intensity model can be easily specified for the covariates
as in Kim et al. (2012). Note that the decision of including either type of
covariates is process-specific. The time-homogeneous covariates are
straightforwardly incorporated in the model, and we denote them here by
$\bm{X}^{(h)}_{i}=\left(X^{(h)}_{i1},\dots,X^{(h)}_{ig^{(h)}}\right)$, for the
$i$-th individual and $h$-th process. On the other hand, the time-varying
continuous ones, denoted by
$\bm{Z}^{(h)}_{i}(t)=\left(Z^{(h)}_{i1}(t),\dots,Z^{(h)}_{iq^{(h)}}(t)\right)$,
require additional assumptions. They are usually included by assuming a piece-
wise constant effect over each interval of observations (Andersen et al.,
2012), or by modelling them as longitudinal processes, linking their
distribution to the ones of the multi-state processes via the inclusion of
suitable random effects (Ferrer et al., 2016). The first option has a clear
computational advantage, while the latter has the potential to yield better
inference on the overall disease progression. In the application presented in
Section 3 we are not provided with any time-varying continuous covariates.
However, the code provided with this manuscript allows for the implementation
of the first method. This assumption leads to a piecewise constant model for
the instantaneous transition rates $\lambda^{(h)}_{ij}(r,s)$, with
$r,s\in\mathcal{S}^{(h)}$, and for the matrix
$\bm{Q}^{(h)}_{ij}:=\bm{Q}^{(h)}_{i}(t^{(h)}_{ij})$, for
$j=1,\dots,n^{(h)}_{i}$, $i=1,\dots,N$ and $h=1,\dots,p$. The model for the
instantaneous $\log$-transition rates for $i=1,\dots,N$ is then:

$\log\left(\lambda^{(h)}_{ij}(r,s)\right)=\phi^{(h)}_{i}(r,s)+\bm{X}^{(h)}_{i}\bm{\beta}^{(h)}_{rs}+\bm{Z}^{(h)}_{ij}\bm{\gamma}^{(h)}_{rs},\quad
j=1,\dots,n^{(h)}_{i},h=1,\dots,p$ (1)

where $\phi^{(h)}_{i}(r,s)$ represents the baseline transition rate (on a log
scale) of a transition $r\rightarrow s$. The parameters
$\bm{\beta}^{(h)}_{rs}\in\mathbb{R}^{g^{(h)}}$ and
$\bm{\gamma}^{(h)}_{rs}\in\mathbb{R}^{q^{(h)}}$ are the vectors of regression
coefficients for the $h$-th process and the $r\rightarrow s$ transition.

Let $\epsilon^{(h)}_{ij}=t^{(h)}_{ij}-t^{(h)}_{ij-1}$ indicate the length of
the $j$-th time interval, for $j=2,\dots,n^{(h)}_{i}$, $i=1,\dots,N$ and
$h=1,\dots,p$. Thanks to the piecewise constant assumption, the Chapman-
Kolmogorov equations can be solved to obtain the interval-specific transition
probabilities,
$\bm{p}_{ij}^{(h)}(\bm{\lambda}^{(h)}_{ij},\epsilon^{(h)}_{ij})=\left\\{p_{ij}^{(h)}(r,s;\bm{\lambda}^{(h)}_{ij},\epsilon^{(h)}_{ij}):r,s\in\mathcal{S}^{(h)}\right\\}$,
for the vector of random variables $\bm{Y}^{(h)}_{i}$ (see Ross et al., 1996).
When $d^{(h)}=2$, closed-form solutions are readily available (Cox and Miller,
1977), while problems involving more than three states are usually tackled
numerically (Moler and Van Loan, 2003). It can be shown that for each process
$h$ a unique stationary distribution exists (Grimmet and Sterzaker, 2001), and
we denote it by $\bm{\pi}^{(h)}_{ij}\left(k;\bm{\lambda}^{(h)}_{ij}\right)$,
with $k\in\mathcal{S}^{(h)}$, highlighting the fact that these are functions
of the subject-specific instantaneous transition rates (see Ross et al., 1996,
for details). The stationary distribution can be used as marginal distribution
for modelling the state of the processes at time $j=1$, considering the
vectors of instantaneous transition rates $\bm{\lambda}^{(h)}_{i1}$, in
contrast to the general practice in multi-state modelling of specifying the
model conditionally on the state at the first time of observation. This is
important, as it allows Bayesian imputation of missing observations at time
one, since they are treated as unknown parameters in the model. This aspect is
particularly useful in our application, where the initial time presents a non-
negligible missing rate.

The specification of the process- and subject-specific transition
probabilities, together with the existence of the stationary distribution,
leads to the joint likelihood for the vector of observed states
$\bm{Y}^{(h)}_{i}$ for $i=1,\dots,N$ and $h=1,\dots,p$, as follows:

$p\left(\bm{Y}\mid\bm{\lambda}^{Y}\right)=\prod_{i=1}^{N}\prod_{h=1}^{p}\prod_{j=2}^{n^{(h)}_{i}}\left(p^{(h)}_{ij}\left(Y^{(h)}_{ij-1},Y^{(h)}_{ij};\bm{\lambda}^{(h)}_{ij},\epsilon^{(h)}_{ij}\right)\right)\pi^{(h)}_{i1}\left(Y^{(h)}_{i1};\bm{\lambda}^{(h)}_{i1}\right)$
(2)

where $\bm{Y}$ and $\bm{\lambda}^{Y}$ indicate the multi-dimensional arrays
containing the observation vectors $\bm{Y}^{(h)}_{i}$ and the instantaneous
transition rate vectors $\bm{\lambda}^{(h)}_{i}$, while $p_{ij}^{(h)}$ denotes
the transition probabilities and $\pi^{(h)}_{i1}$ is the (stationary)
distribution at time one.

###  2.2 Random effect distribution and relationship with SUR

Consider the vector of $\log$-baseline transition rates for the $h$-th process
$\bm{\phi}^{(h)}_{i}=\\{\phi^{(h)}_{i}(r,s):r\rightarrow s\\}$, for
$h=1,\dots,p$ and for each subject $i$, and let
$\bm{\phi}_{i}=(\bm{\phi}^{(1)}_{i},\dots,\bm{\phi}^{(p)}_{i})$ be the vector
containing the $\log$-baseline transition rates of all the $p$ processes. To
capture the inter-individual heterogeneity and allow for clustering of the
subjects, we choose as random effect distribution for
$\bm{\phi}_{1},\dots,\bm{\phi}_{N}$ a mixture prior with random number of
components, where the distribution of the weights is given by the
normalization of a finite point process, as proposed by Argiento and De Iorio
(2019). This approach has several advantages, allowing for flexible modelling
of the weights in the mixture as well as efficient posterior computations
(e.g., as compared to traditional reversible jump algorithms for mixture
models). In more details, we assume the following mixture prior:

$\displaystyle\bm{\phi}_{i}$ $\displaystyle=\bm{\phi}^{\star}_{c_{i}},\quad
i=1,\dots,N$
$\displaystyle\bm{\phi}^{\star}_{1},\dots,\bm{\phi}^{\star}_{M}\mid M$
$\displaystyle\stackrel{{\scriptstyle\rm
iid}}{{\sim}}P_{0}(\bm{\phi}^{\star}\mid\bm{\theta})$
$\displaystyle\mathbb{P}(c_{i}=m)$ $\displaystyle\propto S_{m},\quad
i=1,\dots,N$ (3) $\displaystyle S_{1},\dots,S_{M}$
$\displaystyle\stackrel{{\scriptstyle\rm
iid}}{{\sim}}\text{Gamma}(\gamma_{S},1)$ $\displaystyle M-1$
$\displaystyle\sim\text{Poi}(\Lambda)$

where we denote by $\text{Gamma}(a,b)$ the Gamma distribution with mean $a/b$,
and by $\text{Poi}(\Lambda)$ the Poisson distribution with mean $\Lambda$. The
variables $\bm{c}=(c_{1},\dots,c_{N})$ indicate the component allocations of
the subjects and their corresponding prior probabilities are proportional to
the unnormalized weights $\bm{S}=\left(S_{1},\dots,S_{M}\right)$. As shown by
Argiento and De Iorio (2019), posterior computation is greatly simplified via
the introduction of a latent variable, conditionally on which the unnormalized
weights of the mixture become independent. This computational trick is
borrowed from the Bayesian nonparametric literature (James et al., 2009).
Finally, the vectors $\bm{\phi}^{\star}_{1},\dots,\bm{\phi}^{\star}_{M}$ are a
finite sequence of locations for the mixture distribution and are,
conditionally on the number of components $M$, i.i.d. from the base measure
$P_{0}$. The specification of a joint random effect distribution for
$\bm{\phi}_{1},\dots,\bm{\phi}_{N}$ in model (2.2) and the choice of $P_{0}$
are crucial in our modelling strategy, as it will be shown in Section 2.3,
since this allows inference on the shared dependence structure among the
components of the vectors $\bm{\phi}^{\star}_{m}$, for $m=1,\dots,M$ and,
consequently, on the dependence structure among the $p$ different processes.
As an alternative, a Bayesian nonparametric prior could have been specified as
random effect distribution such as the Dirichlet process (De Iorio et al.,
2018) and the beta-Dirichlet process prior (Kim et al., 2012), or, taking a
complete different approach, flexible modelling of the baseline transition
intensities can be achieved using penalised splines (Kneib and Hennerfeind,
2008).

Our modelling approach resembles the one underlying the SUR framework of
Zellner (1963), where $p$ different regression models are linked by specifying
a joint error distribution, usually multivariate normal. The SUR methodology
is one of the main techniques for handling multiple responses and offers a way
to share information between models which are seemingly unrelated, since they
describe different data-generating processes. However, since these are
observed for the same set of subjects and measurements are taken on often
related processes, the study of their interdependency is of great interest in
most applications. For this reason SUR-type models have gained vast popularity
in different fields, such as Phenomics (Houle et al., 2010; Banterle et al.,
2018). In our application, for instance, where we deal with several processes
associated to different aspects of maternal mental health (e.g., depression,
anxiety, sleep quality), it is important to understand the relationships
between such processes in order to have a broader view of the phenomenon under
study. As in the SUR framework, in our context each process is modelled by its
own seemingly unrelated multi-state Markov process, but then they are related
through the joint random effect distribution on
$\bm{\phi}=\left(\bm{\phi}_{1},\dots,\bm{\phi}_{N}\right)$. Motivated by this
parallelism, we name the proposed model as Seemingly Unrelated Multi-State
(SUMS) processes.

###  2.3 Gaussian Graphical model

We use tools from the Gaussian Graphical models literature to describe the
dependence among the $p$ processes. Referring to model (2.2), we assume that
$\bm{\phi}^{\star}_{1},\dots,\bm{\phi}^{\star}_{M}\mid
M\stackrel{{\scriptstyle\rm
iid}}{{\sim}}P_{0}=\text{N}\left(\bm{\mu},\bm{\Omega}_{G}\right)$, where the
key modelling feature is the specification of the prior on the precision
matrix $\bm{\Omega}_{G}$ conditional on a graph $G$, which captures the
conditional dependence structure among the $\log$-baseline transition rates.
The novelty of our modelling strategy is that $G$ is modelled conditionally on
another random graph $G_{0}$, which characterizes the dependence structure
among processes and is one of the main aspects of our inference.

In details, consider the graph $G_{0}=(V_{0},E_{0})$, defined over the set of
nodes $V_{0}=\\{1,\dots,p\\}$, i.e. each node in the graph corresponds to a
multi-state process $Y^{(h)}(t)$. The edge set $E_{0}$ is formed of the pairs
$E_{0}\subseteq\\{(h,k)\in V_{0}\times V_{0}:h<k\\}$ such that an edge exists
between nodes $h,k\in V_{0}$. We consider only simple graphs, i.e. undirected
graphs, without self-loops nor multiple edges. As mentioned earlier, to
introduce dependence among the elements of the vector
$\bm{\phi}^{\star}_{m}=(\bm{\phi}^{\star,(1)}_{m},\dots,\bm{\phi}^{\star,(p)}_{m})$,
with $\bm{\phi}^{\star,(h)}_{m}=\\{\phi^{\star,(h)}_{m}(r,s):r\rightarrow
s\\}$ for $h=1,\dots,p$, we define a second graph $G$ whose structure is
determined by $G_{0}$. In particular, we let $G=(V,E)$ be the graph whose
nodes are the indices of the vector $\bm{\phi}^{\star}_{m}$, i.e.
$V=\\{1,\dots,D_{p}\\}$, with $D_{p}=\sum_{h\in V_{0}}d^{(h)}(d^{(h)}-1)$. $G$
is a deterministic function of $G_{0}$ specified as follows. First, whatever
the form of $G_{0}$, there exists an edge in $G$ between transition rates of
the same process. Therefore, an empty graph $G_{0}$ corresponds to a graph $G$
with $p$ cliques, one for each process. Second, if there is an edge between
nodes $h$ and $k$ in $G_{0}$ (i.e., $(h,k)\in E_{0}$), then there is and edge
between all the possible pairs of elements of $\bm{\phi}^{\star,(h)}_{m}$ with
those of $\bm{\phi}^{\star,(k)}_{m}$. An illustration for the case of three
binary processes is given in Figure 1. We write $G=f(G_{0})$, $f$ being the
transformation described above. Note that $f$ is bijective and, as such, the
specification of a prior on $G_{0}$ implies a prior on $G$. This construction
is advantageous in terms of dimension reduction, as the dimension of the graph
space where $G_{0}$ is defined can be significantly smaller than the one of
$G$, leading to more efficient exploration of the posterior space.

$G_{0}$ $1$ $2$ $3$ $G$$(1,2)$$(2,1)$$(1,2)$$(2,1)$$(1,2)$$(2,1)$$f$ Figure 1:
Example of graphical structures underlying the SUMS processes: the graph
$G_{0}$ describes the conditional dependence between the processes in
$V_{0}=\\{1,2,3\\}$, while the graph $G$ is obtained as a deterministic
function from $G_{0}$ by connecting the edges of the corresponding transitions
between states 1 and 2 of each process.

Following the literature on GGMs, the conditional independence structure of
the multivariate Gaussian vectors
$\bm{\phi}^{\star}_{m}\sim\text{N}(\bm{\mu},\bm{\Omega}_{G})$, for
$m=1,\dots,M$, is described by constraining the elements of the precision
matrix $\bm{\Omega}_{G}$ (Dempster, 1972). Namely, two elements of the vector
$\bm{\phi}^{\star}_{m}$ are, conditionally on the others, independent if and
only if there is a zero in the corresponding entry of the precision matrix
$\bm{\Omega}_{G}$. Since $G$ is a deterministic function of $G_{0}$, it is the
latter that encodes the conditional independence structure of the vectors
$\bm{\phi}^{\star}_{m}$, for $m=1,\dots,M$ (see Figure 1). The standard
conjugate prior for the precision matrix $\bm{\Omega}_{G}$ is the G-Wishart
distribution, specified conditionally on the graph structure $G$ (Roverato,
2002). The last component needed to fully specify this part of the model, is
the prior distribution for the graph $G$. We do not assign this prior
directly, but rather it is inherited by the prior we choose for the graph
$G_{0}$:

$\pi(G_{0}\mid\eta)\propto\eta^{|E_{0}|}(1-\eta)^{\binom{p}{2}-|E_{0}|},\quad\eta\in(0,1)$

where $|E_{0}|$ is the number of edges in graph $G_{0}$ (i.e., the size of
$E_{0}$), while $\binom{p}{2}$ is the number of possible graphs with nodes
$V_{0}=\\{1,\ldots,p\\}$. This prior is equivalent to assuming a Bernoulli
prior with probability of success (here inclusion) $\eta$ on each edge of the
graph $G_{0}$, independently across edges. Small values of $\eta$ favour
sparser graphs (Armstrong et al., 2009). Finally, we point out that, while the
prior for the graph $G_{0}$ is defined over all possible graphs, including the
non-decomposable ones, the resulting prior distribution on $G$ is defined on a
restricted space due to the clique constraints imposed on the transitions of
the same process which need to be fully connected.

###  2.4 Relationship with PPMx models

The SUMS model has wide applicability in biomedical research as some processes
can be regarded as responses and some as covariates. Indeed, time-varying
categorical covariates are very common in the field of medical research, for
instance in association with the monitoring of a patient’s disease status over
time. In the application to disease progression in Section 3, some processes
represent mental health outcomes of interest, while others correspond to
categorical clinical markers, and the goal of the analysis is to model the
joint evolution of outcomes and clinical factors to determine how the symptom
variables influence the disease course. Handling of time-varying multivariate
categorical information can be problematic in several applications and the
SUMS approach provide a natural framework to deal with this problem. In a
Bayesian framework, De Iorio et al. (2018) discuss possible solutions and
propose an approach based on a latent health function borrowing ideas from
Item Response Theory (Thissen and Steinberg, 2009). The latter approach,
although computationally efficient, does not allow for a direct quantification
of the covariate effect on the clinical response of interest and it may lead
to identifiability problems. A simpler and more common approach to deal with
time-varying categorical covariates is to introduce appropriate dummy
variables, considerably increasing the number of parameters to be estimated,
resulting in slower computations and lower effectiveness in high dimensional
problems. Another computational effective solution is to summarize the
covariates into an often arbitrary time-varying score, but at the cost of
losing information and interpretability.

When some of the multi-state processes can be seen as covariates, then the
SUMS model has interesting connections with Product Partition Models with
Covariates (PPMx), very popular in the Bayesian nonparametric literature
Müller et al. (1996, 2011). Indeed, our approach provides a flexible and
robust modelling strategy, which leads to covariate-driven clustering of the
subjects and enables the inclusion of different types of covariates/responses
in a natural and efficient way (Barcella et al., 2017). We now clarify the
relationship between SUMS and PPMx. The main modelling idea behind the PPMx is
to include covariate information into the partition model (e.g., into a
Dirichlet Process Mixture model framework) by treating each covariate as a
random variable (Müller et al., 1996, 2011; Barcella et al., 2017). In the
SUMS approach we can consider a set of processes
$(Y^{(1)}(t),\dots,Y^{(p_{Y})}(t))$ as responses and another set, denoted by
$(H^{(1)}(t),\dots,H^{(p_{H})}(t))$, as explanatory factors. Then Eq.(2)
specifies a suitable probability model for the joint vector of processes
$\\{(Y^{(h)}(t),H^{(l)}(t));h=1,\ldots,p_{Y};l=1,\ldots,p_{H};t\in\mathbb{R}^{+}\\}$,
with $p_{Y}+p_{H}=p$. Let $\bm{c}$ be the set of allocation variables
introduced in (2.2), and let $\rho_{N}$ be the partition of the indices
$\\{1,\dots,N\\}$ induced by $\bm{c}$. We indicate by $C_{j}$ the set of
indices belonging to the $j$-th cluster, i.e.
$C_{j}=\\{i\in\\{1,\dots,N\\}\mid c_{i}=j\\}$, and thus a partition with
$K_{N}$ clusters corresponds to $\rho_{N}=\\{C_{1},\dots,C_{K_{N}}\\}$. In the
PPMx framework, the SUMS model induces a prior on the partition $\rho_{N}$
which depends on the covariates $\bm{H}$:

$p\left(\rho_{N}\mid\bm{H}\right)=V(N,K_{N})\prod_{j=1}^{K_{N}}\mathcal{C}(C_{j})\mathcal{G}(\bm{H}^{\star}_{j})$
(4)

where $\mathcal{C}(C_{j})$ is the cohesion, i.e. a function of the $j$-th
cluster $C_{j}$, $\mathcal{G}(\bm{H}^{\star}_{j})$ is the similarity, i.e. a
function of the array of covariates corresponding to the subjects in cluster
$j$, denoted as $\bm{H}^{\star}_{j}:=\\{H_{i}:i\in C_{j}\\}$, for
$j=1,\dots,K_{N}$, and $V(N,K_{N})$ is a constant depending only on the sample
size $N$ and the number of clusters $K_{N}$. The cohesion function
$\mathcal{C}$ expresses prior information about the partition, such as the
average size of a cluster, while the similarity function $\mathcal{G}$
captures the contribution of the covariates to the clustering structure. The
presence of $\mathcal{G}$ in (4) allows subjects with similar covariates to be
more likely assigned to the same cluster. Under our modelling assumptions it
can be shown that (4) is given by:

$\displaystyle p\left(\rho_{N}\mid\bm{H}\right)\propto
p\left(\rho_{N}\right)p\left(\bm{H}\mid\rho_{N}\right)=$ $\displaystyle
V(N,K_{N})\prod_{j=1}^{K_{N}}\frac{\Gamma(\gamma_{S}+n_{j})}{\Gamma(\gamma_{S})}\int\left(\prod_{i\in
C_{j}}p\left(\bm{H}_{i}\mid\bm{\lambda}^{H}_{i}\right)\right)P_{0}(d\bm{\phi}^{\star,(p_{Y}+1)}_{j},\dots,d\bm{\phi}^{\star,(p)}_{j})$
(5)

See Supplementary Material Section 1 for a proof. The proposed model induces
the similarity function $\mathcal{G}(\bm{H}^{\star}_{j})=\int\prod_{i\in
C_{j}}p\left(\bm{H}_{i}\mid\bm{\lambda}^{H}_{i}\right)P_{0}(d\bm{\phi}^{\star,(p_{Y}+1)}_{j},\dots,d\bm{\phi}^{\star,(p)}_{j})$,
for $j=1,\dots,K_{N}$. The similarity function $\mathcal{G}$ is not known in
closed form, differently from the common PPMx specification, where the
similarity function is usually obtained from a conjugate model for the
covariates vector via marginalization to simplify computations. In the
proposed approach, the evaluation of (5) would require an expensive numerical
approximation. For this reason, we resort to a conditional MCMC algorithm
analogous to the one proposed by Argiento and De Iorio (2019), not requiring
the evaluation of the integral in (5).

###  2.5 MCMC Algorithm

Posterior inference is performed through a MCMC algorithm, described in
details in Supplementary Material Section 2. The numerous non-conjugate
updates required by the proposed model are tackled using adaptive Metropolis-
Hastings sampling schemes (Haario et al., 2001; Atchadé et al., 2005), which
need an additional short burn-in period. Additionally, inference under the
proposed model is challenging given the presence of the graphs $G$ and
$G_{0}$. We adopt the birth-and-death approach of Mohammadi et al. (2015), and
extend their algorithm to accommodate for MCMC moves on cliques instead of
single edges, recalling that each edge in $G_{0}$ corresponds to a clique in
$G$ through the map $f$. Indeed, the original algorithm of Mohammadi et al.
(2015) is based on theoretical results from the GGM literature (see Wang et
al., 2012), which can be extended to our modelling setting. In Supplementary
Material Section 4, we also compare the performance of our model with the
approach of De Iorio et al. (2018) and with two alternative versions of the
proposed model (i.e., DP and parametric versions). The results of the
comparison show that the proposed model outperforms the parametric approach,
as well as the nonparametric competitors in terms of clustering, leading to
comparable results with respect to the estimation of regression coefficients.

##  3 Application to the GUSTO study

The GUSTO study (Growing Up in Singapore Towards healthy Outcomes, Soh et al.,
2014) is a longitudinal birth cohort study started in 2009 and involving
Singaporean mothers and their children. The study is one of the most carefully
phenotyped parent-offspring cohorts, focusing on the roles of foetal,
developmental and epigenetic factors involved in early body composition as
well as neuro-development. In this work we consider data on $N=301$ mothers,
followed during pre- and post-natal periods, starting from three months before
childbirth. The main focus of the analysis is understanding the relationship
among five psychometric indicators obtained from specific questionnaires: the
Beck’s Depression Inventory II (BDI II, Beck et al., 1961); the Edinburgh
Postnatal Depression Scale (EPDS, Matthey et al., 2006); the State-Trait
Anxiety Inventory (STAI, Spielberger et al., 1983) that can be decomposed into
two different scores describing the anxious states (STAI-s), reflecting
characteristics that can vary with time, and the anxiety traits (STAI-t),
reflecting more stable characteristics; and the Pittsburgh Sleep Quality Index
(PSQI, Buysse et al., 1989). The score ranges of these questionnaires are
discretized to obtain clinically relevant categories, and are recorded at
different time points, as reported in Supplementary Material Table 2. These
five processes represent time-varying categorical observations and are
modelled jointly via SUMS, to capture significant relationships between them
(see also Meaney et al., 2018). In our setting, the four mental health
indicators (BDI, EPDS, STAI-s, STAI-t) represent the main clinical responses
of interest, while the sleep quality indicator (PSQI) is treated as a time-
varying categorical covariate. For all processes, we assume missingness at
random and impute missing values at the first time of observation from their
full conditionals (see Section 2.2 of Supplementary Material). We are also
provided with information regarding socio-demographic and clinical
characteristics, as well as scoring obtained from additional questionnaires
measuring personality traits. In particular, we have individual scores for the
Big Five Inventory (BFI, John et al., 1999) (including the scores for
Extraversion, Agreeableness, Conscientiousness, Neuroticism, Openness, and
Liking) and for the Maternal Childhood Adversity (MCA, Bouvette-Turcot et al.,
2015). Many of the remaining covariates are time-homogeneous categorical,
while no time-varying continuous covariates are available. The time-
homogeneous continuous covariates are centred and scaled so that each column
has null mean and unitary standard deviation, thus estimating the
corresponding regression coefficients $\bm{\beta}^{(h)}$ on the same scale
across processes. The full set of covariates ($g^{(h)}=22$, for $h=1,\dots,4$,
including dummy coding for the categorical ones) is described in more details
in Supplementary Material Table 3, and is included in the specification of the
four psychometric processes, but not of PSQI.

#### Full model specification

We describe the full model used in the application presented in this section,
which is the same implemented in the sensitivity analysis on the
hyperparameters $\Lambda$ and $\gamma_{S}$ appearing in Supplementary Section
3. For each $i=1,\ldots,N$:

$\displaystyle\bm{Y}^{(1)}_{i},\dots,\bm{Y}^{(4)}_{i},\bm{H}_{i}\mid\\{\bm{\lambda}^{(h)}_{i},\bm{\beta}^{(h)},h=1,\ldots,4\\},\bm{\lambda}^{H}_{i}\stackrel{{\scriptstyle\rm
ind}}{{\sim}}\prod_{h=1}^{4}p(\bm{Y}^{(h)}_{i}\mid\bm{\lambda}^{(h)}_{i},\bm{\beta}^{(h)})p(\bm{H}_{i}\mid\bm{\lambda}^{H}_{i})$
$\displaystyle\log\left(\lambda^{(h)}_{i}(r,s)\right)=\phi^{\star,h}_{c_{i}}(r,s)+\bm{X}^{(h)}_{i}\bm{\beta}^{(h)}_{rs},\
r\rightarrow s,\quad h=1,\dots,4$
$\displaystyle\log\left(\lambda^{H}_{i}(r,s)\right)=\phi^{\star,H}_{c_{i}}(r,s),\
r\rightarrow s$ $\displaystyle\bm{\beta}^{(h)}\sim\text{MN}_{g^{(h)}\times
d^{(h)}(d^{(h)}-1)}(\bm{0},\bm{U}_{\bm{\beta}^{(h)}},\bm{V}_{\bm{\beta}^{(h)}}),\quad
h=1,\dots,4$
$\displaystyle\bm{\phi}^{\star}_{m}=(\bm{\phi}^{\star,1}_{m},\dots,\bm{\phi}^{\star,4}_{m},\bm{\phi}^{\star,H}_{m})\mid
M,\bm{\mu},\bm{\Omega}_{G}\sim
P_{0}=\text{N}_{D_{p}}(\bm{\mu},\bm{\Omega}_{G}),\quad m=1,\dots,M$
$\displaystyle\bm{\mu},\bm{\Omega}_{G}\mid
G,\bm{m}_{\mu},k_{0}\sim\text{N}_{D_{p}}(\bm{\mu}\mid\bm{m}_{\bm{\mu}},k_{0}\Omega_{G})G\text{-Wishart}_{G}(\bm{\Omega}_{G}\mid\nu,\bm{\Psi})$
(6) $\displaystyle k_{0}\sim\text{Gamma}(a_{k_{0}},b_{k_{0}})$
$\displaystyle\mathbb{P}(c_{i}=m)\propto S_{m},\quad m=1,\dots,M$
$\displaystyle S_{1},\dots,S_{M}\mid M,\gamma_{S}\stackrel{{\scriptstyle\rm
iid}}{{\sim}}\text{Gamma}(\gamma_{S},1)$ $\displaystyle
M-1\mid\Lambda\sim\text{Poi}(\Lambda)$ $\displaystyle G=f(G_{0}),\quad
p(G_{0})\propto\eta^{|E_{0}|}(1-\eta)^{\binom{p}{2}-|E_{0}|}$

where we indicate with $\bm{\phi}_{m}^{\star}$ the vectors of unique
$\log$-baseline transition rates for the $m$-th component in the model and $M$
is the unknown number of components in the mixture. Here
$\bm{c}=\\{c_{i},i=\ldots,N\\}$ represents the allocation vector, i.e. it
specifies to which component the $i$-th observation is assigned to,
characterised by $\bm{\phi}_{i}=\bm{\phi}_{c_{i}}^{\star}$. The probability of
$c_{i}$ being equal to the $m$-th component of the mixture is proportional to
the unnormalized weights $S_{m}$, for $m=1,\dots,M$. Therefore, due to the
discrete property of the mixing measure, the parameters $\bm{\phi}_{i}$ are
assigned to $K_{N}$ different clusters, with $K_{N}\leq M$. We impose a
conditionally conjugate hyper-prior on $k_{0}$, and fix the hyperparameters
$\gamma_{S},\Lambda$. We refer to Argiento and De Iorio (2019) for a thorough
discussion on prior specification in mixture models with unknown number of
components. However, we point out that the mixture component of the model is
specified conditionally to the graph structure $G$. Finally,
$\text{MN}_{n\times p}(\bm{0},\bm{U},\bm{V})$ is the matrix-variate Normal
distribution of dimension $n\times p$ centred on the null matrix $\bm{0}$ and
with covariance matrices $\bm{U}$ and $\bm{V}$ of dimensions $n\times n$ and
$p\times p$, respectively.

#### Hyper-Prior elicitation

We need to specify the hyperparameters for the priors in the three components
of the model: the transition rates, the mixture model with random number of
components and the graphical model.

In order to induce sparsity in the graph structure and identify meaningful
relationship between the SUMS processes, we set the a-priori probability of
edge inclusion to $\eta=0.1$. The hyperparameters of the centring measure
$P_{0}$ are such that $\bm{m}_{\bm{\mu}}=\bm{0}$,
$k_{0}\sim\text{Gamma}(1,1)$, $\nu=D_{p}+2$ and
$\bm{\Psi}=\mathbb{I}_{D_{p}}/\nu$, where $\mathbb{I}_{p}$ is the identity
matrix of size $p$. In the case of a full graph $G$, the latter corresponds to
$\mathbb{E}(\bm{\Omega}_{G}\mid G)=\mathbb{I}_{D_{p}}$. The regression
coefficients $\bm{\beta}^{(h)}$ are a-priori independent and identically
distributed, i.e.
$\bm{U}_{\bm{\beta}^{(h)}}=\bm{V}_{\bm{\beta}^{(h)}}=\mathbb{I}_{g^{(h)}d^{(h)}(d^{(h)}-1)}$,
for $h=1,\ldots,4$. The mixture prior for the $\log$-baseline transition rates
$\bm{\phi}^{\star}_{1},\dots,\bm{\phi}^{\star}_{M}$ is controlled by the
hyperparameters $\Lambda$ and $\gamma_{S}$. These parameters determine the
distribution of the number of components and the corresponding allocation of
the subjects, and are the object of an extensive sensitivity analysis
presented in Supplementary Section 3. In this application, we fix these
parameters to $\Lambda=0.01$ and $\gamma_{S}=0.1$.

#### Posterior inference

We run the MCMC algorithm described in Section 2.5 for 50000 iterations, after
an initial burn-in period of 1000 iterations used to initialise the adaptive
Metropolis-Hastings, discarding 40000 iterations as burn-in and thinning every
2, obtaining a final sample of 5000.

We explore the relationship between the multi-state processes by imposing
dependencies via the graphical model approach described in Section 2.3.
Inference on the posterior distribution of the graphical structure $G_{0}$ is
obtained by reporting the posterior edge inclusion probability for each pair
of nodes. In Figure 2 we report the posterior median graph, obtained by
including only those edges with posterior edge inclusion probability greater
than 0.5 (Barbieri et al., 2004). The four clinical mental health indicators
BDI, EPDS, STAI-s and STAI-t are strongly associated, presenting a clique in
the posterior median graph. Interestingly, the sleep quality index PSQI is
only related to the anxiety indices STAI-s and STAI-t, forming a clique as
well. Links between probable anxiety and sleeping quality have been reported
in previous studies (Swanson et al., 2011; Ibrahim and Foldvary-Schaefer,
2012), and it is confirmed by our findings. Moreover, as previously reported,
poor sleep quality may feed into poor emotional and mental health states
(Ruiz-Robledillo et al., 2015; Osnes et al., 2019).

Figure 2: Posterior median graph of $G_{0}$: each edge included in the graph
has posterior edge inclusion probability greater than 0.5. Each edge of the
median graph is labelled with the corresponding posterior edge inclusion
probability.

Another important aspect of the proposed model is the possibility of including
covariates in the specification of the transition rates via (1). Posterior
inference on the coefficient $\bm{\beta}^{(h)}$, for $h=1,\dots,4$ is not
trivial, due to the high number of parameters involved. The importance of each
covariate can be assessed through Bayes Factors (BF), defined as the ratio of
the marginal contributions derived from the model with the corresponding
regression coefficient set to zero versus the full model (Kass and Raftery,
1995). Closed form expressions for the Bayes Factor under the SUMS model are
not available, and thus we use the Savage-Dickey density ratio method
(Wagenmakers et al., 2010; Verdinelli and Wasserman, 1995). The applicability
of this method is guaranteed by the component-wise assumption of independence
a-priori for the regression coefficients $\bm{\beta}^{(h)}$, for $h=1,\dots,4$
(see the full model specification in (3)). For each process $h$, the values of
$-\log_{10}\left(\text{BF}^{(h)}_{jk}\right)$ are reported in the heatmap of
Figure 3, for $j=1,\dots,g^{(h)}$ and $k=1,\dots,d^{(h)}(d^{(h)}-1)$. The
magnitude of $-\log_{10}\left(\text{BF}^{(h)}_{jk}\right)$ measures the
evidence in favour of the full model (Kass and Raftery, 1995). The majority of
the coefficients is characterized by a low value of
$-\log_{10}\left(\text{BF}^{(h)}_{jk}\right)$, supporting the hypothesis of no
association, particularly in the case of the STAI processes. However, some
coefficients are characterized by $-\log_{10}(BF)$ values above 1 or 2,
indicating strong evidence in support of the inclusion of the corresponding
covariate in the specific process. Of particular interest are the coefficients
relative to the BFI and MCA scores, representing different traits of
personality, trauma and parental relationship. We present the posterior mean
and 95% credible intervals of the regression coefficients relative to BFI and
MCA in details in Figures 5 of Supplementary Material Section 6.

Figure 3: Heatmap of Bayes Factors ($-\log_{10}(\text{BF}^{(h)}_{jk})$) for
the individual regression coefficients $\beta^{(h)}_{jk}$, for
$j=1,\dots,g^{(h)}$, $k=1,\dots,d^{(h)}(d^{(h)}-1)$ and $h=1,\dots,4$. Each
row refers to a different covariate included in the model. Each column is
associated with a possible transition for each process, excluding PSQI which
is modelled as an explanatory factor.

The personality traits of the mothers as described by the BFI scores have been
previously associated with increased likelihood for both antenatal and
postnatal mood disorder traits (Ritter et al., 2000; Leigh and Milgrom, 2008).
Our analysis supports this as BFI traits have a relevant impact on both BDI
and EPDS (95% credible interval does not contain zero). An interesting result
appears through the estimates of the BFI’s Neuroticism dimension, which
characterizes transitions 2 $\rightarrow$ 3 (deterioration, positive
regression coefficient) and 3 $\rightarrow$ 1 (improvement, negative
regression coefficient) in both BDI and EPDS scores, indicating that higher
Neuroticism scores are associated with higher depressive symptoms during the
peripartum period (Kitamura et al., 1993; O’hara and Swain, 1996). On the
other hand, Openness and Conscientiousness in EPDS (see Supplementary Figure
5) positively influence the transition 3 $\rightarrow$ 1 (improvement). We
also notice the effects of BFI’s Extraversion and Agreeableness differ for BDI
and EPDS’s transitions. This could be explained by the fact that the social
behaviors associated with Extraversion and Agreeableness are distinct (Tobin
et al., 2000; Jensen-Campbell and Graziano, 2001). Extraverts tend to actively
seek out social interactions, whereas people scoring high on Agreeableness
prefer harmonious relationships. Maternal history of developmental adversity
is linked to increased risk for depression (Leigh and Milgrom, 2008), of which
childhood abuse is a strong risk factor (Seng et al., 2014), as highlighted by
the importance of the MCA covariate for the transition 2 $\rightarrow$ 3
(deterioration) in BDI (see Figure 3 and Supplementary Figure 5). This result
is also confirmed by Mandelli et al. (2015) who found that women who were
victims of childhood neglect or abuse were at least twice as likely to suffer
from depression. The quality of relationship with the women’s parents may also
contribute to maternal developmental adversity. Mothers who received low
parental care and high control during childhood are at risk for peripartum
anxiety (Grant et al., 2012) and depression (McMahon et al., 2005).

The choice of the mixture prior (2.2) as random effect distribution for the
vector of $\log$-transition rates $\bm{\phi}_{1},\dots,\bm{\phi}_{N}$ allows
for clustering of the subjects. Inference on the random partition is shown in
Supplementary Material Figure 6, where the posterior distributions of the
number of clusters, components and of the co-clustering probabilities are
reported. An estimate of the random partition induced on the subjects under
study is obtained by minimizing the Binder’s loss function (Binder, 1978) with
equal costs. We obtain a partition with three clusters, which also corresponds
to the posterior mode of the number of clusters. The three clusters contain
140, 126 and 35 subjects, respectively, and are labelled according to their
sizes in decreasing order. In Figure 4 we report the posterior distribution of
$\phi^{(h)}(r,s)$ conditional on the Binder’s partition, for $r\rightarrow s$
and $h=1,\dots,p$. Cluster-specific estimates of transition rates differ among
clusters (see Supplementary Material Section 7 for a discussion). For
instance, transition rates corresponding to improvement in the BDI or EPDS
scores are higher in Clusters 1 and 2 rather than Cluster 3. The binary
processes (STAI-s, STAI-t and PSQI) also seem to present differences between
clusters in the same direction, identifying Cluster 3 as the one most prone to
a deterioration of the mental health status of its subjects.

Figure 4: Posterior means and 95% credible intervals of the instantaneous
$\log$-transition rates $\phi^{(h)}(r,s)$ for each process $h=1,\dots,p_{0}$.
The vertical dashed lines represent the value 0, while the horizontal
continuous lines divide the estimates for the five processes. The estimates
are obtained by fixing the partition of the subjects to the Binder’s
partition, and re-running the algorithm for the conditional model. Each sub-
plot refers to one of the clusters in the fixed partition.

##  4 Conclusions

Observations on time-evolving related processes are very common in biomedical
applications and beyond. In this work, we present a Bayesian semiparametric
approach for joint modelling of several multi-state Markov processes
describing an individual’s transitions between different states in continuous
time. The proposed model builds on the multi-state Markov models, GGM and PPMx
literature. The different multi-state processes are linked by imposing a
flexible prior distribution for the instantaneous transition rates, which
allows for data-driven clustering of the subjects. The dependence among the
processes is captured by a graph and posterior inference is performed through
a tailored MCMC algorithm.

The proposed model finds wide applicability, due to its flexibility,
interpretability and relative ease of computations. In this work, we analyse
data from the GUSTO cohort study with the aim of understanding the evolution
and relationships between mental health indicators over time. Our findings are
in agreement with existing medical literature and shed more light on the
influence of childhood and parental factors on mental health progression.
Potential extensions include higher order Markov dependency and joint
modelling of multi-state processes and continuous longitudinal trajectories.

A possible alternative to our approach is to represent the categorical
covariates with continuous Gaussian latent variables linked to the categorical
outcome by thresholding (Albert and Chib, 1993), allowing for the inclusion of
a time component through auto-regressive terms in the likelihood (e.g.
Barcella et al., 2018). To the best of our knowledge, this strategy has not
been employed in the context of multi-state models, and it represents an
interesting direction for future developments. However, this formulation could
suffer from limited interpretability (García-Zattera et al., 2007) and could
induce further computational challenges (Zhang et al., 2006).

## References

  * Albert and Chib [1993] James H Albert and Siddhartha Chib.  Bayesian analysis of binary and polychotomous response data.  _Journal of the American statistical Association_ , 88(422):669–679, 1993.
  * Andersen et al. [2012] Per K Andersen, Ornulf Borgan, Richard D Gill, and Niels Keiding.  _Statistical models based on counting processes_.  Springer Science & Business Media, 2012.
  * Argiento and De Iorio [2019] R. Argiento and M. De Iorio.  Is infinity that far? a bayesian nonparametric perspective of finite mixture models.  _arXiv preprint arXiv:1904.09733_ , 2019.
  * Armstrong et al. [2009] H. Armstrong, C. K Carter, K. F. K. Wong, and R. Kohn.  Bayesian covariance matrix estimation using a mixture of decomposable graphical models.  _Statistics and Computing_ , 19(3):303–316, 2009\.
  * Atchadé et al. [2005] Yves F Atchadé, Jeffrey S Rosenthal, et al.  On adaptive markov chain monte carlo algorithms.  _Bernoulli_ , 11(5):815–828, 2005.
  * Banterle et al. [2018] Marco Banterle, Leonardo Bottolo, Sylvia Richardson, Mika Ala-Korpela, MR Järvelin, and Alex Lewin.  Sparse variable and covariance selection for high-dimensional seemingly unrelated bayesian regression.  _bioRxiv_ , page 467019, 2018.
  * Barbieri et al. [2004] Maria Maddalena Barbieri, James O Berger, et al.  Optimal predictive model selection.  _The annals of statistics_ , 32(3):870–897, 2004\.
  * Barcella et al. [2017] William Barcella, Maria De Iorio, and Gianluca Baio.  A comparative review of variable selection techniques for covariate dependent dirichlet process mixture models.  _Canadian Journal of Statistics_ , 45(3):254–273, 2017.
  * Barcella et al. [2018] William Barcella, Maria De Iorio, and James Malone-Lee.  Modelling correlated binary variables: an application to lower urinary tract symptoms.  _Journal of the Royal Statistical Society: Series C (Applied Statistics)_ , 67(4):1083–1100, 2018.
  * Beck et al. [1961] A. T. Beck, C. H. Ward, M. Mendelson, J. Mock, and J. Erbaugh.  An Inventory for Measuring Depression.  _Archives of General Psychiatry_ , 4(6):561–571, 06 1961.  ISSN 0003-990X.  doi: 10.1001/archpsyc.1961.01710120031004.  URL https://doi.org/10.1001/archpsyc.1961.01710120031004.
  * Binder [1978] David A Binder.  Bayesian cluster analysis.  _Biometrika_ , 65(1):31–38, 1978.
  * Bouvette-Turcot et al. [2015] Andrée-Anne Bouvette-Turcot, AS Fleming, A Wazana, MB Sokolowski, H Gaudreau, A Gonzalez, J Deslauriers, JL Kennedy, M Steiner, MJ Meaney, et al.  Maternal childhood adversity and child temperament: An association moderated by child 5-httlpr genotype.  _Genes, Brain and Behavior_ , 14(3):229–237, 2015\.
  * Buysse et al. [1989] Daniel J Buysse, Charles F Reynolds, Timothy H Monk, Susan R Berman, David J Kupfer, et al.  The pittsburgh sleep quality index: a new instrument for psychiatric practice and research.  _Psychiatry res_ , 28(2):193–213, 1989.
  * Cook [1999] Richard J Cook.  A mixed model for two-state markov processes under panel observation.  _Biometrics_ , 55(3):915–920, 1999.
  * Cox and Miller [1977] David Roxbee Cox and Hilton David Miller.  _The theory of stochastic processes_ , volume 134.  CRC press, 1977.
  * De Iorio et al. [2018] Maria De Iorio, Natacha Gallot, Beatriz Valcarcel, and Lucy Wedderburn.  A bayesian semiparametric markov regression model for juvenile dermatomyositis.  _Statistics in medicine_ , 37(10):1711–1731, 2018\.
  * Dempster [1972] A. Dempster.  Covariance selection.  _Biometrics_ , 28:157–175, 1972.
  * Ferrer et al. [2016] Loïc Ferrer, Virginie Rondeau, James Dignam, Tom Pickles, Hélène Jacqmin-Gadda, and Cécile Proust-Lima.  Joint modelling of longitudinal and multi-state processes: application to clinical progressions in prostate cancer.  _Statistics in medicine_ , 35(22):3933–3948, 2016\.
  * García-Zattera et al. [2007] María José García-Zattera, Alejandro Jara, Emmanuel Lesaffre, and Dominique Declerck.  Conditional independence of multivariate binary data with an application in caries research.  _Computational statistics & data analysis_, 51(6):3223–3234, 2007.
  * Grant et al. [2012] Kerry-Ann Grant, Alison Bautovich, Catherine McMahon, Nicole Reilly, Leo Leader, and Marie-Paule Austin.  Parental care and control during childhood: associations with maternal perinatal mood disturbance and parenting stress.  _Archives of women’s mental health_ , 15(4):297–305, 2012.
  * Grimmet and Sterzaker [2001] Geoffrey Grimmet and David Sterzaker.  _Probability and Random Processes._ Oxford University press, 2001.
  * Haario et al. [2001] Heikki Haario, Eero Saksman, Johanna Tamminen, et al.  An adaptive metropolis algorithm.  _Bernoulli_ , 7(2):223–242, 2001.
  * Houle et al. [2010] David Houle, Diddahally R Govindaraju, and Stig Omholt.  Phenomics: the next challenge.  _Nature reviews genetics_ , 11(12):855–866, 2010\.
  * Ibrahim and Foldvary-Schaefer [2012] Sally Ibrahim and Nancy Foldvary-Schaefer.  Sleep disorders in pregnancy: implications, evaluation, and treatment.  _Neurologic clinics_ , 30(3):925–936, 2012.
  * Jackson et al. [2011] Christopher H Jackson et al.  Multi-state models for panel data: the msm package for r.  _Journal of statistical software_ , 38(8):1–29, 2011.
  * James et al. [2009] L. F. James, A. Lijoi, and I. Prünster.  Posterior analysis for normalized random measures with independent increments.  _Scandinavian Journal of Statistics_ , 36(1):76–97, 2009.
  * Jensen-Campbell and Graziano [2001] Lauri A Jensen-Campbell and William G Graziano.  Agreeableness as a moderator of interpersonal conflict.  _Journal of personality_ , 69(2):323–362, 2001\.
  * John et al. [1999] Oliver P John, Sanjay Srivastava, et al.  The big five trait taxonomy: History, measurement, and theoretical perspectives.  _Handbook of personality: Theory and research_ , 2(1999):102–138, 1999.
  * Kass and Raftery [1995] Robert E Kass and Adrian E Raftery.  Bayes factors.  _Journal of the american statistical association_ , 90(430):773–795, 1995.
  * Kim et al. [2012] Yongdai Kim, Lancelot James, and Rafael Weissbach.  Bayesian analysis of multistate event history data: beta-dirichlet process prior.  _Biometrika_ , 99(1):127–140, 2012.
  * Kitamura et al. [1993] T Kitamura, S Shima, M Sugawara, and MA Toda.  Psychological and social correlates of the onset of affective disorders among pregnant women.  _Psychological medicine_ , 23(4):967–975, 1993\.
  * Kneib and Hennerfeind [2008] Thomas Kneib and Andrea Hennerfeind.  Bayesian semi parametric multi-state models.  _Statistical Modelling_ , 8(2):169–198, 2008\.
  * Leigh and Milgrom [2008] Bronwyn Leigh and Jeannette Milgrom.  Risk factors for antenatal depression, postnatal depression and parenting stress.  _BMC psychiatry_ , 8(1):24, 2008.
  * Mandelli et al. [2015] L Mandelli, C Petrelli, and A Serretti.  The role of specific early trauma in adult depression: a meta-analysis of published literature. childhood trauma and adult depression.  _European psychiatry_ , 30(6):665–680, 2015.
  * Matthey et al. [2006] S Matthey, C Henshaw, S Elliott, and B Barnett.  Variability in use of cut-off scores and formats on the edinburgh postnatal depression scale–implications for clinical and research practice.  _Archives of women’s mental health_ , 9(6):309–315, 2006.
  * McMahon et al. [2005] Catherine McMahon, Bryanne Barnett, Nicholas Kowalenko, and Christopher Tennant.  Psychological factors associated with persistent postnatal depression: past and current relationships, defence styles and the mediating role of insecure attachment style.  _Journal of affective disorders_ , 84(1):15–24, 2005.
  * Meaney et al. [2018] Michael Meaney, Linde van Lee, Shirong Cai, See Ling Loy, Elaine KH Tham, Fabian KP Yap, Keith M Godfrey, Peter D Gluckman, Lynette PC Shek, Oon Hoe Teoh, Daniel YT Goh, et al.  Relation of plasma tryptophan concentrations during pregnancy to maternal sleep and mental well-being: The gusto cohort.  _Journal of affective disorders_ , 225:523–529, 2018.
  * Mohammadi et al. [2015] Abdolreza Mohammadi, Ernst C Wit, et al.  Bayesian structure learning in sparse gaussian graphical models.  _Bayesian Analysis_ , 10(1):109–138, 2015.
  * Moler and Van Loan [2003] Cleve Moler and Charles Van Loan.  Nineteen dubious ways to compute the exponential of a matrix, twenty-five years later.  _SIAM review_ , 45(1):3–49, 2003.
  * Müller et al. [1996] Peter Müller, Alaattin Erkanli, and Mike West.  Bayesian curve fitting using multivariate normal mixtures.  _Biometrika_ , 83(1):67–79, 1996.
  * Müller et al. [2011] Peter Müller, Fernando Quintana, and Gary L Rosner.  A product partition model with regression on covariates.  _Journal of Computational and Graphical Statistics_ , 20(1):260–278, 2011.
  * O’hara and Swain [1996] Michael W O’hara and Annette M Swain.  Rates and risk of postpartum depression—a meta-analysis.  _International review of psychiatry_ , 8(1):37–54, 1996.
  * Osnes et al. [2019] Rannveig S Osnes, John Olav Roaldset, Turid Follestad, and Malin Eberhard-Gran.  Insomnia late in pregnancy is associated with perinatal anxiety: a longitudinal cohort study.  _Journal of affective disorders_ , 248:155–165, 2019.
  * Ritter et al. [2000] Christian Ritter, Stevan E Hobfoll, Justin Lavin, Rebecca P Cameron, and Michael R Hulsizer.  Stress, psychosocial resources, and depressive symptomatology during pregnancy in low-income, inner-city women.  _Health Psychology_ , 19(6):576, 2000.
  * Ross et al. [1996] Sheldon M Ross, John J Kelly, Roger J Sullivan, William James Perry, Donald Mercer, Ruth M Davis, Thomas Dell Washburn, Earl V Sager, Joseph B Boyce, and Vincent L Bristow.  _Stochastic processes_ , volume 2.  Wiley New York, 1996.
  * Roverato [2002] A. Roverato.  Hyper inverse Wishart distribution for non-decomposable graphs and its application to Bayesian inference for Gaussian graphical models.  _Scandinavian Journal of Statistics_ , 29(3):391–411, 2002.
  * Ruiz-Robledillo et al. [2015] N Ruiz-Robledillo, Catarina Canário, CC Dias, L Moya-Albiol, and Bárbara Figueiredo.  Sleep during the third trimester of pregnancy: the role of depression and anxiety.  _Psychology, health & medicine_, 20(8):927–932, 2015.
  * Seng et al. [2014] Julia S Seng, Wendy D’Andrea, and Julian D Ford.  Complex mental health sequelae of psychological trauma among women in prenatal care.  _Psychological Trauma: Theory, Research, Practice, and Policy_ , 6(1):41, 2014.
  * Soh et al. [2014] Shu-E Soh, Mya Thway Tint, Peter D Gluckman, Keith M Godfrey, Anne Rifkin-Graboi, Yiong Huak Chan, Walter Stünkel, Joanna D Holbrook, Kenneth Kwek, Yap-Seng Chong, et al.  Cohort profile: Growing up in singapore towards healthy outcomes (gusto) birth cohort study.  _International journal of epidemiology_ , 43(5):1401–1409, 2014.
  * Spielberger et al. [1983] Charles D Spielberger, Richard L Gorsuch, R Lushene, Peter R Vagg, and Gerard A Jacobs.  _Manual for the state-trait anxiety inventory_.  Consulting Psychologists Press, Palo Alto, CA, 1983.
  * Swanson et al. [2011] Leslie M Swanson, Scott M Pickett, Heather Flynn, and Roseanne Armitage.  Relationships among depression, anxiety, and insomnia symptoms in perinatal women seeking mental health treatment.  _Journal of Women’s Health_ , 20(4):553–558, 2011\.
  * Thissen and Steinberg [2009] David Thissen and Lynne Steinberg.  Item response theory.  _The Sage handbook of quantitative methods in psychology_ , pages 148–177, 2009.
  * Tobin et al. [2000] Renée M Tobin, William G Graziano, Eric J Vanman, and Louis G Tassinary.  Personality, emotional experience, and efforts to control emotions.  _Journal of personality and social psychology_ , 79(4):656, 2000.
  * van den Hout et al. [2015] Ardo van den Hout, Jean-Paul Fox, and Rinke H Klein Entink.  Bayesian inference for an illness-death model for stroke with cognition as a latent time-dependent risk factor.  _Statistical methods in medical research_ , 24(6):769–787, 2015.
  * Verdinelli and Wasserman [1995] Isabella Verdinelli and Larry Wasserman.  Computing bayes factors using a generalization of the savage-dickey density ratio.  _Journal of the American Statistical Association_ , 90(430):614–618, 1995.
  * Wagenmakers et al. [2010] Eric-Jan Wagenmakers, Tom Lodewyckx, Himanshu Kuriyal, and Raoul Grasman.  Bayesian hypothesis testing for psychologists: A tutorial on the savage–dickey method.  _Cognitive psychology_ , 60(3):158–189, 2010\.
  * Wang [2010] Hao Wang.  Sparse seemingly unrelated regression modelling: Applications in finance and econometrics.  _Computational Statistics & Data Analysis_, 54(11):2866–2877, 2010.
  * Wang et al. [2012] Hao Wang, Sophia Zhengzi Li, et al.  Efficient gaussian graphical model determination under g-wishart prior distributions.  _Electronic Journal of Statistics_ , 6:168–198, 2012.
  * Zellner [1963] Arnold Zellner.  Estimators for seemingly unrelated regression equations: Some exact finite sample results.  _Journal of the American Statistical Association_ , 58(304):977–992, 1963.
  * Zhang et al. [2006] Xiao Zhang, W John Boscardin, and Thomas R Belin.  Sampling correlation matrices in bayesian models with correlated latent variables.  _Journal of Computational and Graphical Statistics_ , 15(4):880–896, 2006.