Liangyuan Hu, Department of Biostatistics and Epidemiology, Rutgers University
School of Public Health, 683 Hoes Lane West, Piscataway, NJ, USA 08854.

# A new tool for clustered survival data and multiple treatments: Estimation
of treatment effect heterogeneity and variable selection

Liangyuan Hu11affiliationmark:  11affiliationmark:  Department of
Biostatistics and Epidemiology, Rutgers University School of Public Health,
USA<EMAIL_ADDRESS>
###### Abstract

A new tool, riAFT-BART, was recently developed to draw causal inferences about
population treatment effect on patient survival from clustered and censored
survival data while accounting for the multilevel data structure. The
practical utility of this tool goes beyond the estimation of population
average treatment effect. In this work, we exposit how riAFT-BART can be used
to solve two important statistical questions with clustered survival data:
estimating the treatment effect heterogeneity and variable selection.
Leveraging the likelihood-based machine learning, we describe a way in which
we can draw posterior samples of the individual survival treatment effect from
riAFT-BART model runs, and use the drawn posterior samples to perform an
exploratory treatment effect heterogeneity analysis to identify subpopulations
who may experience differential treatment effects than population average
effects. There is sparse literature on methods for variable selection among
clustered and censored survival data, particularly ones using flexible
modeling techniques. We propose a permutation based approach using the
predictor’s variable inclusion proportion supplied by the riAFT-BART model for
variable selection. To address the missing data issue frequently encountered
in health databases, we propose a strategy to combine bootstrap imputation and
riAFT-BART for variable selection among incomplete clustered survival data. We
conduct an expansive simulation study to examine the practical operating
characteristics of our proposed methods, and provide empirical evidence that
our proposed methods perform better than several existing methods across a
wide range of data scenarios. Finally, we demonstrate the methods via a case
study of predictors for in-hospital mortality among severe COVID-19 patients
and estimating the heterogeneous treatment effects of three COVID-specific
medications. The methods developed in this work are readily available in the R
package riAFTBART.

###### keywords:

Missing at random, Variable importance, Tree ensemble, Bootstrap
imputation,Variable selection

##  1 Introduction

Large-scale healthcare datasets garnered from multiple facilities are
increasingly available, offering fertile ground for innovative investigation
into emerging medical research questions. The inherent hierarchical data
structure, in conjunction with the censored survival outcome, which is of
strong clinical interest, poses considerable challenges for statistical
analysis. Hu et al. hu2022flexible  recently developed a flexible approach,
riAFT-BART to estimate the population average treatment effect on patient
survival, while accounting for the clustering structure of the observations.
The method riAFT-BART stands for random-intercept accelerated failure time
model with Bayesian additive regression trees (BART).(chipman2010bart, ) The
utility of this new tool, which leverages Bayesian machine learning, goes
beyond the estimation of population average treatment effect. In this article,
we exposit how riAFT-BART can be used to solve two important statistical
questions with clustered survival data: estimating the treatment effect
heterogeneity and variable selection.

An immediate extension of the new tool, riAFT-BART, beyond the population
average treatment effect is the estimation of heterogeneous treatment effect.
Prior work (hu2021estimating, ) has shown that compared to algorithm-based
machine learning, leveraging the likelihood-based machine learning technique
BART can provide more accurate estimation of the individual survival treatment
effect, which in turn can be used to assess the treatment effect
heterogeneity. To the best of our knowledge, no work has yet discussed how
BART can be used to estimate the heterogeneous effect of multiple treatments
on patient survival using clustered survival observations.

Another important utility of riAFT-BART is for variable selection with
clustered survival data, on which there is sparse literature. Fan et
al.(fan2002variable, ) proposed a family of variable selection methods for the
gamma frailty model based on a nonconcave penalized likelihood approach using
the Newton-Raphson method. Androulakis et al.(androulakis2012estimation, )
extended the penalized Cox or Gamma frailty model to a class of frailty models
and proposed a generalized form of the likelihood function to allow the use of
different continuous distributions for the frailty term. Utazirubanda et
al.(utazirubanda2021variable, ) proposed group LASSO with gamma-distributed
frailty for variable selection in Cox regression. However, no corresponding
software is available for these methods. Bender et al. (bender2018generalized,
) proposed a piece-wise exponential additive mixed model (PEAMM) for clustered
survival data, and Marra and Wood (marra2011practical, ) described a practical
approach by which shrinkage penalties can be added to the smoothing parameters
in PEAMM and variable selection can be performed by using restricted maximum
likelihood estimation. Rondeau et al. (rondeau2003maximum, ) proposed to
conduct variable selection by maximizing the penalized likelihood using the
Marquardt algorithm, (marquardt1963algorithm, ) which is a combination of the
Newton-Raphson algorithm and the steepest descent algorithm. Ha et
al.(ha2014variable, ) proposed a penalized h-likelihood estimation approach
for variable selection of fixed effects using semiparametric frailty models.
These methods rely on parametric assumptions about the exact relationships
between covariates and survival outcomes. Misspecifying the parametric forms
of covariate-outcome relationships can reverberate through the variable
selection procedure and yield unsatisfactory results. (hu2021variable, ;
bleich2014variable, )

Flexible machine learning modeling techniques can help relax the parametric
assumptions and improve variable selection results. Ishwaran et
al.(ishwaran2010high, ) derived the distribution of the minimal depth of a
maximal subtree, which measures the predictiveness of a variable in a survival
tree, and used it for high-dimensional variable selection using random
survival forests. Lagani et al.(lagani2010structure, ) proposed an algorithm
based on the theory of Bayesian networks and the Markov blanket for variable
selection suitable for high-dimensional and right-censored survival data.
However, neither of these two machine learning based methods accommodates
clustered survival observations.

In this work, we describe ways in which riAFT-BART can be used to (i) detect
and estimate the treatment effect heterogeneity, and to (ii) select important
predictors in the presence of missing data with clustered survival
observations. The rest of the paper is organized as follows. Section 2
describes the use of riAFT-BART for estimating individual survival treatment
effects and exploring subgroups who may experience enhanced or reduced
treatment effect than population average treatment effect. Section 3 proposes
a variable selection method using riAFT-BART when covariates are subject to
missing data. In Section 4, we conduct an expansive simulation to assess the
practical operating characteristics of riAFT-BART in estimating treatment
effect heterogeneity and performing variable selection. Section 5 illustrates
the methods using a COVID-19 dataset drawn from the Mount Sinai Medical
Center, and Section 6 provides a discussion.

##  2 Estimation of Treatment Effect Heterogeneity

###  2.1 Notation

We maintain notation used in Hu et al. (hu2022flexible, ) Consider a
nonexperimental study with $K$ clusters, each having treated $n_{k}$
individuals, indexed by $i=1,\ldots,n_{k},k=1,\ldots,K$. There are $J$
possible treatment options, denoted by
$A\in\mathscr{A}=\\{a_{1},\ldots,a_{J}\\}$. Denote pre-treatment measured,
$L$-dimensional covariates by $\bm{X}_{ik}=\\{X_{ik1},\ldots,X_{ikL}\\}$ for
each individual $i$ in cluster $k$, whose failure time $T_{ik}$ may be right
censored at $C_{ik}$. The observed outcome consists of
$Y_{ik}=\min(T_{ik},C_{ik})$ and the censoring indicator
$\Delta_{ik}=I(T_{ik}<C_{ik})$. Let $V_{k}$ be the cluster indicators.
Proceeding with the counterfactual framework, the counterfactual failure time
under treatment $a_{j}$ for individual $i$ in cluster $k$ is defined as
$T_{ik}(a_{j})$, $\forall a_{j}\in\mathscr{A}$, and the counterfactual
censoring time under treatment $a_{j}$ is defined as $C_{ik}(a_{j})$.
Throughout, we maintained the standard assumptions: (hu2022flexible, ;
hu2021estimating, ; hu2019causal, ; chen2001causal, ; arpino2016propensity, )
consistency, weak unconfoundedness, positivity and covariate-dependent
censoring, for drawing causal inference with observational clustered survival
data. Detailed assumptions and implications of assumptions can be found in Hu
et al. (hu2022flexible, )

###  2.2 Method

Hu et al. (hu2022flexible, ) adapted BART into a random-intercept accelerated
failure time model,

$\displaystyle\begin{split}&\log
T_{ik}=f(A_{ik},\bm{X}_{ik})+b_{k}+\epsilon_{ik},\\\
b_{k}\mathbin{\overset{i.i.d}{\kern
0.0pt\sim}}&N(0,\alpha_{k}\tau^{2}),\quad\epsilon_{ik}\mathbin{\overset{i.i.d}{\kern
0.0pt\sim}}N(0,\sigma^{2}),\quad b_{k}\perp\epsilon_{ik},\end{split}$ (1)

where $f(A_{ik},\bm{X}_{ik})$ is an unspecified function relating treatment
assignment $A_{ik}$ and pretreatment confounders $\bm{X}_{ik}$ to survival
times $T_{ik}$, $b_{k}$’s are the random intercepts capturing cluster-specific
main effects, and $\epsilon_{ik}$ is the residual term. The unknown function
$f$ is flexibly modeled by BART via a sum of shallow trees. (hu2020estimation,
; hu2021estimation, ) A mean-zero normal distribution and independence are
assumed for $b_{k}$ and $\epsilon_{ik}$. A redundant parameter $\alpha_{k}$
was used for the variance of $b_{k}$ as a parameter expansion technique
(gelman2008using, ) to improve computational performance. Regularizing priors
are placed on the tree parameters to keep the impact of each individual tree
on the overall fit small and thus prevent overfitting. (chipman2010bart, )
Sensible priors are used for parameters $\sigma^{2}$, $\tau^{2}$ and
$\alpha_{k}$. (hu2022flexible, ; henderson2020individualized, ) To deal with
right censoring, data augmentation is used. That is, the unobserved survival
times are imputed from a truncated normal distribution in each Gibbs
iteration. A Metropolis within Gibbs procedure (hu2022flexible, ) was proposed
to draw posterior inferences about the population average treatment effect on
patient survival, which is defined as

$\theta_{a_{j},a_{j^{\prime}}}=\frac{1}{N}\sum_{k=1}^{K}\sum_{i=1}^{n_{k}}E\left[T_{ik}(a_{j})-T_{ik}(a_{j^{\prime}}){\,|\,}\bm{X}_{ik}=\bm{x}_{ik},V_{k}=v_{k}\right],\;\forall
a_{j},a_{j^{\prime}}\in\mathscr{A}.$

Via outcome modeling, the pairwise population average treatment effect between
$a_{j}$ and $a_{j^{\prime}}$ can be estimated as

$\displaystyle\hat{\theta}_{a_{j},a_{j^{\prime}}}=\frac{1}{D}\sum_{d=1}^{D}\frac{1}{n_{k}K}\sum_{k=1}^{K}\sum_{i=1}^{n_{k}}\left\\{\tilde{f}^{d}\left(a_{j},\bm{x}_{ik}\right)-\tilde{f}^{d}\left(a_{j^{\prime}},\bm{x}_{ik}\right)\right\\},$
(2)

where $\tilde{f}^{d}$ is the $d$th draw from the posterior distribution of
$f$. Inferences can be obtained based on the $D$ posterior average treatment
effects,
$(1/n_{k}K)\sum_{k=1}^{K}\sum_{i=1}^{n_{k}}\left\\{\tilde{f}^{d}\left(a_{j},\bm{x}_{ik}\right)-\tilde{f}^{d}\left(a_{j^{\prime}},\bm{x}_{ik}\right)\right\\},d=1,\ldots,D$.

Turning to the estimation of treatment effect heterogeneity. We first define
the _individual survival treatment effect_ (ISTE) for individual $i$ in
cluster $k$ as

$\displaystyle\zeta_{a_{j},a_{j^{\prime}}}(\bm{x}_{ik})=E\left[\log
T_{ik}(a_{j})-\log T_{ik}(a_{j^{\prime}}){\,|\,}\bm{x}_{ik},v_{k}\right],$ (3)

which will be used as the building block for estimating the heterogeneous
treatment effect. From the fitted Model (LABEL:eq:_riAFT-BART), we can
estimate $\zeta_{a_{j},a_{j^{\prime}}}(\bm{x}_{ik})$ as

$\displaystyle\hat{\zeta}_{a_{j},a_{j^{\prime}}}(\bm{x}_{ik})=\tilde{f}(a_{j},\bm{x}_{ik})-\tilde{f}(a_{j^{\prime}},\bm{x}_{ik}),$

where $\tilde{f}$ is the mean of $D$ draws from the posterior predictive
distribution of $f$. The ISTE estimates
$\\{\hat{\zeta}_{a_{j},a_{j^{\prime}}}(\bm{x}_{ik}),i=1,\ldots,n_{k},k=1,\ldots,K\\}$
can then be used to explore the degree of treatment effect heterogeneity. By
the “virtual twins” approach, (foster2011subgroup, ) the
$\\{\hat{\zeta}_{a_{j},a_{j^{\prime}}}(\bm{x}_{ik})\\}$ would be regressed on
the predictors via a tree diagram to identify the predictor subspaces that
explain away differences in $\\{\hat{\zeta}(\bm{x}_{ik})\\}$. In a similar
spirit, the “fit-the-fit” strategy has been used to identify possible
subgroups defined by combination rules of covariates that have differential
treatment effects. (logan2019decision, ; hu2021estimating, ;
hu2021estimatingh, ) Mathematically, we wish to find a covariate subspace
$\mathcal{S}_{a_{j},a_{j^{\prime}}}=\\{S^{1}_{a_{j},a_{j^{\prime}}},S^{2}_{a_{j},a_{j^{\prime}}},\ldots,S^{H}_{a_{j},a_{j^{\prime}}}\\}$,
where $h\in\\{1,\ldots,H\\}$ is the number of combination rules of covariates,
such that the average treatment effects across individuals in this region are
considerably different than the population average treatment effect. Denote
the measure of subgroup treatment effect between $a_{j}$ and $a_{j^{\prime}}$
as
$\mathcal{M}(\mathcal{S}_{a_{j},a_{j^{\prime}}})=\\{M(S^{1}_{a_{j},a_{j^{\prime}}}),\ldots,M(S^{H}_{a_{j},a_{j^{\prime}}})\\}$,
and

$\displaystyle
M(S^{h}_{a_{j},a_{j^{\prime}}})=E\left[\zeta_{a_{j},a_{j^{\prime}}}(\bm{x}_{ik}){\,|\,}\bm{x}_{ik}\in
S^{h}_{a_{j},a_{j^{\prime}}}\right].$ (4)

We adopt the “fit-the-fit” strategy to estimate
$\mathcal{S}_{a_{j},a_{j^{\prime}}}$ and
$\mathcal{M}(\mathcal{S}_{a_{j},a_{j^{\prime}}})$. The algorithm proceeds with
the following steps:

  1. 1.

Add each candidate predictor separately to the random forests
model(breiman2001random, ) using
$\\{\hat{\zeta}_{a_{j},a_{j^{\prime}}}(\bm{x}_{ik})\\}$ as the outcome, and
record the model fit measured by $R^{2}$; include the predictor corresponding
to the largest $R^{2}$ in the model;

  2. 2.

Add each of the remaining candidate predictors separately to the previous
step’s model and retain the predictor corresponding to the largest $R^{2}$
improvement;

  3. 3.

Repeat step 2 until the percent improvement in $R^{2}$ is less than 1%.

Different from prior work, we use random forests instead of a single tree to
regress the ISTE estimates on predictors for improved modeling accuracy and
stability. To interpret the tree ensembles in the final random forests model
yielded in Step 3, we apply the inTrees algorithms (deng2019interpreting, ) to
translate the composite rules into a simple rule model. The inTrees algorithms
extract, measure and rank rules, prune redundant rules, detect frequent
variable interactions, and summarize rules into a simple tree diagram model
that can be used for predicting new data. The branch decision rules in the
final tree diagram model sending individuals to the terminal nodes suggest the
predictor subspace $\hat{\mathcal{S}}_{a_{j},a_{j^{\prime}}}$; and subgroup
treatment effects that are substantially different from the population effect
$\theta_{a_{j},a_{j^{\prime}}}$ can then be estimated by averaging the ISTE
estimates among individuals falling into each of terminal nodes,

$\hat{M}(\hat{S}^{h}_{a_{j},a_{j^{\prime}}})=\frac{1}{{\,|\,}\hat{S}^{h}_{a_{j},a_{j^{\prime}}}{\,|\,}}\sum\limits_{i,k:\bm{x}_{ik}\in\hat{S}^{h}_{a_{j},a_{j^{\prime}}}}\hat{\zeta}_{a_{j},a_{j^{\prime}}}(\bm{x}_{ik}),$

where $\mid\hat{S}^{h}_{a_{j},a_{j^{\prime}}}\mid$ is the size of
$\hat{S}^{h}_{a_{j},a_{j^{\prime}}}$.

##  3 Variable Selection

Variable selection has been a longstanding statistical problem and studied
from both the classical and Bayesian perspective. In the field of causal
inference, variable selection methods can aid in confounder selection.
(vanderweele2019principles, ) Variable selection methods are largely based on
parametric models, in which the exact relationships between the response and
predictor variables need to be explicitly specified. However, incorrectly
specifying the parametric models may produce unsatisfactory variable selection
results. (bleich2014variable, ; hu2021variable, ; lin2022a, ) For example,
noise variables may be selected and important variables ignored. Flexible
machine learning models can better represent functional forms of arbitrary
complexity, thus mitigate parametric assumptions and attendant errors in
variable selection. We develop a way in which the new tool riAFT-BART can be
used for variable selection among clustered survival data. Notation defined in
Section 2.1 is largely retained with the exception of treatment assignment
$A_{ik}$, as the focus now is variable selection.

###  3.1 Random-intercept accelerated failure time model with Bayesian
additive regression trees

We propose a permutation-based approach for variable selection. Following
prior work, bleich2014variable ; hu2021variable ; lin2022a  we use the
variable inclusion proportion (VIP) of each predictor variable as the
foundation to determine how important a predictor variable is. The VIP score
is supplied by the BART model, and is computed as the proportion of times each
predictor is chosen as a splitting rule divided by the total number of
splitting rules appearing in the model. (hu2020ranking, ) We note that in the
riAFT-BART model (LABEL:eq:_riAFT-BART), the predictor variables
$\bm{X}_{ik}$, along with treatment assignment $A_{ik}$, is included in the
unknown function $f(A_{ik},\bm{X}_{ik})$, which is being modeled by BART.
Thus, in each iteration of the riAFT-BART sampling algorithm, (hu2022flexible,
) we can extract the VIP for each predictor variable in the step where a BART
model is fit to regress $\log y^{cent,c}_{ik}-b_{k}$ on
$\\{A_{ik},\bm{X}_{ik}\\}$, where $y^{cent,c}_{ik}$ denotes the centered and
complete survival time. Note that right censoring is addressed via data
augmentation. Detailed riAFT-BART sampling algorithm appears in Web Section 1.
For each predictor variable $X_{l},\;l\in\\{1,\ldots,L\\}$ (we suppress
subscripts $i$ and $k$ for brevity of notation), we compute its VIP by
averaging the VIP scores across $D$ posterior draws, (bleich2014variable, )
$\text{VIP}_{x_{l}}=\frac{1}{D}\sum_{d=1}^{D}\widetilde{\text{VIP}}^{d}_{x_{l}},d=1,\ldots,D$,
where $\widetilde{\text{VIP}}^{d}_{x_{l}}$ is the $d$th draw from the
posterior distribution of the VIP for $X_{l}$. In our simulation (Section 4),
we set $D=4500$ with the first 1000 discarded as burn-in. Trace plot was used
to check the convergence of the posterior distribution of the VIP for $X_{l}$.

Building on prior literature, (bleich2014variable, ; huang2010variable, ) we
propose a non-parametric and permutation-based approach for variable
selection, by which we first permute the event times together with the
censoring indicators and then establish the thresholds for variable selection
using the VIP from observed data and the permutated data. Specifically, we
first create $P$ permutations of the vectors of observed survival times and
event indicators
$(\bm{y},\bm{\Delta})=\\{Y_{11},\Delta_{11},\ldots,Y_{n_{k}K},\Delta_{n_{k}K}\\}$:
$(\bm{y}_{1}^{*},\bm{\Delta}_{1}^{*}),(\bm{y}_{2}^{*},\bm{\Delta}_{2}^{*})\ldots,(\bm{y}_{P}^{*},\bm{\Delta}_{P}^{*})$.
Then for each of the permuted outcome vectors
$(\bm{y}_{p}^{*},\bm{\Delta}_{p}^{*})$, we fit an riAFT-BART model with
$(\bm{y}_{p}^{*},\bm{\Delta}_{p}^{*})$ as the response and the original
$\\{\bm{X}_{ik}\\}$ as predictor variables and $V_{k}$ as cluster indicator.
This permutation strategy removes any dependency between the predictors
$\\{\bm{X}_{ik}\\}$ and the outcomes $\\{(Y_{ik},\Delta_{ik})\\}$ while
keeping possible dependencies among the predictor variables. From the riAFT-
BART run using each permuted outcome vector
$(\bm{y}_{p}^{*},\bm{\Delta}_{p}^{*})$, we retain the VIP for each predictor
variable. We use $\text{VIP}_{x_{l},p}^{*}$ to denote the VIP score from
riAFT-BART for predictor $X_{l}$ from the $p$th permuted outcomes, and use the
notation $\text{VIP}_{\bm{X},p}^{*}$ for the vector of all $L$ VIPs from the
$p$th permuted outcomes,
$(\text{VIP}_{x_{1},p}^{*},\ldots,\text{VIP}_{x_{L},p}^{*})$. As suggested by
Bleich et al. (bleich2014variable, ), we use the VIPs across all $P$
permutations $\text{VIP}^{*}_{\bm{X},1},\ldots,\text{VIP}^{*}_{\bm{X},P}$ as
the “null” distribution for the VIPs from the unpermuted outcomes
$(\bm{y},\bm{\Delta})$, $\text{VIP}_{x_{1}},\ldots,\text{VIP}_{x_{L}}$. The
permutation null distribution will be used to determine which variables should
be selected or deemed as important. We adopt the “local” threshold strategy,
which has been shown in several works (bleich2014variable, ; hu2020machine, ;
hu2020tree, ; hu2021variable, ; lin2022a, ) to have better performance than
more stringent strategies such as “global max” for large-scale healthcare data
in which the ratio of sample size and number of predictors is relatively
large. More stringent strategies compare the VIP of a predictor
$\text{VIP}_{x_{l}}$ to the permutation distribution across all predictor
variables $\text{VIP}_{\bm{X},1}^{*},\ldots,\text{VIP}_{\bm{X},P}^{*}$. By
contrast, a local threshold compares $\text{VIP}_{x_{l}}$ with its own
permutation null distribution
$\text{VIP}_{x_{l},1}^{*},\ldots,\text{VIP}_{x_{l},P}^{*}$. Specifically, we
only select predictor $X_{l}$ if $\text{VIP}_{x_{l}}$ is greater than the
$1-\alpha$ quantile of the permutation distribution
$\text{VIP}_{x_{l},1}^{*},\ldots,\text{VIP}_{x_{l},P}^{*}$. Following prior
work, (hu2021variable, ; lin2022a, ) we set $\alpha=0.05$ in our simulation
and case study, because this threshold value has been shown to yield good
variable selection results with a balance between selecting useful and
selecting noise predictors.

###  3.2 Comparison methods

####  3.2.1 Piece-wise Exponential Additive Mixed Model

For a PEAMM, (bender2018generalized, ) the hazard rate at time $t$ for
individual $i$ in cluster $k$ with covariates $\bm{X_{ik}}$ and cluster
indicator $V_{k}$ is given by

$\displaystyle\lambda\left(t{\,|\,}\bm{X}_{ik},V_{k}\right)=\exp\left\\{g_{0}\left(t_{q}\right)+\sum_{m=1}^{M}g_{m}\left(\bm{x}_{ik,m}\right)+b_{k}\right\\}\forall
t\in(\kappa_{q-1},\kappa_{q}),$ (5)

where $g_{0}\left(t_{q}\right)$ represents the log-baseline hazard rate,
$g_{m}\left(\bm{x}_{ik,m}\right)$ is the basis functions of natural cubic
splines associated with the covariates $\bm{x}_{ik}$, $b_{k}$ is the random
intercept term capturing the cluster-specific effect, and $\kappa_{q-1}$ and
$\kappa_{q}$ are the start and stop times of the censored survival dataset
represented in the (start, stop] counting process form. For variable selection
using PEAMM, following the practical approach by Marra and Wood et
al.,(marra2011practical, ) a shrinkage penalty term can be added to each basis
term $g_{m}\left(\bm{x}_{ik,m}\right)$ of model (5):

$\displaystyle\lambda\left(t{\,|\,}\bm{X}_{ik},V_{k}\right)=\exp\left\\{g_{0}\left(t_{q}\right)+\sum_{m=1}^{M}\left(g_{m}\left(\bm{x}_{ik,m}\right)+\int
g^{\prime}_{m}(\bm{x}_{ik,m})^{2}d\bm{x}_{ik,m}\right)+b_{k}\right\\}\forall
t\in(\kappa_{q-1},\kappa_{q}),$

so that when the whole spline basis terms for a predictor $x_{ikl}$ are shrunk
to zero, the variable $x_{ikl}$ will be selected out of the model. The
shrinkage penalties, along with all other model coefficients, can be estimated
using restricted maximum likelihood estimation.

####  3.2.2 Frailty Models by Regularization Methods

Given the common unobserved frailty $u_{k}$ for cluster $k$, the conditional
hazard function at time $t$ for individual $i$ in cluster $k$ is

$\displaystyle\lambda(t{\,|\,}\bm{X}_{ik},V_{k})=u_{k}\lambda_{0}(t)\exp(\bm{x_{ik}}^{T}\bm{\beta})$
(6)

where $\lambda_{0}(t)$ is an unspecified baseline hazard function and
$\bm{\beta}=(\beta_{1},\beta_{2},...,\beta_{L})^{T}$ is a vector of regression
parameters for $\bm{x_{ik}}$. The common unobserved frailty $u_{k}$ is assumed
to follow a gamma distribution because of its mathematical convenience.
Rondeau et al. (rondeau2003maximum, ) proposed to conduct variable selection
by maximizing the penalized likelihood,

$l(\lambda(t{\,|\,}\bm{X}_{ik},V_{k}))-\alpha\int_{0}^{\infty}\lambda_{0}^{\prime\prime
2}(t)dt,$

where $l(\lambda(t{\,|\,}\bm{X}_{ik},V_{k}))$ is the full marginal log-
likelihood for the frailty model (6) and $\alpha$ is the smoothing parameter
which controls the trade-off between the data fit and the smoothness of the
baseline hazard function. The regression parameters can be estimated by the
robust Marquardt algorithm(marquardt1963algorithm, ) – a combination of the
Newton Raphson algorithm and the steepest descent algorithm. Do Ha et
al.(ha2014variable, ) instead proposed to perform variable selection by
maximizing the penalized profile h-likelihood,

$l^{*}(\lambda(t{\,|\,}\bm{X}_{ik},V_{k}))-\sum_{l=1}^{L}B_{\gamma}(|\beta_{l}|),$

where $l^{*}(\lambda(t{\,|\,}\bm{X}_{ik},V_{k}))$ is the profiled h-likelihood
of the frailty model, $B_{\gamma}(|\beta_{l}|)$ is a penalty function that
controls model complexity using the tuning parameter $\gamma$, with the
optimal value determined by the Bayesian information criterion (BIC). Three
penalty functions were considered for $B_{\gamma}(|\beta_{l}|)$: least
absolute shrinkage and selection operator, (tibshirani1996regression, )
smoothly clipped absolute deviation, (fan2001variable, ) and h-likelihood.
(lee2014new, )

####  3.2.3 Backward stepwise selection

We also implement a classical variable selection approach through a recursive
procedure of backward stepwise selection based on statistical testing, as
described in Wood et al. (wood2008should, ) The backward stepwise selection
starts with a random-intercept Cox regression model (cortinas2005version, )
with all potential predictors and removes from the model the predictor that
has the least impact on the fit determined by the $\alpha$ significance level.
The forward selection process then checks whether removed variables should be
added back into the model at $\alpha(1-\epsilon)$ for $\epsilon$ small
significance level. We use $\alpha=0.05$ as suggested by Wood et al.
(wood2008should, )

###  3.3 Variable selection with missing data

Large-scale healthcare data are susceptible to the missing data issue. Missing
values in predictor variables pose challenges for variable selection. Long and
Johnson (long2015variable, ) proposed a general resampling approach that
combines bootstrap imputation and randomized lasso based variable selection
procedure, and used simulations to demonstrate that this approach had better
performance compared with several existing methods for variable selection in
the presence of missing data. Recent work by Hu et al. (hu2021variable, )
shows that combining bootstrap imputation with flexible machine learning based
variable selection methods can substantially improve variable selection
results over methods that rely on parametric modeling of the covariate-outcome
relationship. In this work, we combine bootstrap imputation and variable
selection, and compare each of the methods considered on the ability to select
predictor variables that are truly associated with the outcome, or useful
predictor variables, among _incomplete_ clustered survival data. To stay
consistent with the literature wood2008should ; bleich2014variable ;
hu2021variable  so that our methods and results can be compared with previous
work in similar contexts, we use the following four performance metrics
commonly used in the variable selection literature. A good method should have
high values of precision, recall and $F_{1}$ , and low Type I error.

  1. (1)

$\text{Precision}=\text{TP}/(\text{TP}+\text{FP})$, where TP and FP are
respectively the number of true positive and false positive selections. A true
positive selection is when a useful predictor was selected; and a false
positive selection is when a noise predictor was selected. The _precision_ of
a variable selection method is the proportion of truly useful predictors among
all selected predictors.

  2. (2)

$\text{Recall}=\text{TP}/(\text{TP}+\text{FN})$, where FN is the number of
false negative selections. A false negative selection is when a truly useful
predictor was not selected. The _recall_ of a variable selection method is the
proportion of truly useful variables selected among all useful variables. This
is sometimes referred to as the _power_ in the literature. wood2008should

  3. (3)

$F_{1}=2\text{
Precision}\times\text{Recall}/(\text{Precision}+\text{Recall}).$ The $F_{1}$
score measures a method’s ability to avoid selecting irrelevant predictors
(precision) with its ability to identify the full set of useful predictors
(recall).

  4. (4)

Type I error. The _Type I error_ measures the overall probability of a method
incorrectly selecting the noise predictor – the mean of the probabilities that
the method incorrectly selects each of the noise predictors.

Among incomplete data, we first conduct bootstrap imputation through two
steps: (i) generate $B$ bootstrap datasets of the original data; (ii) conduct
a single imputation for each bootstrap dataset using an imputation method of
choice. Following Hu et al., (hu2021variable, ) in our simulations and case
study, we used $B=100$ and performed standard imputation program using the R
package mice. Variable selection via each of the methods considered was
performed for all $B$ bootstrap imputed datasets, and the final set of
predictors were selected if they were selected in at least $\pi B$ imputed
datasets, where $\pi\in(0,1)$ is a fraction threshold for selecting a
predictor.

##  4 Simulation

###  4.1 Treatment effect heterogeneity

We carried out a contextualized simulation to evaluate the performance of our
proposed method in estimating the heterogeneous treatment effect using
observational data with clustered survival outcomes and multiple treatments.
For methods comparison, we adapted the popularly used inverse probability
weighting method into the setting of clustered and censored survival data to
form two comparison methods: inverse probability of treatment weighting with
the random-intercept Cox regression model (IPW-riCox) and doubly robust
random-intercept additive hazards model (DR-riAH).(cai2011additive, ) In
addition, we considered the PEAMM (bender2018generalized, ) – another outcome
modeling based method flexible at capturing nonlinear relationships. Because
not all methods directly model the survival times $T_{ik}$, to objectively
compare methods, we use the individual-specific counterfactual survival curve
as the basis. That is, the ISTE $\zeta_{a_{j},a_{j^{\prime}}}(\bm{x}_{ik})$
will now be defined in a similar way as in equation (3) but using functionals
of the counterfactual survival curve – either the survival probability at a
fixed point in time or the conditional restricted mean survival time (RMST).
(royston2013restricted, ) We will compare each method’s ability to accurately
estimate the ISTE $\zeta_{a_{j},a_{j^{\prime}}}(\bm{x}_{ik})$ and the average
subgroup treatment effect $M(S^{h}_{a_{j},a_{j^{\prime}}})$. The performance
metrics used to evaluate the methods were the relative bias (or percent bias),
root-mean-squared-error (RMSE) and precision in the estimation of
heterogeneous effects (PEHE). (hu2021estimating, ) The relative bias evaluates
the bias in the estimation of treatment effect on the relative scale with
respect to the true treatment effect (which is constant across different
methods). The PEHE measures the difference between the true and estimated
survival probabilities or RMST across all data points. A smaller value of PEHE
indicates a higher estimation accuracy.

For weighting based methods IPW-riCox and DR-riAH, we used Super Leaner
(van2007super, ) to estimate the stabilized inverse probability of treatment
weights. Super Learner was implemented using the R package SuperLearner and
the library argument was set to SL.library = c("SL.xgboost",
"SL.bartMachine","SL.gbm"). We used R package coxme to implement IPW-riCox and
timereg to implement DR-riAH. For PEAMM, the gam function from R package mgcv
and two helper functions, (as_ped and add_surv_prob), from R package pammtools
were used. For all methods, all confounding variables available to the analyst
in their original forms were included in the corresponding models.

####  4.1.1 Simulation design

We simulated datasets with structures similar to the COVID-19 data used in our
case study (Section 5). We generated $K=10$ clusters, each with a sample size
of $n_{k}=200$. The total sample size is $N=2000$. We simulated 7 confounding
variables for each individual $i$ in cluster $k$: five continuous variables
$X_{ikl}\sim N(0,1),l\in\\{1,\ldots,5\\}$; and two categorical variables
$X_{ikl}\sim Multinomial(1,.3,.3,.4),l=6,7$. We generated $J=3$ treatment
groups with unequal sample sizes; the ratio of individuals across treatment
groups was 5:3:1, which is similar to the ratio of individuals across the
three treatment groups in our case study. The treatment assignment mechanism
follows a random-intercept multinomial logistic regression model,

$\begin{split}\ln\dfrac{P(A_{ik}=1)}{P(A_{ik}=3)}&=1.5+.1X_{ik1}+.1X_{ik2}+.1X_{ik3}+.5X_{ik4}+.4X_{ik5}+.2X_{ik6}+.3X_{ik7}+.4X_{ik2}^{2}+.4X_{ik2}^{2}X_{ik5}+\tau_{k}\\\
\ln\dfrac{P(A_{ik}=2)}{P(A_{ik}=3)}&=.7+.1X_{ik1}+.3X_{ik2}+.2X_{ik3}+.2X_{ik4}+.1X_{ik5}+.4X_{ik6}+.5X_{ik7}-.3X_{ik2}X_{ik4}+.7X_{ik2}^{2}X_{ik4}+\tau_{k},\end{split}$

where the random intercept $\tau_{k}\sim N(0,1^{2})$. The parameter setup
induces a moderate level of confounding, leading to moderate covariate
overlap. The sparsity of overlap can be visualized by the distributions of
true generalized propensity scores, shown in Web Figure 1. This moderate
covariate overlap mimics the overlap level observed in the COVID-19 dataset
(Web Figure xx). Following strategies (hu2021estimating, ; lu2018estimating, )
recommended for assessing the performance of methods in estimating
heterogeneous treatment effects, we subclassified the simulated individuals
into 40 subgroups based on the distribution of true generalized propensity
scores, (feng2012generalized, ) and calculated the biases and RMSE for each
subgroup across 250 data replications. Note that the generalized propensity
scores for each individual sum to one,
$P(A_{ik}=1\mid\bm{X}_{ik},\tau_{k})+P(A_{ik}=2\mid\bm{X}_{ik},\tau_{k})+P(A_{ik}=3\mid\bm{X}_{ik},\tau_{k})=1$.
The subclassification is based on intervals of the true generalized propensity
scores for treatment 1 and treatment 2, representing a full range of
assignment propensity to each treatment group. Subclasses of generalized
propensity scores are presented in Web Table 1.

We simulated three sets of true counterfactual survival times from the Weibull
survival distribution, which can be parameterised as both the accelerated
failure time model and Cox proportional hazards model, (hu2022flexible, )

$\displaystyle
S_{ik}(t)=\exp\left[-\left\\{\lambda_{a_{j}}\exp\left(q_{a_{j}}\left(\bm{X}_{ik},b_{k}\right)\right)t\right\\}^{\eta}\right],$

where $\lambda_{a_{j}}$ is a treatment group specific parameter, $\eta$ is the
shape parameter and $q_{a_{j}}\left(\bm{X}_{ik},b_{k}\right)$ represents a
general functional form of covariates and the random-intercept. Using the
inverse transform sampling, we generated counterfactual survival times by

$\displaystyle T_{ik}\left(a_{j}\right)=\left\\{\frac{-\log
U}{\lambda_{a_{j}}\exp\left(q_{a_{j}}\left(\bm{X}_{ik},b_{k}\right)\right)}\right\\}^{1/\eta}$
(7)

for $a_{j}\in\\{1,2,3\\}$, where $U$ is a random variable following the
uniform distribution on the interval $[0,1]$ and
$\lambda_{a_{j}}=\\{5000,800,1200\\}$ for $a_{j}=1,2,3$. Observed and
uncensored survival times are generated as

$T_{ik}=\sum\limits_{a_{j}\in\\{1,2,3\\}}T_{ik}(a_{j})I(A_{ik}=a_{j}).$

We considered three treatment effect heterogeneity settings with an increasing
level of complexity:

  1. (a)

$q_{1}\left(\bm{X}_{ik},b_{k}\right)=.1X_{ik1}+.3X_{ik2}+\sin(\pi
X_{ik3})+.6X_{ik4}+.5X_{ik5}+1.2X_{ik6}+.4X_{ik7}+.3X_{ik2}^{2}+.5X_{ik4}X_{ik5}+b_{k}-1$

$q_{2}\left(\bm{X}_{ik},b_{k}\right)=.4X_{ik1}+1.2\sin(\pi
X_{ik2})+.4X_{ik3}+.3X_{ik4}+1.0X_{ik5}+.8X_{ik6}+.2X_{ik7}+.7X_{ik1}^{2}+.4X_{ik1}X_{ik4}+b_{k}$

$q_{3}\left(\bm{X}_{ik},b_{k}\right)=.4X_{ik1}+.9X_{ik2}+.4X_{ik3}+.9X_{ik4}+.4X_{ik5}+.4X_{ik6}+.3X_{ik7}+b_{k}-2$

  2. (b)

$q_{1}\left(\bm{X}_{ik},b_{k}\right)=.1X_{ik1}+.3X_{ik2}+\sin(\pi
X_{ik3})+.6X_{ik4}+.5X_{ik5}+1.2X_{ik6}+.3X_{ik7}+.3X_{ik2}^{2}+.5X_{ik4}X_{ik5}+b_{k}-1$

$q_{2}\left(\bm{X}_{ik},b_{k}\right)=.4X_{ik1}+1.2\sin(\pi
X_{ik2})+.4X_{ik3}+.3X_{ik4}+1.0X_{ik5}+.8X_{ik6}+.1X_{ik7}+.7X_{ik1}^{2}+.4X_{ik1}X_{ik4}+b_{k}$

$q_{3}\left(\bm{X}_{ik},b_{k}\right)=.4\sin(\pi
X_{ik1})+.9X_{ik2}+.9X_{ik3}+.4X_{ik4}+.4X_{ik5}+.9X_{ik6}+.3X_{ik7}+.4X_{ik4}^{2}-.3X_{ik2}X_{ik3}+b_{k}-3$

  3. (c)

$q_{1}\left(\bm{X}_{ik},b_{k}\right)=.1X_{ik1}+.3X_{ik2}+\sin(\pi
X_{ik3})+.6X_{ik4}+.5X_{ik5}+1.2X_{ik6}+.3X_{ik2}^{2}+.5X_{ik4}X_{ik5}+b_{k}-1$

$q_{2}\left(\bm{X}_{ik},b_{k}\right)=.4X_{ik1}+1.2\sin(\pi
X_{ik3})+.4X_{ik4}+.3X_{ik5}+1.0X_{ik6}+.8X_{ik7}+.7X_{ik1}^{2}+.4X_{ik1}X_{ik4}+b_{k}$

$q_{3}\left(\bm{X}_{ik},b_{k}\right)=.4\sin(\pi
X_{ik2})+.9X_{ik3}+.9X_{ik4}+.4X_{ik5}+.4X_{ik6}+.9X_{ik7}+.4X_{ik4}^{2}-.3X_{ik2}X_{ik3}+b_{k}-3$,

where the random intercept $b_{k}\sim N(0,4^{2})$. Across all three
heterogeneity settings, $X_{ik3},X_{ik4},X_{ik5},X_{ik6}$ were confounders
related to both the treatment assignment mechanism and the potential outcome
generating process. In scenario (a), nonlinear covariate effects are present
in treatment groups 1 and 2. In scenario (b), there are nonlinear covariate
effects in all three treatment groups. Scenario (c) is similar to scenario (b)
except that a nonoverlapping set of covariates are included in three treatment
groups. This additional complexity represents an additional source of
treatment effect heterogeneity. The parameter $\eta$ was set to 2 and
$\exp(0.7+0.5x_{ik1})$ to respectively induce proportional hazards and
nonproportional hazards. We further generate the censoring time $C$
independently from an exponential distribution with the rate parameter chosen
to induce a censoring proportion of 50%, which is close to the censoring
proportion observed in the COVID-19 data in our case study. The Kaplan-Meier
survival curves by treatment group for all three heterogeneity settings under
both proportional and nonprotional hazards are presented in Web Figure 2.

####  4.1.2 Simulation results

Table 1 summarizes the overall precision measure PEHE (based on 3-week
survival probability) for each of four methods under three heterogeneity
settings and both proportional and nonproportional hazards. First, flexible
methods riAFT-BART and PEAMM and doubly robust weighting method DR-riAH
yielded substantially smaller PEHE than the parametric single robust method
IPW-riCox across all simulation scenarios. Second, as the heterogeneity
setting became more complex, flexible methods riAFT-BART and PEAMM tended to
have better performance evidenced by decreasing mean and variation of the
PEHE. All methods show decreasing precision in estimating ISTE, demonstrated
by increasing PEHE value, when the proportional hazards assumption fails to
hold. This is possibly due to more complex data generating processes, in which
the rate parameter $\eta$ depends on covariate $x_{ik1}$. The decrease in
estimation precision is the smallest for riAFT-BART. Overall, across all
scenarios, the proposed method riAFT-BART had the smallest PEHE, or the
highest accuracy in estimating the ISTE. Similar observations can be made from
the PEHE measures based on 3-week RMST, which are provided in Web Table 2.

Figure 1 displays the relative biases and RMSE results among 40 generalized
propensity score subgroups, each averaged across 250 simulation runs. Three
pairwise treatment effects were estimated by averaging the ISTE
$\hat{\zeta}_{a_{j},a_{j^{\prime}}}(\bm{x}_{ik})$ (based on 3-week survival
probability) across individuals in each subgroup. Under both settings of
proportional hazards and nonproportional hazards, the proposed method riAFT-
BART boasts the smallest biases and RMSE across all three treatment effects,
particularly for subpopulation near the centroid of the generalized propensity
score distribution. In sharp contrast, IPW-riCox, which encodes parametric and
linearity assumptions and assumes proportional hazards, yielded the largest
biases and RMSE. The increasing complexity level of treatment effect
heterogeneity did not impact the performance of more flexible methods PEAMM
and riAFT-BART, but induced larger biases and RMSE for parametric methods IPW-
riCox and DR-riAH. The violation of the proportional hazards assumption, which
introduced an additional data complexity (the rate parameter $\eta$ depends on
the covariate $x_{ik1}$), had a negligible effect on the performance of riAFT-
BART, but all three comparison methods had a deteriorated performance. It is
worth noting that all methods produced large biases and RMSE for subpopulation
at the tail regions of the generalized propensity score distribution,
suggesting that causal inferences drawn about this subpopulation may not be
reliable. Results based on 3-week RMST, shown in Web Figure 3, convey the same
messages.

Figure 1: Relative biases (Panel A) and root-mean-squared-errors (RMSE) (Panel B) among 40 generalized propensity score subgroups under 6 data configurations: (heterogeneity settings a, b, c) × (proportional hazards (PH) and nonproportional hazards (nPH)) for each of four methods, IPW-riCox, DR-riAH, PEAMM and riAFT-BART. Three pairwise treatment effects were estimated by averaging the individual survival treatment effect (based on 3-week survival probability) across individuals in each subgroup. Each boxplot visualizes the distribution of relative biases or the distribution of RMSE for 40 subgroups, each averaged across 250 simulation runs. Table 1: Mean and (standard deviation) of precision in the estimation of heterogeneity effects (PEHE) across 250 data replications for each of the 4 methods based on 3-week survival probability under six configurations: (proportional hazards [PH] vs. nonproporitonal hazards [nPH]) $\times$ (heterogeneity setting [HS] (a) vs. (b) vs. (c)). |  | PH |  | nPH
---|---|---|---|---
| Methods | HS(a) | HS(b) | HS(c) |  | HS(a) | HS(b) | HS(c)
Trt 1 vs. 2 | IPW-riCox | .061 (.028) | .067 (.029) | .070 (.030) |  | .082 (.028) | .089 (.029) | .092 (.030)
DR-riAH | .018 (.021) | .020 (.022) | .026 (.023) |  | .028 (.020) | .033 (.021) | .035 (.021)
PEAMM | .028 (.023) | .024 (.022) | .020 (.021) |  | .038 (.024) | .034 (.023) | .031 (.023)
riAFT-BART | .012 (.010) | .009 (.009) | .006 (.009) |  | .017 (.012) | .014 (.011) | .011 (.011)
Trt 1 vs. 3 | IPW-riCox | .065 (.027) | .070 (.028) | .073 (.031) |  | .083 (.029) | .090 (.030) | .093 (.031)
DR-riAH | .021 (.018) | .026 (.020) | .029 (.022) |  | .031 (.019) | .036 (.020) | .039 (.022)
PEAMM | .035 (.024) | .030 (.023) | .026 (.022) |  | .049 (.025) | .044 (.024) | .037 (.023)
riAFT-BART | .016 (.010) | .013 (.009) | .010 (.009) |  | .022 (.014) | .019 (.013) | .015 (.012)
Trt 2 vs. 3 | IPW-riCox | .088 (.030) | .095 (.031) | .098 (.032) |  | .105 (.031) | .115 (.032) | .118 (.033)
DR-riAH | .038 (.023) | 043 (.023) | .045 (.024) |  | .045 (.021) | .050 (.022) | .052 (.022)
PEAMM | .049 (.023) | .045 (.022) | .041 (.022) |  | .058 (.024) | .054 (.023) | .049 (.023)
riAFT-BART | .017 (.012) | .015 (.011) | .013 (.011) |  | .022 (.015) | .020 (.014) | .018 (.014)

###  4.2 Variable selection

We adopt the simulation settings used in Hu et al. hu2021variable  to evaluate
the performance of variable selection methods among incomplete, clustered
survival data. We used the multivariate amputation approach
schouten2018generating  to generate missing data scenarios with desired
missingness percentages under the missing at random mechanism. The
multivariate amputation proceeds through the following three steps: (i)
randomly divide the complete data into a certain number of subsets, each
allowing the specification of any missing data pattern; (ii) calculate the
weighted sum scores for individuals in each subset to amputate the data. The
weighted sum score relates the variable to be amputated to the variables on
which the missingness in the amputated variable depends; (iii) apply a
logistic distribution function van2018flexible  on the weighted sum scores to
compute the missingness probability, which is used to determine whether the
data point becomes missing or not. We investigated the practical operating
characteristics of our proposed variable selection method using riAFT-BART,
and compare it to PEAMM, frailty models with penalized likelihood
(FrailtyPenal), frality models with penalized profile h-likelihood (FrailtyHL)
and backward stepwise selection with random-intercept Cox regression model
(riCox). Four performance metrics were used to evaluate the variable selection
results: precision, recall, $F_{1}$ and Type I error (Section 3.3).

The multivariate amputation procedure was implemented via the ampute function
of the R package mice. (schouten2018generating, ) Imputation was performed
using the R package mice. (buuren2010mice, ) To implement our proposed method
using riAFT-BART, we used 4500 posterior draws with the first 1000 discarded
as burn-in. The permutation distributions of VIP scores were constructed from
100 riAFT-BART model runs. To implement the variable selection procedure based
on PEAMM, we used the helper function (as_ped) from R package pammtools to
reformat the data and used the gam function from R package mgcv for variable
selection. Variable selection by frailty models was implemented using the
frailtyPenal function of R package frailtypack for penalized marginal
likelihood; and using the function frailty.vs from R package frailtyHL for
penalized profile h-likelihood. For backward stepwise selection, the random-
intercept Cox regression model was fitted using the coxme function from R
package coxme.

####  4.2.1 Simulation design

The simulation scenarios are motivated by the data structures observed in the
COVID-19 data set used in our case study. We considered $K=10$ clusters, each
with a sample size of $n_{k}=200$. We generated 8 _useful_ predictors that are
truly related to the survival outcomes, $X_{ik1},\ldots,X_{ik8}$, and 20 noise
predictors, $X_{ik9},\ldots,X_{ik28}$. We simulated $X_{ik1}$ and $X_{ik2}$
independently from $\text{Bern}(0.5)$, $X_{ik3}$, $X_{ik4}$ from the standard
normal distribution $N(0,1)$, and $X_{ik5},X_{ik6},X_{ik7},X_{ik8}$ were
designed to have missing values under the missing at random mechanism. For a
missing at random predictor variable, we specify the true forms in which the
predictor depends on the other predictors as the following:

$\displaystyle x_{ik5}{\,|\,}x_{ik2},x_{ik3}\sim N(0.3x_{ik2}-0.2x_{ik3},1)$
$\displaystyle x_{ik6}{\,|\,}x_{ik3},x_{ik4}\sim
N(-0.4x_{ik3}+0.4x_{ik4}+0.3x_{ik3}x_{ik4},1)$ $\displaystyle
x_{ik7}{\,|\,}x_{ik4},x_{ik5},x_{ik6}\sim
N(0.1x_{ik4}(x_{ik5}-2)^{2}-0.1x_{ik6}^{2},1)$ $\displaystyle
x_{ik8}{\,|\,}x_{ik5},x_{ik6},x_{ik7}\sim
N(-0.3x_{ik5}^{2}+0.5x_{ik6}+0.3x_{ik7}+0.2x_{ik6}x_{ik7},1).$

Among the 20 noise predictors, 10 were generated as continuous variables
$X_{ik9},\ldots,X_{ik18}\stackrel{{\scriptstyle i.i.d}}{{\sim}}N(0,1)$, and
the other 10 were generated as binary variables
$X_{ik19},\ldots,X_{ik28}\stackrel{{\scriptstyle
i.i.d}}{{\sim}}\text{Bern}(0.5)$.

Similar to the outcome generating process described in Section 4.1.1, we
generate the observed and uncensored survival times from the Weibull survival
distribution,

$\displaystyle T_{ik}=\left[\frac{-\log
U}{\lambda\exp\left\\{q(\bm{X}_{ik},b_{k})\right\\}}\right]^{1/\eta}$

where $\lambda=3000$, $U$ is a random variable following the uniform
distribution on the interval $[0,1]$, $b_{k}\sim N(0,4^{2})$, and

$\displaystyle q(\bm{X}_{ik},b_{k})$ $\displaystyle=$ $\displaystyle
1.8x_{ik1}+0.5x_{ik2}+1.1x_{ik3}-0.4e^{x_{ik5}}+0.4(x_{ik6}-1.5)^{2}+0.1(x_{ik7}-0.1)^{3}$
$\displaystyle-5\sin(0.1\pi x_{ik4}x_{ik8})-0.4x_{ik5}x_{ik7}+b_{k}$

We further generated the censoring time $C$ independently from an exponential
distribution with the rate parameter selected to induce the censoring
proportion of 50%. The parameter $\eta$ was set to 2 and
$\exp{(0.7+0.5x_{ik1})}$ to respectively produce proportional hazards and
nonproportional hazards. This data generating process was designed to make it
difficult for any method to accurately model the survival outcome. We
considered the outcome model with arbitrary data complexity that reflects
common situations in health datasets: (i) discrete predictors with strong
($X_{ik1}$) and moderate ($X_{ik2}$) associations; (ii) both linear and
nonlinear forms of continuous predictors ($X_{ik6}$ and $X_{ik7}$); (iii)
nonadditive effects ($X_{ik4}X_{ik8}$).

After generating the fully observed data, we amputated the four predictors
designed to have missing values $X_{ik5},X_{ik6},X_{ik7},X_{ik8}$ under the
missing at random mechanism using the multivariate amputation approach.
Additional technical details for the amputation procedure are available in
Schouten et al.; (schouten2018generating, ) and additional amputation details
for the simulation appear in Web Section 2. The resultant simulation data have
an overall missingness proportion of 40%.

####  4.2.2 Simulation results

Table 2 summarizes variable selection results by each method. The optimal
performance among incomplete data is shown for the best selection threshold
value of $\pi$ – chosen on the basis of the $F_{1}$ score (see Web Figures
4-5); and is compared against the performance on fully observed data (before
amputation) and among complete cases only. The proposed method riAFT-BART has
the best performance across all three data panels (fully observed, with
missing data, complete cases), evidenced by the highest precision, recall and
$F_{1}$ score, and the lowest type I error. The second best performing method
is PEAMM, for its modeling flexibility over parametric methods. In the
presence of missing covariate data, all methods had a slightly deteriorated
performance when performing variable selection in combination with bootstrap
imputation, but the performance of variable selection among only complete
cases was considerably lower.

Table 2: Simulation results for each variable selection approach performed on the fully observed data, and among incomplete data with 40% overall missingness in covariates under proportional hazards and non-proportional hazards. Imputation was conducted via mice. For each of five methods, we show results corresponding to the best threshold values of $\pi$ (based on $F_{1}$). | Proportional hazards |  |  | Non-proportional hazards
---|---|---|---|---
| Precision | Recall | $F_{1}$ | Type I error |  |  | Precision | Recall | $F_{1}$ | Type I error
Fully observed data
riAFT-BART | 0.95 | 0.88 | 0.91 | 0.01 |  | riAFT-BART | 0.93 | 0.87 | 0.89 | 0.02
PEAMM | 0.83 | 0.82 | 0.82 | 0.02 |  | PEAMM | 0.81 | 0.80 | 0.80 | 0.03
FrailtyHL | 0.81 | 0.78 | 0.79 | 0.03 |  | FrailtyHL | 0.78 | 0.71 | 0.74 | 0.04
FrailtyPenal | 0.79 | 0.76 | 0.77 | 0.04 |  | FrailtyPenal | 0.75 | 0.69 | 0.72 | 0.05
riCox | 0.78 | 0.72 | 0.75 | 0.05 |  | riCox | 0.72 | 0.67 | 0.69 | 0.05
40% overall missing
riAFT-BART $\pi=0.1$ | 0.93 | 0.85 | 0.89 | 0.02 |  | riAFT-BART $\pi=0.1$ | 0.89 | 0.84 | 0.86 | 0.03
PEAMM $\pi=0.7$ | 0.88 | 0.74 | 0.80 | 0.03 |  | PEAMM $\pi=0.6$ | 0.84 | 0.72 | 0.78 | 0.04
FrailtyHL $\pi=0.6$ | 0.79 | 0.73 | 0.76 | 0.05 |  | FrailtyHL $\pi=0.7$ | 0.75 | 0.70 | 0.71 | 0.04
FrailtyPenal $\pi=0.6$ | 0.76 | 0.71 | 0.74 | 0.04 |  | FrailtyPenal $\pi=0.7$ | 0.72 | 0.67 | 0.69 | 0.05
riCox $\pi=0.8$ | 0.85 | 0.67 | 0.73 | 0.04 |  | riCox $\pi=0.7$ | 0.71 | 0.63 | 0.68 | 0.07
Complete cases
riAFT-BART | 0.85 | 0.81 | 0.83 | 0.03 |  | riAFT-BART | 0.83 | 0.79 | 0.81 | 0.04
PEAMM | 0.77 | 0.73 | 0.75 | 0.02 |  | PEAMM | 0.73 | 0.71 | 0.72 | 0.03
FrailtyHL | 0.73 | 0.68 | 0.71 | 0.04 |  | FrailtyHL | 0.70 | 0.65 | 0.67 | 0.05
FrailtyPenal | 0.71 | 0.65 | 0.69 | 0.05 |  | FrailtyPenal | 0.69 | 0.62 | 0.65 | 0.04
riCox | 0.71 | 0.66 | 0.68 | 0.04 |  | riCox | 0.66 | 0.62 | 0.64 | 0.03
Figure 2: The distribution, across 250 data replications, of the numbers of
selected noise predictors and useful predictors for each of five methods:
riAFT-BART, PEAMM, FrailtyHL, FrailtyPenal and riCox, with clustered survival
data generated under both proportional hazards (PH) and non-proportional
hazards (nPH). The total number of useful predictors is 8 and the total number
of noise predictors is 20. There are $K=10$ clusters, each with a size of 200;
the total sample size is 2000. The overall proportion of missingness is 40%.

Demonstrated by Figure 2, the proposed method riAFT-BART tended to select the
most useful predictors and the least noise predictors, followed by another
flexible method PEAMM, under both proportional hazards and non-proportional
hazards. Among parametric methods, penalized likelihood methods for frailty
models had slightly better performance than riCox. All methods tended to
select more noise predictors and less useful predictors under non-proportional
hazards; the decrease in the performance due to elevated data complexity
associated with nonproportional hazards, was only marginal for riAFT-BART but
considerable for less flexible methods, particularly FrailtyHL and riCox.

Figure 3: Power of each of five methods: riAFT-BART, PEAMM, FrailtyHL,
FrailtyPenal and riCox, for selecting each of 8 useful predictors with
clustered survival data generated under proportional hazards (PH) and non-
proportional hazards (nPH), based on 250 data replications. There are $K=10$
clusters, each with a size of 200; the total sample size is 2000. The overall
proportion of missingness is 40%. Filled symbols represent the PH setting, and
open symbols correspond to the nPH setting.

A perusal of each method’s ability to select each of useful predictors is
visualized in Figure 3. All methods had a good power ($>0.8$) for selecting a
strong discrete predictor ($X_{ik1}$), but unsatisfactory power ($<0.5$) for
detecting a discrete predictor ($X_{ik2}$) that only has a moderate
relationship with the outcome. For continuous variables, all methods performed
well in selecting $X_{ik3}$, which has only a linear relationship with the log
survival time, but riAFT-BART produced a significantly higher power for
selecting the predictor $X_{ik4}$, which has a complex nonadditive form in
model (4.2.1). Among the four predictors with missing data
$X_{ik4},\ldots,X_{ik8}$, the advantage of riAFT-BART is demonstrated by the
higher power of selecting predictors that have more complex functional forms
that are difficult to be captured by other methods considered. A direct
comparison of the proportional hazards setting with nonproportional hazards
shows that all methods performed better under proportional hazards, but the
impact of the violation of proportional hazards is the smallest for riAFT-
BART.

As our proposed method riAFT-BART pivots on the predictor’s VIP, we examined
the convergence of the Markov chain by plotting 4500 posterior draws, with the
first 1000 discarded as burn-in, of the VIP for four useful predictors and
four noise predictors in Web Figure 6. The method converged well.

##  5 Case study: Predictors for in-hospital mortality and heterogeneous
survival treatment effects among post-ICU COVID patients

In this case study, we first applied the proposed variable selection methods
to select important predictors for in-hospital mortality using a comprehensive
COVID-19 dataset, and then evaluated the treatment effect heterogeneity on
patient survival using the selected predictors by implementing our proposed
method riAFT-BART. The COVID-19 dataset was drawn from six hospitals
(clusters) of the Epic electronic medical records system at the Mount Sinai
Medical Center. Included in this dataset were 1955 patients who were diagnosed
with COVID-19 between March 10, 2020 to February 26, 2021 and were in
intensive care (ICU) due to COVID-19 during their hospital stay. (hu2021covid,
) Post ICU admission, each patient received COVID-specific medications: either
dexamethasone, or remdesivir, or the combination of dexamethasone and
remdesivir (dexamethasone+remdesivir). The cluster sizes and number of
patients in each treatment group are provided in Web Table 4.

We used the following 28 baseline (time of ICU admission) variables: age, sex,
self-reported race and ethnicity, smoking status, binary comorbidities
including hypertension, coronary artery disease, cancer, diabetes, asthma, and
chronic obstructive pulmonary disease, vital signs (measured within 24 hours
of ICU admission) including temperature, systolic blood pressure, diastolic
blood pressure, patient oxygen level (see Web Table 3 for definition), heart
rate, the fraction of inspired oxygen, body mass index, oxygen saturation,
risk score (sequential organ failure assessment score and glasgow coma scale),
laboratory test results including partial pressure of oxygen, D dimer value
and inflammatory markers such as lactate dehydrogenase, ferritin, c-reactive
protein, creatinine and white blood cell count. In cases where multiple vital
signs/risk score/lab results were recorded during the first 24 hours of ICU
admission, we used the mean value of the multiple measurements. Among the 28
candidate predictor variables, 17 variables have missing data and the
missingness proportion ranges from 1.4% to 28.6%. Only 1187 (60.7%) patients
have fully observed data for all covariates, which motivated the overall
missingness proportion of 40% used in our simulation study (Section 4.2).
Detailed summary statistics for baseline variables appear in Web Table 4.

The Kaplan-Meier survival curves by treatment group are displayed in Web
Figure 7, with crossovers suggesting nonproportional hazards. Assessed by the
distribution of the generalized propensity scores shown in Web Figure 8, there
is moderate covariate overlap across three treatment groups in the COVID-19
data.

Table 3: Variable selection results by each of 5 methods, with the best
selection threshold value of $\pi$ suggested in simulations under non-
proportional hazards. riAFT-BART and riCox both selected 8 variables, PEAMM
selected 7 variables, FrailtyPenal selected 6 variables and FrailtyHL selected
9 variables.

Variables | riAFT-BART | PEAMM | FrailtyPenal | FrailtyHL | riCox
---|---|---|---|---|---
Oxygen level | Yes | Yes | No | Yes | Yes
Age | Yes | Yes | Yes | Yes | No
Creatinine | Yes | Yes | No | Yes | Yes
White blood cell | Yes | Yes | No | No | Yes
SOFA score | Yes | Yes | Yes | Yes | No
Race | Yes | No | Yes | Yes | Yes
D dimer | Yes | No | No | Yes | No
LDH | Yes | Yes | Yes | Yes | Yes
C-reactive protein | No | Yes | No | Yes | Yes
Ferritin | No | No | No | Yes | Yes
Oxygen saturation | No | No | Yes | No | No
Glasgow coma scale | No | No | Yes | No | Yes

Abbreviations: LDH = Lactate dehydrogenase; SOFA = Sequential organ failure
assessment

Figure 4: The distribution of cross-validated concordance statistics across
250 data replications for each of five methods using the COVID-19 dataset.

We first implemented all five variable selection methods. All 28 candidate
predictors were included in the imputation models, and the optimal threshold
value of $\pi$ for each method, which produced the best $F_{1}$ score in the
simulation study, was used to select the final set of predictors. As
summarized in Table 3, FrailtyHL selected the most (9) predictors, riAFT-BART
and riCox both selected 8 predictors, PEAMM selected 7 preditors and
FrailtyPenal selected the least (6) predictors. Age and lactate dehydrogenase
were the only two predictors selected by all five methods. Variables selected
by the two flexible methods riAFT-BART and PEAMM were largely in agreement
with each other, with the exception that riAFT-BART but not PEAMM selected
Race and D dimer, and PEAMM but not riAFT-BART selected C-reactive protein.

To further distinguish between methods based on their ability to select the
most relevant predictors, we computed the predictive performance of the five
methods that modeled the outcome variable on their respectively selected
predictor variables. For survival outcomes, we used the concordance statistic
as a metric to assess the predictive performance. (harrell1996multivariable, )
The concordance statistic is the fraction of concordant pairs among all pairs
of individuals. A pair of patients is called concordant if the patient with
lower (higher) predicted survival probability at a given time point had a
longer (shorter) time to death. In this case study, we chose 21 days post ICU
admission as the target time point, because it is a sufficiently long time to
evaluate the effects of COVID-19 treatments on in-hospital mortality for
severe patients. As the COVID-19 dataset is subject to missing covariate data,
to calculate the cross validated concordance statistic, we follow the strategy
described in Lin et al. (lin2022a, ) We first split the data into two halves
and implemented each of the five variable selection methods on half of the
data. Then we imputed the other half of the data with a single imputation and
recorded the concordance statistic. Finally we repeated the previous two steps
250 times to get the distribution of the cross-validated concordance
statistics. Figure 4 shows that the proposed method riAFT-BART yielded the
highest predictive performance with a median concordance statistic $>0.7$.

Supported by both the simulation study and case study, riAFT-BART provided the
best variable selection results. Thus we used the 8 predictors selected by
riAFT-BART as confounding variables, and estimated the heterogeneous treatment
effects on patient survival using methodology described in Section 2. To
address the missing data issue when drawing causal inferences about treatment
effects, we fitted an riAFT-BART model (LABEL:eq:_riAFT-BART) on each of the
100 bootstrap and imputed datasets (previously conducted for variable
selection), and drew statistical inferences about the ISTE
$\zeta_{a_{j},a_{j^{\prime}}}(\bm{x}_{ik})$, defined in equation (3), by
pooling posterior samples of the ISTE,
$\tilde{f}^{d}(a_{j},\bm{x}_{ik})-\tilde{f}^{d}(a_{j^{\prime}},\bm{x}_{ik})$
as in equation (2), across model fits arising from the multiple datasets.
(zhou2010note, ; hu2022aflexible, ) Inferences about the population average
treatment effects $\theta_{a_{j},a_{j^{\prime}}}$ can be obtained from the
posterior samples of
$\tilde{f}^{d}(a_{j},\bm{x}_{ik})-\tilde{f}^{d}(a_{j^{\prime}},\bm{x}_{ik})$
across the sample population. (hu2022flexible, ) The population average
treatment effects, summarized in Web Table 5, suggest that there was no
statistically significant difference in treatment effect when comparing
dexamethasone with remdesivir or with a combined use of dexamethasone and
remdesivir. However, there was a statistically significant treatment benefit
associated with dexamethasone+remdesivir compared to using remdesivir alone.
Web Figure 9 demonstrates that substantial variability in the institutional
effect was captured by the riAFT-BART model: the Mount Sinai main hospital had
considerably better outcomes than the Mount Sinai Queens hospital.

Figure 5: Final Random Forests model fit to the posterior mean of the
individual survival treatment effect comparing remdesivir and dexamethasone +
remdesivir. Values in each node correspond to the posterior mean, in terms of
difference in log survival days, for the subgroup of individuals represented
in that node. Uncertainty intervals were obtained by pooling the posterior
samples arising from the multiple imputed data sets. WBC: White blood cell.

Turning to the treatment effect heterogeneity. Figure 5 summarizes the
combination rules of covariates $\mathcal{S}_{a_{j},a_{j^{\prime}}}$ that have
differential subgroup treatment effect
$\mathcal{M}(\mathcal{S}_{a_{j},a_{j^{\prime}}})$, defined in equation (4),
than the sample population treatment effect. When comparing remdesivir and
dexamethasone + remdesivir across the sample population, there was significant
treatment benefit, suggested by a longer expected survival time, for the
combined use of dexamethasone and remdesivir. Our analysis suggests that
combining dexamethasone and remdesivir offered enhanced treatment benefit for
relatively healthier COVID-19 patients with lower oxygen saturation
($<95.5\%$) and a normal white blood cell count ($<11.4K/\mu L$), demonstrated
by the leftmost branch of the tree diagram in Figure 5. By contrast, the
rightmost branch suggests that there was no statistically significant
treatment benefit associated with either treatment choice for comparatively
unhealthier patients with a higher level of oxygen saturation ($>95.5\%$).

##  6 Discussion

We describe two important and practical utilities of a new Bayesian machine
learning tool riAFT-BART for analyzing clustered and censored survival data.
To address the implications of complex data structures of large-scale clinical
datasets generated from multiple institutions for causal inferences about
population treatment effect, riAFT-BART was developed to accurately estimate
the population average treatment effect on patient survival while accounting
for the variation in institutional effects. (hu2022flexible, ) Moving from
broad applicability of a population study to personalized medicine requires
the understanding of treatment effect heterogeneity. The treatment effect
heterogeneity was traditionally identified by first enumerating potential
effect modifiers with subject-matter experts and then estimating the average
treatment effect within each subgroup.kent2020predictive  This approach is
particularly suitable to confirmatory treatment effect heterogeneity analysis
in randomized clinical trials. In observational data with a large number of
pre-treatment covariates, such _a priori_ specification that separates the
issues of confounding and treatment effect heterogeneity is often practically
infeasible. We described a way in which the new tool riAFT-BART can be used to
perform an exploratory treatment effect heterogeneity analysis to generate
scientifically meaningful hypotheses and treatment effect discovery. Through a
comprehensive simulation representative of complex confounding and
heterogeneity settings with right-censored survival data, we provided new
empirical evidence for the better performance of our proposed method in
comparison with three existing methods. In the application to a large COVID-19
dataset, we found that oxygen saturation and white blood cell count were two
key patient factors that modulated the comparative treatment effect between
remdesivir and dexamethasone + remdesivir. We further exploited the posterior
samples of the ISTE from riAFT-BART and the random forests model to explore
clinically meaningful relationships between COVID-19 medications, patient
profile and expected survival time. The results could facilitate treatment
effect discovery in subpopulations and may inform personalized treatment
strategy and the planning of future confirmatory randomized trials.

Drawing on prior work (hu2021variable, ) that combines bootstrap imputation
and BART for variable selection among incomplete data, we develop a way in
which riAFT-BART can be used to select important predictors among incomplete,
clustered and censored survival observations. This strategy is general enough
to accommodate any missing data pattern and is highly flexible as it leverages
BART to model the relationship between survival times and covariates. As
demonstrated by both simulations and the case study, this elevated modeling
flexibility gave rise to substantially better variable selection results,
particularly in terms of identifying predictors of complex functional forms,
compared to existing variable selection methods suitable for clustered
survival data.

There are several important avenues for further research. First, the riAFT-
BART model could be extended to include the random slopes and accommodate the
cluster-level covariates. Correspondingly, a variable selection method can be
developed to identify important predictors at both the individual level and
the cluster level. Second, using flexible modeling techniques alone will not
address violations of the no unmeasured confounding assumption required for
causal inference. Developing a formal sensitivity analysis approach
(hu2022aflexible, ) for unmeasured confounding would be a worthwhile and
important contribution. Finally, our bootstrap imputation-based variable
selection procedure is can be computationally demanding on a alrge dataset.
Although the bootstrap resampling can be computed in parallel on multiple
cores when such resources are available, it would be worthwhile to develop
strategies that could offer substantial computational savings.

## Data Availability Statement

R codes to implement our proposed methods and the comparison methods, and to
replicate our simulation studies are provided in the GitHub page of the author
https://github.com/liangyuanhu/TEH-VS-riAFTBART, and are also available in the
R package riAFTBART. Access to the COVID-19 data used in the case study needs
to be requested and approved by the Icahn School of Medicine at Mount Sinai.
The authors declare that they have no conflict of interest.

This work was supported in part by award ME_2017C3_9041 from the Patient-
Centered Outcomes Research Institute, and by grants R21CA245855 and
1R01HL159077-01A1 from the National Institute of Health.

## References

  * (1) Hu L, Ji J, Ennis RD et al.  A flexible approach for causal inference with multiple treatments and clustered survival outcomes.  _arXiv preprint_ 2022; ArXiv:2202.08318.
  * (2) Chipman HA, George EI and McCulloch RE.  BART: Bayesian additive regression trees.  _The Annals of Applied Statistics_ 2010; 4(1): 266–298.
  * (3) Hu L, Ji J and Li F.  Estimating heterogeneous survival treatment effect in observational data using machine learning.  _Statistics in Medicine_ 2021; 40(21): 4691–4713.
  * (4) Fan J and Li R.  Variable selection for Cox’s proportional hazards model and frailty model.  _The Annals of Statistics_ 2002; 30(1): 74–99.
  * (5) Androulakis E, Koukouvinos C and Vonta F.  Estimation and variable selection via frailty models with penalized likelihood.  _Statistics in Medicine_ 2012; 31(20): 2223–2239.
  * (6) Utazirubanda JC, M León T and Ngom P.  Variable selection with Group LASSO approach: Application to Cox regression with frailty model.  _Communications in Statistics-Simulation and Computation_ 2021; 50(3): 881–901.
  * (7) Bender A, Groll A and Scheipl F.  A generalized additive model approach to time-to-event analysis.  _Statistical Modelling_ 2018; 18(3-4): 299–321.
  * (8) Marra G and Wood SN.  Practical variable selection for generalized additive models.  _Computational Statistics & Data Analysis_ 2011; 55(7): 2372–2387.
  * (9) Rondeau V, Commenges D and Joly P.  Maximum penalized likelihood estimation in a gamma-frailty model.  _Lifetime Data Analysis_ 2003; 9(2): 139–153.
  * (10) Marquardt DW.  An algorithm for least-squares estimation of nonlinear parameters.  _Journal of the Society for Industrial and Applied Mathematics_ 1963; 11(2): 431–441.
  * (11) Ha ID, Pan J, Oh S et al.  Variable selection in general frailty models using penalized h-likelihood.  _Journal of Computational and Graphical Statistics_ 2014; 23(4): 1044–1060.
  * (12) Hu L, Joyce Lin JY and Ji J.  Variable selection with missing data in both covariates and outcomes: Imputation and machine learning.  _Statistical Methods in Medical Research_ 2021; 30(12): 2651–2671.
  * (13) Bleich J, Kapelner A, George EI et al.  Variable selection for BART: an application to gene regulation.  _The Annals of Applied Statistics_ 2014; 8(3): 1750–1781.
  * (14) Ishwaran H, Kogalur UB, Gorodeski EZ et al.  High-dimensional variable selection for survival data.  _Journal of the American Statistical Association_ 2010; 105(489): 205–217.
  * (15) Lagani V and Tsamardinos I.  Structure-based variable selection for survival data.  _Bioinformatics_ 2010; 26(15): 1887–1894.
  * (16) Hu L and Hogan JW.  Causal comparative effectiveness analysis of dynamic continuous-time treatment initiation rules with sparsely measured outcomes and death.  _Biometrics_ 2019; 75(2): 695–707.
  * (17) Chen PY and Tsiatis AA.  Causal inference on the difference of the restricted mean lifetime between two groups.  _Biometrics_ 2001; 57(4): 1030–1038.
  * (18) Arpino B and Cannas M.  Propensity score matching with clustered data. an application to the estimation of the impact of caesarean section on the apgar score.  _Statistics in Medicine_ 2016; 35(12): 2074–2091.
  * (19) Hu L, Gu C, Lopez M et al.  Estimation of causal effects of multiple treatments in observational studies with a binary outcome.  _Statistical Methods in Medical Research_ 2020; 29: 3218–3214.
  * (20) Hu L and Gu C.  Estimation of causal effects of multiple treatments in healthcare database studies with rare outcomes.  _Health Services and Outcomes Research Methodology_ 2021; 21(3): 287–308.
  * (21) Gelman A, Van Dyk DA, Huang Z et al.  Using redundant parameterizations to fit hierarchical models.  _Journal of Computational and Graphical Statistics_ 2008; 17(1): 95–122.
  * (22) Henderson NC, Louis TA, Rosner GL et al.  Individualized treatment effects with censored data via fully nonparametric bayesian accelerated failure time models.  _Biostatistics_ 2020; 21(1): 50–68.
  * (23) Foster JC, Taylor JM and Ruberg SJ.  Subgroup identification from randomized clinical trial data.  _Statistics in Medicine_ 2011; 30(24): 2867–2880.
  * (24) Logan BR, Sparapani R, McCulloch RE et al.  Decision making and uncertainty quantification for individualized treatments using Bayesian Additive Regression Trees.  _Statistical Methods in Medical Research_ 2019; 28(4): 1079–1093.
  * (25) Hu L, Lin JYJ, Sigel K et al.  Estimating heterogeneous survival treatment effects of lung cancer screening approaches: A causal machine learning analysis.  _Annals of Epidemiology_ 2021; 62: 36–42.
  * (26) Breiman L.  Random forests.  _Machine Learning_ 2001; 45(1): 5–32.
  * (27) Deng H.  Interpreting tree ensembles with intrees.  _International Journal of Data Science and Analytics_ 2019; 7(4): 277–287.
  * (28) VanderWeele TJ.  Principles of confounder selection.  _European Journal of Epidemiology_ 2019; 34(3): 211–219.
  * (29) Lin JYJ, Hu L, Huang C et al.  A flexible approach for variable selection in large-scale healthcare database studies with missing covariate and outcome data.  _BMC Medical Research Methodology_ 2022; 22: 132.
  * (30) Hu L, Liu B and Li Y.  Ranking sociodemographic, health behavior, prevention, and environmental factors in predicting neighborhood cardiovascular health: a Bayesian machine learning approach.  _Preventive Medicine_ 2020; 141: 106240.
  * (31) Huang J and Ma S.  Variable selection in the accelerated failure time model via the bridge method.  _Lifetime Data Analysis_ 2010; 16(2): 176–195.
  * (32) Hu L, Li L and Ji J.  Machine learning to identify and understand key factors for provider-patient discussions about smoking.  _Preventive Medicine Reports_ 2020; 20: 101238.
  * (33) Hu L, Liu B, Ji J et al.  Tree-based machine learning to identify and understand major determinants for stroke at the neighborhood level.  _Journal of the American Heart Association_ 2020; 9(22): e016745.
  * (34) Tibshirani R.  Regression shrinkage and selection via the Lasso.  _Journal of the Royal Statistical Society: Series B (Methodological)_ 1996; 58(1): 267–288.
  * (35) Fan J and Li R.  Variable selection via nonconcave penalized likelihood and its oracle properties.  _Journal of the American Statistical Association_ 2001; 96(456): 1348–1360.
  * (36) Lee Y and Oh HS.  A new sparse variable selection via random-effect model.  _Journal of Multivariate Analysis_ 2014; 125: 89–99.
  * (37) Wood AM, White IR and Royston P.  How should variable selection be performed with multiply imputed data?  _Statistics in Medicine_ 2008; 27(17): 3227–3246.
  * (38) Cortinas Abrahantes J and Burzykowski T.  A version of the EM algorithm for proportional hazard model with random effects.  _Biometrical Journal_ 2005; 47(6): 847–862.
  * (39) Long Q and Johnson BA.  Variable selection in the presence of missing data: resampling and imputation.  _Biostatistics_ 2015; 16(3): 596–610.
  * (40) Cai J and Zeng D.  Additive mixed effect model for clustered failure time data.  _Biometrics_ 2011; 67(4): 1340–1351.
  * (41) Royston P and Parmar MK.  Restricted mean survival time: an alternative to the hazard ratio for the design and analysis of randomized trials with a time-to-event outcome.  _BMC Medical Research Methodology_ 2013; 13(1): 1–15.
  * (42) Van der Laan MJ, Polley EC and Hubbard AE.  Super learner.  _Statistical Applications in Genetics and Molecular Biology_ 2007; 6(1).
  * (43) Lu M, Sadiq S, Feaster DJ et al.  Estimating individual treatment effect in observational data using random forest methods.  _Journal of Computational and Graphical Statistics_ 2018; 27(1): 209–219.
  * (44) Feng P, Zhou XH, Zou QM et al.  Generalized propensity score for estimating the average treatment effect of multiple treatments.  _Statistics in Medicine_ 2012; 31(7): 681–697.
  * (45) Schouten RM, Lugtig P and Vink G.  Generating missing values for simulation purposes: a multivariate amputation procedure.  _Journal of Statistical Computation and Simulation_ 2018; 88(15): 2909–2930.
  * (46) Van Buuren S.  _Flexible imputation of missing data_.  2nd ed. Boca Raton, FL: Chapman & Hall, 2018.
  * (47) Van Buuren, Stef and Groothuis-Oudshoorn, Karin.  mice: Multivariate imputation by chained equations in R.  _Journal of Statistical Software_ 2010; 45(3): 1–67.
  * (48) Hu L, Li F, Ji J et al.  Estimating the causal effects of multiple intermittent treatments with application to COVID-19.  _arXiv preprint_ 2021; ArXiv:2109.13368.
  * (49) Harrell Jr FE, Lee KL and Mark DB.  Multivariable prognostic models: Issues in developing models, evaluating assumptions and adequacy, and measuring and reducing errors.  _Statistics in Medicine_ 1996; 15(4): 361–387.
  * (50) Zhou X and Reiter JP.  A note on Bayesian inference after multiple imputation.  _The American Statistician_ 2010; 64(2): 159–163.
  * (51) Hu L, Zou J, Gu C et al.  A flexible sensitivity analysis approach for unmeasured confounding with multiple treatments and a binary outcome with application to SEER-Medicare lung cancer data.  _Annals of Applied Statistics_ 2022; 16(2): 1014–1037.
  * (52) Kent DM, Paulus JK, Van Klaveren D et al.  The predictive approaches to treatment effect heterogeneity (path) statement.  _Annals of Internal Medicine_ 2020; 172(1): 35–45.