Assessing importance of biomarkers: A Bayesian joint modelling approach of longitudinal and survival data with semi-competing risks

Abstract

Longitudinal biomarkers such as patient-reported outcomes (PROs) and quality of life (QOL) are routinely collected in cancer clinical trials or other studies. Joint modelling of PRO/QOL and survival data can provide a comparative assessment of patient-reported changes in specific symptoms or global measures that correspond to changes in survival. Motivated by a head and neck cancer clinical trial, we develop a class of trajectory-based models for longitudinal and survival data with disease progression. Specifically, we propose a class of mixed effects regression models for longitudinal measures, a cure rate model for the disease progression time ( $T_P$ ) and a Cox proportional hazards model with time-varying covariates for the overall survival time ( $T_D$ ) to account for $T_P$ and treatment switching. Under the semi-competing risks framework, the disease progression is the non-terminal event, the occurrence of which is subject to a terminal event of death. The properties of the proposed models are examined in detail. Within the Bayesian paradigm, we derive the decompositions of the deviance information criterion (DIC) and the logarithm of the pseudo-marginal likelihood (LPML) to assess the fit of the longitudinal component of the model and the fit of each survival component, separately. We further develop $\Delta$ DIC as well as $\Delta$ LPML to determine the importance and contribution of the longitudinal data to the model fit of the $T_P$ and $T_D$ data.

1 Introduction

A patient-reported outcome (PRO) is any report of the status of a patient's health condition that comes directly from the patient, without interpretation of the patient's response by a clinician or anyone else (U.S. Food and Drug Administration Guidance for Industry, 2009). PROs can help identify symptoms or problems that may be missed during clinical queries, guide policy and health care delivery, and help patients choose between two equally efficacious therapies. Joint modelling of PRO and survival data can provide a comparative assessment of patient-reported changes in specific symptoms or global measures that correspond to changes in survival.

In joint modelling literature, there are usually a longitudinal component, a survival component and a dependence structure connecting those two components. There has been much work on joint modelling of one longitudinal measure and a survival measure, including Schluchter (1992), Hogan and Laird (1997), Law et al. (2002), Brown and Ibrahim (2003), Chen et al. (2004), Ibrahim et al. (2004), Brown et al. (2005), Chi and Ibrahim (2006), Chi and Ibrahim (2007), Ibrahim et al. (2010), Ye et al. (2008), Brown (2009), Rizopoulos (2012)), Pawitan and Self (1993), De Gruttola and Tu (1994), Lavalley and Degruttola (1996), De Gruttola and Tu (1994), Xu and Zeger (2001b), Xu and Zeger (2001a), Jacqmin-Gadda et al. (2010)) and (Proust-Lima et al. (2014)). (Tudur-Smith et al. (2016)) provided a recent review on joint modelling. In the recent development, multiple longitudinal processes and spline-based methods have been used in the joint modelling (Rizopoulos and Ghosh (2011)). Joint modelling has also been applied to more complex data structures such as competing risks data (Williamson et al., 2008; Elashoff et al., 2008; Huang et al., 2011; Li et al., 2012; Proust-Lima et al., 2016) and semi-competing risks data (Li and Su, 2017). In this article, we propose a mixed effects regression model for longitudinal measures, a cure rate model for the disease progression time ( $T_P$ ) and a Cox proportional hazards model with time-varying covariates for the overall survival time ( $T_D$ ) to account for $T_P$ and treatment switching. Note that Chen et al. (2020) provide a detailed elaboration of treatment switching, and they also develop a method for estimating the intended treatment effect accounting for treatment switching. In addition, the trajectory function from the longitudinal model is substituted into the hazard function of each survival submodel, thereby serving as another time-varying covariate in the survival submodel. Thus, our model is the trajectory model (TM). Although Li and Su (2017) address a very similar context as semi-competing risks within the joint modelling framework, they assume a probit model for the discrete time hazard of dropout and another probit model for the discrete time hazard of death. Li and Su (2017) further assume the shared parameter model (SPM), where the longitudinal model and the survival submodels share common random effects, which then induce correlation between the longitudinal and survival components.

Zhang et al. (2017)) developed a novel decomposition of DIC and LPML to assess the fit of the longitudinal and survival components of the joint model, separately. Based on this decomposition, they then proposed new Bayesian model assessment criteria, namely $\Delta$ DIC and $\Delta$ LPML, to determine the importance and contribution of the longitudinal data to the model fit of the survival data. In this article, we extend the decomposition of Zhang et al. (2017) to our joint model with semi-competing risks. We derive the decompositions of DIC and LPML (i.e., DIC $=$ DIC ${}_{\mathit{Long}}$ $+$ DIC ${}_{\mathit{Pg}\mid \mathit{Long}}$ $+$ DIC ${}_{\mathit{Surv}\mid \mathit{Pg},\mathit{Long}}$ and LPML $=$ LPML ${}_{\mathit{Long}}$ $+$ LPML ${}_{\mathit{Pg}\mid \mathit{Long}}$ $+$ LPML ${}_{\mathit{Surv}\mid \mathit{Pg},\mathit{Long}}$ , respectively) to assess the fit of the longitudinal component of the model and the fit of each survival component, separately. Let $\Delta {\mathit{DIC}}_{\mathit{Surv}}$ be the difference of ${\mathit{DIC}}_{\mathit{Surv}|\mathit{Pg},\mathit{Long}}$ and DIC of the overall survival model alone without joint modelling. $\Delta {\mathit{DIC}}_{\mathit{Surv}}$ measures the additional information gained by adding disease progression time and longitudinal components into modelling of overall survival data. Similarly, $\Delta {\mathit{DIC}}_{\mathit{Pg}}$ measures the additional information gained by adding longitudinal components into modelling of progression data. Similar ideas apply to $\Delta {\mathit{LPML}}_{\mathit{Surv}}$ and $\Delta {\mathit{LPML}}_{\mathit{Pg}}$ . We also develop a simple but more attractive prior for the covariance matrix based on its Cholesky decomposition borrowing the idea of the Lewandowski, Kurowicka, and Joe (LKJ) prior for the correlation matrix $C$ which has a density function of $f(C|\tau )⋑|C{|}^{\tau -1}$ (Lewandowski et al., 2009).

The rest of the article is organized as follows. Section 2 gives a detailed description of the motivating longitudinal and survival data from a head and neck clinical cancer trial. In Section 3, we propose a mixed effects regression model for the longitudinal measure, a cure rate model for $T_P$ and a time-varying covariates model for $T_D$ . The properties of the proposed models are discussed. The likelihood functions under the proposed models are derived in Section 4. In Section 4, we also present the prior and posterior distributions, develop an efficient Markov chain Monte Carlo sampling algorithm for carrying out the posterior computations, and provide the formulation of the decomposition of the LPML and $\Delta$ LPML's for model comparison. A comprehensive analysis of the longitudinal and survival data from a head and neck cancer trial is given in Section 5. We conclude the article with a discussion including some extensions in Section 6.

2 Motivating longitudinal and survival data from a head and neck cancer clinical trial

The motivating study was an open-label, randomized, active-controlled phase III trial comparing the efficacy and safety of afatinib versus methotrexate in patients with recurrent/metastatic squamous carcinoma cancer of head and neck (HNSCC Machiels et al., 2015). Patients were treated on either afatinib or methotrexate, according to randomization, until the occurrence of disease progression or intolerable adverse events. Upon discontinuation of randomized treatment, a patient was free to take other anti-cancer therapies per protocol (except that patients discontinued methotrexate cannot start afatinib). It was reported that more than half of the patients who discontinued study medication started another anti-cancer treatment.

Figure 1:

Graphical representation of the disease progression and survival data from the head and neck cancer clinical trial

Of the total 483 patients, 332 patients were randomized to the treatment arm (methotrexate) and 161 patients were randomized to the control arm (afatinib). Furthermore, 246 patients switched the treatments, 358 died and 125 were censored. The median follow-up of 6.7 months (interquartile range (IQR): 3.1 –9.0). Figure 1 shows the counts of disease progression, treatment switch and death in each treatment arm. Three patient-reported measures: global health status (GH), head and neck pain scale (HNPA), and head and neck swallowing condition (HNSW) were collected during this phase III trial. The number of patients with at least one post-baseline values of the longitudinal measures were 481 for GH, 480 for HNPS and 473 for HNSW as shown in Table 1. In this article, we investigate the associations between changes in the longitudinal measures and risks of disease progression ( $T_P$ ) and death ( $T_D$ ), in which $T_D$ was the key secondary endpoint of the trial.

Table 1:

A summary of the observed counts of longitudinal measures GH, HNPA and HNSW

Number of observed longitudinal measures	Count (percentage)
	GH	HNPA	HNSW
0	6 (1.24%)	7 (1.45%)	14 (2.90%)
1	78 (16.15%)	77 (15.94%)	82 (16.98%)
2	86 (17.81%)	87 (18.01%)	94 (19.46%)
3	111 (22.98%)	113 (23.40%)	104 (21.53%)
4	100 (20.70%)	96 (19.88%)	91 (18.84%)
5	52 (10.77%)	52 (10.77%)	50 (10.35%)
$\ge 6$	50 (10.35%)	51 (10.56%)	48 (9.94%)

We consider the following baseline covariates: patients’ prior treatment with an epidermal growth factor receptor (EGFR)’targeted antibody, cetuximab (PEGFRA: 1 if and 0 if no), age more than 65 years (AGE: 1 if yes and 0 if no), gender (SEX: male/female), prior treatment with plantinum-based chemotherapy and radiotherapy (PCRTDC: 1 if yes and 0 if no), primary tumour site (TMSITED01: 1 if oral cavity, oropharynx, hypopharynx, larynx and 0 if others), metastatic sites of lymph node (METLOC2: 1 if yes and 0 if no), metastatic sites of lung (METLOC3: 1 if yes and 0 if no), liver (METLOC4: 1 if yes and 0 if no), or bone (METLOC5: 1 if yes and 0 if no), Eastern Cooperative Oncology Group (ECOG) performance status at baseline (ECOG: 0/1) and BSDTARG (Baseline target lesion: mm). We exclude two patients from the analysis because of missing BSDTARG measures.

Median survival time in days and its 95% confidence interval are 434 and (386, 470), median switch time, and its 95% confidence interval are 166 and (147, 190), and median progression time and its 95% confidence interval are 84 and (81, 86).

3 The proposed models

Let $Y(t)$ denote the longitudinal measure observed at time $t$ . Also, let $T_P$ denote the time from the study entry until disease progression and let $T_D$ denote the time from the study entry until death. In addition, we let $\mathit{x}$ be the $p$ -dimensional vector of baseline covariates including the randomized treatment but not the intercept.

3.1 Longitudinal component of the joint model

We first assume a mixed effects regression model for the longitudinal measure $Y(t)$ , which is given by

$$Y(t)=({\mathbf{\theta }}_0^{*}+\mathbf{\theta }{)}^{\prime }\mathit{g}(t)+{\mathbf{\gamma }}^{\prime }\mathit{x}+\epsilon (t),$$

(3.1)

where $\mathit{g}(t)$ is a $(q+1)$ -dimensional vector of functions of $t$ including the intercept, ${\mathbf{\theta }}_0^{*}$ is a $(q+1)$ -dimensional vector of subject-dependent random effects, $\mathbf{\theta }$ is the $(q+1)$ -dimensional vector of the overall effects, $\mathbf{\gamma }$ is a $p$ -dimensional vector of regression coefficients, and $\epsilon (t)$ is a measurement error. We further assume that ${\mathbf{\theta }}_0^{*}=({\theta }_{00}^{*},\dots ,{\theta }_{0q}^{*}{)}^{\prime }$ and $\epsilon (t)$ for $t\ge 0$ are independent;

$${\mathbf{\theta }}_0^{*}\sim N(0,\Omega ),$$

(3.2)

where $\Omega$ is a $(q+1)\times (q+1)$ positive definite covariance matrix; and for $t\ge 0$ ,

$$\epsilon (t)\sim N(0,{\sigma }^2),$$

(3.3)

where ${\sigma }^2>0$ is a variance parameter.

We note that in (3.1), if $q=1$ , $\mathit{g}(t)=(1,t{)}^{\prime }$ and $({\mathbf{\theta }}_0^{*}+\mathbf{\theta }{)}^{\prime }\mathit{g}(t)$ represents a linear trajectory, and if $q=2$ , $\mathit{g}(t)=(1,t,t^2{)}^{\prime }$ and $({\mathbf{\theta }}_0^{*}+\mathbf{\theta }{)}^{\prime }\mathit{g}(t)$ leads to a quadratic trajectory. A more general form of $\mathit{g}(t)$ is given by $\mathit{g}(t)=(1,b_1(t),\dots ,b_q(t){)}^{\prime }$ , where $b_1(t)$ , $\dots$ , $b_q(t)$ are a set of base functions. We also note that in (3.1), ${\mathbf{\theta }}_0^{*}$ captures the patient-level dependence of the $Y(t)$ over time $t$ . However, the covariance matrix $\Omega$ depends on the data only via random effects ${\mathbf{\theta }}_0^{*}$ , which can lead to slow convergence of Gibbs sampling. To overcome this issue, we consider the reparametrization: ${\mathbf{\theta }}_0^{*}=\Gamma {\mathbf{\theta }}^{*}$ , where $\Gamma$ is the Cholesky decomposition of $\Omega$ with all diagonal elements being positive such that $\Omega =\Gamma {\Gamma }^{\prime }$ . Under this reparametrization, we have

$${\mathbf{\theta }}^{*}\sim N(0,I_{q+1})$$

(3.4)

and a reparametrized version of (3.1) is written as

$$Y(t)=(\Gamma {\mathbf{\theta }}^{*}+\mathbf{\theta }{)}^{\prime }\mathit{g}(t)+{\mathbf{\gamma }}^{\prime }\mathit{x}+\epsilon (t).$$

(3.5)

Remark 3.1: Write

$$\Gamma =\left[\begin{array}{ccccc}\text{exp}(b_1)& 0& 0& \dots & 0\\ a_{21}& \text{exp}(b_2)& 0& \dots & 0\\ 5\\ a_{q+1,1}& a_{q+1,2}& a_{q+1,3}& \dots & \text{exp}(b_{q+1})\end{array}\right].$$

Then for any values of $b_1,\dots ,b_{q+1}$ and $a_{\mathit{ij}}$ , $i=2,\dots ,q+1$ , $j=2,\dots ,i-1$ , $\Gamma$ is a full-rank matrix. Thus, there are no constraints on $b_1,\dots ,b_{q+1}$ and $a_{\mathit{ij}}$ , $i=2,\dots ,q+1$ , $j=2,\dots ,i-1$ . Consequently, a multivariate normal prior can be specified for these parameters. Therefore, in this sense, our parametrization approach may be more attractive than the one discussed in Lewandowski et al. (2009). In addition, we have $\mathit{g}(t{)}^{\prime }\Gamma {\mathbf{\theta }}^{*}={\sum }_{j=1}^{q+1}{\sum }_{i=j}^{q+1}{g}_i(t)a_{\mathit{ij}}{\theta }_j^{*}$ , where $a_{\mathit{jj}}=\text{exp}(b_j)$ for $j=1,2,\dots ,q+1$ . We may view ${\theta }_j^{*}$ ’s as ‘known data points’ and $a_{\mathit{ij}}$ ’s as the ‘regression coefficients’. Under some mild conditions of the longitudinal data, the parameters $a_{\mathit{ij}}$ ’s are stochastically identifiable.

3.2 Survival component of the joint model

Let $\delta$ denote the disease progression status of subjects such that $\delta =1$ if the subject has disease progression and $\delta =0$ if the subject never has disease progression. We assume a logistic regression model for $\delta$ as

$$\mathit{logit}(P(\delta =1|\mathbf{\varphi },\mathit{x}))=\text{log}\{\frac{P(\delta =1|\mathbf{\varphi },\mathit{x})}{1-P(\delta =1|\mathbf{\varphi },\mathit{x})}\}={\varphi }_0+{\mathbf{\varphi }}_1^{\prime }\mathit{x},$$

(3.6)

where $\mathbf{\varphi }=({\varphi }_0,{\mathbf{\varphi }}_1^{\prime }{)}^{\prime }$ is a ( $p+1$ )’dimensional vector of regression coefficients. Following Berkson and Gage (1952) and Chen et al. (1999)), we assume a cure rate model for $T_p$ with

$$P(T_P=\infty |\mathit{x},{\mathbf{\theta }}^{*},\delta =0)=1$$

(3.7)

and a Cox-type model when $\delta =1$ with

$$\begin{array}{cc}\hfill & P(T_P>t|{\lambda }_{10},{\mathbf{\beta }}_1,{\mathbf{\alpha }}_1,\mathit{x},{\mathbf{\theta }}^{*},\mathbf{\theta },\delta =1)\hfill \\ =\hfill & \text{exp}[-{\int }_0^t{\lambda }_{10}(u)\text{exp}\{{\beta }_1(\Gamma {\mathbf{\theta }}^{*}+\mathbf{\theta }{)}^{\prime }\mathit{g}(u)+{\mathbf{\alpha }}_1^{\prime }\mathit{x}\}\mathit{du}],\hfill \end{array}$$

(3.8)

where ${\lambda }_{10}(t)$ is a proper baseline hazard function such that ${\int }_0^{\infty }{\lambda }_{10}(t)\mathit{dt}=\infty$ , ${\beta }_1$ is a regression coefficient corresponding to the longitudinal trajectory function $(\Gamma {\mathbf{\theta }}^{*}+\mathbf{\theta }{)}^{\prime }\mathit{g}(t)$ and ${\mathbf{\alpha }}_1$ is a $p$ -dimensional vector of regression coefficients corresponding to the baseline covariates. From (3.7) and (3.8), we see that $T_P$ can be written as

$$T_P=\infty \times 1\{\delta =0\}+T_P^{*}\times 1\{\delta =1\},$$

(3.9)

where $\infty \times 0=0$ , $T_P^{*}$ is a proper random survival time that follows the distribution specified in (3.8) and $1\left\{A\right\}$ denotes the indicator function, which takes a value of 1 if $A$ is true and 0 otherwise.

We assume a time-varying covariates Cox model for $T_D$ with hazard function defined as

$$\begin{array}{cc}{\lambda }_2(t|{\lambda }_{20},{\beta }_2,{\mathbf{\theta }}^{*},\mathbf{\theta },{\mathbf{\alpha }}_2,\mathbf{\eta },\mathit{x},\mathit{z}(t))=\hfill & {\lambda }_{20}(t)\text{exp}\{{\beta }_2(\Gamma {\mathbf{\theta }}^{*}+\mathbf{\theta }{)}^{\prime }\mathit{g}(t)+{\mathbf{\alpha }}_2^{\prime }\mathit{x}+{\mathbf{\eta }}^{\prime }\mathit{z}(t)\},\hfill \end{array}$$

(3.10)

where ${\beta }_2$ is a regression coefficient corresponding to the longitudinal trajectory function $(\Gamma {\mathbf{\theta }}^{*}+\mathbf{\theta }{)}^{\prime }\mathit{g}(t)$ , ${\mathbf{\alpha }}_2$ is a $p$ -dimensional vector of regression coefficients and $\mathit{z}(t)$ is the history of time-dependent covariate up to $t$ with $\mathbf{\eta }$ being a vector of the corresponding regression coefficients.

Remark 3.2: We note that although we use $\mathit{x}$ to denote a vector of baseline covariates, different subsets of baseline covariates may be used in (3.1, 3.6, 3.8) and (3.10), respectively. Also $\mathit{x}$ can be different in the longitudinal and survival models. For the notational simplicity, we use the same $\mathit{x}$ for baseline covariates throughout the article.

Remark 3.3: Let $V_{T_P}(t)=1\{T_P\le t\}$ . The variable $V_{T_P}(t)$ essentially tracks the status of disease progression so that $V_{T_P}(t)=0$ if the patient has not had disease progression at time $t$ and $V_{T_P}(t)=1$ if the patient has had disease progression before time $t$ . Similarly, we denote $V_{T_{\mathit{SW}}}(t)=1$ if the subject has switched treatment at time $t$ , and $V_{T_{\mathit{SW}}}(t)=0$ if the subject has switched treatment at time $t$ . Suppose that $R$ denotes the indicator for the randomized treatment such that $R=1$ if the subject is on the treatment arm and $R=0$ if the subject is on the control arm. In (3.10), we take $\mathit{z}(t)=(V_{T_P}(t),V_{T_{\mathit{SW}}}(t),R\times V_{T_{\mathit{SW}}}(t))$ .

Remark 3.4: From the proposed joint model for the longitudinal data and survival data, we see that ${\mathbf{\theta }}^{*}$ captures the dependence between the $Y(t)$ and $(T_P,T_D)$ , while the joint distribution of $T_P$ and $T_D$ is specified by the product of the marginal distribution $[T_P]$ and the conditional distribution $[T_D|T_P]$ . Since the event ‘ $\delta =0$ ’ is not observed, $\mathbf{\varphi }$ is only partially identifiable as shown in Chen et al. (1999). Thus, in our analysis of the head and neck cancer data and the simulation study, we specify a moderate informative prior for $\mathbf{\varphi }$ .

Remark 3.5: The regression coefficient ${\beta }_1$ in (3.8) ( ${\beta }_2$ in (3.10)) captures the association between the longitudinal measure $Y(t)$ and $T_P$ ( $T_D$ ). If ${\beta }_1=0$ and ${\beta }_2=0$ , $Y(t)$ is independent of $T_P$ and $T_D$ . If ${\beta }_1>0$ , an increase in the longitudinal measure $Y(t)$ is associated with a shorter disease programme time $T_P$ . Similarly, if ${\beta }_2>0$ , an increase in $Y(t)$ is associated with a worse overall survival mortality. In addition, the longitudinal trajectory function $(\Gamma {\mathbf{\theta }}^{*}+\mathbf{\theta }{)}^{\prime }\mathit{g}(t)$ is directly included in the hazard function for $T_P$ in (3.8) and the hazard function for $T_D$ in (3.10) as a time-dependent covariate. Therefore, the proposed joint model is the TM.

3.3 The likelihood function

Suppose that there are $n$ subjects. For the $i\mathrm{th}$ subject, let $y_i(t)$ denote the longitudinal measure, which is observed at time $t\in \{a_{i1},a_{i2},\dots ,a_{i{m}_i}\}$ , where $0\le a_{i1}<a_{i2}<\dots <a_{{\mathit{im}}_i}$ and $m_i\ge 0$ . Note that $y_i(0)$ corresponds to the baseline value. Also let ${\mathit{x}}_i$ denote a $p$ -dimensional vector of covariates, which may include the randomized treatment. Let $t_i$ and ${\nu }_i$ denote the overall survival time and the right-censored indicator, respectively, where ${\nu }_i=1$ if $t_i$ is a death time and 0 if $t_i$ is right-censored for the $i\mathrm{th}$ subject. Also let $t_{\mathit{Pi}}$ be the disease progression time and ${\delta }_i$ be the status of disease progression. If disease progression has been observed, then $t_{\mathit{Pi}}<t_i$ and ${\delta }_i=1$ . However, if disease progression has not yet been observed at time $t_i$ , $t_{\mathit{Pi}}$ and ${\delta }_i$ are not observed. Let $t_{\mathit{SWi}}$ denote treatment switch time, and if treatment switch has not been observed yet at time $t_i$ , $t_{\mathit{SWi}}$ is not observed. Let ${\mathit{z}}_i(t)=(V_{t_{\mathit{Pi}}}(t),V_{t_{\mathit{SWi}}}(t),R\times V_{t_{\mathit{SWi}}}(t))$ , where $t\le t_i$ .

Write ${\mathit{y}}_i=(y_i(a_{i1}),\dots ,y_i(a_{{\mathit{im}}_i}){)}^{\prime }$ and $W_i=(W_{i1},W_{i2})$ , which is an $m_i\times (q+1+p)$ matrix with $m_i\times (q+1)$ matrix $W_{i1}=(\mathit{g}(a_{\mathit{ij}}{)}^{\prime },j=1,\dots ,m_i)$ and $m_i\times p$ matrix $W_{i2}=({\mathit{x}}_i^{\prime },j=1,\dots ,m_i)$ . Then, the complete data likelihood function for the longitudinal data ${\mathit{y}}_i$ is given by

$$\begin{array}{cc}f({\mathit{y}}_i|\mathbf{\gamma },{\sigma }^2,\Gamma ,{\mathbf{\theta }}_i^{*},\mathbf{\theta },{\mathit{x}}_i)=\hfill & \frac{1}{(2\pi {\sigma }^2{)}^{\frac{m_i}{2}}}\text{exp}\{-\frac{1}{2{\sigma }^2}({\mathit{y}}_i-W_{i1}^{*}{\mathbf{\theta }}_i^{*}-W_i({\mathbf{\theta }}^{\prime },{\mathbf{\gamma }}^{\prime }{)}^{\prime }{)}^{\prime }\hfill \\ \hfill & ({\mathit{y}}_i-W_{i1}^{*}{\mathbf{\theta }}_i^{*}-W_i({\mathbf{\theta }}^{\prime },{\mathbf{\gamma }}^{\prime }{)}^{\prime })\},\hfill \end{array}$$

(3.11)

where ${\mathbf{\theta }}_i^{*}$ is a $(q+1)$ -dimensional vector of random effects, $W_{i1}^{*}=W_{i1}\Gamma$ and, $W_i({\mathbf{\theta }}^{\prime },{\mathbf{\gamma }}^{\prime }{)}^{\prime }=W_{i1}\mathbf{\theta }+W_{i2}\mathbf{\gamma }$ .

Using (3.2), the density of ${\mathbf{\theta }}_i^{*}$ is given by

$$f({\mathbf{\theta }}_i^{*})=\frac{1}{(2\pi {)}^{\frac{q+1}{2}}}\text{exp}\{-\frac{1}{2}({\mathbf{\theta }}_i^{*}{)}^{\prime }({\mathbf{\theta }}_i^{*})\},$$

(3.12)

for $i=1,\dots ,n$ .

Based on the definition of $\mathit{z}(t)$ , we see that $\mathit{z}(t)$ depends on $t_{\mathit{Pi}}$ and $t$ only through $V_{t_{\mathit{Pi}}}(t)=1\{t_{\mathit{Pi}}\le t\}$ . Write $p(\mathbf{\varphi },{\mathit{x}}_i)=\frac{\text{exp}({\varphi }_0+{\mathbf{\varphi }}_1^{\prime }{\mathit{x}}_i)}{1+\text{exp}({\varphi }_0+{\mathbf{\varphi }}_1^{\prime }{\mathit{x}}_i)}$ . We also let ${\zeta }_i$ denote the group indicator.

Define the survival function for $T_D$ $P(T_D>t_i|{\mathbf{\varphi }}_3,{\mathit{x}}_i,{\mathit{z}}_i(t),\Gamma ,{\mathbf{\theta }}_i^{*},\mathbf{\theta })=\text{exp}[-{\int }_0^{t_i}{\lambda }_{20}(u)$ $\text{exp}\{{\beta }_2((\Gamma {\mathbf{\theta }}_i^{*}+\mathbf{\theta }{)}^{\prime }\mathit{g}(u))+{\mathbf{\alpha }}_2^{\prime }{\mathit{x}}_i+{\mathbf{\eta }}^{\prime }{\mathit{z}}_i(u)\}\mathit{du}]$ , and the density function $f_2(t_i|{\mathbf{\varphi }}_3,{\mathit{x}}_i,{\mathit{z}}_i(t),\Gamma ,{\mathbf{\theta }}_i^{*},\mathbf{\theta })={\lambda }_{20}(t_i)\text{exp}\{{\beta }_2((\Gamma {\mathbf{\theta }}_i^{*}+\mathbf{\theta }{)}^{\prime }\mathit{g}(t_i))+{\mathbf{\alpha }}_2^{\prime }{\mathit{x}}_i+{\mathbf{\eta }}^{\prime }{\mathit{z}}_i(t_i)\}P(T_D>t_i|{\mathbf{\varphi }}_3,{\mathit{x}}_i,{\mathit{z}}_i(t),\Gamma ,{\mathbf{\theta }}_i^{*},\mathbf{\theta })$ . Define the survival function for $T_P$ $P(T_P>t_i|{\mathbf{\varphi }}_{21},{\mathit{x}}_i,\Gamma ,{\mathbf{\theta }}_i^{*},\mathbf{\theta },\delta =1)=\text{exp}[-{\int }_0^{t_i}{\lambda }_{10}(u)\text{exp}\{{\beta }_1((\Gamma {\mathbf{\theta }}_i^{*}+\mathbf{\theta }{)}^{\prime }\mathit{g}(u))+{\mathbf{\alpha }}_1^{\prime }{\mathit{x}}_i\}\mathit{du}]$ , and the density function $f_1(t_{\mathit{Pi}}|{\mathbf{\varphi }}_{21},{\mathit{x}}_i,\Gamma ,{\mathbf{\theta }}_i^{*},\mathbf{\theta },\delta =1)={\lambda }_{10}(t_{\mathit{Pi}})\text{exp}\{{\beta }_1((\Gamma {\mathbf{\theta }}_i^{*}+\mathbf{\theta }{)}^{\prime }\mathit{g}(t_{\mathit{Pi}}))+{\mathbf{\alpha }}_1^{\prime }{\mathit{x}}_i\}P(T_P>t_{\mathit{Pi}}|{\mathbf{\varphi }}_{21},{\mathit{x}}_i,\Gamma ,{\mathbf{\theta }}_i^{*},\mathbf{\theta },\delta =1)$ .

Under the proposed cure rate and time-varying covariates model, the observed data can be divided into four groups of observations. We derive the likelihood function for each of the four groups as follows.

Group 1. Subjects are observed to die at time $t_i$ and no disease progressions have been incurred. The survival data consists of $D_{i1}=\{t_i,{\nu }_i=1,{\mathit{x}}_i,{\mathit{z}}_i(t)\}$ , group indicator ${\zeta }_i=1$ and ${\delta }_i$ is not observed. Since $t_{\mathit{Pi}}=\infty$ when ${\delta }_i=0$ and $t_{\mathit{Pi}}\ge t_i$ when ${\delta }_i=1$ . The likelihood function is given by

$$\begin{array}{cc}L(\mathbf{\varphi },{\lambda }_{10},\hfill & {\mathbf{\beta }}_1,{\mathbf{\alpha }}_1,{\lambda }_{20},{\beta }_2,{\mathbf{\alpha }}_2,\mathbf{\eta }|D_{i1},{\mathbf{\theta }}_i^{*})=f_2(t_i|{\mathbf{\varphi }}_3,{\mathit{x}}_i,{\mathit{z}}_i(t),\Gamma ,{\mathbf{\theta }}_i^{*},\mathbf{\theta })\hfill \\ \hfill & \times (1-p(\mathbf{\varphi },{\mathit{x}}_i)+p(\mathbf{\varphi },{\mathit{x}}_i)P(T_P>t_i|{\mathbf{\varphi }}_{21},{\mathit{x}}_i,\Gamma ,{\mathbf{\theta }}_i^{*},\mathbf{\theta },\delta =1)).\hfill \end{array}$$

(3.13)

Group 2. Subjects are observed to have disease progression at time $t_{\mathit{Pi}}$ and then die at $t_i$ . The survival data consists of $D_{i2}=\{t_i,{\nu }_i=1,t_{\mathit{Pi}},{\delta }_i=1,{\mathit{x}}_i,{\mathit{z}}_i(t)\}$ , $t_i>t_{\mathit{Pi}}$ , and ${\zeta }_i=2$ . The likelihood function is given by

$$\begin{array}{cc}L(\mathbf{\varphi },{\lambda }_{10},\hfill & {\mathbf{\beta }}_1,{\mathbf{\alpha }}_1,{\lambda }_{20},{\mathbf{\beta }}_2,{\mathbf{\alpha }}_2,\mathbf{\eta }|D_{i2},{\mathbf{\theta }}_i^{*})=f_2(t_i|{\mathbf{\varphi }}_3,{\mathit{x}}_i,{\mathit{z}}_i(t),\Gamma ,{\mathbf{\theta }}_i^{*},\mathbf{\theta })p(\mathbf{\varphi },{\mathit{x}}_i)\hfill \\ \hfill & \times f_1(t_{\mathit{Pi}}|{\mathbf{\varphi }}_{21},{\mathit{x}}_i,\Gamma ,{\mathbf{\theta }}_i^{*},\mathbf{\theta },\delta =1)\hfill \end{array}$$

(3.14)

Group 3. Subjects are observed to have disease progression at time $t_{\mathit{Pi}}$ and then censor at $t_i$ . The survival data consists of $D_{i3}=\{t_i,{\nu }_i=0,t_{\mathit{Pi}},{\delta }_i=1,{\mathit{x}}_i,{\mathit{z}}_i(t)\}$ , $t_i>t_{\mathit{Pi}}$ , and ${\zeta }_i=3$ . The likelihood function takes the form

$$\begin{array}{cc}L(\mathbf{\varphi },{\lambda }_{10},\hfill & {\mathbf{\beta }}_1,{\mathbf{\alpha }}_1,{\lambda }_{20},{\mathbf{\beta }}_2,{\mathbf{\alpha }}_2,\mathbf{\eta }|D_{i3},{\mathbf{\theta }}_i^{*})=P(T_D>t_i|{\mathbf{\varphi }}_3,{\mathit{x}}_i,{\mathit{z}}_i(t),\Gamma ,{\mathbf{\theta }}_i^{*},\mathbf{\theta })\hfill \\ \hfill & \times p(\mathbf{\varphi },{\mathit{x}}_i)f_1(t_{\mathit{Pi}}|{\mathbf{\varphi }}_{21},{\mathit{x}}_i,\Gamma ,{\mathbf{\theta }}_i^{*},\mathbf{\theta },\delta =1).\hfill \end{array}$$

(3.15)

Group 4. Subjects are censored at $t_i$ and no disease progression occurs before $t_i$ . The survival data consists of $D_{i4}=\{t_i,{\nu }_i=0,{\mathit{x}}_i,{\mathit{z}}_i(t)\}$ , ${\zeta }_i=4$ and ${\delta }_i$ is not observed. The likelihood function is written as

$$\begin{array}{cc}L(\mathbf{\varphi },{\lambda }_{10},\hfill & {\mathbf{\beta }}_1,{\mathbf{\alpha }}_1,{\lambda }_{20},{\mathbf{\beta }}_2,{\mathbf{\alpha }}_2,\mathbf{\eta }|D_{i4},{\mathbf{\theta }}_i^{*})=P(T_D>t_i|{\mathbf{\varphi }}_3,{\mathit{x}}_i,{\mathit{z}}_i(t),\Gamma ,{\mathbf{\theta }}_i^{*},\mathbf{\theta })\hfill \\ \hfill & \times (1-p(\mathbf{\varphi },{\mathit{x}}_i)+p(\mathbf{\varphi },{\mathit{x}}_i)P(T_P>t_i|{\mathbf{\varphi }}_{21},{\mathit{x}}_i,\Gamma ,{\mathbf{\theta }}_i^{*},\mathbf{\theta },\delta =1)).\hfill \end{array}$$

(3.16)

Remark 3.6: The disease progression is a non-terminal event while the death is a terminal event. The overall survival time $T_D$ censors the disease progression time $T_P$ , but not vice versa. Thus, the disease progress and the death constitute the semicompeting risks. Under the proposed model, these semi-competing risks are factored into the constructions of the likelihood functions in Groups 1 and 4.

4 Bayesian inference

4.1 Prior and posterior distributions

In (3.8) and (3.10), we assume a piecewise constant hazard function for ${\lambda }_{j0}(t)$ defined as

$${\lambda }_{j0}(t)={\lambda }_{\mathit{jk}}, t\in (s_{j,k-1},s_{\mathit{jk}}] \mathit{for}k=1,\dots ,K_j,$$

(4.1)

where $0=s_{j0}<s_{j1}<s_{j2}<\dots <s_{j,K_j-1}<s_{{\mathit{jK}}_j}=\infty$ is a finite partition of the time axis, where $K_j\ge 1$ for $j=1,2$ . Write ${\mathbf{\lambda }}_j=({\lambda }_{j1},\dots ,{\lambda }_{{\mathit{jK}}_j}{)}^{\prime }$ for $j=1,2$ . We assume independent normal priors for $\mathbf{\gamma }$ , $\mathbf{\theta }$ , $\mathbf{\varphi }$ , ${\mathbf{\beta }}_1$ , ${\mathbf{\alpha }}_1$ , ${\mathbf{\beta }}_2$ , ${\mathbf{\alpha }}_2$ and $\mathbf{\eta }$ with $\mathbf{\gamma }\sim N({\mathbf{\mu }}_{01},V_{01})$ , $\mathbf{\theta }\sim N({\mathbf{\mu }}_{02},V_{02})$ $\mathbf{\varphi }\sim N({\mathbf{\mu }}_{03},V_{03})$ $({\mathbf{\beta }}_1^{\prime },{\mathbf{\alpha }}_1^{\prime }{)}^{\prime }\sim N({\mathbf{\mu }}_{04},V_{04})$ , and $({\mathbf{\beta }}_2^{\prime },{\mathbf{\alpha }}_2^{\prime },{\mathbf{\eta }}^{\prime }{)}^{\prime }\sim N({\mathbf{\mu }}_{05},V_{05})$ , where ${\mathbf{\mu }}_{0j}$ and $V_{0j}$ (a positive definite matrix) for $j=1,\dots ,5$ are pre-specified hyperparameters. An inverse gamma (IG) prior is specified for ${\sigma }^2$ , that is, ${\sigma }^2\sim \mathit{IG}(a_0,b_0)$ with mean $b_0/(a_0-1)$ , where $a_0\ge 0$ and $b_0\ge 0$ are pre-specified. When $q=1$ , we write $\Gamma =\left[\begin{array}{cc}a& 0\\ c& b\end{array}\right]$ where $a=\text{exp}(a^{*})$ and $b=\text{exp}(b^{*})$ . Then, we assume normal prior $N({\mathbf{\mu }}_{06},V_{06})$ for $(a^{*},b^{*},c)$ . We further assume independent gamma priors for ${\mathbf{\lambda }}_1$ and ${\mathbf{\lambda }}_2$ with ${\lambda }_{\mathit{jk}}\sim \mathit{Gamma}(a_j,b_j)$ , where $a_j\ge 0$ and $b_j\ge 0$ for $k=1,\dots ,K_j$ and $j=1,2$ . In Section 5 (the real data analysis), the pre-specified hyperparameters are given by ${\mathbf{\mu }}_{01}={\mathbf{\mu }}_{02}={\mathbf{\mu }}_{03}={\mathbf{\mu }}_{04}={\mathbf{\mu }}_{05}={\mathbf{\mu }}_{06}=0$ , $V_{01}=1000I_p$ , $V_{02}=1000I_{q+1}$ , $V_{03}=10I_{p+1}$ , $V_{04}=10000I_{p_3}$ with $p_3=\mathit{dim}(({\mathbf{\beta }}_1^{\prime },{\mathbf{\alpha }}_1^{\prime }{)}^{\prime })$ , $V_{05}=10000I_{p_4}$ with $p_4=\mathit{dim}(({\mathbf{\beta }}_2^{\prime },{\mathbf{\alpha }}_2^{\prime },{\mathbf{\eta }}^{\prime }{)}^{\prime })$ , $V_{06}=10I_2$ $a_0=0.001$ , $b_0=0.001$ , $d_0=q+2$ , ${\Omega }_0=0.001I_{q+1}$ , and $a_j=0.001$ and $b_j=0.001$ for $j=1,2$ , where $I_{p+q+1}$ , $I_{p+1}$ , $I_{p_3}$ , $I_{p_4}$ and $I_2$ are the identity matrices. These choices of the hyperparameters yield noninformative priors for all the model parameters except for $\mathbf{\varphi }$ .

Let ${\mathbf{\varphi }}_1=(\mathbf{\gamma },{\sigma }^2,\mathbf{\theta },\Gamma )$ , ${\mathbf{\varphi }}_2=(\mathbf{\varphi },{\mathbf{\varphi }}_{21})$ , where ${\mathbf{\varphi }}_{21}=({\mathbf{\lambda }}_1,{\mathbf{\beta }}_1,{\mathbf{\alpha }}_1)$ , and ${\mathbf{\varphi }}_3=({\mathbf{\lambda }}_2,{\mathbf{\beta }}_2,{\mathbf{\alpha }}_2,\mathbf{\eta })$ denote the collections of the model parameters for the longitudinal component, the survival component for the disease progression time and the survival component for the overall survival time, respectively, under the joint model. We further let $\mathbf{\varphi }=({\mathbf{\varphi }}_1,{\mathbf{\varphi }}_2,{\mathbf{\varphi }}_3)$ be the collection of all the model parameters. Then, the joint prior of $\mathbf{\varphi }$ is given by

$$\begin{array}{cc}\pi (\mathbf{\varphi })=\hfill & \pi (\mathbf{\gamma })\pi (\mathbf{\theta })\pi ({\sigma }^2)\pi (a^{*},b^{*},c)\pi (\mathbf{\varphi })\pi ({\mathbf{\beta }}_1,{\mathbf{\alpha }}_1)\pi ({\mathbf{\beta }}_2,{\mathbf{\alpha }}_2,\mathbf{\eta })\pi ({\mathbf{\lambda }}_1)\pi ({\mathbf{\lambda }}_2).\hfill \end{array}$$

(4.2)

Let $D_{\mathit{obs}}=({\mathit{y}}_i,D_{i1}1\{{\zeta }_i=1\},D_{i2}1\{{\zeta }_i=2\},D_{i3}1\{{\zeta }_i=3\},D_{i3}1\{{\zeta }_i=4\}, i=1,\dots ,n)$ denote the observed data, where the indictor function $1\{{\zeta }_i=j\}=1$ if ${\zeta }_i=j$ and 0 if ${\zeta }_i\ne j$ for $j=1,\dots ,4$ . The observed likelihood function is

$$\begin{array}{cc}\hfill & L(\mathbf{\varphi }|D_{\mathit{obs}})=\prod _{i=1}^n\int [f({\mathit{y}}_i|\mathbf{\gamma },{\sigma }^2,\mathbf{\theta },\Gamma ,{\mathbf{\theta }}_i^{*},{\mathit{x}}_i)\prod _{j=1}^4\{L({\mathbf{\varphi }}_1,{\mathbf{\varphi }}_2|D_{\mathit{ij}},{\mathbf{\theta }}_i^{*}){\}}^{1\{{\zeta }_i=j\}}]f({\mathbf{\theta }}_i^{*})d{\mathbf{\theta }}_i^{*},\hfill \end{array}$$

(4.3)

where $f({\mathit{y}}_i|\mathbf{\gamma },{\sigma }^2,\mathbf{\theta },\Gamma ,{\mathbf{\theta }}_i^{*},{\mathit{x}}_i)$ , $f({\mathbf{\theta }}_i^{*})$ , and $L({\mathbf{\varphi }}_1,{\mathbf{\varphi }}_2|D_{\mathit{ij}},{\mathbf{\theta }}_i^{*})$ , $j=1,\dots ,4$ , are given in (3.11, 3.12, 3.13) to (3.15) and (3.16), respectively, and $\pi (\mathbf{\varphi })$ is defined by (4.2). The joint posterior distribution of $\mathbf{\varphi }$ given the observed data is written as

$$\begin{array}{cc}\hfill & \pi (\mathbf{\varphi }|D_{\mathit{obs}})=\frac{L(\mathbf{\varphi }|D_{\mathit{obs}})\pi (\mathbf{\varphi })}{c(D_{\mathit{obs}})},\hfill \end{array}$$

(4.4)

where

$$\begin{array}{cc}\hfill & c(D_{\mathit{obs}})=\int L(\mathbf{\varphi }|D_{\mathit{obs}})\pi (\mathbf{\varphi })d\mathbf{\varphi }.\hfill \end{array}$$

(4.5)

4.2 Computational development

To ease exposition, we present the computational algorithm only for ${\beta }_2(\Gamma {\mathbf{\theta }}_i^{*}+\mathbf{\theta }{)}^{\prime }\mathit{g}(t)$ and ${\beta }_1(\Gamma {\mathbf{\theta }}_i^{*}+\mathbf{\theta }{)}^{\prime }\mathit{g}(t)$ in the survival functions. We can derive a similar function when assuming the term ${\beta }_1((\Gamma {\mathbf{\theta }}_i^{*}+\mathbf{\theta }{)}^{\prime }\mathit{g}({\mathit{Tmax}}_i))$ and ${\beta }_2(\Gamma {\mathbf{\theta }}_i^{*}+\mathbf{\theta }{)}^{\prime }\mathit{g}({\mathit{Tmax}}_i)$ after certain time point ${\mathit{Tmax}}_i$ . Models as a similar algorithm can be developed for the SPM. We introduce a latent vector ${\mathbf{\delta }}^{*}=({\delta }_i^{*}:{\zeta }_i=1\mathit{or}{\zeta }_i=4, 1\le i\le n{)}^{\prime }$ . $P({\delta }_i^{*}=1)=\frac{p(\mathbf{\varphi },{\mathit{x}}_i)*P(T_P>t_i|{\mathbf{\varphi }}_{21},{\mathit{x}}_i,\Gamma ,{\mathbf{\theta }}_i^{*},\mathbf{\theta },\delta =1)}{1-p(\mathbf{\varphi },{\mathit{x}}_i)+p(\mathbf{\varphi },{\mathit{x}}_i)*P(T_P>t_i|{\mathbf{\varphi }}_{21},{\mathit{x}}_i,\Gamma ,{\mathbf{\theta }}_i^{*},\mathbf{\theta },\delta =1)}$ and $P({\delta }_i^{*}=0)=1-P({\delta }_i^{*}=1)$ . Let ${{\mathbf{\theta }}^{*}}^R=({{\mathbf{\theta }}_1^{*}}^{\prime },\dots ,{{\mathbf{\theta }}_n^{*}}^{\prime }{)}^{\prime }$ denote the vector of the random effects. Then, the augmented posterior distribution of $\pi (\mathbf{\varphi }|D_{\mathit{obs}})$ in (4.4) can be written as

$$\begin{array}{cc}\hfill & \pi (\mathbf{\varphi },{\mathbf{\delta }}^{*},{{\mathbf{\theta }}^{*}}^R|D_{\mathit{obs}})\hfill \\ \hfill & \propto {\displaystyle \prod _{i=1}^n(\{f_2\left(t_i\mid {\phi }_3,x_i,z_i(t),\Gamma ,{\theta }_i^{*},\theta \right){\left[1-p\left(\varphi ,x_i\right)\right]}^{\left(1-{\delta }_i^{*}\right)}[p\left(\varphi ,x_i\right)}\hfill \\ \hfill & \times \{f_2(t_i|{\mathbf{\varphi }}_3,{\mathit{x}}_i,{\mathit{z}}_i(t),\Gamma ,{\mathbf{\theta }}_i^{*},\mathbf{\theta })p(\mathbf{\varphi },{\mathit{x}}_i)f_1(t_{\mathit{Pi}}|{\mathbf{\varphi }}_{21},{\mathit{x}}_i,\Gamma ,{\mathbf{\theta }}_i^{*},\mathbf{\theta },\delta =1){\}}^{1\{{\zeta }_i=2\}}\hfill \\ \hfill & \times \{P(T_D>t_i|{\mathbf{\varphi }}_3,{\mathit{x}}_i,{\mathit{z}}_i(t),\Gamma ,{\mathbf{\theta }}_i^{*},\mathbf{\theta })p(\mathbf{\varphi },{\mathit{x}}_i)f_1(t_{\mathit{Pi}}|{\mathbf{\varphi }}_{21},{\mathit{x}}_i,\Gamma ,{\mathbf{\theta }}_i^{*},\mathbf{\theta },\delta =1){\}}^{1\{{\zeta }_i=3\}}\hfill \\ \hfill & \times \{P(T_D>t_i|{\mathbf{\varphi }}_3,{\mathit{x}}_i,{\mathit{z}}_i(t),\Gamma ,{\mathbf{\theta }}_i^{*},\mathbf{\theta })[1-p(\mathbf{\varphi },{\mathit{x}}_i){]}^{(1-{\delta }_i^{*})}\hfill \\ \hfill & p(\mathbf{\varphi },{\mathit{x}}_i)P(T_P>t_i|{\mathbf{\varphi }}_{21},{\mathit{x}}_i,\Gamma ,{\mathbf{\theta }}_i^{*},\mathbf{\theta },\delta =1){]}^{{\delta }_i^{*}}{\}}^{1\{{\zeta }_i=4\}}\hfill \\ \hfill & \times f({\mathit{y}}_i|\mathbf{\gamma },{\sigma }^2,\mathbf{\theta },\Gamma ,{\mathbf{\theta }}_i^{*},{\mathit{x}}_i)f({\mathbf{\theta }}_i^{*}))\times \pi (\mathbf{\varphi }),\hfill \end{array}$$

(4.6)

where $f_2(t_i|{\mathbf{\varphi }}_3,{\mathit{x}}_i,{\mathit{z}}_i(t),\Gamma ,{\mathbf{\theta }}_i^{*},\mathbf{\theta })$ , $P(T_D>t_i|{\mathbf{\varphi }}_3,{\mathit{x}}_i,{\mathit{z}}_i(t),\Gamma ,{\mathbf{\theta }}_i^{*},\mathbf{\theta })$ , $f_1(t_{\mathit{Pi}}|{\mathbf{\varphi }}_{21},{\mathit{x}}_i,\Gamma ,{\mathbf{\theta }}_i^{*},\mathbf{\theta },\delta =1)$ and $P(T_P>t_i|{\mathbf{\varphi }}_{21},{\mathit{x}}_i,\Gamma ,{\mathbf{\theta }}_i^{*},\mathbf{\theta },\delta =1)$ are defined in Section 3.3. $f({\mathit{y}}_i|\mathbf{\gamma },{\sigma }^2,\mathbf{\theta },\Gamma ,{\mathbf{\theta }}_i^{*},{\mathit{x}}_i)$ , $f({\mathbf{\theta }}_i^{*})$ , $\pi (\mathbf{\varphi })$ are given in (3.11, 3.12) and (4.2). We note that $\int {\sum }_{{\mathbf{\delta }}^{*}}\pi (\mathbf{\varphi },{\mathbf{\delta }}^{*},{\mathbf{\theta }}^R|D_{\mathit{obs}})d{\mathbf{\theta }}^R=\pi (\mathbf{\varphi }|D_{\mathit{obs}})$ .

Although an analytical evaluation of the posterior distribution in (4.4) is not possible, a convenient Gibbs sampling algorithm can be developed to sample from the augmented posterior distribution in (4.6). Specifically, we sample from the following full conditional distributions in turn: (i) $[\mathbf{\gamma }|{\sigma }^2,\Gamma ,{\mathbf{\theta }}^{*R},\mathbf{\theta },D_{\mathit{obs}}]$ ; (ii) $[{\sigma }^2|{\mathbf{\theta }}^{*R},\mathbf{\theta },\Gamma ,D_{\mathit{obs}},\mathbf{\gamma }]$ ; (iii) $[(a^{*},b^{*},c)|\mathbf{\gamma },{\sigma }^2,\mathbf{\theta },{\mathbf{\theta }}^{*R},D_{\mathit{obs}}]$ ; (iv) $[{\mathbf{\delta }}^{*}|{\mathbf{\varphi }}_2,\Gamma ,{\mathbf{\theta }}_i^{*},\mathbf{\theta },D_{\mathit{obs}}]$ ; (v) $[\mathbf{\varphi }|{\mathbf{\delta }}^{*},D_{\mathit{obs}}]$ ; (vi) $[{\mathbf{\beta }}_1,{\mathbf{\alpha }}_1|{\mathbf{\lambda }}_1,\Gamma ,{\mathbf{\theta }}^{*R},\mathbf{\theta },{\mathbf{\delta }}^{*},D_{\mathit{obs}}]$ ; (vii) $[{\mathbf{\lambda }}_1|{\mathbf{\beta }}_1,{\mathbf{\alpha }}_1,\Gamma ,{\mathbf{\theta }}^{*R},\mathbf{\theta },{\mathbf{\delta }}^{*},D_{\mathit{obs}}]$ ; (viii) $[{\mathbf{\beta }}_2,{\mathbf{\alpha }}_2,\mathbf{\eta }|{\mathbf{\lambda }}_2,\Gamma ,{\mathbf{\theta }}^{*R},\mathbf{\theta },D_{\mathit{obs}}]$ ; (ix) $[{\mathbf{\lambda }}_2|{\mathbf{\alpha }}_2,{\mathbf{\beta }}_2,\Gamma ,{\mathbf{\theta }}^{*R},\mathbf{\theta },D_{\mathit{obs}}]$ ; (x) $[{\mathbf{\theta }}^{*R}|{\mathbf{\varphi }}_1$ , ${\mathbf{\varphi }}_{21},{\mathbf{\varphi }}_3,{\mathbf{\delta }}^{*},D_{\mathit{obs}}]$ ; and (xi) $[\mathbf{\theta }|{\mathbf{\varphi }}_1$ , ${\mathbf{\varphi }}_{21},{\mathbf{\varphi }}_3,{\mathbf{\delta }}^{*},{\mathbf{\theta }}^{*R},D_{\mathit{obs}}]$ . We briefly discuss the full conditional distributions in Section S.1 in the supplementary materials.

4.3 Bayesian model assessment

Let $f({\mathit{y}}_i|{\mathbf{\varphi }}_1,{\mathit{x}}_i)=\int f({\mathit{y}}_i|{\mathbf{\varphi }}_1,{\mathbf{\theta }}_i^{*},{\mathit{x}}_i)f({\mathbf{\theta }}_i^{*})d{\mathbf{\theta }}_i^{*}$ , where marginal distribution of ${\mathit{y}}_i$ follows ${\mathit{y}}_i|{\mathbf{\varphi }}_1^{*},{\mathit{x}}_i\sim N(W_i({\mathbf{\theta }}^{\prime },{\mathbf{\gamma }}^{\prime }{)}^{\prime },{\sigma }^2{I}_{m_i}+W_{i1}\Gamma {\Gamma }^{\prime }{W}_{i1}^{\prime })$ . Let $f(t_{\mathit{Pi}},{\delta }_i|{\mathbf{\varphi }}_1,{\mathbf{\varphi }}_2,{\mathbf{\theta }}_i^{*},{\mathit{x}}_i)=\{1-p(\mathbf{\varphi },{\mathit{x}}_i)+p(\mathbf{\varphi },{\mathit{x}}_i)P(T_P>t_{\mathit{Pi}}|{\mathbf{\varphi }}_{21},{\mathit{x}}_i,\Gamma ,{\mathbf{\theta }}_i^{*},\mathbf{\theta },\delta =1){\}}^{1\{{\zeta }_i=1\}+1\{{\zeta }_i=4\}}\{p(\mathbf{\varphi },{\mathit{x}}_i)f_1(t_{\mathit{Pi}}|{\mathbf{\varphi }}_{21},{\mathit{x}}_i,\Gamma ,{\mathbf{\theta }}_i^{*},\mathbf{\theta },\delta =1){\}}^{1\{{\zeta }_i=2\}+1\{{\zeta }_i=3\}}$ , $f(t_{\mathit{Pi}},{\delta }_i,{\mathit{y}}_i,{\mathbf{\theta }}_i^{*}|{\mathbf{\varphi }}_2,{\mathbf{\varphi }}_1,{\mathit{x}}_i)=f(t_{\mathit{Pi}},{\delta }_i|{\mathbf{\varphi }}_2,\Gamma ,{\mathbf{\theta }}_i^{*},\mathbf{\theta },{\mathit{x}}_i)f({\mathit{y}}_i|{\mathbf{\varphi }}_1,{\mathbf{\theta }}_i^{*},{\mathit{x}}_i)f({\mathbf{\theta }}_i^{*})$ , and $f(t_{\mathit{Pi}},{\delta }_i,{\mathit{y}}_i|{\mathbf{\varphi }}_1,{\mathbf{\varphi }}_2,{\mathit{x}}_i)=\int f(t_{\mathit{Pi}},{\delta }_i,{\mathit{y}}_i,{\mathbf{\theta }}_i^{*}|{\mathbf{\varphi }}_2,{\mathbf{\varphi }}_1,{\mathit{x}}_i)d{\mathbf{\theta }}_i^{*}$ , $f(t_i|{\mathbf{\varphi }}_1,{\mathbf{\varphi }}_3,{\mathbf{\theta }}_i^{*},{\mathit{x}}_i,{\mathit{z}}_i(t),{\delta }_i,{\nu }_i)=\{f_2(t_i|{\mathbf{\varphi }}_3,{\mathit{x}}_i,{\mathit{z}}_i(t),\Gamma ,{\mathbf{\theta }}_i^{*},\mathbf{\theta }){\}}^{1\{{\zeta }_i=1\}+1\{{\zeta }_i=2\}}\{P(T_D>t_i|{\mathbf{\varphi }}_3,{\mathit{x}}_i,{\mathit{z}}_i(t),\Gamma ,{\mathbf{\theta }}_i^{*},\mathbf{\theta }){\}}^{1\{{\zeta }_i=3\}+1\{{\zeta }_i=4\}}$ .

Define the joint distribution of of $({\mathit{y}}_i,t_{\mathit{Pi}},{\delta }_i,t_i,{\mathbf{\theta }}_i^{*})$ to be

$$\begin{array}{cc}f({\mathit{y}}_i,t_{\mathit{Pi}},\hfill & {\delta }_i,t_i,{\mathbf{\theta }}_i^{*}|\mathbf{\varphi },{\mathit{x}}_i,{\mathit{z}}_i(t),{\nu }_i)=f(t_i|{\mathbf{\varphi }}_1,{\mathbf{\varphi }}_3,{\mathbf{\theta }}_i^{*},{\mathit{x}}_i,{\mathit{z}}_i(t),{\delta }_i,{\nu }_i)\hfill \\ \hfill & f(t_{\mathit{Pi}},{\delta }_i|{\mathbf{\varphi }}_2,\Gamma ,{\mathbf{\theta }}_i^{*},\mathbf{\theta },{\mathit{x}}_i)f({\mathit{y}}_i|{\mathbf{\varphi }}_1,{\mathbf{\theta }}_i^{*},{\mathit{x}}_i)f({\mathbf{\theta }}_i^{*}),\hfill \end{array}$$

(4.7)

and the marginal distribution $({\mathit{y}}_i,t_{\mathit{Pi}},{\delta }_i,t_i)$

$$\begin{array}{cc}f({\mathit{y}}_i,t_{\mathit{Pi}},{\delta }_i,t_i|\mathbf{\varphi },{\mathit{x}}_i,{\mathit{z}}_i(t),{\nu }_i)=\int f({\mathit{y}}_i,t_{\mathit{Pi}},{\delta }_i,t_i,{\mathbf{\theta }}_i^{*}|\mathbf{\varphi },{\mathit{x}}_i,{\mathit{z}}_i(t),{\nu }_i)d{\mathbf{\theta }}_i^{*}.\hfill & \hfill \end{array}$$

(4.8)

Based on the joint distributions, define the following conditional distributions: $f(t_{\mathit{Pi}},{\delta }_i|{\mathbf{\varphi }}_2$ , ${\mathbf{\varphi }}_1,{\mathit{x}}_i,{\mathit{y}}_i)=f(t_{\mathit{Pi}},{\delta }_i,{\mathit{y}}_i|{\mathbf{\varphi }}_1,{\mathbf{\varphi }}_2,{\mathit{x}}_i)/f({\mathit{y}}_i|{\mathbf{\varphi }}_1,{\mathit{x}}_i)$ and $f(t_i|\mathbf{\varphi },{\mathit{x}}_i,{\mathit{z}}_i(t),{\mathit{y}}_i,{\delta }_i,{\nu }_i)=f({\mathit{y}}_i,t_{\mathit{Pi}},{\delta }_i,t_i|\mathbf{\varphi },{\mathit{x}}_i,{\mathit{z}}_i(t),{\nu }_i)/f(t_{\mathit{Pi}},{\delta }_i,{\mathit{y}}_i|{\mathbf{\varphi }}_1,{\mathbf{\varphi }}_2,{\mathit{x}}_i)$ .

4.3.1 DIC decomposition

The deviance information criterion (DIC; Spiegelhalter et al. (2002)) for the model is defined as

$$\mathit{DIC}=\mathit{Dev}(\bar{\mathbf{\varphi }})+2p_D,$$

(4.9)

where $\mathit{Dev}(\mathbf{\varphi })=-2\text{log}L(\mathbf{\varphi }|D_{\mathit{obs}})$ , $p_D=\overline{\mathit{Dev}}(\mathbf{\varphi })-\mathit{Dev}(\bar{\mathbf{\varphi }})$ , where $\bar{\mathbf{\varphi }}$ and $\overline{\mathit{Dev}}(\mathbf{\varphi })$ are the posterior means of $\mathbf{\varphi }$ and $\mathit{Dev}(\mathbf{\varphi })$ , respectively, with respect to the posterior distribution in (4.4).

To assess the fit of the longitudinal component of the model and the fit of each survival component, we decompose DIC into three parts:

$$\begin{array}{cc}\mathit{DIC}={\mathit{DIC}}_{\mathit{Long}}+{\mathit{DIC}}_{\mathit{Pg}|\mathit{Long}}+{\mathit{DIC}}_{\mathit{Surv}|\mathit{Pg},\mathit{Long}},\hfill & \hfill \end{array}$$

where ${\mathit{DIC}}_{\mathit{Long}}={\mathit{Dev}}_{\mathit{Long}}(\bar{\mathbf{\varphi }})+2p_{D[\mathit{Long}]}$ , ${\mathit{DIC}}_{\mathit{Pg}|\mathit{Long}}={\mathit{Dev}}_{\mathit{Pg}|\mathit{Long}}(\bar{\mathbf{\varphi }})+2p_{D[\mathit{Pg}|\mathit{Long}]}$ , and ${\mathit{DIC}}_{\mathit{Surv}|\mathit{Pg},\mathit{Long}}={\mathit{Dev}}_{\mathit{Surv}|\mathit{Pg},\mathit{Long}}(\bar{\mathbf{\varphi }})+2p_{D[\mathit{Surv}|\mathit{Pg},\mathit{Long}]}$ .

Write ${\mathit{Dev}}_{\mathit{Long}}(\bar{\mathbf{\varphi }})=-2{\sum }_{i=1}^n\text{log}f({\mathit{y}}_i|\overline{{\mathbf{\varphi }}_1},{\mathit{x}}_i)$ , $p_{D[\mathit{Long}]}=E[-2{\sum }_{i=1}^n\text{log}f({\mathit{y}}_i|{\mathbf{\varphi }}_1,{\mathit{x}}_i)|D_{\mathit{obs}}]+2{\sum }_{i=1}^n\text{log}f({\mathit{y}}_i|\overline{{\mathbf{\varphi }}_1},{\mathit{x}}_i)$ , ${\mathit{Dev}}_{\mathit{Long}}(\bar{\mathbf{\varphi }})+{\mathit{Dev}}_{\mathit{Pg}|\mathit{Long}}(\bar{\mathbf{\varphi }})=-2{\sum }_{i=1}^n\text{log}\int [f(t_{\mathit{Pi}},{\delta }_i|\overline{{\mathbf{\varphi }}_1},\overline{{\mathbf{\varphi }}_2},{\mathbf{\theta }}_i^{*},{\mathit{x}}_i)$ $\times f({\mathit{y}}_i|\overline{{\mathbf{\varphi }}_1},{\mathbf{\theta }}_i^{*},{\mathit{x}}_i)f({\mathbf{\theta }}_i^{*})]d{\mathbf{\theta }}_i^{*}$ , and $p_{D[\mathit{Long}]}+p_{D[\mathit{Pg}|\mathit{Long}]}=E[-2{\sum }_{i=1}^n\text{log}\int f(t_{\mathit{Pi}},{\delta }_i|{\mathbf{\varphi }}_1,{\mathbf{\varphi }}_2,{\mathbf{\theta }}_i^{*}$ , ${\mathit{x}}_i)f({\mathit{y}}_i|{\mathbf{\varphi }}_1,{\mathbf{\theta }}_i^{*},{\mathit{x}}_i)f({\mathbf{\theta }}_i^{*})d{\mathbf{\theta }}_i^{*}|D_{\mathit{obs}}]+2{\sum }_{i=1}^n\text{log}\int f(t_{\mathit{Pi}},{\delta }_i|\overline{{\mathbf{\varphi }}_1},\overline{{\mathbf{\varphi }}_2},{\mathbf{\theta }}_i^{*},{\mathit{x}}_i)f({\mathit{y}}_i|\overline{{\mathbf{\varphi }}_1},{\mathbf{\theta }}_i^{*},{\mathit{x}}_i)f({\mathbf{\theta }}_i^{*})d{\mathbf{\theta }}_i^{*}$ . Then, we have ${\mathit{Dev}}_{\mathit{Long}}(\bar{\mathbf{\varphi }})+{\mathit{Dev}}_{\mathit{Pg}|\mathit{Long}}(\bar{\mathbf{\varphi }})+{\mathit{Dev}}_{\mathit{Surv}|\mathit{Pg},\mathit{Long}}(\bar{\mathbf{\varphi }})=\overline{\mathit{Dev}}(\mathbf{\varphi })$ and $p_{D[\mathit{Long}]}+p_{D[\mathit{Pg}|\mathit{Long}]}+p_{D[\mathit{Surv}|\mathit{Pg},\mathit{Long}]}=p_D$ .

When we fit the survival data alone, assuming ${\beta }_2=0$ and $\mathbf{\eta }=0$ . Define ${\lambda }_2^{*}(t|{\lambda }_{20},{\mathbf{\alpha }}_2,\mathit{x})={\lambda }_{20}(t)\text{exp}\{{\mathbf{\alpha }}_2^{\prime }\mathit{x}\}$ , and $f_0(t_i|{\mathbf{\lambda }}_2,{\mathbf{\alpha }}_2,{\mathit{x}}_i,{\nu }_i)=[{\lambda }_2^{*}(t_i|{\lambda }_{20},{\mathbf{\alpha }}_2,\mathit{x}){]}^{{\nu }_i}\text{exp}[-{\int }_0^{t_i}{\lambda }_2^{*}(u|{\lambda }_{20},{\mathbf{\alpha }}_2,\mathit{x})\mathit{du}]$ .

$$\begin{array}{cc}{\mathit{DIC}}_{\mathit{Surv}0}={\mathit{Dev}}_{\mathit{Surv}0}(\overline{{\mathbf{\lambda }}_2},\overline{{\mathbf{\alpha }}_2})+2p_{D[\mathit{Surv}0]},\hfill & \hfill \end{array}$$

where ${\mathit{Dev}}_{\mathit{Surv}0}(\overline{{\mathbf{\lambda }}_2},\overline{{\mathbf{\alpha }}_2})=-2{\sum }_{i=1}^n\text{log}{f}_0(t_i|\overline{{\mathbf{\lambda }}_2},\overline{{\mathbf{\alpha }}_2},{\mathit{x}}_i,{\nu }_i)$ and $p_{D[\mathit{Surv}0]}=E[-2{\sum }_{i=1}^n\text{log}\{$ $f_0(t_i|{\mathbf{\lambda }}_2,{\mathbf{\alpha }}_2,{\mathit{x}}_i,{\nu }_i)|D_{\mathit{Surv},\mathit{obs}}\}]+2{\sum }_{i=1}^n\text{log}{f}_0(t_i|\overline{{\mathbf{\lambda }}_2},\overline{{\mathbf{\alpha }}_2},{\mathit{x}}_i,{\nu }_i)$ . Define model assessment criterion:

$$\begin{array}{cc}\Delta {\mathit{DIC}}_{\mathit{Surv}}={\mathit{DIC}}_{\mathit{Surv}0}-{\mathit{DIC}}_{\mathit{Surv}|\mathit{Pg},\mathit{Long}}.\hfill & \hfill \end{array}$$

$\Delta {\mathit{DIC}}_{\mathit{Surv}}$ measures the additional information gained by adding time to disease progression and longitudinal components into modelling of overall survival data. Large values of $\Delta {\mathit{DIC}}_{\mathit{Surv}}$ correspond to more gain in the fit of the overall survival data. We can also fit survival data alone with time-dependent covariates, assuming ${\beta }_2=0$ .

Similarly, denote

$$\begin{array}{cc}\Delta {\mathit{DIC}}_{\mathit{Pg}}={\mathit{DIC}}_{\mathit{Pg}0}-{\mathit{DIC}}_{\mathit{Pg}|\mathit{Long}},\hfill & \hfill \end{array}$$

where ${\mathit{DIC}}_{\mathit{Pg}0}$ is obtained by assuming ${\beta }_1=0$ (fitting the time to progression model alone). Define ${\lambda }_1^{*}(t|{\lambda }_{10},{\mathbf{\alpha }}_1,\mathit{x})={\lambda }_{10}(t)\text{exp}\{{\mathbf{\alpha }}_1^{\prime }\mathit{x}\}$ , $P^{*}(T_P>t_{\mathit{Pi}}|{\lambda }_{10},{\mathbf{\alpha }}_1,{\mathit{x}}_i,\delta =1)=\text{exp}[-{\int }_0^{t_{\mathit{Pi}}}{\lambda }_1^{*}(u|{\lambda }_{10},{\mathbf{\alpha }}_1,\mathit{x})\mathit{du}]$ , and $f_0(t_{\mathit{Pi}},{\delta }_i|{\mathbf{\lambda }}_1,{\mathbf{\alpha }}_1,\mathbf{\varphi },{\mathit{x}}_i)=\{1-p(\mathbf{\varphi },{\mathit{x}}_i)+p(\mathbf{\varphi },{\mathit{x}}_i)P^{*}(T_P>t_{\mathit{Pi}}|{\lambda }_{10},{\mathbf{\alpha }}_1,{\mathit{x}}_i,\delta =1){\}}^{1\{{\zeta }_i=1\}+1\{{\zeta }_i=4\}}\{p(\mathbf{\varphi },{\mathit{x}}_i){\lambda }_1^{*}(t|{\lambda }_{10},{\mathbf{\alpha }}_1,\mathit{x})P^{*}(T_P>t_{\mathit{Pi}}|{\lambda }_{10},{\mathbf{\alpha }}_1,{\mathit{x}}_i,\delta =1){\}}^{1\{{\zeta }_i=2\}+1\{{\zeta }_i=3\}}$ . We have

$${\mathit{DIC}}_{\mathit{Pg}0}={\mathit{Dev}}_{\mathit{Pg}0}(\overline{{\mathbf{\lambda }}_1},\overline{{\mathbf{\alpha }}_1},\bar{\mathbf{\varphi }})+2p_{D[\mathit{Pg}0]},$$

where ${\mathit{Dev}}_{\mathit{Pg}0}({\bar{\mathbf{\lambda }}}_1,{\bar{\mathbf{\alpha }}}_1,\bar{\mathbf{\varphi }})=-2{\sum }_{i=1}^n\text{log}{f}_0(t_{\mathit{Pi}},{\delta }_i|\overline{{\mathbf{\lambda }}_1},\overline{{\mathbf{\alpha }}_1},\bar{\mathbf{\varphi }},{\mathit{x}}_i)$ and

$$p_{D[\mathit{Pg}0]}=E[-2\sum _{i=1}^n\text{log}{f}_0(t_{\mathit{Pi}},{\delta }_i|{\mathbf{\lambda }}_1,{\mathbf{\alpha }}_1,\mathbf{\varphi },{\mathit{x}}_i)|D_{\mathit{Pg},\mathit{obs}}]+2\sum _{i=1}^n\text{log}{f}_0(t_{\mathit{Pi}},{\delta }_i|{\bar{\mathbf{\lambda }}}_1,{\bar{\mathbf{\alpha }}}_1,\bar{\mathbf{\varphi }},{\mathit{x}}_i),$$

and $D_{\mathit{Pg},\mathit{obs}}=(t_{\mathit{Pi}},{\delta }_i,{\mathit{x}}_i i=1,\dots ,n)$ . $\Delta {\mathit{DIC}}_{\mathit{Pg}}$ measures the additional information gained by longitudinal components into modelling of progression data. Large values of $\Delta {\mathit{DIC}}_{\mathit{Pg}}$ correspond to more gain in the fit of progression data.

As pointed out in Zhang et al. (2017), when the penalty for the additional parameters in the overall survival model (3.10) or the time to progression model (3.8) outweighs the gain of the fit in the respective survival component, $\Delta {\mathit{DIC}}_{\mathit{Surv}}$ or $\Delta {\mathit{DIC}}_{\mathit{Pg}}$ can be negative. When $\Delta {\mathit{DIC}}_{\mathit{Surv}}$ ( $\Delta {\mathit{DIC}}_{\mathit{Pg}}$ ) is zero or negative, the longitudinal measure is not associated with the overall survival (the time to disease progression). This phenomenon is also applied to the $\Delta$ measures developed based on the CPO statistics below.

4.3.2 CPO decomposition and LPML decomposition

In this subsection, we first establish the conditional predictive ordinate (CPO) identity which will lead to the decomposition of CPO. CPO (Geisser and Eddy, 1979; Chen et al., 2000; Zhang et al., 2017) for the $i$ th subject has the following equivalent definition and identities:

$$\begin{array}{cc}\hfill & {\mathit{CPO}}_i=\int f({\mathit{y}}_i,t_{\mathit{Pi}},{\delta }_i,t_i|\mathbf{\varphi },{\mathit{x}}_i,{\mathit{z}}_i(t),{\nu }_i)\pi (\mathbf{\varphi }|D_{\mathit{obs}}^{(-i)})d\mathbf{\varphi },\hfill \\ \hfill & {\mathit{CPO}}_i=\frac{c(D_{\mathit{obs}})}{c(D_{\mathit{obs}}^{(-i)})}=\frac{f({\mathit{y}}_i,t_{\mathit{Pi}},{\delta }_i,t_i|\mathbf{\varphi },{\mathit{x}}_i,{\mathit{z}}_i(t),{\nu }_i)\pi (\mathbf{\varphi }|D_{\mathit{obs}}^{(-i)})}{\pi (\mathbf{\varphi }|D_{\mathit{obs}})},\hfill \\ \hfill & {\mathit{CPO}}_i=\frac{1}{\int \frac{1}{f({\mathit{y}}_i,t_{\mathit{Pi}},{\delta }_i,t_i|\mathbf{\varphi },{\mathit{x}}_i,{\mathit{z}}_i(t),{\nu }_i)}\pi (\mathbf{\varphi }|D_{\mathit{obs}})d\mathbf{\varphi }},\hfill \end{array}$$

(4.13)

where $c(D_{\mathit{obs}})$ is defined in (4.5), $c(D_{\mathit{obs}}^{(-i)})=\int {\Pi }_{j\ne i}^{n}f({\mathit{y}}_j,t_{\mathit{Pj}},{\delta }_j,t_j|\mathbf{\varphi },{\mathit{x}}_j,{\nu }_j)\pi (\mathbf{\varphi })d\mathbf{\varphi }$ , and $f({\mathit{y}}_i,t_{\mathit{Pi}},{\delta }_i,t_i|\mathbf{\varphi },{\mathit{x}}_i,{\nu }_i)$ is defined in (4.8).

Let ${\mathbf{\varphi }}_1^{*}$ , ${\mathbf{\varphi }}_2^{*}$ and ${\mathbf{\varphi }}_3^{*}$ denote the posterior means ${\mathbf{\varphi }}_1$ , ${\mathbf{\varphi }}_2$ and ${\mathbf{\varphi }}_3$ , and ${\mathbf{\varphi }}^{*}=({\mathbf{\varphi }}_1^{*},{\mathbf{\varphi }}_2^{*},{\mathbf{\varphi }}_3^{*})$ . We have

$$\begin{array}{cc}f({\mathit{y}}_i,t_{\mathit{Pi}},\hfill & {\delta }_i,t_i|{\mathbf{\varphi }}^{*},{\mathit{x}}_i,{\mathit{z}}_i(t),{\nu }_i)=f(t_i|{\mathbf{\varphi }}^{*},{\mathit{x}}_i,{\mathit{z}}_i(t),{\delta }_i,{\nu }_i,{\mathit{y}}_i)f(t_{\mathit{Pi}},{\delta }_i|{\mathit{y}}_i,{\mathbf{\varphi }}_1^{*},{\mathbf{\varphi }}_2^{*},{\mathit{x}}_i)\hfill \\ \hfill & \times f({\mathit{y}}_i|{\mathbf{\varphi }}_1^{*},{\mathit{x}}_i).\hfill \end{array}$$

(4.14)

We also observe that $\pi ({\mathbf{\varphi }}^{*}|D_{\mathit{obs}}^{(-i)})=\pi ({\mathbf{\varphi }}_1^{*}|D_{\mathit{obs}}^{(-i)})\pi ({\mathbf{\varphi }}_2^{*}|{\mathbf{\varphi }}_1^{*},D_{\mathit{obs}}^{(-i)})\pi ({\mathbf{\varphi }}_3^{*}|{\mathbf{\varphi }}_2^{*},{\mathbf{\varphi }}_1^{*},D_{\mathit{obs}}^{(-i)})$ , $\pi ({\mathbf{\varphi }}^{*}|D_{\mathit{obs}})=\pi ({\mathbf{\varphi }}_1^{*}|D_{\mathit{obs}})\pi ({\mathbf{\varphi }}_2^{*}|{\mathbf{\varphi }}_1^{*},D_{\mathit{obs}})\pi ({\mathbf{\varphi }}_3^{*}|{\mathbf{\varphi }}_2^{*},{\mathbf{\varphi }}_1^{*},D_{\mathit{obs}})$ . We can propose the CPO decomposition:

$$\begin{array}{cc}{\mathit{CPO}}_i={\mathit{CPO}}_{i,\mathit{Long}}*{\mathit{CPO}}_{i,\mathit{Pg}|\mathit{Long}}*{\mathit{CPO}}_{i,\mathit{Surv}|\mathit{Pg},\mathit{Long}},\hfill & \hfill \end{array}$$

(4.15)

where ${\mathit{CPO}}_{i,\mathit{Long}}=\frac{f({\mathit{y}}_i|{\mathbf{\varphi }}_1^{*},{\mathit{x}}_i)\pi ({\mathbf{\varphi }}_1^{*}|D_{\mathit{obs}}^{(-i)})}{\pi ({\mathbf{\varphi }}_1^{*}|D_{\mathit{obs}})}$ , ${\mathit{CPO}}_{i,\mathit{Pg}|\mathit{Long}}=\frac{f(t_{\mathit{Pi}},{\delta }_i|{\mathit{y}}_i,{\mathbf{\varphi }}_1^{*},{\mathbf{\varphi }}_2^{*},{\mathit{x}}_i)\pi ({\mathbf{\varphi }}_2^{*}|{\mathbf{\varphi }}_1^{*},D_{\mathit{obs}}^{(-i)})}{\pi ({\mathbf{\varphi }}_2^{*}|{\mathbf{\varphi }}_1^{*},D_{\mathit{obs}})}$ , ${\mathit{CPO}}_{i,\mathit{Surv}|\mathit{Pg},\mathit{Long}}=\frac{f(t_i|{\mathbf{\varphi }}_1^{*},{\mathbf{\varphi }}_2^{*},{\mathbf{\varphi }}_3^{*},{\mathit{x}}_i,t_{\mathit{Pi}},{\delta }_i,{\nu }_i)\pi ({\mathbf{\varphi }}_3^{*}|{\mathbf{\varphi }}_2^{*},{\mathbf{\varphi }}_1^{*},D_{\mathit{obs}}^{(-i)})}{\pi ({\mathbf{\varphi }}_3^{*}|{\mathbf{\varphi }}_2^{*},{\mathbf{\varphi }}_1^{*},D_{\mathit{obs}})}$ .

Let $D_{\mathit{Long},\mathit{obs}}^{(-i)}=\{({\mathit{y}}_j,{\mathit{x}}_j),j=1,...,i-1,...,n\}$ , $D_{\mathit{Pg},\mathit{obs}}^{(-i)}=\{(t_{\mathit{Pj}},{\delta }_j,{\mathit{x}}_j),j=1,...,i-1,...,n\}$ , and $D_{\mathit{Surv},\mathit{obs}}^{(-i)}=\{(t_j,{\nu }_j,{\mathit{x}}_j),j=1,...,i-1,...,n\}$ . When we assume that $D_{\mathit{Long},\mathit{obs}}$ , $D_{\mathit{Pg},\mathit{obs}}$ and $D_{\mathit{Surv},\mathit{obs}}$ are mutually independent. Assume $\pi (\mathbf{\varphi })=\pi ({\mathbf{\varphi }}_1)\pi ({\mathbf{\varphi }}_2)\pi ({\mathbf{\varphi }}_3)$ . We have

$$\begin{array}{cc}\hfill & {\mathit{CPO}}_{i,\mathit{Surv}|\mathit{Pg},\mathit{Long}}={\mathit{CPO}}_{i,\mathit{Surv}0}=\int f_0(t_i|{\mathbf{\lambda }}_2,{\mathbf{\alpha }}_2,{\mathit{x}}_i,{\nu }_i)\pi ({\mathbf{\varphi }}_3|D_{\mathit{Surv},\mathit{obs}}^{(-i)})d{\mathbf{\varphi }}_3,\hfill \\ \hfill & {\mathit{CPO}}_{i,\mathit{Pg}|\mathit{Long}}={\mathit{CPO}}_{i,\mathit{Pg}0}=\int f_0(t_{\mathit{Pi}},{\delta }_i|{\mathbf{\lambda }}_1,{\mathbf{\alpha }}_1,\mathbf{\varphi },{\mathit{x}}_i)\pi ({\mathbf{\varphi }}_2|D_{\mathit{Pg},\mathit{obs}}^{(-i)})d{\mathbf{\varphi }}_2.\hfill \end{array}$$

From Section S.2 in the supplementary materials, we have the following identities for ${\mathit{CPO}}_i$ , ${\mathit{CPO}}_{i,\mathit{Long}}$ , ${\mathit{CPO}}_{i,\mathit{Pg}|\mathit{Long}}$ and ${\mathit{CPO}}_{i,\mathit{Surv}|\mathit{Pg},\mathit{Long}}$ :

$$\begin{array}{cc}\hfill & {\mathit{CPO}}_{i,\mathit{Long}}={\mathit{CPO}}_{i}f({\mathit{y}}_i|{\mathbf{\varphi }}_1^{*},{\mathit{x}}_i)\hfill \\ \hfill & \times \int f^{-1}({\mathit{y}}_i,t_{\mathit{Pi}},{\delta }_i,t_i|{\mathbf{\varphi }}_1^{*},{\mathbf{\varphi }}_2,{\mathbf{\varphi }}_3,{\mathit{x}}_i,{\mathit{z}}_i(t),{\nu }_i)\pi ({\mathbf{\varphi }}_2,{\mathbf{\varphi }}_3|{\mathbf{\varphi }}_1^{*},D_{\mathit{obs}})d{\mathbf{\varphi }}_{2}d{\mathbf{\varphi }}_3,\hfill \\ \hfill & {\mathit{CPO}}_{i,\mathit{Pg}|\mathit{Long}}=f(t_{\mathit{Pi}},{\delta }_i|{\mathit{y}}_i,{\mathbf{\varphi }}_1^{*},{\mathbf{\varphi }}_2^{*},{\mathit{x}}_i)\hfill \\ \hfill & \times \frac{[\int f^{-1}({\mathit{y}}_i,t_{\mathit{Pi}},{\delta }_i,t_i|{\mathbf{\varphi }}_1^{*},{\mathbf{\varphi }}_2^{*},{\mathbf{\varphi }}_3,{\mathit{x}}_i,{\mathit{z}}_i(t),{\nu }_i)\pi ({\mathbf{\varphi }}_3|{\mathbf{\varphi }}_1^{*},{\mathbf{\varphi }}_2^{*},D_{\mathit{obs}})d{\mathbf{\varphi }}_3]}{[\int f^{-1}({\mathit{y}}_i,t_{\mathit{Pi}},{\delta }_i,t_i|{\mathbf{\varphi }}_1^{*},{\mathbf{\varphi }}_2,{\mathbf{\varphi }}_3,{\mathit{x}}_i,{\mathit{z}}_i(t),{\nu }_i)\pi ({\mathbf{\varphi }}_2,{\mathbf{\varphi }}_3|{\mathbf{\varphi }}_1^{*},D_{\mathit{obs}})d{\mathbf{\varphi }}_{2}d{\mathbf{\varphi }}_3]},\hfill \\ \hfill & {\mathit{CPO}}_{i,\mathit{Surv}|\mathit{Pg},\mathit{Long}}=[\int f^{-1}(t_i|{\mathbf{\varphi }}_1^{*},{\mathbf{\varphi }}_2^{*},{\mathbf{\varphi }}_3,{\mathit{x}}_i,t_{\mathit{Pi}},{\delta }_i,{\nu }_i,{\mathit{y}}_i)\pi ({\mathbf{\varphi }}_3|{\mathbf{\varphi }}_1^{*},{\mathbf{\varphi }}_2^{*},D_{\mathit{obs}})d{\mathbf{\varphi }}_3{]}^{-1}.\hfill \end{array}$$

The logarithm of the pseudo-marginal likelihood (LPML; Ibrahim et al., 2001) is defined as

$$\begin{array}{cc}\mathit{LPML}=\sum _{i=1}^n\text{log}({\mathit{CPO}}_i).\hfill & \hfill \end{array}$$

Then LPML can be decomposed as

$$\begin{array}{cc}\mathit{LPML}={\mathit{LPML}}_{\mathit{Long}}+{\mathit{LPML}}_{\mathit{Pg}|\mathit{Long}}+{\mathit{LPML}}_{\mathit{Surv}|\mathit{Pg},\mathit{Long}},\hfill & \hfill \end{array}$$

where ${\mathit{LPML}}_{\mathit{Long}}={\sum }_{i=1}^n\text{log}({\mathit{CPO}}_{i,\mathit{Long}})$ , ${\mathit{LPML}}_{\mathit{Pg}|\mathit{Long}}={\sum }_{i=1}^n\text{log}({\mathit{CPO}}_{i,\mathit{Pg}|\mathit{Long}})$ and ${\mathit{LPML}}_{\mathit{Surv}|\mathit{Pg},\mathit{Long}}={\sum }_{i=1}^n\text{log}({\mathit{CPO}}_{i,\mathit{Surv}|\mathit{Pg},\mathit{Long}})$ .

Under the same assumptions for the calculation of ${\mathit{DIC}}_{\mathit{Surv}0}$ , define ${\mathit{LPML}}_{\mathit{Surv}0}={\sum }_{i=1}^n\text{log}{\mathit{CPO}}_{i,\mathit{Surv}0}$ . We propose the model assessment criterion:

$$\begin{array}{cc}\Delta {\mathit{LPML}}_{\mathit{Surv}}={\mathit{LPML}}_{\mathit{Surv}|\mathit{Pg},\mathit{Long}}-{\mathit{LPML}}_{\mathit{Surv}0}.\hfill & \hfill \end{array}$$

$\Delta {\mathit{LPML}}_{\mathit{Surv}}$ measures the additional information gained by adding time to disease progression and longitudinal components into modelling of overall survival data. Large values of $\Delta {\mathit{LPML}}_{\mathit{Surv}}$ correspond to more gain in the fit of the overall survival data. Similarly, denote $\Delta {\mathit{LPML}}_{\mathit{Pg}}={\mathit{LPML}}_{\mathit{Pg}|\mathit{Long}}-{\mathit{LPML}}_{\mathit{Pg}0}$ . ${\mathit{LPML}}_{\mathit{Pg}0}$ is from assuming ${\beta }_1=0$ and fit time to progression model alone. $\Delta {\mathit{LPML}}_{\mathit{Pg}}$ measures the additional information gained by longitudinal components into modelling of progression data. Large values of $\Delta {\mathit{LPML}}_{\mathit{Pg}}$ correspond to more gain in the fit of progression data.

Algorithm I. For $f({\mathit{y}}_i,t_{\mathit{Pi}},{\delta }_i,t_i|{\mathbf{\varphi }}_1,{\mathbf{\varphi }}_2,{\mathbf{\varphi }}_3,{\mathit{x}}_i,{\mathit{z}}_i(t),{\nu }_i)$ , $f(t_{\mathit{Pi}},{\delta }_i,{\mathit{y}}_i{\mathbf{\varphi }}_1,{\mathbf{\varphi }}_2,{\mathit{x}}_i)$ can be approximated by adaptive Gaussian quadrature (AGQ).

Now, let $\{{\mathbf{\varphi }}_b,b=1,\dots ,B\}$ denote a Gibbs sample of $\mathbf{\varphi }$ . An estimate of ${\mathit{CPO}}_{\mathit{i}}^{-1}$ is given by

$${\widehat{\mathit{CPO}}}_i^{-1}=\frac{1}{B}\sum _{b=1}^B{f}^{-1}({\mathit{y}}_i,t_{\mathit{Pi}},{\delta }_i,t_i|{\mathbf{\varphi }}_b,{\mathit{x}}_i,{\mathit{z}}_i(t),{\nu }_i).$$

Next, fix ${\mathbf{\varphi }}_1={\mathbf{\varphi }}_1^{*}$ , and let $\{{\mathbf{\varphi }}_{2b},{\mathbf{\varphi }}_{3b},b=1,\dots ,B\}$ denote a Gibbs sample of $({\mathbf{\varphi }}_2,{\mathbf{\varphi }}_3)$ .

$$\begin{array}{cc}{\widehat{\mathrm{CPO}}}_{\mathrm{i},\mathit{Long}}\hfill & =\widehat{{\mathit{CPO}}_{\mathit{i}}}f({\mathit{y}}_i|{\mathbf{\varphi }}_1^{*},{\mathit{x}}_i)\frac{1}{B}\sum _{b=1}^B{f}^{-1}({\mathit{y}}_i,t_{\mathit{Si}},{\delta }_i,t_i|{\mathbf{\varphi }}_1^{*},{\mathbf{\varphi }}_{2b},{\mathbf{\varphi }}_{3b},{\mathit{x}}_i,{\mathit{z}}_i(t),{\nu }_i).\hfill \end{array}$$

Next, fix ${\mathbf{\varphi }}_1={\mathbf{\varphi }}_1^{*}$ and ${\mathbf{\varphi }}_2={\mathbf{\varphi }}_2^{*}$ , and let $\{{\mathbf{\varphi }}_{3b},b=1,\dots ,B\}$ denote a Gibbs sample of $({\mathbf{\varphi }}_3)$ . We can calculate ${\sum }_{b=1}^B{f}^{-1}({\mathit{y}}_i,t_{\mathit{Si}},{\delta }_i,t_i|{\mathbf{\varphi }}_1^{*},{\mathbf{\varphi }}_2^{*},{\mathbf{\varphi }}_{3b},{\mathit{x}}_i,{\mathit{z}}_i(t),{\nu }_i)$ . Then we get

$$\begin{array}{cc}{\widehat{\mathrm{CPO}}}_{\mathrm{i},\mathit{Pg}|\mathit{Long}}\hfill & =\frac{f(t_{\mathit{Si}},{\delta }_i,{\mathit{y}}_i|{\mathbf{\varphi }}_1^{*},{\mathbf{\varphi }}_2^{*},{\mathit{x}}_i)}{f({\mathit{y}}_i|{\mathbf{\varphi }}_1^{*},{\mathit{x}}_i)}\frac{{\sum }_{b=1}^B{f}^{-1}({\mathit{y}}_i,t_{\mathit{Si}},{\delta }_i,t_i|{\mathbf{\varphi }}_1^{*},{\mathbf{\varphi }}_2^{*},{\mathbf{\varphi }}_{3b},{\mathit{x}}_i,{\mathit{z}}_i(t),{\nu }_i)}{{\sum }_{b=1}^B{f}^{-1}({\mathit{y}}_i,t_{\mathit{Si}},{\delta }_i,t_i|{\mathbf{\varphi }}_1^{*},{\mathbf{\varphi }}_{2b},{\mathbf{\varphi }}_{3b},{\mathit{x}}_i,{\mathit{z}}_i(t),{\nu }_i)}.\hfill \end{array}$$

Finally, ${\widehat{\mathrm{CPO}}}_{\mathrm{i},\mathit{Surv}|\mathit{Pg},\mathit{Long}}=\frac{{\widehat{\mathrm{CPO}}}_{\mathrm{i}}}{{\widehat{\mathrm{CPO}}}_{\mathrm{i},\mathit{Long}}{\widehat{\mathrm{CPO}}}_{\mathrm{i},\mathit{Pg}|\mathit{Long}}}$ .

Algorithm II. Following Section S.2 of the supplementary materials, Let $w_i({\mathbf{\theta }}_i^{*})$ be a normalized weight function such that $\int w_i({\mathbf{\theta }}_i^{*})d{\mathbf{\theta }}_i^{*}=1$ . $\widehat{{{\mathit{CPO}}_{\mathit{i}}}^{-1}}=\int w_i({\mathbf{\theta }}_i^{*})f^{-1}({\mathit{y}}_i,t_{\mathit{Si}},{\delta }_i,t_i,{\mathbf{\theta }}_i^{*}|\mathbf{\varphi },{\mathit{x}}_i,{\mathit{z}}_i(t),{\nu }_i)\pi (\mathbf{\varphi },{\mathbf{\theta }}^{*R}|D_{\mathit{obs}})d{\mathbf{\theta }}^{*R}d\mathbf{\varphi }$ . The opti- mal choice for $w_i({\mathbf{\theta }}_i^{*})$ in minimizing the variance of the Monte Carlo estimator for ${\widehat{\mathit{CPO}}}_{\mathrm{i}}$ is $\frac{f({\mathit{y}}_i,t_{\mathit{Si}},{\delta }_i,t_i,{\mathbf{\theta }}_i^{*}|\mathbf{\varphi },{\mathit{x}}_i,{\mathit{z}}_i(t),{\nu }_i)}{f({\mathit{y}}_i,t_{\mathit{Si}},{\delta }_i,t_i|\mathbf{\varphi },{\mathit{x}}_i,{\mathit{z}}_i(t),{\nu }_i)}$ . Now, let $\{{\mathbf{\varphi }}_b,{\mathbf{\theta }}_b^{*R},b=1,\dots ,B\}$ denote a Gibbs sample of $(\mathbf{\varphi },{\mathbf{\theta }}^{*R})$ . A estimate of ${\mathit{CPO}}_{\mathit{i}}^{-1}$ is given by

$${{\widehat{\mathrm{CPO}}}_{\mathrm{i}}}^{-1}=\frac{1}{B}\sum _{b=1}^B{w}_i({\mathbf{\theta }}_{\mathit{ib}}^{*})f^{-1}({\mathit{y}}_i,t_{\mathit{Pi}},{\delta }_i,t_i,{\mathbf{\theta }}_{\mathit{ib}}^{*}|{\mathbf{\varphi }}_b,{\mathit{x}}_i,{\mathit{z}}_i(t),{\nu }_i).$$

(4.10)

Similarly, $\int f^{-1}({\mathit{y}}_i,t_{\mathit{Pi}},{\delta }_i,t_i|{\mathbf{\varphi }}_1^{*},{\mathbf{\varphi }}_2^{*},{\mathbf{\varphi }}_3,{\mathit{x}}_i,{\mathit{z}}_i(t),{\nu }_i)\pi ({\mathbf{\varphi }}_3|{\mathbf{\varphi }}_1^{*},{\mathbf{\varphi }}_2^{*},D_{\mathit{obs}})d{\mathbf{\varphi }}_3=\int w_i({\mathbf{\theta }}_i^{*})f^{-1}({\mathit{y}}_i,t_{\mathit{Pi}},{\delta }_i,t_i,{\mathbf{\theta }}_i^{*}|{\mathbf{\varphi }}_1^{*},{\mathbf{\varphi }}_2^{*},{\mathbf{\varphi }}_3,{\mathit{x}}_i,{\mathit{z}}_i(t),{\nu }_i)\pi ({\mathbf{\varphi }}_3,{\mathbf{\theta }}^{*R}|{\mathbf{\varphi }}_1^{*},{\mathbf{\varphi }}_2^{*},D_{\mathit{obs}})d{\mathbf{\theta }}^{*R}d{\mathbf{\varphi }}_3$ . The optimal choice for $w_i({\mathbf{\theta }}_i^{*})$ in minimizing the variance of the Monte Carlo estimator is $f({\mathit{y}}_i,t_{\mathit{Pi}},{\delta }_i,t_i,{\mathbf{\theta }}_i^{*}|{\mathbf{\varphi }}_1^{*},{\mathbf{\varphi }}_2^{*},{\mathbf{\varphi }}_3,{\mathit{x}}_i,{\mathit{z}}_i(t),{\nu }_i)f^{-1}({\mathit{y}}_i,t_{\mathit{Pi}},{\delta }_i,t_i|{\mathbf{\varphi }}_1^{*},{\mathbf{\varphi }}_2^{*},{\mathbf{\varphi }}_3,{\mathit{x}}_i,{\mathit{z}}_i(t),{\nu }_i)$ . $\int f^{-1}({\mathit{y}}_i,t_{\mathit{Pi}},{\delta }_i,t_i|{\mathbf{\varphi }}_1^{*},{\mathbf{\varphi }}_2,{\mathbf{\varphi }}_3,{\mathit{x}}_i,{\mathit{z}}_i(t),{\nu }_i)\pi ({\mathbf{\varphi }}_2,{\mathbf{\varphi }}_3|{\mathbf{\varphi }}_1^{*},D_{\mathit{obs}})d{\mathbf{\varphi }}_{2}d{\mathbf{\varphi }}_3=\int w_i({\mathbf{\theta }}_i^{*})f^{-1}({\mathit{y}}_i,t_{\mathit{Pi}},{\delta }_i,t_i,{\mathbf{\theta }}_i^{*}|{\mathbf{\varphi }}_1^{*},{\mathbf{\varphi }}_2,{\mathbf{\varphi }}_3,{\mathit{x}}_i,{\mathit{z}}_i(t),{\nu }_i)\pi ({\mathbf{\varphi }}_2,{\mathbf{\varphi }}_3,{\mathbf{\theta }}^{*R}|{\mathbf{\varphi }}_1^{*},D_{\mathit{obs}})d{\mathbf{\theta }}^{*R}d{\mathbf{\varphi }}_{2}d{\mathbf{\varphi }}_3$ . The optimal choice for $w_i({\mathbf{\theta }}_i^{*})$ in minimizing the variance of the Monte Carlo estimator is $f({\mathit{y}}_i,t_{\mathit{Pi}},{\delta }_i,t_i,{\mathbf{\theta }}_i^{*}|{\mathbf{\varphi }}_1^{*},{\mathbf{\varphi }}_2,{\mathbf{\varphi }}_3,{\mathit{x}}_i,{\mathit{z}}_i(t),{\nu }_i)f^{-1}({\mathit{y}}_i,t_{\mathit{Pi}},{\delta }_i,t_i|{\mathbf{\varphi }}_1^{*},{\mathbf{\varphi }}_2,{\mathbf{\varphi }}_3,{\mathit{x}}_i,{\mathit{z}}_i(t),{\nu }_i)$ .

Next, fix ${\mathbf{\varphi }}_1={\mathbf{\varphi }}_1^{*}$ , and let $\{{\mathbf{\varphi }}_{2b},{\mathbf{\varphi }}_{3b},{\mathbf{\theta }}_b^{*R},b=1,\dots ,B\}$ denote a Gibbs sample of $({\mathbf{\varphi }}_2,{\mathbf{\varphi }}_3,{\mathbf{\theta }}^{*R})$ . Set the weight function to be $w_{i1}({\mathbf{\theta }}_i^{*})$ .

$$\begin{array}{cc}{\widehat{\mathrm{CPO}}}_{\mathrm{i},\mathit{Long}}=\hfill & \widehat{{\mathit{CPO}}_{\mathit{i}}}f({\mathit{y}}_i|{\mathbf{\varphi }}_1^{*},{\mathit{x}}_i)\hfill \\ \hfill & \times \frac{1}{B}\sum _{b=1}^B{w}_{i1}({\mathbf{\theta }}_{\mathit{ib}}^{*})f^{-1}({\mathit{y}}_i,t_{\mathit{Si}},{\delta }_i,t_i,{\mathbf{\theta }}_{\mathit{ib}}^{*}|{\mathbf{\varphi }}_1^{*},{\mathbf{\varphi }}_{2b},{\mathbf{\varphi }}_{3b},{\mathit{x}}_i,{\mathit{z}}_i(t),{\nu }_i).\hfill \end{array}$$

Next, fix ${\mathbf{\varphi }}_1={\mathbf{\varphi }}_1^{*}$ and ${\mathbf{\varphi }}_2={\mathbf{\varphi }}_2^{*}$ , and let $\{{\mathbf{\varphi }}_{3b},{\mathbf{\theta }}_b^{*R},b=1,\dots ,B\}$ denote a Gibbs sample of $({\mathbf{\varphi }}_3,{\mathbf{\theta }}^{*R})$ . Set the weight function to be $w_{i2}({\mathbf{\theta }}_i^{*})$ . We can calculate ${\sum }_{b=1}^B{w}_{i2}({\mathbf{\theta }}_{\mathit{ib}}^{*})f^{-1}({\mathit{y}}_i,t_{\mathit{Si}},{\delta }_i,t_i,{\mathbf{\theta }}_{\mathit{ib}}^{*}|{\mathbf{\varphi }}_1^{*},{\mathbf{\varphi }}_2^{*},{\mathbf{\varphi }}_{3b},{\mathit{x}}_i,{\mathit{z}}_i(t),{\nu }_i)$ . Then we get

$$\begin{array}{cc}{\widehat{\mathrm{CPO}}}_{\mathrm{i},\mathit{Pg}|\mathit{Long}}=\hfill & \frac{f(t_{\mathit{Si}},{\delta }_i,{\mathit{y}}_i|{\mathbf{\varphi }}_1^{*},{\mathbf{\varphi }}_2^{*},{\mathit{x}}_i)}{f({\mathit{y}}_i|{\mathbf{\varphi }}_1^{*},{\mathit{x}}_i)}\hfill \\ \hfill & \times \frac{{\sum }_{b=1}^B{w}_{i2}({\mathbf{\theta }}_{\mathit{ib}}^{*})f^{-1}({\mathit{y}}_i,t_{\mathit{Si}},{\delta }_i,t_i,{\mathbf{\theta }}_{\mathit{ib}}^{*}|{\mathbf{\varphi }}_1^{*},{\mathbf{\varphi }}_2^{*},{\mathbf{\varphi }}_{3b},{\mathit{x}}_i,{\mathit{z}}_i(t),{\nu }_i)}{{\sum }_{b=1}^B{w}_{i1}({\mathbf{\theta }}_{\mathit{ib}}^{*})f^{-1}({\mathit{y}}_i,t_{\mathit{Si}},{\delta }_i,t_i,{\mathbf{\theta }}_{\mathit{ib}}^{*}|{\mathbf{\varphi }}_1^{*},{\mathbf{\varphi }}_{2b},{\mathbf{\varphi }}_{3b},{\mathit{x}}_i,{\mathit{z}}_i(t),{\nu }_i)}.\hfill \end{array}$$

Finally, ${\widehat{\mathrm{CPO}}}_{\mathrm{i},\mathit{Surv}|\mathit{Pg},\mathit{Long}}=\frac{{\widehat{\mathrm{CPO}}}_{\mathrm{i}}}{{\widehat{\mathrm{CPO}}}_{\mathrm{i},\mathit{Long}}{\widehat{\mathrm{CPO}}}_{\mathrm{i},\mathit{Pg}|\mathit{Long}}}$ .

Due to the computational difficulty of the optimal weight functions, another possible choice is $w_{i,\mathit{cond}}({\mathbf{\theta }}_i^{*})=f({\mathbf{\theta }}_i^{*}|{\mathbf{\varphi }}_1,{\mathit{x}}_i,{\mathit{y}}_i)$ . The conditional distribution of ${\mathbf{\theta }}_i^{*}|{\mathbf{\varphi }}_1,{\mathit{x}}_i,{\mathit{y}}_i$ is ${\mathbf{\theta }}_i^{*}|{\mathbf{\varphi }}_1,{\mathit{x}}_i,{\mathit{y}}_i\sim N({\mathbf{\mu }}^{*},{\Sigma }^{*})$ , where ${\mathbf{\mu }}^{*}={\Gamma }^{\prime }{W}_{i1}^{\prime }({\sigma }^2{I}_{m_i}+W_{i1}\Gamma {\Gamma }^{\prime }{W}_{i1}^{\prime }{)}^{-1}({\mathit{y}}_i-W_i({\mathbf{\theta }}^{\prime },{\mathbf{\gamma }}^{\prime }{)}^{\prime })$ and ${\Sigma }^{*}=I_p-{\Gamma }^{\prime }{W}_{i1}^{\prime }({\sigma }^2{I}_{m_i}+W_{i1}\Gamma {\Gamma }^{\prime }{W}_{i1}^{\prime }{)}^{-1}{W}_{i1}\Gamma$ .

5 Analysis of the head and neck cancer data

We apply our proposed method in Section 3 to analyse the head and neck cancer data introduced in Section 2. Out of the total 483 patients, we exclude two patients from the analysis because of missing BSDTARG0 measures. We are interested in identifying which of the three patient-reported measures: GH, HNPA and HNSW contributes most to the model fit of the survival data. Each of GH, HNPA and HNSW was measured on a continuous scale from 0 to 100. Following Ibrahim et al. (2001), we apply transformation to the original longitudinal measures: $\sqrt{100-\mathit{GH}}$ , $\sqrt{100-\mathit{HNPA}}$ and $\sqrt{100-\mathit{HNSW}}$ , respectively, to achieve a better normal approximation. For each transformed longitudinal measure, we fit the joint model for the longitudinal measure, disease progression time and the overall survival, with a time-dependent covariate accounting for the model fit of the survival data. We include 11 covariates described in Section 2, namely ECOG, AGE, PEGFRA, FEMALE PCRTDC, TMSITED01, METLOC2, METLOC3, METLOC4, ETLOC5, BSDTARG, and treatment indicator R in the time to progression model (3.8) and the overall survival model (3.10). We include five covariates, namely ECOG, AGE, PEGFRA, FEMALE, and PCRTDC, and treatment indicator R in the cure rate model (3.6). In the time to progression model the survival model, and the cure rate model, all of the covariates including those binary endpoints, except for the treatment indicator R, are standardized to have mean 0 and standard division 1. We use a linear trajectory in the longitudinal model. No standardization is done for any of these covariates in the longitudinal model. In the overall survival model, we also include time-dependent indicator for treatment switch ${\mathit{V}}_{{\mathit{T}}_{\mathit{SW}}}$ , time-dependent indicator for progression status ${\mathit{V}}_{{\mathit{T}}_{\mathit{P}}}$ and an interaction term $\mathit{R}*{\mathit{V}}_{{\mathit{T}}_{\mathit{P}}}$ .

We used DIC and LPML under the survival model alone for overall survival and the cure rate model alone for time to progression to determine $K_j$ , the number of pieces for the piecewise constant hazard function for $j=1,2$ . As reported in the footnotes under Tables S.5 and S.6 in the supplementary materials, the ‘best choices’ were $K_1=25$ for time to progression and $K_2=7$ for overall survival according to both DIC and LPML. Under these choices of $K_1$ and $K_2$ , we fit the models of the survival data jointly with each of GH, HNPA and HNSW longitudinal data and calculated the decompositions of DIC and LPML. The values of DIC, $\Delta {\mathit{DIC}}_{\mathit{Pg}}$ , $\Delta {\mathit{DIC}}_{\mathit{Surv},1}$ $\Delta {\mathit{DIC}}_{\mathit{Surv},2}$ , LPML, $\Delta {\mathit{LPML}}_{\mathit{Pg}}$ , $\Delta {\mathit{LPML}}_{\mathit{Surv},1}$ , and $\Delta {\mathit{LPML}}_{\mathit{Surv},2}$ are shown in Figure 2 for GH and reported in Table 2 for all three longitudinal measures.

Figure 2:

Diagram of the decompositions of DIC and LPML for the disease progression and survival data jointly with GH

Larger values of $\Delta {\mathit{DIC}}_{\mathit{Surv}}$ and $\Delta {\mathit{LPML}}_{\mathit{Surv}}$ correspond to more gain in the fit of the overall survival data by adding time to disease progression and different longitudinal components. All values of $\Delta {\mathit{DIC}}_{\mathit{Surv},1}$ and $\Delta {\mathit{LPML}}_{\mathit{Surv},1}$ in Table 2 were positive, implying that all three longitudinal measures together with the time-dependent covariates helped in the fit of the overall survival data. All values of $\Delta {\mathit{DIC}}_{\mathit{Surv},2}$ and $\Delta {\mathit{LPML}}_{\mathit{Surv},2}$ were positive, indicating that all three longitudinal measures helped in the fit of the overall survival data in addition to the time-dependent covariates. Comparing the magnitudes of $\Delta {\mathit{DIC}}_{\mathit{Surv},2}$ and $\Delta {\mathit{LPML}}_{\mathit{Surv},2}$ in Table 2, in terms of gain in the fit of overall survival, GH is more associated with the overall survival time than HNSW and HNPA, while HNSW has a stronger association with the overall survival time than HNPA. Similarly, larger values of $\Delta {\mathit{DIC}}_{\mathit{Pg}}$ and $\Delta {\mathit{LPML}}_{\mathit{Pg}}$ correspond to more gain in the fit of disease progression data by adding different longitudinal components. The results in Table 2 show that GH is helpful in the fit of the time to progression data, while there were weak associations between HNPA or HNSW and the time to progression data. After comparing the values of overall DIC for the models jointly with different longitudinal measures in Table 2, we see that HNPA $<$ GH $<$ HNSW in terms of the corresponding DIC values, while based on ${\mathit{DIC}}_{\mathit{Pg}|\mathit{Long}}$ and ${\mathit{DIC}}_{\mathit{Pg}|\mathit{Long}}$ , GH has the smallest values. The reason for this is that there were different numbers of observations (1 589 for GH, 1 587 for HNPA and 1 525 for HNSW) for the three longitudinal measures. Thus, the values of the overall DIC were not comparable. The similar results of LPML and LPML decompositions were observed.

Table 2:

The Decompositions of DIC and LPML for three longitudinal measures

	GH	HNPA	HNSW		GH	HNPA	HNSW
DIC	9 348.07	9 226.44	9 430.50	LPML	$-$ 4 687.24	$-$ 4 631.66	$-$ 4 805.36
${\mathit{DIC}}_{\mathit{Long}}$	6 111.25	5 964.54	6 336.95	${\mathit{LPML}}_{\mathit{Long}}$	$-$ 3 060.59	$-$ 2 991.56	$-$ 3 173.07
${\mathit{DIC}}_{\mathit{Pg}|\mathit{Long}}$	1 090.54	1 096.82	1 098.59	${\mathit{LPML}}_{\mathit{Pg}|\mathit{Long}}$	$-$ 552.44	$-$ 556.68	$-$ 555.80
${\mathit{DIC}}_{\mathit{Surv}|\mathit{Pg},\mathit{Long}}$	2 146.27	2 165.09	2 151.79	${\mathit{LPML}}_{\mathit{Surv}|\mathit{Pg},\mathit{Long}}$	$-$ 1 074.21	$-$ 1 083.42	$-$ 1 076.50
$\Delta {\mathit{DIC}}_{\mathit{Surv},1}$	78.53	59.71	73.00	$\Delta {\mathit{LPML}}_{\mathit{Surv},1}$	38.93	29.72	36.64
$\Delta {\mathit{DIC}}_{\mathit{Surv},2}$	34.75	15.93	29.23	$\Delta {\mathit{LPML}}_{\mathit{Surv},2}$	17.17	7.96	14.88
$\Delta {\mathit{DIC}}_{\mathit{Pg}}$	7.48	1.21	$-$ 0.56	$\Delta {\mathit{LPML}}_{\mathit{Pg}}$	2.57	$-$ 1.66	-0.78

Source: ${\mathit{DIC}}_{\mathit{Surv}0,1}=2224.80$ , ${\mathit{DIC}}_{\mathit{Surv}0,2}=2181.02$ , ${\mathit{DIC}}_{\mathit{Pg}0}=1098.03$ , ${\mathit{LPML}}_{\mathit{Surv}0,1}=-1113.14$ , ${\mathit{LPML}}_{\mathit{Surv}0,2}=-1091.38$ , ${\mathit{LPML}}_{\mathit{Pg}0}=-555.01$

Table 3 shows the posterior estimates and 95% HPD intervals of the hazard ratios of variables/parameters for the overall survival model (3.10) and the time to progression model (3.8), and the odds ratios of variables for the cure rate logistic regression model (3.6) and the posterior estimates and the 95% HPD intervals of parameters for the longitudinal model (3.1) under the joint model with GH. Posterior means, posterior standard deviations, and 95% HPD intervals of the parameters under the joint models with GH, HNPA and HNSW are reported in Tables S.2 to S.4 in the supplementary materials. We see from Table S.2 that both ${\beta }_1$ and ${\beta }_2$ with 95% HPD intervals (0.052, 0.246) and (0.195, 0.393), respectively, were significant, in the sense that their corresponding 95% HPD intervals do not contain zero. From Table 3, we also see that the 95% HPD intervals of the hazard ratios corresponding to ${\beta }_1$ and ${\beta }_2$ are (1.053, 1.279) and (1.215, 1.482), respectively, which do not contain 1. Thus, a covariate is significant if its 95% HPD interval of the hazard ratio does not contain 1. These results imply that an increase in the time-dependent longitudinal trajectory function is associated with a shorter disease progression time as well as a worse overall survival mortality, which is clinically appealing since a large value of the longitudinal trajectory function corresponds to a small GH value. Under the proposed joint model, ${\beta }_1$ constitutes relationship between the longitudinal measure and the time to progression, while ${\beta }_2$ induces correlation between the longitudinal measure and the overall survival time. In particular, ${\beta }_1=0$ or ${\beta }_2=0$ implies that the longitudinal measure does not influence the risk of disease progression or overall survival. From Tables S.3 and S.4, we see that only ${\beta }_2$ was significant for both HNPA with respective 95% HPD intervals ( $-$ 0.237, $-$ 0.087) and ( $-$ 0.193, $-$ 0.095). These results were consistent with our decompositions of DIC and LPML as the values of $\Delta {\mathit{DIC}}_{\mathit{Surv},2}$ and $\Delta {\mathit{LPML}}_{\mathit{Surv},2}$ were 34.75 and 17.17 for GH, 15.93 and 7.96 for HNPA, and 29.23 and 14.88 for HNSW, while the values of $\Delta {\mathit{DIC}}_{\mathit{Pg}}$ and $\Delta {\mathit{LPML}}_{\mathit{Pg}}$ were 7.48 and 2.57 for GH, 1.21 and $-$ 1.66 for HNPA, and $-$ 0.56 and $-$ 0.78 for HNSW. Similarly, from Tables S.2 to S.4, we see that the 95%HPD for the coefficient of time-dependent covariate ${\mathit{V}}_{{\mathit{T}}_{\mathit{P}}}$ was significant for all of GH, HNPA and HNSW, which is consistent with the results shown in Table 2, in which $\Delta {\mathit{DIC}}_{\mathit{Surv},1}>\Delta {\mathit{DIC}}_{\mathit{Surv},2}$ and $\Delta {\mathit{LPML}}_{\mathit{Surv},1}>\Delta {\mathit{LPML}}_{\mathit{Surv},2}$ under all of GH, HNPA and HNSW. Interestingly, the longitudinal measures could have a minor impact on the significance of covariates under the model for the overall survival. Specifically, METLOC5 was not significant under the joint model with GH, while METLOC5 was significant under the joint models with HNPA and HNSW. Also, the treatment effect (R) without treatment switch was significant under the joint models with GH and HNSW, while it was not significant under the joint model with HNPA. Finally, we note that (a) the treatment effect $R$ was not significant under the cure rate model alone for time to progression (Table S.5); (b) the survival model alone for overall survival without time-dependent covariates (Table S.6); and (c) the treatment effect (R) without treatment switch was significant under the survival model alone for overall survival with time-dependent covariates (Table S.7).

Table 3:

Parameter estimates under the joint model with GH

Parameter	EST	95%HPD	Parameter	EST	95%HPD
Survival (HR)			Progression (HR)
ECOG	2.221	(1.669, 2.949)	ECOG	1.062	(0.780, 1.442)
AGE	0.937	(0.707, 1.207)	AGE	0.833	(0.577, 1.190)
PEGFRAD	1.174	(0.916, 1.500)	PEGFRAD	1.548	(1.170, 2.043)
Female	0.992	(0.707, 1.375)	Female	0.964	(0.642, 1.446)
PCRTDC	1.242	(0.978, 1.566)	PCRTDC	1.028	(0.750, 1.349)
TMSITED01	1.507	(1.150, 1.947)	TMSITED01	0.853	(0.624, 1.123)
METLOC2	1.072	(0.797, 1.408)	METLOC2	0.942	(0.686, 1.229)
METLOC3	0.755	(0.575, 0.949)	METLOC3	1.046	(0.800, 1.389)
METLOC4	1.699	(1.137, 2.499)	METLOC4	1.983	(1.225, 3.007)
METLOC5	1.397	(0.910, 2.055)	METLOC5	1.127	(0.696, 1.758)
BSDTARG	1.009	(1.006, 1.012)	BSDTARG	1.003	(1.000, 1.006)
R	0.705	(0.501, 0.945)	R	0.973	(0.723, 1.288)
${\mathit{V}}_{{\mathit{T}}_{\mathit{P}}}$	2.490	(1.904, 3.283)	${\beta }_1$	1.161	(1.053, 1.279)
${\mathit{V}}_{{\mathit{T}}_{\mathit{SW}}}$	0.782	(0.515, 1.129)	Longitudinal
$\mathit{R}*{\mathit{V}}_{{\mathit{T}}_{\mathit{SW}}}$	1.410	(0.861, 2.118)	${\theta }_1$	4.863	(4.367, 5.397)
${\beta }_2$	1.346	(1.215, 1.482)	${\theta }_2$	0.126	(0.080, 0.173)
Logistic (OR)			ECOG	1.285	(0.940, 1.626)
ECOG	0.893	(0.133, 7.473)	AGE	0.187	(-0.144, 0.539)
AGE	0.569	(0.104, 16.586)	PEGFRAD	0.263	(-0.021, 0.578)
PEGFRAD	0.679	(0.056, 4.709)	Female	0.113	(-0.326, 0.534)
Female	1.270	(0.114, 95.950)	PCRTDC	0.102	(-0.189, 0.415)
PCRTDC	2.715	(0.447, 24.536)	TMSITED01	0.204	(-0.131, 0.520)
R	1.796	(0.457, 9.381)	METLOC2	0.169	(-0.196, 0.537)
			METLOC3	-0.028	(-0.342, 0.267)
			METLOC4	-0.604	(-1.182, -0.093)
			METLOC5	0.339	(-0.197, 0.898)
			BSDTARG	0.002	(-0.002, 0.006)
			R	-0.029	(-0.331, 0.295)
			${\Omega }_{11}$	2.036	(1.644, 2.463)
			${\Omega }_{21}$	-0.021	(-0.097, 0.059)
			${\Omega }_{22}$	0.028	(0.000, 0.057)

In all the Bayesian computations, we used 5 000 MCMC samples, which were taken from every third iteration, after a burn-in of 1 500 iterations for each model to compute all posterior estimates, including posterior means, posterior SDs, 95% HPD intervals, DIC and LPML. The convergence of the MCMC sampling algorithm was checked using several diagnostic procedures discussed in (Chen et al., 2000). The trace and autocorrelation plots of MCMC chains for the model parameters are shown in Figures S.1 to S.18. We see from these figures that the autocorrelations disappear even after lag 2 for some parameters and after lag 20 for most of the parameters. These trace and autocorrelation plots show good convergence and mixing of these MCMC chains. The HPD intervals were computed via the Monte Carlo method developed by (Chen and Shao, 1999). Computer code was written for the FORTRAN 95 compiler, and we used IMSL subroutines with double precision accuracy. The FORTRAN code is available from the authors upon request. In addition, an R interface of our Fortran code is developed, which can be downloaded from the journal website (http://www.statmod.org/smij/archive.html).

6 Discussion

In this article, we proposed a class of mixed effects regression models for longitudinal measures, a cure rate model for the disease progression time and a time-varying covariates model for the overall survival. In Section 5, we carried out a Bayesian analysis of the head and neck cancer data using the computational algorithm discussed in Section 4.2. To examine the empirical performance of the posterior estimates under the proposed joint model, we conducted a simulation study. The design and results of the simulation study are presented in Section S.3 of the supplementary materials. As shown in Table S.1, the posterior estimates are close to the true values for almost all of the parameters except for $\mathbf{\varphi }$ , and the CP's are close to 95% for most of the parameters. The simulation results empirically confirm convergence of the proposed Markov chain Monte Carlo sampling algorithm and identifiability of the model parameters except for $\mathbf{\varphi }$ under the proposed joint model.

Within the semi-competing risk setting, we developed the decompositions of DIC and LPML to determine which PRO is most associated with the survival outcome. As discussed in Section 4.3, the DIC and DIC decompositions require a numerical approximation of an integral over random effects ${\theta }^{*}$ , while the LPML and LPML decompositions can be computed using the Monte Carlo methods. Thus, the LPML and LPML decompositions are more advantageous when the dimension of random effects becomes large. The proposed joint models and the corresponding decompositions of Bayesian model selection criteria can be extended to discrete longitudinal biomarkers such as Eastern Cooperative Oncology Group (ECOG) performance status and other non-normal longitudinal biomarkers. Other extensions include the joint models of the multiple longitudinal measures and the survival data with semi-competing risks. These extensions are currently under investigation.