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.

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 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 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 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 where

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:

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.

Simulation II.

Simulation III. 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 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.

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: 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.

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.