A Two-stage Approach for Variable Selection in Joint Modeling of Multiple Longitudinal Markers and Competing Risk Outcomes

A multivariate linear mixed-effects model is considered for the multiple longitudinal markers as follows:

where $\varepsilon_{ijk}\overset{\text{iid}}{\sim}N(0,\sigma_{k}^{2})$. We assume that

where ${X}_{ik}(s)$ and ${Z}_{ik}(s)$ are the design vectors at time $s$ for the $p^{\beta}_{k}$-dimensional fixed effects ${\beta}_{k}$ and $q_{k}$-dimensional random effects ${b}_{ik}$, respectively. We assume that ${b}_{i}=({b}_{i1}^{\top},\ldots,{b}_{iK}^{\top})^{\top}$ takes into account the association within and between the multiple longitudinal markers such that ${b}_{i}\sim N({0},{\Sigma})$, and ${\Sigma}$ is a $q\times q$ covariance matrix, where $q=\sum_{k=1}^{K}q_{k}$. It is assumed that ${\varepsilon}_{ijk}$ is independent of ${b}_{ik}$, $i=1,\ldots,K$.

For each of the $L$ causes, a cause-specific proportional hazard model is postulated. The risk of an event is assumed to be a function of the subject-specific linear predictors $\eta_{ik}(t|{b}_{ik})$ and of time-invariant covariates as follows [40]:

where $\lambda_{0l}(t)$ denotes the baseline hazard function for the $l$th cause, ${\omega}_{il}$ is a $p^{\gamma}_{l}$ vector of covariates with corresponding regression coefficients ${\gamma}_{l}$, and ${\alpha}_{l}=({\alpha}_{l1},\ldots,{\alpha}_{lK})^{\top}$ represents the vector of association parameters between the longitudinal markers and cause $l$. There are different approaches for specifying the baseline hazard function in joint modeling [40]. The methodology is not restricted to a specific form of the baseline hazard, and we consider a piecewise constant baseline hazard [27]. Therefore, a finite partition of the time axis is constructed as $0=\xi_{0}<\xi_{1}<\xi_{2}<\ldots,\xi_{\mathcal{K}}$, such that $max(T_{1},\ldots,T_{N})<\xi_{\mathcal{K}}$. Thus, there are $\mathcal{K}$ intervals such as $(0,\xi_{1}],(\xi_{1},\xi_{2}],\ldots,(\xi_{\mathcal{K}-1},\xi_{\mathcal{K}}]$. We consider a constant baseline hazard for each interval as $\lambda_{0l}(t)=\lambda_{\kappa l},~{}t\in(\xi_{\kappa-1},\xi_{\kappa}],~{}\kappa=1,\ldots,\mathcal{K}$.

Let ${\psi}_{k}$ and ${\Upsilon}$ represent the unknown parameters for the $k$th longitudinal marker and the time to event outcomes, respectively, excluding the covariance matrix ${\Sigma}$. Let ${\psi}=\{{\psi}_{1},\ldots,{\psi}_{K}\}$. If ${y}=\{{y}_{1},\ldots,{y}_{N}\}$, ${t}=(t_{1},\ldots,t_{N})^{\top}$, and ${\delta}=(\delta_{1},\ldots,\delta_{N})^{\top}$, the joint likelihood function for models (1) and (3) is as follows:

where $\Lambda_{il}(t|{\omega}_{i},{{b}}_{i})$ represents the cumulative hazard function of Equation (3), and the densities of the univariate normal distribution with mean $\mu$ and variance $\sigma^{2}$, and the $q$-variate normal distribution with mean ${\mu}$ and covariance matrix ${\Sigma}$, are denoted as $\phi(x;\mu,\sigma^{2})$ and $\phi_{q}({x};{\mu},{\Sigma})$, respectively.
In the frequentist paradigm, parameter estimation can be computed by maximizing (2.3), which involves the use of numerical integration [14, 48]. This integration will be more challenging when the number of outcomes increases. In the Bayesian paradigm [57, 42, 4, 30], the following joint posterior distribution can be used for parameter estimation:

where $p({\psi}_{k})$, $p({\Upsilon})$ and $p({\Sigma})$ are the corresponding prior distributions for the unknown parameters. In this paper, we consider zero-mean normal distributions with a large variance as prior distributions for the regression coefficients and association parameters. A vague inverse gamma distribution is considered as the prior distribution for the variance of errors. An inverse Wishart distribution with degrees of freedom equal to the dimension of the random effects and an identity scale matrix is utilized as the prior for the variance of the random effects. The parameters of the baseline hazard are assumed to follow a gamma distribution with a mean of 1 and a large variance. In the next section, we introduce our proposed two-stage approach for variable selection of longitudinal markers and covariates in a Bayesian paradigm.

The proposed method adopts a two-stage approach for variable selection in joint modeling, targeting the selection of markers (or time-varying covariates) and time-invariant covariates. This approach is specifically designed for handling multiple longitudinal markers and time-to-event outcomes, with the primary goal of avoiding bias caused by informative dropouts. This relies on the idea of TSJM approach of [5] but considering competing risks and replacing the estimation of a Cox model by partial likelihood in the second stage by a fully Bayesian approach allowing variable selection.
The first stage involves estimating $K$ one-marker joint models for the competing risks and each marker. In the second stage, we estimate cause-specific proportional hazard models including the time-invariant explanatory variables and the predicted current values of the linear predictor of each longitudinal markers obtained from stage 1. Therefore, the approach consists of two stages, outlined as follows:

Stage 1

For each $k$, $k=1,\ldots,K$, estimate the joint model for the longitudinal marker $k$ and competing risks. The one-marker JMs are estimated using the joint posterior distribution given by equation 2.3 for $K=1$. Then, predict the random effects $\hat{{b}}_{ik}$ and the linear predictor $\eta_{ik}(t|\hat{{\beta}}_{k},\hat{{b}}_{ik})$, for marker $k$ using the empirical mean of the posterior distribution $({b}_{ik}|{y}_{ik},t_{i},\delta_{i},\hat{{\theta}}_{k})$, where $\hat{{\theta}}_{k}$ represents the estimated parameters for the one-marker joint model for the $k$th longitudinal marker. This first stage allow the estimation of $\eta_{ik}(t|\hat{{\beta}}_{k},\hat{{b}}_{ik})$ without the bias due to informative dropout but neglecting the correlation between the markers..

Stage 2

In the second stage, the values of $\eta_{ik}(t|\hat{{\beta}}_{k},\hat{{b}}_{ik})$ are included as time-varying covariates in a proportional hazard model to estimate parameters ${\gamma}$ and ${\alpha}_{k}$, as well as the baseline hazard from (3). The cumulative hazard function of (3) can be approximated using Gaussian quadrature [9, 54]. In the Bayesian paradigm, the following joint posterior distribution is considered:

	$\displaystyle\pi({\lambda},\bm{\alpha}\mid\bm{t},\bm{\delta},\hat{\bm{b}},\bm{\omega})$	$\displaystyle\propto$	$\displaystyle\prod_{i=1}^{N}\prod_{l=1}^{L}\bigr{[}\lambda_{il}(t_{i}|{\omega}_{i},\hat{{b}}_{i},{\Upsilon})\bigr{]}^{I(\delta_{i}=l)}$
		$\displaystyle\times$	$\displaystyle\exp\bigr{(}-\sum_{l=1}^{L}\Lambda_{il}(t_{i}|{\omega}_{i},\hat{{b}}_{i},{\Upsilon})\bigr{)}\times p(\bm{\Upsilon}),$

where $\bm{\Upsilon}=({\lambda}_{l},{\alpha}_{l})$, such that, ${\lambda}_{l}=\{\lambda_{l1},\ldots,\lambda_{l\mathcal{K}}\}$ denotes the vector of the piecewise constant baseline hazard for cause $l$. The prior distributions for unknown parameters are the same as those considered for Equation 2.3.

As above, we consider the prior distributions ${\beta}_{k}\sim N_{p_{k}}(0,c\rm{I}_{p_{k}})$, where $c$ is a large value, $\sigma_{k}^{2}\sim{\rm{I\Gamma}}(a_{\sigma},b_{\sigma})$, ${\Sigma}\sim\rm{IWishart}({\mathcal{V}},\nu)$ and $\lambda_{l\kappa}\sim{\rm{\Gamma}}(a_{l\kappa},b_{l\kappa}),~{}{\kappa}=1,\ldots,\mathcal{K},~{}l=1,\ldots,L$. Our purpose is to implement variable selection for the regression coefficients ${\gamma}_{l},~{}l=1,\ldots,L$ of the cause-specific hazard model (3) and the association parameters ${\alpha}_{l},~{}l=1,\ldots,L$ of (3). Therefore, we consider two different prior distributions, including continuous spike and Dirac spike, for them.

To select significant variables in the covariates of the cause-specific hazard model and by considering Equation (11) or (20), the following hypothesis tests are conducted:

Also, to select significant longitudinal markers based on the association parameters and by considering Equation (16) or (24), the following hypothesis tests are considered:

In practice, when implementing the hypotheses (27) and (30) for CS, it is possible to consider the following equivalent hypotheses, respectively:

We consider two different Bayesian testing approaches for evaluating these hypotheses: a local Bayesian false discovery rate [18] and a Bayes factor [28].

The joint modeling of longitudinal measurements and time-to-event outcomes is a popular methodology for calculating dynamic prediction (DP). The DP represents the risk predictions for subject $i$ who has survived up to time point $s$, taking into account longitudinal measurements and other information available until time $s$. Consider the notation described in Section 2. Let the proposed joint models be fitted to a random sample, which is referred to as training data. For validation purposes, we are interested in predicting the risk of subjects in a new set of data called the validation set. The longitudinal measurements until time $s$ are denoted by ${\mathcal{Y}}_{i}(s)$ and given by $\{Y_{ik}(s_{ijk});0\leq s_{ijk}\leq s,j=1,\ldots,n_{i},k=1,\ldots,K\}$. Additionally, the vector of available explanatory variables until time $s$ is ${\mathcal{X}}_{i}(s)$ and given by $\{{w}_{i},{X}_{ik}(s_{ijk});0\leq s_{ijk}\leq s,j=1,\ldots,n_{i},k=1,\ldots,K\}$. Let $\delta_{i}=1$ denote the event of interest, then, for the landmark time $s$ and the prediction window $t$, the risk probability is as follows:

where ${{\theta}}$ is a vector of all unknown parameters of the model, including ${{\theta}}_{k}$ which is estimated at stage 1, and ${\bm{\Upsilon}}$ which is estimated at stage 2, and $T_{i}^{*}$ represents the true event time for the $i$th subject. The predicted probability (38) can be computed as:

where $S(\cdot)$ represents the survival probability to survive to both events given random effects. The first order approximation for estimating this quantity [39, 5] is given as follows:

where $\hat{{\theta}}=\{\hat{{\theta}}_{1},\ldots,\hat{{\theta}}_{K},\hat{\bm{\Upsilon}}\}$ are estimated based on the training data. To predict $\hat{{b}}_{i}$, we consider the one-marker joint models for each marker and predict the $p_{k}$-dimensional random effect ${b}_{ik}$ given $(T^{*}_{i}>s,{\mathcal{Y}}_{i}(s),{\mathcal{X}}_{i}(s),\hat{{\theta}}_{k})$. This is denoted by $\hat{{b}}_{ik}$, where $\hat{{\theta}}_{k}$ represents the estimated parameters for the one-marker joint model in the first stage of VSJM. Merge the predicted random effects to obtain the $p=\sum_{k=1}^{K}p_{k}$-dimensional predicted random effect $\hat{{b}}_{i}=(\hat{{b}}_{i1},\ldots,\hat{{b}}_{iK})$.
For Monte-Carlo approximation of (39), we follow the same strategy as those discussed in [62]. For this aim, consider a realization of the generated MCMC sample for estimating the unknown parameters from the one-marker JM for marker $k$ in the training data, denoted ${\theta}_{k}^{(\ell)},~{}\ell=1,\ldots,\mathcal{L}$. Draw $M$ values ${b}_{ik}^{(\ell,m)}$ from the posterior distribution $({b}_{ik}|T^{*}_{i}>s,{\mathcal{Y}}_{i}(s),{\mathcal{X}}_{i}(s),{\theta}_{k}^{(\ell)})$. For $~{}\ell=1,\ldots,\mathcal{L},~{}m=1,\ldots,M,$

Calculate any desired quantity from the generated chains, such as percentiles, to determine the credible interval.

To reduce estimation biases resulting from variable selection, we propose incorporating an additional stage to calculate dynamic predictions. After variable selection using CS or DS prior, we recommend re-estimating the proportional hazard model by substituting CS or DS with non-informative normal priors for the association parameters of the selected markers and the regression coefficients of the selected covariates. Then calculate the DP as described in Section 3.1.

We use the Area Under the Receiver Operating Characteristics Curve (AUC) and the Brier Score (BS) to assess the predictive performance of prediction models [50]. AUC represents the probability that the model correctly ranks a randomly chosen individual with the disease higher than a randomly chosen individual without the disease based on their predicted risks Additionally, BS measures the preduction accuracy by quantifying the mean squared error of prediction. For censored and competing time-to-event and time-dependent markers, a non-parametric inverse probability of censoring weighted estimator of AUC and BS was proposed by [8] and utilized in this paper.

A simulation study was conducted to assess the effectiveness of the proposed VSJM approach for variable selection and risk prediction. To evaluate the effectiveness of the proposed approaches for variable selection, we examine three distinct criteria: the true positive rate (TPR), which measures the proportion of correctly selected variables among all significant variables; the false positive rate (FPR), which measures the proportion of incorrectly selected variables among all irrelevant variables; and the Matthews correlation coefficient [33] (MCC). The latter is defined as follows:

where TP, TN, FP, and FN represent the number of true positives, true negatives, false positives, and false negatives, respectively. The MCC and TPR are expected to reach 1, while the FPR is expected to be close to zero for optimal performance. BF and lBFDR are used for variable selection, with the results based on a threshold of 1 for BF and a threshold of 0.05 for lBFDR.
The simulation scenarios varied based on the strength of the association between the markers and the events, the residual variances for the markers that quantify measurement error, and the correlation between the markers that is neglected in the first stage of the approach. The performance of the approaches is also evaluated based on their predictive abilities using AUC and BS. The predictive performance of the approach is compared with that of the true model using the true values of the parameters and random effects (referred to as ”Real”) as well as the oracle model (true model, where parameters are estimated by accounting for the multi-marker joint modeling).
Data was generated using multi-marker joint models with $K=10$ longitudinal markers and a competing risks model with two causes. The longitudinal data was generated from the following model:

where $\varepsilon_{ikt}\sim N(0,\sigma_{k}^{2})$, ${b}_{i}=(b_{01i},b_{11i},b_{02i},b_{12i},\ldots,b_{0,10,i},b_{1,10,i})^{\top}\sim N({0},{\Sigma})$, $i=1,\ldots,n$. Let $k=1,2,\ldots,K=10$, the measurement times $t=0,0.1,0.2,\ldots,2$ and ${\beta}_{k}=(-0.5,0.5)^{\top},k=1,2,..,10$.
The time-to-event outcome was simulated using a cause-specific hazard model that depends on the current value of the $K$ markers and on 24 time-invariant covariates as follows:

where ${\omega}_{il}=(x_{1i},w_{1i},\ldots,w_{24,i})^{T}$. Half of the covariates were independently generated from a standard normal distribution, while the others were independently generated from a binary distribution with a success probability of 0.5. To simulate time-to-event data, we apply the inverse of the cumulative hazard function [3] using the inverse probability integral transform method [43].
The covariance matrix of the random effects a structure is considered as follows:

For each scenario, $100$ data sets of $n=1000$ subjects were generated and split into a training sample of 500 subjects and a validation of 500 subjects. The investigation focused on estimating parameters on the training sample using various approaches. The performance of these approaches for risk prediction was then evaluated using the validation samples. To achieve this goal, we consider four landmark times at $s=0,0.25,0.5,0.75$, using $t=0.25$ as window for prediction.

Among the 10 generated longitudinal markers and the 24 generated time-invariant covariates, only 4 markers and 4 covariates (2 Gaussian and 2 binary) were truly associated with the two event risks, sharing identical regression parameters for both events. To evaluate the impact of increasing measurement errors, association parameters, and marker correlation on the performance of the three methods, we considered different values for these parameters, as described below:

1. 1.

   $\sigma_{k}^{2}=0.5,1$ for $k=1,2,\ldots,10$,

2. 2.

   ${\gamma}_{l}=(\gamma^{*},-\gamma^{*},\underbrace{0,\ldots,0}_{10},\gamma^{*},-\gamma^{*},\underbrace{0,\ldots,0}_{10})^{T},~{}l=1,2$ and ${\alpha}_{l}=(-\alpha^{*},-\alpha^{*},\alpha^{*},\alpha^{*},0,0,0,0,0,0)^{\top}$, where $\gamma^{*}=\alpha^{*}=0.5,1$.

3. 3.

   $\rho=0.1,~{}0.7$.

Based on the insights obtained from both Table 1, which pertains to variable selection, and Table 2, which focuses on risk prediction, regarding the performance of the VSJM, we can conclude:

• •

  For all scenarios, AUC and BS computed on the validation set after Bayesian variable selection on the training set were close to those of the oracle model.

• •

  Increasing the between markers correlation values tends to decrease the performance of the approach in terms of TPR, FPR, and MCC, albeit to an acceptable extent.

• •

  Increasing the residual variance of the longitudinal markers ($\sigma^{2}$) leads to a decrease in the performance of the approach.

• •

  Enhanced values of association parameters and regression coefficients tend to improve the approach’s performance.

• •

  Although the predictive performance of BF and lBFDR is similar, BF generally outperforms lBFDR in terms of total variable selection, except in the easiest scenario characterized by high association and low measurement error.