Joint modeling of longitudinal and competing-risk data using cumulative incidence functions for the failure submodels accounting for potential failure cause misclassification through double sampling

Summary Most of the literature on joint modeling of longitudinal and competing-risk data is based on cause-specific hazards, although modeling of the cumulative incidence function (CIF) is an easier and more direct approach to evaluate the prognosis of an event. We propose a flexible class of shared parameter models to jointly model a normally distributed marker over time and multiple causes of failure using CIFs for the survival submodels, with CIFs depending on the “true” marker value over time (i.e., removing the measurement error). The generalized odds rate transformation is applied, thus a proportional subdistribution hazards model is a special case. The requirement that the all-cause CIF should be bounded by 1 is formally considered. The proposed models are extended to account for potential failure cause misclassification, where the true failure causes are available in a small random sample of individuals. We also provide a multistate representation of the whole population by defining mutually exclusive states based on the marker values and the competing risks. Based solely on the assumed joint model, we derive fully Bayesian posterior samples for state occupation and transition probabilities. The proposed approach is evaluated in a simulation study and, as an illustration, it is fitted to real data from people with HIV.


INTRODUCTION
A special feature of epidemiological cohort studies is that surrogate markers (e.g., markers related to disease progression) are usually collected over time along with multiple time-to-event outcomes.Such outcomes are often mutually exclusive events, referred to as competing risks (Beyersmann and others, 2011).In our motivating example from the epidemiology of human immunodeficiency virus (HIV) infection, patients receiving antiretroviral treatment (ART) can die while in care or disengage from care, which are competing risks; the number of CD4 cells is a longitudinal marker typically collected over time to keep track of HIV progression.If a patient dies or disengages from care (two competing risks), data on CD4 counts are no longer available.
Joint modeling of marker and time-to-event data has been an active research area (e.g., Rizopoulos, 2012;Hickey and others, 2018).The aim of joint modeling is 2-fold: to estimate the risk for an event conditional on an endogenous time-updated covariate (Wulfsohn and Tsiatis, 1997) and to adjust for notat-random missingness (MNAR), as most joint modeling approaches assume that missing marker data after the event are MNAR (Rizopoulos, 2012).However, the distinction between MNAR and missing at random (MAR) marker data is complex and requires further consideration (Thomadakis and others, 2019).In this article, we focus on the possibly most frequent case in which a linear mixed model (LMM) is assumed for the marker values, with the "true" marker value (i.e., predicted by both the fixed and the random effects) being included in the survival model.That is, the two processes are linked through common parameters, hence the term shared-parameter models (SPMs) that are frequently used in the literature.
Most of the research in joint modeling assumes that there is a single cause of failure.However, joint modeling of longitudinal data and competing-risk survival data has also gained attention (e.g., Hickey and others, 2018;Proust-Lima and others, 2016).In principle, competing-risk data can be analyzed through either cause-specific hazards or cumulative incidence functions (CIFs), with the latter being more direct for evaluating the prognosis of a disease.In most cases, though, the competing-risk submodels are specified in terms of the cause-specific hazards under the SPM framework (Elashoff and others, 2008;Dantan and others, 2011;Rizopoulos, 2012;Andrinopoulou and others, 2014), probably due to easier implementation, with the exception of the proportional subdistribution hazard joint model by Deslandes and Chevret (2010).In theory, a cause-specific CIF can be obtained from cause-specific hazards by integrating the product of the respective cause-specific hazard and the overall survival function.However, SPMs require additional integration over the random effects, and the overall survival function is also often approximated by numerical integration (Rizopoulos, 2012).Thus, an SPM in terms of the CIFs (or some function of them) would be more natural and could substantially reduce the computational burden of formally deriving CIF estimates based on cause-specific hazard estimates.
The literature on regression modeling of CIFs has expanded since the seminal paper by Fine and Gray (1999).An issue in such models is that the all-cause CIF should be bounded by 1, which is ignored by some approaches (Fine and Gray, 1999;Jeong and Fine, 2006;Mozumder and others, 2018).This can be dealt with by modeling the baseline asymptote for one cause-specific CIF (Shi and others, 2013), by adding a small positive number to force the survival function to be positive (Mao and Lin, 2017), or by incorporating a formal boundedness constraint in the maximization process (Bakoyannis and others, 2017).However, how to impose such a constraint in SPMs is not so clear as SPMs are defined conditionally on the random effects, and integration over the prior distribution of the random effects is required to obtain the observed data likelihood.Under the Bayesian paradigm, Gelfand and others (1992) suggested that when the constraints involve the data (as it is the case in CIF modeling), it is more natural to build the constraints into the likelihood function rather than into the prior distribution.
In HIV studies, especially in those from resource-constrained countries, under-reporting of deaths is often a major issue; patients who have actually died may have been incorrectly classified as disengaged from care, which is failure cause misclassification.To deal with this issue, a double sampling design can be used (Bakoyannis and others, 2019); that is, the true vital status of a small random sample drawn from patients reported to have disengaged from care is actively ascertained.This is performed for typically 10-20% of the patients, thus the true failure cause for the remaining patients is missing.Various methods to adjust for outcome misclassification using double sampling data have been proposed (e.g., Bakoyannis and others, 2019;Daniel Paulino and others, 2003).
In applied medical research, the progression of cohorts over time is often monitored by using mutually exclusive states defined jointly by competing-risk data and discretized continuous marker data.For example, the United Nations (UN) Joint Programme in HIV/AIDS (UNAIDS) produces various projections for the HIV epidemic using clinically relevant discrete CD4 states (i.e., [0,50), [50,100), [100,200), [200,250), [250,350), [350,500), and [500,] cells/μL) and time-to-event outcomes through the Spectrum software (Stover and others, 2019).This definition of states is clinically meaningful in HIV infection as ART was initiated in the past based on some of the previous CD4 cutoff points.In chronic kidney disease longitudinal studies, glomerular filtration rate (GFR) is a key surrogate of kidney function, which is typically collected over time.Survival states (death and end-stage kidney disease) together with discrete states based on the GFR levels have been used in this context (Hu and others, 2012).Hu and others (2012) proposed an interesting approach estimating multistate probabilities, but they relied on a two-stage approach without accounting for potential failure cause misclassification.
In this article, we propose a unified and flexible approach to jointly model a continuous marker over time and competing-risk data using CIFs for the failure submodels, accounting for failure cause misclassification using double sampling.Inference is obtained through a Bayesian procedure, and based on the postulated model, we also derive posterior samples for multistate probabilities defined jointly by marker and competing-risk data.In Section 2, we describe the structure of the proposed models, and we describe the extensions required to account for potential failure cause misclassification in Section 3. In Section 4, we derive a procedure to obtain posterior samples for multistate probabilities.In Section 5, a simulation study is performed to evaluate the performance of the proposed methodology, while, in Section 6, we fit the examined models to real data.Finally, in Section 7, we present concluding remarks and discuss limitations and possible extensions.

Marker model
For the longitudinal marker model, we use an LMM model of the form where x i (t) and z i (t) denote the fixed-and random-effects design matrices at time t, respectively, and β and b i ∼ N (0, D) denote the fixed and random effects, respectively.Furthermore, i (t) ∼ N (0, ω −1 ) denotes the within-subject residuals with ω being the within-subject precision.As usually assumed in the joint modeling literature (Rizopoulos, 2012), the "true" marker value at time t is defined as m i (t) = x i (t)β + z i (t)b i and the history of "true" values up to time t are denoted by M i (t) = {m i (s) : 0 ≤ s ≤ t}.The vector of the observed marker values on the ith subject is denoted by Y i = {y i (t i1 ), y i (t i2 ), . . ., y i (t in i )}, where t i1 , . . ., t in i are the observation times and n i is the number of observed marker values on subject i.The corresponding design matrices for the fixed and random effects at times t i1 , . . ., t in i are denoted by X i and Z i , respectively.The marker model parameters are θ L = (β , vech(D) , ω) and b = (b 1 , . . ., b N ).

Competing-risk survival models
Let T i be the time to the first occurring event for the ith individual and K i ∈ {1, 2, . . ., K} be the corresponding failure cause.We propose to model the CIFs for all causes simultaneously conditional on the "true" marker values, that is, where w ik is a vector of baseline covariates for cause k and individual i and θ sk is the parameter vector for cause k.Note that (2.1) depends on β and b i through the history of the "true" marker values, M i (t).Since all CIFs are modeled simultaneously, the all-cause CIFs should be bounded by 1 at each failure time.To account for that, we assume that all CIFs increase over time up to a common point τ i at which the all-cause CIF approaches 1 and they reach a plateau thereafter, where F M ik {t|M i (t), w ik ; θ sk } is a certain parametric model for the CIF of cause k conditional on the "true" marker values, M i (t), and some baseline covariates, w ik .To be more formal, τ i depends on the values of (β , θ s , b i ), that is, where θ s = (θ s1 , θ s2 , . . ., θ sK ) .In other words, the support of the distribution of the survival time The motivation for assuming (2.2) becomes clearer when considering the survival likelihood conditionally on the random effects under noninformative right censoring (Jeong and Fine, 2006); we can only observe T i = min(T i , C i ) and δ ik = I (K i = k), where C i is the hypothetical censoring time for the ith individual, δ ik is the corresponding failure indicator for cause k, k = 1, 2, . . ., K, and δ i = K k=1 δ ik is the overall failure indicator.As a convention, we assume that K i = 0 denotes right censoring.In this case, the survival likelihood is equal to where , w ik ; θ sk } is the overall survival function, and w i = (w i1 , w i2 , . . ., w iK ) .
For some specific set of parameter values, (β, θ s , b i ), suppose that the assumed model yields an all-cause CIF evaluated at the observed survival time, T i , greater than or equal to 1, that is, . By the definition of τ i , it is implied that T i does not lie within the support of T i given the parameter values, (β, are equal to zero, ensuring that the likelihood function is equal to zero.Thus, (2.3) is equivalent to including the model-based CIF, F M ik {T i |M i (T i ), w ik ; θ sk }, and its derivative along with the indicator function 3).We propose to model the CIFs using the class of generalized odds rate transformation models (e.g., Bakoyannis and others, 2017;Jeong and Fine, 2006) where model SPM-1 is a proportional subdistribution hazard model (Deslandes and Chevret, 2010;Fine and Gray, 1999), since λ M ik {t|m i (t), w ik ; , where λ M ik {t|m i (t), w ik ; θ sk } is the assumed subdistribution hazard function, whereas SPM-2 reduces to SPM-1 as c k → 0 (Jeong and Fine, 2006).B k (s) denotes a B-splines basis matrix at time s, with associated parameter ψ k , and γ k measures the effect of the baseline covariates on the kth CIF.The parameters α k , k = 1, . . ., K, correspond to the effects of the "true" marker values on the CIFs and indicate the level of association between the marker and the survival submodels (referred to as the association parameters).Under SPM-1, exp(α k ) denotes the kth subdistribution hazard ratio at time t, resulting from one unit increase in m i (t) at the same time point.The interpretation of the parameters of SPM-2 does not seem so appealing in general, but assuming c k = 1 a proportional rate of odds increase model is implied, that is, . Thus, SPM-2 would be useful when the SPM-1 model does not provide good fit to the data.We assume that the parameters c k are known, as trying to estimate them can lead to nonidentifiability issues (Bakoyannis and others, 2017).When this is not true, a feasible, although not optimal approach to select the link function parameters is to perform a grid search over plausible values of (c 1 , c 2 , . . ., c K ) and select the one that optimizes some model comparison criterion.

Bayesian inferential procedures
The observed data likelihood of the model requires multidimensional integration over the random effects.Since such integration is challenging, especially given the constrained space due to (2.2), we rely on a Bayesian inferential procedure based on a Markov chain Monte Carlo (MCMC) algorithm.Letting θ = (θ L , θ s ) be the entire parameter vector of the models and where f (θ ) is the prior distribution of the parameters.The integrals involved in the CIFs can be accurately approximated by Gauss-Legendre rules with 30 nodes.A Normal prior distribution, N (μ 0 , C 0 ), is used for β, a Gamma(λ 1 , λ 2 ) for ω and a Normal, N (μ s 0 , C s 0 ), distribution for θ s .For the covariance matrix of the random effects, D, we assumed the Inverse-Wishart IW (A, df ) distribution.To update the parameter values, we used Gibbs steps wherever possible and Metropolis-Hastings steps for the remaining parameters.Further details on the MCMC algorithm are presented in Section S1 of the Supplementary material available at Biostatistics online.It needs to be emphasized that any proposed value leading to an all-cause model-based CIF greater than or equal to 1, that is, K k=1 F M ik {t|M i (t), w ik ; θ sk } ≥ 1, is rejected as the posterior ratio is equal to zero.Thus, calculation of τ i (β, θ s , b i ) is not required within the MCMC algorithm.

INFERENCE UNDER POTENTIALLY MISCLASSIFIED CAUSES OF FAILURE
When the true failure cause, K i , is not observed for all individuals, we assume that a cause of failure, Ki , is always reported although potentially misclassified (i.e., K i = Ki ).Let π jk (D misc,i ) = Pr( Ki = j|K i = k, D misc,i ; θ misc ) be the probability of observing failure cause j given the kth true failure cause and the history of the observed information up to T i (including but not limited to the observed marker values), D misc,i , with θ misc being the associated parameter vector.Note also that π kk (D misc,i ) is the probability of correctly classifying cause k, whereas K j=1 π jk (D misc,i ) = 1, for any k ∈ {1, 2, . . ., K}.Moreover, we assume that (i) noninformative right censoring (e.g., administrative censoring) is correctly classified, that is, K i = 0 ⇔ Ki = 0 and (ii) the true failure cause is known in a small random sample of individuals, leading to a double sampling design (e.g., Bakoyannis and others, 2019).In this context, the observed data are where R i is an indicator function of the ith individual being doubly sampled.We make the MAR assumption for the missing failure causes, that is, the probability of being in the double sample depends on the observed data, but not on the missing true failure cause and the random effects.This implies that the true failure cause can be validly predicted based on the observed data.As shown in Section S2 of the Supplementary material available at Biostatistics online, the true failure cause probabilities conditionally on the observed data, For individuals that are not doubly sampled, we observe only Ki , which may be different from K i .To deal with this issue, one can use data augmentation (Tanner and Wong, 1987), augmenting the observed likelihood for individuals who have failed from any event but have not been included in the double sampling by the unobserved true failure causes, K i .Letting I mis be the indices of individuals that have failed from any event but are not doubly sampled (I mis ≡ {i : Ki = 0 & R i = 0}), the augmented posterior distribution of all unknown quantities is equal to where the full survival likelihood, f {T i , K i , Ki |M i (T i ), w i , D misc,i ; θ s , θ misc }, has been factorized as the product of f {T i , K i |M i (T i ), w i ; θ s } and Pr( Ki |K i , D misc,i ; θ misc ), with δij = I ( Ki = j).Note also that Pr( Ki = 0|K i = 0, D misc,i ; θ misc ) = 1, for the right censored individuals.The following algorithm outlines the modified MCMC procedure to account for misclassification: choosing appropriate initial values θ (0) , b (0) , {K (0) i : i ∈ I mis }, θ (0) misc , meeting the likelihood constraints for all individuals, for l = 1, 2, . . ., L: • Update (θ (l−1) , b (l−1) ) to (θ (l) , b (l) ) according to the posterior distribution f (θ , b|{K (l−1)   i , that is, the posterior distribution of the parameters of main interest, with the missing failure causes being equal to their current values.The MCMC algorithm for fully observed failure causes can be used.
• Update θ (l−1)  misc to θ (l)  misc according to f (θ misc |{K (l−1) Thus, data augmentation results in a simple MCMC scheme as, conditionally on the true failure causes, K i , the posterior distribution of (θ, b) is independent of θ misc and it has the same form as with the case of fully observed failure causes.The simplicity, though, of data augmentation comes at a price as repeatedly sampling K i may slowly converge towards its limit distribution.
Although the misclassification probabilities do not directly depend on the random effects, they are conditional on the true failure cause, with the missing failure causes imputed based on the CIFs, which in turn depend on the random effects.The simplest approach is to not include any covariate information in π jk (D misc,i ), that is, D misc,i assumed to be an empty set.Then, a natural choice for the prior distributions of the π jk 's would be the Dirichlet (Beta for K = 2) distributions, leading to conditional conjugacy.More generally, to model π jk (D misc,i ) conditional on the observed information (i.e., D misc,i ), multinomial logistic regression could be used.More information is provided in Section S3 of the Supplementary material available at Biostatistics online.
To compare the fit of the models, we adopted the marginalized version (Quintero and Lesaffre, 2018) of the deviance information criterion (DIC), which requires integration of the random effects.Integration was performed through Monte Carlo integration or importance sampling when the former fails (Section S4 of the Supplementary material available at Biostatistics online).

POSTERIOR INFERENCES FOR POPULATION-AVERAGED CIFs AND MARKER STATES
To describe the cohort evolution over time, states defined by marker data and clinical outcomes are often used.A pragmatic approach to do so is to discretize the marker values into nonoverlapping intervals {[s 0 , s 1 ), . . ., [s J −1 , s J )} and define mutually exclusive states based on clinical events and (discretized) marker data.If the focus of the analysis lies in describing the "true" biological process, as often is the case in the joint modeling literature, states may be defined in terms of the "true" marker values, that is, for Progression of the whole cohort can be easily monitored by a series of estimated multistate probabilities Pr{m i (t) ∈ S h , T i > t|w i ; θ }, h = 1, . . ., J and Pr(T i ≤ t, K i = k|w ik ; θ ), k = 1, . . ., K, through a multistate probability plot.The first quantity, referred to as latent marker state probability, expresses the probability of being event free and having true marker values in S h .The second expression, that is, Pr(T i ≤ t, K i = k|w ik ; θ ), is the population-averaged CIF for a particular cause.To get better insight into the dynamics of the processes, one may be also interested in transitions between states.In real-life applications (e.g., Stover and others, 2019), for simplicity, transitions are often defined from baseline states.Letting p g (0) = Pr{m i (0) ∈ S g ; θ L }, it can be easily shown that Inference on (4.6) and (4.7) involves two distinct problems (i) approximation of the integral over the random effects and (ii) accounting for the variability in θ .

Estimation procedure
We initially describe the estimation of (4.6) and (4.7) for any given θ .Specifically, (4.6) can be approximated by drawing samples {b for b i from the N (0, D) distribution under the linear constraint m i (0) ∈ S g , which can be carried out, among many other options (e.g., Gibbs sampling), through Hamiltonian Monte Carlo (Pakman, 2015).Specifically, where m (j) Similarly, it follows that, after multiplying and dividing (4.7) by Pr{m i (t) ∈ S h , m i (0) ∈ S g ; θ }, (4.7) can be approximated using samples {b where m igh and M (j) and Pr{m i (0) ∈ S g ; θ } can be easily computed using the cumulative distribution function of the (bivariate) Normal distribution.Due to (2.2), it should be noted that if ig ), thus calculation of the upper bound is required only for the random draws that do not fulfill the boundedness constraint.

CIF estimates conditional on observed marker states
In a clinical application, estimating the population-averaged CIF conditional on the observed marker state could be valuable for making projections about the future cohort evolution.Thus, CIFs given observed baseline state, Pr{T i ≤ t, K i = k|y i (0) ∈ S g , w ik ; θ }, g = 1, . . ., J , could be of interest.Similarly, one may be also interested in CIFs conditional on being in certain observed states at specific time points.In this case, it would be reasonable to also condition on survival up to the last time point and the baseline state, that is, Pr{T i ≤ t, K i = k|T i > s, y i (0) ∈ S g , y i (s) ∈ S h , w i ; θ }, for 0 ≤ s < t and g, h ∈ {1, 2, . . ., J }.Such estimates could be useful for identifying certain subsets of the population who are event free and at high risk for developing any of the events.Estimation of these probabilities is outlined in Section S5 of the Supplementary material available at Biostatistics online.

SIMULATION STUDIES
A simulation study was carried out to evaluate the performance of the proposed methodology under certain conditions.Marker data were generated using a piece-wise linear LMM, and i (t) ∼ N (0, ω −1 ).Thus, the population slopes are β 1 , β 2 , and β 3 when t ∈ [0, 1), t ∈ [1, 5), and t > 5, respectively.Measurements were assumed to be collected biannually and the maximum study duration was assumed to be 10 years.We assumed two competing risks with K i = 1, 2 corresponding to death and disengagement from care, respectively.Two scenarios regarding the competing-risk data were considered: survival data were simulated based on both the SPM-1 model and the SPM-2 model with c 1 = c 2 = 1, with respective equations , where k = 1, 2 and w i is a binary baseline covariate following the Bernoulli distribution with probability 0.5.An independent right censoring mechanism was also applied using C i ∼ min(U i , 10), where U i ∼ Exp(0.025)(i.e., the exponential distribution with rate = 0.025), leading to around 50% censoring rate.The reported failure cause was generated conditional on the first and the second true failure cause with probabilities π 11 = 0.75 and π 22 = 0.90, respectively, whereas 15% of the noncensored observations was included in the double sample.For each simulation scenario, we simulated 500 data sets including N = 2000 individuals per data set.
Under each of the two scenarios for the survival submodels, we fitted the proposed model using both the SPM-1 and SPM-2 parameterizations, with the marker model being correctly specified in the fitted models.We also examined the performance of the proposed approach in deriving inferences on the quantities described in Section 4. These estimates were produced at times 0, 2, 4, 6, 8, and 10 years.For each fitted model, we calculated the marginal DIC criterion.Its performance in correctly identifying the true model was assessed by recording the proportion of time the true model was chosen under both scenarios.Further details on the simulation study are provided in Section S5 of the Supplementary material available at Biostatistics online.
To assess model performance, we present the bias, the Monte Carlo standard deviation, the average model-based standard error, and the empirical coverage probability of credible intervals.Since parameters γ k and α k do not have the same interpretation under SPM-1 and SPM-2, we did not provide bias and coverage probability results for γ k and α k when the fitted model was misspecified.The results under the SPM-1 and SPM-2 scenarios are presented in Tables 1 and 2, respectively.The fixed-effect estimates were approximately unbiased for both models under the two scenarios (bias from −0.010 to 0.007), while the coverage probabilities were close to the nominal level (from 93.2% to 95.8%).The estimates of γ k and α k were nearly unbiased along with approximately 95% coverage probabilities when the fitted model coincided with the true data generating mechanism.For the misclassification parameters (i.e., π 11 and π 22 ), the bias was small, with decent coverage rates.The DIC criterion had a moderate ability to identify the correct model as it selected the true model 75.0% and 58.2% of the time under the SPM-1 and SPM-2 scenarios, respectively.Its ability to correctly identify the true model substantially increased to 88.0% and 82.0%, respectively, when a simulation study including 8000 individuals with 50 replications was performed (further results not shown).Focusing on the population CIF estimates, both models yielded estimates with negligible bias along with adequate empirical coverage probabilities, while the Monte Carlo standard deviation was close to the average model-based standard error.Results for the remaining quantities presented in Section 4 were roughly similar and are presented in detail in Section S7 of the Supplementary material available at Biostatistics online.
To investigate the added value of double sampling, we performed an additional simulation study in which only the reported failure causes are used in the model.The main conclusion was that the estimated association parameters can be seriously biased (Tables S9 and S10 of the Supplementary material available at Biostatistics online).As expected, when the true failure cause was not misclassified (Tables S11 and S12 of the Supplementary material available at Biostatistics online), the model performance was satisfactory.Results from 500 replications with each data set including 2000 individuals.The "true" marker evolution was based on linear splines with knots at 1 and 5 years since baseline and it was correctly specified in the fitted SPM-1 and SPM-2 models."True" denotes the true parameter values; "Median" the mean of posterior medians over the 500 replications; "Bias" the mean bias for posterior median estimates; "ASD" the average posterior standard deviation, "MCSD" the empirical Monte Carlo deviation of estimates and "Coverage" the empirical coverage probability (%) of posterior credible intervals.Results from 500 replications with each data set including 2000 individuals.The "true" marker evolution was based on linear splines with knots at 1 and 5 years since baseline and it was correctly specified in the fitted SPM-1 and SPM-2 models.True, denotes the true parameter values; Median, the mean of posterior medians over the 500 replications; Bias, the mean bias for posterior median estimates; ASD, the average posterior standard deviation; MCSD, the empirical Monte Carlo deviation of estimates; Coverage, the empirical coverage probability (%) of posterior credible intervals.

APPLICATION
The proposed methodology was applied to data from the EastAfrica International Epidemiologic Databases to Evaluate AIDS (IeDEA) Regional Consortium.We aimed to jointly estimate the CD4 evolution after ART initiation and the CIFs for death and disengagement from care.To adjust for potential death underreporting, we incorporated information from double sampling, that is, a random sample from disengaged patients whose true vital status was ascertained by tracing these patients in the community.In this application, misclassification can be safely assumed to be unidirectional, that is, a true death can be incorrectly classified as a disengagement from care but a true disengagement from care cannot be misclassified as an observed death.Thus, assuming nondifferential misclassification, there is actually only one misclassification parameter π 11 , with K i = 1 and K i = 2 corresponding to death and disengagement from care, respectively.
To illustrate the proposed methodology, we used a 60% random sample of women aged from 35 to 45 years (the most frequent covariate pattern in the data), leading to 8005 patients included.3275 (40.9%) and 273 (3.4%) disengagements from care and deaths were reported, respectively, whereas 4457 (55.7%) were event free at the end of the follow-up (noninformative right censoring).In 443 (13.5%) patients out of those who were reported as disengaged from care, the true vital status was ascertained through double sampling, and among them, there were 80 (18.1%) hidden deaths.As mortality among doubly sampled individuals was very high, disengagement from care should not be treated as a noninformative censoring event.The median (interquartile range) number of CD4 observations per individual, CD4 count at ART initiation, and follow-up time were 2.0 (1.0, 5.0), 163 (80, 264) cells/μL, and 1.3 (0.4, 2.7) years, respectively.
To model the CD4 evolution, we used an LMM on the square root scale using natural cubic splines of time with knots at 0.55, 1.25, and 2.35 years since ART initiation for both the fixed and random effects.The psition of knots was based on the quartiles of the measurement times.Baseline CIF levels for death and care disengagement were modeled through cubic B-splines with 2 and 3 knots, respectively.The model for each cause was selected among SPM-1 and SPM-2 with c k ∈ {0.25, 0.5, 0.75, . . ., 2} by optimizing the DIC, assuming nondifferential misclassification, π 11 .The DIC criterion was optimized at c 1 = 0.5 for death and the subdistribution hazard model (SPM-1) for disengagement from care.We also examined whether π 11 depends on the time since ART initiation and the last observed √ CD4, but the effects were nonsignificant, thus, the final model assumed nondifferential misclassification.For comparison, we also fitted the corresponding model without misclassification.Details on the prior distributions used and the specification of the MCMC algorithm are presented in Section S8 of the Supplementary material available at Biostatistics online.
The main results are presented in Table 3.There is substantial underestimation of mortality, as only 29.79% of the estimated "true" deaths were reported.The fixed-effect estimates and the effects of the "true" marker value on mortality were roughly similar between the two models.Of note, the corresponding effects on the CIF of disengagement from care were discordant, that is, an increase of m i (t) was associated with a significantly greater subdistribution hazard for disengagement from care when death misclassification was accounted for through double sampling, whereas the model using the observed failure causes implied no effect.This finding may be attributed to the considerable proportion of hidden deaths among the patients flagged as lost to clinic, as deceased patients had lower "true" marker values on average.From now on, we focus only on the model accounting for misclassification.The estimated CD4 evolution is presented in Figure 1A 1 .After an initial very rapid increase from ART initiation to around 6 months, CD4 cell counts continued increasing but at a lower rate, until reaching a plateau at about 6 years since ART initiation.The estimates from the proposed model were similar to those obtained by a corresponding LMM (Figure 1A 1 ).The mean evolution was based on natural cubic splines of time with knots at 0.55, 1.25, and 2.35 years since ART initiation while the random-effects specification was based on a random intercept and slope structure."Median," "SD," "LB," and "UB" denote the posterior median, standard deviation, 2.5% and 97.5% quantiles, respectively.sHR, denotes the subdistribution hazard ratio; π 11 , probability of correctly classifying a death.
The results for the estimated population-averaged CIFs are presented in Figure 1A 2 .Also shown is the corresponding CIFs ignoring potential misclassification, that is, Pr(T i ≤ t, Ki = 1; θ , θ misc ) = Pr(T i ≤ t, K i = 1; θ )π 11 and Pr(T i ≤ t, Ki = 2; θ , θ misc ) = Pr(T i ≤ t, K i = 1; θ )(1−π 11 )+Pr(T i ≤ t, K i = 2; θ).These estimates are in close agreement with the correspondingAalen-Johansen estimates, implying that the model is flexible enough to model the observed patterns of events.As expected, ignoring misclassification led to underestimation of mortality and overestimation of the risk for disengagement from care.In Figure 2, we present multistate probabilities for all states simultaneously (CD4 states, death, and disengagement from care).By 7 years since ART initiation, we estimated that 14.2% (95% CI 12.6-16.0)had died and 58.3% (95% CI 55.5-61.1)had disengaged from care.The corresponding results by baseline CD4 state are presented in Figure S5 of the Supplementary material available at Biostatistics online.The CIF of death at 7 years for those starting at m i (0) < √ 50 was 30.9% (95% CI 27.6-34.3),remarkably higher than that of the remaining baseline CD4 states.Given m i (0) < √ 50, the transition probability to m i (5) ≥ √ 500 at 5 years while being event free was low 6.7% (95% CI 5.6-7.9),whereas the corresponding probability for those with √ 50 ≤ m i (0) < √ 100 was 12.5% (95% CI 11.3-13.6).Among those with 100-200 observed baseline CD4 cells/μL who were alive and progressed to 200-250 or > 500 CD4 cells/μL at 1 year since ART initiation, the conditional probabilities of dying within the next year were 2.8% (95% CI 2.1-3.5) and 1.1% (95% CI 0.8-1.4),respectively.

DISCUSSION
In this article, we proposed a flexible and unified class of models to jointly model a normally distributed marker over time and competing risks using CIFs for the survival submodels, with inference on model parameters obtained through a hybrid MCMC algorithm.The proposed models assume that the CIFs depend on the "true" marker value over time, m i (t), thus the association between the marker and survival processes is induced via the random effects.Hence, the proposed models lie within the family of SPMs.Most competing-risk SPMs rely on cause-specific hazards; but CIF estimates may be of particular interest when the focus lies on prognosis.Though it is feasible to derive CIF estimates based on estimated causespecific hazards, it requires complex integration, being particularly challenging in joint models.In contrast, under our proposed approach, the effects on the CIFs are described in a direct and straightforward way.To model the link functions, we used the generalized odds rate transformation, with the proportional subdistribution hazards model (Deslandes and Chevret, 2010;Fine and Gray, 1999) being a special case.Due to potential failure cause misclassification in our motivating example, we extended our methodology by incorporating information from doubly sampled patients, that is, a random sample from patients to whom a gold standard diagnostic procedure was performed.Accounting for misclassification, based solely on the joint model, we also estimated multistate probabilities jointly defined by marker data and competing risks.A simulation study was carried out to examine the performance of the methodology, indicating that the model performance is satisfactory when the marker trajectory and the association structure between the marker and the CIFs are correctly specified but the link function is misspecified (using two scenarios for c k ).The proposed models were also fitted to data from the IeDEA study using CD4 count data from ART initiation until the occurrence of death or disengagement from care.To reduce computation time, a 60% random sample was used.Ignoring double sampling led to seriously underestimated mortality, whereas it implied no effect of the "true" CD4 count on the risk for disengagement from care, but after adjusting for misclassification, moderate, but statistically significant, positive correlation was found.We suppose that the latter discrepancy could be at least partly explained by the considerable proportion of deaths among those observed to disengage from care.
One important issue when specifying models for CIFs is that the all-cause CIF should be bounded by 1 at each failure time.When there are no random effects, this can be dealt with in the maximization process, but how to address this issue in the presence of random effects has not been resolved in the literature.Our model assumed τ i (β, θ s , b i ) as the upper bound of the survival time, which mathematically led to zero likelihood when the constraint is violated.This was equivalent to introducing an indicator function in the likelihood in the parameter estimation process.However, to further derive population-averaged CIFs and marker state probabilities, integration over the random effects is required, and thus, CIFs should be evaluable at any random-effect value drawn from its prior N (0, D).Thus, having an explicitly defined model for the CIFs accounting for the constraints, population-averaged quantities can be estimated directly.
Our approach of multistate modeling differs from standard approaches (e.g., Putter and others, 2007) in which states are assumed to be directly observed and usually rely on the Markov assumption.In contrast, under our approach, marker states were not assumed to be directly observed, with the computations being solely based on the assumed joint model by formally deriving posterior samples for multistate probabilities.Our approach has also substantial differences from the work by Hu and others (2012); they proposed a two-stage approach where marker trends are first estimated using subject-specific regression models and then marker states are evaluated by averaging over individuals.Thus, the effects of the marker on the competing risks are not modeled explicitly, whereas individuals with highly irregular visit times or just one marker measurement may cause additional difficulties.
In this work, we have adjusted for failure cause misclassification through double sampling.The issue of missing failure cause in a joint modeling setting has been addressed before (Sheikh and others, 2021).However, there are major differences with our approach as Sheikh and others (2021) did not consider reported failure causes and used cause-specific hazards.
The proposed methodology relies on parametric assumptions, thus, as always with parametric approaches, certain assumptions may not hold.Among all model assumptions, though, the ones that are most difficult to verify are perhaps those related to missing data mechanisms, for example, the proposed models, lying within the general class of SPMs, assume that missing marker data after the first occurring event are MNAR.However, the question of missing data being MAR or MNAR is complex and probably depends on the richness of the design.Death is often considered to cause MNAR marker data as it corresponds to underlying disease progression that is unlikely to be reflected in the observed marker data, measured at prespecified time points; the nature of the dropout mechanism due to disengagement is less clear though.It has been shown that if missing marker data are MAR, fixed-effect estimates from specific SPMs are susceptible to bias (Thomadakis and others, 2019).In our application, we feel that this is unlikely as the fixed-effects estimates from the proposed model were in line with the corresponding ones from the LMM.
There are some extensions that could be incorporated into the proposed methodology: for example, it would be interesting to provide dynamic survival predictions (Rizopoulos, 2012) or use more flexible forms for the association between the marker and the CIFs.Moreover, some aspects of the proposed methodology may require further consideration.For example, it would be interesting to evaluate the model performance under misspecification of the misclassification model, different true parameter values, and different percentages of doubly sampled individuals.
To sum up, we have proposed a flexible class of SPMs to jointly model a normally distributed marker and competing risks using CIFs in the survival submodels, extended to account for potential failure cause misclassification.As most approaches in the literature rely on cause-specific hazards, our proposed approach can be a useful alternative when the focus is on identifying risk factors for the risk of an event.
SUPPLEMENTARY MATERIAL Supplementary material is available online at Biostatistics online and software in the form of R code is available at https://github.com/cthomadak/JointModelCumInc.

Fig. 1 .
Fig.1.Estimated CD4 evolution and population-averaged CIFs based on the proposed model assuming c 1 = 0.5 for death and the proportional subdistribution hazard model (SPM-1) for disengagement from care, taking into account the double sampling data, applied to East Africa IeDEA data.Shades in gray show pointwise credible intervals.A 1 : estimated CD4 evolution over time since ART initiation (CD4 counts up to 800 cells/μL are shown).A 2 : populationaveraged CIFs for death and disengagement from care along with the corresponding CIFs for an observed death and disengagement from care.

Fig. 2 .
Fig.2.Stacked multistate probability plot of latent marker states and competing risks of death and disengagement from care over time since ART initiation based on the proposed model assuming c 1 = 0.5 for death and the proportional subdistribution hazard model (SPM-1) for disengagement from care, taking into account the double sampling data, applied to East Africa IeDEA data.The corresponding state occupancy probabilities are visualized through the difference between two adjacent curves with different shades of gray.

Table 2 .
Simulation study results from fitted SPM-1 and SPM-2 models when the data are generated by the SPM-2 model under failure cause misclassification

Table 3 .
Results from the proposed model assuming c 1 = 0.5 for death and the proportional subdistribution hazard model (SPM-1) for disengagement from care applied to East Africa IeDEA data