Dynamic predictions from longitudinal CD4 count measures and time to death of HIV/AIDS patients using a Bayesian joint model

A Bayesian joint modeling approach to dynamic prediction of HIV progression and mortality allows individualized predictions to be made for HIV patients, based on monitoring of their CD4 counts. This study aims to provide predictions of patient-specific trajectories of HIV disease progression and survival. Longitudinal data on 254 HIV/AIDS patients who received ART between 2009 and 2014, and who had at least one CD4 count observed, were employed in a Bayesian joint model of disease progression. Different forms of association structure that relate the longitudinal CD4 biomarker and time to death were assessed; and predictions were averaged over the different models using Bayesian model averaging. The individual follow-up times ranged from 1 to 120 months, with a median of 22 months and IQR 7–39 months. The estimates of the association structure parameters from two of the three models considered indicated that the HIV mortality hazard at any time point is associated with the rate of change in the underlying value of the CD4 count. Model averaging the dynamic predictions resulted in only one of the hypothesized association structures having non-zero weight in almost all time points for each individual, with the exception of twelve patients, for whom other association structures were preferred at a few time points. The predictions were found to be different when we averaged them over models than when we derived them from the highest posterior weight model alone. The model with highest posterior weight for almost all time points for each individual gave an estimate of the association parameter of –0.02 implying that for a unit increase in the CD4 count, the hazard of HIV mortality decreases by a factor (hazard ratio) of 0.98. Functional status and alcohol intake are important contributing factors that affect the mean square root of CD4 measurements.


Introduction
Globally new HIV infections have started to decrease steadily over the last 10 years. International effort s to improve HIV prevention and treatment services also minimize HIV transmission. As with AIDS-related mortality, the decrease in new HIV infections between 2010 and 2018 was the highest in Eastern and Southern Africa (28% decline) [ 1 , 2 ]. Global decreases in AIDS-related deaths have largely been seen in sub-Saharan Africa, especially in Eastern and Southern Africa, which consists 53% of the world's population live with HIV. AIDS-related mortality has decreased by 42% between 2010 and 2017 in Eastern and Southern Africa, indicating the rapid rate of increase in treatment in the region [ 3 , 4 ].
Ethiopia has been severely affected by the HIV epidemic for the last three decades and there are an increasing number of people living with HIV/AIDS and taking antiretroviral therapy (ART). Although overall estimates of HIV prevalence in the general population remain low, with the 2016 Ethiopian Demographic and Health Survey reporting a prevalence of adult HIV of 0.9%, there are major differences by region (4.8% in Gambella, 3.4% in Addis Ababa, and 0.4% in SNNPR), and by type of region (2.9% urban versus 0.4% rural). According to the latest UNAIDS Spectrum-derived estimates, 613,803 people were living with HIV (PLHIV) in Ethiopia by the end of 2017, of whom 443,213 (72%) were on treatment [5] .
There has been an increasing interest in personalized medical research in recent years, including in the HIV field. Monitoring of a biomarker of disease progression can allow doctors to customize treatment decisions to the characteristics of the patient, to improve medical care and survival [6] . For example, monitoring a HIV patient's CD4 counts allows a clinician to decide when to start and adjust anti-retroviral treatment [7] . However, analysing a longitudinal biomarker measuring disease progression, such as CD4 count, and a survival outcome separately can lead to biased estimates of both the biomarker and survival processes, as such an analysis ignores the dependence between the repeated longitudinal measurements and survival. Joint modelling of longitudinal and survival data is preferred to separate analyses both to optimally use the available information and to obtain unbiased estimates of parameters describing both processes [ 8 , 9 ]. Joint models are versatile methods for deriving probabilities of survival and forecasts for future levels of biomarkers [10] . There is a long history of the use of joint models of CD4 counts and survival in the HIV literature [ 7 , 11-17 ],among these studies, the most common structure assumed for the association between the longitudinal CD4 count process and survival is a random effects structure, with random intercepts and random slopes for the longitudinal CD4 counts, where survival is regressed on these random effects in the joint model [ 11 , 13 ]. A few have used more flexible structures, including integrated Ornstein-Uhlenbeck processes [18] or splines [15] . Here, following Rizopoulos et al. [19] , we consider that different association structures may be better for predicting biomarker and survival processes at different times for different individuals.
A Bayesian approach to joint modeling can be used to derive dynamic individualized predictions of disease progression and survival [6] , and to account for the potential for multiple competing association structures to be preferred for prediction at different time points and individuals, Rizopoulos et al. [19] used Bayesian Model Averaging (BMA) to average across dynamic predictions from the competing models. Motivated by the work of Rizopoulos et al. [19] , we aim to investigate whether different association structures predict better the longitudinal CD4 count and survival processes for different individuals at different times in their follow-up, for a cohort of people living with HIV who are on anti-retroviral treatment. We derive individual-level predictions of both CD4 count progression and HIV survival, based on a collection of possible models simultaneously, and consider whether combining them using Bayesian model averaging is valuable for this cohort or not.

Data, setting and participants
This study used data from HIV patients from Jimma University Specialized Hospital, Ethiopia, in order to predict the probabilities of HIV-related survival over time by jointly modeling the longitudinal CD4 count representing HIV disease progression and time-to-death processes. While the cohort comprised a total of 854 patients who first received ART between 2009 and 2014, the specific HIV data for this analysis came from a longitudinal study of the subset of 254 patients aged at least 18 years who had at least one measurement of CD4 count. 96 subjects (38%) had only a single CD4 count: given this large proportion of a small sample size, these patients were retained in the analysis population, to maintain power. These patients were followed up for a maximum of 48 months, with visit times at which CD4 counts were recorded occurring approximately every 6 months starting from the first time each of these patients received ART. A histogram of observed CD4 counts in Fig. C1 (Appendix C) shows the data fail to fulfil a normality assumption: logarithmic and square root transformations were considered to achieve normality, but the square-root-transformed CD4 counts appeared more Gaussian than the log-transformed counts ( Figs. C2 and C3 , Appendix C). In what follows, we therefore use the square root transformation.

Outcome
The two outcome variables considered for this study were the survival outcome, i.e., time from attendance date to death, measured in months, and the observed CD4 counts, measured in cells/mm 3 of blood. Of the study population, 13% had died, 16% had transferred their care to another facility, 22% were lost to follow-up and 49% were still under active follow-up, by the end of the study. The World Health Organization (WHO)'s clinical stage of disease [20] , the functional status of patients (ambulatory, bedridden or working), age, sex, weight, alcohol use, smoking, drug use and marital status were the covariates available for analysis. The covariate weight is time-dependent and other remaining covariates are fixed at baseline values. The alcohol, smoking and drug use covariates are binary variables for use or not. Further details about the data are found in Temesgen et al. [16] .

Statistical analysis
Since HIV survival is known to be dependent on disease progression, such as measured by CD4 count [ 7 , 14 , 16 ], careful consideration of the statistical method used to relate the two outcome variables is important. A key characteristic of HIV disease progression is its dynamic nature: the rate of progression is not only different from patient to patient, but also dynamically changes in time for the same patient. Thus, the true potential of the CD4 count biomarker in describing disease progression and its association with survival can only be exploited when repeated measurements of CD4 count are considered in the analysis. The structure of the dependence between the two outcomes is not fully known, and may vary between populations, and even over time within a single patient. To address research questions involving characterization of the association structures between repeated measures and event times, a class of statistical models has been developed known as joint models for longitudinal and time-to-event data [ 6 , 21 ]. Briefly, a mixed effects model is proposed for the longitudinal biomarker observations, of a general form y i (t) = m i (t) + ε i (t) where the mean m i ( t ) is the (linear) predictor comprising both fixed and random effects, and ε i ( t ) is a normally distributed error term. Simultaneously, a standard survival model is posited for the time-to-event data, regressed both on covariates and on the mean m i ( t ) of the longitudinal process. Different forms of the regression on m i ( t ) are possible, including regressing only on the current value of the mean, regressing on both the current value and rate of change, or instead regressing on the random effects that are included in m i ( t ), for example. Given the random effects in m i ( t ), both the longitudinal and survival processes are assumed independent, as are the longitudinal responses of each individual. The random effects therefore account both for the association between the longitudinal and the survival outcomes and the correlation between the repeated measurements in the longitudinal process.
Extensions of joint models such as dynamic predictions and accuracy measures have also been implemented [ 15 , 19 , 22 ]. Dynamic prediction is a method for updating predictions ahead in time of both the longitudinal and survival processes, whenever a new measurement of the longitudinal biomarker is taken. Here three joint models of the evolution of the CD4 count process and HIV survival are fitted to the data in a Bayesian framework, each with a different association structure. Dynamic predictions are derived from each of the three models and are combined using Bayesian model averaging [ 19 , 23 ]. The Bayesian approach to joint modelling [ 24 , 25 ] was implemented using the JMbayes package in R version 1.2.5033 that implements both the joint models and the Bayesian model averaging. Full details on the formulation of these joint models and dynamic prediction can be found in Appendix A, but are briefly summarized below.
In a preliminary step, covariate selection for the sub-model for the longitudinal square-root CD4 process was carried out, fitting linear mixed models in a frequentist framework and using likelihood ratio tests, Bayesian Information Criterion (BIC) and Akaike's Information Criterion (AIC) to choose between models (results not shown). The chosen sub-model included a flexible specification of the subject-specific square-root CD4 trajectories, using natural cubic splines of time to capture the non-linear trajectories ( Figs. 1 ; 2 ). In addition to the spline of time, significant covariates retained in the linear mixed model were the functional status of patients (ambulatory, bedridden or working) and alcohol consumption, with an interaction between functional status and the spline of time: where y i (t) is the square-root CD4 count for individual i at time t , B n (t) for n = 1, 2 denotes the B-spline basis for a natural cubic spline with 2 degrees of freedom; FNS and Alcohol denote the variables functional status and alcohol intake respectively; and the term VT * FNS denotes the interaction between the 2df spline of patient visit times and functional status.
The β's are fixed effects whereas the b ' s are random effects. For the survival process, after an initial covariate selection step which resulted in functional status, alcohol intake, marital status and weight as significant covariates for survival, we consider three relative risk models, each positing a different association structure between the two processes, namely: where MS and wt are the marital status and weight variables; the baseline hazard h 0 (t) is approximated with splines (see Eq. (A3) of Appendix A), m i (t) is the current true value of the CD4 count trajectory, m i (t) is the slope of the trajectories at time t (rate of change in CD4 count), γ are regression parameters of the survival model and α are parameters describing the strength of the association between the CD4 count and survival processes.

Individualized dynamic predictions
Based on each joint model, prediction of survival probabilities and future CD4 counts for a new individual j who has a set of longitudinal square-root CD4 counts Y j (t) = { y j (s ) ; 0 ≤ s ≤ t } and a vector of baseline covariates w j is required. For any time u > t , the focus of interest is in predicting both the conditional probability π j (u | t ) that subject j will survive at least up to u and his/her predicted CD4 count at u . At each time of interest (e.g. a clinic visit time) t' , t < t < u , these predictions are dynamically updated, as extra information is recorded for the patient. That is, the prediction ω j (u | t ) of the square-root CD4 count y j (u ) that is based on the information available up to time t , can be updated at time t , to produce a new prediction ω j (u | t ) that uses the additional longitudinal information up to the latter time point t' . Under the Bayesian joint modelling framework, both predictions π j ( u | t ) and ω j (u | t ) are based on the posterior predictive distribution, as given in Appendix A.
Standard model selection techniques for choosing between the three association structures may prefer different models, depending which model selection criteria is used [ 26 , 27 ], particularly in contexts where different association structures may produce better predictions for different individuals at different time points. BMA [ 19 , 22 ] explicitly takes into account model uncertainty by applying Bayesian inference to model selection. Each model is given a prior weight, in this case assuming each association structure is equally likely, and the resulting posterior model weights are used to average over the estimates.
Here, following Rizopoulos et al., instead of averaging estimates over the association structures, the predictions π j (u | t ) and ω j ( u | t ) are averaged over the different association structures. In some contexts, this BMA approach can produce less risky predictions via a straightforward model choice criteria [28] .
Parameter estimates from the three joint models Table 1 shows posterior mean estimates and their corresponding 95% credible intervals for the parameters in the longitudinal and survival sub-models for each of the three different association structures. In the models including current CD4 count in the association structure (Models 1 and 2), the association parameters indicate that HIV survival at any time t is not significantly associated with the current underlying CD4 count at the same time point, whereas the rate of change in CD4 count over time is significantly associated with HIV survival (Models 2 and 3). The parameter estimates in the relative risk models and in the linear mixed models show slight variability between the posited association structures. Model 3, having an association structure which assumes only the random effects are shared between the two processes, has the lowest Deviance Information Criterion (DIC) where a lower DIC implies a better fitting model relative to models with higher DIC. In contrast, in the Bayesian model-averaged predictions, model 1, where the subject-specific linear predictor from the mixed model (current value of CD4 count) is included as a time-varying covariate in the relative risk model is the only model with almost all non-zero posterior weights for all time points for each individual ( Table B2 , Appendix B). The following table ( Table 2 ) shows the 12 patients who have non-zero posterior weights for Model 2 and Model 3 at some time points. Given the preference of the model averaging approach for Model 1 in predicting CD4 counts and HIV survival, for the majority of patients, we interpret here only the parameter estimates for Model 1.
Convergence of the Markov chain Monte Carlo algorithm for this selected Model 1 is demonstrated in Appendix C ( Fig. C4 ). The Model 1 regression coefficient estimates for the CD4 count process suggest that subjects whose functional status was working were associated with a higher CD4 count ( β = 1.84; 95% CrI: 0.10, 3.51) whereas those whose functional status was bedridden were associated with lower CD4 count ( β = −3.09; 95% CI: −5.07, −1.05), compared with the baseline ambulatory functional status group. The cubic spline of visit times had positive effects on the square-root CD4 measurements, indicating increasing CD4 counts with time, and interacted significantly with functional status. Functional status, alcohol consumption and weight also had substantial effects on the hazard of HIV mortality: hazard ratio 0.47 for patients able to work and 2.11 for patients who are bedridden compared to those with ambulatory functional status; 1.64 times higher for patients who consume alcohol compared to the non-alcohol drinking group; and hazard ratio 0.98 for each unit increase in weight. A unit increment in the patient-specific current square-root CD4 count decreases the hazard of HIV mortality 0.98 times ( β = −0.02; 95% CrI: −0.07, 0.03). The posterior probability which is close to zero for the association between current CD4 count and mortality may be due to the inclusion of functional status in both the CD4 count process model and the survival model: perhaps this covariate, since it is associated with CD4 count and has interaction with time, already adequately accounts for the time-varying association of health with mortality Table 2 shows model-averaged predictions for twelve individuals, who have non-zero posterior weights for Model 2 and Model 3 at some time points. The table also shows the time-dependent subject-specific BMA weights for the three joint models. The dynamic predictions are updated at each visit time for each patient, using the CD4 counts observed up to the current visit time to predict both future CD4 count trajectories and HIV survival. As noted, only Model 1 (which considers the current estimated value of the CD4 count trajectory) has non-zero posterior weight at almost all time points. Figs. 3 and 4 display these averaged predictions for one example patient (id 2127), who has non-zero posterior weights for model 1. This patient's observed and expected CD4 counts increase over the study period.

Model-averaged predictions
The vertical dotted lines in both Figs. 3 and 4 represent the time point of the last square-root CD4 count observation. The left sides of the panels present the observed square root of CD4 counts and the right shows the predicted survival ( Fig. 3 ) and predicted longitudinal trajectory of CD4 count for subject 2127. The dashed lines represent the corresponding 95% point wise posterior predictive intervals. The dynamic survival probabilities and dynamic predictions of CD4 counts for patients 2769 and 3790, who have non-zero posterior weights for model 3 and model 2 respectively, are given in Appendix C ( Fig. C5 -C8 ).

Discussion
We have fitted three joint models incorporating different association structures that relate HIV survival to the longitudinal trajectories of patient CD4 counts; then used Bayesian model averaging techniques to base dynamic predictions on the best model for each individual and time point respectively. The model averaging allows for a precision medicine approach to Table 1 Parameter estimates (posterior mean), 95% credible intervals (CrI, in brackets), two-sided posterior probability (tail area) that the parameter is more extreme than 0 (i.e. p = 2 × min{Pr( θ > 0), Pr( θ < 0)} for each parameter θ, in square brackets) and model comparison results (Deviance Information Criterion, DIC) from the three joint models fitted to the HIV data. Parameter Survival process    individual prediction, by allowing for robustness to model misspecification. For most patients in our study, Model 1, which considers the current estimated value of the CD4 count trajectory, was the only model with non-zero posterior weight at most time points. But a notable proportion (12/254 = 0.047) of patients had visit times at which either Model 2 or Model 3 was preferred for predicting future CD4 counts and survival. Despite this preference for other models a different time points, there is not much difference in prediction using individual models versus the averaging method in our data set. This small magnitude of differences suggests robustness of the predictions to alternative association structures ( Table B3 in the Appendix) for this study. The estimates of the association parameter in Model 1, together with the dynamic survival probability plots, showed that a unit increase in the CD4 count decreases the hazard of HIV mortality, indicating that those with higher CD4 count survive longer. The association parameter has close to zero posterior probability, however, suggesting that other covariates in the survival model already adequately explained the variation in mortality. In particular, since the posterior probability for the association of functional status with both CD4 count and survival is a value close to one, as well as the interaction with time in the CD4 count sub-model, perhaps there is some collinearity between functional status and CD4 count which lead to the lack of significant association with survival. Note that we found an increase in the width of the prediction intervals for the future CD4 counts as time progressed. Note also that Rizopoulos et al. [22] state that differences in prediction performance between different specifications of the association structure might be due to the dependence structure not being correctly specified. Rizopoulos et al. [19] compared dynamic prediction results from five joint models with different parameterizations in the survival sub-model and found that the predicted conditional survival probabilities showed considerable variability between the six parameterizations. As a response to the challenge of different results from different parameterizations, Barrett & Su [15] use a different approach, using a flexible non-or semi-parametric specification of the association structure, rather than a parametric specification, to avoid mis-specifying the dependence structure. Although BMA has some interesting features, it has some potential drawbacks as well, including being computationally intensive, since it requires fitting all the models we wish to average. Other approaches to individualized dynamic prediction include stacking and pseudo-BMA [29] and land marking [30] . A number of other papers have also compared prediction models for dynamic prediction [31][32][33] . This study is limited to a single university hospital; it will of interest to apply the method more widely to data from different hospitals throughout the country.
The main aim of this study was to investigate how dynamic predictions from joint models incorporating different association structures can be combined using Bayesian model averaging, and whether BMA is useful for our particular cohort. The result from this study showed that there is some variability between the three association parameterizations we used in terms of parameter estimates. Many studies use standard likelihood information criteria to decide on the selection of appropriate joint model and make predictions. However, we found that in terms of model fit (balanced by a measure of model complexity), model 3 was preferred; whereas in terms of dynamic prediction, the BMA approach showed model 1 was preferred for the majority of patients. For the subset of 12 patients for whom there was non-zero posterior probability of either model 2 or 3 at different time points, we found that BMA gave almost the same dynamic predictions to those obtained from a single model (Model 1) alone. Even if the averaged predictions are not that different from model-specific predictions in general in this particular dataset, as also recognized in other applications [19] , a single prognostic model may not be adequate for all patients at all times, in which case BMA may provide a good solution to sensitivity of predictions to model specification, accounting for model uncertainty. Investigation of different parameterizations and sensitivity of parameter estimates and predictions to model specification is essential to the provision of individual-level predictions of patient progression and survival, i.e. to a precision medicine approach.

Authors' information
FKM is a Ph.D student at Hawasa University, Ethiopia; AMP is a Senior Investigator Statistician at the MRC Biostatistics Unit, University of Cambridge, United Kingdom, DBB is an assistant professor at Department of Statistics, Bahir Dar University, Ethiopia and ATS is an assistant professor at Department of Statistics, Haramaya University, Dire Dawa, Ethiopia.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Availability of data and materials
The datasets analysed in this study are not publicly available due to patient confidentiality, but are available from the corresponding author on reasonable request.

Ethics approval and consent to participate
Ethical clearance was obtained from Department of Statistics of Jimma University. Personal information was kept confidentially without disclosing to others during data collection from patient cards.

Consent for publication
Not applicable.

Model Specification
We denote T * i as a true event time for the i th subject, T i the observed event time, defined as the minimum of the potential censoring time C i and T * i , and by δ i = I( T * i < C i ) the event indicator. A standard approach is the Cox model (Cox, 1972) [34] , postulating a relative risk model of the form where ν 0 < ν 1 < .... < ν Q denotes knots splitting the time scale, with ν Q being larger than the largest observed time, and ξ q denotes the value of the hazard in the interval (ν q −1 , ν q ] . As the number of the knots increases the specification of the baseline hazard becomes more flexible. The log baseline risk function log h 0 ( t ) is extended to B-spline basis functions for cubic splines for the regression splines model as follows: Where k T = ( k 0 , k 1 , . . . , k m ) ar e the spline coefficients, q denotes the degree of the B-splines basis functions B(.), and m = m + q − 1 ,with m denoting the number of interior knots. According to the piecewise-constant model, increasing the number of knots increases the approximate h 0 (. ) flexibility. However, in both approaches, we need to keep a balance between bias and variance and avoid over fitting. Laird and Ware, 1982 [35] ; Harville, 1997 [36] ; and Verbeke and Molenberghs, 20 0 0 [37] used linear mixed models to model the longitudinal process. In linear mixed models, the observed longitudinal data consist of the measurements: where y i (t) denotes the value of the longitudinal outcome at any particular time point t, x i ( t ) and z i ( t ) denote the timedependent design vectors for the fixed-effects β and for the random effects b i , respectively. C is the variance-covariance matrix for the random effects and ε i (t) is the corresponding error terms that are assumed independent of the random effects with mean zero and variance σ 2 .

Estimation of joint models
The main estimation method proposed for joint models is (semi-parametric) maximum likelihood (Wolfsohn and Tsiatis, 1997 [38] ; Henderson et al., 20 0 0 [21] ; Hsieh et al., 2006 [39] ). Bayesian estimation of joint models using MCMC techniques has been employed by Rizopoulos et al. (2014) [19] whereas Barrett and Su (2017) [15] used maximum penalized likelihood approach that maximizes the likelihood function corresponding to the joint distribution of the longitudinal and time-toevent outcomes. Under the Bayesian approach, estimation of the joint model's parameters proceeds using Markov chain Monte Carlo (MCMC) algorithms. The expression for the posterior distribution of the model parameters is derived under the assumptions that given the random effects ( b i ), both the longitudinal and event time process are assumed independent, and the longitudinal responses of each subject are assumed independent. This means that the random effects account both for the association between the longitudinal and the event outcomes, and the correlation between the repeated measurements in the longitudinal process. Formally, we have T denotes the full parameter vector, with θ t denoting the parameters for the event time outcome, θ y the parameters for the longitudinal outcomes and θ b the unique parameters of the random effects covariance matrix, and y i is the n i x 1 vector of longitudinal responses of the i th subject. And p(.) denotes an appropriate probability density function. Under these assumptions the posterior distribution is analogous to: With ψ il (b i ) and φ denoting the natural and dispersion parameters in the exponential family, respectively, c( ·), a( ·), and d( ·) are known functions specifying the member of the exponential family, and for the survival part With h i (.) given by (1) .The integral in the definition of survival function does not have a closed-form solution, and thus a numerical method must be employed for its evaluation. Standard options are the Gauss-Kronrod and Gauss-Legendre quadrature rules. For the parameters θ we take standard prior distributions. In particular, for the vector of fixed effects of the longitudinal submodel β, for the regression parameters of the survival model γ , for the vector of spline coefficients for the baseline hazard γ h 0 , and for the association parameter α we use independent univariate diffuse normal priors. The penalized version of the B-spline approximation to the baseline hazard can be fitted by specifying for γ h 0 the prior (Lang and Brezger 2004) [40] : where τ h is the smoothing parameter that takes a Gamma(1, 0.005) hyper-prior in order to ensure a proper posterior for γ h 0 , K = T r r + 10 −6 I, where r denotes the r th difference penalty matrix, and ρ(K) denotes the rank of K. For the covariance matrix of the random effects we assume an inverse Wishart prior, and when fitting a joint model with a normally distributed longitudinal outcome, we take an inverse-Gamma prior for the variance of the error terms σ 2 .

Model Selection Techniques
When interest is in comparing non-nested models, information criteria are typically used. The two most commonly used information criteria are the Akaike's Information Criteria (AIC; Akaike, 1974) [26] and the Bayesian Information Criteria (BIC; Schwarz, 1978) [27] ( Figs. C6 , C7 ).
where p denotes the number of parameters in the model and n denote the number of observations.
We focus here on dynamic BMA predictions of survival probabilities. BMA predictions for the longitudinal outcome can be produced with similar methodology. Here we assume that we have available data D n = { T i , δ i , y i ; i = 1 , ..., n } based on which we fit M 1 , ..., M K joint models with different association structures. Interest is in calculating predictions for a new subject j from the same population who has provided a set of longitudinal measurements Y j (t) , and has a vector of baseline covariates w j . We let D j (t) = { T * J > t, Y j (t) , w t } denote the available data for this subject. The model-averaged probability of subject j surviving time u > t , given survival up to t is given by the expression: The first term in the right-hand side of Eq. (14) denotes the model-specific survival, while the second denotes the posterior weights of each competing joint models. For calculating the model weights, we observe that these are written as (Rizopoulos et al. 2014) [19] : and p(D n | M k ) is defined analogously. The likelihood part p(D n | θ k ) is based on (3) , and Thus, the subject-specific information in the model weights at time t comes from the available longitudinal measurements p(Y j (t) but also from the fact that this subject has survived up till this time.

Table B3
Subject-specific predictions for the survival and longitudinal outcomes at the final visit time and next three time points for the three joint models with non-zero posterior weights for model 2 and model 3.