Dynamic predictions of kidney graft survival in the presence of longitudinal outliers

TL;DR

This paper introduces a robust Bayesian joint model with t-distributed random effects and errors to improve dynamic survival predictions for kidney transplant graft failure, especially in the presence of outliers.

Contribution

It extends traditional Gaussian joint models by incorporating t-distributions for robustness and evaluates prediction accuracy on real transplant data.

Findings

Robust joint models outperform Gaussian models in outlier scenarios.

Enhanced calibration and discrimination in survival predictions.

Application to kidney transplant data demonstrates practical utility.

Abstract

Dynamic predictions of survival outcomes are of great interest to physicians and patients, since such predictions are useful elements of clinical decision-making. Joint modelling of longitudinal and survival data has been increasingly used to obtain dynamic predictions. A common assumption of joint modelling is that random-effects and error terms in the longitudinal sub-model are Gaussian. However, this assumption may be too restrictive, e.g. in the presence of outliers as commonly encountered in many real-life applications. A natural extension is to robustify the joint models by assuming more flexible distributions than Gaussian for the random-effects and/or error terms. Previous research reported improved performance of robust joint models compared to the Gaussian version in terms of parameter estimation, but dynamic prediction accuracy obtained from such approach has not been yet evaluated. In this study, we define a general robust joint model with t-distributed random-effects and error terms under a Bayesian paradigm. Dynamic predictions of graft failure were obtained for kidney transplant recipients from the French transplant cohort, DIVAT. Calibration and discrimination performances of Gaussian and robust joint models were compared for a validation sample. Dynamic predictions for two individuals are presented.

Figures

Equations

$$Y_{ij}=Y_i^{\ast }(t_{ij})+Z_{ij}={\mathrm{x}}_{ij}^{\top }\alpha +{\mathrm{d}}_{ij}^{\top }{\mathrm{B}}_i+Z_{ij},$$

$$h_i(t)=h_0(t)\exp ({\mathrm{c}}_i^{\top }\omega +g\left({\mathcal{Y}}_i^{\ast }(t),\eta \right)).$$

$$f(\alpha ,\mathrm{\Sigma },\varphi ,\sigma ,\delta ,\zeta ,\omega ,\eta ,\mathrm{B},\mathrm{V},\mathrm{W}|\mathrm{Y},\mathrm{T},\mathrm{E},\mathrm{x},\mathrm{d},\mathrm{c})\propto$$

$$f(\mathrm{T},\mathrm{E}|\zeta ,\omega ,\eta ,\alpha ,\mathrm{B},\mathrm{x},\mathrm{d},\mathrm{c})\times$$

$$f(\mathrm{B}|\mathrm{V},\mathrm{\Sigma })f(\mathrm{V}|\varphi )f(\mathrm{W}|\delta )\times$$

$$f(\alpha ,\mathrm{\Sigma },\varphi ,\sigma ,\delta ,\zeta ,\omega ,\eta ),$$

$$Y_{ij}|\alpha ,\sigma ,{\mathrm{B}}_i,W_{ij},{\mathrm{x}}_{ij},{\mathrm{d}}_{ij}\sim \mathcal{N}({\mathrm{x}}_{ij}^{\top }\alpha +{\mathrm{d}}_{ij}^{\top }{\mathrm{B}}_i,W_{ij}{\sigma }^2).$$

$$h(T_i{)}^{E_i}S(T_i),$$

$$f(\alpha ,\mathrm{\Sigma },\varphi ,\sigma ,\delta ,\zeta ,\omega ,\eta )=f(\alpha )f(\mathrm{\Sigma })f(\varphi )f(\sigma )f(\delta )f(\zeta )f(\omega )f(\eta ).$$

$${\pi }_k(s,u)$$

$$=1-\mathbb{P}\left(T_k^{\ast }>s+u|T_k^{\ast }>s,{\mathrm{Y}}_k,{\mathrm{x}}_k,{\mathrm{d}}_k,{\mathrm{c}}_k\right),$$

$$\mathbb{P}\left(T_k^{\ast }>s+u|T_k^{\ast }>s,{\mathrm{Y}}_k,{\mathrm{x}}_k,{\mathrm{d}}_k,{\mathrm{c}}_k\right)$$

$${\int }_{\theta }{\int }_{{\mathrm{W}}_k}$$

$$f({\mathrm{B}}_k,V_k,{\mathrm{W}}_k|T_k^{\ast }>s,{\mathrm{Y}}_k,{\mathrm{x}}_k,{\mathrm{d}}_k,{\mathrm{c}}_k,\theta )\times$$

$$f(\theta |\mathrm{Y},\mathrm{x},\mathrm{d},\mathrm{c},\mathrm{T},\mathrm{E})d{\mathrm{B}}_{k}d{V}_{k}d{\mathrm{W}}_{k}d\theta ,$$

$$f(V_k,{\mathrm{W}}_k,{\mathrm{B}}_k|T_k^{\ast }>s,{\mathrm{Y}}_k,{\mathrm{x}}_k,{\mathrm{d}}_k,{\mathrm{c}}_k,\theta )\propto$$

$$f(T_k,E_k|{\mathrm{B}}_k,{\mathrm{x}}_k,{\mathrm{d}}_k,{\mathrm{c}}_k,\theta )\times$$

$$f({\mathrm{B}}_k|V_k,\theta )f(V_k|\theta )f({\mathrm{W}}_k|\theta ),$$

$$\displaystyle\mbox{AUC}_{\pi}(s,u)=\mathbb{P}\left(\pi_{k}(s,u)>\pi_{k^{\prime}}(s,u)\big{|}s<T_{k}^{*}\leq s+u,T_{k^{\prime}}^{*}>s+u\right),$$

$$\text{mbox}{R}_{\pi }^2(s,u)=1-\frac{\text{mbox}{BS}_{\pi }(s,u)}{\text{mbox}{BS}_0(s,u)},$$

$$Y_{ij}=Y_i^{\ast }(t_{ij})+Z_{ij}={\alpha }_1+t_{ij}{\alpha }_2+B_{1i}+t_{ij}{B}_{2i}+Z_{ij},$$

$$\displaystyle\begin{array}[]{l}\displaystyle\smash[b]{h}_{i}(t)=\lambda\nu t^{\nu-1}\exp\Big{(}SCr3_{i}\omega_{1}+AR_{i}\omega_{2}+Age_{i}\omega_{3}+CV_{i}\omega_{4}+ACI_{i}\omega_{5}+\\ \qquad\qquad\qquad\qquad\ \ P2008_{i}\omega_{6}+GR_{i}\omega_{7}+\eta_{1}Y_{i}^{*}(t)+\eta_{2}\frac{\partial Y_{i}^{*}(t)}{\partial t}\Big{)}\end{array}$$

Full text

Keywords: Dynamic prediction; kidney transplantation; longitudinal outliers; predictive accuracy; repeated measures; time-to-event

1 Introduction

In the context of chronic diseases, prediction scores of clinical events have become increasingly popular. These scores may help patients and physicians in a shared decision-making and facilitate the implementation of the P4-medicine (predictive, personalized, preventive and participative) (Flores et al., 2013) in clinical practice. Often measured to assess patients’ health evolution during the follow-up, longitudinal markers can be used to improve time-fixed (static) predictions obtained using only baseline information. Dynamic predictions are therefore defined as updated predictions, whenever any new data become available along the follow-up (Rizopoulos, 2012; Proust-Lima and Blanche, 2014).

In the statistical literature, there is a growing interest in methods to compute dynamic predictions. Among them, joint modelling of longitudinal and survival data is one of the most popular (Proust-Lima and Taylor, 2009; Rizopoulos, 2011; Taylor et al., 2013; Asar et al., 2015). Classically, survival outcomes are modelled using a Cox model with a time-varying frailty term, whereas the repeated measures are modelled using a mixed-effect model with Normally distributed random-effects and error terms. In longitudinal clinical studies, some observations may be highly apart from the others, and two types of outliers may be defined: i) at population level, outlying subjects who do not follow the typical population trend, ii) at individual level, outlying observations that do not follow the typical trajectory of an individual (Pinheiro et al., 2001; Wu, 2009). Gaussian assumption would not give appropriate weights to these individuals or observations, i.e. it would not be robust against outliers (Sutradhar and Ali, 1986; Lange et al., 1989).

In the context of kidney transplantation, a systematic review on predictive models emphasized the need for dynamic predictions (Kaboré et al., 2017). Since longitudinal measures of serum creatinine (SCr) have been demonstrated as associated with kidney graft failure (Fournier et al., 2016), we recently proposed to use them to obtain dynamic predictions of long-term kidney graft failure (Fournier et al., 2019). For the French kidney transplant cohort DIVAT (www.divat.fr), we obtained dynamic predictions of graft survival based on a joint model with Gaussian random-effects and error terms. Nevertheless, such dynamic predictions based on this joint model may be sub-optimal in presence of longitudinal outliers.

The objective of the current paper is to investigate possible impacts of longitudinal outliers on dynamic predictions of long-term kidney graft survival. We compare prognostic accuracies (discrimination and calibration) of the Gaussian and robust joint modelling approaches for the DIVAT patients. For the robust joint modelling, we relax the Gaussian assumption by postulating $t$ distributions for the random-effects and error terms. The robust model we consider is novel in the sense that it postulates independent $t$ distributions for the random-effects and error terms, and error terms are allowed to be independent within a subject. To the best of our knowledge, there is no work in the literature that considered such a joint model. For inference and dynamic predictions, we take a Bayesian paradigm and sample from the joint posterior densities. While several authors reported better performance of robust joint models compared to the Gaussian version in terms of parameter estimation, especially in terms of standard error estimation (Li et al., 2009; Huang et al., 2010; Sungduk and Albert, 2016; Baghfalaki et al., 2013, 2014), the advantages of robust joint modelling regarding individual dynamic predictions are yet to be explored.

The rest of the paper is organised as follows. In Section 2, we give details of the DIVAT data that motivates this work. Sections 3 introduces the general modelling framework, distributional assumptions and inferential procedures. Section 4 describes dynamic predictions for newcomer subjects and accuracy measures to evaluate prognostic capabilities of the models. In Section 5, we present a concrete application to the DIVAT data-set. Section 6 closes the paper.

2 DIVAT data

Data were extracted from the multicentric French kidney transplant cohort DIVAT (French Research Ministry: RC12_0452, last agreement No 13 334, No CNIL for the cohort: 891735). The inclusion criteria are: adult recipients who received a first or second renal graft from a living or heart-beating deceased donor, alive with a functioning graft at 1 year post-transplantation, and maintained under Tacrolimus and Mycofenolate. The aim is to predict graft failure risk in the chronic phase of transplantation. Therefore, time origin is selected as 1 year after transplantation. Graft failure is defined as the occurrence of any of the following events: return to dialysis, pre-emptive re-transplantation, or death with a functioning graft. Serum creatinine (SCr, measured in $\mu$mol/L) is the longitudinal biomarker, with high values indicating worse kidney health. Following Fournier et al. (2019), the following set baseline covariates are considered: recipient age (Age: in years), history of cardiovascular disease (CV: yes/no), 3-month SCr (SCr3), occurrence of an acute rejection episode in the first year post-transplantation (AR: yes/no), pre-transplantation anti-class I immunisation (ACI: yes/no) and graft rank (GR: second/first). Follow-up time $(t_{ij})$ is the time-varying covariate.

The learning and validation samples are constructed based on two data extracts. While the first extract covers the period of 1 January 2000 – 31 August 2013, the second covers 1 September 2013 – 31 October 2016. Randomly chosen two-third of the first extract constitutes the learning sample (2,584 patients), whereas remaining portion of the first extract plus the whole data from the second extract constitute the validation sample (1,577+946 = 2,523 patients). Number of repeated measures ranged between 1 and 14 (with median of 4) for the learning sample. Spagetti-plots SCr (in log-scale) for a random sample of 60 patients from the learning sample are shown in Figure 1. For more details of the data, reader is referred to Fournier et al. (2019). The quantile-quantile plot displayed in Figure 2 indicates that there are considerable departures from the Gaussian assumption for the DIVAT data-set.

3 Joint modelling of longitudinal and survival outcomes

From Figure 2, we are not able to discern the source of heavy-tailedness, i.e. whether it is due to outlying individuals or outlying observations or both. Therefore, in what follows we consider a general modelling framework which postulates that both types of outliers may be present. Based on the posterior summaries for the degrees-of-freedom parameters for random-effects and error terms, we can then drop the heavy-tailedness assumption for the term for which Normality is indicated.

3.1 Notations and framework

To make inference, we assume to observe a sample of $n$ independent and identically distributed subjects with the data of $\left\{T_{i},E_{i},Y_{i},t_{i},a_{i};i=1,\ldots,n\right\}$. Here, $Y_{i}=\left\{Y_{ij};j=1,\ldots,m_{i}\right\}$ is the set of longitudinal marker with the set of corresponding timings, $t_{i}=\left\{t_{ij};j=1,\ldots,m_{i}\right\}$; $a_{i}=\left\{a_{ik};k=1,\ldots,f\right\}$ the baseline explanatory variables; and $T_{i}$ the time elapsed between the origin and occurrence of the survival event. $T_{i}$ is typically right-censored, i.e. $T_{i}=\min(T_{i}^{*},C_{i})$, where $T_{i}^{*}$ being the true survival time and $C_{i}$ the censoring time for subject $i$. Therefore, an event indicator, $E_{i}=\mathbbm{1}\{T_{i}^{*}\leq C_{i}\}$, where $\mathbbm{1}\{\cdot\}$ denotes the indicator function, completes the survival information.

In this study, we consider the so-called shared-parameter version of the joint model Wulfson and Tsiatis (1997), based on two linked sub-models. The modelling framework can be defined as follows:

[TABLE]

Equation (1) corresponds to a linear mixed-effects model that defines the longitudinal process. The observed longitudinal measure, $Y_{ij}$, is assumed to be a noisy version of the underlying continuous-time signal, $Y_{i}^{*}(t)$, at time $t_{ij}$, with $Z_{ij}$ being the noise, or as often called measurement error. ${{\boldsymbol{\mathrm{x}}}}_{ij}=[x_{ij1}\ldots x_{ijp}]^{\top}$ is a design matrix that consists of elements from $a_{i}$ and $t_{i}$, and ${{\boldsymbol{\mathrm{d}}}}_{ij}=[d_{ij1}\ldots d_{ijq}]^{\top}$ another design matrix that is typically structured as a subset of ${{\boldsymbol{\mathrm{x}}}}_{ij}$. ${{\boldsymbol{\mathrm{\alpha}}}}=[\alpha_{1}\ldots\alpha_{p}]^{\top}$ and ${{\boldsymbol{\mathrm{B}}}}_{i}=[B_{i1}\ldots B_{iq}]^{\top}$ are population-averaged parameters and subject-specific random-effects (latent variables), respectively.

Equation (2) corresponds to a Cox model with time-varying frailty for the time-to-event. $h_{i}(t)$ is the instantaneous risk of experiencing the event at time $t$, and $h_{0}(t)$ the baseline hazard. $h_{0}(t)$ can be left un-specified as in Cox (1972), or specified parametrically using hazard function of a life-time distribution, e.g. for Weibull $h_{0}(t)=\lambda\nu t^{(\nu-1)}$, or using piecewise-constants or splines; for details see Rizopoulos (2012). ${{\boldsymbol{\mathrm{c}}}}_{i}=[c_{i1}\ldots c_{ir}]^{\top}$ is a design matrix that is composed of elements from $a_{i}$, and ${{\boldsymbol{\mathrm{\omega}}}}$ are the associated parameters. $g\left(\mathcal{Y}_{i}^{*}(t),{{\boldsymbol{\mathrm{\eta}}}}\right)$ is the term that links features of underlying signal up to and including time $t$, $\mathcal{Y}_{i}^{*}(t)$, and hazard function at time $t$ through a known link function $g(\cdot)$ and parameters ${{\boldsymbol{\mathrm{\eta}}}}$. There are a number of choices for $g\left(\mathcal{Y}_{i}^{*}(t),{{\boldsymbol{\mathrm{\eta}}}}\right)$. Widely used examples include, among others, the current value parametrization, $g\left(\mathcal{Y}_{i}^{*}(t),{{\boldsymbol{\mathrm{\eta}}}}\right)=\eta Y_{i}^{*}(t)$, or current value and rate of change parametrization, $g\left(\mathcal{Y}_{i}^{*}(t),{{\boldsymbol{\mathrm{\eta}}}}\right)=\eta_{1}Y_{i}^{*}(t)+\eta_{2}\frac{\partial Y_{i}^{*}(t)}{\partial t}$; for other parametrizations, see (Rizopoulos, 2012).

3.2 Distributional assumptions

Standard joint models assume that ${{\boldsymbol{\mathrm{B}}}}_{i}$ and $Z_{ij}$ are both zero-mean Gaussian, such that ${{\boldsymbol{\mathrm{B}}}_{i}}|{{\boldsymbol{\mathrm{\Sigma}}}}\sim\mathcal{MVN}\left({{\boldsymbol{\mathrm{0}}}},{{\boldsymbol{\mathrm{\Sigma}}}}\right)$ and $Z_{ij}|\sigma\sim\mathcal{N}(0,\sigma^{2})$. The terms are further assumed to have the following properties, ${{\boldsymbol{\mathrm{B}}}}_{i}\perp Z_{ij}$ and $Z_{ij}\perp Z_{ij^{\prime}}$ for $j\neq j^{\prime}$. Gaussian assumption might be too restrictive for some real-life applications, because the data-sets typically consist of subjects that exhibit outlying behaviours.

Pinheiro et al. (2001) discuss two types of longitudinal outliers:

1. outlying individuals, and 2) outlying observations within individuals. Examples of outlying individuals include subjects with very high/low health status at baseline, and/or with rapid progression. In other words, outlying individuals correspond to outliers in ${{\boldsymbol{\mathrm{B}}}}_{i}$. Examples of outlying observations might be a few observations that are quite different than the rest of the observations for a given individual. In other words, outlying observations correspond to outliers in $Z_{ij}$. Figure 3 displays simulated realisations from a robust joint model. As expected, most of the subjects are homogeneous, while a few subjects seem to have extreme trajectory or extreme observations. For instance, subjects 10, 26 and 27 seem to have relatively higher slopes compared to the rest. Fourth observation for subject 24 seems to deviate more around the individual line compared to the other observations for the same subject. Similar features can also be seen for the DIVAT data-set: Figure 1 illustrates subjects with high level of SCr at baseline, with rapid progression, and with a few observations that seem different from the rest of the observations for the same subjects.

To accommodate the aforementioned outliers, the Gaussian assumption for ${{\boldsymbol{\mathrm{B}}}}_{i}$ and $Z_{ij}$ can be relaxed using $t$ distribution that would give lower weights to outliers. $t$ distribution can be specified using the variance mixtures as ${{\boldsymbol{\mathrm{B}}}}_{i}=\sqrt{V_{i}}{{\boldsymbol{\mathrm{B}}}}_{i}^{*}$ and $Z_{ij}=\sqrt{W_{ij}}\sigma Z_{ij}^{*}$, where $V_{i}$ and $W_{ij}$ are inverse-Gamma random variables, such that $V_{i}|\phi\sim\mathcal{IG}(\phi/2,\phi/2)$ and $W_{ij}|\delta\sim\mathcal{IG}(\delta/2,\delta/2)$, with the following properties, $V_{i}\perp W_{ij}$, $W_{ij}\perp W_{ij^{\prime}}$ for $j\neq j^{\prime}$, ${{\boldsymbol{\mathrm{B}}}_{i}^{*}}|{{\boldsymbol{\mathrm{\Sigma}}}}\sim\mathcal{MVN}\left({{\boldsymbol{\mathrm{0}}}},{{\boldsymbol{\mathrm{\Sigma}}}}\right)$, and $Z_{ij}^{*}\sim\mathcal{N}(0,1)$. With these specifications, one would obtain ${{\boldsymbol{\mathrm{B}}}_{i}}|{{\boldsymbol{\mathrm{\Sigma}}}},\phi\sim\mathcal{MV}t\left({{\boldsymbol{\mathrm{0}}}},{{\boldsymbol{\mathrm{\Sigma}}}},\phi\right)$ and $Z_{ij}|\sigma,\delta\sim t(0,\sigma^{2},\delta)$, with the following properties, ${{\boldsymbol{\mathrm{B}}}}_{i}\perp Z_{ij}$ and $Z_{ij}\perp Z_{ij^{\prime}}$, for $j\neq j^{\prime}$. Note that the conditionals on the mixing variates are still Normal such that ${{\boldsymbol{\mathrm{B}}}}_{i}|V_{i},{{\boldsymbol{\mathrm{\Sigma}}}}\sim\mathcal{MVN}({{\boldsymbol{\mathrm{0}}}},V_{i}{{\boldsymbol{\mathrm{\Sigma}}}})$ and $Z_{ij}|\sigma,W_{ij}\sim\mathcal{N}(0,W_{ij}\sigma^{2})$. $\phi$ and $\delta$ are the so-called degree-of-freedom parameters, and it is well known that if such a parameter converges to infinity, $t$ distribution converges to the Gaussian.

To the best of our knowledge, there is no work in the literature that considered robust joint modelling with the above properties. Taylor et al. (2013), Li et al. (2009) and Huang et al. (2010) considered robust joint modelling with Normally distributed ${{\boldsymbol{\mathrm{B}}}}_{i}$ and $t$-distributed $Z_{ij}$ with fixed degree-of-freedom parameters. Sungduk and Albert (2016) considered Normally distributed ${{\boldsymbol{\mathrm{B}}}}_{i}$, and generalised $t$-distributed $Z_{ij}$. Baghfalaki et al. (2013, 2014) considered our model with $W_{ij}=W_{i}$. Note that under their specification, the property, $Z_{ij}\perp Z_{ik}$ for $j\neq k$, does not hold, as $W_{i}$ is shared across $Z_{ij}$’s. Among these works, only Taylor et al. (2013) considered dynamic predictions, whereas the rest focused on parameter estimation.

3.3 Bayesian inference

In this section, we present inference for the joint model with $t$-distributed ${{\boldsymbol{\mathrm{B}}}}_{i}$ and $Z_{ij}$ terms. Inference for the model with at least one of these terms being Gaussian, or $Z_{ij}$ being $t$-distributed based on $W_{ij}=W_{i}$, are just special cases.

Let ${{\boldsymbol{\mathrm{Y}}}}=[Y_{1}^{\top}\ldots Y_{n}^{\top}]^{\top}$ with ${{\boldsymbol{\mathrm{Y}}}}_{i}=[Y_{i1}\ldots Y_{im_{i}}]^{\top}$; ${{\boldsymbol{\mathrm{T}}}}=[T_{1}\ldots T_{n}]^{\top}$; ${{\boldsymbol{\mathrm{E}}}}=[E_{1}\ldots E_{n}]^{\top}$; ${{\boldsymbol{\mathrm{x}}}}=[{{\boldsymbol{\mathrm{x}}}}_{1}^{\top}\ldots{{\boldsymbol{\mathrm{x}}}}_{n}^{\top}]^{\top}$ with ${{\boldsymbol{\mathrm{x}}}}_{i}=[{{\boldsymbol{\mathrm{x}}}}_{i1}\ldots{{\boldsymbol{\mathrm{x}}}}_{im_{i}}]^{\top}$ and ${{\boldsymbol{\mathrm{x}}}}_{ij}$ as before; ${{\boldsymbol{\mathrm{d}}}}=[{{\boldsymbol{\mathrm{d}}}}_{1}^{\top}\ldots{{\boldsymbol{\mathrm{d}}}}_{n}^{\top}]^{\top}$ with ${{\boldsymbol{\mathrm{d}}}}_{i}=[{{\boldsymbol{\mathrm{d}}}}_{i1}\ldots{{\boldsymbol{\mathrm{d}}}}_{im_{i}}]^{\top}$ and ${{\boldsymbol{\mathrm{d}}}}_{ij}$ as before; ${{\boldsymbol{\mathrm{B}}}}=[{{\boldsymbol{\mathrm{B}}}}_{1}^{\top}\ldots{{\boldsymbol{\mathrm{B}}}}_{n}^{\top}]^{\top}$ with ${{\boldsymbol{\mathrm{B}}}}_{i}$ as before; ${{\boldsymbol{\mathrm{V}}}}=[V_{1}\ldots V_{n}]^{\top}$; ${{\boldsymbol{\mathrm{W}}}}=[{{\boldsymbol{\mathrm{W}}}}_{1}^{\top}\ldots{{\boldsymbol{\mathrm{W}}}}_{n}^{\top}]^{\top}$ with ${{\boldsymbol{\mathrm{W}}}}_{i}=[W_{i1}\ldots W_{im_{i}}]^{\top}$; ${{\boldsymbol{\mathrm{\zeta}}}}$ the parameters of $h_{0}(t)$; ${{\boldsymbol{\mathrm{c}}}}=[{{\boldsymbol{\mathrm{c}}}}_{1}^{\top}\ldots{{\boldsymbol{\mathrm{c}}}}_{n}^{\top}]^{\top}$ with ${{\boldsymbol{\mathrm{c}}}}_{i}$ as before. The joint posterior density of the parameters and latent variables can be written as

[TABLE]

with $f(\cdot)$ being a general notation for probability density function. The first distribution on the right-hand side of (3.3), $f({{\boldsymbol{\mathrm{Y}}}}|{{\boldsymbol{\mathrm{\alpha}}}},\sigma,{{\boldsymbol{\mathrm{B}}}},{{\boldsymbol{\mathrm{W}}}},{{\boldsymbol{\mathrm{x}}}},{{\boldsymbol{\mathrm{d}}}})$, is based on the longitudinal sub-model (1), and is constructed based on univariate Normal distributions such that

[TABLE]

The second term is based on the survival sub-model (2), and constructed by

[TABLE]

where $S(\cdot)$ being the survival function, defined as $S(t)=\mathbb{P}(T>t)=\exp\left(-\int_{0}^{T_{i}}h_{i}(v)dv\right)$. As mentioned in Section 3.2, $f({{\boldsymbol{\mathrm{B}}}}|{{\boldsymbol{\mathrm{V}}}},{{\boldsymbol{\mathrm{\Sigma}}}})$ is constructed based on $\mathcal{MVN}({{\boldsymbol{\mathrm{0}}},V_{i}{{\boldsymbol{\mathrm{\Sigma}}}}})$, and $f({{\boldsymbol{\mathrm{V}}}}|\phi)$ and $f({{\boldsymbol{\mathrm{W}}}}|\delta)$ are based on inverse-Gamma distributions, $\mathcal{IG}(\phi/2,\phi/2)$ and $\mathcal{IG}(\delta/2,\delta/2)$, respectively. $f({{\boldsymbol{\mathrm{\alpha}}}},{{\boldsymbol{\mathrm{\Sigma}}}},\phi,\sigma,\delta,{{\boldsymbol{\mathrm{\zeta}}}},{{\boldsymbol{\mathrm{\omega}}}},{{\boldsymbol{\mathrm{\eta}}}})$ corresponds to the joint prior distribution of the parameters. We assume independent priors for the parameters such that

[TABLE]

$\alpha_{h}$ $(h=2,\ldots,p)$ are given zero-mean Cauchy prior with scale parameter of 5, $\mathcal{C}(0,5)$, whereas $\alpha_{1}$ is given $\mathcal{C}(0,20)$. ${{\boldsymbol{\mathrm{\Sigma}}}}$ is decomposed as ${{\boldsymbol{\mathrm{R}}}}{{\boldsymbol{\mathrm{\Omega}}}}{{\boldsymbol{\mathrm{R}}}}$, where ${{\boldsymbol{\mathrm{R}}}}$ is diagonal matrix of scale parameters that are specific to $B_{h}$ $(h=1,\ldots,q)$, and ${{\boldsymbol{\mathrm{\Omega}}}}$ is in the form of a correlation matrix. Elements of ${{\boldsymbol{\mathrm{R}}}}$ are given half-Cauchy priors with scale of 5, $\mathcal{C}_{+}(0,5)$, whereas elements of ${{\boldsymbol{\mathrm{\Omega}}}}$ are given LKJ prior with the parameter of 2, $\mathcal{LKJ}(2)$. $\phi$ and $\delta$ are given uniform priors, between 2 and 100. $\sigma$ is given $\mathcal{C}_{+}(0,5)$. Log-transformed elements of ${{\boldsymbol{\mathrm{\zeta}}}}$ and elements of ${{\boldsymbol{\mathrm{\omega}}}}$ and ${{\boldsymbol{\mathrm{\eta}}}}$ are given $\mathcal{C}(0,5)$.

Samples from the joint posterior (3.3) are drawn using HMC (Neal, 2011), specifically using the NUTS algorithm (Hoffman and Gelman, 2014), that is an adaptive version of HMC. Methods are implemented in the R package robjm (github.com/ozgurasarstat/robjm) (R Core Team, 2018), that internally uses the so-called HMC engine Stan (Carpenter et al., 2017) through the RStan package (Stan Development Team, 2018).

4 Dynamic predictions

4.1 Definition

Our target of inference for a newcomer subject $k$ is the prediction of subject-specific conditional failure probability between time points $s$ and $s+u$ given that the patient did not experience the event until time point $s$, i.e. $T_{k}^{*}>s$, and subject-specific data recorded up to and including time $s$ $\left({{\boldsymbol{\mathrm{Y}}}}_{k}=[Y_{k1}\ldots Y_{km_{k}}]^{\top},{{\boldsymbol{\mathrm{x}}}}_{k}=[x_{k1}\ldots x_{km_{k}}]^{\top},{{\boldsymbol{\mathrm{d}}}}_{k}=[d_{k1}\ldots d_{km_{k}}]^{\top},{{\boldsymbol{\mathrm{c}}}}_{k}=[c_{i1}\ldots c_{ir}]^{\top}\right)$:

[TABLE]

where $u>0$ is called the lead-time or forecast horizon. The conditional survival probability in (4.1) could be obtained as

[TABLE]

where ${{\boldsymbol{\mathrm{\theta}}}}$ consists of all the parameters. For the a new subject, we would not have samples from $f({{\boldsymbol{\mathrm{B}}}}_{k},V_{k},{{\boldsymbol{\mathrm{W}}}}_{k}|T_{k}^{*}>s,{{\boldsymbol{\mathrm{Y}}}}_{k},{{\boldsymbol{\mathrm{x}}}}_{k},{{\boldsymbol{\mathrm{d}}}}_{k},{{\boldsymbol{\mathrm{c}}}}_{k},{{\boldsymbol{\mathrm{\theta}}}})$, that are readily obtained during the inference step. We can draw the latent variables, $V_{k},{{\boldsymbol{\mathrm{W}}}}_{k},{{\boldsymbol{\mathrm{B}}}}_{k}$, from

[TABLE]

using the MC samples of the parameters obtained from $[{{\boldsymbol{\mathrm{\theta}}}}|{{\boldsymbol{\mathrm{Y}}}},{{\boldsymbol{\mathrm{x}}}},{{\boldsymbol{\mathrm{d}}}},{{\boldsymbol{\mathrm{c}}}},{{\boldsymbol{\mathrm{T}}}},{{\boldsymbol{\mathrm{E}}}}]$; see (3.3). One would obtain dynamic predictions by updating $\pi_{k}(s,u)$ to $\pi_{k}(s^{\prime},u)$, if subject $k$ is still at risk at time $s^{\prime}$ $(s^{\prime}>s)$ (such that $T_{k}^{*}>s^{\prime}$) by incorporating any new data recorded for her/him in the time interval $(s,s^{\prime}]$.

4.2 Accuracy measures for dynamic predictions

A prediction score, whether dynamic or not, requires good properties of discrimination and calibration in order for having a practical use in personalized medicine (Steyerberg et al., 2010). Methods to assess these properties have already been extensively published (Graf et al., 1999; Gerds and Schumacher, 2006; Heagerty et al., 2000) and have recently been extended to dynamic predictions (Schoop et al., 2008; Blanche et al., 2015; Fournier et al., 2018).

A widely used discrimination measure is the Area Under the Receiver Operating Characteristics Curve (AUC) that aims to assess how the predictions distinguish between a patient who has the event from another patient who does not. In a dynamic prediction context, AUC is calculated for each landmark time-point $s$, with a forecast horizon of $u$, such that

[TABLE]

where $k$ and ${k^{\prime}}$ are indices for two randomly selected subjects. We use the estimator of Blanche et al. (2015) for taking into account right-censoring. Higher AUC values indicate better discrimination.

Brier score is one of the widely used measure to globally assess the prognostic performances. A disadvantage of the Brier score is that it depends on the marginal failure probabilities that could potentially take different values at different landmark times, and thus can be misleading in the dynamic prediction context. Fournier et al. (2018) proposed an $\mbox{R}^{2}$-type criterion that builds on the Brier score by adjusting it with the marginal failure probability such that

[TABLE]

where $\mbox{BS}_{\pi}(s,u)$ is the Brier score and $\mbox{BS}_{0}(s,u)$ the Brier score for the reference model that does not use any subject-specific information. Higher $\mbox{R}^{2}$ indicates better performance.

As shown in Fournier et al. (2018), both Brier score and the $\mbox{R}^{2}$-type criterion measure calibration and discrimination simultaneously. We therefore also use calibration plots to solely check the calibration properties of the predictions. The calibration is described as comparing predicted values within subgroups, e.g. defined from deciles of predictions for the DIVAT application, to observed event survival that is computed using the Kaplan-Meier method. From a Bayesian perspective, we have calculated deciles for each of element of the MC samples of dynamic predictions.

5 Application to the kidney transplantation data

5.1 Joint modelling for the learning sample

Following Fournier et al. (2019), the following joint model is fitted to the learning data:

[TABLE]

where $Y=\log(SCr)$, and $P2008$ stands for transplantation period: before 31 December 2007 versus after January 2008. No baseline covariates were included in (8) since SCr evolution is on the causal pathway between baseline factors and graft failure risk as discussed in Fournier et al. (2019). Besides, it was confirmed in the web supplementary material of Fournier et al. (2019) that inclusion of baseline covariates in the longitudinal sub-model did not improve predictive performances. Note that the current value and rate of change parametrization is suggested by Figure 1: increase in the mean SCr level is higher for the graft failure group compared to the censored group, and for the former group there is an acceleration in the SCr increase towards the event.

In the following, we will distinguish between four joint models for the DIVAT data-set based on the following distributional assumptions: Normally distributed ${{\boldsymbol{\mathrm{B}}}}_{i}$ and Normally distributed $Z_{ij}$ ($N-N$ model), $t$-distributed ${{\boldsymbol{\mathrm{B}}}}_{i}$ and $t$-distributed $Z_{ij}$ ($t-t$), Normally distributed ${{\boldsymbol{\mathrm{B}}}}_{i}$ and $t$-distributed $Z_{ij}$ ($N-t$), $t$-distributed ${{\boldsymbol{\mathrm{B}}}}_{i}$ and Normally distributed $Z_{ij}$ ($t-N$). Our general strategy would be to fit the $t-t$ model when there is evidence against Normal, and check the posterior summaries of the degrees-of-freedom parameters. As mentioned before, this approach is due to the fact that we do not know the source of heavy-tailedness in the standardised marginal residuals (see Figure 2). If results for any of the degrees-of-freedom parameters indicates that Normal assumption is reasonable, we then switch to $N-t$ or $t-N$ model. As will be discussed in the next paragraph, for the DIVAT application, $t-t$ model does not indicate Normality to either of the ${{\boldsymbol{\mathrm{B}}}}_{i}$ or $Z_{ij}$ terms. Therefore, we include the $N-t$ and $t-N$ models for the sake of seeing the likely effects of wrongly assuming any of ${{\boldsymbol{\mathrm{B}}}}_{i}$ or $Z_{ij}$ as Normal. For each of the models, 4 chains with lengths of 2,000 were started from random initials. For each chain, first halves were considered as warm-up that results MC samples of 4,000 for each model. Convergence of the chains were checked using trace-plots, density plots for the chains, and R-hat statistic of Brooks and Gelman (1997).

Table 1 presents posterior summaries for each of the 4 joint models, specifically the 2.5%th, 50%th and 97.5%th percentiles, of the MC samples. Results for $\phi$ and $\delta$ under the $t-t$ model imply that both ${{\boldsymbol{\mathrm{B}}}}_{i}$ and $Z_{ij}$ reflect heavier tails than the Gaussian. When we let both ${{\boldsymbol{\mathrm{B}}}}_{i}$ and $Z_{ij}$ being Normally distributed (i.e. the $N-N$ model), the population-averaged slope, $\alpha_{2}$, is largely over-estimated, compared to the $t-t$ model. When we only let the error terms being $t$-distributed (i.e. the $N-t$ model), this over-estimation gets milder. The results of $\alpha_{2}$ under the $t-N$ model are the most similar to those of the $t-t$ model. These can be explained as the following: profiles for subjects who had extreme progression (high slopes) are better captured by the $t$-distributed ${{\boldsymbol{\mathrm{B}}}}_{i}$ terms. However, when the error term is assumed to be Normal ($t-N$ model), all the outlying behaviours are forced to be in ${{\boldsymbol{\mathrm{B}}}}_{i}$, hence we obtained somewhat over-estimated $\alpha_{2}$. Additionally, $\eta_{2}$ is largely over-estimated by the models with Normally distributed ${{\boldsymbol{\mathrm{B}}}}_{i}$ regardless of $Z$ being Normal or $t$. This can also be explained by better capturing subject-specific slopes that in turn would results better predictions for $\frac{\partial Y^{*}_{i}(t)}{\partial t}$. Smoothed posterior densities for $\alpha_{2}$ and $\eta_{2}$ are presented in Figure 4. Scatter-plots of the posterior quantiles of $B_{1i}$ and $B_{2i}$ to compare the models are presented in Figures $1-6$ of the supplementary material. Note that the differences between the models are more apparent for $B_{2i}$.

5.2 Dynamic predictions for the validation sample

5.2.1 Accuracy measures

Dynamic predictions are calculated at six post-transplantation landmark times: $s=1,2,3,4,5,6$. A clinically meaningful forecast window of 5 years is considered, i.e. $u=5$. Table 2 is the frequency table regarding subjects who were at risk at the landmark times, subjects who had the event, and who were censored within the forecast horizons, and subjects who survived beyond the forecast horizons.

Posterior summaries of the AUC and R2 are displayed in Figure 5. Calibration plots are dislayed in Figures $7-10$ of the supplementary material, whereas the estimated slopes fit to predicted risk versus observed risk are reported in Table 3. Note that as mentioned in Fournier et al. (2018), $\mbox{R}^{2}$ can take negative values. Globally, discrimination is good for all the models, AUC’s increase with landmark times, and calibration performances appear reasonable. More interestingly, comparing the models, we observe the best results under the $t-N$ model in terms of AUC and $R^{2}$. It is followed by the $t-t$ model and the $N-t$ model. The $N-N$ has the worst performance. In terms of calibration, again the $t-N$ and $t-t$ models appeared the best two. Whereas in landmark times $1-3$, $t-t$ model had higher slope estimates compared to the $t-N$ model, for the later times the estimates are almost the same.

We also calculated 2.5th, 50th and 97.5th percentiles of the dynamic predictions for each individual for each landmark times. Scatter-plots of these statistics are displayed in Figures $11-16$ of the supplementary material. It can be seen, e.g. based on the medians, that there are some subjects for whom the models quite disagree, e.g. see the dots that are above the $x=y$ lines. One can deduce that there are some subjects (e.g. outlying subjects) whose dynamic predictions benefit from the robust joint model with $t$ distributions, especially $t$-distributed ${{\boldsymbol{\mathrm{B}}}}_{i}$ terms.

5.2.2 Results for two patients

We present 5-year dynamic predictions for two patients from the validation sample. They are the same subjects presented in Fournier et al. (2019). The first patient (Figure 6) was a female, aged 51 years, transplanted in 2005 for the first time, without history of a cardiovascular disease, immunized against HLA class I, with a SCr measurement at 3 months post-transplantation of 88 $\mu$mol/L, and without acute rejection episode in the first year post-transplantation. The recipient was returned to dialysis at 9.34 years after transplantation. The second patient (Figure 7) was a female, aged 60 years at transplantation, transplanted in 2007 for the second time, without history of a cardiovascular disease, immunized against HLA class I, and a SCr measurement at 3 months post-transplantation of 100 $\mu$mol/L, and at least one acute rejection episode in the first year after transplantation. The recipient was still alive with a functioning graft at 10.14 years post-transplantation. Figures 6 & 7 only present the results under the $N-N$ and $t-t$ models, whereas results under all the four models are presented in Figures 17 & 18 of the supplementary material. For these patients, dynamic predictions obtained under the models generally agree. There is a notable difference for the first patient in the fourth landmark time-point; see mid-right panel of Figure 6. The difference is due to the fourth SCr measurement that is a typical example of outlying observations: the fourth $SCr$ was 209, whereas the first three were 93, 110, and 115. $N-N$ and $t-N$ assumptions overreact and produce considerably lower survival probabilities compared to the $N-t$ and $t-t$. Median (2.5% and 97.5% percentiles of the MC samples) of the survival probabilities for 8.97 years after transplantation, i.e. $\mathbb{P}(T_{k}^{*}>8.97|T_{k}^{*}>3.97,Y_{k}=\{\log(93),\log(110),\log(115),\log(209)\},{\boldsymbol{\mathrm{x}}}_{k},{\boldsymbol{\mathrm{d}}}_{k},{\boldsymbol{\mathrm{c}}}_{k})$, were 0.35 (0.01, 0.78), 0.35 (0.00002, 0.82), 0.74 (0.27, 0.92), 0.77 (0.08, 0.92) under $N-N$, $t-N$, $N-t$ and $t-t$ models, respectively. Note that the patient had graft failure at year 9.3, i.e. $T_{k}^{*}=9.3$. The fifth $SCr$ measurement for the same subject is 195 that is closer in magnitude to the fourth measurement compared to the first three. $N-N$, $t-N$ and $t-t$ models correctly react to the fifth measurement and update the predictions accordingly. However, it $N-t$ model seems to produce somewhat higher median probabilities compared to these three models.

6 Discussion

In this work, we considered dynamic predictions of kidney graft survival using joint modelling of longitudinal and survival outcomes. We mainly focused on how the distributional assumptions would impact the predictions: widely used Normal distribution was compared against $t$. The proposed joint model with $t$ distributional assumptions is novel, as no work in the literature considered such a general model. Bayesian methods were considered for estimation and dynamic predictions. The proposed methods are implemented in the R package robjm. Methods are applied to data for kidney transplant patients from the French cohort DIVAT. Impacts of distributional assumptions on dynamic prediction performances were inspected through accuracy measures and predictions on two individuals are illustrated.

Regarding the DIVAT data-set, degree-of-freedom results indicated that there are both outlying individuals and outlying observations. We observed important differences between the standard and the robust joint models. The population averaged slope in the mixed-effects sub-model and the association parameter for the individual rate of change in the survival sub-model were largely over-estimated by the Gaussian joint model. In terms of dynamic predictions, we observed better calibration and discrimination for the robust joint models. Intermediate joint models, i.e. the $N-t$ and $t-N$ assumptions, suggest that outlying subjects have greater impact on dynamic predictions if they are not considered compared to outlying observations. Regarding individual dynamic predictions, robust models produced better results for the patient with outlying observation, whereas the two models produced similar results for the subject with no outlying observation.

Our current work voluntarily does not include a simulation study. Indeed, the objective is to present a case-study on predicting kidney graft failure risk, specifically for the patients in the DIVAT cohort. We clearly illustrated that considering longitudinal outliers have an impact on prognostic accuracy in the kidney transplantation context. Our modelling strategies may be applicable to others clinical contexts. When comparing the four joint models, we also voluntarily do not consider metrics for model selection. Indeed, since we are mainly interested in dynamic predictions, we compare the models using predictive accuracy measures obtained from a validation data-set.

In our robust joint model proposal, we considered symmetric $t$ distribution as an alternative to the Gaussian. Given that distributional assumptions might have considerable impacts on dynamic predictions, it would be worth to study other distributions than symmetric $t$, e.g. skew-$t$. We considered joint modelling framework with the shared-parameter formulation to link the longitudinal and survival sub-models. Performances of other methods, e.g. latent class joint modelling (Proust-Lima et al., 2014), Cox model with time-varying covariates, or landmarking methods (van Houwelingen, 2007), in the presence of longitudinal outliers would also be interesting to investigate. Diagnostics tools for checking the appropriateness of distributional assumptions for non-Gaussian joint models would be interesting. In this study, this was secondary to us, since our aim was to inspect if we can improve the dynamic predictions obtained from the widely used Normal assumption, and checked this using accuracy measures for the validation sample.

Fournier et al. (2019) previously presented dynamic predictions for the DIVAT patients using the $N-N$ joint model with a frequentist point of view. We are able to re-produce their results using our $N-N$ model. The AUC and $R^{2}$ values presented in (Fournier et al., 2019) are higher than those presented in the current paper. The reason for this is that they used median of the MC samples to obtain the point estimates and frequentist methods to obtain the associated intervals rather than obtaining MC samples of the accuracy measures. Yet, following their approach, we are able to obtain the values they presented.

In conclusion, this study presents improved dynamic predictions of kidney graft failure based on robust joint modelling framework. The most important component to obtain improved predictions seems the random-effects terms. Nonetheless, even letting only the error term being $t$-distributed improves predictions compared to the Gaussian model. We prefer working with the $t-t$ model, since i) it consists of the $N-t$ and $t-N$ models as special cases, ii) it is one of the best in terms of accuracy measures, iii) it provides a good compromise between the two special cases as presented for the first patient. The predictions based on the $t-t$ model will be deployed into the DynPG Shiny application of Fournier et al. (2019) (available at https://shiny.idbc.fr/DynPG) that is currently based on Gaussian joint modelling. This would allow physicians and patients get benefit from our methods.

Acknowledgements

Dr. Özgür Asar was funded by the French Embassy to Turkey for a visit to The University of Nantes, which enabled great progress for the current work. Helpful discussions with Dr. Jonas Wallin (Lund University), Prof. Peter Diggle (Lancaster University) and Dr. Yohann Foucher (Nantes University) are greatfully acknowledged. The analysis and interpretation of data collected from the French DIVAT (Données Informatisées et VAlidées en Transplantation) cohort* (www.divat.fr, No. CNIL 914184) are the responsibility of the authors. We wish to thank members of the data manager of the DIVAT cohort (Clarisse Kerleau) and clinical research assistant team (S. Le Floch, A. Petit, J. Posson, C. Scellier, V. Eschbach, K. Zurbonsen, C. Dagot, F. M’Raiagh, V. Godel, X. Longy, P. Przednowed). We are also grateful to Roche Pharma, Novartis and Sanofi laboratories for supporting the DIVAT cohort as the CENTAURE foundation (www.fondation-centaure.org).

• DIVAT cohort collaborators (Medical Doctors, Surgeons, HLA Biologists): Nantes: G. Blancho, J. Branchereau, D. Cantarovich, A. Chapelet, J. Dantal, C. Deltombe, L. Figueres, C. Garandeau, M. Giral, C. Gourraud-Vercel, M. Hourmant, G. Karam, C. Kerleau, A. Meurette, S. Ville, C. Kandell, A. Moreau, K. Renaudin, A. Cesbron, F. Delbos, A. Walencik, A. Devis; Paris-Necker: L. Amrouche, D. Anglicheau, O. Aubert, L. Bererhi, C. Legendre, A. Loupy, F. Martinez, R. Sberro-Soussan, A. Scemla, C. Tinel, J. Zuber; Nancy: P. Eschwege, L. Frimat, S. Girerd, J. Hubert, M. Ladriere, E. Laurain, L. Leblanc, P. Lecoanet, J-L. Lemelle; Lyon E. Hériot: L. Badet, M. Brunet, F. Buron, R. Cahen, S. Daoud, C. Fournie, A. Grégoire, A. Koenig, C. Lévi, E. Morelon, C. Pouteil-Noble, T. Rimmelé, O. Thaunat; Montpellier: S. Delmas, V. Garrigue, M. Le Quintrec, V. Pernin, J-E. Serre.