Bayesian Joint Modeling Analysis of Longitudinal Proportion and Survival Data

1. Introduction

Joint analysis of longitudinal and survival data have been widely applied to cancer and HIV/AIDS clinical studies where time to event outcomes, such as disease free and overall survivals, are usually the primary endpoints. In particular, after Faucett and Thomas [1] and Wulfsohn and Tsiatis [2] introduced what is the standard joint model, the joint model has been extensively explored. The advantages of joint models have been discussed by researchers [3,4,5,6,7,8]. However, in clinical trials, some patients with poor quality of life may drop out of the study due to disease recurrence or death. In this case, the absence of quality of life measures resulting from the withdrawal of these patients is informative. In order to provide strong evidence for the tradeoff between intensive treatment and the associated poor quality of life for patients, we conducted joint modeling of longitudinal life measures and survival data to investigate the relationship between them. For longitudinal quality of life and survival data, Henderson et al. [9] and Zeng and Cai [10] considered the use of shared normal distribution random effects to jointly analyze the relationship between longitudinal quality of life and survival time. Tang et al. [11] considered a novel semiparametric joint model for multivariate longitudinal and survival data to analyze data from the International Breast Cancer Study. Longitudinal life quality measurement data can be linearly converted into longitudinal proportion data whose value range is in the unit interval (0,1) [12]. For the longitudinal component, there are two methods to deal with it. The first method applied the classic linear mixed model to the longitudinal scale data after logit transformation, and the second method directly used the simplex distribution to model the longitudinal scale data. The models established by the two methods both used the EM algorithm and the Laplace approximation to estimate the unknown parameters. In order to be more flexible and practical, this paper will use a partial linear mixed effect model to the logit transformed longitudinal scale data, and use the B-splines method to model the nonlinear function in the model. Meanwhile, to promote the feasibility of our proposed model, we use CDPMM method to model random effects.

In addition, variables selection in the joint model is also considered. In traditional regression models, methods of variable selection include the forward selection method, backward elimination method, stepwise selection method or some information criterion such as the Akaike information criterion. But these approaches are computationally expensive and unstable for the complicated models with a large number of covariates. Therefore, some penalized likelihood methods are gradually proposed, especially Lasso of Tibshirani [13] estimating linear regression coefficients by $L_1$-constrained least squares. Tibshirani [13] proposed that Lasso estimates can be interpreted as posterior norm estimates when the regression parameters have independent and identical Laplacian priors (ie, double exponential priors). Motivated by this connection, under the Bayesian framework, Park and Casella [14] proposed the BLasso variable selection method by imposing the double exponential prior on the regression coefficients and the gamma distribution on the shrinkage parameter. The BLasso method has been applied to linear regression model [15], semi-parametric structural equation model [16] and joint models of longitudinal and survival data [17]. Therefore, this paper extended that BLasso variable selection method to the joint model of longitudinal proportion data and survival data and a BLasso approach is developed to simultaneously estimate unknown parameters and identify the significant effect of the important covariates.

The rest of this paper is organized as follows. In Section 2, the joint model of longitudinal proportion and survival data is introduced. In Section 3, the Bayesian estimations of joint model are proposed. In Section 4, three numerical simulations are presented to evaluate the performance of the proposed methods. In Section 5, we apply the proposed approach to the MA.5 research experiment data and some concluding remarks are given in Section 6. Technical details are given in Appendix A.

2. Model and Notation

Consider a data set from n individuals. Let $y_{ij}$ be a longitudinal proportional measurement for the i-th individual ($i=1,2,\dots ,n$) at observation time point $t_{ij}$ for $j=1,2,\dots ,n_i$, and $y_{ij}\in (0,1)$. We assume $y_{ij}^{*}$ is the logit-transformation of $y_{ij}$, and $y_{ij}^{*}=\mathrm{logit}\left(y_{ij}\right)=log\left(\frac{1-y_{ij}}{y_{ij}}\right)$. Further, $T_i^{*}$ and $C_i$ are the true survival time and censoring time, respectively. Let $T_i=min(T_i^{*},C_i)$ denotes the corresponding observed event time, and ${\delta }_i=\mathit{1}(T_i^{*}\le C_i)$ denotes the failure indicator, where $\mathit{1}(·)$ is a indicator function.

Denote ${\mathit{y}}^{*}=\left\{{\mathit{y}}_1^{*},{\mathit{y}}_2^{*},\dots ,{\mathit{y}}_n^{*}\right\}$, ${\mathit{y}}_i^{*}=\left\{y_{i1}^{*},y_{i2}^{*},\dots ,y_{i{n}_i}^{*}\right\}$, $\mathit{T}=\left\{T_1,T_2,\dots ,T_n\right\}$, $\Delta =\left\{{\delta }_1,{\delta }_2,\dots ,{\delta }_n\right\}$. Let $\mathit{b}=\{{\mathit{b}}_1,{\mathit{b}}_2,\dots ,{\mathit{b}}_n\}$ be time-independent random effects underlying both the longitudinal and survival processes for the i-th individual. Given the random effects ${\mathit{b}}_i$, assume that $y_{ij}^{*}$ follows partially linear mixed-effect model

$$\begin{array}{c}\hfill y_{ij}^{*}|{\mathit{b}}_i={\mathit{X}}_{ij}^{\top }\mathsf{\beta }+g\left(t_{ij}\right)+{\mathit{Z}}_{ij}^{\top }{\mathit{b}}_i+{\epsilon }_{ij},\end{array}$$

(1)

where ${\mathit{X}}_{ij}$ and ${\mathit{Z}}_{ij}$ represent the time-independent design vectors of fixed and random effects associated with $y_{ij}^{*}$, respectively. And $\mathsf{\beta }$ is a $p_1\times 1$ fixed effects regression parameters, ${\mathit{b}}_i$ is a $q\times 1$ random effects vector, $g\left(t\right)$ is a twice continuous differentiable unknown function, and ${\epsilon }_{ij}$ is a white noise process with variance ${\sigma }^2$. Also, we assume that ${\epsilon }_{ij}$’s are independent of ${\mathit{b}}_i$. To facilitate the feasibility of our proposed model, instead of the traditional normality assumption, which may be violated in some applications [18], it is specified the random effects by Dirichlet process (DP) mixture of normals.

For event time $T_i$, given random effects ${\mathit{b}}_i$, we assume that $T_i$ follows the hazard model

$$\begin{array}{c}\hfill {\lambda }_i\left(t\right|{\mathit{b}}_i)={\lambda }_0\left(t\right)exp({\mathit{W}}_i^{\top }\mathsf{\gamma }+{\mathsf{\varphi }}^{\top }{\mathit{b}}_i),\end{array}$$

(2)

where ${\mathit{W}}_i$ is the known fixed effects design matrix linking the unknown $p_2\times 1$ parameter vector $\mathsf{\gamma }$ to ${\lambda }_i\left(t\right|{\mathit{b}}_i)$, $\mathsf{\varphi }$ is an unknown $q\times 1$ parameter vector linking ${\mathit{b}}_i$ to ${\lambda }_i\left(t\right|{\mathit{b}}_i)$, and ${\lambda }_0\left(t\right)$ is an unknown basic hazard function.

From the above discussion, it is suggested to link models (1) and (2) through shared random effects, called a shared random effects joint model (JMSRE). Parameter $\mathsf{\varphi }$ in model JMSRE reflects the correlation between transformed longitudinal proportional data and survival data, given random effects. When $\mathsf{\varphi }=0_{q\times 1}$, it means that the longitudinal index is not necessarily related to the event time, that is, longitudinal proportional data and survival data can be modeled separately. So in this case joint modeling is not necessary and longitudinal indicators can be ignored for modeling survival data.

Further, to make Bayesian inference on $\mathsf{\beta }$ based on model (1), we approximate $g\left(t\right)$ through a B-splines method

$$\begin{array}{c}\hfill g\left(t\right)\approx B_1\left(t\right){\phi }_1+B_2\left(t\right){\phi }_2+\dots +B_L\left(t\right){\phi }_L={\mathit{B}}^{\top }\left(t\right)\mathsf{\phi },\end{array}$$

where $L=d+K+1$, d is the degree of B-splines, K is the number of knots, $\mathsf{\phi }={\left({\phi }_1,{\phi }_2,\dots ,{\phi }_L\right)}^{\top }$ is a $L\times 1$ unknown coefficient vector, and $\mathit{B}\left(t\right)={\left(B_1\left(t\right),B_2\left(t\right),\dots ,B_L\left(t\right)\right)}^{\top }$.

Denote ${\mathsf{\theta }}_y=\left\{\mathsf{\beta },\mathsf{\phi },{\sigma }^2\right\}$ be unknown parameters associated with model (1), and ${\mathsf{\theta }}_T=\left\{\mathsf{\gamma },\mathsf{\varphi }\right\}$ be unknown parameters associated with model (2). Thus, given $({\mathsf{\theta }}_y,{\mathsf{\theta }}_T,\mathit{b})$, a joint likelihood function of $({\mathit{y}}^{*},\mathit{T},\mathsf{\Delta })$ can be written as

$$p({\mathit{y}}^{*},\mathit{T},\mathsf{\Delta }|{\mathsf{\theta }}_y,{\mathsf{\theta }}_T,\mathit{b})=\prod _{i=1}^{n}p\left({\mathit{y}}_i^{*}|{\mathit{b}}_i;{\mathsf{\theta }}_y\right)p\left(\mathit{T},\mathsf{\Delta }|{\mathit{b}}_i;{\mathsf{\theta }}_T\right),$$

(3)

where

$$p\left({\mathit{y}}_i^{*}|{\mathit{b}}_i;{\mathsf{\theta }}_y\right)=\prod _{j=1}^{n_i}\frac{1}{\sqrt{2\pi {\sigma }^2}}exp\left\{-\frac{{\left(y_{ij}^{*}-{\mathit{X}}_{ij}^{\top }\mathsf{\beta }-{\mathit{B}}^{\top }\left(t_{ij}\right)\mathsf{\phi }-{\mathit{Z}}_{ij}^{\top }{\mathit{b}}_i\right)}^2}{2{\sigma }^2}\right\},$$

$$p\left(\mathit{T},\mathsf{\Delta }|{\mathit{b}}_i;{\mathsf{\theta }}_T\right)=\prod _{i=1}^n{\left[\frac{exp({\mathit{W}}_i^{\top }\mathsf{\gamma }+{\mathsf{\varphi }}^{\top }{\mathit{b}}_i)}{{\sum }_{j\in {\mathcal{R}}_i}exp({\mathit{W}}_j^{\top }\mathsf{\gamma }+{\mathsf{\varphi }}^{\top }{\mathit{b}}_j)}\right]}^{{\delta }_i},{\mathcal{R}}_i=\left\{j:T_j\ge T_i\right\}.$$

3. Bayesian estimation of joint model

3.1. Prior specification

To develop Bayesian inference on the considered models, we need specifying the prior distributions for ${\sigma }^2,\mathsf{\phi },\mathsf{\beta },\mathsf{\gamma }\mathrm{and}\mathsf{\varphi }$. For conjugation, we consider the following priors for ${\sigma }^2,\mathsf{\phi }$:

$${\sigma }^{-2}\sim \mathsf{\Gamma }\left(a_0,b_0\right),\mathsf{\phi }\sim N_L\left(\mathit{0},{\mathit{H}}_{\phi }^0\right),$$

(4)

where $a_0,b_0$, and ${\mathit{H}}_{\phi }^0$ are pre-given hyperparameters, $\Gamma \left(a_0,b_0\right)$ denotes the Gamma distribution with parameter $a_0$ and the shape parameter $b_0$.

Following Tang et al. [17], we specify the DP mixture of normals for the random effects ${\mathit{b}}_i$, that is, ${\mathit{b}}_i\stackrel{\mathrm{i}.\mathrm{i}.\mathrm{d}.}{\sim }{\sum }_{g=1}^{\infty }{\pi }_g{N}_q({\mathsf{\mu }}_g,{\mathsf{\Omega }}_g)$ with $({\mathsf{\mu }}_g,{\mathsf{\Omega }}_g)\sim \mathcal{P}$, where $\mathcal{P}$ is a random probability with an unknown form. one common approach is to specify DP prior to approximate $\mathcal{P}$, i.e., $\mathcal{P}\sim \mathrm{DP}\left(\tau F_0\right)$ in which $F_0$ is a base distribution that serves as a starting-point for constructing the nonparametric distribution, and $\tau$ is a weight that indicates the researcher’s certainty of $F_0$ as the distribution of $\mathcal{P}$. To gain the form of $\mathcal{F}$, motivated by Ohlssen et al. [18] and Yang et al. [19], we consider the following truncated approximate centered Dirichlet Process Mixture Model (CDPMM), which can specify the prior distribution of ${\mathit{b}}_i$ as follows:

$$b_i\stackrel{\mathrm{i}.\mathrm{i}.\mathrm{d}.}{\sim }\sum _{g=1}^G{\pi }_g{N}_q({\mathsf{\mu }}_g,{\mathsf{\Omega }}_g)\mathrm{with}{\mathsf{\mu }}_g={\mathsf{\mu }}_g^{*}-\sum _{g=1}^G{\pi }_g{\mathsf{\mu }}_g^{*}\mathrm{and}({\mathsf{\mu }}_g^{*},{\mathsf{\Omega }}_g)\stackrel{\mathrm{i}.\mathrm{i}.\mathrm{d}.}{\sim }{F}_0,$$

where $1\le G<\infty$, and ${\pi }_g$ is a random probability weight chosen to be independent of $({\mathsf{\mu }}_g^{*},{\mathsf{\Omega }}_g)$ such that $0\le {\pi }_g\le 1$ and ${\sum }_{g=1}^{\infty }{\pi }_g=1$. In order to ease the computational intensity, according to take moderate value proposed by Ishwaran and Zarepour (2000) [20], we consider $G=25$, and ${\pi }_g$ is given by the following stick-breaking procedure:

$${\pi }_1={\vartheta }_1\mathrm{and}{\pi }_g={\vartheta }_g\prod _{\iota =1}^{g-1}(1-{\vartheta }_{\iota })\mathrm{for}g=2,\dots ,G,$$

(5)

where ${\vartheta }_g\stackrel{\mathrm{i}.\mathrm{i}.\mathrm{d}.}{\sim }\mathrm{Beta}(1,\tau )$ for $g=1,2,\dots ,G-1$, and ${\vartheta }_G=1$ so that ${\sum }_{g=1}^G{\pi }_g=1$. The prior distribution for the unknown parameter $\tau$ is given by $\tau \sim \Gamma (a_1,a_2)$ with prespecified hyperparameters $a_1$ and $a_2$ [21]. Here, we took the hyperparameters $a_1$ and $a_2$ to be 25 and 5, respectively.

Sampling from posterior distributions of ${\mathit{b}}_i$ with the above specified DP prior is quite difficult and inefficient. An efficient and flexible method for solving the above question is to represent ${\mathit{b}}_i$ in terms of a latent variable $L_i\in \{1,2,\dots ,G\}$, which records each ${\mathit{b}}_i$’s cluster membership and conveys its parametric value to the distribution of ${\mathit{b}}_i$. Let $\mathit{L}=\{L_1,L_n,\dots ,L_n\}$, $\mathsf{\pi }=\{{\pi }_1,{\pi }_2,\dots ,{\pi }_G\}$, ${\mathsf{\mu }}^{*}=\{{\mathsf{\mu }}_1^{*},{\mathsf{\mu }}_2^{*},\dots ,{\mathsf{\mu }}_G^{*}\}$ and $\mathsf{\Omega }=\{{\mathsf{\Omega }}_1,{\mathsf{\Omega }}_2,\dots ,{\mathsf{\Omega }}_G\}$ in which ${\mathsf{\Omega }}_g=\mathrm{diag}({\omega }_{g1},{\omega }_{g2},\dots ,{\omega }_{gq})$, and they can be reformulated by

$$L_i|\mathsf{\pi }\stackrel{\mathrm{i}.\mathrm{i}.\mathrm{d}}{\sim }\sum _{g=1}^G{\pi }_g{\delta }_g(·)\mathrm{and}(\mathsf{\pi },{\mathsf{\mu }}^{*},\mathsf{\Omega })\sim f_1\left(\mathsf{\pi }\right)f_2\left({\mathsf{\mu }}^{*}\right)f_3\left(\mathsf{\Omega }\right),$$

where ${\delta }_g(·)$ denotes a discrete probability measure concentrated at g, $f_1\left(\mathsf{\pi }\right)$ is specified by the stick-breaking prior as given in Equation (5), $f_2\left({\mathsf{\mu }}^{*}\right)={\prod }_{g=1}^G{f}_2\left({\mathsf{\mu }}_g^{*}\right)$ and $f_3\left(\mathsf{\Omega }\right)={\prod }_{g=1}^G{\prod }_{j=1}^q{f}_3\left({\omega }_{g_j}\right)$. Here, ${\mathsf{\mu }}_g^{*}$ and ${\omega }_{g_j}$ can be specified by

$${\mathsf{\mu }}_g^{*}|\mathsf{\xi },\mathsf{\Psi }\stackrel{\mathrm{i}.\mathrm{i}.\mathrm{d}}{\sim }{N}_q(\mathsf{\xi },\mathsf{\Psi }),\mathsf{\xi }|{\mathsf{\xi }}^0,{\mathsf{\Psi }}^0\sim N_q({\mathsf{\xi }}^0,{\mathsf{\Psi }}^0),{\psi }_j^{-1}|c_1,c_2\sim \Gamma (c_1,c_2)\mathrm{for}j=1,2,\dots ,q,$$

(6)

$${\omega }_{g_j}^{-1}|{\omega }_j^a,{\varpi }_j\sim \Gamma ({\omega }_j^a,{\varpi }_j)\mathrm{and}{\varpi }_j|{\varpi }_j^a,{\varpi }_j^b\sim \Gamma ({\varpi }_j^a,{\varpi }_j^b),$$

respectively, where $\mathsf{\Psi }=\mathrm{diag}({\psi }_1,{\psi }_2,\dots ,{\psi }_q)$, $\Gamma (c_1,c_2)$ denotes the Gamma distribution with parameters $c_1$ and $c_2$, and ${\mathsf{\xi }}^0,{\mathsf{\Psi }}^0,c_1,c_2,{\omega }_j^a,{\varpi }_j^a$ and ${\varpi }_j^b$ are prespecified hyperparameters [18]. We consider ${\mathsf{\xi }}^0=0_{q\times 1}$, ${\mathsf{\Psi }}^0=I_q$, $c_1=11$, $c_2=2.5$, ${\omega }_j^a=3$, ${\varpi }_j^a=n$ and ${\varpi }_j^b=10$ in the paper. Given the values of $L_i$, ${\mathsf{\mu }}^{*}$ and $\mathsf{\Omega }$, we can sample ${\mathit{b}}_i$ from $N_q({\mathsf{\mu }}_{L_i},{\mathsf{\Omega }}_{L_i})$ with ${\mathsf{\mu }}_{L_i}={\mathsf{\mu }}_{L_i}^{*}-{\Sigma }_{g=1}^G{\mathsf{\mu }}_g^{*}$.

Hereafter we mainly introduce the variable selection principle of BLasso method [14,16] for the proposed joint model JMSRE. We need to identify not only the important variables in models (1) and (2). In particular, we also need to identify whether the parameter $\mathsf{\varphi }$ is $0_{q\times 1}$, and our proposed BLasso method does just that. In general, the prior distribution of the regression parameters is set to a multivariate normal distribution. Following the idea of Bayesian Lasso inference [14], we consider the following hierarchical priors for $\mathsf{\beta }$, $\mathsf{\gamma }$ and $\mathsf{\varphi }$:

$$\begin{array}{cc}\hfill \mathsf{\beta }|{\mathit{H}}_{\beta }\sim N_{p_1}(0,{\mathit{H}}_{\beta }),& \mathrm{with}{\mathit{H}}_{\beta }=diag(h_{{\beta }_1}^2,h_{{\beta }_2}^2,\dots ,h_{{\beta }_{p_1}}^2),\hfill \\ \hfill f(h_{{\beta }_1}^2,h_{{\beta }_2}^2,\dots ,& h_{{\beta }_{p_1}}^2)=\prod _{j=1}^{p_1}\frac{{\vartheta }_{{\beta }_j}^2}{2}exp\left(-\frac{{\vartheta }_{{\beta }_j}^2}{2}{h}_{{\beta }_j}^2\right),\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill \mathsf{\gamma }|{\mathit{H}}_{\gamma }\sim N_{p_2}(0,{\mathit{H}}_{\gamma }),& \mathrm{with}{\mathit{H}}_{\gamma }=diag(h_{{\gamma }_1}^2,h_{{\gamma }_2}^2,\dots ,h_{{\gamma }_{p_2}}^2),\hfill \\ \hfill f(h_{{\gamma }_1}^2,h_{{\gamma }_2}^2,\dots ,& h_{{\gamma }_{p_2}}^2)=\prod _{j=1}^{p_2}\frac{{\vartheta }_{{\gamma }_j}^2}{2}exp\left(-\frac{{\vartheta }_{{\gamma }_j}^2}{2}{h}_{{\gamma }_j}^2\right),\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill \mathsf{\varphi }|{\mathit{H}}_{\varphi }\sim N_q(0,{\mathit{H}}_{\varphi }),& \mathrm{with}{\mathit{H}}_{\varphi }=diag(h_{{\varphi }_1}^2,h_{{\varphi }_2}^2,\dots ,h_{{\varphi }_q}^2),\hfill \\ \hfill f(h_{{\varphi }_1}^2,h_{{\varphi }_2}^2,\dots ,& h_{{\varphi }_q}^2)=\prod _{j=1}^q\frac{{\vartheta }_{{\varphi }_j}^2}{2}exp\left(-\frac{{\vartheta }_{{\varphi }_j}^2}{2}{h}_{{\varphi }_j}^2\right),\hfill \end{array}$$

(9)

where ${\mathsf{\vartheta }}_{\beta }=\left\{{\vartheta }_{{\beta }_1},{\vartheta }_{{\beta }_2},\dots ,{\vartheta }_{{\beta }_{p_1}}\right\}$, ${\mathsf{\vartheta }}_{\gamma }=\left\{{\vartheta }_{{\gamma }_1},{\vartheta }_{{\gamma }_2},\dots ,{\vartheta }_{{\gamma }_{p_2}}\right\}$ and ${\mathsf{\vartheta }}_{\varphi }=\left\{{\vartheta }_{{\varphi }_1},{\vartheta }_{{\varphi }_2},\dots ,{\vartheta }_{{\varphi }_q}\right\}$ are the regularization parameters that control the tail decay. In particular, in order to better control the effect of tail decay, this paper sets different regularization parameters for different components of the same parameter. Inspired by Park and Casella [14], we further consider the following super-priorities for these tuning parameters:

$$\begin{array}{cc}\hfill {\vartheta }_{{\beta }_j}^2& \sim \Gamma \left(a_{{\vartheta }_{\beta }},b_{{\vartheta }_{\beta }}\right),j=1,2,\dots ,p_1,\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill {\vartheta }_{{\gamma }_j}^2& \sim \Gamma \left(a_{{\vartheta }_{\gamma }},b_{{\vartheta }_{\gamma }}\right),j=1,2,\dots ,p_2,\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill {\vartheta }_{{\varphi }_j}^2& \sim \Gamma \left(a_{{\vartheta }_{\varphi }},b_{{\vartheta }_{\varphi }}\right),j=1,2,\dots ,q.\hfill \end{array}$$

(12)

3.2. Bayesian analysis of joint model

To obtain Bayesian estimates of unknown parameters $\mathsf{\beta },\mathsf{\phi },{\sigma }^2,\mathit{b},\mathsf{\gamma }$ and $\mathsf{\varphi }$, a hybrid algorithm combining the block Gibbs sampler and the Metropolis-Hastings algorithm is employed to iteratively draw $\{\mathsf{\beta },\mathsf{\phi },{\sigma }^2,\mathit{b},\mathsf{\gamma },\mathsf{\varphi }\}$ as follows.

(A) Conditional distribution of $\mathsf{\beta }$

It follows from Equation (3) and Equation (7) that the conditional posterior distribution $p\left(\mathsf{\beta }\right|\mathsf{\phi },{\sigma }^2,\mathit{b},{\mathit{y}}^{*})$ is expressed as

$$\begin{array}{c}\hfill p\left(\mathsf{\beta }\right|\mathsf{\phi },{\sigma }^2,\mathit{b},{\mathit{y}}^{*})\propto exp\left\{-\frac{1}{2}\left[\sum _{i=1}^n\sum _{j=1}^{n_i}\frac{1}{{\sigma }^2}\left(y_{ij}^{*}-{\mathit{X}}_{ij}^{\top }\mathsf{\beta }-{\mathit{B}}^{\top }\left(t_{ij}\right)\mathsf{\phi }-{\mathit{Z}}_{ij}^{\top }{\mathit{b}}_i\right)+{\mathsf{\beta }}^{\top }{\mathit{H}}_{\beta }^{-1}\mathsf{\beta }\right]\right\},\end{array}$$

which yields

$$\begin{array}{c}\hfill \mathsf{\beta }|\mathsf{\phi },{\sigma }^2,\mathit{b},{\mathit{y}}^{*}\sim N_{p_1}\left({\mathcal{A}}_{\beta },{\mathcal{V}}_{\beta }\right),\end{array}$$

where ${\mathcal{V}}_{\beta }=\sum _{i=1}^n\sum _{j=1}^{n_i}{\sigma }^{-2}{\mathit{X}}_{ij}^{\top }{\mathit{X}}_{ij}+{\mathsf{\Sigma }}_{\beta }^{-1}$, ${\mathcal{A}}_{\beta }={\mathcal{V}}_{\beta }\left(\sum _{i=1}^n\sum _{j=1}^{n_i}{\mathit{X}}_{ij}\left(y_{ij}^{*}-{\mathit{Z}}_{ij}^{\top }{\mathit{b}}_i-{\mathit{B}}^{\top }\left(t_{ij}\right)\mathsf{\phi }\right){\sigma }^{-2}\right)$.

(B) Conditional distribution of $\mathsf{\phi }$ According to Equation (3) and the prior of $\mathsf{\phi }$ in Equation (4), the conditional distribution $p\left(\mathsf{\phi }\right|\mathsf{\beta },{\sigma }^2,\mathit{b},{\mathit{y}}^{*})$ is given by

$$\begin{array}{c}\hfill p\left(\mathsf{\phi }\right|\mathsf{\beta },{\sigma }^2,\mathit{b},{\mathit{y}}^{*})\propto exp\left\{-\frac{1}{2}\left[\sum _{i=1}^n\sum _{j=1}^{n_i}\frac{1}{{\sigma }^2}\left(y_{ij}^{*}-{\mathit{X}}_{ij}^{\top }\mathsf{\beta }-{\mathit{B}}^{\top }\left(t_{ij}\right)\mathsf{\phi }-{\mathit{Z}}_{ij}^{\top }{\mathit{b}}_i\right)+{\mathsf{\phi }}^{\top }{\left({\mathit{H}}_{\phi }^0\right)}^{-1}\mathsf{\phi }\right]\right\},\end{array}$$

which yields

$$\begin{array}{c}\hfill \mathsf{\phi }|\mathsf{\beta },{\sigma }^2,\mathit{b},{\mathit{y}}^{*}\sim N_L\left({\mathcal{A}}_{\phi },\mathcal{V}\phi \right),\end{array}$$

where $\mathcal{V}\phi =\sum _{i=1}^n\sum _{j=1}^{n_i}{\sigma }^{-2}B{\left(t_{ij}\right)}^{\top }B\left(t_{ij}\right)+{\Sigma }_{\phi }^{-1}$, $\mathcal{A}\phi =\mathcal{V}\phi \left(\sum _{i=1}^n\sum _{j=1}^{n_i}B\left(t_{ij}\right)\left(y_{ij}^{*}-X_{ij}^{\top }\mathsf{\beta }-{\mathit{Z}}_{ij}^{\top }{\mathit{b}}_i\right){\sigma }^{-2}\right)$.

(C) Conditional distribution of ${\sigma }^{-2}$

It follows from Equation (3) and the prior of ${\sigma }^2$ in Equation (4) that the conditional distribution $p\left({\sigma }^{-2}\right|\mathsf{\beta },\mathsf{\phi },b,y^{*})$ is given by

$$\begin{array}{cc}\hfill p\left({\sigma }^{-2}\right|\mathsf{\beta },\mathsf{\phi },b,y^{*})& \propto exp\left\{-\left(\frac{1}{2}\sum _{i=1}^n\sum _{j=1}^{n_i}\left(y_{ij}^{*}-{\mathit{X}}_{ij}^{\top }\mathsf{\beta }-{\mathit{B}}^{\top }\left(t_{ij}\right)\phi -{\mathit{Z}}_{ij}^{\top }{\mathit{b}}_i\right)+b_0\right){\sigma }^{-2}\right\}\hfill \\ \hfill & \times {\left({\sigma }^{-2}\right)}^{\frac{1}{2}\sum _{i=1}^n{n}_i+a_0-1},\hfill \end{array}$$

which yields

$$\begin{array}{c}\hfill {\sigma }^{-2}|\mathsf{\beta },\mathsf{\phi },b,y^{*}\sim \Gamma \left(\frac{1}{2}\sum _{i=1}^n{n}_i+a_0,\frac{1}{2}\left(y_{ij}^{*}-{\mathit{X}}_{ij}^{\top }\mathsf{\beta }-{\mathit{B}}^{\top }\left(t_{ij}\right)\phi -{\mathit{Z}}_{ij}^{\top }{\mathit{b}}_i\right)+b_0\right).\end{array}$$

(D) Conditional distribution of $b_i$

For reasons of space, the sampling of $b_i,i=1,2,\dots ,n$ follows the steps in Appendix A, which can also be seen in Tang et al. [24]

(E) Conditional distribution of $\mathsf{\gamma }$

It follows from Equation (3) and Equation (8) that the conditional distribution $p\left(\mathsf{\gamma }\right|b,\mathsf{\varphi },T,\Delta )$ is proportional to

$$\begin{array}{c}\hfill \prod _{i=1}^n{\left[\frac{exp({\mathit{W}}_i^{\top }\mathsf{\gamma }+{\mathsf{\varphi }}^{\top }{\mathit{b}}_i)}{{\sum }_{j\in {\mathcal{R}}_i}exp({\mathit{W}}_j^{\top }\mathsf{\gamma }+{\mathsf{\varphi }}^{\top }{\mathit{b}}_j)}\right]}^{{\delta }_i}exp\left\{-\frac{1}{2}{\mathsf{\gamma }}^T{H}_{\mathsf{\gamma }}^{-1}\mathsf{\gamma }\right\},\end{array}$$

which is not a familiar distribution. Therefore, the well-known Metropolis-Hastings (MH) algorithm is adopted to simulate observations from the above given conditional distribution, which is implemented as follows. Given the current value ${\mathsf{\gamma }}^{\left(m\right)}$, new candidates $\mathsf{\gamma }$ are generated from $N_{p_2}({\mathsf{\gamma }}^{\left(m\right)},{\sigma }_{\mathsf{\gamma }}^2{\Sigma }_{\mathsf{\gamma }})$, where ${\Sigma }_{\mathsf{\gamma }}={\left(-{\partial }^2{(ln\left(p\left(\mathsf{\gamma }\right|b,\mathsf{\varphi },T,\Delta )\right)/\partial \mathsf{\gamma }\partial {\mathsf{\gamma }}^{\top }|}_{\mathsf{\gamma }={\mathsf{\gamma }}^{\left(m\right)}}\right)}^{-1}$. The new ${\mathsf{\gamma }}^{\left(m\right)}$ is accepted with probability

$$\begin{array}{c}\hfill min\left\{1,\frac{p\left(\mathsf{\gamma }\right|b,\mathsf{\varphi },T,\Delta )}{p\left({\mathsf{\gamma }}^{\left(m\right)}\right|b,\mathsf{\varphi },T,\Delta )}\right\},\end{array}$$

where

$$\begin{array}{c}\hfill {\Sigma }_{\mathsf{\gamma }}={\left(\sum _{i=1}^n\frac{{\sum }_{j\in {\mathcal{R}}_i}exp(·){\mathit{W}}_j{\mathit{W}}_j^{\top }{\sum }_{j\in {\mathcal{R}}_i}exp(·)+{\sum }_{j\in {\mathcal{R}}_i}exp(·){\mathit{W}}_j{\sum }_{j\in {\mathcal{R}}_i}exp(·){\mathit{W}}_j^{\top }}{{\left({\sum }_{j\in {\mathcal{R}}_i}exp(·)\right)}^2}+H_{\mathsf{\gamma }}^{-1}\right)}^{-1}.\end{array}$$

with $exp(·)=exp({\mathit{W}}_j^{\top }{\mathsf{\gamma }}^{\left(m\right)}+{\mathsf{\varphi }}^{\top }{\mathit{b}}_j).$ The variance coefficient ${\sigma }_{\mathsf{\gamma }}^2$ can be chosen such that the average acceptance rate are approximately 0.25 or more.

(F) Conditional distribution of $\mathsf{\varphi }$

From Equation (3) and Equation (9), the conditional distribution $p\left(\mathsf{\varphi }\right|\mathsf{\gamma },b,T,\Delta )$ is proportional to

$$\begin{array}{c}\hfill \prod _{i=1}^n{\left[\frac{exp({\mathit{W}}_i^{\top }\mathsf{\gamma }+{\mathsf{\varphi }}^{\top }{\mathit{b}}_i)}{{\sum }_{j\in {\mathcal{R}}_i}exp({\mathit{W}}_j^{\top }\mathsf{\gamma }+{\mathsf{\varphi }}^{\top }{\mathit{b}}_j)}\right]}^{{\delta }_i}exp\left\{-\frac{1}{2}{\mathsf{\gamma }}^T{H}_{\mathsf{\varphi }}^{-1}\mathsf{\gamma }\right\},\end{array}$$

which is not a familiar distribution. Similar to above (E), given the current value ${\mathsf{\varphi }}^{\left(m\right)}$, new candidates $\mathsf{\varphi }$ are generated from $N_q({\mathsf{\varphi }}^{\left(m\right)},{\sigma }_{\mathsf{\varphi }}^2{\Sigma }_{\mathsf{\varphi }})$, where ${\Sigma }_{\mathsf{\varphi }}={\left(-{\partial }^2{(ln\left(p\left(\mathsf{\varphi }\right|b,\mathsf{\gamma },T,\Delta )\right)/\partial \mathsf{\varphi }\partial {\mathsf{\varphi }}^{\top }|}_{\mathsf{\varphi }={\mathsf{\varphi }}^{\left(m\right)}}\right)}^{-1}$. The new ${\mathsf{\varphi }}^{\left(m\right)}$ is accepted with probability

$$\begin{array}{c}\hfill min\left\{1,\frac{p\left(\mathsf{\varphi }\right|\mathsf{\gamma },b,T,\Delta )}{p\left({\mathsf{\varphi }}^{\left(m\right)}\right|\mathsf{\gamma },b,T,\Delta )}\right\},\end{array}$$

where

$$\begin{array}{c}\hfill {\Sigma }_{\mathsf{\varphi }}={\left(\sum _{i=1}^n\frac{{\sum }_{j\in {\mathcal{R}}_i}exp(·){\mathit{b}}_j{\mathit{b}}_j^{\top }{\sum }_{j\in {\mathcal{R}}_i}exp(·)+{\sum }_{j\in {\mathcal{R}}_i}exp(·){\mathit{b}}_j{\sum }_{j\in {\mathcal{R}}_i}exp(·){\mathit{b}}_j^{\top }}{{\left({\sum }_{j\in {\mathcal{R}}_i}exp(·)\right)}^2}+H_{\mathsf{\varphi }}^{-1}\right)}^{-1},\end{array}$$

and $exp(·)=exp({\mathit{W}}_j^{\top }\mathsf{\gamma }+{{\mathsf{\varphi }}^{\left(m\right)}}^{\top }{\mathit{b}}_j).$ The variance coefficient ${\sigma }_{\mathsf{\varphi }}^2$ can be chosen such that the average acceptance rate are approximately 0.25 or more.

By the above iterative process, we can obtain a series of sample $\{({\mathsf{\beta }}^{\left(m\right)},{\mathsf{\phi }}^{\left(m\right)},{{\sigma }^2}^{\left(m\right)},b_i^{\left(m\right)},{\mathsf{\gamma }}^{\left(m\right)},{\mathsf{\varphi }}^{\left(m\right)}):m=1,2,\dots ,M\}$. Then, Bayesian estimates of $\mathsf{\beta }$, $\mathsf{\phi }$, ${\sigma }^2$, ${\mathit{b}}_i$, $\mathsf{\gamma }$ and $\mathsf{\varphi }$ can be obtained by

$$\begin{array}{c}\hfill \widehat{\mathsf{\beta }}=\frac{1}{M}\sum _{m=1}^M{\mathsf{\beta }}^{\left(m\right)},\widehat{\mathsf{\phi }}=\frac{1}{M}\sum _{m=1}^M{\mathsf{\phi }}^{\left(m\right)},\widehat{{\sigma }^2}=\frac{1}{M}\sum _{m=1}^M{{\sigma }^2}^{\left(m\right)},\end{array}$$

$$\begin{array}{c}\hfill {\widehat{b}}_i=\frac{1}{M}\sum _{m=1}^M{b}_i^{\left(m\right)}\widehat{\mathsf{\gamma }}=\frac{1}{M}\sum _{m=1}^M{\mathsf{\gamma }}^{\left(m\right)},\widehat{\mathsf{\varphi }}=\frac{1}{M}\sum _{m=1}^M{\mathsf{\varphi }}^{\left(m\right)}.\end{array}$$

Similarly, the consistent estimates of the posterior covariance matrices of $\mathit{var}\left(\mathsf{\beta }\right|{\mathit{y}}^{*},\mathit{X},\mathit{Z})$, $\mathit{var}\left(\mathsf{\phi }\right|{\mathit{y}}^{*},\mathit{X},\mathit{Z})$, $\mathit{var}\left({\mathsf{\sigma }}^2\right|{\mathit{y}}^{*},\mathit{X},\mathit{Z})$, $\mathit{var}\left(\mathsf{\gamma }\right|\mathit{W},\mathit{T},\Delta )$, and $\mathit{var}\left(\mathsf{\varphi }\right|{\mathit{y}}^{*},\mathit{X},\mathit{Z},\mathit{W},\mathit{T},\Delta )$ can be obtained via the sample covariance matrices. For example

$$\begin{array}{c}\hfill \widehat{\mathit{var}}\left(\mathsf{\beta }\right|{\mathit{y}}^{*},\mathit{X},\mathit{Z})=\frac{1}{M-1}\sum _{m=1}^M({\mathsf{\beta }}^{\left(m\right)}-\widehat{\mathsf{\beta }}){({\mathsf{\beta }}^{\left(m\right)}-\widehat{\mathsf{\beta }})}^{\top }.\end{array}$$

Therefore, the variance of the corresponding parameter can be obtained through the diagonal elements of the sample covariance matrix of the random sample sequence.

4. Simulation Studies

In this section, we conducted three simulation studies to investigate the finite performance of the above proposed methods.

We considered the model defined in models (1) and (2) with 200 individuals. Here, the specific model is as follows:

$$y_{ij}^{*}|{\mathit{b}}_i=X_{1ij}{\mathsf{\beta }}_1+X_{2ij}{\mathsf{\beta }}_2+X_{3ij}{\mathsf{\beta }}_3+X_{4ij}{\mathsf{\beta }}_4+X_{5ij}{\mathsf{\beta }}_5+X_{6ij}{\mathsf{\beta }}_6+sin\left(2\pi t_{ij}\right)+{\mathit{b}}_i+{\epsilon }_{ij},$$

$${\lambda }_i\left(t|b_i\right)={\lambda }_0\left(t\right)exp\left(W_{1i}{\mathsf{\gamma }}_1+W_{2i}{\mathsf{\gamma }}_2+W_{3i}{\mathsf{\gamma }}_3+W_{4i}{\mathsf{\gamma }}_4+\mathsf{\varphi }{b}_i\right).$$

The data were generated as follows: observation time $t_{ij}$ was generated randomly between 0 and 3, the covariates $X_{1ij}$ and $X_{6ij}$ followed the Bernoulli distribution with the success probability 0.5 and 0.3, respectively. The covariates $X_{2ij},X_{3ij},X_{4ij}$ and $X_{5ij}$ were generated form the multivariate normal distribution $N_4\left(0,\Sigma \right)$ with mean vector 0 and covariance matrix $\Sigma$. And the covariance matrix $\Sigma$ is a symmetric positive definite matrix whose diagonal elements are 1 and the rest are $0.5$. The random error ${\epsilon }_{ij}$ is generated from the normal distribution with mean 0 and variance ${\sigma }^2=0.6^2$. And $W_i={(W_{i1},W_{i2},W_{i3},W_{i4})}^{\top }={(X_{3i1},X_{4i1},X_{5i1},X_{6i1})}^{\top }$. The baseline hazard function ${\lambda }_0\left(t\right)=0.7$ and $\mathsf{\varphi }=0.6$. The censoring time $C_i$ was generated from the Uniform distribution $U\left[0,3\right]$, and $T_i^{*}$ was generated from the Exponential distribution with mean $1/{\lambda }_i\left(t|b_i\right)$, $T_i=min(T_i^{*},C_i)$. Our main purpose is to use the proposed approaches to identify the unimportant covariates and estimate nonzero coefficients. Bayesian results were obtained from 200 replications.

To show that our proposed method above can both accurately estimate parameters of interest and identify unimportant variables, but also flexibly capture the feature of nonlinear function $g\left(t\right)$ and random effects $b_i$, we consider three different simulation studies as follows.

Simulation I.

$$\begin{array}{c}\hfill \mathsf{\beta }={\left({\mathsf{\beta }}_1,\dots ,{\mathsf{\beta }}_6\right)}^{\top }={\left(1,0,0,-0.5,0.5,-1\right)}^{\top },\mathsf{\gamma }={\left({\mathsf{\gamma }}_1,\dots ,{\mathsf{\gamma }}_4\right)}^{\top }={\left(0,1,-0.5,0\right)}^{\top },\end{array}$$

$$\begin{array}{c}\hfill g\left(t\right)=sin\left(\frac{3}{4}\pi t\right),b_i\stackrel{\mathrm{i}.\mathrm{i}.\mathrm{d}}{\sim }0.6N(-0.8,0.1^2)+0.4N(1.2,0.5^2),\end{array}$$

Simulation II.

$$\begin{array}{c}\hfill \mathsf{\beta }={\left({\mathsf{\beta }}_1,\dots ,{\mathsf{\beta }}_6\right)}^{\top }={\left(1,0,0,-0.5,0.5,-1\right)}^{\top },\mathsf{\gamma }={\left({\mathsf{\gamma }}_1,\dots ,{\mathsf{\gamma }}_4\right)}^{\top }={\left(0,1,-0.5,0\right)}^{\top },\end{array}$$

$$\begin{array}{c}\hfill g\left(t\right)=t,b_i\stackrel{\mathrm{i}.\mathrm{i}.\mathrm{d}}{\sim }0.4N(0,0.3^2)+0.3N(-1.5,0.1^2)+0.3N(1.5,0.1^2),\end{array}$$

Simulation III.

$$\begin{array}{c}\hfill \mathsf{\beta }={\left({\mathsf{\beta }}_1,\dots ,{\mathsf{\beta }}_6\right)}^{\top }={\left(1,0.5,-0.5,-0.5,0.5,-1\right)}^{\top },\mathsf{\gamma }={\left({\mathsf{\gamma }}_1,\dots ,{\mathsf{\gamma }}_4\right)}^{\top }={\left(-0.5,1,-0.5,1\right)}^{\top },\end{array}$$

$$\begin{array}{c}\hfill g\left(t\right)=t^2,b_i\stackrel{\mathrm{i}.\mathrm{i}.\mathrm{d}}{\sim }N(0,0.8^2).\end{array}$$

which the mean censoring rates of the survival times in the three simulation studies are 44%, 45%, and 37%, respectively.

For each of the above generated data sets, the proposed semiparametric Bayesian procedure was used to simultaneously evaluate Bayesian estimates of unknown parameters and select the important covariates. The prior hyperparameters $a_{{\vartheta }_{\mathsf{\beta }}}=a_{{\vartheta }_{\mathsf{\gamma }}}=a_{{\vartheta }_{\mathsf{\varphi }}}=1,b_{{\vartheta }_{\mathsf{\gamma }}}=b_{{\vartheta }_{\mathsf{\beta }}}=b_{{\vartheta }_{\mathsf{\varphi }}}=0.1$, corresponding to the hyperparameters of the hyper priors for the adjustment coefficients in Equations (10)-(). $a_0=1,b_0=1$ and $H_{\phi }^0=100I_4$, corresponding to the prior parameters of ${\sigma }^2,\phi$. We set the degree of B-splines $d=3$, the number of konts $K=4$, and $G=25$. To investigate the convergence of the proposed algorithm, we calculated the estimated potential scale reduction (EPSR) values of parameters. In addition, we also need to test whether the nonlinear function fitted by the B-splines method is convergent. It can be seen from Figure 1 that the EPSR values were less than 1.2 after about 3000 iterations in three simulation studies. Thus, 3000 ($M=3000$) observations were collected to calculated Bayesian estimates of parameters after 3000 iterations in producing Bayesian results for each of 200 replications. Results obtained under three simulation studies were reported in Table 1, where “Bias” is the difference between the true value and the mean of the estimates based on 200 replications, “RMS” is the root mean square between the estimates based on 200 replications and its true value, and “F0” is the proportion that parameter is identified to be zero in 200 replications in terms of the criterion that a parameter is identified to be 0 if its $95\%$ confidence interval contains 0.

Results of Table 1 indicated that (i) Bayesian estimates of parameters are reasonably accurate because their absolute biases are less than 0.10 and the root mean square of their RMS values are less than 0.20; (ii) BLasso could identify the correct models in most cases regardless of prior inputs of parameters because the F0 values corresponding to the important covariates were less than $10\%$, but the F0 values corresponding to unimportant covariates were more than $90\%$. The performance of the proposed approach to recover the nonlinear function $g\left(t\right)$ can be measured by the root mean square error

$$\begin{array}{c}\hfill \mathrm{RASE}\left(g^{\left(r\right)}\right)=\sqrt{\frac{1}{300}{\sum }_{l=1}^{300}\left(g\left(u_l\right)-{\widehat{g}}^{\left(r\right)}\left(u_l\right)\right)},r=1,2,\dots ,200,\end{array}$$

where ${\widehat{g}}^{\left(r\right)}\left(t\right)=B^{\top }\left(t\right){\widehat{\mathsf{\phi }}}^{\left(r\right)}$, ${\widehat{\mathsf{\phi }}}^{\left(r\right)}$ represents the Bayesian estimated value of the parameter vector $\phi$ in the r-th replication. Similar to the root mean square error of nonlinear function, we also calculate the root mean square error of random effects. Figure 2 plots the estimated curve and estimated density of the nonlinear function $\widehat{g}\left(t\right)$ and the random effects $b_i$ of the replication in which the mean of the RASE of the nonlinear function and the random effects is in the middle of the 200 replications, and against true curves and true density in three simulation studies, respectively. Inspection of Figure 2 showed that (i) the Bayesian B-splines method proposed in this paper can be flexible enough fit the true curve of the nonlinear function $g\left(t\right)$; (ii) the CDPMM method proposed is flexible enough to capture the general shapes of our considered three distribution assumptions for $b_i$; (iii) the estimated 95% confidence region of the nonlinear function can cover its true curve with a fairly narrow area. From Table 2, based on the results of 200 replications in three simulation studies, the estimated means and standard deviation (SD) of random effects $b_i$ are enough near to their corresponding true values, and 25%, 50% and 75% quantile of the mean of the RASE of the nonlinear function and the random effects are enough small and close to illustrate that it is robust to apply CDPMM to estimate random effects. All these findings indicated that (i) our proposed Bayesian B-splines can be flexible enough fit the true curve of the nonlinear; (ii) Bayesian procedure could well capture the true information of $b_i$ regardless of their true distributions and forms; (iii) BLasso could identify the true model with a high probability.

5. An Example

In this section, we apply the method proposed in the previous sections to the MA.5 research experiment conducted by the Clinical Trial Group of the National Cancer Institute of Canada. The data contains 716 women with early-stage breast cancer before menopause. 356 patients were randomly selected to receive cyclophosphamide, epirubicin, and fluorouracil (CEF) adjuvant chemotherapy as the experimental group. The remaining 360 patients received cyclophosphamide, methotrexate, and fluorouracil (CMF) adjuvant chemotherapy as the control group of the trial. In clinical trials, visits were made before the start of treatment, each of the 6 treatment cycles, and every three months after treatment. At each visit, medical history and physical examination are accepted, and BCQ (Breast Cancer Questionnaire) is used to assess the patient’s quality of life. The dataset has a total of 7807 observations. By the end of the study, a total of 366 patients had died, and the censorship rate was about $49\%$. For a detailed study of this data, please refer to the literature of Song et al. [22] and Levine et al. [24]. We linearly convert the evaluated BCQ score into a unit interval $(0,1)$, and the longitudinal data constrained on the interval $(0,1)$ is the longitudinal proportion data that we are interested in. The survival time of interest in the trial is the recurrence-free survival time (RFS), which is defined as the time from randomization to disease recurrence. Different treatment options, age, and the number of tumor-positive lymph nodes may directly affect RFS and the quality of life of patients. Similar to Tang et al. [17], we fitted the MA.5 research experiment data set to the following model:

$$\begin{array}{c}\hfill y_{ij}^{*}|b_i={\mathsf{\beta }}_1{\mathrm{EM}}_i+{\mathsf{\beta }}_2{{\mathrm{NODE}}_{\_}\mathrm{POS}}_i+{\mathsf{\beta }}_3{\mathrm{AGE}}_i+g\left(t_{ij}\right)+b_i+{\epsilon }_{ij},\end{array}$$

$$\begin{array}{c}\hfill {\lambda }_i\left(t|b_i\right)={\lambda }_0\left(t\right)exp\left({\mathsf{\gamma }}_1{\mathrm{EM}}_i+{\mathsf{\gamma }}_2{{\mathrm{NODE}}_{\_}\mathrm{POS}}_i+{\mathsf{\gamma }}_3{\mathrm{AGE}}_i+\mathsf{\varphi }{b}_i\right),\end{array}$$

where $y_{ij}^{*}$ represents the BCQ score after logit function transformation, ${\mathrm{EM}}_i$ is a two-class treatment index, ${\mathrm{EM}}_i=1$ indicates that the i-th patient underwent CEF treatment, ${\mathrm{EM}}_i=0$ indicates that the i-th patient underwent CMF treatment. Age and the number of lymph node metastases are binary variables. Patients who are less than or equal to 40 years old belong to the younger group, that is, $\mathrm{AGE}=1$. Patients who are older than 40 years old belong to the elderly group, that is, $\mathrm{AGE}=0$. When the number of lymph node metastases is 0-3, $\mathrm{NODE}\_\mathrm{POS}=0$, otherwise it is 1. The nonlinear $g\left(t\right)$ is estimated by using cubic B-splines function, and the domain of cubic B-splines function is $[min\left(t_{ij}\right),max\left(t_{ij}\right)]$. The prior distributions and the values of all hyperparameters in the case study are the same as those set in the simulation study above. Based on the above settings, we calculated EPSR values for all parameters and are shown in the left panel in Figure 3, which shows that after about 3000 iterations, all EPSR values are less than 1.2. Therefore, we use the 3000 iterations after the 3000th time to calculate the Bayesian estimation. The example analysis results are shown in Table 3.

Inspection of Figure 3 showed that (i) the estimated density of random effects is approximately normally distributed. The distribution settings for random effects in our Simulation III correspond to. This shows that our CDPMM method is effective and reasonable; (ii) the estimated nonlinear function first decreases, then increases and finally stabilizes with time t. Although there are some fluctuations in the tail, this is reasonable and can be ignored. From Table 3, we can see that (i) the risk ratio of randomly receiving CEF and CMF treatment is $\mathrm{HR}=exp\left({\mathsf{\gamma }}_1\right)=71.605\%$, which means patients who randomly receive CEF chemotherapy have a lower risk; (ii) the credible interval $(0.171,0.321)$ of ${\mathsf{\beta }}_1$ does not contain 0, indicating that different adjuvant chemotherapy regimens have a significant impact on the QOL of patients and CEF chemotherapy is more toxic than CMF chemotherapy; (iii) the risk ratio of the number of lymph node metastases greater than or equal to 4 to less than 4 is $\mathrm{HR}=exp\left({\mathsf{\gamma }}_2\right)=205.238\%$, which means that the number of lymph node metastases is greater than or equal to 4. The risk of breast cancer recurrence is greater in patients, and the RFS is shorter; (iv)the regression coefficient ${\mathsf{\beta }}_2$ of the patient’s lymph node metastasis number greater than or equal to 4 is $0.309$, and its credible interval does not contain 0, which is highly significant. This suggests that patients with more lymph nodes have a lower quality of life, which is also consistent with clinical experience; (v) the risk ratio between the young group and the old group is $\mathrm{HR}=exp\left({\mathsf{\gamma }}_3\right)=182.212\%$, which means the risk of breast cancer recurrence is higher and the RFS is shorter in the young group; (vii) the credible interval $(0.168,0.344)$ of ${\mathsf{\beta }}_3=0.263$ does not contain 0, indicating that different ages have a significant impact on the QOL of patients. It indicates that the quality of life of the olderly group is better than that of the young group. (viii) $\mathsf{\varphi }=0.259$ and the credible interval $(0.003,0.542)$ of $\mathsf{\varphi }$ does not contain 0, which means that $\mathsf{\varphi }$ is significantly different from 0. This shows that there is a significant correlation between the longitudinal proportional data and survival data, and the JMSRE model proposed in this paper is applicable and reasonable for the analysis of the MA.5 research experiment data.

6. Concluding Remarks

In this paper, a semiparametric joint model for longitudinal proportional data and survival data is proposed by relaxing the normality assumption of random effects and without specifying a nonlinear function that affects longitudinal responses. The advantages of the proposed model include: (i) the proposed method improves the flexibility of jointly modeling longitudinal proportional data and survival data; (ii) the proposed B-splines method can flexibly capture different nonlinear functions; (iii) the proposed CDPMM method can flexibly capture the unimodal, bimodal and multimodal features of random effects; (iv) the computational burden is not heavy, e.g. it takes about 4 minutes to run the replication in the above simulation study and about 2 hours to run the breast cancer dataset.

Our simulation studies and example analysis show that the Bayesian estimation proposed based on the joint model are quite accurate and relatively robust. The nonlinear function curve estimated by Bayesian B-splines is more flexible and can capture the true characteristics of the nonlinear function. And CDPMM method could also well capture the true information of $b_i$. BLasso method could identify the true model with a high probability. Compared with the method proposed by Song et al. [22] to jointly model longitudinal proportional data and survival data, the joint model proposed in this paper is more flexible.

The joint model of longitudinal proportional data and survival data proposed in this paper still has many unsolved problems and we need to solve the following problems in the future: (i) does not make any constraints on the form of the basic hazard function; (ii) the distribution of random effects does not depend on the normality assumption and still has robust inferences; (iii) consider more complex spline models, for example nodes can be automatically selected to improve the performance of the proposed model; (iv) consider a joint model for the variable longitudinal proportional outcome and the multivariate survival outcome.