A latent class approach for joint modeling of a time-to-event outcome and multiple longitudinal biomarkers subject to limits of detection

Abstract

Multiple biomarkers on different biological pathways are often measured over time to investigate the complex mechanism of disease development and progression. Identification of informative subpopulation patterns of longitudinal biomarkers and clinical endpoint may assist in risk stratification and provide insights into new therapeutic targets. Motivated by a multicenter study to assess the inflammatory markers of sepsis in patients with community-acquired pneumonia, we propose a joint latent class analysis of multiple biomarkers and a time-to-event outcome while accounting for censored biomarker measurements due to detection limits. The interrelationship between biomarker trajectories and clinical endpoint is fully captured by a latent class structure, which reveals the subpopulation profiles of biomarkers and clinical outcome. The estimation of joint latent class models becomes more complicated when biomarkers are subject to detection limits. Based on a Metropolis–Hastings method, we develop a Monte Carlo Expectation–Maximization (MCEM) algorithm to estimate model parameters. We demonstrate the satisfactory performance of our MCEM algorithm using simulation studies, and apply our method to the motivating study to examine the heterogeneous patterns of cytokine responses to pneumonia and associated mortality risks.

Keywords

Joint modeling latent class limit of detection Monte Carlo EM multiple biomarkers

1 Introduction

In biomedical and clinical research, multiple biomarkers are often measured over time to investigate the complex mechanism of disease development and progression. For example, CD4+ cell counts and plasma HIV RNA viral load have been studied as prognostic markers of HIV infection.^1,2 Comprehensive neuropsychological tests were used to evaluate the underlying cognitive level in dementia study.³ It is of great interest to determine the interrelationship between multiple biomarkers and their joint effects on the clinical outcomes. Our work was motivated by the Genetic and Inflammatory Markers of Sepsis (GenIMS) study,⁴ a multicenter cohort study of patients admitted to the emergency departments with community acquired pneumonia (CAP). Multiple biomarkers were measured on the patients admitted to the hospital daily in the first week, and weekly thereafter while they were still in the hospital. One of the primary goals of GenIMS study is to examine the activation of pro- and anti-inflammatory biomarkers after pneumonia and their roles in predicting subsequent adverse events. The biomarker analysis of GenIMS study is challenging because some biomarkers are censored at detection limits due to the low sensitivity of the bioassays used in the study. In particular, two inflammatory biomarkers, interleukin-6 (IL-6) and interleukin-10 (IL-10), were censored at lower detection limits with moderate to heavy censoring rates that increase over time, necessitating the use of appropriate statistical methods to accommodate this level of censoring.

Joint models have been widely used to explore the association between biomarker trajectories and a clinical endpoint. The relationship between the longitudinal biomarkers and primary endpoint can be modeled through shared random effects⁵ or latent classes.⁶ The shared random effects approach assumes the population is homogeneous and links the two processes using a specific structure of random effects, while the joint latent class approach assumes a heterogeneous population and uses latent classes to fully capture the associations between two processes. For the GenIMS study, the investigators hypothesized that patients with CAP may have different patterns of cytokine responses to pneumonia that are associated with different mortality risks. In the previous report of GenIMS study,⁴ the biomarker analysis was conducted in two stages. First, inflammatory biomarkers IL-6 and IL-10 were analyzed separately to identify the informative subgroups and then the two group memberships were used as covariates to predict the clinical outcome. This approach may not be efficient since the correlation between biomarkers was ignored in the subgroup identification. The results may also be biased because the group memberships are treated as if they are observed. Furthermore, the clinical outcome is not taken into account to determine the subgroups. In order to incorporate both longitudinal biomarkers and clinical outcome into subgroup identification while accounting for the correlation between biomarkers, we take a joint latent class modeling approach which allows exploring the heterogeneity of cytokine responses as well as the clinical outcome.

Several methods have been developed to account for the associations between multiple biomarkers in the joint modeling framework. In the joint models with shared random effects, the interrelationship between multiple biomarkers is often reflected by between-marker covariances of random effects in a multivariate mixed model.^7–10 However, this approach has rarely been used for the joint latent class models to our best knowledge. To handle multiple cognitive tests in a joint latent class model, model, Proust-Lima et al.³ considered a common longitudinal latent process, where each test was linked to this process through a flexible nonlinear measurement model. This idea is well suited to the study of cognitive decline and dementia when comprehensive tests were designed to estimate the underlying cognitive level, but it may not be very applicable for other studies such as GenIMS study when multiple biomarker profiles actually represent different pathways and may not share an underlying latent process. In this paper, we develop a joint latent class model to accommodate the correlations between biomarkers via a covariance structure of marker-specific random effects, similarly to the joint models with shared random effects.

As mentioned previously, the censored biomarker data in the GenIMS study requires appropriate handling. Statistical methods based on maximum likelihood estimation have been proposed to handle the censored longitudinal biomarkers in mixed models^2,11–13 and joint models with shared random effects.^10,14 Nevertheless, there is still no readily available algorithm and software to accommodate two censored biomarkers and a survival outcome in the shared random effects models. The common strategy to address censored data is to incorporate the contribution of left-censored observations by using a cumulative distribution function in the likelihood function. We adopt this strategy to handle left-censored biomarkers in joint latent class models. Since the random effects for longitudinal biomarkers are not shared between submodels, most computational methods for joint latent class models used simple Expectation–Maximization (EM) algorithm or direct optimization of the joint likelihood function. However, these methods cannot be readily applied to the case of censored data. Recently, we developed a Monte Carlo EM (MCEM) method based on Metropolis–Hastings algorithm for joint latent class analysis of a censored biomarker and a binary outcome.¹⁵ We now extend this approach to multiple censored biomarkers and a time to event outcome.

The remainder of the paper is organized as follows. In Section 2, we lay out the framework of the joint latent class model and present the estimation procedure via the MCEM algorithm. A simulation study is presented in Section 3. In Section 4, we apply our proposed method to the GenIMS study and demonstrate the utility of joint latent class models for subgroup identification. We conclude the paper with discussions in Section 5.

2 Joint latent class modeling of longitudinal and survival outcomes

We consider a joint latent class model, with three submodels corresponding to latent class membership, survival outcome, and longitudinal biomarkers, respectively. Given the latent class, longitudinal biomarkers and survival outcome are assumed to be independent. In the following, we first specify the models for each component, and then provide the joint likelihood function to account for the censored data in biomarker measurements. Next, we present a MCEM algorithm to estimate the joint latent class model. In the end, we describe the posterior classification method, prediction of survival probability using biomarkers, and a score test to check the conditional independence assumption.

2.1 Model description

2.1.1 Class membership model

Suppose a study cohort of n subjects consists of K heterogeneous latent classes (subpopulations). Let $C_{i} = (C_{i 1}, \dots, C_{iK}) T$ be a vector indicating class membership for subject i, where C_ik = 1 if the ith subject belongs to the kth class $(k = 1, 2, \dots, K)$ and $\sum_{k = 1}^{K} C_{ik} = 1$ We can model the probability of belonging to the kth latent class with a multinomial logistic regression model given by

\begin{matrix} P (C_{ik} = 1) = exp (X_{1 i}^{T} γ_{k}) / [{\sum_{k = 1}}^{K} exp (X_{1 i}^{T} γ_{k})] \end{matrix}

(1)

where γ_k is the coefficient vector associated with a vector of time-independent covariates,

X_{1 i}

. We set

γ_{K} = 0

to satisfy the identifiability criterion.

2.1.2 Class-specific survival model

Let T_i be the true survival time and U_i be the potential right censoring time for the ith subject. Define $V_{i} = min (T_{i}, U_{i})$ as the observed survival time and $Δ_{i} = I (T_{i} < U_{i})$ as the censoring indicator. Given latent class membership, the survival time T_i can be described by a Cox model with hazard function

λ_{i} (t | C_{ik} = 1) = λ_{0 k} (t) exp [X_{2 i}^{T} (t) β_{k}]

(2)

where

λ_{0 k} (t)

is a class-specific baseline hazard function and

β_{k}

is a vector of class-specific regression coefficients for a set of covariates

X_{2 i} (t)

. To reduce the complexity of estimating

λ_{0 k} (t)

, we assume a parametric model for

λ_{0 k} (t)

such as Weibull hazard

λ_{0 k} (t) = λ_{k} ρ_{k} t^{(ρ_{k} - 1)},

where λ_k and ρ_k are the class-specific scale and shape parameters of a Weibull distribution.

2.1.3 Class-specific multivariate linear mixed model

Suppose we have R biomarkers in total and let $W_{ri} = {W_{ri} (t_{i 1}), \dots, W_{ri} (t_{{im}_{i}})}^{T}$ be a vector of measurements for the rth biomarker $(r = 1, \dots, R)$ , recorded at time $t_{ij} (j = 1, \dots, m_{i})$ . The biomarker trajectory, conditional on the latent class, is described by a standard linear mixed model

\begin{matrix} W_{ri} (t_{ij} | b_{ri}, C_{ik} = 1) = X_{3 ri}^{T} (t_{ij}) τ_{rk} + Z_{ri} (t_{ij})^{T} b_{ri} + ɛ_{ri} (t_{ij}) \end{matrix}

(3)

The fixed effect τ_rk is a common or class-specific vector of coefficients for biomarker r, associated with a covariate vector $X_{3 ri} (t)$ , which includes baseline and possibly time-dependent covariates. The vector $Z_{ri} (t_{ij})$ , including polynomial terms of time, is associated with a $p_{r} \times 1$ vector of random effects b_ri, with $b_{ri} = b_{ri} | C_{ik} = 1 \sim N_{p_{r}} (η_{rk}, Σ_{r})$ . The measurement errors, $ɛ_{ri} (t_{ij}) \sim iid N (0, δ_{r}^{2})$ and are independent of random effects.

The interrelationship among multiple biomarkers is characterized by the correlations of marker-specific random effects. Let $b_{i} = (b_{1 i}^{T}, \dots, b_{Ri}^{T}) T$ be a p-vector consisting of random effects for all biomarkers. We assume that b_i follows a multivariate normal distribution $N_{p} (B_{k}, G)$ , with mean vector $B_{k} = (η_{1 k}^{T}, \dots, η_{Rk}^{T}) T$ and a variance-covariance matrix G. The diagonal element of G, Σ _r , is the covariance matrix of random effects for the rth biomarker. The off-diagonal element is the between-marker covariance matrix $cov (b_{ri}, b_{r' i})$ $(r \neq r')$ , which reflects the correlation between two biomarkers.

When there are detection limits due to the low sensitivity of the bioassay used, biomarker level $W_{ri} (t_{ij})$ is not observable if it is out of a detectable range. Suppose the rth biomarker is subject to a lower detection limit, d_ijr, for the jth observation of ith subject, we define the censoring indicator, q_ijr = 1 when the biomarker level is below d_ijr and $q_{ijr} = 0$ when the biomarker level is above d_ijr.

2.2 Joint likelihood function

Let $O_{i} = {v_{i}, Δ_{i}, w_{ri} (t_{ij})^{I (q_{ijr} = 0)}, d_{ijr}^{I (q_{ijr} = 1)}, r = 1, \dots, R, j = 1, \dots, m_{i}}$ be the observed data of survival outcome and biomarker profiles for the ith subject. In order to handle the left-censored biomarker measurements, we incorporate the contribution of censored measurements to the joint likelihood function by using a cumulative distribution function. Given latent class C_i and random effects b_i, multiple biomarkers are independent of each other and they are also independent of the survival outcome. Thus, the joint likelihood function of the observed data is given by

\begin{matrix} Π_{i = 1}^{n} \sum_{k = 1}^{K} P (C_{ik} = 1) f (v_{i}, Δ_{i} | C_{ik} = 1) \int Π_{r = 1}^{R} Π_{j = 1}^{m_{i}} [f (w_{ri} (t_{ij} | b_{i}, C_{ik} = 1)) I (q_{ijr} = 0) \\ \times F (d_{ijr} | b_{i}, C_{ik} = 1) I (q_{ijr} = 1)] f (b_{i} | C_{ik} = 1) d b_{i} \end{matrix}

(4)

where

P (C_{ik} = 1)

is given in equation (1)

\begin{matrix} f (v_{i}, Δ_{i} | C_{ik} = 1) = [λ_{0 k} (v_{i}) exp {X_{2 i}^{T} (v_{i}) β_{k}}]^{Δ_{i}} exp [- \int_{0}^{v_{i}} λ_{0 k} (u) exp {X_{2 i}^{T} (u) β_{k}} d u], \\ f (b_{i} | C_{ik} = 1) = (2 π) - p / 2 | G | - 1 / 2 exp [- \frac{(b_{i} - B_{k}) T G^{- 1} (b_{i} - B_{k})}{2}] \end{matrix}

(5)

f (w_{ri} (t_{ij} | b_{i}, C_{ik} = 1)) = \frac{1}{\sqrt{2 π δ_{r}^{2}}} exp {- \frac{[w_{ri} (t_{ij}) - Z_{ri} (t_{ij})^{T} b_{ri} - X_{3 ri}^{T} (t_{ij}) τ_{rk}]^{2}}{2 δ_{r}^{2}}}

(6)

which is a normal density function with mean

Z_{ri} (t_{ij}) T b_{ri} + X_{3 ri}^{T} (t_{ij}) τ_{rk}

and variance

δ_{r}^{2}

, and

F (d_{ijr} | b_{i}, C_{ik} = 1)

is the corresponding cumulative distribution function evaluated at d_ijr.

2.3 Estimation procedure

When there are no censored data in the biomarker measurements, the estimates of parameters in joint latent class models can be obtained using a standard EM procedure as shown in Lin et al.¹⁶ However, the use of cumulative distribution function for the censored observations makes the calculation of conditional expectations in the E-step intractable. We recently developed an MCEM algorithm for joint latent class modeling of a binary outcome and one longitudinal biomarker with censored observations.¹⁵ This approach can be extended to the case of multiple biomarkers and survival outcome as illustrated below.

Let θ represent all the parameters in the joint latent class model. The complete-data log-likelihood of θ based on the observed data O_i, unobserved random effects b_i and latent class vector C_i is defined as

\begin{matrix} {logL}_{c} (θ) = {\sum_{i = 1}}^{n} log L_{ic} (θ) = {\sum_{i = 1}}^{n} {log f (v_{i}, Δ_{i} | C_{i}) + {\sum_{r = 1}}^{R} {\sum_{j = 1}}^{m_{i}} [I (q_{ijr} = 0) log f (w_{ri} (t_{ij} | b_{i}, C_{i})) \\ + I (q_{ijr} = 1) log F (d_{ijr} | b_{i}, C_{i})] + log f (b_{i} | C_{i}) + log f (C_{i})} \end{matrix}

where

\begin{matrix} log f (v_{i}, Δ_{i} | C_{i}) = Δ_{i} [log {\sum_{k = 1}}^{K} C_{ik} λ_{0 k} (v_{i}) + {\sum_{k = 1}}^{K} C_{ik} X_{2 i}^{T} (v_{i}) β_{k}] \\ - \int_{0}^{v_{i}} [{\sum_{k = 1}}^{K} C_{ik} λ_{0 k} (u)] exp ({\sum_{k = 1}}^{K} C_{ik} X_{2 i}^{T} (u) β_{k}) d u \end{matrix}

(7)

\begin{matrix} log f (w_{ri} (t_{ij} | b_{i}, C_{i})) = - (1 / 2) {log (2 π δ_{r}^{2}) + {w_{ri} (t_{ij}) - {\sum_{k = 1}}^{K} C_{ik} [Z_{ri} (t_{ij})^{T} b_{ri} + X_{3 ri}^{T} (t_{ij}) τ_{rk}]}^{2} / δ_{r}^{2}} \end{matrix}

(8)

log f (b_{i} | C_{i}) = - \frac{p}{2} log (2 π) - \frac{1}{2} log | G | - \frac{(b_{i} - \sum_{k = 1}^{K} C_{ik} B_{k})^{T} G^{- 1} (b_{i} - \sum_{k = 1}^{K} C_{ik} B_{k})}{2}

(9)

\begin{matrix} log f (C_{i}) = \sum_{k = 1}^{K} C_{ik} X_{1 i}^{T} γ_{k} - log [{\sum_{k = 1}}^{K} exp (X_{1 i}^{T} γ_{k})] \end{matrix}

(10)

In order to calculate the conditional expectations in E-step given the observed data O_i, we generate posterior samples of b_i and C_i using Metropolis–Hastings algorithm, as described in our previous work.¹⁵

In the M-step of the $(t + 1)$ th iteration, we update the distribution parameters for random effects as follows

\begin{matrix} {\overset{\land}{B}}_{k}^{(t + 1)} = {\sum_{i = 1}}^{n} E (C_{ik} b_{i} | θ^{(t)}, O_{i}) / n \end{matrix}

(11)

\begin{matrix} {\overset{\land}{G}}^{(t + 1)} = {\sum_{i = 1}}^{n} [E (b_{i} - B_{k}) (b_{i} - B_{k})^{T} | θ^{(t)}, O_{i})] / n \end{matrix}

(12)

When there are no censored data in the biomarker measurements, the updates of ${\overset{\land}{τ}}_{rk}$ and ${\overset{\land}{δ}}_{r}^{2}$ have a closed form as

\begin{matrix} {\overset{\land}{τ}}_{rk}^{(t + 1)} = [{\sum_{i = 1}}^{n} E (C_{ik} X_{3 ri}^{T} X_{3 ri} | θ^{(t)}, O_{i})]^{- 1} {\sum_{i = 1}}^{n} E (C_{ik}) X_{3 ri}^{T} E [w_{ri} (t_{ij}) - Z_{ri} (t_{ij})^{T} b_{ri} - {\sum_{s \neq k}}^{K} C_{is} X_{3 ri}^{T} (t_{ij}) τ_{rs} | θ^{(t)}, O_{i}] \end{matrix}

(13)

\begin{matrix} {\overset{\land}{δ_{r}^{2}}}^{(t + 1)} = {\sum_{i = 1}}^{n} {\sum_{j = 1}}^{m_{i}} E [w_{ri} (t_{ij}) - Z_{ri} (t_{ij})^{T} b_{ri} - {\sum_{k = 1}}^{K} C_{ik} X_{3 ri}^{T} (t_{ij}) τ_{rk} | θ^{(t)}, O_{i}]^{2} / {\sum_{i = 1}}^{n} m_{i} \end{matrix}

(14)

However, there are no closed form for ${\overset{\land}{τ}}_{rk}$ and ${\overset{\land}{δ}}_{r}^{2}$ in the presence of censored data in biomarkers because the cumulative distribution function of biomarker is involved in the likelihood function, so ${\overset{\land}{τ}}_{rk}^{(t + 1)}$ and ${\overset{\land}{δ_{r}^{2}}}^{(t + 1)}$ are estimated by Newton-Raphson iteration method. As the models for latent class membership and survival outcome are nonlinear, we also used Newton-Raphson method to update the parameters in these two models.

The choices of initial values are critical for joint latent class models. We first replace the censored biomarker measurements with half of the detection limits, and then obtain the initial values for longitudinal model and class membership model by using the R package (LCMM) developed by Proust-Lima et al.¹⁷ The LCMM function fits a latent class linear mixed model using only biomarker information. The initial values of parameters in the survival model can be chosen by fitting a parametric survival model. As the mixture distribution under joint latent class models can have multiple local maxima, we fit our model with a set of values around the initials in order to find the global maxima. The MCEM algorithms are implemented with a given number of latent classes. In practice, Bayesian information criterion¹⁸ (BIC) can be used to select the optimal number of classes. Additionally, clinical insights should be taken into account to ensure the latent classes are clinically meaningful.

As shown in Lin et al.,¹⁶ the standard errors of parameters estimates can be obtained through the inverse of the observed information matrix. Under the EM algorithm, the observed information matrix can be approximated by the empirical information matrix

\begin{matrix} I_{e} (θ) = {\sum_{i = 1}}^{n} s_{i} (θ) s_{i} (θ) T - \frac{1}{n} [{\sum_{i = 1}}^{n} s_{i} (θ)] [{\sum_{i = 1}}^{n} s_{i} (θ) T] \end{matrix}

(15)

where

s_{i} (θ)

is the observed score function corresponding to the complete-data log-likelihood, i.e.,

s_{i} (θ) = E (\frac{\partial {logL}_{ic} (θ)}{\partial θ} | O_{i})

. Because

s_{i} (θ) = 0

θ = \overset{\land}{θ},

the second term of

I_{e} (\overset{\land}{θ})

is 0. Therefore, the variance matrix of

\overset{\land}{θ}, Var (\overset{\land}{θ})

, can be estimated by

I_{e} (\overset{\land}{θ}) - 1

with

I_{e} (\overset{\land}{θ}) = (\sum_{i = 1}^{n} s_{i} (θ) s_{i} (θ) T) | θ = \overset{\land}{θ}

2.4 Posterior classification

Given the parameter estimates $\overset{\land}{θ}$ and observed data O_i, the posterior probability of belonging to class k for subject i can be computed using Bayes rule as follows:

\begin{matrix} {\overset{\land}{π}}_{i, k} = p (C_{ik} = 1 | O_{i}; \overset{\land}{θ}) = \frac{p (C_{ik} = 1 | \overset{\land}{θ}) f (v_{i}, Δ_{i} | C_{ik} = 1; \overset{\land}{θ}) f (w_{i} | C_{ik} = 1; \overset{\land}{θ})}{\sum_{k = 1}^{K} p (C_{ik} = 1 | \overset{\land}{θ}) f (v_{i}, Δ_{i} | C_{ik} = 1; \overset{\land}{θ}) f (w_{i} | C_{ik} = 1; \overset{\land}{θ})} \end{matrix}

(16)

where

p (C_{ik} = 1 | \overset{\land}{θ})

is the prior probability defined in the equation (1) and

\begin{matrix} f (w_{i} | C_{ik} = 1; \overset{\land}{θ}) = \int Π_{r = 1}^{R} Π_{j = 1}^{m_{i}} [f (w_{ri} (t_{ij} | b_{i}, C_{ik} = 1)) I (q_{ijr} = 0) F (d_{ijr} | b_{i}, C_{ik} = 1) I (q_{ijr} = 1)] \\ \times f (b_{i} | C_{ik} = 1) P (C_{ik} = 1) d b_{i} \end{matrix}

We can classify subjects based on the derived posterior probabilities, with each subject being classified to the class which has the largest posterior probability.

2.5 Prediction of event probability using biomarkers

Let $w_{i}^{(s)}$ be the biomarker measurements up to time s. We can predict the probability of having an event by time t (t > s) given $w_{i}^{(s)}$ as follows

\begin{matrix} P (T_{i} \leq t | w_{i}^{(s)}, \overset{\land}{θ}) = {\sum_{k = 1}}^{K} P (T_{i} \leq t | C_{ik} = 1, \overset{\land}{θ}) P (C_{ik} = 1 | w_{i}^{(s)}; \overset{\land}{θ}) \end{matrix}

(17)

where

\begin{matrix} P (T_{i} \leq t | C_{ik} = 1, \overset{\land}{θ}) = 1 - exp [- \int_{0}^{t} λ_{0 k} (u) exp (X_{2 i}^{T} (u) β_{k}) d u] \end{matrix}

and

\begin{matrix} P (C_{ik} = 1 | w_{i}^{(s)}; \overset{\land}{θ}) = \frac{P (C_{ik} = 1 | \overset{\land}{θ}) f (w_{i}^{(s)} | C_{ik} = 1; \overset{\land}{θ})}{\sum_{k = 1}^{K} P (C_{ik} = 1 | \overset{\land}{θ}) f (w_{i}^{(s)} | C_{ik} = 1; \overset{\land}{θ})} \end{matrix}

which is the posterior probability of subject i belonging to class k given only longitudinal biomarkers up to time s, $w_{i}^{(s)}$ .

2.6 Score test for conditional independence

The key assumption of a joint latent class model is the conditional independence between biomarkers and the survival outcome given latent classes and correct specification of covariates for biomarker and survival models. We evaluate this assumption by extending the score test in Jacqmin-Gadda et al.¹⁹ to censored biomarker data. We test the null hypothesis of conditional independence given latent classes versus the alternative hypothesis that the survival outcome and the biomarkers are conditionally independent given both latent classes and random effects in the biomarker models. Therefore, the joint latent class model described in equations (1) to (3) is under null hypothesis H₀. Under the alternative hypothesis $H_{a}$ , we retain equations (1) and (3), but include shared random effects b_i in the survival model as follows

\begin{matrix} λ_{a} (t | b_{i}, C_{ik} = 1) = λ_{0 k} (t) exp [X_{2 i}^{T} (t) β_{k} + b_{i}^{T} φ] \end{matrix}

(18)

where the

p \times 1

vector of parameters φ captures the additional association between survival outcome and biomarker profiles within the latent classes. Then a hypothesis test for conditional independence is equivalent to test

H_{0} ⋮ φ = 0

versus

H_{a} ⋮ φ \neq 0

As shown in Jacqmin-Gadda et al.,¹⁹ the score function under H₀ can be expressed as

\begin{matrix} U (φ, θ) | φ = 0 = U (0, θ) = \sum_{i = 1}^{n} U_{i} (0, θ) = \sum_{i = 1}^{n} \sum_{k = 1}^{K} π_{i, k} [Δ_{i} - Λ_{k} (v_{i} | C_{ik} = 1)] E (b_{i} | w_{i}, C_{ik} = 1) \end{matrix}

(19)

where

π_{i, k}

is the posterior probability, which can be estimated by equation (16).

Λ_{k} (.)

is the class-specific cumulative hazard function.

E (b_{i} | w_{i}, C_{ik} = 1)

can be estimated using the posterior samples obtained from the MCEM algorithm. The score statistic

U (0, \overset{\land}{θ})

is asymptotically normally distributed with mean 0. Freedman²⁰ proposed a robust variance of the score statistic using the empirical covariance matrix:

\begin{matrix} Var (U) = \sum_{i = 1}^{n} U_{i} (0, \overset{\land}{θ}) U_{i} (0, \overset{\land}{θ})^{T} - U_{i} (0, \overset{\land}{θ}) U_{i} (0, \overset{\land}{θ})^{T} / n \end{matrix}

It has been shown that the robust variance has similar performance as the asymptotic variance.¹⁹ Therefore, we use the robust variance to construct the score test. Under the null hypothesis, $U (0, \overset{\land}{θ}) T Var (U) - 1 U (0, \overset{\land}{θ}) T$ follows a $χ_{p}^{2}$ distribution asymptotically.

3 Simulation study

In the simulation study, we generated a survival outcome and repeated measurements of two longitudinal biomarkers from the following joint latent class model with two latent classes (k = 1, 2).

Latent class model: $P (C_{i 1} = 1 | x_{i}) = exp (γ_{1} + γ_{2} x_{i}) / [1 + exp (γ_{1} + γ_{2} x_{i})]$ , where the covariate x_i is generated from a normal distribution N(0, 1).

Survival model: $λ (t | C_{ik} = 1) = λ_{0 k} (t) exp (β_{k} x_{i})$ , where the baseline hazard function follows a Weibull distribution, $λ_{0 k} (t) = λ_{k} ρ_{k} t^{(ρ_{k} - 1)}$ . The censoring time is generated from a uniform distribution $7 + U (14, 90)$ .

Bivariate mixed model:

Biomarker 1 ⋮ W_{1 i} (t_{ij} | b_{1 i}, C_{ik} = 1) = b_{0 i}^{(1)} + τ_{1 k}^{x} x_{i} + b_{1 i}^{(1)} t_{ij} + ɛ_{1 i} (t_{ij})

Biomarker 2 ⋮ W_{2 i} (t_{ij} | b_{2 i}, C_{ik} = 1) = b_{0 i}^{(2)} + τ_{2 k}^{x} x_{i} + b_{1 i}^{(2)} t_{ij} + ɛ_{2 i} (t_{ij})

where class-specific random intercept

b_{0 i}

, and slope

b_{1 i}

are considered for each biomarker. We assume that the joint distribution of random intercepts and slopes follow a multivariate normal distribution with class-specific mean and a common variance-covariance matrix G, i.e.,

(b_{0 i}^{(1)}, b_{1 i}^{(1)}, b_{0 i}^{(2)}, b_{1 i}^{(2)}) \sim N ((η_{01 k}, η_{11 k}, η_{02 k}, η_{12 k}), G)

, with

\begin{matrix} G = (\begin{matrix} 1 & 0.05 & 0.1 & 0.1 \\ 0.05 & 1 & 0.1 & 0.1 \\ 0.1 & 0.1 & 1 & 0.05 \\ 0.1 & 0.1 & 0.05 & 1 \end{matrix}) \end{matrix}

The measurement errors $ɛ_{1 i} (t_{ij}) \sim i . i . d N (0, δ_{1}^{2}), ɛ_{2 i} (t_{ij}) \sim i . i . d N (0, δ_{2}^{2})$ .

To mimic the biomarker sampling schedule in the GenIMS study, each subject has at most seven repeated measurements on the biomarkers, observed at days

t_{ij} = 1, 2, \dots, 7

before the occurrence of the event of interest or censoring. The coefficients,

τ_{1 k}^{x}

and

τ_{2 k}^{x}

, denote class-specific covariate effects of x_i on two biomarkers. The true parameter values of the above joint latent class model are given in Table 1 along with the simulation results. The overall censoring rate of survival outcome is about 27.3% with the median event time at 17.5. The censoring percentages for two longitudinal biomarkers are set at 20% and 25%, respectively. The corresponding lower detection limits are selected empirically for each simulated dataset.

Table 1.

Simulation results for a joint latent class model with two biomarkers and two latent classes.

		No censoring				20% + 25% censoring
Model	Parameter	Bias	SD	SE	CP	Bias	SD	SE	CP
Latent class	$γ_{1} = 0.4$	−0.0087	0.0880	0.0887	0.95	−0.0188	0.0908	0.0877	0.96
	$γ_{2} = 0.2$	0.0091	0.0723	0.0744	0.94	0.0254	0.0803	0.0812	0.94
Biomarker 1
Class 1	$η_{011} = 5$	0.0003	0.0453	0.0471	0.97	0.0017	0.0481	0.0451	0.97
	$η_{111} = - 0.05$	0.0009	0.0115	0.0121	0.96	0.0016	0.0137	0.0129	0.95
	$τ_{11}^{x} = 0.8$	0.0015	0.0218	0.0255	0.94	0.0085	0.0318	0.0327	0.96
Class 2	$η_{012} = 3$	0.0070	0.0371	0.0388	0.94	0.0090	0.0452	0.0461	0.95
	$η_{112} = - 0.1$	−0.0005	0.0135	0.0143	0.96	−0.0058	0.0136	0.0137	0.94
	$τ_{12}^{x} = 0.7$	0.0082	0.0215	0.0221	0.96	0.0093	0.0325	0.0361	0.97
Variance	$δ_{1}^{2} = 0.1$	−0.0022	0.0013	0.0021	0.94	−0.0042	0.0032	0.0040	0.94
Biomarker 2
Class 1	$η_{021} = 6$	−0.0020	0.0481	0.0467	0.96	−0.0081	0.0518	0.0512	0.95
	$η_{121} = - 0.1$	0.0004	0.0127	0.0134	0.97	0.0077	0.0143	0.0141	0.96
	$τ_{21}^{x} = 0.5$	−0.0063	0.0134	0.0121	0.94	−0.0069	0.0236	0.0241	0.95
Class 2	$η_{022} = 4$	0.0014	0.0393	0.0427	0.94	0.0094	0.0445	0.0425	0.95
	$η_{122} = - 0.15$	0.0015	0.0139	0.0137	0.96	0.0085	0.0149	0.0151	0.94
	$τ_{22}^{x} = 0.6$	0.0063	0.0231	0.0199	0.95	0.0096	0.0312	0.0299	0.96
Variance	$δ_{2}^{2} = 0.2$	−0.0011	0.0022	0.0035	0.95	−0.0037	0.0041	0.0053	0.95
Survival response
Class 1	$λ_{1} = 0.05$	−0.0008	0.0072	0.0070	0.95	−0.0032	0.0080	0.0077	0.94
	$ρ_{1} = 0.7$	0.0033	0.0165	0.0159	0.95	0.0033	0.0355	0.0316	0.96
	$β_{1} = - 1$	−0.0028	0.0344	0.0312	0.94	−0.0068	0.0604	0.0599	0.94
Class 2	$λ_{2} = 0.1$	−0.0018	0.0113	0.0132	0.94	−0.0051	0.0134	0.0147	0.95
	$ρ_{2} = 0.5$	−0.0067	0.0274	0.0288	0.96	−0.0093	0.0416	0.0431	0.96
	$β_{2} = - 0.7$	0.0054	0.0521	0.0529	0.95	0.0083	0.0861	0.0872	0.96

SD: standard deviation; SE: standard error; CP: coverage probability.

A total of 100 datasets with sample size n = 1000 are simulated. The model parameters are estimated by the proposed MCEM Algorithm. The estimating procedure was implemented in R and the source codes for simulation studies are provided in the Supplementary Material. The average computation time for a single simulated dataset was approximately 1 h when the R program was run on a computer with an Intel Core i7 7200 U @ 2.5 GHz processor and 8 GB of RAM. We summarize the bias, sample standard deviation (SD), mean standard error (SE), and empirical 95 % coverage probability (CP) for each parameter in Table 1. For comparison, we also provided the estimation results when the biomarkers are not censored at the detection limits. The estimates for all the parameters are approximately unbiased and the empirical coverage probabilities are close to 0.95 for both scenarios, with and without censoring. The means of the estimated standard errors are in good agreement with the empirical standard deviations of the estimates. As expected, the standard errors of regression parameter estimates are larger in general when the biomarkers are censored due to detection limits.

4 Application to GenIMS study

4.1 Data overview

The GenIMS study is a multicenter cohort study that enrolled 2320 subjects from the emergency department with CAP in 28 US academic and community hospitals between 2001 and 2003.⁴ Multiple biomarkers on different biological pathways were measured on the patients admitted to the hospital daily in the first week, and weekly thereafter while they were still in the hospital. To illustrate our proposed methods, we focus on two cytokines, IL-6, a pro-inflammatory biomarker, and IL-10, an anti-inflammatory biomarker. Both IL-6 and IL-10 are subject to lower detection limits, with thresholds 2 or 5 pg/mL for IL-6 (depending on the bioassay used) and 5 pg/mL for IL-10, respectively.

The objective of our analysis is to determine the patterns of IL-6 and IL-10 profiles in the first week of hospitalization, and identify the subgroups of patients with different patterns of cytokine responses and associated risks of 90-day mortality. In the analysis cohort of 1882 subjects who had confirmed CAP and had at least one IL-6 measurement and one IL-10 measurement, the mean (SD) age was 67.3 (16.75), 52.1% were male, 80.8% were white, and 72.5% had Charlson Comorbidity Index (CCI) greater than zero. The overall 90-day mortality rate is 11.4%. The censoring percentage in biomarkers increased over time, ranging from 13.5% to 35.5% for IL-6, and from 46.9% to 78.9% for IL-10.

4.2 Model fitting

We fit a joint latent class model of survival time and two biomarkers: IL-6 and IL-10, with the number of latent classes varying from 1 to 3. The biomarker measurements were transformed to natural log scales in the analysis. Age, gender (1 for male, and 0 for female), and CCI (>0 vs. 0) are included as covariates in each component of the joint latent class model, and their effects on biomarkers’ trajectories and survival time are assumed to be common across the latent classes. Random intercept and slope are included in the mixed model for each biomarker. An unstructured covariance matrix is specified for the random effects to reflect the correlation between IL-6 and IL-10. The maximized log-likelihood values and BIC values are calculated for each of the three joint models and are listed in Table 2. The computation times are approximately 3, 10, and 28 h under K = 1, 2, and 3, respectively. A two-class joint model is optimal with the lowest BIC value, and the corresponding fitted model is summarized as follows. The standard errors and p-values based on Wald tests are also provided in Table 3.

Class membership model (k = 1, 2)

logit [P (C_{i 2} = 1)] = 2.136 - 0.004 \times Age - 0.149 \times Gender + 0.106 \times I (CCI > 0)

Survival model for 90-day mortality

λ (t | C_{ik} = 1) = λ_{0 k} (t) exp [0.033 \times Age + 0.305 \times Gender + 0.669 \times I (CCI > 0)]

where Weibull hazard function

\begin{matrix} λ_{0 k} (t) = λ_{k} ρ_{k} t^{(ρ_{k} - 1)} = {\begin{matrix} 0.0010 \times 0.552 \times t^{(0.552 - 1)}, & class 1 \\ 0.0001 \times 0.796 \times t^{(0.796 - 1)}, class 2 \end{matrix} \end{matrix}

Bivariate mixed model for IL-6 and IL-10 in log scale

Table 2.

GenIMS study: BIC for different numbers of latent classes.

K	Negative Log-Likelihood	Bayesian information criterion
1	78,501	156,102
2	50,156	100,622
3	65,450	133,278

Table 3.

GenIMS study: Parameter estimates for the joint latent class model with two latent classes.

Model	Parameter	Estimate	SE	p-value
Latent class (probability of belonging to class 2)	Intercept	2.136	0.191	<0.001
	Age	−0.004	0.013	0.380
	Male	−0.149	0.072	0.005
	CCI > 0	0.106	0.053	0.090
Longitudinal IL-6
Class 1	Intercept	4.801	0.158	<0.001
	Linear	−0.754	0.026	<0.001
	Quadratic	0.052	0.007	<0.001
Class 2	Intercept	3.010	0.151	<0.001
	Linear	−0.574	0.015	<0.001
	Quadratic	0.043	0.006	<0.001
Covariates
	Age	0.006	0.001	<0.001
	Male	0.451	0.152	0.001
	CCI > 0	0.328	0.136	0.021
Longitudinal IL-10
Class 1	Intercept	2.448	0.153	<0.001
	Linear	−0.039	0.014	0.007
	Quadratic	−0.021	0.003	<0.001
Class 2	Intercept	1.122	0.173	<0.001
	Linear	−0.302	0.043	<0.001
	Quadratic	0.029	0.009	0.020
Covariates
	Age	0.007	0.002	<0.001
	Male	0.572	0.199	0.003
	CCI > 0	0.368	0.163	0.024
Survival
Class 1	Scale	0.001	0.00005	<0.001
Class 1	Shape	0.552	0.034	<0.001
Class 2	Scale	0.0001	0.00006	<0.001
Class 2	Shape	0.796	0.057	<0.001
Covariates
	Age	0.033	0.002	<0.001
	Male	0.305	0.018	<0.001
	CCI > 0	0.669	0.016	<0.001

SE: standard error; IL: interleukin; CCI: Charlson Comorbidity Index.

Mean IL-6:

\begin{matrix} {\begin{matrix} 4.801 - 0.754 \times t_{ij} + 0.052 \times t_{ij}^{2} + 0.006 \times Age + 0.451 \times Gender + 0.328 \times I (CCI > 0), & class 1 \\ 3.010 - 0.574 \times t_{ij} + 0.043 \times t_{ij}^{2} + 0.006 \times Age + 0.451 \times Gender + 0.328 \times I (CCI > 0), & class 2 \end{matrix} \end{matrix}

Mean IL-10:

\begin{matrix} {\begin{matrix} 2.448 - 0.039 \times t_{ij} - 0.021 \times t_{ij}^{2} + 0.007 \times Age + 0.572 \times Gender + 0.368 \times I (CCI > 0), & class 1 \\ 1.122 - 0.302 \times t_{ij} + 0.029 \times t_{ij}^{2} + 0.007 \times Age + 0.572 \times Gender + 0.368 \times I (CCI > 0), & class 2 \end{matrix} \end{matrix}

The variance estimates of measurement errors are 0.014 for IL-6 and 0.013 for IL-10. Denote the random intercept and slope by $b_{0 i}^{(1)}, b_{1 i}^{(1)}$ for IL-6 and $b_{0 i}^{(2)}, b_{1 i}^{(2)}$ for IL-10. The estimate of covariance matrix of random effect vector $b_{i} = (b_{0 i}^{(1)}, b_{1 i}^{(1)}, b_{0 i}^{(2)}, b_{1 i}^{(2)})$ is

\begin{matrix} \overset{\land}{G} = (\begin{matrix} 3.312 & 0.178 & - 0.234 & - 0.052 \\ 0.178 & 0.525 & 0.000 & - 0.004 \\ - 0.234 & 0.000 & 0.023 & 0.003 \\ - 0.052 & - 0.004 & 0.003 & 0.018 \end{matrix}) \end{matrix}

The between-subject variations are larger in IL-6 than in IL-10, especially in the baseline level of IL-6. The correlation between random slope and random intercept are 0.14 and 0.15 for IL-6 and IL-10, respectively. The correlation between random intercepts of two biomarkers is −0.84, indicating a strong correlation between two biomarkers.

4.3 Latent classes characteristics

Figures 1 shows the class-specific estimated biomarker curves and survival functions for a male with age 67 and CCI score >0, respectively. We refer to class 1 with high IL-6 and high IL-10 values as “high cytokine” class and class 2 with low IL-6 and low IL-10 values as “low cytokine” class. The cytokine responses to pneumonia appear to be balanced between pro-inflammatory IL-6 and anti-inflammatory IL-10. After posterior classification, we found that about 89% of patients belong to “low cytokine” class. For patients with the same age, gender, and CCI status, those in “high cytokine” class have a consistently lower survival probability. On average, IL-6 decreases more rapidly in the first half of the week in “high cytokine” class compared to “low cytokine” class, while IL-10 declines relatively faster at the beginning in “low cytokine” class, but shows faster decline post day 4 in “high cytokine” class.

For subjects in the same latent class, being older, male, and having CCI score >0 are significantly associated with higher cytokine trajectories and increased risk of 90-day mortality. The results from latent class membership model indicate that gender is significantly associated with latent class membership. There are more males in “high cytokine” and low survival rate class.

4.4 Goodness of fit

We evaluate the goodness-of-fit of our joint model with respect to posterior classification, comparison of fitted values versus observed values, and the conditional independence assumption. In posterior classification, subjects were classified to the classes which have the highest posterior probabilities. Higher posterior probabilities close to 1 indicate better discriminatory power of the model. The mean maximal posterior probabilities of our model are 0.96 and 0.99 for class 1 and 2, respectively, indicating unambiguous classification.

Next, we compare the predicted values with observed values for both biomarker measurements and survival outcome. The predicted values for rth biomarker at time t can be calculated using class-specific predictions averaged over all subjects who had observations at time t, i.e., ${\overset{\land}{W}}_{r}^{k} (t) = \sum_{i = 1}^{N (t)} {\overset{\land}{π}}_{i, k} {\overset{\land}{W}}_{ri} (t) | C_{ik} = 1,$ where N(t) is the number of subjects who had observations at time t and ${\overset{\land}{π}}_{i, k}$ is the posterior probability from equation (16). Then the predicted values ${\overset{\land}{W}}_{r}^{k} (t)$ are compared with the class-specific mean of observed values ${\bar{W}}_{r}^{k} (t) = \sum_{i = 1}^{N (t)} {\overset{\land}{π}}_{i, k} W_{ri} (t)$ , with censored values of biomarker measurements replaced by half of detection limits. As shown in Figure 2, the fitted curves track the observed curves very well for both biomarkers. In a similar way, we compare the predicted survival probability with the Kaplan–Meier estimate. The predicted class-specific survival function is ${\overset{\land}{S}}_{k} (t) = \sum_{i = 1}^{n} {\overset{\land}{π}}_{i, k} {\overset{\land}{S}}_{ik} (t)$ , where ${\overset{\land}{S}}_{ik} (t)$ can be derived from survival model (2). Then we compare ${\overset{\land}{S}}_{k} (t)$ with class-specific weighted Kaplan–Meier estimates. We can see that for both classes, the fitted curves track the Kaplan–Meier curves very well (Figure 3).

The key assumption for our joint latent class model is the conditional independence between biomarkers and the survival outcome given latent classes. We apply the score test described in Section 2.6 to evaluate the conditional independence assumption. We consider an alternative hypothesis that the survival outcome depends on two biomarkers through shared random intercepts and slopes in addition to the latent classes. Therefore, the score test statistic asymptotically follows a $χ^{2}$ distribution with 4 degrees of freedom under null hypothesis. Our two-class joint model has a score test statistic equal to 5.80, which is smaller than the critical value (i.e., 9.48) of $χ_{4}^{2}$ at 0.05 significance level, thus the null hypothesis of conditional independence between biomarkers and survival outcome given latent classes is not rejected.

4.5 Comparison with one biomarker model on posterior classification and prediction

For comparison, we fit joint latent class models with only IL-6 or IL-10 biomarker. The optimal number of classes is also two when fitting with only IL-6 or IL-10 biomarker based on BIC scores. In order to compare the discrimination ability between the two-marker model and one-marker models, we calculate the mean of maximal posterior probabilities for subjects classified in each latent class. One-marker model based on IL-6 yields mean values 0.86 for class 1, and 0.95 for class 2. The model based on IL-10 results in 0.94 and 0.98 for classes 1 and 2. Using a two-marker model, the values increase to 0.95 and 0.99 for two classes. Since higher posterior probabilities stand for more unambiguous classification, the model using both IL-6 and IL-10 biomarkers shows clear improvement over the single-marker models.

We also compare the prediction accuracy between the two-marker model and one-marker models. Mortality rates at 30 days and 90 days are estimated using equation (17) and the corresponding receiver operating characteristic (ROC) curves are provided in Figure 4. The two-marker model shows better discrimination power than single-marker models for both 30-day and 90-day mortality. The areas under the ROC curves (AUCs) for 30-day mortality are 0.66, 0.72, and 0.75, respectively, for the models based on IL-6 only, IL-10 only, and two biomarkers together. The corresponding AUCs for 90-day mortality are 0.70, 0.73, and 0.75. Incorporating the information of two biomarkers into the model improves prediction accuracy in terms of both specificity and sensitivity.

5 Discussion

Biomarker data are often collected over time to understand the development and progression of diseases. However, the presence of limits of detection in bioassays can hamper the statistical evaluation of biomarkers in biomedical and clinical studies. In this paper, we focused on addressing left-censoring issue in the biomarker measurements for joint latent class models. Our research work facilitates the utilization of this joint modeling approach to explore the heterogeneity of biomarker profiles and assist in risk stratification. We developed an estimation procedure based on likelihood methods to incorporate the censored biomarker measurements and correlations between multiple biomarkers. Simulation studies showed that the proposed MCEM algorithm for parameter estimation had satisfactory performance. As illustrated in the application to GenIMS study, cytokine biomarker responses to pneumonia were not homogeneous, indicating that different treatment strategies may be needed to effectively manipulate cytokine activation. The advantage of our model compared to the previous analysis of GenIMS cytokine data⁴ is that we identified subgroups based on longitudinal biomarkers and survival outcome simultaneously instead of using a two-stage analysis. In addition, our model demonstrated better explanatory ability by considering subgroups associated with different risks. It is worth noting that the latent class modeling is a data-driven approach. The identified subgroups should be interpreted in an exploratory fashion and their utility for prediction and risk stratification should be validated in new data. There are several limitations with our proposed method. First, the MCEM algorithm has high computational cost. The required number of posterior samples increases exponentially as the number of biomarkers increases. In order to reduce the computational burden, Bernhardt et al.¹⁰ proposed a fast, approximate EM algorithm to reduce the dimension in the E-step of the algorithm. Also, Chen et al.²¹ developed a fully Bayesian method for joint modeling of longitudinal and survival data with censored covariates, where censored values were treated as extra parameters in the Markov chain Monte Carlo algorithm. These methods were both developed under the framework of shared random-effects models. Adoption of these methods to joint latent class models in the presence of censored biomarkers requires further consideration. Second, we assumed conditional independence between biomarkers and clinical outcomes given latent classes. In some scenarios, this assumption may not be adequate to capture the correlation between the two processes. A more general setting to let the biomarkers and clinical outcomes depend on both latent classes and shared random effects may be considered.²² Finally, the parametric distribution assumption for the survival outcome can be restrictive. For semiparametric survival models such as Cox models, a Breslow-type estimator or piecewise constant estimator can be considered for baseline hazard functions. Recently, joint latent class models have been increasingly used as a prognostic tool for dynamic prediction of a clinical outcome based on the biomarker profiles. Several prediction accuracy measures have been developed for joint latent class models.²³ Extension of these methods to censored biomarkers warrants further study.

Supplemental Material

Supplemental material for A latent class approach for joint modeling of a time-to-event outcome and multiple longitudinal biomarkers subject to limits of detection

Supplemental Material for A latent class approach for joint modeling of a time-to-event outcome and multiple longitudinal biomarkers subject to limits of detection by Menghan Li, Ching-Wen Lee and Lan Kong in Statistical Methods in Medical Research

Footnotes

Acknowledgements

We thank Dr. Derek Angus and the Clinical Research, Investigation, and Systems Modeling of Acute Illness (CRISMA) Center at the Department of Critical Care Medicine, University of Pittsburgh, for access to the GenIMS data. GenIMS was funded by NIGMS R01 GM61992 with additional support from GlaxoSmithKline for enrolment and clinical data collection, and Diagnostic Products Corporation for the cytokine assays. The authors thank all the investigators who contributed to the generation of the data, a full list of whom was published in Critical Care Medicine 2007; 35: 1061–1067.

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) received no financial support for the research, authorship, and/or publication of this article.

Supplemental material

Supplemental material for this article is available online.

References

Mellors

Munoz

Giorgi

, et al. Plasma viral load and CD4+ lymphocytes as prognostic markers of HIV-1 infection. Ann Intern Med 1997; 126: 946–954.

Thiebaut

Jacqmin-Gadda

Babiker

. Joint modeling of bivariate longitudinal data with informative dropout and left-censoring, with application to the evolution of CD4+ cell count and HIV RNA viral load in response to treatment of HIV infection. Stat Med 2005; 24: 65–82.

Proust-Lima

Jacqmin-Gadda

Joly

. Joint modelling of multivariate longitudinal outcomes and a time-to-event: a nonlinear latent class approach. Comput Stat Data Anal 2009; 53: 1142–1154.

Kellum

Kong

Fink

, et al. Understanding the inflammatory cytokine response in pneumonia and sepsis. Arch Intern Med 2007; 167: 1655–1663.

Tsiatis

Davidian

. Joint modeling of longitudinal and time-to-event data: an overview. Stat Sin 2004; 14: 809–834.

Proust-Lima

Sene

Jeremy

, et al. Joint latent class models for longitudinal and time-to-event data: a review. Stat Meth Med Res 2014; 23: 74–90.

Lin

McCulloch

Mayne

. Maximum likelihood estimation in the joint analysis of time-to-event and multiple longitudinal variables. Stat Med 2002; 21: 2369–2382.

Ibrahim

Chen

Sinha

. Bayesian methods for joint modeling of longitudinal and survival data with applications to cancer vaccine trials. Stat Sin 2004; 14: 863–883.

Zeger

. The evaluation of multiple surrogate endpoints. Biometrics 2001; 57: 81–87.

10.

Bernhardt

Zhang

Wang

. A fast EM algorithm for fitting joint models of a binary response and multiple longitudinal covariates subject to detection limits. Comput Stat Data Anal 2015; 114: 961–968.

11.

Hughes

. Mixed effects models with censored data with application to HIV RNA levels. Biometrics 1999; 55: 625–629.

12.

Lyles

Taylor

. Random regression models for human immunodeficiency virus ribonucleic acid data subject to left censoring and informative drop-outs. Appl Stat 2000; 49: 485–498.

13.

Jacqmin-Gadda

Thiebaut

Chene

, et al. Analysis of left-censored longitudinal data with application to viral load in HIV infection. Biostatistics 2000; 1: 355–368.

14.

Wannemuehler

Lyles

Manatunga

, et al. Likelihood-based methods for estimating the association between a health outcome and left-or interval-censored longitudinal exposure data. Stat Med 2010; 29: 1661–1672.

15.

Kong

. Joint analysis of left-censored longitudinal biomarker and binary outcome via latent class modeling. Stat Med. Epub ahead of print 2 April 2018. doi: 10.1002/sim.7642.

16.

Lin

Turnbull

McCulloch

, et al. Latent class models for joint analysis of longitudinal biomarker and event process data: application to longitudinal prostate-specific antigen readings and prostate cancer. J Am Stat Assoc 2002; 97: 53–65.

17.

Proust-Lima

Philipps

Liquet

. Estimation of extended mixed models using latent classes and latent processes: the R package LCMM. J Stat Software 2017; 78: 1–56.

18.

Schwarz

. Estimating the dimension of a model. Ann Statist 1978; 6: 461–464.

19.

Jacqmin-Gadda

Proust-Lima

Taylor

, et al. Score test for conditional independence between longitudinal outcome and time to event given the classes in the joint latent class model. Biometrics 2010; 66: 11–19.

20.

Freedman

. How can the score test be inconsistent? Am Stat 2007; 61: 291–295.

21.

Chen

May

Ibrahim

, et al. Joint modeling of longitudinal and survival data with missing and left-censored time-varying covariates. Stat Med 2014; 33: 4560–4576.

22.

Liu

Zhou

. Joint latent class model of survival and longitudinal data: an application to CPCRA study. Comput Stat Data Anal 2015; 91: 40–50.

23.

Blanche

Proust-Lima

Loubre

, et al. Quantifying and comparing dynamic predictive accuracy of joint models for longitudinal marker and time-to-event in presence of censoring and competing risks. Biometrics 2015; 71: 102–113.

Supplementary Material

Please find the following supplemental material available below.

For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.

For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.

0.01 MB

0.00 MB