# Combining chains of Bayesian models with Markov melding

Andrew A. Manderson111MRC Biostatistics Unit, University of Cambridge, United
Kingdom, and The Alan Turing Institute, The British Library, London.
<EMAIL_ADDRESS>Robert J. B. Goudie222MRC Biostatistics
Unit, University of Cambridge, United Kingdom<EMAIL_ADDRESS>
(November, 2021)

###### Abstract

A challenge for practitioners of Bayesian inference is specifying a model that
incorporates multiple relevant, heterogeneous data. It may be easier to
instead specify distinct submodels for each source of data, then join the
submodels together. We consider chains of submodels, where submodels directly
relate to their neighbours via common quantities which may be parameters or
deterministic functions thereof. We propose _chained Markov melding_ , an
extension of Markov melding, a generic method to combine chains of submodels
into a joint model. One challenge we address is appropriately capturing the
prior dependence between common quantities within a submodel, whilst also
reconciling differences in priors for the same common quantity between two
adjacent submodels. Estimating the posterior of the resulting overall joint
model is also challenging, so we describe a sampler that uses the chain
structure to incorporate information contained in the submodels in multiple
stages, possibly in parallel. We demonstrate our methodology using two
examples. The first example considers an ecological integrated population
model, where multiple data are required to accurately estimate population
immigration and reproduction rates. We also consider a joint longitudinal and
time-to-event model with uncertain, submodel-derived event times. Chained
Markov melding is a conceptually appealing approach to integrating submodels
in these settings.

##  1 Introduction

The Bayesian philosophy is appealing in part because the posterior
distribution quantifies all sources of uncertainty. However, a joint model for
all data and parameters is a prerequisite to posterior inference, and in
situations where multiple, heterogeneous sources of data are available,
specifying such a joint model is a formidable task. Models that consider such
data are necessary to describe complex phenomena at a useful precision. One
possible approach begins by specifying individual submodels for each source of
data. These submodels could guide the statistician when directly specifying
the joint model, but to use the submodels only informally seems wasteful.
Instead, it may be preferable to construct a joint model by formally joining
the individual submodels together.

Some specific forms of combining data are well established. Meta-analyses and
evidence synthesis methods are widely used to summarise data, often using
hierarchical models (Ades and Sutton, 2006; Presanis _et al._ , 2014). Outside
of the statistical literature, a common name for combining multiple data is
_data fusion_ (Kedem _et al._ , 2017; Lahat _et al._ , 2015), though there are
many distinct methods that fall under this general name. Interest in
integrating data is not just methodological; applied researchers often collect
multiple disparate data sets, or data of different modalities, and wish to
combine them. For example, to estimate SARS-CoV-2 positivity Donnat _et al._
(2020) build an intricate hierarchical model that integrates both testing data
and self-reported questionnaire data, and Parsons _et al._ (2021) specify a
hierarchical model of similar complexity to estimate the number of injecting
drug users in Ukraine. Both applications specify Bayesian models with data-
specific components, which are united in a hierarchical manner. In
conservation ecology, _integrated population models_ (IPMs) (Zipkin and
Saunders, 2018) are used to estimate population dynamics, e.g. reproduction
and immigration rates, using multiple data on the same population. Such data
have standard models associated with them, such as the Cormack-Jolly-Seber
model (Lebreton _et al._ , 1992) for capture-recapture data, and the IPM
serves as the framework in which the standard models are combined. More
generally, the applications we list illustrate the importance of generic,
flexible methods for combining data to applied researchers.

_Markov melding_ (Goudie _et al._ , 2019) is a general statistical methodology
for combining submodels. Specifically, it considers $M$ submodels that share
some common quantity $\phi$, with each of the $m=1,\ldots,M$ submodels
possessing distinct parameters $\psi_{m}$, data $Y_{m}$, and form
$\text{p}_{m}(\phi,\psi_{m},Y_{m})$. Goudie _et al._ (2019) then propose to
combine the submodels into a joint model, denoted
$\text{p}_{\text{meld}}(\phi,\psi_{1},\ldots,\psi_{M},Y_{1},\ldots,Y_{M})$.
However, it is unclear how to integrate models where there is no single
quantity $\phi$ common to all submodels, such as for submodels that are linked
in a chain structure.

We propose an extension to Markov melding, which we call _chained Markov
melding_ 333 _“Chained graphs”_ were considered by Lauritzen and Richardson
(2002), however they are unrelated to our proposed model. We use “chained” to
emphasise the nature of the relationships between submodels., which
facilitates the combination of $M$ submodels that are in a chain structure.
For example, when $M=3$ we address the case in which submodel 1 and 2 share a
common quantity $\phi_{1\mathrel{{\cap}}2}$, and submodel 2 and 3 share a
different quantity $\phi_{2\mathrel{{\cap}}3}$. Our extension addresses
previously unconsidered complications including the distinct domains (and
possibly supports) of the common quantities, and the desire to capture
possible prior correlation between them. Two examples serve to illustrate our
methodology, which we introduce in the following section. The computational
effort required to fit a complex, multi-response model is a burden to the
model development process. We propose a multi-stage posterior estimation
method that exploits the properties of our chained melded model to reduce this
burden. We can parallelise aspects of the computation across the submodels,
using less computationally expensive techniques for some submodels. Reusing
existing software implementations of submodels, and subposterior samples where
available, is also possible. Multi-stage samplers can aid in understanding the
contribution of each submodel to the final posterior, and are used in many
applied settings, including hierarchical modelling (Lunn _et al._ , 2013) and
joint models (Mauff _et al._ , 2020).

One contribution of our work is to clarify the informal process commonly used
in applied analyses of summarising and/or approximating submodels for use in
subsequent analyses. The two most common approximation strategies seem to be
(i) approximating the subposterior of the common quantity with a normal
distribution for use in subsequent models (see, e.g. Jackson and White, 2018;
Nicholson _et al._ , 2021) and (ii) taking only a point estimate of the
subposterior, and treating it as a known value in further models. These
strategies may, but not always, produce acceptable approximations to the
chained melded model. Both the chained melded model and these approximation
strategies are examples of ‘multi-phase’ and ‘multi-source’ inference (Meng,
2014), with the melding approach most comprehensively accounting for
uncertainty.

###  1.1 Example introduction

In this section we provide a high-level overview of two applications that
require integrating a chain of submodels, with more details in Sections 4 and
5. Our first example decomposes a joint model into its constituent submodels
and rejoins them. This simple situation allows us to compare the output from
the chained melding process to the complete joint model, and is meant to
illustrate both the ‘chain-of-submodels’ notion and the mechanics of chained
melding. The second example is a realistic and complex setting in which the
combining of submodels without chained Markov melding is non-obvious. Our
comparator is the common technique of summarising previously considered
submodels with point estimates, and demonstrates the importance of fully
accounting for uncertainty.

####  1.1.1 An integrated population model for little owls

Integrated population models (IPMs) (Zipkin and Saunders, 2018) combine
multiple data to estimate key quantities governing the dynamics of a specific
population. Schaub _et al._ (2006) and Abadi _et al._ (2010) used an IPM to
estimate fecundity, immigration, and yearly survival rates for a population of
little owls. These authors collect and model three types of data, illustrated
in Figure 1. Capture-recapture data $Y_{1}$, and associated capture-recapture
submodel $\text{p}_{1}(\phi_{1\mathrel{{\cap}}2},\psi_{1},Y_{1})$, are
acquired by capturing and tagging owls each year, and then counting the number
of tagged individuals recaptured in subsequent years. Population counts
$Y_{2}$ are obtained by observing the number of occupied nesting sites, and
are modelled in
$\text{p}_{2}(\phi_{1\mathrel{{\cap}}2},\phi_{2\mathrel{{\cap}}3},\psi_{2},Y_{2})$.
Finally, nest-record data $Y_{3}$ counts both the number of reproductive
successes and possible breading pairs, and is associated with a submodel for
fecundity $\text{p}_{3}(\phi_{2\mathrel{{\cap}}3},\psi_{3},Y_{3})$. The
population count model $\text{p}_{2}$ shares the parameter
$\phi_{1\mathrel{{\cap}}2}$ with the capture-recapture model $\text{p}_{1}$,
and the parameter $\phi_{2\mathrel{{\cap}}3}$ with the fecundity model
$\text{p}_{3}$; each of the $m=1,2,3$ submodels has distinct, submodel-
specific parameters $\psi_{m}$. No single source of data is sufficient to
estimate all quantities of interest, so it is necessary to integrate the three
submodels into a single joint model to produce acceptably precise estimates of
fecundity and immigration rates. We will show that the chained Markov melding
framework developed in Section 2 encapsulates the process of integrating these
submodels, producing results that are concordant with the original joint IPM.

$Y_{2}$$\psi_{2}$$\phi_{2\mathrel{{\cap}}3}$$\phi_{1\mathrel{{\cap}}2}$$\phi_{2\mathrel{{\cap}}3}$$\phi_{1\mathrel{{\cap}}2}$$\psi_{3}$$Y_{3}$$\psi_{1}$$Y_{1}$$\text{p}_{1}$$\text{p}_{2}$$\text{p}_{3}$
Figure 1: A simplified DAG of the integrated population model (IPM) for the
little owls. The capture-recapture submodel ($\text{p}_{1}$) is surrounded by
the blue line, the population count submodel ($\text{p}_{2}$) by the black
line, and the fecundity submodel ($\text{p}_{3}$) by the red line. The
capture-recapture and population count submodels share parameters affecting
the juvenile and adult survival rate ($\phi_{1\mathrel{{\cap}}2}$), whilst the
parameter for fecundity is common to both the population count and fecundity
submodels ($\phi_{2\mathrel{{\cap}}3}$). The combination of all the submodels
forms the IPM.

####  1.1.2 Survival analysis with time varying covariates and uncertain event
times

Our second example considers the time to onset of respiratory failure (RF)
amongst patients in intensive care units, and factors that influence the onset
of RF. A patient can be said to be experiencing RF if the ratio of the partial
pressure of arterial blood oxygen (PaO2) to the faction of inspired oxygen
(FiO2) is less than 300mmHg (The ARDS Definition Task Force, 2012), though
this is not the only definition of RF. Patients’ PaO2/FiO2 (P/F) ratios are
typically measured only a few times a day. This relatively infrequency of P/F
ratio data, when combined with the intrinsic variability in each individual’s
blood oxygen level, results in significant uncertainty in about the time of
onset of RF.

Factors that influence the time to onset of RF are both longitudinal and time
invariant. Both types of data can be considered in _joint models_ (Rizopoulos,
2012), which are composed of two distinct submodels, one for each data type.
However, existing joint models are not able to incorporate the uncertainty
surrounding the event time, which may result in overconfident and/or biased
estimates of the parameters in the joint model.

Chained Markov melding offers a conceptually straightforward, Bayesian
approach to incorporating uncertain event times into joint models.
Specifically, we consider the event time as a submodel-derived quantity from a
hierarchical regression model akin to Lu and Meeker (1993). We call this
submodel the _uncertain event time_ submodel and denote it
$\text{p}_{1}(\phi_{1\mathrel{{\cap}}2},\psi_{1},Y_{1})$, where
$\phi_{1\mathrel{{\cap}}2}$ incorporates the event time. The survival submodel
$\text{p}_{2}(\phi_{1\mathrel{{\cap}}2},\phi_{2\mathrel{{\cap}}3},\psi_{2},Y_{2})$
uses the event time within $\phi_{1\mathrel{{\cap}}2}$, the common quantity,
as the response. We treat the longitudinal submodel,
$\text{p}_{3}(\phi_{2\mathrel{{\cap}}3},\psi_{3},Y_{3})$, separately from the
survival submodel, as is common in two-stage joint modelling (Mauff _et al._ ,
2020), and denote the subject-specific parameters that also appear in the
survival model as $\phi_{2\mathrel{{\cap}}3}$. Each of the $m=1,2,3$ has
submodel-specific data $Y_{m}$ and parameters $\psi_{m}$. The high level
submodel relationships are displayed as a DAG in Figure 2.

$Y_{1}$$\psi_{1}$$\phi_{1\mathrel{{\cap}}2}$$\phi_{1\mathrel{{\cap}}2}$$\psi_{2}$$Y_{2}$$\phi_{2\mathrel{{\cap}}3}$$\phi_{2\mathrel{{\cap}}3}$$\psi_{3}$$Y_{3}$$\text{p}_{1}$$\text{p}_{2}$$\text{p}_{3}$
Figure 2: A simplified DAG of the submodels consider in the survival analysis
example. The event time submodel $\text{p}_{1}$ defines the event time
$\phi_{1\mathrel{{\cap}}2}$ as noninvertible function of the other model
parameters (denoted by the dotted line), whilst the survival submodel
$\text{p}_{2}$ considers $\phi_{1\mathrel{{\cap}}2}$ as the response. The
longitudinal submodel $\text{p}_{3}$ has parameters
$\phi_{2\mathrel{{\cap}}3}$ in common with the survival submodel.

###  1.2 Markov melding

We now review Markov melding (Goudie _et al._ , 2019) before detailing our
proposed extension. As noted in the introduction, Markov melding is a method
for combining $M$ submodels
$\text{p}_{1}(\phi,\psi_{1},Y_{1}),\ldots,\text{p}_{M}(\phi,\psi_{M},Y_{M})$
which share the same $\phi$. When the submodel prior marginals
$\text{p}_{m}(\phi)$ are identical, i.e. $\text{p}_{m}(\phi)=\text{p}(\phi)$
for all $m$, it is possible to combine the submodels using _Markov
combination_ (Dawid and Lauritzen, 1993; Massa and Lauritzen, 2010)

$\displaystyle\text{p}_{\text{comb}}\left(\phi,\psi_{1},\ldots,\psi_{M},Y_{1},\ldots,Y_{M}\right)$
$\displaystyle=\text{p}(\phi)\prod_{m=1}^{M}\text{p}_{m}\left(\psi_{m},Y_{m}\mid\phi\right)$
(1)
$\displaystyle=\frac{\prod_{m=1}^{M}\text{p}_{m}\left(\phi,\psi_{m},Y_{m}\right)}{\text{p}(\phi)^{M-1}}.$

Markov combination is not immediately applicable when submodel prior marginals
are distinct, so Goudie _et al._ define a _marginal replacement_ procedure,
where individual submodel prior marginals are replaced with a common marginal
$\text{p}_{\text{pool}}(\phi)=h(\text{p}_{1}(\phi),\ldots,\text{p}_{M}(\phi))$
which is the result of a pooling function $h$ that appropriately summarises
all prior marginals (the choice of which is described below). The result of
marginal replacement is

$\text{p}_{\text{repl},m}(\phi,\psi_{m},Y_{m})=\text{p}_{\text{pool}}(\phi)\frac{\text{p}_{m}(\phi,\psi_{m},Y_{m})}{\text{p}_{m}(\phi)}.$
(2)

Goudie _et al._ show that $\text{p}_{\text{repl},m}(\phi,\psi_{m},Y_{m})$
minimises the Kullback–Leibler (KL) divergence between a distribution
$\text{q}(\phi,\psi_{m},Y_{m})$ and $\text{p}_{m}(\phi,\psi_{m},Y_{m})$ under
the constraint that $\text{q}(\phi)=\text{p}_{\text{pool}}(\phi)$, and that
marginal replacement is valid when $\phi$ is a deterministic function of the
other parameters in submodel $m$. Markov melding joins the submodels via the
Markov combination of the marginally replaced submodels

$\displaystyle\text{p}_{\text{meld}}(\phi,\psi_{1},\ldots,\psi_{M},Y_{1},\ldots,Y_{M})$
$\displaystyle=\text{p}_{\text{pool}}(\phi)\prod_{m=1}^{M}\text{p}_{\text{repl},m}(\psi_{m},Y_{m}\mid\phi),$
(3)
$\displaystyle=\text{p}_{\text{pool}}(\phi)\prod_{m=1}^{M}\frac{\text{p}_{m}(\phi,\psi_{m},Y_{m})}{\text{p}_{m}(\phi)}.$

####  1.2.1 Pooled prior

Goudie _et al._ proposed forming $\text{p}_{\text{pool}}(\phi)$ using linear
or logarithmic prior pooling (O’Hagan _et al._ , 2006)

$\displaystyle\text{p}_{\text{pool, lin}}(\phi)$
$\displaystyle=\frac{1}{K_{\text{lin}}(\lambda)}\sum_{m=1}^{M}\lambda_{m}\text{p}_{m}(\phi),\quad
K_{\text{lin}}(\lambda)=\int\sum_{m=1}^{M}\lambda_{m}\text{p}_{m}(\phi)\text{d}\phi$
(4) $\displaystyle\text{p}_{\text{pool, log}}(\phi)$
$\displaystyle=\frac{1}{K_{\text{log}}(\lambda)}\prod_{m=1}^{M}\text{p}_{m}(\phi)^{\lambda_{m}},\quad
K_{\text{log}}(\lambda)=\int\prod_{m=1}^{M}\text{p}_{m}(\phi)^{\lambda_{m}}\text{d}\phi$
(5)

where $\lambda=(\lambda_{1},\ldots,\lambda_{M})$ are nonnegative weights,
which are chosen subjectively to ensure $\text{p}_{\text{pool}}(\phi)$
appropriately represents prior knowledge about the common quantity. Two
special cases of pooling are of particular interest. _Product of experts (PoE)
pooling_ (Hinton, 2002) is a special case of logarithmic pooling that occurs
when $\lambda_{m}=1$ for all $m$. _Dictatorial pooling_ is a special case of
either pooling method when $\lambda_{m^{\prime}}=1$ and, for all $m\neq
m^{\prime}$, $\lambda_{m}=0$.

##  2 Chained model specification

Consider $m=1,\ldots,M$ submodels each with data $Y_{m}$ and parameters
$\theta_{m}$ denoted $\text{p}_{m}(\theta_{m},Y_{m})$, with $M\geq 3$. We
assume that the submodels are connected in a manner akin to a chain and so can
be ordered such that only ‘adjacent’ submodels in the chain have parameters in
common. Specifically we assume that submodels $m$ and $m+1$ have some
parameter $\phi_{m\mathrel{{\cap}}m+1}$ in common for $m=1,\ldots,M-1$. For
notational convenience define
$\phi_{1}=\phi_{1\mathrel{{\cap}}2},\phi_{M}=\phi_{M-1\mathrel{{\cap}}M}$ and
$\phi_{m}=(\phi_{m-1\mathrel{{\cap}}m},\,\phi_{m\mathrel{{\cap}}m+1})$ for
$m=2,\ldots,M-1$, so that $\phi_{m}\subseteq\theta_{m}$ denotes the parameters
in model $m$ shared with another submodel. The submodel-specific parameters of
submodel $m$ are thus $\psi_{m}=\theta_{m}\setminus\phi_{m}$. Define the
vector of all common quantities
$\boldsymbol{\phi}=(\phi_{1},\ldots,\phi_{M})=(\phi_{1\mathrel{{\cap}}2},\ldots,\phi_{M-1\mathrel{{\cap}}M})$,
and denote by $\boldsymbol{\phi}_{-m}$ the subvector of $\boldsymbol{\phi}$
excluding the $m$th element. It will also be convenient to define
$\boldsymbol{\psi}=(\psi_{1},\ldots,\psi_{M})$ and likewise
$\boldsymbol{Y}=(Y_{1},\ldots,Y_{M})$. Note that all components of
$\boldsymbol{\phi},\boldsymbol{\psi}$ and $\boldsymbol{Y}$ may themselves be
multivariate. Additionally, because $\phi_{m\mathrel{{\cap}}m+1}$ may be a
deterministic function of either $\theta_{m}$ or $\theta_{m+1}$ we refer to
$\phi_{m\mathrel{{\cap}}m+1}$ as a common parameter or a common quantity as
appropriate.

All submodels, and marginal and conditional distributions thereof, have
density functions that are assumed to exist and integrate to one. When
considering conditional distributions we assume that the parameter being
conditioned on has support in the relevant region. We define the $m$th
_subposterior_ as $\text{p}_{m}(\phi_{m},\psi_{m}\mid Y_{m})$.

###  2.1 Extending marginal replacement

We now define the chained melded model by extending the marginal replacement
procedure to submodels linked in a chain-like way. The $m$th marginally
replaced submodel is

$\text{p}_{\text{repl},m}(\boldsymbol{\phi},\psi_{m},Y_{m})=\text{p}_{\text{pool}}(\boldsymbol{\phi})\text{p}_{m}(\psi_{m},Y_{m}\mid\boldsymbol{\phi})=\text{p}_{\text{pool}}(\boldsymbol{\phi})\frac{\text{p}_{m}(\phi_{m},\psi_{m},Y_{m})}{\text{p}_{m}(\phi_{m})},$
(6)

where
$\text{p}_{\text{pool}}(\boldsymbol{\phi})=g(\text{p}_{1}(\phi_{1}),\text{p}_{2}(\phi_{2}),\ldots,\text{p}_{m}(\phi_{m}))$
is a pooling function that appropriately summarises all submodel prior
marginals. The second equality in Equation (6) is because of the conditional
independence
$(\psi_{m},Y_{m}\perp\\!\\!\\!\perp\boldsymbol{\phi}_{-m})\mid\phi_{m}$ that
exists due to the chained relationship between submodels. It is important to
note that $\text{p}_{\text{repl},m}(\boldsymbol{\phi},\psi_{m},Y_{m})$ is
defined on a larger parameter space than
$\text{p}_{m}(\phi_{m},\psi_{m},Y_{m})$, as it includes
$\boldsymbol{\phi}_{-m}$.

Define
$\text{p}_{\text{repl},m}(\phi_{m},\psi_{m},Y_{m})=\int\text{p}_{\text{repl},m}(\boldsymbol{\phi},\psi_{m},Y_{m})\text{d}\boldsymbol{\phi}_{-m}$.
Each marginally replaced submodel, as defined in Equation (6), minimises the
following KL divergence444This is shown in Appendix B of the online supplement
to Goudie _et al._ (2019).

$\text{p}_{\text{repl},m}(\phi_{m},\psi_{m},Y_{m})=\operatorname*{arg\,min}_{\text{q}}\left\\{D_{\mathrm{KL}}\left(\text{q}\,\|\,\text{p}_{m}\right)\mid\text{q}(\phi_{m})=\text{p}_{\text{pool}}(\phi_{m})\text{
for all }\phi_{m}\right\\},$ (7)

where
$\text{p}_{\text{pool}}(\phi_{m})=\int\text{p}_{\text{pool}}(\boldsymbol{\phi})\text{d}\boldsymbol{\phi}_{-m}$.
We can thus interpret $\text{p}_{\text{repl},m}(\phi_{m},\psi_{m},Y_{m})$ as a
minimally modified $\text{p}_{m}(\phi_{m},\psi_{m},Y_{m})$ which admits
$\text{p}_{\text{pool}}(\phi_{m})$ as a marginal. Note that it is the
combination of $\text{p}_{\text{repl},m}(\phi_{m},\psi_{m},Y_{m})$ and
$\text{p}_{\text{pool}}(\boldsymbol{\phi}_{-m}\mid\phi_{m})$ that uniquely
determine (6).

We form the chained melded model by taking the Markov combination of the
marginally replaced submodels

$\displaystyle\text{p}_{\text{meld}}(\boldsymbol{\phi},\boldsymbol{\psi},\boldsymbol{Y})$
$\displaystyle=\text{p}_{\text{pool}}(\boldsymbol{\phi})\prod_{m=1}^{M}\text{p}_{\text{repl},m}(\psi_{m},Y_{m}\mid\boldsymbol{\phi})$
(8)
$\displaystyle=\text{p}_{\text{pool}}(\boldsymbol{\phi})\prod_{m=1}^{M}\frac{\text{p}_{m}(\phi_{m},\psi_{m},Y_{m})}{\text{p}_{m}(\phi_{m})}.$
(9)

Rewriting (9) in terms of $\phi_{m\mathrel{{\cap}}m+1}$ for $m=1,\ldots,M-1$
yields

$\displaystyle\text{p}_{\text{meld}}(\boldsymbol{\phi},\boldsymbol{\psi},\boldsymbol{Y})=\text{p}_{\text{pool}}$
$\displaystyle(\boldsymbol{\phi})\frac{\text{p}_{1}(\phi_{1\mathrel{{\cap}}2},\psi_{1},Y_{1})}{\text{p}_{1}(\phi_{1\mathrel{{\cap}}2})}\frac{\text{p}_{M}(\phi_{M-1\mathrel{{\cap}}M},\psi_{M},Y_{M})}{\text{p}_{M}(\phi_{M-1\mathrel{{\cap}}M})}$
(10)
$\displaystyle\times\prod_{m=2}^{M-1}\left(\frac{\text{p}_{m}(\phi_{m-1\mathrel{{\cap}}m},\phi_{m\mathrel{{\cap}}m+1},\psi_{m},Y_{m})}{\text{p}_{m}(\phi_{m-1\mathrel{{\cap}}m},\phi_{m\mathrel{{\cap}}m+1})}\right).$

Finally, we use _chained melded posterior_
$\text{p}_{\text{meld}}(\boldsymbol{\phi},\boldsymbol{\psi}\mid\boldsymbol{Y})\propto\text{p}_{\text{meld}}(\boldsymbol{\phi},\boldsymbol{\psi},\boldsymbol{Y})$
to refer to posterior of the chained melded model conditioned on all data.

###  2.2 Pooled prior

Specifying (9) requires a joint prior for $\boldsymbol{\phi}$. As in Markov
melding we form the joint prior by pooling the marginal priors, selecting a
pooling function $g$ that appropriately represents prior knowledge about the
common quantities. We define $\text{p}_{\text{pool}}(\boldsymbol{\phi})$ as a
generic function of all prior marginals

$\displaystyle\text{p}_{\text{pool}}(\boldsymbol{\phi})$
$\displaystyle=g(\text{p}_{1}(\phi_{1}),\text{p}_{2}(\phi_{2}),\ldots,\text{p}_{M}(\phi_{M}))$
(11)
$\displaystyle=g(\text{p}_{1}(\phi_{1\mathrel{{\cap}}2}),\text{p}_{2}(\phi_{1\mathrel{{\cap}}2},\phi_{2\mathrel{{\cap}}3}),\ldots,\text{p}_{M}(\phi_{M-1\mathrel{{\cap}}M})),$
(12)

because we do not always wish to assume independence between the components of
$\boldsymbol{\phi}$.

Two special cases of Equation (12) are noteworthy. Firstly, if all components
of $\boldsymbol{\phi}$ are independent, then we can form
$\text{p}_{\text{pool}}(\boldsymbol{\phi})$ as the product of $M-1$ standard
pooling functions $h_{m}$ defined in Section 1.2.1

$\displaystyle\text{p}_{\text{pool}}(\boldsymbol{\phi})=\prod_{m=1}^{M-1}\text{p}_{\text{pool},m}(\phi_{m\mathrel{{\cap}}m+1}),$
(13)
$\displaystyle\text{p}_{\text{pool},m}(\phi_{m\mathrel{{\cap}}m+1})=h_{m}(\text{p}_{m}(\phi_{m\mathrel{{\cap}}m+1}),\text{p}_{m+1}(\phi_{m\mathrel{{\cap}}m+1})).$
(14)

A second case, in between complete dependence (12) and independence (14), is
that if
$\text{p}_{m}(\phi_{m-1\mathrel{{\cap}}m},\phi_{m\mathrel{{\cap}}m+1})=\text{p}_{m}(\phi_{m-1\mathrel{{\cap}}m})\text{p}_{m}(\phi_{m\mathrel{{\cap}}m+1})$
then we can define

$\text{p}_{\text{pool}}(\boldsymbol{\phi})=g_{1}(\text{p}_{1}(\phi_{1\mathrel{{\cap}}2}),\ldots,\text{p}_{m}(\phi_{m-1\mathrel{{\cap}}m}))g_{2}(\text{p}_{m}(\phi_{m\mathrel{{\cap}}m+1}),\ldots,\text{p}_{M}(\phi_{M}))$
(15)

without any additional assumptions. That is, if any two consecutive components
of $\boldsymbol{\phi}$ are independent in the submodel containing both of
them, we can divide the pooled prior specification problem into two pooling
functions. The smaller number of arguments to $g_{1}$ and $g_{2}$ make it is
easier to choose appropriate forms for those functions.

Selecting a specific form of $g$ is not trivial given the many choices of
functional form and pooling weights (the latter we discuss momentarily). One
complication is that standard linear and logarithmic pooling, as defined in
Equations (4) and (5), are not immediately applicable when the submodel
marginal distributions consider different quantities. We now propose
extensions to logarithmic, linear, and dictatorial pooling for use in the case
of chained melding.

####  2.2.1 Chained logarithmic pooling

Extending logarithmic pooling for chained Markov melding is straightforward.
We define the logarithmically pooled prior to be

$\begin{gathered}\text{p}_{\text{pool,
log}}(\boldsymbol{\phi})=\frac{1}{K_{\text{log}}(\lambda)}\prod_{m=1}^{M}\text{p}_{m}(\phi_{m})^{\lambda_{m}},\end{gathered}$
(16)

with
$K_{\text{log}}(\lambda)=\int\prod_{m=1}^{M}\text{p}_{m}(\phi_{m})^{\lambda_{m}}\text{d}\boldsymbol{\phi}$
for nonnegative weight vector $\lambda=(\lambda_{1},\ldots,\lambda_{M})$. Note
that (16) does not imply independence between the elements of
$\boldsymbol{\phi}$ because

$\prod_{m=1}^{M}\text{p}_{m}(\phi_{m})^{\lambda_{m}}=\text{p}_{1}(\phi_{1\mathrel{{\cap}}2})^{\lambda_{1}}\prod_{m=2}^{M-1}\left(\text{p}_{m}(\phi_{m-1\mathrel{{\cap}}m},\phi_{m\mathrel{{\cap}}m+1})^{\lambda_{m}}\right)\text{p}_{M}(\phi_{M-1\mathrel{{\cap}}M})^{\lambda_{M}}.$
(17)

When $\lambda_{1}=\lambda_{2}=\ldots=\lambda_{M}=1$ we obtain a special case
which we call product-of-experts (PoE) pooling (Hinton, 2002).

####  2.2.2 Chained linear pooling

Our generalisation of linear pooling to handle marginals of different
quantities is a two step procedure. The first step forms intermediary pooling
densities via standard linear pooling, using appropriate marginals of the
relevant quantity

$\displaystyle\text{p}_{\text{pool},m}(\phi_{m\mathrel{{\cap}}m+1})\,\propto\,\lambda_{m,1}\text{p}_{m}(\phi_{m\mathrel{{\cap}}m+1})+\lambda_{m,2}\text{p}_{m+1}(\phi_{m\mathrel{{\cap}}m+1}),$
(18)

where $\lambda_{m}=(\lambda_{m,1},\lambda_{m,2})$ are nonnegative pooling
weights, and for $m=2,\ldots,M-1$

$\displaystyle\text{p}_{m}(\phi_{m\mathrel{{\cap}}m+1})=\int\text{p}_{m}(\phi_{m-1\mathrel{{\cap}}m},\phi_{m\mathrel{{\cap}}m+1})\text{d}\phi_{m-1\mathrel{{\cap}}m}.$
(19)

For $m=1$ and $m=M$ the relevant marginals are
$\text{p}_{1}(\phi_{1\mathrel{{\cap}}2})$ and
$\text{p}_{M}(\phi_{M-1\mathrel{{\cap}}M})$. In step two we form the pooled
prior as the product of the intermediaries

$\text{p}_{\text{pool,
lin}}(\boldsymbol{\phi})=\frac{1}{K_{\text{lin}}(\lambda)}\prod_{m=1}^{M-1}\text{p}_{\text{pool},m}(\phi_{m\mathrel{{\cap}}m+1}),$
(20)

with
$K_{\text{lin}}(\lambda)=\int\prod_{m=1}^{M-1}\text{p}_{\text{pool},m}(\phi_{m\mathrel{{\cap}}m+1})\text{d}\boldsymbol{\phi}$,
for $\lambda=(\lambda_{1},\ldots,\lambda_{M})$. Clearly, this assumes prior
independence amongst all components of $\boldsymbol{\phi}$ which may be
undesirable, particularly if this independence was not present under one or
more of the submodel priors. We discuss extensions to linear pooling that
enable prior dependence between the components of $\boldsymbol{\phi}$ in
Section 6.

####  2.2.3 Dictatorial pooling

Chained Markov melding does not admit a direct analogue to dictatorial pooling
as defined in Section 1.2.1 because not all submodel prior marginals contain
all common quantities. For example, consider the logarithmically pooled prior
of Equation (16) with, say, the $m^{\prime}$th entry in $\lambda$ set to $1$
and all others set to $0$. This choice of $\lambda$ results in
$\text{p}_{\text{pool}}(\boldsymbol{\phi})=\text{p}(\phi_{m^{\prime}})$, which
is flat for $\boldsymbol{\phi}_{-m^{\prime}}$. It seems reasonable to require
any generalisation of dictatorial pooling to result in a reasonable prior for
all components in $\boldsymbol{\phi}$. Such a generalisation should also
retain the original intention of dictatorial pooling, i.e. ‘ _the
authoritative prior for_ $\phi_{m}$ _is_ $\text{p}_{m}(\phi_{m})$’.

We propose two possible forms of dictatorial pooling that satisfy the
aforementioned criteria. _Partial dictatorial pooling_ enforces a single
submodel prior for the relevant components of $\boldsymbol{\phi}$, with no
restrictions on the pooling of the remaining components; and _complete
dictatorial pooling_ which requires selecting one of the two possible submodel
priors for each component of $\boldsymbol{\phi}$.

Partial dictatorial pooling considers $\text{p}_{m}(\phi_{m})$ as the
authoritative prior for
$\phi_{m}=(\phi_{m-1\mathrel{{\cap}}m},\phi_{m\mathrel{{\cap}}m+1})$. This
results in,

$\displaystyle\text{p}_{\text{pool, dict}}(\boldsymbol{\phi})=$ $\displaystyle
g_{1}\\!\left(\text{p}_{1}(\phi_{1\mathrel{{\cap}}2}),\ldots,\text{p}_{m-1}(\phi_{m-2\mathrel{{\cap}}m-1})\right)$
(21)
$\displaystyle\times\text{p}_{m}(\phi_{m-1\mathrel{{\cap}}m},\phi_{m\mathrel{{\cap}}m+1})$
$\displaystyle\times
g_{2}\\!\left(\text{p}_{m+1}(\phi_{m+1\mathrel{{\cap}}m+2}),\ldots,\text{p}_{M}(\phi_{M-1\mathrel{{\cap}}M})\right),$

where $g_{1}$ and $g_{2}$ are linear or logarithmic pooling functions as
desired555Some care is required if the authoritative submodel is
$\text{p}_{m}$ for $m\in\\{1,2,M-1,M\\}$. If it is taken to be
$m\in\\{1,2\\}$, then $g_{1}$ does not exist, and additionally in the $m=1$
case
$\text{p}_{1}(\phi_{0\mathrel{{\cap}}1},\phi_{1\mathrel{{\cap}}2})\coloneqq\text{p}_{1}(\phi_{1\mathrel{{\cap}}2})$.
The $m\in\\{M-1,M\\}$ cases have analogous definitions..

Complete dictatorial pooling requires the marginal pooled prior for each
component in $\boldsymbol{\phi}$ to be chosen solely on the basis of only one
of the two priors specified for it under the submodels. For $m=1,\ldots,M-1$,
the $m$th marginal of the pooled prior is either

$\text{p}_{\text{pool,
dict}}(\phi_{m\mathrel{{\cap}}m+1})\coloneqq\begin{cases}\text{p}_{m}(\phi_{m\mathrel{{\cap}}m+1})&\text{or}\\\
\text{p}_{m+1}(\phi_{m\mathrel{{\cap}}m+1}).&\end{cases}$ (22)

If two consecutive marginals are chosen to have the same submodel prior then
we instead define

$\displaystyle\text{p}_{\text{pool,
dict}}(\phi_{m-1\mathrel{{\cap}}m})\text{p}_{\text{pool,
dict}}(\phi_{m\mathrel{{\cap}}m+1})$
$\displaystyle=\text{p}_{m}(\phi_{m-1\mathrel{{\cap}}m})\text{p}_{m}(\phi_{m\mathrel{{\cap}}m+1})$
(23)
$\displaystyle=\text{p}_{m}(\phi_{m-1\mathrel{{\cap}}m},\phi_{m\mathrel{{\cap}}m+1}),$

to preserve any dependence between $\phi_{m-1\mathrel{{\cap}}m}$ and
$\phi_{m\mathrel{{\cap}}m+1}$ that may be present under $\text{p}_{m}$. The
complete dictatorially pooled prior is thus

$\text{p}_{\text{pool,
dict}}(\boldsymbol{\phi})=\prod_{m=1}^{M-1}\text{p}_{\text{pool,
dict}}(\phi_{m\mathrel{{\cap}}m+1}),$ (24)

where, subject to the potential modification in Equation (23), the terms in
the product are as defined in Equation (22). For example, if $M=5$ and we wish
to ignore $\text{p}_{2}$ and $\text{p}_{4}$ when constructing the pooled prior
and instead associate $\phi_{1\mathrel{{\cap}}2}$ with $\text{p}_{1}$, both
$\phi_{2\mathrel{{\cap}}3}$ and $\phi_{3\mathrel{{\cap}}4}$ with
$\text{p}_{3}$, and $\phi_{4\mathrel{{\cap}}5}$ with $\text{p}_{5}$, then

$\displaystyle\text{p}_{\text{pool, dict}}(\boldsymbol{\phi})$
$\displaystyle=\text{p}_{1}(\phi_{1\mathrel{{\cap}}2})\text{p}_{3}(\phi_{2\mathrel{{\cap}}3})\text{p}_{3}(\phi_{3\mathrel{{\cap}}4})\text{p}_{5}(\phi_{4\mathrel{{\cap}}5})$
(25)
$\displaystyle=\text{p}_{1}(\phi_{1\mathrel{{\cap}}2})\text{p}_{3}(\phi_{2\mathrel{{\cap}}3},\phi_{3\mathrel{{\cap}}4})\text{p}_{5}(\phi_{4\mathrel{{\cap}}5}).$

####  2.2.4 Pooling weights

Choosing values for the pooling weights is an important step in specifying the
pooled prior. Because appropriate values for the weights depend on the
submodels being pooled and the information available a priori, universal
recommendations are impossible, so we illustrate the impact of different
choices in a straightforward example. It is important that prior predictive
visualisations of the pooled prior are produced (Gabry _et al._ , 2019; Gelman
_et al._ , 2020) to guide the choice of pooling weights and ensure that the
result suitably represents the available information. Figure 3 illustrates how
$\lambda$ and the choice of pooling method impacts
$\text{p}_{\text{pool}}(\boldsymbol{\phi})$ when pooling normal distributions.
Specifically, we consider $M=3$ submodels and pool

$\begin{gathered}\text{p}_{1}(\phi_{1\mathrel{{\cap}}2})=\text{N}(\phi_{1\mathrel{{\cap}}2};\mu_{1},\sigma_{1}^{2}),\quad\text{p}_{3}(\phi_{2\mathrel{{\cap}}3})=\text{N}(\phi_{2\mathrel{{\cap}}3};\mu_{3},\sigma_{3}^{2}),\\\\[3.00003pt]
\text{p}_{2}(\phi_{1\mathrel{{\cap}}2},\phi_{2\mathrel{{\cap}}3})=\text{N}\left(\begin{bmatrix}\phi_{1\mathrel{{\cap}}2}\\\
\phi_{2\mathrel{{\cap}}3}\end{bmatrix}\\!\\!;\begin{bmatrix}\mu_{2,1}\\\
\mu_{2,2}\end{bmatrix}\\!\\!,\begin{bmatrix}\sigma_{2}^{2}&\rho\sigma_{2}^{2}\\\
\rho\sigma_{2}^{2}&\sigma_{2}^{2}\end{bmatrix}\right),\end{gathered}$ (26)

where $\text{N}(\phi;\mu,\sigma^{2})$ is the normal density function with mean
$\mu$ and variance $\sigma^{2}$ (or covariance matrix where appropriate). The
two dimensional density function $\text{p}_{2}$ has an additional parameter
$\rho$, which controls the intra-submodel marginal correlation. We set
$\mu_{1}=-2.5,\mu_{2}=\left[\mu_{2,1}\,\,\mu_{2,2}\right]^{\prime}=\left[0\,\,0\right]^{\prime},\mu_{3}=2.5,\sigma_{1}^{2}=\sigma_{2}^{2}=\sigma_{3}^{2}=1$
and $\rho=0.8$. In the logarithmic case we set $\lambda_{1}=\lambda_{3}$ and
parameterise $\lambda_{2}=1-2\lambda_{1}$, so that
$\lambda_{1}+\lambda_{2}+\lambda_{3}=1$ whilst limiting ourselves to varying
only $\lambda_{1}$. Similarly, in the linear case we set
$\lambda_{1,1}=\lambda_{2,2}=\lambda_{1}$ and
$\lambda_{1,2}=\lambda_{2,1}=1-2\lambda_{1}$. We consider 5 evenly spaced
values of $\lambda_{1}\in[0,0.5]$.

Figure 3: Contour plots of $\text{p}_{\text{pool}}(\boldsymbol{\phi})$ (red)
under logarithmic and linear pooling (left and right column respectively). The
three original densities $\text{p}_{1}(\phi_{1\mathrel{{\cap}}2})$,
$\text{p}_{3}(\phi_{2\mathrel{{\cap}}3})$ and
$\text{p}_{2}(\phi_{1\mathrel{{\cap}}2},\phi_{2\mathrel{{\cap}}3})$ are shown
in blue, with the univariate densities shown on the appropriate axis. The
pooling weight parameter $\lambda_{1}$ is indicated in the plot titles.

For both pooling methods, as the weight $\lambda_{1}$ associated with models
$\text{p}_{1}$ and $\text{p}_{3}$ increases, the relative contributions of
$\text{p}_{1}(\phi_{1\mathrel{{\cap}}2})$ and
$\text{p}_{3}(\phi_{2\mathrel{{\cap}}3})$ increase. Note the lack of
correlation in $\text{p}_{\text{pool}}$ under linear pooling (right column of
Figure 3) due to Equation (20). A large, near-flat plateau is visible in the
$\lambda_{1}=0.25$ and $\lambda_{1}=0.375$ cases, which is a result of the
mixture of four, 2-D normal distributions that linear pooling produces in this
example. The logarithmic pooling process produces a more concentrated prior
for small values of $\lambda_{1}$, and does not result in a priori
independence between $\phi_{1\mathrel{{\cap}}2}$ and
$\phi_{2\mathrel{{\cap}}3}$. Appendix A shows analytically that $\lambda_{2}$
controls the quantity of correlation present in $\text{p}_{\text{pool}}$ in
this setting.

##  3 Posterior estimation

We now present a multi-stage MCMC method for generating samples from the
melded posterior. Whilst the melded posterior is a standard Bayesian posterior
and so can, in principle, be targetted using any suitable Monte Carlo method,
in practice this may be cumbersome or infeasible. More specifically, it may be
feasible to fit each submodel separately using standard methods, but when the
submodels are combined – either through Markov melding, or by expanding the
definition of one submodel to include another – the computation required to
estimate the posterior in a single step poses an insurmountable barrier. In
such settings we can employ multi-stage posterior estimation methods including
Tom _et al._ (2010), Lunn _et al._ (2013), Hooten _et al._ (2019), and Mauff
_et al._ (2020). We propose a multi-stage strategy that uses the chain-like
relationship to both avoid evaluating all submodels simultaneously, and
parallelise the computation required in the first stage to produce posterior
samples in less time than an equivalent sequential method666Appendix B
describes such a sequential MCMC sampler..

We also describe an approximate method, where stage one submodels are
summarised by normal distributions for use in stage two.

We consider the $M=3$ case, as this setting includes both of our examples. Our
approach can be extended to $M>3$ settings, although we anticipate that it is
unlikely to be suitable for large $M$ settings. We discuss some of
difficulties associated with generic, parallel methodology for efficient
posterior sampling in Section 6.

###  3.1 Parallel sampler

Our proposed strategy involves obtaining in stage one samples from submodels 1
and 3 in parallel. Stage two reuses these samples in a Metropolis-within-Gibbs
sampler, which targets the full melded posterior. The stage specific targets
are displayed in Figure 4.

$\text{p}_{1}$$\text{p}_{2}$$\text{p}_{3}$$\phi_{1\mathrel{{\cap}}2}$$\phi_{2\mathrel{{\cap}}3}$$s_{1}$$s_{1}$$s_{2}$
Figure 4: A graphical depiction of the submodels and their shared quantities,
with the parallel sampling strategy overlaid. The stage one ($s_{1}$) targets
are surrounded by blue dashed lines, with the stage two ($s_{2}$) target in
red.

The parallel sampler assumes that the pooled prior decomposes such that

$\text{p}_{\text{pool}}(\boldsymbol{\phi})=\text{p}_{\text{pool},1}(\phi_{1\mathrel{{\cap}}2})\text{p}_{\text{pool},2}(\phi_{1\mathrel{{\cap}}2},\phi_{2\mathrel{{\cap}}3})\text{p}_{\text{pool},3}(\phi_{2\mathrel{{\cap}}3}).$
(27)

All pooled priors trivially satisfy (27) by assuming
$\text{p}_{\text{pool},1}(\phi_{1\mathrel{{\cap}}2})$ and
$\text{p}_{\text{pool},3}(\phi_{2\mathrel{{\cap}}3})$ are improper, flat
distributions. Alternatively we may choose
$\text{p}_{\text{pool},1}(\phi_{1\mathrel{{\cap}}2})=\text{p}_{1}(\phi_{1\mathrel{{\cap}}2})$
and
$\text{p}_{\text{pool},3}(\phi_{2\mathrel{{\cap}}3})=\text{p}_{3}(\phi_{2\mathrel{{\cap}}3})$,
with appropriate adjustments to
$\text{p}_{\text{pool},2}(\phi_{1\mathrel{{\cap}}2},\phi_{2\mathrel{{\cap}}3})$.
This choice targets, in stage one, the subposteriors of $\text{p}_{1}$ and
$\text{p}_{3}$ under their original prior for $\phi_{1\mathrel{{\cap}}2}$ and
$\phi_{2\mathrel{{\cap}}3}$ respectively.

##### Stage one

Two independent, parallel sampling processes occur in stage one. Terms from
the melded model that pertain to $\text{p}_{1}$ and $\text{p}_{3}$ are
isolated

$\displaystyle\text{p}_{\text{meld},1}(\phi_{1\mathrel{{\cap}}2},\psi_{1}\mid
Y_{1})\propto\text{p}_{\text{pool},1}(\phi_{1\mathrel{{\cap}}2})\frac{\text{p}_{1}(\phi_{1\mathrel{{\cap}}2},\psi_{1},Y_{1})}{\text{p}_{1}(\phi_{1\mathrel{{\cap}}2})}$
(28)
$\displaystyle\text{p}_{\text{meld},3}(\phi_{2\mathrel{{\cap}}3},\psi_{3}\mid
Y_{3})\propto\text{p}_{\text{pool},3}(\phi_{2\mathrel{{\cap}}3})\frac{\text{p}_{3}(\phi_{2\mathrel{{\cap}}3},\psi_{3},Y_{3})}{\text{p}_{3}(\phi_{2\mathrel{{\cap}}3})},$
(29)

and targeted using standard MCMC methodology. This produces $N_{1}$ samples
$\\{(\phi_{1\mathrel{{\cap}}2},\psi_{1})_{n}\\}_{n=1}^{N_{1}}$ from
$\text{p}_{\text{meld},1}(\phi_{1\mathrel{{\cap}}2},\psi_{2}\mid Y_{1})$ and
$N_{3}$ samples $\\{(\phi_{2\mathrel{{\cap}}3},\psi_{3})_{n}\\}_{n=1}^{N_{3}}$
from $\text{p}_{\text{meld},3}(\phi_{2\mathrel{{\cap}}3},\psi_{3}\mid Y_{3})$.

##### Stage two

Stage two targets the melded posterior of Equation (9) using a Metropolis-
within-Gibbs sampler, where the proposal distributions are

$\displaystyle\phi_{1\mathrel{{\cap}}2}^{*},\psi_{1}^{*}\mid\phi_{2\mathrel{{\cap}}3},\psi_{2},\psi_{3}$
$\displaystyle\sim\text{p}_{\text{meld},1}(\phi_{1\mathrel{{\cap}}2}^{*},\psi_{1}^{*}\mid
Y_{1})$ (30)
$\displaystyle\phi_{2\mathrel{{\cap}}3}^{*},\psi_{3}^{*}\mid\phi_{1\mathrel{{\cap}}2},\psi_{1},\psi_{2}$
$\displaystyle\sim\text{p}_{\text{meld},3}(\phi_{2\mathrel{{\cap}}3}^{*},\psi_{3}^{*}\mid
Y_{3})$ (31)
$\displaystyle\psi_{2}^{*}\mid\phi_{1\mathrel{{\cap}}2},\phi_{2\mathrel{{\cap}}3},\psi_{1},\psi_{3}$
$\displaystyle\sim\text{q}(\psi_{2}^{*}\mid\psi_{2}),$ (32)

where $\text{q}(\psi_{2}^{*}\mid\psi_{2})$ is a generic proposal distribution
for $\psi_{2}$. In practice, we draw an index $n_{1}^{*}$ uniformly from
$\\{1,\ldots,N_{1}\\}$ and use the corresponding value
$(\phi_{1\mathrel{{\cap}}2}^{*},\psi_{1}^{*})_{n_{1}^{*}}$ as the proposal,
doing likewise for $n_{3}^{*}$ and
$(\phi_{2\mathrel{{\cap}}3}^{*},\psi_{3}^{*})_{n_{3}^{*}}$. The acceptance
probabilities for these updates are

$\displaystyle\alpha((\phi_{1\mathrel{{\cap}}2}^{*},\psi_{1}^{*})_{n_{1}^{*}},(\phi_{1\mathrel{{\cap}}2},\psi_{1})_{n_{1}})$
$\displaystyle=\frac{\text{p}_{\text{pool},2}(\phi_{1\mathrel{{\cap}}2}^{*},\phi_{2\mathrel{{\cap}}3})}{\text{p}_{\text{pool},2}(\phi_{1\mathrel{{\cap}}2},\phi_{2\mathrel{{\cap}}3})}\frac{\text{p}_{2}(\phi_{1\mathrel{{\cap}}2}^{*},\phi_{2\mathrel{{\cap}}3},\psi_{2},Y_{2})}{\text{p}_{2}(\phi_{1\mathrel{{\cap}}2},\phi_{2\mathrel{{\cap}}3},\psi_{2},Y_{2})}\frac{\text{p}_{2}(\phi_{1\mathrel{{\cap}}2},\phi_{2\mathrel{{\cap}}3})}{\text{p}_{2}(\phi_{1\mathrel{{\cap}}2}^{*},\phi_{2\mathrel{{\cap}}3})}$
(33)
$\displaystyle\alpha((\phi_{2\mathrel{{\cap}}3}^{*},\psi_{3}^{*})_{n_{3}^{*}},(\phi_{2\mathrel{{\cap}}3},\psi_{3})_{n_{3}})$
$\displaystyle=\frac{\text{p}_{\text{pool},2}(\phi_{1\mathrel{{\cap}}2},\phi_{2\mathrel{{\cap}}3}^{*})}{\text{p}_{\text{pool},2}(\phi_{1\mathrel{{\cap}}2},\phi_{2\mathrel{{\cap}}3})}\frac{\text{p}_{2}(\phi_{1\mathrel{{\cap}}2},\phi_{2\mathrel{{\cap}}3}^{*},\psi_{2},Y_{2})}{\text{p}_{2}(\phi_{1\mathrel{{\cap}}2},\phi_{2\mathrel{{\cap}}3},\psi_{2},Y_{2})}\frac{\text{p}_{2}(\phi_{1\mathrel{{\cap}}2},\phi_{2\mathrel{{\cap}}3})}{\text{p}_{2}(\phi_{1\mathrel{{\cap}}2},\phi_{2\mathrel{{\cap}}3}^{*})}$
(34) $\displaystyle\alpha(\psi_{2}^{*},\psi_{2})$
$\displaystyle=\frac{\text{p}_{2}(\phi_{1\mathrel{{\cap}}2},\phi_{2\mathrel{{\cap}}3},\psi_{2}^{*},Y_{2})}{\text{p}_{2}(\phi_{1\mathrel{{\cap}}2},\phi_{2\mathrel{{\cap}}3},\psi_{2},Y_{2})}\frac{\text{q}(\psi_{2}\mid\psi_{2}^{*})}{\text{q}(\psi_{2}^{*}\mid\psi_{2})},$
(35)

where $\alpha(x,z)$ denotes the probability associated with a move from $z$ to
$x$. Note that all stage two acceptance probabilities only contain terms from
the second submodel and the pooled prior, and thus do not depend on $\psi_{1}$
or $\psi_{3}$. If a move is accepted then we also store the index, i.e.
$n_{1}^{*}$ or $n_{3}^{*}$, associated with the move, otherwise we store the
current value of the index. The stored indices are used to appropriately
resample $\psi_{1}$ and $\psi_{3}$ from the stage one samples.

###  3.2 Normal approximations to submodel components

Normal approximations are commonly employed to summarise submodels for
subsequent use in more complex models. For example, two-stage meta-analyses
often use a normal distribution centred on each studies’ effect estimate
(Burke _et al._ , 2017). Suppose we employ such an approximation to summarise
the prior and posterior of $\phi_{1\mathrel{{\cap}}2}$ and
$\phi_{2\mathrel{{\cap}}3}$ under $\text{p}_{1}$ and $\text{p}_{3}$
respectively. In addition, assume that (a) such approximations are appropriate
for
$\text{p}_{1}(\phi_{1\mathrel{{\cap}}2}),\text{p}_{1}(\phi_{1\mathrel{{\cap}}2}\mid
Y_{1}),\text{p}_{3}(\phi_{2\mathrel{{\cap}}3})$, and
$\text{p}_{3}(\phi_{2\mathrel{{\cap}}3}\mid Y_{3})$, (b) we are not interested
in $\psi_{1}$ and $\psi_{3}$, and can integrate them out of all relevant
densities, and (c) we employ our second form of dictatorial pooling and choose
$\text{p}_{2}(\phi_{1\mathrel{{\cap}}2},\phi_{2\mathrel{{\cap}}3})$ as the
authoritative prior. The latter two assumptions imply that the melded
posterior of interest is proportional to

$\text{p}_{\text{meld}}(\phi_{1\mathrel{{\cap}}2},\phi_{2\mathrel{{\cap}}3},\psi_{2}\mid\boldsymbol{Y})\,\propto\,\frac{\text{p}_{1}(\phi_{1\mathrel{{\cap}}2}\mid
Y_{1})}{\text{p}_{1}(\phi_{1\mathrel{{\cap}}2})}\text{p}_{2}(\phi_{1\mathrel{{\cap}}2},\phi_{2\mathrel{{\cap}}3},\psi_{2}\mid
Y_{2})\frac{\text{p}_{3}(\phi_{2\mathrel{{\cap}}3}\mid
Y_{3})}{\text{p}_{3}(\phi_{2\mathrel{{\cap}}3})}.$ (36)

Denote the normal approximation of $\text{p}_{1}(\phi_{1\mathrel{{\cap}}2}\mid
Y_{1})$ as
$\widehat{\text{p}}_{1}(\phi_{1\mathrel{{\cap}}2}\mid\widehat{\mu}_{1},\widehat{\Sigma}_{1})$,
which is a normal distribution with mean $\widehat{\mu}_{1}$ and covariance
matrix $\widehat{\Sigma}_{1}$. The corresponding normal approximation of the
prior $\text{p}_{1}(\phi_{1\mathrel{{\cap}}2})$ is
$\widehat{\text{p}}_{1}(\phi_{1\mathrel{{\cap}}2}\mid\widehat{\mu}_{1,0},\widehat{\Sigma}_{1,0})$.
The equivalent approximations for the subposterior and prior of $\text{p}_{3}$
are
$\widehat{\text{p}}_{3}(\phi_{2\mathrel{{\cap}}3}\mid\widehat{\mu}_{3},\widehat{\Sigma}_{3})$
and
$\widehat{\text{p}}_{3}(\phi_{2\mathrel{{\cap}}3}\mid\widehat{\mu}_{3,0},\widehat{\Sigma}_{3,0})$
respectively. Substituting in the approximations and using standard results
for Gaussian density functions (see Bromiley (2003) and Appendix C) results in

$\begin{gathered}\widehat{\text{p}}_{\text{meld}}(\phi_{1\mathrel{{\cap}}2},\phi_{2\mathrel{{\cap}}3},\psi_{2}\mid\boldsymbol{Y})\,\propto\,\widehat{\text{p}}\left((\phi_{1\mathrel{{\cap}}2},\phi_{2\mathrel{{\cap}}3})\mid\widehat{\mu},\,\widehat{\Sigma}\right)\text{p}_{2}(\phi_{1\mathrel{{\cap}}2},\phi_{2\mathrel{{\cap}}3},\psi_{2}\mid
Y_{2}),\end{gathered}$ (37)

where

$\begin{gathered}\widehat{\mu}_{\text{nu}}=\begin{bmatrix}\widehat{\mu}_{1}\\\
\widehat{\mu}_{3}\end{bmatrix}\\!\\!,\quad\widehat{\Sigma}_{\text{nu}}=\begin{bmatrix}\widehat{\Sigma}_{1}&0\\\
0&\widehat{\Sigma}_{3}\end{bmatrix}\\!\\!,\quad\widehat{\mu}_{\text{de}}=\begin{bmatrix}\widehat{\mu}_{1,0}\\\
\widehat{\mu}_{3,0}\end{bmatrix}\\!\\!,\quad\widehat{\Sigma}_{\text{de}}=\begin{bmatrix}\widehat{\Sigma}_{1,0}&0\\\
0&\widehat{\Sigma}_{3,0}\end{bmatrix}\\!\\!,\\\\[6.45831pt]
\widehat{\Sigma}=\left(\widehat{\Sigma}_{\text{nu}}^{-1}-\widehat{\Sigma}_{\text{de}}^{-1}\right)^{-1},\quad\widehat{\mu}=\widehat{\Sigma}\left(\widehat{\Sigma}_{\text{nu}}^{-1}\widehat{\mu}_{\text{nu}}-\widehat{\Sigma}_{\text{de}}^{-1}\widehat{\mu}_{\text{de}}\right).\end{gathered}$
(38)

Standard MCMC methods can be used to sample from the approximate melded
posterior. If instead we opt for product-of-experts pooling, all
$\widehat{\mu}_{\text{de}}$ and $\widehat{\Sigma}_{\text{de}}$ terms disappear
from the parameter definitions in Equation (38).

##  4 An integrated population model for little owls

We now return to the integrated population model (IPM) for the little owls
introduced in Section 1.1.1. Finke _et al._ (2019) consider a number of
variations on the original model of Schaub _et al._ (2006) and Abadi _et al._
(2010): here we consider only the variant from Finke _et al._ (2019) with the
highest marginal likelihood (Model 4 of their online supplement).

Before we detail the specifics of each submodel, we must introduce some
notation. Data and parameters are stratified into two age-groups
$a\in\\{J,A\\}$ where $J$ denotes juvenile owls (less than one year old) and
$A$ adults, two sexes $s\in\\{M,F\\}$, and observations occur annually at
times $t\in\\{1,\ldots,T\\}$, with $T=25$. The sex- and age-specific
probability of an owl surviving from time $t$ to $t+1$ is $\delta_{a,s,t}$,
and the sex-specific probability of a previously captured owl being recaptured
at time $t+1$ is $\pi_{s,t+1}$ so long as the owl is alive at time $t+1$.

###  4.1 Capture recapture: $\text{p}_{1}$

Capture-recapture data pertain to owls that are released at time $t$ (having
been captured and tagged), and then recaptured at time $u=t+1,\dots,T$, or not
recaptured before the conclusion of the study, in which case $u=T+1$. Define
$M_{a,s,t,u}$ as the number of owls of age-group $a$ and sex $s$ last observed
at time $t$, recaptured at time $u$. We aggregate these observations into age-
and sex-specific matrices $\boldsymbol{M}_{a,s}$, with $T$ rows, corresponding
to release times, and $T+1$ columns, corresponding to recapture times. Let
$R_{a,s,t}=\sum_{u=1}^{T+1}M_{a,s,t,u}$ be the number of owls released at time
$t$, i.e. a vector containing the row-wise sum of the entries in
$\boldsymbol{M}_{a,s}$. The recapture times for owls released at time $t$
follow an age- and sex-specific multinomial likelihood

$(M_{a,s,t,1},\ldots,M_{a,s,t,T+1})\sim\text{Multinomial}(R_{a,s,t},\boldsymbol{Q}_{a,s,t}),$
(39)

with probabilities $\boldsymbol{Q}_{a,s,t}=(Q_{a,s,t,1},\ldots,Q_{a,s,t,T+1})$
such that

$Q_{a,s,t,u}=\left\\{\begin{array}[]{ll}0,&\text{for}\,\,u=1,\,\ldots\,,t\\\
\delta_{a,s,t}\pi_{s,u}\prod_{r=t+1}^{u-1}\delta_{a,s,r}\left(1-\pi_{s,r}\right),&\text{for}\,\,u=t+1,\,\ldots\,,T\\\
1-\sum_{r=1}^{T}Q_{a,s,t,r},&\text{if}\,\,u=T+1.\end{array}\right.$ (40)

###  4.2 Count data model: $\text{p}_{2}$

To estimate population abundance, a two level model is used: the first level
models the observed (counted) number of females at each point in time denoted
$y_{t}$, with a second, latent process modelling the total number of females
in population. The observation model is

$y_{t}\mid x_{t}\sim\text{Poisson}\left(x_{t}\right),$ (41)

where we denote the number of juvenile and adult females in the population at
time $t$ as $x_{J,t}$ and $x_{A,t}$ respectively, with
$x_{t}=x_{J,t}+x_{A,t}$. If $\text{sur}_{t}$ adult females survive from time
$t-1$ to time $t$, and $\text{imm}_{t}$ adult females immigrate over the same
time period, then the latent, population level model is

$\displaystyle x_{J,t}\mid x_{t-1},\rho,\delta_{J,F,t-1}$
$\displaystyle\sim\,\,\text{Poisson}\left(x_{t-1}\,\frac{\rho}{2}\,\delta_{J,F,t-1}\right),$
(42) $\displaystyle\text{sur}_{t}\mid x_{t-1},\delta_{A,F,t-1}$
$\displaystyle\sim\,\,\text{Binomial}\left(x_{t-1},\,\delta_{A,F,t-1}\right),$
$\displaystyle\text{imm}_{t}\mid x_{t-1},\eta_{t}$
$\displaystyle\sim\,\,\text{Poisson}\left(x_{t-1}\eta_{t}\right),$
$\displaystyle x_{A,t}$ $\displaystyle=\text{sur}_{t}+\text{imm}_{t},$

where $\eta_{t}$ is the immigration rate. The initial population sizes
$x_{J,1}$ and $x_{A,1}$ have independent discrete uniform priors on
$\\{0,1,\ldots,50\\}$.

###  4.3 Fecundity: $\text{p}_{3}$

The fecundity submodel considers the number of breeding females at time $t$
denoted $N_{t}$, and the number of chicks produced that survive and leave the
nest denoted $n_{t}$. A Poisson model is employed to estimate fecundity
(reproductive rate) $\rho$

$n_{t}\sim\text{Poisson}(N_{t}\rho).$ (43)

###  4.4 Parameterisation and melding quantities

Abadi _et al._ (2010) parameterise the time dependent quantities via linear
predictors, to minimise the number of parameters in the submodels. The
specific parameterisation of Finke _et al._ (2019) we employ is

$\begin{gathered}\text{logit}(\delta_{a,s,t})=\alpha_{0}+\alpha_{1}\mathbb{I}(s=M)+\alpha_{2}\mathbb{I}(a=A),\quad\text{log}(\eta_{t})=\alpha_{6},\\\
\text{logit}(\pi_{s,u})=\alpha_{4}\mathbb{I}(s=M)+\alpha_{5,u},\,\,\text{for}\,\,u=2,\ldots
T,\end{gathered}$ (44)

thus the quantities in common between the submodels are
$\phi_{1\mathrel{{\cap}}2}=(\alpha_{0},\alpha_{2})$ and
$\phi_{2\mathrel{{\cap}}3}=\rho$. To align the notation of this example with
our chained melding notation we define, for all permitted values of $a,s$ and
$t$, $Y_{1}=(\boldsymbol{M}_{a,s})$,
$\psi_{1}=\left(\alpha_{1},\alpha_{4},(\alpha_{5,u})_{u=2}^{T}\right)$;
$Y_{2}=(y_{t})$,
$\psi_{2}=(x_{J,t},\alpha_{6},\text{sur}_{t},\text{imm}_{t})$; and
$Y_{3}=(N_{t},n_{t})$, $\psi_{3}=\varnothing$. Note that the definition of
$\phi_{1\mathrel{{\cap}}2}$ does not include $\alpha_{1}$ as it is male
specific and does not exist in $\text{p}_{2}$. The model variant of Finke _et
al._ (2019) we consider does not include $\alpha_{3}$, and for comparability
we keep the other parameter indices the same.

###  4.5 Priors

We use the priors of Finke _et al._ (2019) for the parameters in each
submodel. Denote
$\boldsymbol{\alpha}=(\alpha_{0},\alpha_{1},\alpha_{2},\alpha_{4},\alpha_{6})$.
In both $\text{p}_{1}$ and $\text{p}_{2}$ the elements of
$\boldsymbol{\alpha}$ are assigned independent $\text{Normal}(0,2^{2})$ priors
truncated to $[-10,10]$. The time varying recapture probabilities
$\alpha_{5,u}$ also have $\text{Normal}(0,2^{2})$ priors truncated to
$[-10,10]$. A $\text{Uniform}(0,10)$ prior is assigned to $\rho$ in
$\text{p}_{2}$ and $\text{p}_{3}$.

To completely specify $\text{p}_{\text{meld}}$ we must choose how to form
$\text{p}_{\text{pool}}(\phi_{1\mathrel{{\cap}}2},\phi_{2\mathrel{{\cap}}3})$.
We form
$\text{p}_{\text{pool}}(\phi_{1\mathrel{{\cap}}2},\phi_{2\mathrel{{\cap}}3})$
using three different pooling methods and estimate the melded posterior in
each case. The first pooling method is product-of-experts (PoE) pooling, which
is logarithmic pooling with $\lambda=(1,1,1)$, and we denote the melded
posterior as $\text{p}_{\text{meld, PoE}}$. We also use logarithmic pooling
with $\lambda=(\frac{1}{2},\frac{1}{2},\frac{1}{2})$, which is denoted
$\text{p}_{\text{meld, log}}$. The final pooling method is linear pooling with
$\lambda=(\frac{1}{2},\frac{1}{2},\frac{1}{2},\frac{1}{2})$, denoted
$\text{p}_{\text{meld, lin}}$.

###  4.6 Posterior estimation

We estimate the melded posterior –
$\text{p}_{\text{meld}}(\boldsymbol{\phi},\boldsymbol{\psi}\mid\boldsymbol{Y})$,
proportional to Equation (9) – using both the parallel sampler (Section 3.1)
and normal approximation (Section 3.2). This allows us to use pre-existing
implementations of the submodels. Specifically, the capture-recapture submodel
is written in BUGS (Lunn _et al._ , 2009) and sampled via rjags (Plummer,
2019). The fecundity submodel is written in Stan (Carpenter _et al._ , 2017)
and sampled via rstan (Stan Development Team, 2021). The count data submodel
is also written in BUGS, and we reuse this implementation in stage two of the
multi-stage sampler via Nimble (de Valpine _et al._ , 2017) and its R
interface (NIMBLE Development Team, 2019). The approximate melded posterior
obtained by Section 3.2 is sampled using rjags. Code and data for this
example, as well as trace plots and numerical convergence measures (Vehtari
_et al._ , 2020) for both stages of the parallel sampling process, are
available in the accompanying online
repository777https://doi.org/10.5281/zenodo.5717455.

###  4.7 Results

Figure 5: Top row: credible intervals for
$\phi_{1\mathrel{{\cap}}2}=(\alpha_{0},\alpha_{2})$ and
$\phi_{2\mathrel{{\cap}}3}=\rho$ from the posterior of the original integrated
population model $\text{p}_{\text{ipm}}$, and the individual subposteriors
from submodels $\text{p}_{1},\text{p}_{2}$, and $\text{p}_{3}$. Bottom row:
credible intervals for the same quantities, but with a different x-axis scale,
from the original IPM (repeated from top row); the chained melded posteriors
using product-of-experts pooling, logarithmic pooling, and linear pooling
denoted $\text{p}_{\text{meld}}$, $\text{p}_{\text{meld, log}}$ and
$\text{p}_{\text{meld, lin}}$; and the melded posterior using the normal
approximation $\widehat{\text{p}}_{\text{meld}}$. Intervals are 50%, 80%, 95%,
and 99% wide.

We empirically validate our methodology and sampler by comparing the melded
posterior samples to a large sample – 6 chains, each containing $1\times
10^{5}$ post-warmup iterations – from the original IPM posterior. Similarity
in the posteriors is expected as the IPM is effectively the joint model we
wish to approximate with the chained melded model. It is simply fortunate,
from a modelling standpoint, that this example’s joint model is easy to
construct and computationally feasible with standard tools. Note that under
logarithmic pooling with $\lambda=(\frac{1}{2},\frac{1}{2},\frac{1}{2})$ the
melded posterior is identical to the original IPM, so any differences are
attributable to the multi-stage sampler. Figure 5 depicts the posterior
credible intervals (Gabry _et al._ , 2021; Kay, 2020) for the common
quantities from the individual submodels, the melded models, and the original
IPM. The top row in Figure 5 indicates that the count data alone
($\text{p}_{2}$) contain minimal information about $\alpha_{0},\alpha_{2}$ and
$\rho$; incorporating the data from the other submodels is essential for
precise estimates.

The multi-stage sampler works well by producing melded posterior estimates
generally similar to the original IPM estimate, and are near identical for
logarithmic pooling. PoE pooling produces the posterior most different from
the original IPM, as it yields a prior for $(\alpha_{0},\alpha_{2})$ that is
more concentrated around zero than the other pooling methods. The lack of
large differences between the melded posteriors that use different pooled
priors indicates that the prior has almost no effect on the posterior. The
similarity of the approximate approach ($\widehat{\text{p}}_{\text{meld}}$ \-
bottom row of Figure 5) to the melding approaches suggests that the normal
approximations are good summaries of the subposteriors, and that the
approximate melding procedure of Section 3.2 is suitable for this example.

##  5 Survival analysis with time varying covariates and uncertain event times

We return now to the respiratory failure example introduced in Section 1.1.2.
Our intention is to illustrate the application of chained Markov melding to an
example of realistic complexity, and explore empirically the importance of
accounting for all sources of uncertainty by comparing chained Markov melding
to equivalent analyses which use only a point estimate summary of the
uncertainty.

There are $i=1,\ldots,N$ individuals in the data set. Each individual is
admitted to the ICU at time $0$, and is discharged or dies at time $C_{i}$.
See Appendix I for information on how the $N=37$ individuals were selected
from MIMIC-III (Johnson _et al._ , 2016).

###  5.1 P/F ratio submodel (B-spline): $\text{p}_{1}$

The first submodel fits a B-spline to the PaO2/FiO2 data to calculate if and
when an individual experiences respiratory failure. Each individual has
PaO2/FiO2 ratio observations $z_{i,j}$ (in units of mmHg) at times $t_{i,j}$,
with $j=1,\ldots,J_{i}$. For each individual denote the vector of observations
$\boldsymbol{z}_{i}=(z_{i,1},\ldots,z_{i,J_{i}})$ and observation times
$\boldsymbol{t}_{i}=(t_{i,1},\ldots,t_{i,J_{i}})$. To improve computational
performance, we standardise the P/F ratio data for each individual such that
$z_{i,j}=\frac{\tilde{z}_{i,j}-\overline{z}_{i}}{\hat{s}_{i}}$, where
$\tilde{z}_{i,j}$ is the underlying unstandardised observation with mean
$\overline{z}_{i}$ and standard deviation $\hat{s}_{i}$. Similarly we rescale
the threshold for respiratory failure:
$\tau_{i}=\frac{300-\overline{z}_{i}}{\hat{s}_{i}}$.

We choose to model the P/F ratio using cubic B-splines and 7 internal knots,
and do not include an intercept column in the spline basis (for background on
B-splines see: Chapter 2 in Hastie and Tibshirani, 1999; and the supplementary
material of Wang and Yan, 2021). The internal knots are evenly spaced between
two additional boundary knots at $\min(\boldsymbol{t}_{i})$ and
$\max(\boldsymbol{t}_{i})$. These choices result in $k=1,\ldots,10$ spline
basis terms per individual, with coefficients $\zeta_{i,k}$ where
$\boldsymbol{\zeta}_{i}=(\zeta_{i,1},\ldots,\zeta_{i,10})$. We denote the
individual specific B-spline basis evaluated at time $t_{i,j}$ as
$B_{i}(t_{i,j})\in[0,\infty)^{10}$ so that the submodel can be written as

$z_{i,j}=\beta_{0,i}+B_{i}(t_{i,j})^{\top}\boldsymbol{\zeta}_{i}+\varepsilon_{i,j}.$
(45)

We employ a weakly informative prior for the intercept
$\beta_{0,i}\sim\text{N}(0,1^{2})$, a heavy tailed distribution for the error
term888P/F data contain many outliers for, amongst many possible reasons,
arterial/venous blood sample mislabelling; incorrectly recorded oxygenation
support information; and differences between sample collection time, lab
result time, and the observation time as recorded in the EHR.
$\varepsilon_{i,j}\sim t_{5}(0,\omega_{i})$, and a weakly informative half-
normal prior for the unknown scale parameter
$\omega_{i}\sim\text{N}_{+}(0,1^{2})$. For the spline basis coefficients we
set $\zeta_{i,1}\sim\text{N}(0,0.1^{2})$, and for $k=2,\ldots,10$ we employ
the random-walk prior $\zeta_{i,k}\sim\text{N}(\zeta_{i,k-1},0.1^{2})$ from
Kharratzadeh (2017).

We identify that a respiratory failure event occurred (which we denote by
$d_{i}=1$) at event time $T_{i}$ if a solution to the following optimisation
problem exists

$T_{i}=\min_{t}\left\\{\tau_{i}=\beta_{0,i}+B_{i}(t)\boldsymbol{\zeta}_{i}\mid
t\in[\max(0,\min(\boldsymbol{t}_{i})),\,\max(\boldsymbol{t}_{i})]\right\\},$
(46)

We attempt to solve Equation 46 using a standard multiple root finder
(Soetaert _et al._ , 2020). If there are no roots then the individual died or
was discharged before respiratory failure occurred so we set $T_{i}=C_{i}$ and
$d_{i}=0$. The relationship between $T_{i}$ and other model coefficients is
displayed in the left hand panel of Figure 6.

2503003504000510152025$t$P/F
ratio$\beta_{0,i}$$\genfrac{}{}{0.0pt}{}{T_{i}}{(d_{i}=1)}$${B_{i}(t)}^{\top}\boldsymbol{\zeta}_{i}$

010200510152025$t$Cumulative
Fluid$\eta_{1,i}^{b}$$\eta_{1,i}^{a}$$\eta_{0,i}$$\kappa_{i}$

Figure 6: Parameters and form for the P/F ratio submodel ($\text{p}_{1}$,
left) and cumulative fluid submodel ($\text{p}_{3}$, right).

###  5.2 Cumulative fluid submodel (piecewise linear) $\text{p}_{3}$

The rate of fluid administration reflects the clinical management of patients
by ICU staff, and hence changes to the rate reflect decisions to change
treatment strategy. We employ a breakpoint regression model to capture the
effect of such decisions, and consider only one breakpoint as this appears
sufficient to fit the observed data. Specifically, we model the 8-hourly
cumulative fluid balance data $x_{i,l}$ (in litres) at times $u_{i,l}$,
$l=1,\ldots,L_{i}$. The cumulative data are derived from the raw fluid
input/output observations, which we detail in Appendix D. We denote the
complete vector of observations by
$\boldsymbol{x}_{i}=(x_{i,1},\ldots,x_{i,L_{i}})$ and times by
$\boldsymbol{u}_{i}=(u_{i,1},\ldots,u_{i,L_{i}})$.

We assume a piecewise linear model with $\eta_{0,i}$ as the value at the
breakpoint at time $\kappa_{i}$, slope $\eta_{1,i}^{b}$ before the breakpoint,
and slope $\eta_{1,i}^{a}$ after the breakpoint. We write this submodel as

$\begin{gathered}x_{i,l}=m_{i}(u_{i,l})+\epsilon_{i,l},\\\
m_{i}(u_{i,l})=\eta_{0,i}+\eta^{b}_{1,i}(u_{i,l}-\kappa_{i})\boldsymbol{1}_{\\{u_{i,l}<\kappa_{i}\\}}+\eta^{a}_{1,i}(u_{i,l}-\kappa_{i})\boldsymbol{1}_{\\{u_{i,l}\geq\kappa_{i}\\}}.\end{gathered}$
(47)

It will be useful to refer to the fitted value of this submodel at arbitrary
time as $m_{i}(t)$. We assume a weakly informative prior for the error term
$\epsilon_{i,l}\sim\text{N}(0,\sigma^{2}_{x,i})$, with individual-specific
error variances $\sigma_{x,i}\sim\text{N}_{+}(0,5^{2})$, and specific,
informative priors for the slope before the breakpoint
$\eta^{b}_{1,i}\sim\text{Gamma}(1.53,0.24)$ and after
$\eta^{a}_{1,i}\sim\text{Gamma}(1.53,0.24)$. An appropriate prior for
$\kappa_{i}$ and $\eta_{0,i}$ is challenging to specify due to the
relationship between the two parameters and the individual-specific support
for $\kappa_{i}$. We address both challenges by reparameterisation, resulting
in a prior for $\kappa_{i}$ that, in the absence of other information, places
the breakpoint in the middle of an individuals ICU stay, and a prior for
$\eta_{0,i}$ that captures the diverse pathways into ICU that an individual
can experience. Details and justifications for all the informative priors are
available in Appendix E. Figure 6 displays the parameters and their
relationship to the fitted regression line.

###  5.3 Survival submodel $\text{p}_{2}$

The rate at which fluid is administered is thought to influence the time to
respiratory failure (Seethala _et al._ , 2017), so we explore this
relationship using a survival model. Individuals experience respiratory
failure ($d_{i}=1$) at time $0<t<C_{i}$, or are censored $(d_{i}=0,t=C_{i})$.
We assume a Weibull hazard with shape parameter $\gamma$ for the event times.
All individuals have baseline (time invariant) covariates $w_{i,a}$,
$a=1,\ldots,A$, with $\boldsymbol{w}_{i}=(1,w_{i,1},\ldots,w_{i,A})$ (i.e.
including an intercept term), and common coefficients
$\boldsymbol{\theta}=(\theta_{0},\ldots,\theta_{A})$. The hazard is assumed to
be influenced by these covariates and the rate of increase
$\frac{\partial}{\partial t}m_{i}(t)$ in the cumulative fluid balance. The
strength of the latter relationship is captured by $\alpha$. Hence, the hazard
is

$\displaystyle h_{i}(t)=\gamma
t^{\gamma-1}\exp\left\\{\boldsymbol{w}_{i}^{\top}\boldsymbol{\theta}+\alpha\frac{\partial}{\partial
t}m_{i}(t)\right\\}$ (48) $\displaystyle\frac{\partial}{\partial
t}m_{i}(t)=\eta^{b}_{1,i}\boldsymbol{1}_{\\{t<\kappa_{i}\\}}+\eta^{a}_{1,i}\boldsymbol{1}_{\\{t\geq\kappa_{i}\\}},$
(49)

The survival function at an individual’s observed event time and status,
$(T_{i},d_{i})$, denoted
$S_{i}(T_{i})=\exp\\{-\int_{0}^{T_{i}}h_{i}(u)\text{d}u\\}$, has an analytic
form which we derive in Appendix F. Hence the likelihood for individual $i$ is

$\text{p}(T_{i},d_{i}\mid\gamma,\boldsymbol{\theta},\alpha,\kappa_{i},\eta_{1,i}^{b},\eta_{1,i}^{a},\boldsymbol{w}_{i})=h_{i}(T_{i})^{d_{i}}S_{i}(T_{i}),\\\
$ (50)

where we suppress the dependence on the parameters on the right hand side for
brevity.

Our priors, which we justify in Appendix G, for the submodel specific
parameters are $\gamma\sim\text{Gamma}(9.05,8.72)$,
$\alpha\sim\text{SkewNormal}(0,0.5,-2)$,
$\theta_{a}\sim\text{SkewNormal}(0,0.5,-1)$, and
$\theta_{0}\sim\text{N}(\hat{E},0.5^{2})$ where $\hat{E}$ is the log of the
crude event rate (Brilleman _et al._ , 2020). We adopt the same priors as the
cumulative fluid balance submodel for $\kappa_{i},\eta_{1,i}^{b}$, and
$\eta_{1,i}^{a}$.

###  5.4 Chained Markov melding details

Figure 7: The P/F ratio data ($Y_{1}$, top row); cumulative fluid data
($Y_{3}$, bottom row); subposterior means and 95% credible intervals for each
of the submodels (black solid lines and grey intervals); and stage one event
times ($T_{i}$, red rug in the top row) for individuals $i=17$ and $29$.

To combine the submodels with chained Markov melding we must define the common
quantities $\phi_{1\mathrel{{\cap}}2}$ and $\phi_{2\mathrel{{\cap}}3}$. We
meld $\text{p}_{1}$ and $\text{p}_{2}$ by treating the derived event times and
indicators $\\{(T_{i},d_{i})\\}_{i=1}^{N}$ under $\text{p}_{1}$ as the
“response”, i.e. event times, in $\text{p}_{2}$. Care is required when
defining $\phi_{1\mathrel{{\cap}}2}$ under $\text{p}_{1}$ as it is a
deterministic function of $\beta_{0,i}$ and $\boldsymbol{\zeta}_{i}$. Define
$\chi_{1,i}=(\beta_{0,i},\boldsymbol{\zeta}_{i})$ and
$\phi_{1\mathrel{{\cap}}2,i}=f(\chi_{1,i})=(T_{i},d_{i})$, where $f$ is the
output from attempting to solve Equation (46), so that
$\phi_{1\mathrel{{\cap}}2}=(f(\chi_{1,i}),\ldots,f(\chi_{1,N}))$. The
parameters shared by Equations (47) and (49) constitute
$\phi_{2\mathrel{{\cap}}3}=(\eta^{b}_{1,i},\eta^{a}_{1,i},\kappa_{i})_{i=1}^{N}$.

To completely align with our chained melding notation we also define, for the
P/F submodel, $Y_{1}=(\boldsymbol{z}_{i},\boldsymbol{t}_{i})_{i=1}^{N}$ and
$\psi_{1}=(\omega_{i})_{i=1}^{N}$, noting that $\psi_{1}$ and
$(\chi_{1,i},\ldots,\chi_{1,N})$ have no components in common. For the
cumulative fluid submodel we define
$Y_{3}=(\boldsymbol{x}_{i},\boldsymbol{u}_{i})_{i=1}^{N}$, and
$\psi_{3}=(\eta_{0,i},\sigma^{2}_{x,i})_{i=1}^{N}$. Finally, for the survival
submodel we define $Y_{2}=(\boldsymbol{w}_{i})_{i=1}^{N}$ and
$\psi_{2}=(\gamma,\boldsymbol{\theta},\alpha)$.

###  5.5 Pooling and estimation

We consider logarithmic pooling with
$\lambda=(\frac{4}{5},\frac{4}{5},\frac{4}{5})$ (any smaller value of
$\lambda$ results in a prior that is so uninformative that it causes
computational problems) and with $\lambda=(1,1,1)$ (Product-of-Experts).
Because the correlation between $\phi_{1\mathrel{{\cap}}2}$ and
$\phi_{2\mathrel{{\cap}}3}$ in
$\text{p}_{2}(\phi_{1\mathrel{{\cap}}2},\phi_{2\mathrel{{\cap}}3})$ is
important, we do not consider linear pooling in this example. Logarithmic
pooling requires estimates of $\text{p}_{1}(\phi_{1\mathrel{{\cap}}2})$ and
$\text{p}_{2}(\phi_{1\mathrel{{\cap}}2},\phi_{2\mathrel{{\cap}}3})$. Because
these are mixed distributions, with both discrete and continuous components,
standard kernel density estimation, as suggested by Goudie _et al._ (2019), is
inappropriate. Instead, we fit appropriate parametric mixture distributions to
the unknown prior marginal distributions. Details for all the parametric
mixture distribution estimates are contained in Appendix H.

We use the parallel multi-stage sampler with
$\text{p}_{\text{pool},1}(\phi_{1\mathrel{{\cap}}2})=\text{p}_{1}(\phi_{1\mathrel{{\cap}}2}),\,\,\text{p}_{\text{pool},3}(\phi_{2\mathrel{{\cap}}3})=\text{p}_{3}(\phi_{2\mathrel{{\cap}}3})$
and
$\text{p}_{\text{pool},2}(\phi_{1\mathrel{{\cap}}2},\phi_{2\mathrel{{\cap}}3})=\text{p}_{\text{pool}}(\boldsymbol{\phi})\mathop{/}\left(\text{p}_{1}(\phi_{1\mathrel{{\cap}}2})\text{p}_{3}(\phi_{2\mathrel{{\cap}}3})\right)$.
That is, in stage one we target the subposteriors
$\text{p}_{1}(\phi_{1\mathrel{{\cap}}2},\psi_{1}\mid Y_{1})$ and
$\text{p}_{3}(\phi_{2\mathrel{{\cap}}3},\psi_{3}\mid Y_{3})$; in stage two we
target the full melded model. It is crucial for the convergence of our multi-
stage sampler that $\phi_{1\mathrel{{\cap}}2}$ and $\phi_{2\mathrel{{\cap}}3}$
are updated _individual-at-a-time_. Appendix K contains the details of this
scheme. The stage one subposteriors are sampled using Stan, using 5 chains
with $10^{3}$ warm-up iterations and $10^{4}$ post warm-up iterations. We use
Stan to sample $\psi_{2}$ where, in every MH-within-Gibbs step, we run Stan
for 9 warm-up iterations and 1 post warm-up iteration999We also initialise
Stan at the previous value of $\psi_{2}$, and disable all adaptive procedures
as the default (identity) mass matrix and step size are suitable for this
example.. We run 5 chains of $10^{4}$ iterations for all stage two targets.
Visual and numerical diagnostics (Vehtari _et al._ , 2020) are assessed and
are available in the repository accompanying this
paper101010https://doi.org/10.5281/zenodo.5717455.

###  5.6 Results

We first inspect the subposterior fitted values for $\text{p}_{1}$ and
$\text{p}_{3}$. The top row of Figure 7 displays the P/F data, the fitted
submodel, and derived event times for individuals $i=17$ and $29$. The spline
appears to fit the raw P/F data well, with the heavy tailed error term
accounting for the larger deviations away from the fitted value. It is
interesting to see the relatively wide, multimodal distribution for
$(T_{29},d_{29})$ (there is a second mode at $(T_{29}=C_{29},d_{29}=0)$ and
for other individuals not shown here). The bottom row of Figure 7 displays the
cumulative fluid data and the fitted submodel, with the little noise in the
data resulting in minimal uncertainty about the fitted value and a
concentrated subposterior distribution.

To assess the importance of fully accounting for the uncertainty in
$\phi_{1\mathrel{{\cap}}2}$ and $\phi_{2\mathrel{{\cap}}3}$, we compare the
posterior for $\psi_{2}$ obtained using the chained melding approach with the
posterior obtained by fixing $\phi_{1\mathrel{{\cap}}2}$ and
$\phi_{2\mathrel{{\cap}}3}$. Plugging in a point estimate reflects common
applied statistical practice when combining submodels, particularly when a
distributional approximation is difficult to obtain (as it is for
$\text{p}_{1}(\phi_{1\mathrel{{\cap}}2}\mid Y_{1})$). Additionally, standard
survival models and software typically do not permit uncertainty in event
times and indicators, rendering such a plug-in approach necessary.

Specifically, we fix $\phi_{1\mathrel{{\cap}}2}$ to the median value111111For
each individual the samples of $(T_{i},d_{i})_{i=1}^{N}$-pairs are sorted by
$T_{i}$, and the $\lfloor\frac{N}{2}\rfloor$th tuple
$(\widehat{T}_{i},\widehat{d}_{i})$ is chosen as the median. for each
individual under $\text{p}_{1}(\phi_{1\mathrel{{\cap}}2}\mid Y_{1})$ and
denote it $\widehat{\phi}_{1\mathrel{{\cap}}2}$, and use the subposterior mean
of $\text{p}_{3}(\phi_{2\mathrel{{\cap}}3}\mid Y_{3})$ denoted
$\widehat{\phi}_{2\mathrel{{\cap}}3}$. With these fixed values we sample
$\text{p}(\psi_{2}\mid\widehat{\phi}_{1\mathrel{{\cap}}2},\widehat{\phi}_{2\mathrel{{\cap}}3},Y_{2})$.
We also compare the melded posterior to the submodel marginal prior
$\text{p}_{2}(\psi_{2})$, but we note that this comparison is difficult to
interpret, as the melding process alters the prior for $\psi_{2}$. Figure 8
displays the aforementioned densities for
$(\theta_{3},\theta_{17},\gamma,\alpha)\subset\psi_{2}$, with
$(\theta_{3},\theta_{17})$ chosen as they exhibit the greatest sensitivity to
the fixing of $\phi_{1\mathrel{{\cap}}2}$ and $\phi_{2\mathrel{{\cap}}3}$. For
the baseline coefficients ($\theta_{3},\theta_{17}$) the chained melding
posterior differs slightly in location from
$\text{p}(\psi_{2}\mid\widehat{\phi}_{1\mathrel{{\cap}}2},\widehat{\phi}_{2\mathrel{{\cap}}3},Y_{2})$,
with a small increase in uncertainty. A more pronounced change is visible for
$\alpha$, where the melding process has added a notable degree of uncertainty
and shifted the posterior leftwards.

Figure 8: Density estimates for a subset of $\psi_{2}$. The submodel marginal
prior $\text{p}_{2}(\psi_{2})$ is shown as the grey dotted line (note that
this is not the marginal prior under the melded model). The figure also
contains the subposteriors obtained from chained melding using PoE pooling
(red, solid line) and logarithmic pooling (blue, solid line), as well as the
posterior using the fixed values
$\text{p}(\psi_{2}\mid\widehat{\phi}_{1\mathrel{{\cap}}2},\widehat{\phi}_{2\mathrel{{\cap}}3},Y_{2})$
(black, dashed line).

To investigate which part of the melding process causes this change in the
posterior of $\alpha$, we consider fixing either one of
$\phi_{1\mathrel{{\cap}}2}$ and $\phi_{2\mathrel{{\cap}}3}$ to their
respective point estimates. That is, we employ Markov melding as described in
Section 1.2, using either logarithmic or PoE pooling, to obtain
$\text{p}_{\text{meld}}(\alpha\mid\widehat{\phi}_{1\mathrel{{\cap}}2},Y_{2},Y_{3})$
and
$\text{p}_{\text{meld}}(\alpha\mid\widehat{\phi}_{2\mathrel{{\cap}}3},Y_{1},Y_{2})$.
Figure 9 displays the same distributions for $\alpha$ as Figure 8, and adds
the posteriors obtained using one fixed value
($\widehat{\phi}_{1\mathrel{{\cap}}2}$ or
$\widehat{\phi}_{2\mathrel{{\cap}}3}$) whilst melding the other non-fixed
parameter.

Evident for both choices of pooling is the importance of incorporating the
uncertainty in $\phi_{1\mathrel{{\cap}}2}$. This is expected given the large
uncertainty and multimodal nature of $\phi_{1\mathrel{{\cap}}2}$ compared to
$\phi_{2\mathrel{{\cap}}3}$ (see Figure 7). We suspect that it is the
multimodality in $\text{p}_{1}(\phi_{1\mathrel{{\cap}}2}\mid Y_{1})$ that
produces the shift in posterior mode of $\phi_{1\mathrel{{\cap}}2}$, with the
width of $\text{p}_{1}(\phi_{1\mathrel{{\cap}}2}\mid Y_{1})$ affecting the
increase in uncertainty. Because we prefer the chained melded posterior, under
either pooling method, for its full accounting of uncertainty we conclude that
$\text{p}(\alpha\mid\widehat{\phi}_{1\mathrel{{\cap}}2},\widehat{\phi}_{2\mathrel{{\cap}}3},Y_{2})$
is both overconfident and biased.

Figure 9: Median (vertical line), 50%, 80%, 95%, and 99% credible intervals
(least transparent to most transparent) for $\alpha$. The marginal prior
(grey, top row) and posterior using fixed
$\widehat{\phi}_{1\mathrel{{\cap}}2}$ and
$\widehat{\phi}_{2\mathrel{{\cap}}3}$ (black, bottom row) are as in Figure 8.
For the chained melded posteriors (red and blue, rows 2 and 3) and the melded
posteriors (red and blue, rows 4 – 7), the tick label on the y-axis denotes
the type of pooling used, and which of $\phi_{1\mathrel{{\cap}}2}$ and/or
$\phi_{2\mathrel{{\cap}}3}$ are fixed.

The marginal changes to the components of $\psi_{2}$ visible in Figure 8
appear small, however the cumulative effect of such changes becomes apparent
when inspecting the posterior of the survival function. Figure 10 displays the
model-based, mean survival function under the melded posterior (using PoE
pooling), and corresponding draws of $\phi_{1\mathrel{{\cap}}2}$ converted
into survival curves using the Kaplan-Meier estimator. Also shown are the
Kaplan-Meier estimate of $\widehat{\phi}_{1\mathrel{{\cap}}2}$ and the mean
survival function computed using
$\text{p}(\psi_{2}\mid\widehat{\phi}_{1\mathrel{{\cap}}2},\widehat{\phi}_{2\mathrel{{\cap}}3},Y_{2})$.
The posterior survival functions differ markedly, with the 95% intervals
overlapping only for small values of time. It is also interesting to see that
$\widehat{\phi}_{1\mathrel{{\cap}}2}$, despite being a reasonable point
estimate of $\text{p}_{1}(\phi_{1\mathrel{{\cap}}2}\mid Y_{1})$, is not very
likely under the melded posterior. Figure 10 also suggests that the Weibull
hazard is insufficiently flexible for this example. We discuss the
complexities of other hazards in Section 6.

Figure 10: Survival curves and mean survival function at time $t$. The red,
stepped lines are draws of $\phi_{1\mathrel{{\cap}}2}$ from the melded
posterior using PoE pooling, converted into survival curves via the Kaplan-
Meier estimator. The smooth red line and interval (posterior mean and 95%
credible interval) denote the model-based, mean survival function obtained
from the melded posterior (PoE pooling) values of $\psi_{2}$ and
$\phi_{2\mathrel{{\cap}}3}$. The blue dashed line is the Kaplan-Meier estimate
of $\widehat{\phi}_{1\mathrel{{\cap}}2}$, and the blue solid line and interval
are the corresponding model-based estimate from
$\text{p}(\psi_{2}\mid\widehat{\phi}_{1\mathrel{{\cap}}2},\widehat{\phi}_{2\mathrel{{\cap}}3},Y_{2})$.

##  6 Conclusion

This paper introduces the chained Markov melded model. In doing so we make
explicit the notion of submodels related in a chain-like way, describe a
generic methodology for joining together any number of such submodels and
illustrate its application with our examples. Our examples also demonstrate
the importance of quantifying the uncertainty when joining submodels; not
doing so can produce biased, over-confident inference.

We have introduced extensions to linear and logarithmic pooling to marginals
of different but overlapping quantities. Linear pooling, introduced in Section
1.2.1, could be extended to induce dependence between the components of
$\boldsymbol{\phi}$ using multivariate or vine copulas (Kurowicka and Joe,
2011; Nelsen, 2006), or other techniques (Lin _et al._ , 2014). Copula methods
are particularly appealing as, depending on the choice of copula, they yield
computationally cheap to evaluate expressions for the density function, are
easy to sample, and induce correlation between an arbitrary number of
marginals.

Our parallel multi-stage sampler currently only considers $M=3$ submodels,
rather than the fully generic definition of chained Markov melding in Equation
(10). Whilst we anticipate needing more complex methods in large $M$ settings,
the value of $M$ at which the performance of our multi-stage sampler becomes
unacceptable will depend on the specific submodels and data under
consideration. A general method would consider a large and arbitrary number of
submodels in a chain, and initially split the chain into more/fewer pieces
depending on the computational resources available. Designing such a method is
complex, as it would have to:

  * •

avoid requiring the inverse of any component of $\boldsymbol{\phi}$ with a
noninvertible definition,

  * •

estimate the relative cost of sampling each submodel’s subposterior, to split
the chain of submodels into steps/jobs of approximately the same computational
cost,

  * •

decide the order in which pieces of the chain are combined.

These are substantial challenges. It may be possible to use combine the ideas
in Lindsten _et al._ (2017) and Kuntz _et al._ (2021), who propose a parallel
Sequential Monte Carlo method, with the aforementioned constraints to obtain a
generic methodology. Ideally we would retain the ability to use existing
implementations of the submodels, however the need to recompute the weights of
the particles, and hence reevaluate previously considered submodels, may
preclude this requirement. Our currently sampler is also sensitive to large
differences in location or scale of the target distribution between the
stages. The impact of these differences can be ameliorated using Manderson and
Goudie (2021), and, more generally, Sequential Monte Carlo samplers are likely
to perform better in these settings.

Our chained Markov melding methodology is general and permits any form of
uncertainty in the common quantities. In Section 5 we use our chained melded
model to incorporate uncertainty in the event times and indicators into a
survival submodel. Some specific forms of uncertainty in the event times have
been considered in previous work. These include Wang _et al._ (2020), who
consider uncertain event times arising from record linkage, where the event
time is assumed to be one of a finite number of event times arising from the
record linkage; and Giganti _et al._ (2020), Oh _et al._ (2018), and Oh _et
al._ (2021), who leverage external validation data to account for measurement
error in the event time. However, the general and Bayesian nature of our
methodology readily facilitates any form of uncertainty in the event times and
the event indicators; uncertainty in the latter is not considered the cited
papers.

The example in Section 5 has three more interesting aspects to discuss.
Firstly, the P/F ratio data used in the first submodel is obtained by finding
all blood gas measurements from arterial blood samples. Approximately $20\%$
of the venous/arterial labels are missing. In these instances a logistic
regression model, fit by the MIMIC team121212The coefficients, classification
threshold, and the imputation used in the case of missing data are supplied in
the blood-gasses.sql file in the GitHub repository accompanying this paper. No
other information is available about this model (the data used to produce the
coefficients, and the performance of the fitted model)., is used to predict
the missing label based on other covariates. It is theoretically possible to
refit the model in a Bayesian framework and use the chained melded model to
incorporate the uncertainty in the predicted sample label – adding another
‘link’ to the chain.

Secondly, the application of our multi-stage sampler to this example is
similar to the two-stage approach used for joint longitudinal and time-to-
event models (see Mauff _et al._ , 2020 for a description of this approach).
In the two-stage approach, the longitudinal model is fit using standard MCMC
methods in stage one, and the samples are reused in stage two when considering
the time-to-event data. This can significantly reduce the computational effort
required to fit the joint model. However, unlike our multi-stage sampler, the
typical two-stage approach does not target the full posterior distribution,
which can lead to biased estimates (though Mauff _et al._ (2020) extend the
typical two-stage approach to reduce this bias).

Thirdly, we observe a lack of flexibility the baseline hazard, visible in
Figure 10. More complex hazards could be employed, e.g. modelling the
(log-)hazard using a (penalised) B-spline (Rosenberg, 1995; Royston and
Parmar, 2002; Rutherford _et al._ , 2015). However, this increased flexibility
precludes an analytic form for the survival function. Whilst numerical
integration is possible it is not trivial, particularly when the hazard is
discontinuous, as our hazard is at the breakpoint. Splines also have more
coefficients than the single parameter of the Weibull hazard. Identifiability
issues arise with a small number of individuals, many of whom are censored,
and are compounded when there are a relatively large number of other
parameters $(\alpha,\boldsymbol{\theta})$. Whilst we do not believe these
costs are worth incurring for our example, for settings with a larger number
of patients and more complicated longitudinal submodels the increased
flexibility may be vital.

### Acknowledgements

We thank Sarah L Cowan for assistance in understanding respiratory failure,
and Anne Presanis and Brian Tom for many helpful discussions about the
methodological aspects of this work. This work was supported by The Alan
Turing Institute under the UK Engineering and Physical Sciences Research
Council (EPSRC) [EP/N510129/1] and the UK Medical Research Council [programme
code MC_UU_00002/2].

## Appendices

###  A Log pooling Gaussian densities

We can exactly compute $\text{p}_{\text{pool}}$ when logarithmically pooling
Gaussian densities. Noting that, in the one dimensional case,
$\text{N}(\phi;\mu,\sigma^{2})^{\lambda_{m}}=\text{N}(\phi;\mu,\frac{\sigma^{2}}{\lambda_{m}})$,
we use the results of Bromiley (2003) and write

$\displaystyle\Omega_{1}=\left(\begin{bmatrix}\lambda_{1}^{-1}&0\\\
0&\lambda_{3}^{-1}\end{bmatrix}\begin{bmatrix}\sigma_{1}^{2}&0\\\
0&\sigma_{3}^{2}\end{bmatrix}\right)^{-1},\;\;\Omega_{2}=\left(\begin{bmatrix}\lambda_{2}^{-1}&0\\\
0&\lambda_{2}^{-1}\end{bmatrix}\begin{bmatrix}\sigma_{2}^{2}&\rho\sigma_{2}^{2}\\\
\rho\sigma_{2}^{2}&\sigma_{2}^{2}\end{bmatrix}\right)^{-1}$ (51)
$\displaystyle\Sigma_{\text{log}}=\left(\Omega_{1}+\Omega_{2}\right)^{-1},\quad\mu_{\text{log}}=\Sigma_{\text{log}}\left(\Omega_{1}\begin{bmatrix}\mu_{1}\\\
\mu_{3}\end{bmatrix}+\Omega_{2}\begin{bmatrix}\mu_{2,1}\\\
\mu_{2,1}\end{bmatrix}\right),$ (52)

hence
$\text{p}_{\text{pool}}(\phi_{1\mathrel{{\cap}}2},\phi_{2\mathrel{{\cap}}3})=\text{N}((\phi_{1\mathrel{{\cap}}2},\phi_{2\mathrel{{\cap}}3});\,\mu_{\text{log}},\,\Sigma_{\text{log}})$.
The choice of $\lambda_{2}$ is critical; by controlling the contribution of
$\text{p}_{2}$ to $\text{p}_{\text{pool}}$, $\lambda_{2}$ controls the degree
of correlation present in the latter. The left hand column of Figure 3
illustrates this phenomena. When
$\lambda_{1}=\lambda_{3}=0\implies\lambda_{2}=1$, all correlation in
$\text{p}_{2}$ is present in $\text{p}_{\text{pool}}$. The correlation
decreases for increasing values of $\lambda_{1}$ until
$\lambda_{1}=\lambda_{3}=0.5\implies\lambda_{2}=0$, where no correlation
persists.

###  B Sequential sampler

$\text{p}_{1}$$\text{p}_{2}$$\text{p}_{3}$$\phi_{1\mathrel{{\cap}}2}$$\phi_{2\mathrel{{\cap}}3}$$s_{1}$$s_{2}$$s_{3}$
Figure 11: A DAG of the submodels and their common quantities, with the
sequential sampling strategy overlaid. The stage one target ($s_{1}$) is
encapsulated in the light blue dashed line; stages two and three
($s_{2},s_{3}$) are in dark blue and red respectively.

Figure 11 depicts graphically the strategy employed by the sequential sampler.
The sequential sampler assumes that the pooled prior decomposes such that

$\text{p}_{\text{pool}}(\boldsymbol{\phi})=\text{p}_{\text{pool},1}(\phi_{1\mathrel{{\cap}}2})\text{p}_{\text{pool},2}(\phi_{1\mathrel{{\cap}}2},\phi_{2\mathrel{{\cap}}3})\text{p}_{\text{pool},3}(\phi_{1\mathrel{{\cap}}2},\phi_{2\mathrel{{\cap}}3}).$
(53)

This is necessary to avoid sampling all components of $\boldsymbol{\phi}$ in
the first stage. All pooled priors trivially satisfy (53), as we can assume
all but
$\text{p}_{\text{pool},3}(\phi_{1\mathrel{{\cap}}2},\phi_{2\mathrel{{\cap}}3})$
are improper, flat distributions. However, including some portion of the
pooled prior in each stage of the sampler can improve performance, and
eliminate computational instabilities when submodel likelihoods contain little
information.

####  B.1 Stage one

Stage one of the sequential sampler targets

$\text{p}_{\text{meld},1}(\phi_{1\mathrel{{\cap}}2},\psi_{1}\mid
Y_{1})\propto\text{p}_{\text{pool},1}(\phi_{1\mathrel{{\cap}}2})\frac{\text{p}_{1}(\phi_{1\mathrel{{\cap}}2},\psi_{1},Y_{1})}{\text{p}_{1}(\phi_{1\mathrel{{\cap}}2})},$
(54)

using a generic proposal kernel for both $\phi_{1\mathrel{{\cap}}2}$ and
$\psi_{1}$. The corresponding acceptance probability for a proposed update
from $(\phi_{1\mathrel{{\cap}}2},\psi_{1})$ to
$(\phi_{1\mathrel{{\cap}}2}^{*},\psi_{1}^{*})$ is

$\alpha((\phi_{1\mathrel{{\cap}}2}^{*},\psi_{1}^{*}),(\phi_{1\mathrel{{\cap}}2},\psi_{1}))=\\\
\frac{\text{p}_{\text{pool},1}(\phi_{1\mathrel{{\cap}}2}^{*})}{\text{p}_{\text{pool},1}(\phi_{1\mathrel{{\cap}}2})}\frac{\text{p}_{1}(\phi_{1\mathrel{{\cap}}2}^{*},\psi_{1}^{*},Y_{1})\text{p}_{1}(\phi_{1\mathrel{{\cap}}2})}{\text{p}_{1}(\phi_{1\mathrel{{\cap}}2},\psi_{1},Y_{1})\text{p}_{1}(\phi_{1\mathrel{{\cap}}2}^{*})}\frac{\text{q}(\phi_{1\mathrel{{\cap}}2},\psi_{1}\mid\phi_{1\mathrel{{\cap}}2}^{*},\psi_{1}^{*})}{\text{q}(\phi_{1\mathrel{{\cap}}2}^{*},\psi_{1}^{*}\mid\phi_{1\mathrel{{\cap}}2},\psi_{1})}.$
(55)

####  B.2 Stage two

The stage two target augments the stage one target by including the second
submodel, corresponding prior marginal distribution, and an additional pooled
prior term

$\text{p}_{\text{meld},2}(\phi_{1\mathrel{{\cap}}2},\phi_{2\mathrel{{\cap}}3},\psi_{1},\psi_{2}\mid
Y_{1},Y_{2})\propto\\\
\text{p}_{\text{pool},1}(\phi_{1\mathrel{{\cap}}2})\text{p}_{\text{pool},2}(\phi_{1\mathrel{{\cap}}2},\phi_{2\mathrel{{\cap}}3})\frac{\text{p}_{1}(\phi_{1\mathrel{{\cap}}2},\psi_{1},Y_{1})}{\text{p}_{1}(\phi_{1\mathrel{{\cap}}2})}\frac{\text{p}_{2}(\phi_{1\mathrel{{\cap}}2},\phi_{2\mathrel{{\cap}}3},\psi_{2},Y_{2})}{\text{p}_{2}(\phi_{1\mathrel{{\cap}}2},\phi_{2\mathrel{{\cap}}3})}.$
(56)

A Metropolis-within-Gibbs strategy is employed, where the stage one samples
are used as a proposal for $\phi_{1\mathrel{{\cap}}2}$, whilst a generic
proposal kernel is used for $\psi_{2}$ and $\phi_{2\mathrel{{\cap}}3}$. Thus
the proposal distributions for $\phi_{1\mathrel{{\cap}}2}^{*}$ and
$(\phi_{2\mathrel{{\cap}}3}^{*},\psi_{2}^{*})$ are

$\displaystyle\phi_{1\mathrel{{\cap}}2}^{*},\psi_{1}^{*}\mid\phi_{2\mathrel{{\cap}}3},\psi_{2}$
$\displaystyle\sim\text{p}_{\text{meld},1}(\phi_{1\mathrel{{\cap}}2}^{*},\psi_{1}^{*}\mid
Y_{1})$ (57)
$\displaystyle\phi_{2\mathrel{{\cap}}3}^{*},\psi_{2}^{*}\mid\phi_{1\mathrel{{\cap}}2},\psi_{1}$
$\displaystyle\sim\text{q}(\phi_{2\mathrel{{\cap}}3}^{*},\psi_{2}^{*}\mid\phi_{2\mathrel{{\cap}}3},\psi_{2}).$
(58)

The acceptance probability for this proposal strategy is

$\displaystyle\alpha$
$\displaystyle\\!\\!\begin{multlined}\left((\phi_{1\mathrel{{\cap}}2}^{*},\psi_{1}^{*}),(\phi_{1\mathrel{{\cap}}2},\psi_{1})\right)=\\\
\frac{\text{p}_{\text{pool},2}(\phi_{1\mathrel{{\cap}}2}^{*},\phi_{2\mathrel{{\cap}}3})}{\text{p}_{\text{pool},2}(\phi_{1\mathrel{{\cap}}2},\phi_{2\mathrel{{\cap}}3})}\frac{\text{p}_{2}(\phi_{1\mathrel{{\cap}}2}^{*},\phi_{2\mathrel{{\cap}}3},\psi_{2},Y_{2})\text{p}_{2}(\phi_{1\mathrel{{\cap}}2},\phi_{2\mathrel{{\cap}}3})}{\text{p}_{2}(\phi_{1\mathrel{{\cap}}2},\phi_{2\mathrel{{\cap}}3},\psi_{2},Y_{2})\text{p}_{2}(\phi_{1\mathrel{{\cap}}2}^{*},\phi_{2\mathrel{{\cap}}3})}\end{multlined}\left((\phi_{1\mathrel{{\cap}}2}^{*},\psi_{1}^{*}),(\phi_{1\mathrel{{\cap}}2},\psi_{1})\right)=\\\
\frac{\text{p}_{\text{pool},2}(\phi_{1\mathrel{{\cap}}2}^{*},\phi_{2\mathrel{{\cap}}3})}{\text{p}_{\text{pool},2}(\phi_{1\mathrel{{\cap}}2},\phi_{2\mathrel{{\cap}}3})}\frac{\text{p}_{2}(\phi_{1\mathrel{{\cap}}2}^{*},\phi_{2\mathrel{{\cap}}3},\psi_{2},Y_{2})\text{p}_{2}(\phi_{1\mathrel{{\cap}}2},\phi_{2\mathrel{{\cap}}3})}{\text{p}_{2}(\phi_{1\mathrel{{\cap}}2},\phi_{2\mathrel{{\cap}}3},\psi_{2},Y_{2})\text{p}_{2}(\phi_{1\mathrel{{\cap}}2}^{*},\phi_{2\mathrel{{\cap}}3})}$
(61) $\displaystyle\alpha$
$\displaystyle\\!\\!\begin{multlined}\left((\phi_{2\mathrel{{\cap}}3}^{*},\psi_{2}^{*}),(\phi_{2\mathrel{{\cap}}3},\psi_{2})\right)=\\\
\frac{\text{p}_{\text{pool},2}(\phi_{1\mathrel{{\cap}}2},\phi_{2\mathrel{{\cap}}3}^{*})}{\text{p}_{\text{pool},2}(\phi_{1\mathrel{{\cap}}2},\phi_{2\mathrel{{\cap}}3})}\frac{\text{p}_{2}(\phi_{1\mathrel{{\cap}}2},\phi_{2\mathrel{{\cap}}3}^{*},\psi_{2}^{*},Y_{2})\text{p}_{2}(\phi_{1\mathrel{{\cap}}2},\phi_{2\mathrel{{\cap}}3})}{\text{p}_{2}(\phi_{1\mathrel{{\cap}}2},\phi_{2\mathrel{{\cap}}3},\psi_{2},Y_{2})\text{p}_{2}(\phi_{1\mathrel{{\cap}}2},\phi_{2\mathrel{{\cap}}3}^{*})}\frac{\text{q}(\phi_{2\mathrel{{\cap}}3},\psi_{2}\mid\phi_{2\mathrel{{\cap}}3}^{*},\psi_{2}^{*})}{\text{q}(\phi_{2\mathrel{{\cap}}3}^{*},\psi_{2}^{*}\mid\phi_{2\mathrel{{\cap}}3},\psi_{2})}.\end{multlined}\left((\phi_{2\mathrel{{\cap}}3}^{*},\psi_{2}^{*}),(\phi_{2\mathrel{{\cap}}3},\psi_{2})\right)=\\\
\frac{\text{p}_{\text{pool},2}(\phi_{1\mathrel{{\cap}}2},\phi_{2\mathrel{{\cap}}3}^{*})}{\text{p}_{\text{pool},2}(\phi_{1\mathrel{{\cap}}2},\phi_{2\mathrel{{\cap}}3})}\frac{\text{p}_{2}(\phi_{1\mathrel{{\cap}}2},\phi_{2\mathrel{{\cap}}3}^{*},\psi_{2}^{*},Y_{2})\text{p}_{2}(\phi_{1\mathrel{{\cap}}2},\phi_{2\mathrel{{\cap}}3})}{\text{p}_{2}(\phi_{1\mathrel{{\cap}}2},\phi_{2\mathrel{{\cap}}3},\psi_{2},Y_{2})\text{p}_{2}(\phi_{1\mathrel{{\cap}}2},\phi_{2\mathrel{{\cap}}3}^{*})}\frac{\text{q}(\phi_{2\mathrel{{\cap}}3},\psi_{2}\mid\phi_{2\mathrel{{\cap}}3}^{*},\psi_{2}^{*})}{\text{q}(\phi_{2\mathrel{{\cap}}3}^{*},\psi_{2}^{*}\mid\phi_{2\mathrel{{\cap}}3},\psi_{2})}.$
(64)

Our judicious choice of proposal distribution has resulted in a cancellation
in Equation (61) which removes all terms related to $\text{p}_{1}$. Similarly,
all terms related to $\text{p}_{1}$ are constant – hence cancel – in Equation
(64). This eliminates any need to re-evaluate the first submodel.

####  B.3 Stage three

In stage three we target the full melded posterior

$\text{p}_{\text{meld},3}(\phi_{1\mathrel{{\cap}}2},\phi_{2\mathrel{{\cap}}3},\psi_{1},\psi_{2},\psi_{3}\mid
Y_{1},Y_{2},Y_{3})\propto\\\
\text{p}_{\text{meld}}(\phi_{1\mathrel{{\cap}}2},\phi_{2\mathrel{{\cap}}3},\psi_{1},\psi_{2},\psi_{3}\mid
Y_{1},Y_{2},Y_{3}).$ (65)

The target has now been broadened to include terms from the third submodel and
the entirety of the pooled prior. Again, we employ a Metropolis-within-Gibbs
sampler, with proposals drawn such that

$\displaystyle\phi_{1\mathrel{{\cap}}2}^{*},\phi_{2\mathrel{{\cap}}3}^{*},\psi_{1}^{*},\psi_{2}^{*}\mid\psi_{3}$
$\displaystyle\sim\text{p}_{\text{meld},2}(\phi_{1\mathrel{{\cap}}2}^{*},\phi_{2\mathrel{{\cap}}3}^{*},\psi_{1}^{*},\psi_{2}^{*}\mid
Y_{1},Y_{2})$ (66)
$\displaystyle\psi_{3}^{*}\mid\phi_{1\mathrel{{\cap}}2},\phi_{2\mathrel{{\cap}}3},\psi_{1},\psi_{2}$
$\displaystyle\sim\text{q}(\psi_{3}^{*}\mid\psi_{3}),$ (67)

which leads to acceptance probabilities of

$\displaystyle\alpha$
$\displaystyle\\!\\!\begin{multlined}\left((\phi_{1\mathrel{{\cap}}2}^{*},\phi_{2\mathrel{{\cap}}3}^{*},\psi_{1}^{*},\psi_{2}^{*}),(\phi_{1\mathrel{{\cap}}2},\phi_{2\mathrel{{\cap}}3},\psi_{1},\psi_{2})\right)=\\\
\frac{\text{p}_{\text{pool},3}(\phi_{1\mathrel{{\cap}}2}^{*},\phi_{2\mathrel{{\cap}}3}^{*})}{\text{p}_{\text{pool},3}(\phi_{1\mathrel{{\cap}}2},\phi_{2\mathrel{{\cap}}3})}\frac{\text{p}_{3}(\phi_{2\mathrel{{\cap}}3}^{*},\psi_{3},Y_{3})\text{p}_{3}(\phi_{2\mathrel{{\cap}}3})}{\text{p}_{3}(\phi_{2\mathrel{{\cap}}3},\psi_{3},Y_{3})\text{p}_{3}(\phi_{2\mathrel{{\cap}}3}^{*})}\end{multlined}\left((\phi_{1\mathrel{{\cap}}2}^{*},\phi_{2\mathrel{{\cap}}3}^{*},\psi_{1}^{*},\psi_{2}^{*}),(\phi_{1\mathrel{{\cap}}2},\phi_{2\mathrel{{\cap}}3},\psi_{1},\psi_{2})\right)=\\\
\frac{\text{p}_{\text{pool},3}(\phi_{1\mathrel{{\cap}}2}^{*},\phi_{2\mathrel{{\cap}}3}^{*})}{\text{p}_{\text{pool},3}(\phi_{1\mathrel{{\cap}}2},\phi_{2\mathrel{{\cap}}3})}\frac{\text{p}_{3}(\phi_{2\mathrel{{\cap}}3}^{*},\psi_{3},Y_{3})\text{p}_{3}(\phi_{2\mathrel{{\cap}}3})}{\text{p}_{3}(\phi_{2\mathrel{{\cap}}3},\psi_{3},Y_{3})\text{p}_{3}(\phi_{2\mathrel{{\cap}}3}^{*})}$
(70) $\displaystyle\alpha$
$\displaystyle\left(\psi_{3}^{*},\psi_{3}\right)=\frac{\text{p}_{3}(\phi_{2\mathrel{{\cap}}3},\psi_{3}^{*},Y_{3})\text{q}(\psi_{3}\mid\psi_{3}^{*})}{\text{p}_{3}(\phi_{2\mathrel{{\cap}}3},\psi_{3},Y_{3})\text{q}(\psi_{3}^{*}\mid\psi_{3})}.$
(71)

The informed choice of proposal distribution for
($\phi_{1\mathrel{{\cap}}2},\phi_{2\mathrel{{\cap}}3},\psi_{1},\psi_{2}$) has
allowed us to target the full melded posterior without needing to evaluate all
submodels simultaneously.

###  C Normal approximation calculations

Substituting in the approximations of Section 3.2 to Equation (36) yields the
approximate melded posterior

$\widehat{\text{p}}_{\text{meld}}(\phi_{1\mathrel{{\cap}}2},\phi_{2\mathrel{{\cap}}3},\psi_{2}\mid
Y_{1},Y_{2},Y_{3})\,\propto\,\\\
\frac{\widehat{\text{p}}_{1}(\phi_{1\mathrel{{\cap}}2}\mid\widehat{\mu}_{1},\widehat{\Sigma}_{1})}{\widehat{\text{p}}_{1}(\phi_{1\mathrel{{\cap}}2}\mid\widehat{\mu}_{1,0},\widehat{\Sigma}_{1,0})}\frac{\widehat{\text{p}}_{3}(\phi_{2\mathrel{{\cap}}3}\mid\widehat{\mu}_{3},\widehat{\Sigma}_{3})}{\widehat{\text{p}}_{3}(\phi_{2\mathrel{{\cap}}3}\mid\widehat{\mu}_{3,0},\widehat{\Sigma}_{3,0})}\text{p}_{2}(\phi_{1\mathrel{{\cap}}2},\phi_{2\mathrel{{\cap}}3},\psi_{2}\mid
Y_{2}).$ (72)

Noting that the product of independent normal densities is an unnormalised
multivariate normal density with independent components, we rewrite Equation
(72) as

$\begin{gathered}\widehat{\text{p}}_{\text{meld}}(\phi_{1\mathrel{{\cap}}2},\phi_{2\mathrel{{\cap}}3},\psi_{2}\mid
Y_{1},Y_{2},Y_{3})\propto\frac{\widehat{\text{p}}_{\text{nu}}(\phi_{1\mathrel{{\cap}}2},\phi_{2\mathrel{{\cap}}3}\mid\widehat{\mu}_{\text{nu}},\widehat{\Sigma}_{\text{nu}})}{\widehat{\text{p}}_{\text{de}}(\phi_{1\mathrel{{\cap}}2},\phi_{2\mathrel{{\cap}}3}\mid\widehat{\mu}_{\text{de}},\widehat{\Sigma}_{\text{de}})}\text{p}_{2}(\phi_{1\mathrel{{\cap}}2},\phi_{2\mathrel{{\cap}}3},\psi_{2}\mid
Y_{2}),\\\\[6.45831pt]
\widehat{\mu}_{\text{nu}}=\begin{bmatrix}\widehat{\mu}_{1}\\\
\widehat{\mu}_{3}\end{bmatrix}\\!\\!,\quad\widehat{\Sigma}_{\text{nu}}=\begin{bmatrix}\widehat{\Sigma}_{1}&0\\\
0&\widehat{\Sigma}_{3}\end{bmatrix}\\!\\!,\quad\widehat{\mu}_{\text{de}}=\begin{bmatrix}\widehat{\mu}_{1,0}\\\
\widehat{\mu}_{3,0}\end{bmatrix}\\!\\!,\quad\widehat{\Sigma}_{\text{de}}=\begin{bmatrix}\widehat{\Sigma}_{1,0}&0\\\
0&\widehat{\Sigma}_{3,0}\end{bmatrix}\\!\\!.\end{gathered}$ (73)

The ratio of normal densities is also an unnormalised normal density, and
hence Equation (73) simplifies to

$\begin{gathered}\widehat{\text{p}}_{\text{meld}}(\phi_{1\mathrel{{\cap}}2},\phi_{2\mathrel{{\cap}}3},\psi_{2}\mid
Y_{1},Y_{2},Y_{3})\,\propto\,\widehat{\text{p}}(\phi_{1\mathrel{{\cap}}2},\phi_{2\mathrel{{\cap}}3}\mid\widehat{\mu},\,\widehat{\Sigma})\text{p}_{2}(\phi_{1\mathrel{{\cap}}2},\phi_{2\mathrel{{\cap}}3},\psi_{2}\mid
Y_{2}),\\\\[6.45831pt]
\widehat{\Sigma}=\left(\widehat{\Sigma}_{\text{nu}}^{-1}-\widehat{\Sigma}_{\text{de}}^{-1}\right)^{-1},\quad\widehat{\mu}=\widehat{\Sigma}\left(\widehat{\Sigma}_{\text{nu}}^{-1}\widehat{\mu}_{\text{nu}}-\widehat{\Sigma}_{\text{de}}^{-1}\widehat{\mu}_{\text{de}}\right).\end{gathered}$
(74)

###  D Calculating the cumulative fluid balance from the raw fluid data

In the raw fluid data each patient has $\tilde{l}=1,\ldots,\tilde{L}_{i}$
observations. Each observation $\tilde{x}_{i,\tilde{l}}$ is typically a small
fluid administration (e.g. an injection of some medicine in saline solution),
or a fluid discharge (almost always urine excretion). The observations have
corresponding observation times $\tilde{u}_{i,\tilde{l}}$, with
$\tilde{\boldsymbol{u}}_{i}=\\{\tilde{u}_{i,1},\ldots,\tilde{u}_{i,\tilde{L}_{i}}\\}$
and
$\tilde{\boldsymbol{x}}_{i}=\\{\tilde{x}_{i,1},\ldots,\tilde{x}_{i,\tilde{L}_{i}}\\}$.
We code the fluid administrations/inputs as positive values, and the
excretions/outputs as negative values. Each patient has an enormous number of
raw fluid observations $(L_{i}\ll\tilde{L}_{i})$ and it is computationally
infeasible to consider them all at once. We aggregate the raw fluid
observations into 8-hourly changes in fluid balance. From these 8-hourly
changes we calculate the cumulative fluid balance.

Mathematically, we define an ordered vector of boundary values

$\boldsymbol{v}_{i}=(\lfloor\min\\{\tilde{\boldsymbol{u}}_{i}\\}\rfloor,\lfloor\min\\{\tilde{\boldsymbol{u}}_{i}\\}\rfloor+\frac{1}{3},\ldots,\lceil\max\\{\tilde{\boldsymbol{u}}_{i}\\}\rceil),$
(75)

noting that $\dim(\boldsymbol{v}_{i})=L_{i}+1$. Because the observation times
encoded as _days since ICU admission_ and we are interested in the 8-hourly
changes, our floor and ceiling functions round down or up to the appropriate
third respectively. The raw fluid observations are then divided up into
$L_{i}$ subsets of
$\\{\tilde{\boldsymbol{x}}_{i},\tilde{\boldsymbol{u}}_{i}\\}$ based on which
boundary values the observation falls in between:

$V_{i,l}=\left\\{\\{\tilde{\boldsymbol{x}}_{i},\tilde{\boldsymbol{u}}_{i}\\}\mid
v_{i,l}\leq\tilde{\boldsymbol{u}}_{i}<v_{i,l+1}\right\\},$ (76)

for $l=1,\ldots,L_{i}$. Denote $N^{V}_{i,l}=\lvert V_{i,l}\rvert\mathop{/}2$
(dividing by two as $V_{i,l}$ contains both the observation and the
observation time). The $l$th 8-hourly fluid change $\Delta_{i,l}$ and
corresponding observation time $u_{i,l}$ can then be computed as

$\Delta_{i,l}=\sum_{s=1}^{N^{V}_{i,l}}\tilde{x}_{i,s},\,\,\text{s.t.}\,\,\tilde{x}_{i,s}\in
V_{i,l},\qquad
u_{i,l}=\frac{1}{N^{V}_{i,l}}\sum_{s=1}^{N^{V}_{i,l}}\tilde{u}_{i,s},\,\,\text{s.t.}\,\,\tilde{u}_{i,s}\in
V_{i,l}.$ (77)

Finally, the 8-hourly cumulative fluid balance data are computed by
$x_{i,l}=\sum_{s=1}^{l}\Delta_{i,s}$, and we assume they too are observed at
$u_{i,l}$.

###  E Priors and justification for the cumulative fluid submodel

The parameters for the gamma prior for $\eta^{b}_{1,i}$ and $\eta^{a}_{1,i}$
are obtained by assuming that the 2.5-, 50-, and 97.5- percentiles are at 0.5,
5, and 20 (Belgorodski _et al._ , 2017). A slope of $0.5$ (i.e. the change in
cumulative fluid balance per day) is unlikely but possible due to missing
data; a slope of $20$ is also unlikely but possible as extremely unwell
patients can have very high respiratory rates and thus require large fluid
inputs.

The prior for the breakpoint $\kappa_{i}$ is derived as follows. Define
$u_{i,(1)}=\min(\boldsymbol{u}_{i})$ and $u_{i,(n)}=\max(\boldsymbol{u}_{i})$,
with $r_{i}=u_{i,(n)}-u_{i,(1)}$. We reparameterise the breakpoint by noting
that $\kappa_{i}=\kappa^{\text{raw}}_{i}r_{i}+u_{i,(1)}$, where
$\kappa^{\text{raw}}\in[0,1]$. We then set
$\kappa^{\text{raw}}_{i}\sim\text{Beta}(5,5)$ to regularise the breakpoint
towards the middle of each individual’s stay in ICU. This is crucial to ensure
the submodel is identifiable when there is little evidence of a breakpoint in
the data. Note that this results in the following analytic expression for
$\text{p}_{2}(\phi_{2\mathrel{{\cap}}3})$

$\text{p}_{3}(\phi_{2\mathrel{{\cap}}3})=\prod_{i=1}^{N}\text{p}(\eta^{b}_{1,i})\text{p}(\eta^{a}_{1,i})\text{p}(\kappa_{i}),\,\,\text{with}\,\,\,\text{p}(\kappa_{i})=\text{p}_{\kappa^{\text{raw}}_{i}}\left(\frac{\kappa_{i}-u_{i,(1)}}{r_{i}}\right)\frac{1}{r_{i}}$
(78)

by the change of variables formula.

Specifying a prior for $\eta_{0,i}$, the cumulative fluid balance at
$\kappa_{i}$, is difficult because it too depends on the length of stay.
Instead, we reparameterise so that $\eta_{0,i}$ is a function of the
y-intercept $\eta_{0,i}^{\text{raw}}$.

$\eta_{0,i}=(\eta_{0,i}^{\text{raw}}+\eta^{b}_{1,i}\kappa_{i})\boldsymbol{1}_{\\{0<\kappa_{i}\\}}+(\eta_{0,i}^{\text{raw}}+\eta^{a}_{1,i}\kappa_{i})\boldsymbol{1}_{\\{0\geq\kappa_{i}\\}}$
(79)

We place a $\text{LogNormal}(1.61,0.47^{2})$ prior on
$\eta_{0,i}^{\text{raw}}$. These values are obtained assuming that, a priori,
the $2.5\%,50\%$, and $99\%$ percentiles of $\eta_{0,i}^{\text{raw}}$ are
$0.5,5$, and $15$ respectively (Belgorodski _et al._ , 2017). This is a broad
prior that reflects the numerous possible admission routes into the ICU. We
expect those admitted from the wards to have little pre-admission fluid data.
Those admitted from the operating theatre occasionally have their in-theatre
fluid input recorded after admission into the ICU, with no easy way to
distinguish these records in the data.

###  F Analytic form for the survival function

The hazard at arbitrary time $t$ is

$\displaystyle h_{i}(t)=\gamma
t^{\gamma-1}\exp\left\\{\boldsymbol{w}_{i}\theta+\alpha\frac{\partial}{\partial
t}m_{i}(t)\right\\}$ $\displaystyle
m_{i}(t)=\eta_{0,i}+\eta^{b}_{1,i}(t-\kappa_{i})\boldsymbol{1}_{\\{t<\kappa_{i}\\}}+\eta^{a}_{1,i}(t-\kappa_{i})\boldsymbol{1}_{\\{t\geq\kappa_{i}\\}}$
$\displaystyle\frac{\partial}{\partial
t}m_{i}(t)=\eta^{b}_{1,i}\boldsymbol{1}_{\\{t<\kappa_{i}\\}}+\eta^{a}_{1,i}\boldsymbol{1}_{\\{t\geq\kappa_{i}\\}}.$

Then, for $t>\kappa_{i}$, the cumulative hazard is

$\displaystyle\int_{0}^{t}h_{i}(u)\text{d}u$ $\displaystyle=\int_{0}^{t}\gamma
u^{\gamma-1}\exp\left\\{\boldsymbol{w}_{i}\theta+\alpha\eta^{b}_{1,i}\boldsymbol{1}_{\\{u<\kappa_{i}\\}}+\alpha\eta^{a}_{1,i}\boldsymbol{1}_{\\{u\geq\kappa_{i}\\}}\right\\}\text{d}u$
$\displaystyle=\gamma\exp\\{\boldsymbol{w}_{i}\theta\\}\int_{0}^{t}u^{\gamma-1}\exp\left\\{\alpha\eta^{b}_{1,i}\boldsymbol{1}_{\\{u<\kappa_{i}\\}}+\alpha\eta^{a}_{1,i}\boldsymbol{1}_{\\{u\geq\kappa_{i}\\}}\right\\}\text{d}u$
$\displaystyle=\gamma\exp\\{\boldsymbol{w}_{i}\theta\\}\left[\int_{0}^{\kappa_{i}}u^{\gamma-1}\exp\left\\{\alpha\eta^{b}_{1,i}\right\\}\text{d}u+\int_{\kappa_{i}}^{t}u^{\gamma-1}\exp\left\\{\alpha\eta^{a}_{1,i}\right\\}\text{d}u\right]$
$\displaystyle=\exp\\{\boldsymbol{w}_{i}\theta\\}\left[\exp\left\\{\alpha\eta^{b}_{1,i}\right\\}\kappa_{i}^{\gamma}+\exp\left\\{\alpha\eta^{a}_{1,i}\right\\}(t^{\gamma}-\kappa_{i}^{\gamma})\right]$

and for $t<\kappa_{i}$

$\displaystyle\int_{0}^{t}h_{i}(u)\text{d}u$
$\displaystyle=\gamma\exp\\{\boldsymbol{w}_{i}\theta\\}\left[\int_{0}^{t}u^{\gamma-1}\exp\left\\{\alpha\eta^{b}_{1,i}\right\\}\text{d}u\right]$
$\displaystyle=\exp\\{\boldsymbol{w}_{i}\theta\\}\left[\exp\left\\{\alpha\eta^{b}_{1,i}\right\\}t_{i}^{\gamma}\right]$
$\displaystyle=t_{i}^{\gamma}\exp\\{\boldsymbol{w}_{i}\theta+\alpha\eta^{b}_{1,i}\\}.$

The survival functions are defined with appropriately for $t>\kappa_{i}$ and
$t<\kappa_{i}$ as $S_{i}(t)=\exp\\{-\int_{0}^{t}h_{i}(u)\text{d}u\\}$.

###  G Survival submodel prior justification

Our prior for $(\gamma,\alpha,\boldsymbol{\theta})$ must result in a plausible
distribution for $\text{p}_{2,i}(T_{i}\mid d_{i}=1)$, and a reasonable balance
between $d_{i}=1$ and $d_{i}=0$ events. The primary concern is unintentionally
specifying a prior for which the bulk of $\text{p}_{2,i}(T_{i}\mid d_{i}=1)$
is very close to zero. In addition, certain extreme configurations of
$(\gamma,\alpha,\boldsymbol{\theta})$ cause issues for the methodology of
Crowther and Lambert (2013), particularly the numerical root finding and
numerical integration steps. We would like to rule out such extreme
configurations a priori. Ideally we would encode this information a joint
prior for $(\gamma,\alpha,\boldsymbol{\theta})$, but specifying the
appropriate correlation structure for these parameters is prohibitively
challenging. Instead we focus on specifying appropriate marginals for each of
$\gamma,\alpha$, and $\boldsymbol{\theta}$, and create visual prior predictive
checks (Gabry _et al._ , 2019; Gelman _et al._ , 2020) to ensure the induced
prior for $(T_{i},d_{i})$ is acceptable.

Before justifying our chosen marginal prior, we note that the
$\exp\\{\boldsymbol{x}_{i}\theta+\alpha\frac{\partial}{\partial
T_{i}}m_{i}(T_{i})\\}$ term implies that the priors for $\theta$ and $\alpha$
are on the log-scale. Hence the magnitude of these parameters must be small,
otherwise all event times would be very near zero or at infinity. The
nonsymmetric effect of the transformation from the log scale also implies that
symmetric priors are not obviously sensible. From these observations we deduce
that $\theta$ and $\alpha$ must not be too large in magnitude, however if they
are negative then they can be slightly larger. Hence, we specify the skew-
normal priors detailed in Section 5.3, noting that the skewness parameter for
$\alpha$ is smaller, because $\frac{\partial}{\partial T_{i}}m_{i}(T_{i})$ is
strictly positive and typically between 0.5 and 20, whilst
$\boldsymbol{w}_{i}$ is standardised to be approximately standard normal.
Lastly, if $\gamma$ is too far away from $1$ (in either direction), then the
event times are very small either because the hazard increases rapidly
($\gamma\gg 1$), or because almost all of the cumulative hazard is in the
neighbourhood of 0 ($\gamma\ll 1$). We specify a gamma distribution for
$\gamma$ with the $1$th-, $50$th-, and $99$th-percentiles of
$\text{p}_{2}(\gamma)$ at $0.2,1$, and $2$, allowing for a wide range of
hazard shapes, but removing many of the extremes.

###  H Estimating submodel prior marginal distributions

For $\text{p}_{1}$, we note that
$\text{p}_{1}(\phi_{1\mathrel{{\cap}}2})=\prod_{i=1}^{N}\text{p}_{1,i}(T_{i},d_{i})$,
and that $\text{p}_{1,i}(T_{i},d_{i})$ conditions on each individual’s length
of stay (in specifying the location of the knots), as well as the range, mean,
and standard deviation of the P/F data (by standardising $\tilde{z}_{i,j}$).
Simple Monte Carlo samples are drawn from
$\text{p}_{1}(\phi_{1\mathrel{{\cap}}2})$ and used to estimate
$\widehat{\text{p}}_{1}(\phi_{1\mathrel{{\cap}}2})$. Under the second submodel
we obtain samples from
$\text{p}_{2}(\phi_{1\mathrel{{\cap}}2},\phi_{2\mathrel{{\cap}}3})$ using the
methodology of Crowther and Lambert (2013) as implemented in simsurv
(Brilleman, 2021). These samples are use to estimate
$\widehat{\text{p}}_{2}(\phi_{1\mathrel{{\cap}}2},\phi_{2\mathrel{{\cap}}3})$.

####  H.1 P/F submodel

We approximate $\text{p}_{1}(\phi_{1\mathrel{{\cap}}2})$ using a mixture of
discrete and continuous distributions, with a discrete spike at $C_{i}$ for
the censored events and a beta distribution for the (rescaled) event times.
Monte Carlo samples of $T_{i}$ and $d_{i}$ are obtained from
$\text{p}_{1,i}(T_{i},d_{i})$ by drawing $\beta_{0,i}$ and
$\boldsymbol{\zeta}_{i}$ from their respective prior distributions and then
solving (46). Denoting the estimated mixture weight
$\widehat{\pi}_{i}\in[0,1]$, the density estimate is

$\widehat{\text{p}}_{1,i}(T_{i},d_{i})=\widehat{\pi}_{i}\,\text{Beta}\left(\frac{T_{i}}{C_{i}};\widehat{a},\widehat{b}\right)\frac{1}{C_{i}}\boldsymbol{1}_{\\{d_{i}=1\\}}+(1-\widehat{\pi}_{i})\boldsymbol{1}_{\\{d_{i}=0,T_{i}=C_{i}\\}}$
(80)

where $\widehat{\pi}_{i},\widehat{a}_{i}$ and $\widehat{b}_{i}$ are maximum
likelihood estimates obtained using the prior samples. Examples of
$\widehat{\text{p}}_{1,i}(T_{i},d_{i})$ for a subset of individuals are
displayed in Figure 12.

Figure 12: Fitted distribution (curve) and Monte Carlo samples drawn from
$\text{p}_{1}(\phi_{1\mathrel{{\cap}}2})$ (histogram) for a subset of the
individuals in the cohort. The height of the atom at $C_{i}$ (red bar and
point) has been set to $1-\widehat{\pi}_{i}$

####  H.2 Survival submodel

Our estimate of
$\text{p}_{2}(\phi_{1\mathrel{{\cap}}2},\phi_{2\mathrel{{\cap}}3})$ relies on
the fact that
$\text{p}_{2}(\phi_{1\mathrel{{\cap}}2},\phi_{2\mathrel{{\cap}}3})=\prod_{i=1}^{N}\text{p}_{2,i}(T_{i},d_{i},\kappa_{i},\eta^{b}_{1,i},\eta^{a}_{1,i})$.
As such we estimate
$\text{p}_{2,i}(T_{i},d_{i},\kappa_{i},\eta^{b}_{1,i},\eta^{a}_{1,i})$ for
each individual and take the product of these estimates. Drawing samples from
$\text{p}_{2,i}(T_{i},d_{i},\kappa_{i},\eta^{b}_{1,i},\eta^{a}_{1,i})$ is
challenging: we use the approach proposed by Crowther and Lambert (2013) as
implemented in Brilleman (2021). Inspecting the samples reveals correlation
between $(T_{i},d_{i})$ and $(\kappa_{i},\eta^{b}_{1,i},\eta^{a}_{1,i})$ that
we would like to capture in our estimate. To do so, we fit a mixture of
multivariate normal distributions to transformations of the continuous
parameters with support on $\mathbb{R}$,

$\tilde{T}_{i}=\text{Logit}\left(\frac{T_{i}}{C_{i}}\right),\quad\tilde{\kappa}_{i}=\text{Logit}\left(\frac{\kappa_{i}-u_{i,(1)}}{u_{i,(n)}-u_{i,(1)}}\right),\quad\tilde{\eta}^{b}_{1,i}=\log(\eta^{b}_{1,i}),\quad\tilde{\eta}^{a}_{1,i}=\log(\eta^{a}_{1,i}).$
(81)

The resulting density estimate, with estimated mixture weight
$\widehat{\theta}_{i}\in[0,1]$, is

$\displaystyle\widehat{\text{p}}_{2}(T_{i},d_{i},\kappa_{i},\eta^{b}_{1,i},\eta^{a}_{1,i})=\Big{[}$
$\displaystyle(\widehat{\theta}_{i})\text{N}\left(\left[\tilde{T}_{i},\tilde{\kappa}_{i},\tilde{\eta}^{b}_{1,i},\tilde{\eta}^{a}_{1,i}\right]^{\mathsf{T}};\widehat{\mu}_{1,i},\widehat{\Sigma}_{1,i}\right)\boldsymbol{1}_{\\{d_{i}=1\\}}+$
$\displaystyle(1-\widehat{\theta}_{i})\text{N}\left(\left[\tilde{\kappa}_{i},\tilde{\eta}^{b}_{1,i},\tilde{\eta}^{a}_{1,i}\right]^{\mathsf{T}};\widehat{\mu}_{0,i},\widehat{\Sigma}_{0,i}\right)\boldsymbol{1}_{\\{d_{i}=0,T_{i}=C_{i}\\}}\Big{]}J_{i},$

where
$\widehat{\theta}_{i},\widehat{\mu}_{1,i},\widehat{\Sigma}_{1,i},\widehat{\mu}_{0,i}$,
and $\widehat{\Sigma}_{0,i}$ are maximum likelihood estimates, and the
Jacobian correction $J_{i}$ is

$J_{i}=\left[\left(\frac{1}{C_{i}-T_{i}}+\frac{1}{T_{i}}\right)^{d_{i}}\left(\frac{1}{u_{i,(n)}-\kappa_{i}}+\frac{1}{\kappa_{i}-u_{i,(1)}}\right)\left(\frac{1}{\eta^{b}_{1,i}}\right)\left(\frac{1}{\eta^{a}_{1,i}}\right)\right].$
(82)

We assess the fit of this estimate by drawing samples from
$\widehat{\text{p}}_{2,i}(\phi_{1\mathrel{{\cap}}2},\phi_{2\mathrel{{\cap}}3})$
and comparing them to the Monte Carlo samples drawn using simsurv. Our visual
assessment is displayed in Figure 14 for individual $i=19$. The normal
approximation seems to fit the samples well, with the shape of
$\text{p}_{2,i}(T_{i}\mid d_{i}=1)$ closely matching that of the Monte Carlo
samples, and with a similar mix of $d_{i}=0$ and $d_{i}=1$ samples.

We also require an estimate of $\text{p}_{2}(\phi_{1\mathrel{{\cap}}2})$ for
experiments discussed in Section 4.7. This is obtained using the samples
generated under the survival submodel prior and the methodology of Section
H.1. The raw samples and fit are displayed in Figure 13.

Figure 13: Fitted distribution (curve) and Monte Carlo samples drawn from
$\text{p}_{2}(\phi_{1\mathrel{{\cap}}2})$ (histogram) for a subset of the
individuals in the cohort. The height of the atom at $C_{i}$ (red bar and
point) has been set to $1-\widehat{\pi}_{i}$ Figure 14: Monte Carlo (MC)
samples from
$\text{p}_{2,i}(\phi_{1\mathrel{{\cap}}2},\phi_{2\mathrel{{\cap}}3})$ obtained
using simsurv and samples from the fitted normal approximation (NA) for
$i=19$. The panels on the off diagonal elements contain a 2D kernel density
estimate for $d_{i}=1$ and the samples for $d_{i}=0$. Diagonal and lower-
triangular panels are on their original scales, whilst the upper-triangular
panels are on the log scale.

###  I Cohort selection criteria

This appendix details the cohort selection criteria and our rationale for
them. In the text we speak of the $i$th individual. This is because in our
final data set (the data that results from applying the following criteria) we
are dealing with unique individuals, however some individuals in MIMIC have
multiple ICU stays. In this appendix $i$ represents a single stay in ICU.

  1. 1.

Each ICU stay must have at least 12 PaO2/FiO2 observations ($J_{i}\geq 12$),
with the first 6 being greater than 350 ($z_{i,j}>350$ for $j=1,\dots,1$).

     * •

This is to ensure we have enough data to fit a B-spline with 7 internal knots.
The restriction on the first 6 observations is to avoid selecting those who
have already started to experience respiratory failure prior to ICU admission.

  2. 2.

The time between any 2 consecutive PaO2/FiO2 observations cannot exceed 2
days.

     * •

This is because we believe this lack of PaO2/FiO2 observations is likely to be
an error in the data: e.g. what appears as a single ICU stay in the data is
actually two or more separate stays.

  3. 3.

The fluid observations must be after ICU admission (some observations are
entered as ‘Pre-admission intake’) and cannot be associated with fluid
administered in the operating room (OR). Note that this does not mean all OR
fluid administrations are removed, as some are mis-labelled.

  4. 4.

There must be sufficient temporal overlap between the fluid data and the
PaO2/FiO2 data. Specifically,

$\frac{\max\left\\{0,\min\left[\max(\boldsymbol{t}_{i}),\max(\tilde{\boldsymbol{u}}_{i})\right]-\max\left[\min(\boldsymbol{t}_{i}),\min(\tilde{\boldsymbol{u}}_{i})\right]\right\\}}{\max\left[\max(\boldsymbol{t}_{i}),\max(\tilde{\boldsymbol{u}}_{i})\right]-\min\left[\min(\boldsymbol{t}_{i}),\min(\tilde{\boldsymbol{u}}_{i})\right]}>0.9$
(83)

     * •

The numerator of (83) is strictly positive, and the denominator ensures that
the quantity is bounded between 0 and 1.

     * •

We cannot investigate the relationship between the rate of fluid intake and
respiratory failure if the latter occurs without sufficient fluid data
surrounding the event.

###  J Baseline covariate information

The baseline covariate vector $\boldsymbol{w}_{i}$ contains the median of the
measurements taken in the first 24 hours of the ICU stay, which are then
stardardised, of the following covariates: Anion gap, Bicarbonate, Creatinine,
Chloride, Glucose, Hematocrit, Hemoglobin, Platelet, Partial thromboplastin
time, International normalized ratio, Prothrombin time, Sodium, blood Urea
nitrogen, White blood cell count, Age at ICU admission, and Gender.

###  K Updating $\phi_{1\mathrel{{\cap}}2}$ and $\phi_{2\mathrel{{\cap}}3}$ in
stage two individual-at-a-time

The parallel sampler described in Section 3.1 is a MH-within-Gibbs scheme,
with each iteration sampling from the conditionals of the melded model

$\begin{gathered}\text{p}_{\text{meld}}(\phi_{1\mathrel{{\cap}}2},\psi_{1}\mid\boldsymbol{Y},\psi_{2},\phi_{2\mathrel{{\cap}}3},\psi_{3}),\\\
\text{p}_{\text{meld}}(\phi_{2\mathrel{{\cap}}3},\psi_{3}\mid\boldsymbol{Y},\phi_{1\mathrel{{\cap}}2},\psi_{1},\psi_{2}),\\\
\text{p}_{\text{meld}}(\psi_{2}\mid\boldsymbol{Y},\phi_{1\mathrel{{\cap}}2},\psi_{1},\phi_{2\mathrel{{\cap}}3},\psi_{3}).\end{gathered}$
(84)

We would like to update the first two conditionals ‘individual-at-a-time’. For
simplicity we will assume the pooled prior is formed using product-of-experts
pooling, and focus on the first conditional
$\text{p}_{\text{meld}}(\phi_{1\mathrel{{\cap}}2},\psi_{1}\mid\boldsymbol{Y},\psi_{2},\phi_{2\mathrel{{\cap}}3},\psi_{3})$.
Similar arguments apply to the second conditional and other pooling types.
Recall that

$\displaystyle\phi_{1\mathrel{{\cap}}2}$
$\displaystyle=(\phi_{1\mathrel{{\cap}}2,1},\ldots,\phi_{1\mathrel{{\cap}}2,N})$
(85)
$\displaystyle=(\phi_{1\mathrel{{\cap}}2,1}(\chi_{1,1}),\ldots,\phi_{1\mathrel{{\cap}}2,N}(\chi_{1,N}))$
$\displaystyle=\left((T_{1},d_{1}),\ldots,(T_{N},d_{N})\right),$

and $\psi_{1}=(\psi_{1,1},\ldots,\psi_{1,N})=(\omega_{i})_{i=1}^{N}$. We would
like to update $\phi_{1\mathrel{{\cap}}2}$ and $\psi_{1}$ by updating
$(T_{i},d_{i},\omega_{i})$ for one individual at a time, for a total of $N$
‘sub-steps’.

####  K.1 Sub-step 1

Suppose we are at step $t-1$ of the Markov chain, and we are proposing values
for step $t$. Without any loss of generality we assume that we are updating
individual $i=1$ in sub-step 1 – in practice we update the individuals in a
random order for each iteration of the MH-within-Gibbs scheme.

Our target is

$\text{p}_{1}(\chi_{1,1},\psi_{1,1}\mid
Y_{1},(\chi_{1,i},\psi_{1,i})_{i=2}^{N})\text{p}_{2}(\phi_{1\mathrel{{\cap}}2,1}(\chi_{1,1})\mid
Y_{2},(\phi_{1\mathrel{{\cap}}2,i}(\chi_{1,i}))_{i=2}^{N},\psi_{2},\phi_{2\mathrel{{\cap}}3}).$
(86)

The model, detailed in Section 5.1, factorises such that

$\text{p}_{1}(\chi_{1},\psi_{1}\mid
Y_{1})=\prod_{i=1}^{N}\text{p}_{1}(\chi_{1,i},\psi_{1,i}\mid Y_{1}),$ (87)

which implies131313This property is also true of $\text{p}_{3}$.
$(\chi_{1,i},\psi_{1,i}\perp\\!\\!\\!\perp\chi_{1,i^{{}^{\prime}}},\psi_{1,i^{{}^{\prime}}})\mid
Y_{1}$ for $i\neq i^{{}^{\prime}}$, hence

$\text{p}_{1}(\chi_{1,1},\psi_{1,1}\mid
Y_{1},(\chi_{1,i},\psi_{1,i})_{i=2}^{N})=\text{p}_{1}(\chi_{1,1},\psi_{1,1}\mid
Y_{1}).$ (88)

It will be convenient to rewrite (86) as

$\text{p}_{1}(\chi_{1,1},\psi_{1,1}\mid
Y_{1})\frac{\text{p}_{2}(\phi_{1\mathrel{{\cap}}2,1}(\chi_{1,1}),(\phi_{1\mathrel{{\cap}}2,i}(\chi_{1,i}))_{i=2}^{N},\psi_{2},\phi_{2\mathrel{{\cap}}3}\mid
Y_{2})}{\text{p}_{2}((\phi_{1\mathrel{{\cap}}2,i}(\chi_{1,i}))_{i=2}^{N},\psi_{2},\phi_{2\mathrel{{\cap}}3}\mid
Y_{2})}.$ (89)

Suppose there are $K_{1}$ stage one samples from
$\text{p}_{1}(\chi_{1},\psi_{1}\mid Y_{1})$, and each sample has a
corresponding $\phi_{1\mathrel{{\cap}}2}$. We propose $\psi_{1,1}^{*}$ and
$\chi_{1,1}^{*}$ (hence $\phi_{1\mathrel{{\cap}}2,1}^{*})$ by sampling a
random integer $k_{1}^{*}$ from $\\{1,\ldots,K_{1}\\}$, retrieving the
corresponding values of $\psi_{1}$ and $\chi_{1}$, and ignoring
$(\chi_{1,2},\ldots,\chi_{1,N},\psi_{1,2},\ldots,\psi_{1,N})$. This is
equivalent to proposing from $\text{p}_{1}(\chi_{1,1},\psi_{1,1}\mid Y_{1})$
Under this proposal, noting that the denominator term in (89) does not depend
on the proposed parameters, the acceptance probability for this sub-step is

$\displaystyle\alpha$
$\displaystyle\left((\chi_{1,1}^{*},\psi_{1}^{*})_{k_{1}^{*}},(\chi_{1,1},\psi_{1})\right)$
(90) $\displaystyle=\frac{\text{p}_{1}(\chi_{1,1}^{*},\psi_{1,1}^{*}\mid
Y_{1})\text{p}_{2}(\phi_{1\mathrel{{\cap}}2,1}(\chi_{1,1}^{*}),(\phi_{1\mathrel{{\cap}}2,i}(\chi_{1,i}))_{i=2}^{N},\phi_{2\mathrel{{\cap}}3},\psi_{2}\mid
Y_{2})}{\text{p}_{1}(\chi_{1,1},\psi_{1,1}\mid
Y_{1})\text{p}_{2}(\phi_{1\mathrel{{\cap}}2,1}(\chi_{1,1}),(\phi_{1\mathrel{{\cap}}2,i}(\chi_{1,i}))_{i=2}^{N},\phi_{2\mathrel{{\cap}}3},\psi_{2}\mid
Y_{2})}\frac{\text{p}_{1}(\chi_{1,1},\psi_{1}\mid
Y_{1})}{\text{p}_{1}(\chi_{1,1}^{*},\psi_{1}^{*}\mid Y_{1})}$
$\displaystyle=\frac{\text{p}_{2}(\phi_{1\mathrel{{\cap}}2,1}(\chi_{1,1}^{*}),(\phi_{1\mathrel{{\cap}}2,i}(\chi_{1,i}))_{i=2}^{N},\phi_{2\mathrel{{\cap}}3},\psi_{2}\mid
Y_{2})}{\text{p}_{2}(\phi_{1\mathrel{{\cap}}2,1}(\chi_{1,1}),(\phi_{1\mathrel{{\cap}}2,i}(\chi_{1,i}))_{i=2}^{N},\phi_{2\mathrel{{\cap}}3},\psi_{2}\mid
Y_{2})},$

which does not depend on $\text{p}_{1}$. We store the value of $k_{1}^{*}$
associated with the proposal (if it is accepted) in order to resample the
stage one values of $\psi_{1,1}$ to reflect the information in the other
submodels.

####  K.2 Sub-step $n$

At sub-step $n$, for $1<n\leq N$ with $(\chi_{1,i})_{i=N+1}^{N}=\varnothing$,
we have updated $((\chi_{1,i},\psi_{1,i})_{i=1}^{n-1})$ to
$((\chi_{1,i}^{[t]},\psi_{1,i}^{[t]})_{i=1}^{n-1})$. Thus the target is

$\displaystyle\text{p}_{1}$ $\displaystyle\left(\chi_{1,n},\psi_{1,n}\mid
Y_{1},(\chi_{1,i}^{[t]},\psi_{1,i}^{[t]})_{i=1}^{n-1},(\chi_{1,i},\psi_{1,i})_{i=n+1}^{N}\right)$
(91)
$\displaystyle\times\text{p}_{2}\left(\phi_{1\mathrel{{\cap}}2,n}(\chi_{1,n})\mid
Y_{2},\psi_{2},\phi_{2\mathrel{{\cap}}3},(\phi_{1\mathrel{{\cap}}2,i}(\chi_{1,i}^{[t]}))_{i=1}^{n-1},(\phi_{1\mathrel{{\cap}}2,i}(\chi_{1,i}))_{i=n+1}^{N}\right)$
$\displaystyle=$ $\displaystyle\text{p}_{1}\left(\chi_{1,n},\psi_{1,n},\mid
Y_{1}\right)\frac{\text{p}_{2}\left(\phi_{1\mathrel{{\cap}}2,n}(\chi_{1,n}),\psi_{2},\phi_{2\mathrel{{\cap}}3},(\phi_{1\mathrel{{\cap}}2,i}(\chi_{1,i}^{[t]}))_{i=1}^{n-1},(\phi_{1\mathrel{{\cap}}2,i}(\chi_{1,i}))_{i=n+1}^{N}\mid
Y_{2}\right)}{\text{p}_{2}\left(\psi_{2},\phi_{2\mathrel{{\cap}}3},(\phi_{1\mathrel{{\cap}}2,i}(\chi_{1,i}^{[t]}))_{i=1}^{n-1},(\phi_{1\mathrel{{\cap}}2,i}(\chi_{1,i}))_{i=n+1}^{N}\mid
Y_{2}\right)}$ (92)

We propose $\chi_{1,n}^{*},\psi_{1,n}^{*}$ in the same manner as the previous
section. The corresponding acceptance probability is

$\displaystyle\alpha((\chi_{1,n}^{*},\psi_{1,n}^{*}),$
$\displaystyle(\chi_{1,n},\psi_{1,n}))$ (93) $\displaystyle=$
$\displaystyle\frac{\text{p}_{1}\left(\chi_{1,n}^{*},\psi_{1,n}^{*}\mid
Y_{1}\right)}{\text{p}_{1}\left(\chi_{1,n},\psi_{1,n}\mid Y_{1}\right)}$
$\displaystyle\times\frac{\text{p}_{2}(\phi_{1\mathrel{{\cap}}2,n}(\chi_{1,n}^{*}),\psi_{2},\phi_{2\mathrel{{\cap}}3},(\chi_{1,i}^{[t]})_{i=1}^{n-1},(\chi_{1,i})_{i=n+1}^{N}\mid
Y_{2})}{\text{p}_{2}(\phi_{1\mathrel{{\cap}}2,n}(\chi_{1,n}),\psi_{2},\phi_{2\mathrel{{\cap}}3},(\chi_{1,i}^{[t]})_{i=1}^{n-1},(\chi_{1,i})_{i=n+1}^{N}\mid
Y_{2})}$ $\displaystyle\times\frac{\text{p}_{1}(\chi_{1,n},\psi_{1,n}\mid
Y_{1})}{\text{p}_{1}(\chi_{1,n}^{*},\psi_{1,n}^{*}\mid Y_{1})}$
$\displaystyle=$
$\displaystyle\frac{\text{p}_{2}(\phi_{1\mathrel{{\cap}}2,n}(\chi_{1,n}^{*}),\psi_{2},\phi_{2\mathrel{{\cap}}3},(\chi_{1,i}^{[t]})_{i=1}^{n-1},(\chi_{1,i})_{i=n+1}^{N}\mid
Y_{2})}{\text{p}_{2}(\phi_{1\mathrel{{\cap}}2,n}(\chi_{1,n}),\psi_{2},\phi_{2\mathrel{{\cap}}3},(\chi_{1,i}^{[t]})_{i=1}^{n-1},(\chi_{1,i})_{i=n+1}^{N}\mid
Y_{2})}.$

## References

reAbadi, F., Gimenez, O., Ullrich, B., et al. (2010) Estimation of immigration
rate using integrated population models. _Journal of Applied Ecology_ , 47,
393–400. DOI: 10.1111/j.1365-2664.2010.01789.x.

preAdes, A. E. and Sutton, A. J. (2006) Multiparameter evidence synthesis in
epidemiology and medical decision-making: Current approaches. _Journal of the
Royal Statistical Society: Series A (Statistics in Society)_ , 169, 5–35. DOI:
10.1111/j.1467-985X.2005.00377.x.

preBelgorodski, N., Greiner, M., Tolksdorf, K., et al. (2017)
rriskDistributions: Fitting distributions to given data or known quantiles.

preBrilleman, S. (2021) Simsurv: Simulate survival data.

preBrilleman, S. L., Elci, E. M., Novik, J. B., et al. (2020) Bayesian
survival analysis using the rstanarm R package. _arXiv:2002.09633 [stat]_.
Available at: https://arxiv.org/abs/2002.09633.

preBromiley, P. (2003) Products and convolutions of Gaussian probability
density functions. _Tina-Vision Memo_ , 3, 1.

preBurke, D. L., Ensor, J. and Riley, R. D. (2017) Meta-analysis using
individual participant data: One-stage and two-stage approaches, and why they
may differ. _Statistics in Medicine_ , 36, 855–875. DOI: 10.1002/sim.7141.

preCarpenter, B., Gelman, A., Hoffman, M., et al. (2017) Stan: A probabilistic