Model formulation
A shared parameter joint model consists of related submodels which are specified separately for each of the longitudinal and time-to-event outcomes. These are therefore commonly referred to as the longitudinal submodel(s) and the event submodel. The longitudinal and event submodels are linked using shared individual-specific parameters, which can be parameterised in a number of ways. We describe each of these submodels below.
Longitudinal submodel(s)
We assume \(y_{ijm}(t) = y_{im}(t_{ij})\) corresponds to the observed value of the \(m^{th}\) \((m = 1,...,M)\) biomarker for individual \(i\) \((i = 1,...,N)\) taken at time point \(t_{ij}\), \(j = 1,...,n_{im}\). We specify a (multivariate) generalised linear mixed model that assumes \(y_{ijm}(t)\) follows a distribution in the exponential family with mean \(\mu_{ijm}(t)\) and linear predictor
\[
\eta_{ijm}(t) = g_m(\mu_{ijm}(t)) =
\boldsymbol{x}^T_{ijm}(t) \boldsymbol{\beta}_m +
\boldsymbol{z}^T_{ijm}(t) \boldsymbol{b}_{im}
\]
where \(\boldsymbol{x}^T_{ijm}(t)\) and \(\boldsymbol{z}^T_{ijm}(t)\) are both row-vectors of covariates (which likely include some function of time, for example a linear slope, cubic splines, or polynomial terms) with associated vectors of fixed and individual-specific parameters \(\boldsymbol{\beta}_m\) and \(\boldsymbol{b}_{im}\), respectively, and \(g_m\) is some known link function. The distribution and link function are allowed to differ over the \(M\) longitudinal submodels. We let the vector \(\boldsymbol{\beta} = \{ \boldsymbol{\beta}_m ; m = 1,...,M\}\) denote the collection of population-level parameters across the \(M\) longitudinal submodels. We assume that the dependence across the different longitudinal submodels (i.e. the correlation between the different longitudinal biomarkers) is captured through a shared multivariate normal distribution for the individual-specific parameters; that is, we assume
\[
\begin{pmatrix} \boldsymbol{b}_{i1} \\ \vdots \\ \boldsymbol{b}_{iM} \end{pmatrix} =
\boldsymbol{b}_i \sim
\mathsf{Normal} \left( 0 , \boldsymbol{\Sigma} \right)
\]
for some unstructured variance-covariance matrix \(\boldsymbol{\Sigma}\).
Event submodel
We assume that we also observe an event time \(T_i = \mathsf{min} \left( T^*_i , C_i \right)\) where \(T^*_i\) denotes the so-called “true” event time for individual \(i\) (potentially unobserved) and \(C_i\) denotes the censoring time. We define an event indicator \(d_i = I(T^*_i \leq C_i)\). We then model the hazard of the event using a parametric proportional hazards regression model of the form
\[
h_i(t) = h_0(t; \boldsymbol{\omega}) \mathsf{exp}
\left(
\boldsymbol{w}^T_i(t) \boldsymbol{\gamma} +
\sum_{m=1}^M \sum_{q=1}^{Q_m}
f_{mq}(\boldsymbol{\beta}, \boldsymbol{b}_{i}, \alpha_{mq}; t)
\right)
\]
where \(h_i(t)\) is the hazard of the event for individual \(i\) at time \(t\), \(h_0(t; \boldsymbol{\omega})\) is the baseline hazard at time \(t\) given parameters \(\boldsymbol{\omega}\), \(\boldsymbol{w}^T_i(t)\) is a row-vector of individual-specific covariates (possibly time-dependent) with an associated vector of regression coefficients \(\boldsymbol{\gamma}\) (log hazard ratios), \(f_{mq}(.)\) are a set of known functions for \(m=1,...,M\) and \(q=1,...,Q_m\), and the \(\alpha_{mq}\) are regression coefficients (log hazard ratios).
The longitudinal and event submodels are assumed to be related via an “association structure”, which is a set of functions each \(\{ f_{mq} ; m = 1,...,M, q = 1,...,Q_m \}\) that may each be conditional on the population-level parameters from the longitudinal submodel \(\boldsymbol{\beta}\), the individual-specific parameters \(\boldsymbol{b}_{i}\), and the population-level parameters \(\alpha_{mq}\) for \(m=1,...,M\) and \(q=1,...,Q_m\). That is, the association structure of the joint model is captured via the \(\sum_{m=1}^M \sum_{q=1}^{Q_m} f_{mq}(\boldsymbol{\beta}_m, \boldsymbol{b}_{im}, \alpha_{mq}; t)\) term in the linear predictor of the event submodel. The \(\alpha_{mq}\) are referred to as the “association parameters” since they quantify the strength of the association between the longitudinal and event processes. The various ways in which the association structure can be are described in the next section.
The probability of individual \(i\) still being event-free at time \(t\), often referred to as the “survival probability”, is defined as
\[
S_i(t) =
\text{Prob} \Big( T_i^* \geq t \Big) =
\exp \Big( -H_i(t) \Big)
\]
where \(H_i(t) = \int_{s=0}^t h_i(s) ds\) is the cumulative hazard for individual \(i\).
We assume that the baseline hazard \(h_0(t; \boldsymbol{\omega})\) is modelled parametrically. In the stan_jm modelling function the baseline hazard be specified as either: an approximation using B-splines on the log hazard scale (the default); a Weibull distribution; or an approximation using a piecewise constant function on the log hazard scale (sometimes referred to as piecewise exponential). The choice of baseline hazard can be made via the basehaz argument. In the case of the B-splines or piecewise constant baseline hazard, the user can control the flexibility by specifying the knots or degrees of freedom via the basehaz_ops argument. (Note that currently there is slightly limited post-estimation functionality available for models estimated with a piecewise constant baseline hazard, so this is perhaps the least preferable choice).
Association structures
As mentioned in the previous section, the dependence between the longitudinal and event submodels is captured through the association structure, which can be specified in a number of ways. The simplest association structure is likely to be
\[
f_{mq}(\boldsymbol{\beta}, \boldsymbol{b}_{im}, \alpha_{mq}; t) = \alpha_{mq} \eta_{im}(t)
\]
and this is often referred to as a current value association structure since it assumes that the log hazard of the event at time \(t\) is linearly associated with the value of the longitudinal submodel’s linear predictor also evaluated at time \(t\). This is the most common association structure used in the joint modelling literature to date. In the situation where the longitudinal submodel is based on an identity link function and normal error distribution (i.e. a linear mixed model) the current value association structure can be viewed as a method for including the underlying “true” value of the biomarker as a time-varying covariate in the event submodel.
However, other association structures are also possible. For example, we could assume the log hazard of the event is linearly associated with the current slope (i.e. rate of change) of the longitudinal submodel’s linear predictor, that is
\[
f_{mq}(\boldsymbol{\beta}, \boldsymbol{b}_{i}, \alpha_{mq}; t) = \alpha_{mq} \frac{d\eta_{im}(t)}{dt}
\]
There are in fact a whole range of possible association structures, many of which have been discussed in the literature [14-16].
The stan_jm modelling function in the rstanarm package allows for the following association structures, which are specified via the assoc argument:
Current value (of the linear predictor or expected value) \[
f_{mq}(\boldsymbol{\beta}, \boldsymbol{b}_{i}, \alpha_{mq}; t) = \alpha_{mq} \eta_{im}(t) \\
f_{mq}(\boldsymbol{\beta}, \boldsymbol{b}_{i}, \alpha_{mq}; t) = \alpha_{mq} \mu_{im}(t)
\]
Current slope (of the linear predictor or expected value) \[
f_{mq}(\boldsymbol{\beta}, \boldsymbol{b}_{i}, \alpha_{mq}; t) = \alpha_{mq} \frac{d\eta_{im}(t)}{dt} \\
f_{mq}(\boldsymbol{\beta}, \boldsymbol{b}_{i}, \alpha_{mq}; t) = \alpha_{mq} \frac{d\mu_{im}(t)}{dt}
\]
Area under the curve (of the linear predictor or expected value) \[
f_{mq}(\boldsymbol{\beta}, \boldsymbol{b}_{i}, \alpha_{mq}; t) = \alpha_{mq} \int_0^t \eta_{im}(u) du \\
f_{mq}(\boldsymbol{\beta}, \boldsymbol{b}_{i}, \alpha_{mq}; t) = \alpha_{mq} \int_0^t \mu_{im}(u) du
\]
Interactions between different biomarkers \[
f_{mq}(\boldsymbol{\beta}, \boldsymbol{b}_{i}, \alpha_{mq}; t) = \alpha_{mq} \eta_{im}(t) \eta_{im'}(t)
\text{ for some } m = m' \text{ or } m \neq m' \\
f_{mq}(\boldsymbol{\beta}, \boldsymbol{b}_{i}, \alpha_{mq}; t) = \alpha_{mq} \eta_{im}(t) \mu_{im'}(t)
\text{ for some } m = m' \text{ or } m \neq m' \\
f_{mq}(\boldsymbol{\beta}, \boldsymbol{b}_{i}, \alpha_{mq}; t) = \alpha_{mq} \mu_{im}(t) \mu_{im'}(t)
\text{ for some } m = m' \text{ or } m \neq m'
\]
Interactions between the biomarker (or it’s slope) and observed data \[
f_{mq}(\boldsymbol{\beta}, \boldsymbol{b}_{i}, \alpha_{mq}; t) = \alpha_{mq} c_{i}(t) \eta_{im}(t)
\text{ for some covariate value } c_{i}(t) \\
f_{mq}(\boldsymbol{\beta}, \boldsymbol{b}_{i}, \alpha_{mq}; t) = \alpha_{mq} c_{i}(t) \mu_{im}(t)
\text{ for some covariate value } c_{i}(t) \\
f_{mq}(\boldsymbol{\beta}, \boldsymbol{b}_{i}, \alpha_{mq}; t) = \alpha_{mq} c_{i}(t) \frac{d\eta_{im}(t)}{dt}
\text{ for some covariate value } c_{i}(t) \\
f_{mq}(\boldsymbol{\beta}, \boldsymbol{b}_{i}, \alpha_{mq}; t) = \alpha_{mq} c_{i}(t) \frac{d\mu_{im}(t)}{dt}
\text{ for some covariate value } c_{i}(t)
\]
As well as using lagged values for any of the above. That is, replacing \(t\) with \(t-u\) where \(u\) is some lag time, such that the hazard of the event at time \(t\) is assumed to be associated with some function of the longitudinal submodel parameters at time \(t-u\).
Lastly, we can specify some time-independent function of the random effects, possibly including the fixed effect component. For example,
\[
f_{mq}(\boldsymbol{\beta}, \boldsymbol{b}_{i}, \alpha_{mq}; t) = \alpha_{mq} \boldsymbol{b}_{im0}
\]
or
\[
f_{mq}(\boldsymbol{\beta}, \boldsymbol{b}_{i}, \alpha_{mq}; t) =
\alpha_{mq} \Big( \boldsymbol{\beta}_{m0} + \boldsymbol{b}_{im0} \Big)
\]
where \(\boldsymbol{\beta}_{m0}\) is the population-level intercept for the \(m^{th}\) longitudinal submodel and \(\boldsymbol{b}_{im0}\) is the \(i^{th}\) individual’s random deviation from the population-level intercept for the \(m^{th}\) longitudinal submodel.
Note that more than one association structure can be specified, however, not all possible combinations are allowed. Moreover, if you are fitting a multivariate joint model (i.e. more than one longitudinal outcome) then you can optionally choose to use a different association structure(s) for linking each longitudinal submodel to the event submodel. To do this you can pass a list of length \(M\) to the assoc argument.
Assumptions
Here we define a set of assumptions for the multivariate shared parameter joint model.
The so-called conditional independence assumption of the shared parameter joint model postulates
\[
y_{im}(t) \perp y_{im'}(t) \mid \boldsymbol{b}_i, \boldsymbol{\theta} \\
y_{im}(t) \perp y_{im}(t') \mid \boldsymbol{b}_i, \boldsymbol{\theta} \\
y_{im}(t) \perp T_i^* \mid \boldsymbol{b}_i, \boldsymbol{\theta}
\]
for some \(m \neq m'\) and \(t \neq t'\), and where \(\boldsymbol{\theta}\) denotes the combined vector of all remaining population-level parameters in the model. That is, conditional on the individual-specific parameters \(\boldsymbol{b}_i\) and population-level parameters \(\boldsymbol{\theta}\), the following are assumed: (i) any biomarker measurement for individual \(i\) is independent of that individual’s true event time \(T_i^*\); (ii) any two measurements of the \(m^{th}\) biomarker taken on the \(i^{th}\) individual at two distinct time points \(t\) and \(t'\) (i.e. longitudinal or repeated measurements) are independent of one another; and (iii) any two measurements of two different biomarkers, taken on the \(i^{th}\) individual at some time point \(t\) are independent of one another. These conditional independence assumptions allow for a convenient factorisation of the full likelihood for joint model into the likelihoods for each of the component parts (i.e. the likelihood for the longitudinal submodel, the likelihood for the event submodel, and the likelihood for the distribution of the individual-specific parameters), which facilitates the estimation of the model.
Moreover, we require two additional assumptions: (i) that the censoring process for the event outcome is independent of the true event time, that is \(C_i \perp T_i^* \mid \boldsymbol{\theta}\) (i.e. uninformative censoring); and (ii) that the visiting process by which the observation times \(t_{ijm}\) are determined is independent of the true event time \(T_i^*\) and all missing future unobserved longitudinal biomarker measurements.
Log posterior distribution
Under the conditional independence assumption, the log posterior for the \(i^{th}\) individual can be specified as
\[
\log p(\boldsymbol{\theta}, \boldsymbol{b}_{i} \mid \boldsymbol{y}_{i}, T_i, d_i)
\propto
\log \Bigg[
\Bigg(
\prod_{m=1}^M
\prod_{j=1}^{n_i}
p(y_{ijm}(t) \mid \boldsymbol{b}_{i}, \boldsymbol{\theta})
\Bigg)
p(T_i, d_i \mid \boldsymbol{b}_{i}, \boldsymbol{\theta})
p(\boldsymbol{b}_{i} \mid \boldsymbol{\theta})
p(\boldsymbol{\theta})
\Bigg]
\]
where \(\boldsymbol{y}_i = \{ y_{ijm}(t); j = 1,...,n_i, m = 1,...,M \}\) denotes the collection of longitudinal biomarker data for individual \(i\) and \(\boldsymbol{\theta}\) denotes all remaining population-level parameters in the model.
We can rewrite this log posterior as
\[
\log p(\boldsymbol{\theta}, \boldsymbol{b}_{i} \mid \boldsymbol{y}_{i}, T_i, d_i)
\propto
\Bigg(
\sum_{m=1}^M
\sum_{j=1}^{n_i}
\log p(y_{ijm}(t) \mid \boldsymbol{b}_{i}, \boldsymbol{\theta})
\Bigg) +
\log p(T_i, d_i \mid \boldsymbol{b}_{i}, \boldsymbol{\theta}) +
\log p(\boldsymbol{b}_{i} \mid \boldsymbol{\theta}) +
\log p(\boldsymbol{\theta})
\]
where \(\sum_{j=1}^{n_{im}} \log p(y_{ijm} \mid \boldsymbol{b}_{i}, \boldsymbol{\theta})\) is the log likelihood for the \(m^{th}\) longitudinal submodel, \(\log p(T_i, d_i \mid \boldsymbol{b}_{i}, \boldsymbol{\theta})\) is the log likelihood for the event submodel, \(\log p(\boldsymbol{b}_{i} \mid \boldsymbol{\theta})\) is the log likelihood for the distribution of the group-specific parameters (i.e. random effects), and \(\log p(\boldsymbol{\theta})\) represents the log likelihood for the joint prior distribution across all remaining unknown parameters.
We can rewrite the log likelihood for the event submodel as
\[
\log p(T_i, d_i \mid \boldsymbol{b}_{i}, \boldsymbol{\theta}) =
d_i * \log h_i(T_i) - \int_0^{T_i} h_i(s) ds
\]
and then use Gauss-Kronrod quadrature with \(Q\) nodes to approximate \(\int_0^{T_i} h_i(s) ds\), such that
\[
\int_0^{T_i} h_i(s) ds \approx \frac{T_i}{2} \sum_{q=1}^{Q} w_q h_i \bigg( \frac{T_i(1+s_q)}{2} \bigg)
\]
where \(w_q\) and \(s_q\), respectively, are the standardised weights and locations (“abscissa”) for quadrature node \(q\) \((q=1,...,Q)\) [17]. The default for the stan_jm modelling function is to use \(Q=15\) quadrature nodes, however if the user wishes, they can choose between \(Q=15\), \(11\), or \(7\) quadrature nodes (specified via the qnodes argument).
Therefore, once we have an individual’s event time \(T_i\) we can evaluate the design matrices for the event submodel and longitudinal submodels at the \(Q+1\) necessary time points (which are the event time \(T_i\) and the quadrature points \(\frac{T_i(1+s_q)}{2}\) for \(q=1,...,Q\)) and then pass these to Stan’s data block. We can then evaluate the log likelihood for the event submodel by simply calculating the hazard \(h_i(t)\) at those \(Q+1\) time points and summing the quantities appropriately. This calculation will need to be performed each time we iterate through Stan’s model block. A simplified example of the underlying Stan code used to fit the joint model can be found in Brilleman et al. (2018) [12].
Model predictions
Before discussing the methods by which we can generate posterior predictions, first let us define some additional relevant quantities. Let \(\mathcal{D} = \{ \boldsymbol{y}_i, T_i, d_i; i = 1,...,N \}\) be the entire collection of outcome data in the sample. We will refer to this sample as the “training data”. Let \(T_{max} = \max \{ T_i; i = 1,...,N \}\) denote the maximum event or censoring time across the \(i = 1,...,N\) individuals in our training data.
Individual-specific predictions for in-sample individuals (for \(0 \leq t \leq T_i\))
We can generate posterior predictions for the longitudinal and time-to-event outcomes in the following manner. For the \(i^{th}\) individual in our training data, a predicted value for the \(m^{th}\) longitudinal biomarker at time \(t\), denoted \(y^*_{im}(t)\), can be generated from the posterior predictive distribution
\[
p \Big( y^{*}_{im}(t) \mid \mathcal{D} \Big) =
\int
\int
p \Big( y^{*}_{im}(t) \mid \boldsymbol{\theta}, \boldsymbol{b}_i \Big)
p \Big( \boldsymbol{\theta}, \boldsymbol{b}_i \mid \mathcal{D} \Big)
\space d \boldsymbol{b}_i
\space d \boldsymbol{\theta}
\]
and, similarly, a predicted probability of the \(i^{th}\) individual being event-free at time \(t\), denoted \(S^*_i(t)\), can be generated from the posterior predictive distribution
\[
p \Big( S^{*}_{i}(t) \mid \mathcal{D} \Big) =
\int
\int
p \Big( S^{*}_i(t) \mid \boldsymbol{\theta}, \boldsymbol{b}_i \Big)
p \Big( \boldsymbol{\theta}, \boldsymbol{b}_i \mid \mathcal{D} \Big)
d \boldsymbol{b}_i \space d \boldsymbol{\theta}
\]
Note that for simplicity we have ignored the implicit conditioning on covariates; \(\boldsymbol{x}_{im}(t)\) and \(\boldsymbol{z}_{im}(t)\), for \(m = 1,...,M\), and \(\boldsymbol{w}_{i}(t)\). Since individual \(i\) is included in the training data, it is easy for us to approximate these posterior predictive distributions by drawing from \(p(y^{*}_{im}(t) \mid \boldsymbol{\theta}^{(l)}, \boldsymbol{b}_i^{(l)})\) and \(p(S^{*}_i(t) \mid \boldsymbol{\theta}^{(l)}, \boldsymbol{b}_i^{(l)})\) where \(\boldsymbol{\theta}^{(l)}\) and \(\boldsymbol{b}_i^{(l)}\) are the \(l^{th}\) \((l = 1,...,L)\) MCMC draws from the joint posterior distribution \(p(\boldsymbol{\theta}, \boldsymbol{b}_i \mid \mathcal{D})\).
These draws from the posterior predictive distributions can be used for assessing the fit of the model. For example,
the draws from \(p(y^{*}_{im}(t) \mid \mathcal{D})\) for \(0 \leq t \leq T_i\) can be used to evaluate the fit of the longitudinal trajectory for the \(m^{th}\) biomarker for the \(i^{th}\) individual, and
the draws from \(p(S^{*}_{i}(t) \mid \mathcal{D})\) for \(0 \leq t \leq T_{max}\) can be averaged across the \(N\) individuals to obtain a standardised survival curve (discussed in greater detail in later sections) which can then be compared to the observed survival curve, for example, the Kaplan-Meier curve.
Individual-specific predictions for in-sample individuals (for \(t > C_i\))
However, given that we know the event or censoring time for each individual in our training data, it may make more sense to consider what will happen to censored individuals in our study when we look beyond their last known survival time (i.e. extrapolation).
For an individual \(i\), who was in our training data, and who was known to be event-free up until their censoring time \(C_i\), we wish to draw from the conditional posterior predictive distribution for their longitudinal outcome at some time \(t > C_i\), that is
\[
p \Big( y^{*}_{im}(t) \mid \mathcal{D}, t > C_i \Big) =
\int
\int
p \Big( y^{*}_{im}(t) \mid \boldsymbol{\theta}, \boldsymbol{b}_i, t > C_i \Big)
p \Big( \boldsymbol{\theta}, \boldsymbol{b}_i \mid \mathcal{D} \Big)
d \boldsymbol{b}_i \space d \boldsymbol{\theta}
\]
and the conditional posterior predictive distribution for their survival probability at some time \(t > C_i\), that is
\[
\begin{aligned}
p \Big( S^{*}_{i}(t) \mid \mathcal{D}, t > C_i, T_i^* > C_i \Big)
& =
\frac
{p \Big( S^{*}_{i}(t) \mid \mathcal{D} \Big)}
{p \Big( S^{*}_{i}(C_i) \mid \mathcal{D} \Big)} \\
& =
\int
\int
\frac
{p \Big( S^{*}_i(t) \mid \boldsymbol{\theta}, \boldsymbol{b}_i \Big)}
{p \Big( S^{*}_i(C_i) \mid \boldsymbol{\theta}, \boldsymbol{b}_i \Big)}
\space p \Big( \boldsymbol{\theta}, \boldsymbol{b}_i \mid \mathcal{D} \Big)
d \boldsymbol{b}_i \space d \boldsymbol{\theta}
\end{aligned}
\]
These draws from the conditional posterior predictive distributions can be used to extrapolate into the future for individual \(i\), conditional on their longitudinal biomarker data collected between baseline and their censoring time \(C_i\). For example,
the draws from \(p(y^{*}_{im}(t) \mid \mathcal{D}, t > C_i)\) for \(C_i \leq t \leq T_{max}\) can be used to show the forecasted longitudinal trajectory for the \(m^{th}\) biomarker for the \(i^{th}\) individual, and
the draws from \(p(S^{*}_{i}(t) \mid \mathcal{D}, t > C_i, T_i^* > C_i))\) for \(C_i \leq t \leq T_{max}\) can be used to show the estimated conditional probability of individual \(i\) remaining event-free into the future.
Individual-specific predictions for out-of-sample individuals (i.e. dynamic predictions)
TBC. Describe dynamic predictions under the framework of Rizopoulos (2011) [18]. These types of individual-specific predictions can be obtained using the posterior_traj and posterior_survfit functions by providing prediction data and specifying dynamic = TRUE (which is the default); see the examples provided below.
Population-level (i.e. marginal) predictions
We can also generate posterior predictions for the longitudinal and time-to-event outcomes that do not require any conditioning on observed outcome data for a specific individual. Here, we will discuss two ways in which this can be done.
The first way is to “marginalise” over the distribution of the individual-specific parameters. We wish to generate a predicted value for the \(m^{th}\) longitudinal biomarker at time \(t\) for a new individual \(k\) for whom we do not have any observed data. We will denote this prediction \(y^*_{km}(t)\) and note that it can be generated from the posterior predictive distribution for the longitudinal outcome
\[
\begin{aligned}
p \Big( y^{*}_{km}(t) \mid \mathcal{D} \Big)
& =
\int
\int
p \Big( y^{*}_{km}(t) \mid \boldsymbol{\theta}, \boldsymbol{\tilde{b}}_k \Big)
p \Big( \boldsymbol{\theta}, \boldsymbol{\tilde{b}}_k \mid \mathcal{D} \Big)
\space d \boldsymbol{\tilde{b}}_{k}
\space d \boldsymbol{\theta} \\
& =
\int
\int
p \Big( y^{*}_{km}(t) \mid \boldsymbol{\theta}, \boldsymbol{\tilde{b}}_k \Big)
p \Big( \boldsymbol{\tilde{b}}_k \mid \boldsymbol{\theta} \Big)
p \Big( \boldsymbol{\theta} \mid \mathcal{D} \Big)
\space d \boldsymbol{\tilde{b}}_{k}
\space d \boldsymbol{\theta}
\end{aligned}
\]
and similarly for the survival probability
\[
\begin{aligned}
p \Big( S^{*}_{k}(t) \mid \mathcal{D} \Big)
& =
\int
\int
p \Big( S^{*}_k(t) \mid \boldsymbol{\theta}, \boldsymbol{\tilde{b}}_k \Big)
p \Big( \boldsymbol{\theta}, \boldsymbol{\tilde{b}}_k \mid \mathcal{D} \Big)
d \boldsymbol{b}_k \space d \boldsymbol{\theta} \\
& = \int
\int
p \Big( S^{*}_k(t) \mid \boldsymbol{\theta}, \boldsymbol{\tilde{b}}_k \Big)
p \Big( \boldsymbol{\tilde{b}}_k \mid \boldsymbol{\theta} \Big)
p \Big( \boldsymbol{\theta} \mid \mathcal{D} \Big)
d \boldsymbol{b}_k \space d \boldsymbol{\theta} \\
\end{aligned}
\]
We can obtain draws for \(\boldsymbol{\tilde{b}}_k\) in the same manner as for the individual-specific parameters \(\boldsymbol{b}_i\). That is, at the \(l^{th}\) iteration of the MCMC sampler we draw \(\boldsymbol{\tilde{b}}_k^{(l)}\) and store it. However, individual \(k\) did not provide any contribution to the training data and so we are effectively taking random draws from the posterior distribution for the individual-specific parameters. We are therefore effectively marginalising over the distribution of the group-specific coefficients when we obtain predictions using the draws \(\boldsymbol{\tilde{b}}_k^{(l)}\) fro \(l = 1,\dots,L\). In other words, we are predicting for a new individual whom we have no information except that they are drawn from the same population as the \(i = 1,...,N\) individuals in the training data. Because these predictions will incorporate all the uncertainty associated with between-individual variation our 95% credible intervals are likely to be very wide. These types of marginal predictions can be obtained using the posterior_traj and posterior_survfit functions by providing prediction data and specifying dynamic = FALSE; see the examples provided below.
The second way is to effectively ignore the group-level structure in the model. That is, to only predict with only the population-level parameters contributing to the model. For example, under a identity link function and normal error distribution (i.e. linear mixed effect longitudinal submodel), we would obtain draws from the distribution \(y^{(l)}_{km}(t) \sim N \Big( \boldsymbol{x}^T_{km}(t) \boldsymbol{\beta}_m^{(l)}, \sigma_m^{(l)} \Big)\) where \(\boldsymbol{\beta}_m^{(l)}\) and \(\sigma_m^{(l)}\) are the population-level parameters and residual error standard deviation, respectively, for the \(l^{th}\) draw of the MCMC samples. However, referring to this as a “marginal” prediction is somewhat misleading since we are not explicitly conditioning on the individual-specific parameters but we are implicitly assuming that we know they are equal to zero with absolute certainty. That is, we are actually drawing from the posterior predictive distribution for the longitudinal outcome
\[
\begin{aligned}
p \Big( y^{*}_{km}(t) \mid \mathcal{D} \Big)
& =
\int
p \Big( y^{*}_{km}(t) \mid \boldsymbol{\theta}, \boldsymbol{b}_k = 0 \Big)
p \Big( \boldsymbol{\theta} \mid \mathcal{D} \Big)
d \boldsymbol{\theta} \\
\end{aligned}
\] and similarly for the survival probability
\[
p \Big( S^{*}_{k}(t) \mid \mathcal{D} \Big) =
\int
p \Big( S^{*}_k(t) \mid \boldsymbol{\theta}, \boldsymbol{b}_k = 0 \Big)
p \Big( \boldsymbol{\theta} \mid \mathcal{D} \Big)
d \boldsymbol{\theta} \\
\]
These types of so-called “marginal” predictions can not currently be obtained using the posterior_traj and posterior_survfit functions.
Standardised survival probabilities
All of the previously discussed population-level (i.e. marginal) predictions assumed implicit conditioning on some covariate values for the longitudinal submodel, \(\boldsymbol{x}_{im}(t)\) and \(\boldsymbol{z}_{im}(t)\) for \(m = 1,...,M\), and for the event submodel, \(\boldsymbol{w}_{i}(t)\). Even though we marginalise over the distribution of the individual-specific parameters we were still assuming that we obtained predictions for some known values of the covariates. However, sometimes we wish to marginalise (i.e. average) over the observed distribution of covariates as well. Here we discuss a method by which we can do that for the predicted survival probabilities.
At any time \(t\), it is possible to obtain a standardised survival probability by averaging the individual-specific survival probabilities. That is, we can obtain
\[
S^*(t) = \frac{\sum_{i=1}^{N^{pred}} S_i^*(t)}{N^{pred}}
\]
where \(S_i^*(t)\) is the predicted survival probability for individual \(i\) (\(i = 1,\dots,N^{pred}\) at time \(t\), and \(N^{pred}\) is the number of individuals included in the prediction dataset. We refer to these predictions as standardised survival probabilities.
Note however, that if \(N_{pred}\) is not sufficiently large (e.g. we pass new data with just 2 individuals, say) then marginalising over their covariate distribution may not be meaningful and, similarly, their joint random effects distribution may be a poor representation of the random effects distribution for the entire population. It is better to calculate these standardised survival probabilities using where, say, \(N^{pred}\) is equal to the total number of individuals in the training data.
