A joint latent class model of longitudinal and survival data with a time-varying membership probability

Joint latent class modelling has been developed considerably in the past two decades. In some instances, the models are linked by the latent class k (i.e. the number of subgroups), in others they are joined by shared random effects or a heterogeneous random covariance matrix. We propose an extension to the joint latent class model (JLCM) in which probabilities of subjects being in latent class k can be set to vary with time. This can be a more flexible way to analyse the effect of treatments to patients. For example, a patient may be in period I at the first visit time and may move to period II at the second visit time, implying the treatment the patient had before might be noneffective at the following visit time. For a dataset with these particular features, the joint latent class model which allows jumps among different subgroups can potentially provide more information as well as more accurate estimation and prediction results compared to the basic JLCM. A Bayesian approach is used to do the estimation and a DIC criterion is used to decide the optimal number of classes. Simulation results indicate that the proposed model produces accurate results and the time-varying JLCM outperforms the basic JLCM. We also illustrate the performance of our proposed JLCM on the aids data (Goldman et al., 1996).


I. INTRODUCTION
Joint models for longitudinal and survival data are a popular approach for handling the dependencies between these two types of data. They are particularly useful for modelling longitudinal data with informative dropout or performing time to event analysis with time dependent covariates, but they can also be used if the interest is on both processes, while also considering their association. Incorporating all information can help increase the efficiency and decrease the bias of inferences (Ibrahim et al., 2010).
Joint modelling of longitudinal measurements and survival data was first introduced in the early 90s. Schluchter (1992) pointed out the problem of non-ignorable censoring data in longitudinal studies and gave an approach which was based on lognormal survival. Self and Pawitan (1992), De Gruttola and Tu (1994) and Tsiatis et al. (1995) applied joint modelling approaches to AIDS research, which adjusted the inferences on the longitudinal responses with informative missing data by proposing the two-stage method to do the parameter estimation. Henderson et al. (2000) proposed a flexible joint model, where the longitudinal and survival processes were linked via a latent Gaussian process, allowing for both serial correlation and measurement error in the longitudinal model. Rizopoulos (2009) proposed a new computational approach, in the form of a fully exponential Laplace approximation for the joint modelling of survival and longitudinal data, making is feasible to handle high dimensional random effects structures in joint models. Motivated by a primary biliary cirrhosis (PBC) study (Murtaugh et al., 1994), a joint modelling of multi-state event times and longitudinal data with informative observation time points was proposed by Dai and Pan (2018), where the joint modelling was linked by random effects under no distributional assumption. This model extended the corrected score method (Wang, 2006) to cases with longitudinal data which is collected at informative time points.
Typically, joint models are linked by shared parameters/random effects or by random effects sharing the same distribution. The development of specialised software, such as JM (Rizopoulos, 2010), joineR (Philipson et al., 2012), stjm (Crowther et al., 2013) and JMBayes (Rizopoulos, 2016) among others, has also helped popularise this type of approach. Other work in the literature includes copula joint models (Kürüm et  Traditional joint modelling via shared random effects in a homogenous group, though widely used, is limited to describe only one pattern for subjects. In many circumstances, the presence of heterogeneous subgroups can not be ignored and thus potential subtypes of longitudinal and survival outcomes should be considered. In order to deal with the underlying subpopulations with different patterns of responses, joint latent class models can be used. These models are usually applied to analyse data from health and aging studies, e.g. prostate cancer, AIDS or dementia studies. Lin et al. (2002) first proposed the joint latent class model (JLCM) by extracting meaningful subgroups with respect to the joint distribution of the longitudinal data and the event process, which could be used to model informative dropout and help identify patterns of longitudinal measurements over time. They proposed and applied a joint latent class model of longitudinal PSA and the PCA onset, permitting distinct behaviours within each subgroup to model the association between longitudinal and time-to-event data. Garre et al. (2008) proposed a joint latent class changepoint model to improve the prediction of time to graft failure. They divided patients into two different groups, such that subjects in latent class 1 were modelled by an intercept-only randomeffects model and subjects within latent class 2 were modelled by a segmented random changepoint model, achieving better predictions of time to graft failure compared to other joint models.
Chiang (2009) postulated a more flexible joint latent varying-coefficient model to handle population heterogeneity, which generalised the correlation mechanism in the joint model with two latent processes, proposed by Henderson et al. (2000). Additionally, Liu et al. (2015) introduced a joint latent class model by integrating the joint random effects model and applying this model to an AIDS study to analyse the relationship between longitudinal CD4+ measurements and the time to death. This type of proposed model is appropriate when nonignorable heterogeneity exists among subjects. An overview of earlier work on joint latent class models is provided in a review article by Proust-Lima et al. (2014). Proust-Lima et al. (2017) also developed the R package lcmm which can implement various extended mixed models, including joint latent class models, using the maximum likelihood method. More recently, Andrinopoulou et al. (2020) introduced a Bayesian approach to joint latent class models, which incorporated a new latent class selection procedure via a mixture model and applied the proposed approach to cystic fybrosis data.
Reboussin and Anthony (2001) proposed the dynamic latent class regression model for longitudinal data by adding time-varying covariates into the baseline-category logistic regression model, Lin et al. (2014) added random effects to the aforementioned dynamic latent class model for longitudinal data, to jointly model the informative event and a class-specific logistic model with shared random effects. Zhang and Simonoff (2020) discussed both the strengths and weaknesses of the joint latent class modelling and presented an idea for a nonparametric joint modelling approach. They pointed that the joint latent class tree (JLCT) model with a tree-based approach could be a good alternative by addressing the time-invariant limitation of the JLCM. Zhang and Simonoff (2022) proposed a tree-based approach to model time-to-event and longitudinal data. This semiparametric joint latent class tree model uses time varying covariates as the splitting variables for constructing a tree in the first stage. Once the latent class has been identified for each individual at each time-point, in the second stage, it can fit the longitudinal and survival models independently. Despite the advantages in prediction performance and computational cost, the JLCT can not handle the scenario where longitudinal and survival outcomes are associated even after conditioning on the latent class.
Motivated by Garre et al. (2008) and Liu et al. (2015), we want to explore changes in time in subpopulations of subjects. Under this situation, patients may experience pathological changes such that a particular treatment regime might stop being effective for a patient and a different regime might need to be considered. This could be reflected by a change in the latent class for that patient. In order to handle this kind of particular dataset with jumping behaviours, we propose adding time-varying covariates into the latent class membership probability, which allows jumps among different subgroups.
Based on the standard JLCM, the JLCM with a time-varying membership probability is proposed, where the longitudinal and survival processes are joint via shared random effects as well as the latent class k. The proposed time-varying JLCM can relax the time-invariant restrictions of JLCM and also capture jumping behaviors for longitudinal markers. In contrast to the dynamic latent class model proposed by Reboussin and Anthony (2001), we consider a multinomial logistic regression latent class submodel. Compared with the Bayesian JLCM proposed by Andrinopoulou et al. (2020), we compute dynamic membership probabilities for each subject i across time, by introducing time-varying covariates in the latent class submodel.
A Bayesian approach is used for estimation and inference. The detailed structure of the time-varying JLCM is presented in section II. In section III, we conduct a simulation study to compare the performance of the proposed JLCM with the basic JLCM when there is jumping behaviour in the data, while in section IV, we apply our model to analyse the AIDS dataset (Goldman et al., 1996). For both simulation and real data analyses, we also investigate the performance of the proposed model in the dynamic prediction of survival probabilities. Finally, section V summarises the work in this paper and discusses future extensions of the proposed model.

II. JOINT LATENT CLASS MODEL WITH TIME-VARYING PROBABILITY
We define our reference basic JLCM is consist of the latent class submodel without time-varying covariates, the longitudinal and survival submodels which have shared random effects. Based on the basic JLCM, the structure of our proposed JLCM is shown as follows.

A. General latent class submodel with time varying probability
We consider a logistic model to describe the class membership probability, allowing this probability to vary with time t. Suppose there are N subjects, labelled i = 1, 2, . . . N and K latent classes labelled k = 1, 2, . . . K. Assume the design matrices X 1i may have time-dependent covariates or polynomial terms in time measured at time points t i = t ij | j=1,2,...mi = (t i1 , t i2 , · · · , t imi ) T , where t ij denotes the jth visit time of subject i (j = 1, 2, . . . m i ) and m i is the number of obervations for subject i. Define R ij = k as the latent class indicator, which means subject i belongs to class k at time point t ij . The probability π ijk for subject i measured at time t ij , of belonging to class k with K k=1 π ijk = 1 is described using a multinomial logistic regression model with covariate vector X 1i and associated class-specific coefficient vector ξ k , as follows:

B. Longitudinal submodel
Denoting y ijk as a repeated measurement for subject i at the jth visit in class k, the longitudinal submodel for y ij is (Troxel et al., 1998): is the covariate vector of subject i, which might contain time-varying covariates or time polynomial terms, β k is the parameter vector for class k and U ik (q × 1) are the class-specific random effects for class k. Z i (t ij ) denotes the vector of covariates at the jth visit time associated with U ik , again allowing for time-varying covariates or time polynomial terms. ijk is the class-specific error term with variance τ k . We assume ijk is independent of U ik (Liu et al., 2015). The class-specific random effects U ik can also be written as . For simplification, we assume Z k = 1, i.e. we assume the random effects across different latent classes are the same. Then, a simplified form of the longitudinal submodel is:

C. Survival submodel
Let T * i denote the true event time which may or may not be observed (i = 1, 2, . . . N ). Let C i denote the censoring time for subject i. The observed follow-up time is represented by as the censoring indicator, where∆ i takes the value 0 if censoring occurs and 1 otherwise. The time to event in each latent class can be described by a Cox frailty model with class-specific baseline hazard function and parameters (Hougaard, 1995): where λ 0k (t) denotes the baseline hazard function in class k, X 3i (t) is the covariate vector for subject i at time t and ω k is the parameter vector for class k. In addition, δ k denotes the association vector for class k associated with the random effects U i , which is considered as the strength of the association between the longitudinal measurements and survival outcomes. If δ k = 0, it means the longitudinal and survival submodels are only linked via the latent class k. We further assume that the kth baseline hazard is a step function as follows < ∞ is a partition of (0, ∞), k = 1, 2, . . . K and S k denotes the number of steps for the kth baseline hazard, labelled as s = 1, 2, . . . S k .

1) The likelihood :
We denote Ψ as the collection of model parameters, such that Ψ = (ξ, β, ω, δ, τ, λ, Σ u ) where ξ = (ξ 1 , . . . , ξ K ), β = (β 1 , . . . , β K ), ω = (ω 1 , . . . , ω K ), δ = (δ 1 , . . . , δ K ), τ = (τ 1 , . . . , τ K ), λ = (λ 01 , . . . , λ 0K ), and R i = (R i1 , . . . , R imi ). Assuming that the observed data (Y , T ,∆) are independent conditional on the latent class indicator R and the random effects U , and that (T ,∆) only depends on R imi , the likelihood contribution for the ith patient is written as Then, we obtain the joint likelihood function where 2) MCMC sampling procedure: We employ a Bayesian estimation approach to obtain the parameter estimates in our proposed JLCM. Based on the likelihood given in (4), the posterior distribution is given by: is the prior distribution of parameter vector Ψ. Assuming apriori independence, we can write π(Ψ) as the product of the prior for each parameter component. We can draw β, τ, λ, Σ ui directly from their full conditional distributions using Gibbs sampling. We set normal priors for β, ω, δ, U, ξ, obtaining conjugate posterior distributions for β. We also choose a gamma distribution as the prior for λ, an inverse Wishart prior for Σ ui and an inverse gamma distribution as the prior for τ. Details of the calculation of the full conditional densities for β, τ, λ and Σ ui , can be found in Appendix A. We also sample R ij directly from its discrete posterior distribution where the posterior probability of subject i at time j being a member of latent class k is P (R ij = k |R, Y , T ,∆, U ) = P ijk /(P ij1 + · · · + P ijK ).
However, for w, δ, ξ, U the conditional distributions have non-standard densities, making sampling from these distributions more complex. A traditional Metropolis-Hastings type algorithm (Metropolis et al., 1953;Hastings 1970) although widely used and easy to implement, can often have convergence issues. Adaptive MCMC algorithms can help to deal with this problem by automatically learning better parameter values of MCMC algorithms while they run (Roberts & Rosenthal, 2009). Hence, an adaptive MCMC method is preferred to sample the parameters/components with non-standard distributions in our proposed JLCM. Due to convergence issues with sampling τ and λ directly from their conditional posteriors, we sample them via the AM scheme, instead.
According to the Adaptive Metropolis (AM) algorithm in Haario et al. (2001), we consider ad-dimensional target distribution π(θ, τ , λ), where θ = (w, δ, ξ, U ) and define θ p = (w, δ, ξ), withd p denoting the dimension of θ p .d τ is defined as the dimension of τ andd λ is defined as the dimension of λ. We assume that the elements of θ p have individual independent standard normal priors N (0, 1) with the density function denoted as f N (0,1) and the elements of vector τ have individual independent vague inverse Gamma priors IG(0.01, 0.01) with the density function denoted as f IG(0.01,0.01) . We also set that the elements of vector λ to have individual independent uninformative Gamma priors Γ(0.01, 0.01) with the density function denoted as f Γ(0.01,0.01) .
where Ψ −(θp,τ ,λ) is the parameter vector excluding the components in θ p , τ and λ. Then, we can constuct an adaptive Metropolis algorithm with proposal distribution π(θ, τ , λ|R, Y , T ,∆, Ψ −(θp,τ ,λ) ) given at iteration m by where Σ m denotes the current empirical estimate of the covariance structure of the target distribution based on the runs so far, σ 2 is a fixed constant variance which we choose, to determine how well the adaptive algorithm is performing and α prop is a small positive constant. Gelman et al. (1997) and Roberts and Rosenthal (2001) noted that the proposal N (θ, (2.38) 2 Σ m /d) is optimal in a particular large-dimensional context. Here, we found that N (θ, (0.1) 2 Id/d) is a more suitable proposal. According to the usual Metropolis-Hastings formula (Hastings, 1970), a proposal from (θ, τ , λ) to (θ * , τ * , λ * ) is accepted with probability where θ * = (w * , δ * , ξ * , U * ) and θ * , τ * and λ * are vectors of parameters for τ and λ, which are the generated "candidate values" from Q m ((θ, τ , λ), · ). Through the MCMC sampling scheme, we can also calculate the posterior membership probability to perform modelbased classification. The posterior probability of subject i belonging to class k at time j iŝ is computed using the posterior sample of the parameters Ψ and random effects U , and π ijk is the prior probability from equation (1). The greater the value of the posterior probabilityπ ijk (Ψ) is, the more likely the subject is divided into this group.

3) Misclassification rate:
The misclassification rate, also known as the error rate, is calculated as the percentage of miss-classified counts over the whole dataset. For example, if we only have two classes, then the error rate for the JLCM without time-varying probability is computed as (Garre et al., 2008): where n M 1 and n M 2 denote the counts which are mis-classified into class 1 and class 2, respectively and N denotes the total number of subjects.
For the time-varying JLCM, we proposed a similar measure to calculte the error rate.
where N M 1 and N M 2 denote the counts which are mis-classified into class 1 and class 2 across time for all individuals and m i denotes the number of observations for subject i.

4) Predicted time to event:
Following Garre (2008) and Rizopoulos (2011), we can compute dynamic predictions of the time-to-event probabilities based on our proposed JLCM. The survival function for subject i is given by whereΨ andÛ i denote the posterior means of the parameter vector Ψ and random effects U i , and y i (t) contains the longitudinal measurements up to time t to imply survival up to this time point. (5). The posterior predictive survival function S i at t to predict the survival of a patient up to time t is given byŜ which is a mixture of survival distributions for patient i conditional on subgroup k to get survival probabilities for each subject.
To assess the prediction accuracy of the proposed model, we calculate the IPCW estimators of the area under the ROC curve (AUC) within the time interval [t, t + ∆t) for joint models of longitudinal markers and time-to-event outcomes (Blanche et al., 2015) as follows: . Notations of T i and∆ i can be found in the descriptions of the survival submodel and π i (t, ∆t) is given in equation (9).

III. SIMULATION STUDY FOR JOINT LATENT CLASS MODEL WITH A TIME-VARYING PROBABILITY A. Simulation setup
For simplicity, in our simulation study we just consider two latent classes and use baseline covariates in X 3i (fixed effects in the survival model) but a more complex model with time-varying X 3i is also possible. The fixed covariate X 1i which is generated from a standard normal distribution, is used to capture certain baseline characteristics of each subject and X 3i corresponds to a treatment indicator simulated from the bernoulli distribution (1 if Here we set a constant baseline hazard function h 0k (t) = λ 0k (i.e. we sample the baseline hazard λ 0k from an exponential distribution). We simulate survival times using For the membership probability, we set p ijk = exp(ξ k1 X 1i + ξ k2 T ime ij )/( K c=1 (exp(ξ c1 X 1i + ξ c2 T ime ij ))). p ijk is used to classify subgroups based on a multinomial distribution, which is known in advance. The error rate can be calculated by the percentage of the counts of being mis-classified into class 1 and class 2 for the whole dataset (equation (8)).
Each dataset is generated as follows:

B. Simulation results
For each simulation run, we set N as the total number of subjects and m i as the number of repeated longitudinal measurements for each subject i. Here is an example for one of 100 simulations with 5000 MCMC iterations. We have N = 50 individuals with at least 6 repeated measurements for each subject i. The sample size for this particular dataset is   Figure 1.
The mean longitudinal trace and KM plot based on the basic JLCM without time-varying changes are also shown in Figure 1. It is obvious that allowing for classification is important, to distinguish between the two different subgroups in this dataset. However, this kind of static classification can not tell us how individuals might change group membership over time. Based on our model, we select eight individuals to show how changes in the group membership are reflected based on the longitudinal trajectories.
We generate some new data under the same settings of simulation study as we have described in subsection B to illustrate group jumping behaviors and do the predictions for some specific subjects. From Figures 2 and 3, we observe three different jumping behaviours: some subjects, e.g. 28, 34 and 43 remain stable in the same group (no jumping). Others, e.g. 3 and 38, jump from one group to the other and remain there until the event or censoring occurs. Finally, other subjects, such as 7, 9 and 13 experience multiple jumps between groups. Unless we allow for the group membership probability to vary with time, we cannot capture these different behaviours across individuals. For example, for patients undergoing cancer staging (i.e. the process of determining the extent to which a cancer has developed by spreading), different medicines and treatment measures should be considered for patients with moderate or severe symptoms. Therefore, it becomes important to monitor any changes in the spread of cancer during therapy and adjust treatment plans accordingly. Our proposed model can provide this flexibility by allowing changes in group membership with time.
Tables 1 and 2 summarise the results of our simulation study, over 100 simulated datasets and using 5000 MCMC iterations with 2000 iterations as burn-in. The optimal number of classes for each model will be chosen via DIC. From Table  1, we can see an overwhelming support for the correct model, which is the proposed JLCM with time-varying probability and 2 classes, according to both DIC and the error rate. In the small number of cases when the basic JLCM is selected as the optimal model (only according to DIC), we notice that a model with 3 classes is preferred. This is to be expected, as the basic JLCM might be trying to compensate for any jumping behaviours by introducing an additional class.   For subject 7 in Figure 2, we can see that at the time of event, the probability of survival from the proposed model is lower when the event occurred (0.9677501 as opposed 0.9894454 computed via the basic JLCM). This individual exhibits jumping behaviours which implies the proposed model can predict the survival probability for this individual more accurately. For individuals who have a stable trajectory (i.e. subjects 28, 34 and 43), the proposed model provides marginally lower but similar survival probabilities to the basic JLCM, which indicates that our proposed JLCM includes the basic model as a more general JLCM. For individuals who exhibit jumping behaviours as in Figure 3, we can see that our proposed model is more conservative and provides lower predictive survival probabilities. It is also more sensitive to changes in group membership, e.g. for subject 3, the probability is more reactive to the drop in the response at the final two measurements. The dynamic prediction plots for these subjects are shown in figure C1 of Appendix C.
To compare the prediction accuracy of the proposed and basic JLCM, we calculate the IPCW estimator of the AUC within the time interval

IV. REAL DATA APPLICATION A. Data information
The Aids dataset (Goldman et al., 1996) contains both longitudinal and survival data from a randomized clinical trial looking to compare the efficacy and safety of two antiretroviral drugs in treating patients who had failed or were intolerant of zidovudine (AZT) therapy. The aids drug trial data is a data frame with 1408 observations on 9 different variables, i.e. patient, T ime, death, CD4, obstime, drug, gender, prevOI and AZT . There are 467 patients in total, distinguished by patients identifiers. Time denotes the time to death or censoring. Death is the indicator of a numeric vector with 0 denoting censoring and 1 death. The CD4 value is the CD4 cell counts in blood, which was recorded at study entry in 0, 2, 6, 12 and 18 months. CD4 cells are white blood cells which can fight infection as well as used to check the health of the immune system in patients infected with human immunodeficiency virus (HIV). Obstime denotes the time points at which the CD4 cell counts was recorded. Drug is a factor with levels ddC denoting zalcitabine and ddI denoting didanosine. Gender is a factor with levels female and male. prevOI is a factor with levels AIDS, denoting previous opportunistic infection (AIDS diagnosis) at study entry, and noAIDS, denoting no previous infection. AZT is a factor with levels intolerance and failure denoting AZT intolerance and AZT failure, respectively.
We analyse the dataset with the basic joint model and our proposed JLCM with time varying probability. Based on the proposed time-varying JLCM with K subgroups, the longitudinal submodel is as follows: CD4 ijk =β k1 + β k2 T ime ij + β k3 gender ij + β k4 prevOI ij + β k5 AZT ij + U 1i + U 2i T ime ij + ijk the class specific hazards submodel is: and the fitted latent class submodel is: U 1i is the random intercept and U 2i is treated as the random slope for time in the linear mixed-effects model, with a bivariate normal distribution U i ∼ N (q) (0, Σ ui ) and ijk are the measurement errors which follow a normal distribution: . For comparison, we obtain the basic JLCM with MLE approach from "lcmm" package as following.
The longitudinal submodel is: N (0, τ ) and the class-specific hazard submodel is: where the class-specific baseline hazrad is sample from Weinull distribution and the scale parameter λ * k > 0, the shape parameter ν * k > 0. For identification, the last latent class is treated as reference, i.e.: ξ K = 0. The fitted latent class submodel is: U 1i is the random intercept and U 2i is treated as the random slope for time in the linear mixed-effects model, with a bivariate normal distribution U i ∼ N (2) (0, Σ ui ) and ijk are the measurement errors which follow a normal distribution:

B. Estimation results and comparison with results from basic JLCM
For the real data application, we consider seven different model settings, similar to our earlier simulation study, including the proposed JLCM with 1, 2 or 3 classes and the basic JLCM with 1-4 classes. According to Table 3, the proposed JLCM with 2 classes is the optimal model. In terms of the basic JLCM, the optimal model is the one with 3 classes (which matches our findings in the simulation study). This result also agrees with the results obtained via the lcmm package, which fits the basic JLCM under a maximum likelihood approach. More details can be found in Appendix B. Table 4 summarises the jumping behaviours of individuals in the AIDS dataset. We can see that most patients remain stable in the same class and only 16% of patients have jumping behaviours. For patients with jumping behaviours, around 61% jump from the better group (higher CD4 cell count) to the lower CD4 group, which might indicate worsening of their health.  Figure 5 illustrates the differences in the mean trajectories and survival probabilities in the aids data, between a scenario with no classification and one with a static 2 group classification. We notice two distinct groups, with patients with a higher CD4 cell count having no risk of death, compared to those with lower CD4 cell count. The wide confidence bounds around the survival probabilities of the latter group, indicate that perhaps not all patients in this group are associated with lower survival probabilities and allowing for jumps between groups could provide a more accurate prediction of these probabilities. We take the first 2000 iterations as burn-in. The posterior means and sds, as well as 89% credible intervals, for each parameter in the optimal model, are presented in Table 5. Variables T ime, P revOI and intercept provide similar trends in the longitudinal model for both proposed and basic JLCM. Compared with the results of the basic JLCM with optimal K = 3 obtained using a frequentist approach via the lcmm package and shown in Table B2 in Appendix B, we notice the CD4 level decreases with time. P revOI has similar effects as variable T ime in the longitudinal submodel but P revOI accounts for a sharper decrease in the CD4 level. We also notice similar trends in terms of the effect of variables in the survival models compared to the basic JLCM. In particular, both P revOI and AZT can increase risks in both proposed and basic JLCM, however, gender and drug only seem to have a significant effect in our proposed JLCM. The δ k values in the survival function are not zero, which indicates that there exists association between the longitudinal and survival processes. Non-zero ξ k2 values indicate there is evidence of a time-varying membership probability, providing further support for our proposed model. What is more, ξ 12 is positive while ξ 22 is negative, which implies the membership probability of class 1 (2) increases (decreases) with time. In table 6, we compare the posterior predictive probabilities for each group between the basic JLCM (under both our Bayesian approach and the standard likelihood based approach from the literature) and the proposed time varying JLCM. We notice that the majority of individuals are classified in class 2 (for the time varying JLCM, this is the classification at the final time-point), with comparable percentages across the 3 models. Class 1 contains the second largest group of patients, with only very few in class 3 of the basic JLCMs. By allowing group membership to vary with time, our proposed model can capture the variability in the data without the need to introduce an additional class. Similar to the simulation study, we are interested in using the CD4 measurements to estimate the expected survival probabilities and also assess how well we are able to discriminate between patients with high risk and low risk of dying. According to equations (9) and (10), we obtain dynamic predictions based on our proposed JLCM for eight specific subjects, shown in Figures 6 and 7. We incorporate both the longitudinal trajectories and the survival probabilities on the same plot, to assess how various jumping behaviours might affect survival. The green dotted line represents the survival prediction based on the proposed JLCM with 2 classes (the best model), the blue dotted line corresponds to the survival probabilities of the basic JLCM with 3 classes (the best basic model) and the vertical brown dotted (purple dashed) lines indicate the time of event (censoring time). In the longitudinal trajectories, the black (red) circles indicate the subject is classified in class 1 (2) (based on the proposed JLCM) at the particular timepoint.
According to Figure 6, subjects 2, 90 and 60, who exhibit no jumps, have a slowly decreasing survival probability around 80% by the end of the study. For subject 60, the survival probability at the time of event based on the proposed JLCM (0.9341326) is more accurate compared to that based on the basic JLCM (0.9355598), although the estimation could be improved. For subject 188, who has a stable trajectory in class 2, with a large jump to class 1 at the last time point, the proposed JLCM seems to curb the decreasing survival probability slightly to account for the jump.
In Figure 7, we also notice that our proposed model reacts faster to the steeply decreasing trajectory for subjects 68 and 318, adjusting the survival probability accordingly. Subject 75 also exhibits a jump to the lower CD4 group, however the change in the longitudinal response is minimal and thus the survival probability is higher in this instance (but still computed more accurately compared to the basic JLCM). The survival probability for subject 341 is also similar to subject 75, as both exhibit jumps between groups but minimal change in the longitudinal response.
Compared with the basic JLCM, it is obvious that our proposed JLCM can provide more reasonable dynamic predictions as well as improve the accuracy of predicted survival probabilities. We can also see a link between the dynamic predictions from the proposed JLCM and the jumping trends exhibited in the longitudinal trajectories. The dynamic prediction plots of these eight subjects is shown in Figure C2 of Appendix C.

V. CONCLUSION AND FUTURE WORK
In this paper, we propose to extend the basic joint latent class model with a more general JLCM which allows the class membership probability to change with time. In our model, the longitudinal and survival processes are linked via both shared random effects and latent classes, and a Bayesian approach is used for estimation and inference. Through both simulation and real data analysis of the aids dataset, our proposed JLCM with time varying class membership probability, outperforms the basic JLCM, in terms of both DIC and estimation accuracy. In the dynamic prediction of survival probabilities, our model is generally more conservative than the basic JLCM and also more responsive to changes in the class membership over time. In the instances where the time to event is known, our model can predict survival probabilities more accurately, although there is still room for improvement.
In future work, the variance-covariance structure of the random effects can also be modelled simultaneously. For this particular joint model, rather than assuming the longitudinal and survival processes are joint via shared random effects, we would consider different random effects sharing the same distribution. It would then be useful to model the random effects covariance matrix because it might uncover any effects of covariates on the association between the joint models. This general joint latent class model would contain two more regression submodels for the heterogenous random covariance matrix, which can be decomposed using the Cholesky decomposition technique. To deal with outlying longitudinal measurements, a mixture distribution can also be considered for the longitudinal submodel or the normality assumption of the random errors can be relaxed, by considering a skew-normal or t distribution, leading to a more robust JLCM. (1)Sampling β k from normal distribution.
We use InverseWishart as the prior distribution for Σ ui i.e. π(Σ ui ) ∼ InverseW ishart(ν * u , S * So we can see: (4)Sampling τ k from inverse gamma distribution. Further specify the inverse Gamma prior for the variance of the measurement error τ k , i.e. we set the prior distribution for τ k : π(τ k ) ∼ InverseGamma(α * , β * ) and then we can write: π(τ k ) = β * α * /Γ(α * )(τ k ) −α * −1 exp(−β * /τ k ) (α * > 0, β * > 0). Thus, we can get the conditional distribution for τ k : If N i=1 m i I(R i = k) = 0, we sample τ k from InverseGamma(α * , β * ); if not, we sample τ k from InverseGamma( Appendix B: Estimation results for real data application on basic JLCM We use basic JLCM with maximum likelihood estimation method to do the inference using lcmm package in R. We set basic JLCM with group numbers from 1 to 4 and model it by the function "Jointlcmm". From Figure B1 in Appendix B, we can see mean trajectories as well as K-M plot of aids data for the optimal K = 3 from basic JLCM in MLE method. It is obvious that it makes sense to add time-varying membership probability into the basic JLCM for aids dataset. CD4 values are divided into different levels which indicates different health conditions for aids patients. The jumping behaviours can be considered for this specific aids data. Table B1 in Appendix B shows the summary information for basic JLCMs of aids data with various k settings. Aids data with 3 classes gets the smallest DIC value among four model settings. Basic JLCM with three subgroups is picked up as the optimal model here. Estimation results of basic JLCM with K = 3 are shown in Table B2.