Hu et al

*Liangyuan Hu, PhD
683 Hoes Lane West, Piscataway, NJ 08854

# A flexible approach for causal inference with multiple treatments and
clustered survival outcomes

Liangyuan Hu  Jiayi Ji  Ronald D. Ennis  Joseph W. Hogan  Department of
Biostatistics and Epidemiology, Rutgers University, Piscataway, NJ 08854,
U.S.A.  Department of Radiation Oncology, Cancer Institute of New Jersey of
Rutgers University, New Brunswick, NJ 08903, U.S.A.  Department of
Biostatistics, Brown University, Providence, RI 02912, U.S.A.
<EMAIL_ADDRESS>Hu L  Ji J  Ennis D. R  Hogan W. J

(<day> <Month>, <year>; <day> <Month>, <year>; <day> <Month>, <year>)

###### Abstract

[Abstract]When drawing causal inferences about the effects of multiple
treatments on clustered survival outcomes using observational data, we need to
address implications of the multilevel data structure, multiple treatments,
censoring and unmeasured confounding for causal analyses. Few off-the-shelf
causal inference tools are available to simultaneously tackle these issues. We
develop a flexible random-intercept accelerated failure time model, in which
we use Bayesian additive regression trees to capture arbitrarily complex
relationships between censored survival times and pre-treatment covariates and
use the random intercepts to capture cluster-specific main effects. We develop
an efficient Markov chain Monte Carlo algorithm to draw posterior inferences
about the population survival effects of multiple treatments and examine the
variability in cluster-level effects. We further propose an interpretable
sensitivity analysis approach to evaluate the sensitivity of drawn causal
inferences about treatment effect to the potential magnitude of departure from
the causal assumption of no unmeasured confounding. Expansive simulations
empirically validate and demonstrate good practical operating characteristics
of our proposed methods. Applying the proposed methods to a dataset on older
high-risk localized prostate cancer patients drawn from the National Cancer
Database, we evaluate the comparative effects of three treatment approaches on
patient survival, and assess the ramifications of potential unmeasured
confounding. The methods developed in this work are readily available in the R
package riAFTBART.

, , , (2022), A flexible approach for causal inference with multiple
treatments and clustered survival outcomes, Statistics in Medicine,
2022;00:1–18.

###### keywords:

Observational studies; Bayesian machine learning; Sensitivity analysis;
Multilevel survival data

††articletype: Article Type

##  1 Introduction

In cancer research, decision makers are starting to rely more heavily on real
world evidence because clinical trials can be enormously expensive, time
consuming and restrictive. The increasing availability of observational data
sources like large registries and electronic health records provides new
opportunities to obtain real world comparative effectiveness evidence. Recent
efforts have been made to evaluate the comparative effectiveness of multiple
treatment approaches on patient survival for high-risk localized prostate
cancer using the large-scale national cancer database (NCDB). 1, 2, 3, 4 The
complex data structures, however, pose three main challenges for statistical
analyses that have not been fully addressed in the extant literature.

First, multiple (i.e., more than two) active treatments are involved. For
high-risk localized prostate cancer, there are three popular treatment
approaches: radical prostatectomy (RP), external beam radiotherapy (EBRT)
combined with androgen deprivation (AD) (EBRT+AD) and EBRT plus brachytherapy
with or without AD (EBRT+brachy$\pm$AD). Each treatment option has
historically been performed on different types of patients as various health
and demographic factors are closely linked to treatment choice. Second,
treatment and patient information on a large, national population is collected
from various treating facilities. The NCDB hospital registry data are
collected in more than 1,500 Commission on Cancer accredited facilities. The
participating institutions are not selected at random, and there can be
substantial institutional variation in treatment effect. Third, some important
confounders (pre-treatment variables predicting both treatment and outcome)
may not be collected in the observational data. Two known main confounders in
high-risk prostate cancer are the number of positive cores and magnetic
resonance imaging findings, which are often not collected in the NCDB data or
other large cancer registries alike.

Despite numerous recent advances in causal inference, the literature for
handling data structures of this type – which arise frequently in population
cancer research – is sparse. Causal inference techniques traditionally focus
on a binary treatment. There is now a substantial body of research on causal
inference methods with multiple treatments and a continuous outcome 5, 6, 7 or
a binary outcome. 8, 9, 10 Ennis et al. 1 and Zeng et al. 4 described
propensity score weighting based methods for drawing inferences about causal
effects of multiple treatments on censored survival outcomes, but neither work
considered the multilevel data structure. There are two main reasons why it
might be important to account for the cluster-level, or institutional
variation when estimating treatment effect in the general population. First,
if there are substantial institutional effects, then a causal analysis
ignoring institution would be based on an incorrect model, which can lead to
invalid inferences about treatment effects. Second, neither participating
institutions nor patients in the registry data were selected at random. With
substantial institutional variation, it is unclear exactly what treatment
effect would be seen in the general patient population across various
institutions with different outcomes. Finally, inferring causal links from
observational data inevitably involves the untestable assumption of no
unmeasured confounding, which holds that all pre-treatment variables are
sufficient to predict both treatment and outcome. If there are unmeasured
confounders, it is important to evaluate how departures from the no unmeasured
confounding assumption might alter causal conclusions. Sensitivity analysis is
useful to address this causal assumption, and is recommended by the
Strengthening the Reporting of Observational Studies in Epidemiology (STROBE)
guidelines. 11 Existing sensitivity analysis approaches have largely focused
on a binary treatment. Hu et al. 12 recently proposed a flexible Monte Carlo
sensitivity analysis approach in the context of multiple treatments and a
binary outcome. There is still sparse literature on sensitivity analysis
methods that simultaneously accommodate multiple treatments and multilevel
censored survival outcomes.

To fill these research gaps, we propose a flexible approach for drawing causal
inferences about multiple treatments while respecting the multilevel survival
data structure, and develop an interpretable sensitivity analysis approach to
evaluate how the drawn causal conclusions might be altered in response to the
potential magnitude of departure from the no unmeasured confounding
assumption. We propose a random-intercept accelerated failure time (AFT) model
utilizing Bayesian additive regression trees (BART) 13, termed as riAFT-BART.
In this model, we use the random intercepts for cluster-specific main effects
capturing the variation across the institutions, and leverage the flexibility
of the BART model 8, 14, 15, 16 to accurately capture arbitrarily complex
relationships among survival times, treatments and covariates. We then develop
an efficient Markov chain Monte Carlo algorithm to draw posterior inferences
about the population survival effects of multiple treatments. We further
propose an interpretable sensitivity analysis approach leveraging our riAFT-
BART model in the context of multiple treatments and clustered survival
outcomes. Finally, we apply the proposed methods to the NCDB data and
elucidate the causal effects of three treatment approaches (RP, EBRT+AD and
EBRT+brachy$\pm$AD) on patient survival among older and high-risk prostate
cancer patients and the impact of unmeasured confounding, as well as examine
the institutional effects.

The rest of the paper is organized as follows. Section 2 describes notation
and proposes the riAFT-BART method for the estimation of causal effects.
Section 3 describes a corresponding sensitivity analysis approach using the
model developed in Section 2. Section 4 develops a wide variety of simulation
scenarios to examine the practical operating characteristics of our proposed
methods against three alternative methods, and presents findings. In Section
5, we apply our methods to NCDB data to estimate the causal effects of three
treatment approaches on patient survival among high-risk localized prostate
cancer patients, and perform a sensitivity analysis to evaluate how sensitive
the estimated causal effects are to possible unmeasured confounding. Section 6
concludes with a discussion.

##  2 Estimation of causal effects

###  2.1 Notation, definitions and assumptions

Consider an observational study possessing a two-level data structure that has
$K$ clusters (institutions), each having treated $n_{k}$ individuals, indexed
by $i=1,\ldots,n_{k},k=1,\ldots,K$. The total number of individuals in the
study is $N=\sum_{k=1}^{K}n_{k}$. Our goal is to infer the causal effect of
treatment $A\in\mathscr{A}=\\{a_{1},\ldots,a_{J}\\}$ on time to failure $T$,
where $J$ is the total number of treatment options. For each individual $i$ in
cluster $k$, there is a vector of pre-treatment measured covariates
$\bm{X}_{ik}$, and let $T_{ik}$ be the individual’s failure time, which 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. There
are no cluster-level covariates in our study, but our work can be extended to
include them in $V_{k}$. We proceed in 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}$. We
similarly define $C_{ik}(a_{j})$ as the counterfactual censoring time under
treatment $a_{j}$. Throughout, we maintained the standard assumptions for
drawing causal inference with observational clustered survival data: 3, 17, 18

  1. (A1)

Consistency: the observed failure time
$T_{ik}=\sum_{j=1}^{J}T_{ik}(a_{j})I(A_{ik}=a_{j})$ and censoring time
$C_{ik}=\sum_{j=1}^{J}C_{ik}(a_{j})I(A_{ik}=a_{j})$, where $I(\cdot)$ is the
usual indicator function;

  2. (A2)

Weak unconfoundedness: $T_{ik}(a_{j})\amalg A_{ik}\mid\bm{X}_{ik},V_{k}\text{
for }A_{ik}=a_{j},a_{j}\in\mathscr{A}$;

  3. (A3)

Positivity: the generalized propensity score for treatment assignment
$e(\bm{X}_{ik},V_{k})=P(A_{ik}=a_{j}|\bm{X}_{ik},V_{k})$ is bounded away from
0 and 1 for $A_{ik}=a_{j},a_{j}\in\mathscr{A}$;

  4. (A4)

Covariate-dependent censoring: $T_{ik}(a_{j})\amalg
C_{ik}(a_{j})\mid\bm{X}_{ik},V_{k},A_{ik}$, for
$A_{ik}=a_{j},a_{j}\in\mathscr{A}$.

The counterfactual outcomes are linked to the observed outcomes via Assumption
(A1). It allows us to write
$Y_{ik}=\sum_{j=1}^{J}Y_{ik}(a_{j})I(A_{ik}=a_{j})$, where
$Y_{ik}(a_{j})=\min\left(T_{ik}(a_{j}),C_{ik}(a_{j})\right)$. Similarly,
$\Delta_{ik}=\sum_{j=1}^{J}\Delta_{ik}(a_{j})I(A_{ik}=a_{j})$, where
$\Delta_{ik}(a_{j})=I\\{T_{ik}(a_{j})<C_{ik}(a_{j})\\}$. Assumption (A2) is
referred to as the “no unmeasured confounding” assumption. Because this is an
untestable assumption, we will develop a sensitivity analysis approach in
Section 3 to gauge the impact of violations of this assumption. Assumption
(A3) requires that the treatment assignment is not deterministic within each
strata formed by the covariates. 19 This assumption can be directly assessed
by visualizing the distribution of estimated generalized propensity scores.
Finally, Assumption (A4) states that the counterfactual survival time is
independent of the counterfactual censoring time given pre-treatment
covariates, cluster-level variables and treatment variable. This condition
directly implies $T_{ik}\amalg C_{ik}\mid\bm{X}_{ik},V_{k},A_{ik}$, and is
akin to the (conditionally) independent censoring assumption in the
traditional survival analysis literature. 20

We define the causal estimands directly in terms of counterfactual survival
times. Alternatively, one can define causal estimands based on functionals
(e.g., median) of the counterfactual survival curves. 3 In this paper, we
focus on the average treatment effect (ATE) defined either over the sample or
the population. Consider a pairwise comparison between treatments $a_{j}$ and
$a_{j^{\prime}}$. Common sample estimands are the sample average treatment
effect (SATE),

$\frac{1}{N}\sum_{k=1}^{K}\sum_{i=1}^{n_{k}}[T_{ik}(a_{j})-T_{ik}(a_{j^{\prime}})].$

Common population estimands are the population average treatment effect
(PATE),

$E[T(a_{j})-T(a_{j^{\prime}})].$

Conditional average treatment effect (CATE)

$\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},V_{k}\right]$

is another estimand that preserves some of the properties of the previous two.
14 As our methods are developed in a Bayesian framework, CATE is a natural
estimand to use in this paper. 14 We obtain the sample marginal effects by
averaging the individual conditional expectation of the counterfactual
survival times
$E[T_{ik}(a_{j})-T_{ik}(a_{j^{\prime}}){\,|\,}\bm{X}_{ik},V_{k}]$ across the
empirical distribution of $\\{\bm{X}_{ik},V_{k}\\}_{i=1,k=1}^{n_{k},K}$. 8
Another causal estimand of interest is the average treatment effect on the
treated (ATT). By averaging the differenced counterfactual survival times over
those in the reference group, we can define the ATT counterparts of all three
estimands described above. For brevity of exposition, we focus on CATE in this
paper, but our methods can be straightforwardly extended for the ATT effects.
For example, the conditional average treatment effect among those who received
treatment $a_{j}$ $\text{CATT}_{a_{j}{\,|\,}a_{j},a_{j^{\prime}}}$ is

$\frac{1}{N_{j}}\mathop{\sum\sum}_{i,k:A_{ik}=a_{j}}E\left[T_{ik}(a_{j})-T_{ik}(a_{j^{\prime}}){\,|\,}\bm{X}_{ik},V_{k}\right],$

where $N_{j}=\sum_{i=1}^{n_{k}}\sum_{k=1}^{K}I(A_{ik}=a_{j})$ is the size of
the reference group $a_{j}$.

###  2.2 The riAFT-BART model for clustered survival data

For clustered survival data, we propose the following random-intercept AFT
model utilizing the likelihood-based machine learning technique BART:

$\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 treatments
and covariates to survival times $T_{ik}$, $b_{k}$’s are the random intercepts
for cluster-specific main effects capturing the institutional variation, and
$\epsilon_{ik}$ is the residual term. We use BART to flexibly model the
unknown function $f$ by a sum of shallow trees
$f(A_{ik},\bm{X}_{ik})=\sum_{h=1}^{H}g(A_{ik},\bm{X}_{ik};\mathcal{W}_{h},\mathcal{M}_{h})$,
where $\mathcal{W}_{h}$ is the $h$th binary tree structure,
$\mathcal{M}_{h}=(\mu_{1h},\ldots,\mu_{{c_{h}}h})^{T}$ is the set of $c_{h}$
terminal node parameters associated with tree structure $\mathcal{W}_{h}$. For
a given value $A_{ik}$ and $\bm{X}_{ik}$ in the predictor space, the
assignment function $g(A_{ik},\bm{X}_{ik};\mathcal{W}_{h},\mathcal{M}_{h})$
returns the parameter $\mu_{lh},l\in\\{1,\ldots,c_{h}\\}$ associated with the
terminal node of the predictor subspace in which $\\{A_{ik},\bm{X}_{ik}\\}$
falls. We will place regularizing priors on
$\\{\mathcal{W}_{h},\mathcal{M}_{h}\\}$ to keep the impact of each individual
tree on the overall fit small and prevent overfitting. 13, 14, 21 We assume a
mean-zero normal distribution for and independence between $b_{k}$ and
$\epsilon_{ik}$ with variance $\alpha_{k}\tau^{2}$ and $\sigma^{2}$,
respectively. Here we adopt the parameter expansion technique 22 and introduce
a redundant parameter $\alpha_{k}$ for the variance of $b_{k}$ to improve
computational performance of our proposed Markov chain Monte Carlo (MCMC)
algorithm (Section 2.3).

Model (LABEL:eq:_riAFT-BART) has two main advantages. First, it allows for
direct specification of treatment effect on changes in life expectancy, which
substantially facilitates the interpretability of an sensitivity analysis.
Second, unlike the proportional hazards regression, the AFT model formulation
can naturally incorporates BART to not only flexibly capture arbitrarily
complex functional form of $f(A_{ik},\bm{X}_{ik})$ but also provide coherent
inferences based on a probability model and proper representations of
uncertainty intervals via the posterior.

We decompose the joint prior distribution as

$\displaystyle
P\left[\left(\mathcal{W}_{1},\mathcal{M}_{1}\right),\ldots,\left(\mathcal{W}_{H},\mathcal{M}_{H}\right),\sigma,\tau,\alpha_{k}\right]$
$\displaystyle=\left[\prod_{h=1}^{H}P\left(\mathcal{W}_{h},\mathcal{M}_{h}\right)\right]P(\sigma)P(\tau)P(\alpha_{k})$
$\displaystyle=\left[\prod_{h=1}^{H}P\left(\mathcal{M}_{h}{\,|\,}\mathcal{W}_{h}\right)P\left(\mathcal{W}_{h}\right)\right]P(\sigma)P(\tau)P(\alpha_{k})$
$\displaystyle=\left[\prod_{h=1}^{H}\left\\{\prod_{l=1}^{c_{h}}P\left(\mu_{lh}\mid\mathcal{W}_{h}\right)\right\\}P\left(\mathcal{W}_{h}\right)\right]P(\sigma)P(\tau)P(\alpha_{k}).$

Following Chipman et al. 13 and Hendersen et al., 23 we center the observed
responses $y_{ik}$ via the following two steps: (i) fit a parametric
intercept-only AFT model assuming log-normal residuals, and estimate the
intercept $\hat{\mu}_{AFT}$ and the residual scale $\hat{\sigma}_{AFT}$; (ii)
transform the responses as
$y_{ik}^{cent}=y_{ik}\exp\left(-\hat{\mu}_{AFT}\right)$. Then for the terminal
node values $\mu_{lh}$, we place the prior distribution $\mu_{lh}\sim
N\left(0,\xi^{2}/(4Hk^{2})\right)$, where $k$ is the terminal node $\mu_{lh}$
hyperparameter and $\xi=4\hat{\sigma}_{AFT}$. This prior induces a
$N(0,4\hat{\sigma}_{AFT}^{2}/k^{2})$ prior on the regression function
$f(\cdot)$ in model (LABEL:eq:_riAFT-BART), and with a default setting of
$k=2$ assigns 95% prior probability to the interval
$[-2\hat{\sigma}_{AFT},2\hat{\sigma}_{AFT}]$, which is sensible.23 As
suggested in Chipman et al., 13, we use an inverse gamma distribution
$IG(\frac{\nu}{2},\frac{\nu\lambda}{2})$ as the prior for $\sigma^{2}$ and the
default value for $\nu$ ($\nu=3$), and defer to the defauls for other
hyperparameters of the BART trees $\\{\mathcal{W}_{h},\mathcal{M}_{h}\\}$. We
place a prior $IG(1,1)$ on $\tau^{2}$ and $\alpha_{k}$.22 For initial values,
we first set an initial random intercept $\hat{b}_{k}^{(0)}$ to be the mean of
the lognormal residuals from a parametric AFT model with $\bm{X}_{ik}$ as the
predictors for each cluster $k$. We then set the initial value
${\sigma}^{(0)}$ to be the standard deviation of the pooled model residuals
over the $K$ clusters, and initialize $\lambda$ as the value such that
$P\left(\sigma<{\sigma}^{(0)};\nu,\lambda^{(0)}\right)=0.9$.

We use data augmentation to deal with right censoring. 23, 24 Working with the
centered responses $y_{ik}^{cent}$, when $\Delta_{ik}=0$, we impute the
unobserved and centered survival times $z_{ik}$ from a truncated normal
distribution:

$\left[\log Z_{ik}{\,|\,}\log Z_{ik}>\log y^{cent}_{ik}\right]\sim N_{(\log
y_{ik}^{cent},\infty)}\left(f(A_{ik},\bm{X}_{ik})+b_{k},\sigma^{2}\right)$

in each Gibbs iteration, where $\log Z_{ik}\sim
N\left(f(A_{ik},\bm{X}_{ik})+b_{k},\sigma^{2}\right)$. The centered complete-
data survival times are

$y^{cent,c}_{ik}=\begin{cases}y^{cent}_{ik}\text{ \; \; if }\Delta_{ik}=1\\\
z_{ik}\text{ \; \; if }\Delta_{ik}=0\end{cases}.$

###  2.3 Posterior inferences for treatment effects

Here we employ a Metropolis within Gibbs procedure for posterior inferences
about treatment effects on patient survival. Using the centered complete-data
survival times $y_{ik}^{cent,c}$, the joint posterior is

$\displaystyle\hskip
12.0ptP\left(b_{k},\tau^{2},\alpha_{k},\mu_{lh},\sigma^{2}{\,|\,}y^{cent,c}_{ik},\bm{X}_{ik},A_{ik},V_{k},\\{\mathcal{W}_{h},\mathcal{M}_{h}\\}\right)$
$\displaystyle\propto
P\left(y^{cent,c}_{ik}{\,|\,}\bm{X}_{ik},A_{ik},V_{k},b_{k},\tau^{2},\alpha_{k},\sigma^{2},\\{\mathcal{W}_{h},\mathcal{M}_{h}\\}\right)P\left(b_{k}{\,|\,}\tau^{2},\alpha_{k}\right)P\left(\tau^{2}\right)P\left(\alpha_{k}\right)P\left(\mu_{lh}\right)P\left(\sigma^{2}\right).$

We can draw the values of BART sum-of-trees model parameters, $\mu_{lh}$ and
$\sigma^{2}$, directly from the fitted BART model. Their posterior
distributions
$P\left(\mu_{lh}{\,|\,}y^{cent,c}_{ik},\bm{X}_{ik},A_{ik},V_{k},b_{k},\tau^{2},\alpha_{k},\sigma^{2},\\{\mathcal{W}_{h}\\}\right)$
and
$P\left(\sigma^{2}{\,|\,}y^{cent,c}_{ik},\bm{X}_{ik},A_{ik},V_{k},b_{k},\tau^{2},\alpha_{k},\\{\mathcal{W}_{h},\mathcal{M}_{h}\\}\right)$
are presented in Web Section S1. We can show the posterior distribution of the
random intercept $b_{k}$ is

$\left[b_{k}{\,|\,}y^{cent,c}_{ik},\bm{X}_{ik},A_{ik},V_{k},\tau^{2},\alpha_{k},\sigma^{2},\\{\mathcal{W}_{h},\mathcal{M}_{h}\\}\right]\sim
N\left(\dfrac{\tau^{2}\alpha_{k}\sum_{i=1}^{n_{k}}\left(y^{cent,c}_{ik}-\hat{f}(\bm{X}_{ik},A_{ik})\right)}{n_{k}\tau^{2}\alpha_{k}+\sigma^{2}},\dfrac{\sigma^{2}\tau^{2}\alpha_{k}}{n_{k}\tau^{2}\alpha_{k}+\sigma^{2}}\right).$

The posterior of $\alpha_{k}$, used for parameter expansion, is

$\left[\alpha_{k}{\,|\,}y^{cent,c}_{ik},\bm{X}_{ik},A_{ik},V_{k},\tau^{2},b_{k},\sigma^{2},\\{\mathcal{W}_{h},\mathcal{M}_{h}\\}\right]\sim
IG\left(1,1+\frac{\sum_{k=1}^{K}b_{k}^{2}}{2\tau^{2}}\right).$

We obtain the posterior of $\tau^{2}$ as

$\left[\tau^{2}{\,|\,}y^{cent,c}_{ik},\bm{X}_{ik},A_{ik},V_{k},b_{k},\alpha_{k},\sigma^{2},\\{\mathcal{W}_{h},\mathcal{M}_{h}\\}\right]\sim
IG\left(\dfrac{K}{2}+1,\dfrac{\sum_{k=1}^{K}b_{k}^{2}+2\alpha_{k}}{2\alpha_{k}}\right).$

Complete derivation of the posterior distributions are provided in Web Section
S1.

We now describe our Metropolis within Gibbs procedure to draw from the
posterior distribution of our proposed riAFT-BART model (LABEL:eq:_riAFT-
BART). A _single iteration_ of our sampling algorithm proceeds through the
following steps:

Algorithm 1 A single iteration of riAFT-BART sampling algorithm

  1. 1.

Update $b_{k}$, $\tau^{2}$ and $\alpha_{k}$ from their respective posterior
distributions.

  2. 2.

Using $\log y^{cent,c}_{ik}-b_{k}$ as the responses and
$\\{A_{ik},\bm{X}_{ik}\\}$ as the covariates, update BART sum-of-trees model
via parameters $\mu_{lh}$ and $\sigma^{2}$, using the Bayesian backfitting
approach of Chipman et al. 13 Directly update $f(A_{ik},\bm{X}_{ik})$ using
the updated BART model, for $i=1,\ldots,n_{k},k=1,\ldots,K$.

  3. 3.

For each $\\{i,k\\}\in\\{i=1,\ldots,n_{k},k=1,\ldots,K\\}$, update $z_{ik}$ by
sampling

$\log z_{ik}\sim\text{Truncated-
Normal}\left(f(A_{ik},\bm{X}_{ik})+b_{k},\sigma^{2};\log
y^{cent}_{ik}\right).$

Because we use the centered responses
$\log(y^{cent}_{ik})=\log(y_{ik})-\hat{\mu}_{AFT}$ in posterior computation,
we add $\hat{\mu}_{AFT}$ back to the posterior draws of
$f(A_{ik},\bm{X}_{ik})$ in the final output.

To draw posterior inferences about the CATE effects via riAFT-BART, we note
that under the causal assumptions (A1)–(A4),

$\displaystyle\begin{split}&\hskip 13.0ptE\left\\{\log T_{ik}(a_{j})-\log
T_{ik}(a_{j^{\prime}})\right\\}\\\ &=E_{\bm{x}_{ik},b_{k}}\left\\{E\left(\log
T_{ik}{\,|\,}A_{ik}=a_{j},\bm{X}_{ik}=\bm{x}_{ik},b_{k}\right)-E\left(\log
T_{ik}{\,|\,}A_{ik}=a_{j^{\prime}},\bm{X}_{ik}=\bm{x}_{ik},b_{k}\right)\right\\}\\\
&=E_{\bm{x}_{ik},b_{k}}\left[E\left\\{\left(f\left(a_{j},\bm{x}_{ik}\right)+b_{k}+\epsilon_{ik}\right)-\left(f\left(a_{j^{\prime}},\bm{x}_{ik}\right)+b_{k}+\epsilon_{ik}\right)\right\\}\right]\\\
&=E_{\bm{x}_{ik}}\left[E\left\\{f\left(a_{j},\bm{x}_{ik}\right)-f\left(a_{j^{\prime}},\bm{x}_{ik}\right)\right\\}\right].\end{split}$
(2)

This allows us to estimate treatment effect via outcome modeling.
Specifically,

$\displaystyle\widehat{CATE}_{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\\{f^{d}\left(a_{j},\bm{x}_{ik}\right)-f^{d}\left(a_{j^{\prime}},\bm{x}_{ik}\right)\right\\},$
(3)

where $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\\{f^{d}\left(a_{j},\bm{x}_{ik}\right)-f^{d}\left(a_{j^{\prime}},\bm{x}_{ik}\right)\right\\},d=1,\ldots,D$.

##  3 Sensitivity analysis

###  3.1 Overview

The estimation of causal effects with observational data relies on the weak
unconfoundedness assumption (A2), which cannot be verified empirically.
Violations of this assumption can lead to biased treatment effect estimates.
One widely recognized way to address concerns about violations of this
assumption is sensitivity analysis. In fact, the STROBE guidelines recommend
observational studies be accompanied by sensitivity analysis investigating the
ramifications of potential unmeasured confounding. 11 Many sensitivity
analysis methods have been developed, including Rosenbaum’s $\Gamma$, 25
external adjustment, 26 confounding functions 27, 12 and the E value, 28 to
name a few. These methods differ in how unmeasured confounding is formulated
and parameterized. Sensitivity analysis approaches in the context of multiple
treatments and clustered censored survival outcomes are an underdeveloped
area.

With multilevel data, there can be unobserved confounders at both cluster- and
individual-level. It has been shown in the literature that with propensity
score based methods, the fixed-effects model for propensity score estimation
automatically controls for the effects of unmeasured cluster-level
confounders. 29, 30, 31 In situations where the cluster sizes are small, a
random-effects propensity score model may provide more accurate effect
estimates, but is reliant on inclusion of important cluster-level covariates
as regressors. As the cluster size increases, results from the random-effects
model converge to those from a corresponding fixed-effects model. Although
there is sparse literature on whether accounting for the clustered structure
in potential outcome models would protect against misspecification due to
cluster-level confounders, the outcome model is connected to the propensity
score model in that the sufficient statistics (treatment group means of
covariates) that must be balanced to eliminate confounding differences under
both models are the same. 30 Li et al. 30 conducted a simulation to show that
ignoring the clustered data structure in both the propensity score and outcome
models would lead to biased ATE estimates; respecting the structure in at
least one of the models gives consistent estimates. Based on these grounds and
that we are dealing with large clusters of a national database, we believe our
riAFT-BART model (LABEL:eq:_riAFT-BART) – which will converge to a fixed-
effects model with large cluster sizes – will represent heterogeneity in
cluster-level unmeasured confounding by the random effects $\\{b_{k}\\}$. We
then assume by conditioning on $\\{b_{k}\\}$ and $\\{\bm{X}_{ik}\\}$, the
potential outcome and treatment at cluster-level are independent, and
introduce our sensitivity analysis approach for individual-level unmeasured
confounding.

Our sensitivity analysis approach is along the line of work by Hu et al., 12
and is based within the framework of confounding function. 27, 32 The
confounding function based methods have the advantage of avoiding introducing
a hypothetical unmeasured confounder and making an assumption about its
underlying structure, on which there is a lack of consensus, and are preferred
when the primary interest is in understanding the total effect of all
unmeasured confounders.32, 12

###  3.2 Confounding function adjusted treatment effect estimates

For notational brevity, we suppress the $ik$ subscript denoting individual.
Following Brumback et al.32 and Hu et al., 12 we first define the confounding
function for any pair of treatments $(a_{j},a_{j^{\prime}})$ as

$\displaystyle c(a_{j},a_{j^{\prime}},\bm{x},v)=E\left[\log
T(a_{j}){\,|\,}A=a_{j},\bm{X}=\bm{x},V=v\right]-E\left[\log
T(a_{j}){\,|\,}A=a_{j^{\prime}},\bm{X}=\bm{x},V=v\right].$ (4)

This confounding function directly represents the difference in the mean
potential log survival times under treatment $a_{j}$ between those treated
with $a_{j}$ and those treated with $a_{j^{\prime}}$, who have the same level
of $\bm{x}$. Under the assumption of no unmeasured confounding, given measured
individual- and cluster-level covariates $\bm{X}$ and $V$, the potential
outcome is independent of treatment assignment. Had they received the same
treatment $a_{j}$, their mean potential survival times would have been the
same, or $c(a_{j},a_{j^{\prime}},\bm{x},v)=0$,
$\forall\\{a_{j},a_{j^{\prime}}\\}\in\mathscr{A}$. When this assumption is
violated and there exists unmeasured confounding, the causal effect estimates
using measured confounders will be biased. The bias in the estimated treatment
effect $\widehat{CATE}_{a_{j},a_{j^{\prime}}}$ takes the following form:

$\displaystyle\begin{split}\text{Bias}(a_{j},a_{j^{\prime}}{\,|\,}\bm{x},v)=&-p_{j}c(a_{j^{\prime}},a_{j},\bm{x},v)+p_{j^{\prime}}c(a_{j},a_{j^{\prime}},\bm{x},v)\\\
&-\sum\limits_{m:a_{m}\in\mathscr{A}\setminus\\{a_{j},a_{j^{\prime}}\\}}p_{m}\left\\{c(a_{j^{\prime}},a_{m},\bm{x},v)-c(a_{j},a_{m},\bm{x},v)\right\\},\end{split}$
(5)

where $p_{j}=P(A=a_{j}{\,|\,}\bm{X}=\bm{x},V=v)$, $j\neq
j^{\prime}\in\\{1,\ldots,J\\}$. A proof of this result is presented in Web
Section S2.

Given known confounding functions $c$, we can construct the confounding
function adjusted estimators by first modifying the actual survival time $T$,
and then estimating the causal effect by fitting our riAFT-BART model to
modified outcomes. In this way, the bias in equation (5) will be effectively
removed from the adjusted effect estimate. Because the survival times $T$ may
be right censored and we deal with right censoring using data augmentation,
the outcome modification can be implemented on the complete-data survival
times. We propose a Monte Carlo sensitivity analysis approach along the line
of work by Hu et al., 12 which was developed for binary outcomes. We extend
their work to accommodate multilevel censored survival outcomes. Our
sensitivity analysis proceeds with steps listed in Algorithm 2.

Algorithm 2 Sensitivity analysis algorithm

  1. 1.

Fit a multinomial probit BART model
$f^{\text{MBART}}(A_{ik}{\,|\,}\bm{X}_{ik},V_{k})$ to estimate the generalized
propensity scores, $p_{j}\equiv
P(A_{ik}=a_{j}{\,|\,}\bm{X}_{ik}=\bm{x}_{ik},V_{k}=v_{k})\;\forall
a_{j}\in\mathscr{A}$, for each individual.

  2. 2.

For each treatment $a_{j}\in\mathscr{A}$, draw $Q_{1}$ generalized propensity
scores $\tilde{p}_{m1},\ldots,\tilde{p}_{mQ_{1}},\forall m\neq j\wedge
a_{m}\in\mathscr{A}$ from the posterior predictive distribution of
$f^{\text{MBART}}(A_{ik}{\,|\,}\bm{X}_{ik},V_{k})$ for each individual.

  3. 3.

For $q\in\\{1,\ldots,Q_{1}\\}$, draw $Q_{2}$ values
$\gamma^{*}_{mq1},\ldots,\gamma^{*}_{mqQ_{2}}$ from the prior distribution of
each of the confounding functions $c(a_{j},a_{m},\bm{x},v)$, for each $m\neq
j\wedge a_{m}\in\mathscr{A}$.

  4. 4.

For each treatment $a_{j}$, adjust the centered complete-data survival times
(Step 2 of Algorithm 1) as follows:

$\displaystyle\log y_{ik}^{CF}\equiv\log y_{ik}^{cent,c}-\sum_{m\neq
j}^{J}P\left(A_{ik}=a_{m}{\,|\,}\bm{X}_{ik}=\bm{x},V_{k}=v\right)c(a_{j},a_{m},\bm{x},v),$
(6)

for each of $Q_{1}Q_{2}$ draws of
$\\{\tilde{p}_{m1},\gamma^{*}_{m11},\ldots,\gamma^{*}_{m1Q_{2}},\ldots,\tilde{p}_{mQ_{1}},\gamma^{*}_{mQ_{1}1},\ldots,\gamma^{*}_{mQ_{1}Q_{2}};m\neq
j\wedge a_{m}\in\mathscr{A}\\}$.

  5. 5.

Run Algorithm 1 for riAFT-BART on each of $Q_{1}\times Q_{2}$ sets of observed
data with $\log y_{ik}^{CF}$. Estimate the combined adjusted causal effects
and uncertainty intervals by pooling posterior samples across model fits
arising from the $Q_{1}\times Q_{2}$ data sets.

Note that in Step 1 of Algorithm 2, we can fit a flexible fixed-effects
multinomial probit BART model for generalized propensity scores. Steps 2 and 3
constitute a nested multiple imputation, 33 which is used to draw samples for
the product term
$P\left(A_{ik}=a_{m}{\,|\,}\bm{X}_{ik}=\bm{x},V_{k}=v\right)c(a_{j},a_{m},\bm{x},v)$
in equation (6). Step 4 “corrects" the complete-data survival times to adjust
the treatment effect estimate for individual-level unmeasured confounding.
This is because the causal effect is defined as the between-group difference
in mean potential outcomes and is estimated based on the observed outcomes. To
correct the bias in equation (5) due to individual-level unmeasured
confounding, we adjust the actual survival time $T$ of an individual who
received treatment $a_{j}$ as

$\displaystyle\log T^{CF}$ $\displaystyle=$ $\displaystyle\log
T-\left[E\left\\{\log T(a_{j}){\,|\,}a_{j},\bm{x},v\right\\}-E\left\\{\log
T(a_{j}){\,|\,}\bm{x},v\right\\}\right].$

Because the survival time $T$ can be right censored, we can replace the
centered complete-data survival time $\log y_{ik}^{cent,c}$ used in our riAFT-
BART sampling Algorithm 1 for treatment effect estimation with adjusted $\log
y_{ik}^{CF}$ as in equation (6). Web Section S2 provides a detailed
justification of this strategy for obtaining confounding function adjusted
causal effect estimates. Finally, we obtain the overall estimates of the
adjusted causal effect and sampling variance as the posterior mean and
variance of the pooled posterior samples across the $Q_{1}\times Q_{2}$ data
sets. 34 We used $Q_{1}=30$ and $Q_{2}=30$ when implementing our sensitivity
analysis Algorithm 2 in both simulation (Web Section S3) and case study
(Section 5).

The confounding functions are not identifiable from the observe data. We can
assume the form of confounding function to represent our prior beliefs about
the possible direction and magnitude of the effect of unmeasured confounding.
We follow strategies discussed in earlier work 27, 32, 35, 12 to specify the
signs and bounds of the confounding functions. For example, by assigning
$c(a_{j},a_{j^{\prime}},{\,|\,}\bm{x},v)>0$ and
$c(a_{j^{\prime}},a_{j},{\,|\,}\bm{x},v)<0$, we assume the unmeasured factors
tend to lead clinicians to systematically prescribe $a_{j}$ to healthier
patients relative to $a_{j^{\prime}}$, because patients treated with $a_{j}$
will on average have longer potential survival time to both $a_{j}$ and
$a_{j^{\prime}}$ than patients treated with $a_{j^{\prime}}$. Web Table 1
presents interpretations of confounding functions with other specifications of
signs. When setting the bounds, we assume the unmeasured confounding would
account for less than $\omega$ units of the remaining standard deviation
unexplained by measured confounders $\bm{X}_{ik}$. Using the NCDB data in our
case study as an example, a specification of
$c(1,2{\,|\,}\bm{x},v)\leq\hat{\sigma}=0.90$ and
$c(2,1{\,|\,}\bm{x},v)\geq-\hat{\sigma}=-0.90$ assumes that patients assigned
to RP will on average have $\exp(0.90)=2.45$ months longer potential survival
times than patients assigned to EBRT+AD to both treatment options; and
therefore clinicians tend to prescribe RP to healtier patients. This bound of
$\omega=1$ unit of remaining standard deviation is a plausible assumption. In
the NCDB data, the median survival time was 7.7 and 7.8 years for EBRT+AD and
RP group respectively and was not reached for the EBRT+brachy$\pm$AD treatment
group (Web Figure 1). As suggested by Hu et al.12 we draw the values of
confounding functions from the uniform distribution.

##  4 Simulation

###  4.1 Comparison methods

Through a contextualized simulation, we investigate the practical operating
characteristics of our proposed method riAFT-BART. We also adapt 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). In
addition, we consider another outcome modeling based method, random-intercept
generalized additive proportional hazards model (riGAPH), 36, 3 that is
flexible at capturing nonlinear relationships. We use the counterfactual
survival curve as the basis to objectively compare methods. Note that we can
derive the individual survival curve corresponding to our riAFT-BART model
(LABEL:eq:_riAFT-BART) as

$\displaystyle
P\left(T_{ik}>t{\,|\,}A_{ik},\bm{X}_{ik},\sigma,b_{k}\right)=1-\Phi\left(\frac{\log
t-f(A_{ik},\bm{X}_{ik})-b_{k}}{\sigma}\right).$ (7)

We can define the causal estimands by contrasting the conditional survival
probability up to a fixed time $t^{*}$, or by comparing the conditional
restricted mean survival time (RMST). 37 For some arbitrary time bound
$t^{*}$, the RMST can be represented as the area under the survival curve
$S(t)$ up to $t^{*}$, $\text{RMST}=\int_{0}^{t^{*}}S(t)dt$. In our simulation,
we present results based on both metrics, because in our motivating prostate
cancer research question, 5-year survival and RMST are of most clinical
relevance. 1

For the DR-riAH method, following suggestions by Li et al.,30 we obtain the
DR-riAH treatment effect estimator based on the survival probability at
$t^{*}$ as

$\widehat{CATE}_{a_{j},a_{j^{\prime}}}=\frac{1}{n_{k}K}\sum_{k=1}^{K}\sum_{i=1}^{n_{k}}\left(\frac{A_{ik}P(T_{ik}>t^{*})-(A_{ik}-\hat{e}_{ik})P\left(T_{ik}(a_{j})>t^{*}\right)}{\hat{e}_{ik}}-\frac{(1-A_{ik})P(T_{ik}>t^{*})-(A_{ik}+\hat{e}_{ik})P\left(T_{ik}(a_{j^{\prime}})>t^{*}\right)}{1-\hat{e}_{ik}}\right),$

where
$\hat{e}_{ik}=P\left(A_{ik}=a_{m}{\,|\,}\bm{X}_{ik}=\bm{x},V_{k}=v\right)$ is
the estimated generalized propensity score, the observed survival probability
$P(T_{ik}>t^{*})$ is the Kaplan-Meier estimator and the predicted
counterfactual survival probability
$P\left(T_{ik}(a_{j^{\prime}})>t^{*}\right)$ is calculated from the random-
intercept additive hazards model. 38 For CATE effects based on RMST, we
replace the observed (counterfactual) survival probability with the area under
the observed (counterfactual) survival curve.

To assess the performance of each method, we compare the relative bias defined
as

$\frac{\widehat{CATE}_{a_{j},a_{j^{\prime}}}-CATE^{0}_{a_{j},a_{j^{\prime}}}}{CATE^{0}_{a_{j},a_{j^{\prime}}}},$

where $CATE^{0}_{a_{j},a_{j^{\prime}}}$ is the true treatment effect, and the
frequentist coverage probability for $\widehat{CATE}_{a_{j},a_{j^{\prime}}}$
in terms of 5-year survival probability and 5-year RMST among 250 data
replications.

When implementing the comparison methods, for weighting based methods IPW-
riCox and DR-riAH, we used Super Leaner39 to estimate the stabilized inverse
probability of treatment weights for improved modeling flexibility and
accuracy of the estimated weights. Super Learner was implemented via the R
package SuperLearner with SL.library = c("SL.xgboost",
"SL.bartMachine","SL.gbm"). We fitted a weighted random-intercept Cox
regression model using the coxme function of R  package coxme to obtain the
IPW-riCox estimator. For DR-riAH, we fitted the random-intercept additive
hazards model using the aalen function from R  package timereg to compute the
counterfactual outcomes used in the DR-riAH estimator. To implement ri-GAPH,
we used the gam function from R  package mgcv and two helper functions (as_ped
and add_surv_prob) from R  package pammtools. RMST was calculated by the
trapezoidal rule using the R  function rmst of RISCA package. For all methods,
the same confounders available to the analyst were used in the linear forms in
the corresponding models.

###  4.2 Simulation design

Our data generating processes are contextualized in NCDB data settings. We
generate $K=20$ clusters, each with a sample size of $n_{k}=500$, and the
total sample size is $N=10000$. We simulate 10 confounding variables, with
five continuous variables independently generated from the standard normal
distribution $X_{ikj}\sim N(0,1),\;j=1,\ldots,5$, two categorical variables
independently generated from the multinomial distribution
$X_{ikj}\sim\text{Multinomial}(1,.3,.3,.4),\;j=6,7$ and three binary variables
$X_{ik8}\sim\text{Bern}(0.6)$, $X_{ik9}\sim\text{Bern}(0.4)$,
$X_{ik10}\sim\text{Bern}(0.5)$ generated for each individual $i$ in cluster
$k$. Throughout we consider three treatment groups. The treatment assignment
mechanism follows a random intercept multinomial logistic regression model,

$\begin{split}\ln\dfrac{P(A_{ik}=1)}{P(A_{ik}=3)}&=\xi_{01}+\bm{X}_{ik}\xi_{1}^{L}+\bm{G}_{ik}\xi_{1}^{NL}+\tau_{k}\\\
\ln\dfrac{P(A_{ik}=2)}{P(A_{ik}=3)}&=\xi_{02}+\bm{X}_{ik}\xi_{2}^{L}+\bm{G}_{ik}\xi_{2}^{NL}+\tau_{k}\\\
\end{split}$ (8)

where $\tau_{k}\sim N(0,1^{2})$, $\bm{G}_{ik}$ denotes the nonlinear
transformations and higher-order terms of the predictors $\bm{X}_{ik}$, and
$\xi^{L}_{1},\xi^{L}_{2}$ and $\xi^{NL}_{1},\xi^{NL}_{2}$ are respectively
vectors of coefficients for the untransformed versions of the confounders
$\bm{X}_{ik}$ and for the transformed versions of the confounders captured in
$\bm{G}_{ik}$. The intercepts $\xi_{01}$ and $\xi_{02}$ control the ratio of
units across three treatment groups, for which we use 6:3:1 to mimic the ratio
of individuals in the NCDB data.

We generate the potential survival times from a Weibull survival curve,

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

where $\beta^{L}_{a_{j}}$ is a treatment-specific vector of coefficients for
$\bm{X}_{ik}$ and $\beta^{NL}_{a_{j}}$ for $\bm{G}_{ik}$, $\forall
a_{j}\in\\{1,2,3\\}$. We explicitly avoid generating survival times from a
lognormal AFT model to assess the robustness of the assumption of lognormal
residuals in our riAFT-BART model formulation in equation (LABEL:eq:_riAFT-
BART). The parameter $\eta$ is set to 2 and $\exp(0.7+0.5x_{1})$ to
respectively produce proportional hazards (PH) and nonproportional hazards
(nPH). Three sets of non-parallel response surfaces are generated as

$\displaystyle T_{ik}\left(a_{j}\right)=\left\\{\frac{-\log
U}{\lambda_{a_{j}}\exp\left(\bm{X}_{ik}\beta^{L}_{a_{j}}+\bm{G}_{ik}\beta^{NL}_{a_{j}}+b_{k}\right)}\right\\}^{1/\eta}$
(10)

for $a_{j}\in\\{1,2,3\\}$, where $U$ is a random variable following the
uniform distribution on the interval $[0,1]$, $b_{k}\sim
N(0,4^{2}),\lambda_{a_{j}}=\\{3000,1200,2000\\}$ for $a_{j}=1,2,3$. Observed
and uncensored survival times are generated as
$T_{ik}=\sum_{a_{j}\in\\{1,2,3\\}}T_{ik}(a_{j})I(A_{ik}=a_{j})$. We further
generate censoring time $C$ independently from an exponential distribution
with the rate parameter selected to induce two different censoring
proportions: 10% and 40%.

The data generating processes will produce four configurations: (PH vs. nPH)
$\times$ (10% censoring proportion vs. 40% censoring proportion). Detailed
model specifications for treatment assignment, as in model (8), and for
outcomes, as in model (10), are given in Web Table 2. Web Figure 2 presents
the Kaplan-Meier survival curves stratified by treatment for each of four data
configurations in our simulation. The assessment of covariate overlap
displayed in Web Figure 3 suggests there is moderate to strong overlap across
three simulated treatment groups, which represents the overlap in the NCDB
dataset (Web Figure 4).

We conduct an illustrative simulation, using one individual-level binary
measured confounder and one individual-level binary unmeasured confounder, to
examine how our sensitivity analysis approach performs in comparison to two
causal analyses: (i) including unmeasured confounders and (ii) ignoring
unmeasured confounders. Simulation design details are provided in Web Section
S3.

Figure 1: Relative biases among 250 replications for each of four methods,
IPW-riCox, DR-riAH, riGAPH and riAFT-BART, and three treatment effects
$CATE_{1,2}$, $CATE_{1,3}$ and $CATE_{2,3}$ based on 5-year RMST under four
data configurations: (proportional hazards vs. nonproporitonal hazards)
$\times$ (10% censoring proportion vs. 40% censoring proportion). The true
treatment effects under proportional hazards are $CATE^{0,PH}_{1,2}=7.7$
months, $CATE^{0,PH}_{1,3}=3.6$ months and $CATE^{0,PH}_{2,3}=-4.1$ months.
The true treatment effects under nonproportional hazards are
$CATE^{0,nPH}_{1,2}=8.1$ months, $CATE^{0,nPH}_{1,3}=3.9$ months and
$CATE^{0,nPH}_{2,3}=-4.2$ months. Figure 2: Relative biases in the estimates
of three pairwise treatment effects $\text{CATE}_{1,2}$, $\text{CATE}_{1,3}$
and $\text{CATE}_{2,3}$ among 1000 replications using data simulated for
illustrative sensitivity analysis. Three causal analyses were performed: (i)
including unmeasured confounder (Including UMC), (ii) sensitivity analysis
given known confounding functions (SA), and (iii) ignoring unmeasured
confounder (Ignoring UMC).

###  4.3 Results

Figure 1 displays boxplots of relative biases in three treatment effect
estimates $\widehat{CATE}_{1,2}$, $\widehat{CATE}_{1,3}$ and
$\widehat{CATE}_{2,3}$ based on 5-year RMST, among 250 simulations under four
data configurations. Our proposed method riAFT-BART boasts the smallest biases
and variability in all treatment effect estimates across all simulation
scenarios, followed by DR-riAH and riGAPH; while IPW-riCox yields the largest
biases and variability. When the censoring proportion increases, all four
methods show decreased performance. The violation of proportional hazards has
the largest impact on the IPW-riCox method, demonstrated by the largest bias
increase in CATE estiamtes, but only a small impact on riAFT-BART. Even though
our method does not require proportional hazards, the elevated data complexity
in the non-proportional hazards setting – the shape parameter in the outcome
model (10) is covariate-dependent – may have contributed to the slight
increase in the biases. The bias assessment based on 5-year survival
probability is provided in Web Figure 5, which conveys the same messages as
the RMST based results.

Table 1: The coverage probability for three treatment effect estimates $\widehat{CATE}_{1,2}$, $\widehat{CATE}_{1,3}$ and $\widehat{CATE}_{2,3}$ based on 5-year RMST under four data configurations: (proportional hazards vs. nonproporitonal hazards) $\times$ (10% censoring proportion vs. 40% censoring proportion). |  | Proportional hazards |  | Nonproportional hazards
---|---|---|---|---
Censoring % | Methods | $\text{CATE}_{1,2}$ | $\text{CATE}_{1,3}$ | $\text{CATE}_{2,3}$ |  | $\text{CATE}_{1,2}$ | $\text{CATE}_{1,3}$ | $\text{CATE}_{2,3}$
10% | IPW-riCox | 25.2 | 29.2 | 28.8 |  | 17.2 | 22.4 | 19.6
DR-riAH | 81.6 | 81.2 | 82.0 |  | 79.2 | 79.6 | 79.6
ri-GAPH | 84.4 | 85.2 | 84.8 |  | 82.4 | 84.0 | 83.2
riAFT-BART | 94.8 | 95.2 | 95.2 |  | 94.4 | 94.4 | 94.8
40% | IPW-riCox | 22.0 | 26.4 | 24.0 |  | 13.2 | 19.6 | 16.8
DR-riAH | 78.4 | 79.6 | 79.6 |  | 76.2 | 78.8 | 77.2
ri-GAPH | 81.2 | 83.2 | 82.0 |  | 79.2 | 82.0 | 80.4
riAFT-BART | 94.8 | 95.2 | 94.8 |  | 94.0 | 93.2 | 94.0

Table 1 presents the frequentist coverage probability of each of four
estiamtors for each simulation configuration. The proposed method riAFT-BART
provides nominal frequentist coverage probability under proportional hazards
with 10% censoring proportion. Even in the most complex data settings with
nonproportional hazards and 40% censoring proportion, riAFT-BART still
provides close-to-nominal frequetist coverage probability. By comparison, the
IPW-riCox estimator is the least efficient producing unsatisfactory coverage
probability; DR-riAH and ri-GAPH deliver similar coverage probabilities around
0.8. Web Table 3 examines the frequentist coverage probability of the
estimators based on 5-year survival probability, and we observed the same
differences across the estimators demonstrated in the RMST based results.

Web Figure 6 suggests that our riAFT-BART algorithm converges well by plotting
3500 posterior draws of the variance parameters $\tau$ and $\sigma$, and
cluster-specific parameter $\alpha_{k}$ and the random intercepts $b_{k}$ for
clusters $k=1$, $k=10$ and $k=20$.

The illustrative simulation for sensitivity analysis, displayed in Figure 2,
empirically supports our proposed sensitivity analysis Algorithm 2. Given the
known confounding functions, our sensitivity analysis estimators are similar
to the results that could be achieved had the unmeasured confounders been made
available to the analyst. The naive analysis where we ignored the unmeasured
confounding produced substantially biased estimators.

##  5 Application to prostate cancer

We applied the proposed method riAFT-BART to estimate the comparative
effectiveness of three treatment approaches, RP, EBRT+AD and
EBRT+brachy$\pm$AD, on patient survival among high-risk localized prostate
cancer patients who are older than 65 years of age. We then applied the
proposed sensitivity analysis approach to evaluate how the causal conclusions
about treatment effects would change in response to various degrees of
departure from the no unmeasured confounding assumption.

The analysis dataset was drawn from the NCDB and included 23058 high-risk
localized prostate cancer patients who were older than 65 when diagnosed
between 2004 and 2015. Among these patients, 14237 received RP, 6683 undertook
EBRT+AD with at least 7920 cGy EBRT dose 3 and 2138 were treated with
EBRT+brachy$\pm$AD. Included in the dataset are pre-treatment patient
information on age, race and ethnicity, insurance status, income, education
level, clinical T stage, year of diagnosis, prostate-specific antigen and
gleason score, and geographic locations of treating facilities. There are nine
hospital locations, which were considered as the clusters in our analysis.
Detailed descriptions of the individual- and cluster-level variables (hospital
locations) are presented in Web Table 4. Covariates are deemed to have good
overlap across three treatment groups based on the estimated generalized
propensity scores shown in Web Figure 4.

Under the assumption of no unmeasured confounding, the treatment effect
estimates, shown in Table 2, suggest that RP is the most beneficial treatment,
which would on average lead to a ratio of 1.5 (1.3, 1.7) in the expected
survival time compared to EBRT+AD and a ratio of 1.2 (1.1, 1.4) compared to
EBRT+brachy$\pm$AD. Between the two radiotherapy approaches, EBRT +AD leads to
a shorter expected survival time that is 0.9 (0.7, 1) times the expected
survival time for EBRT+brachy$\pm$AD. Figure 3presents the posterior mean of
the predicted counterfactual survival curves for each of three treatment
groups. After correcting for confounding and accounting for the variability in
location effects, our riAFT-BART estimators suggest increased treatment
benefit associated with RP and reduced survival for EBRT+AD compared to the
unadjusted Kaplan-Meier estimators. A perusal of the posterior distribution of
the random intercepts $b_{k}$’s displayed in Figure 4 suggests that there was
substantial variability in the location effect. Hospitals in New England had
significantly better outcomes (longer expected survival times) than hospitals
in East Central area.

Figure 3: The posterior mean of the counterfactual survival curves for each of
three treatment groups in NCDB data. The solid curves are the average by
treatment group of the individual-specific survival curves estimated following
equation (7). The dashed survival curves are the Kaplan-Meier estimates for
each treatment group. Solid gray curves are estimates of individual-specific
survival curves for 15 randomly selected patients from three treatment groups.
Figure 4: The location effect in terms of the log survival time in months
represented by the posterior mean and credible intervals of the random
intercept $b_{k}$, $k=1,\ldots,9$.

It is possible that some important confounders were not collected in the NCDB
data. For example, patient functional status has been shown to be strongly
associated with both treatment choices and survival among men with prostate
cancer. 40, 41 The performance status measured by Eastern Clinical Oncology
Group (ECOG) score may also be a confounder as it both affects the likelihood
of clinicians choosing AD 42 and predicts survival for prostate cancer. 43
Futhermore, magnetic resonance imaging findings and the number of positive
biopsy cores are both likely confounders as they are related to patient
selection for RP or brachytherapy and indicate the degree of agressiveness of
high-risk prostate cancer. 2, 44

To evaluate the sensitivity of the treatment effect estimates to these
unmeasured confounders, we first leverage the subject-area literature to
specify the confounding functions. For the sign, we assume that the unmeasured
factors guiding clinicians to prescribe RP lead them systematically to
prescribe it to relatively healthier patients. This is because magnetic
resonance imaging findings supportive of resectability were used for patient
selection for RP, 2 and RP was recommended to those with lower number of
positive biopsy cores, better functional score and better performance status.
40, 45 Between the two radiotherapy based treatment approaches, on the one
hand, unhealthier patients with lower functional scores or ECOG scores were
not recommended to use AD in treatment as they would not tolerate strong side
effects induced by AD. 42, 46 On the other hand, brachytheray has been
recommended as a boost to EBRT only to relatively unhealthier patients with
multiple positive biopsy cores.44 Based on these pieces of evidence, we assume
clinicians may have a preference to recommend EBRT+AD over EBRT+brachy$\pm$AD
to healthier patients, but other directions of unmeasured confounding may also
be plausible.

We next postulate the bounds of confounding functions based on published
scientific work. Lehtonen et al. 43 shows that ECOG score has a large effect
(Cohen’s $d>1$) and number of positive biopsy cores has a small effect
(Cohen’s $d<0.2$) on prostate cancer patient survival. Loosely translating
Cohen’s $d$ to the proportion of the total variation in the outcome explained
by a given covariate, 47 we assume the unmeasured confounders approximately
account for 25% of the total variance in overall survival. Following Chan et
al., 48 we fitted an AFT log-normal model to the NCDB data with fixed effects
for locations, and computed the $R^{2}$ to be 60%. This suggests that the
unmeasured confounders will explain no more than 40% of the variation in the
outcomes. Based on these grounds, we assume that the unmeasured confounding
would account for $\omega=0.75$ units of the remaining standard deviation
unexplained by measured variables. The remaining standard deviation in the
outcome unexplained by measured covariates was estimated to be $\hat{\sigma}$
= 0.90 (months) via our riAFT-BART model. We hence set the bounds of the
confounding functions to be $\pm\omega\hat{\sigma}=\pm.675$ (months).

Table 2 displays the sensitivity analysis results in comparison to treatment
effect estimates under the weak unconfoundedness assumption. We assume
relatively healthier patients were assigned to RP:
$c(1,2{\,|\,}\bm{x},v)\sim\mathcal{U}(0,.675),\;c(2,1{\,|\,}\bm{x},v)\sim\mathcal{U}(-.675,0)$,
and
$c(1,3{\,|\,}\bm{x},v)\sim\mathcal{U}(0,.675),\;c(3,1{\,|\,}\bm{x},v)\sim\mathcal{U}(-.675,0)$.
Between EBRT+AD and EBRT+brachy$\pm$AD, all four different forms of
confounding functions are deemed plausible (see Web Table 1 for
interpretations). Results show that the significant treatment benefit
associated with RP over EBRT+AD is robust to different magnitudes and
directions of unmeasured confounding. However, the significant gain in the
expected survival time offered by RP in the causal analysis assuming weak
unconfoundedness is negated in the presence of unmeasured confounding. Turning
to the comparative effect between EBRT+AD and EBRT+brachy$\pm$AD, under the
assumption of no unmeasured confounding, EBRT+brachy$\pm$AD had a survival
benefit bordering on being statistically significant over EBRT +AD. Assuming
unmeasured factors guiding clinicians to prescribe EBRT+AD lead them
systematically to prescribe it to relatively healthier patients “tips” the
result over to significant benefit for EBRT+brachy$\pm$AD.

Table 2: Estimation of causal effects of three treatment approaches on patient survival for high-risk localized prostate cancer, based on differences in log survival months, using the NCDB database. Three treatment options are $A=1$: RP, $A=2$: EBRT+AD and $A=3$: EBRT+brachy$\pm$AD. Row 2–5 present confounding function adjusted effect estimates using our sensitivity analysis Algorithm 2. Interval estimates are based on pooled posterior samples across model fits arising from $30\times 30$ data sets. We assume relatively healthier patients were assigned to RP, $c(1,2{\,|\,}\bm{x},v)>0$, $c(2,1{\,|\,}\bm{x},v)<0$; $c(1,3{\,|\,}\bm{x},v)>0$, $c(3,1{\,|\,}\bm{x},v)<0$, and four unmeasured confounding directions between $A=2$: EBRT+AD and $A=3$: EBRT+brachy$\pm$AD. We drew the values of confounding functions from the uniform distributions bounded by $\pm.675$ (months). |  | $CATE_{1,2}$ | $CATE_{1,3}$ | $CATE_{2,3}$
---|---|---|---|---
| Assuming weak unconfoundedness | $.38(.24,.52)$ | $.22(.07,.37)$ | $-.16(-.32,.00)$
Sensitivity analysis | $c(2,3{\,|\,}\bm{x},v)>0,\;c(3,2{\,|\,}\bm{x},v)<0$ | $.27(.09,.45)$ | $.12(-.07,.31)$ | $-.28(-.48,-.08)$
$c(2,3{\,|\,}\bm{x},v)<0,\;c(3,2{\,|\,}\bm{x},v)>0$ | $.31(.13,.49)$ | $.13(-.06,.32)$ | $-.05(-.25,.15)$
$c(2,3{\,|\,}\bm{x},v)>0,\;c(3,2{\,|\,}\bm{x},v)>0$ | $.27(.09,.45)$ | $.10(-.09,.29)$ | $-.13(-.33,.07)$
$c(2,3{\,|\,}\bm{x},v)<0,\;c(3,2{\,|\,}\bm{x},v)<0$ | $.28(.10,.46)$ | $.11(-.08,.30)$ | $-.19(-.39,.01)$

##  6 Summary and discussion

The increased availability of large-scale healthcare databases allows
researchers to conduct comparative effectiveness analysis of modern treatment
approaches for high-risk cancer patients. Recent efforts have been made to
investigate the effects of surgical and radiotherapy based treatments on
patient survival using the national cancer databases. However, the multilevel
data structure presented in these databases and its implications for causal
analyses require special statistical consideration but have not been well
studied. In addition, there is the lack of tools for evaluating the
sensitivity of treatment effect estimates to the presence of unmeasured
confounding in the context of multiple treatments and multilevel survival
data.

Motivated by these research gaps, our work makes two primary contributions to
the causal inference literature. First, we develop a flexible causal modeling
tool and MCMC algorithm for causal inferences about effects of multiple
treatments on patient survival while respecting the multilevel data structure.
The proposed riAFT-BART model flexibly captures the relationship between
survival times and individual-level covariates and cluster-specific main
effects, while providing proper representations of uncertainty intervals via
the posterior based on a probability model. Second, leveraging the flexible
riAFT-BART model, we develop an interpretable sensitivity analysis algorithm
to address the causal assumption of no unmeasured confounding. In our
sensitivity analysis approach, we present the confounding function on the
basis of expected survival time, which can be easily interpreted, and propose
methods to adjust the estimation of causal effects by effectively removing the
bias due to posited levels of unmeasured confounding.

Applying our proposed methods to NCDB data on older high-risk localized
prostate cancer patients, we confirmed the survival benefit of RP relative to
EBRT+AD. Inferences about other two pairwise treatment effects were
inconclusive because they were impacted by the potential unmeasured
confounding. Our causal analysis also demonstrates that there is substantial
variability in the effects of hospital locations, which reinforces the
importance of examining the cluster-level variation when estimating treatment
effect in the general population using data with hierarchical structure.

There are several important avenues for future research. First, our riAFT-BART
model can be extended to include the random slopes and accommodate the
cluster-level covariates. Second, developing a sensitivity analysis for
cluster-level unmeasured confounding could be a worthwhile and important
contribution. Third, although our simulation results suggest that our methods
are robust to the normality assumption for the random intercepts and
residuals, it may be worthwhile to develop nonparametric priors to further
improve the modeling flexibility. Finally, to address the causal assumption of
positivity, we can extend the work by Hill and Su 49 and Hu et al. 8 to
develop a strategy to identify a common support region for inferential units.

## ACKNOWLEDGEMENTS

This work was supported in part by the National Cancer Institute under grant
NIH NCI R21CA245855, and by award ME_2017C3_9041 from the Patient-Centered
Outcomes Research Institute (PCORI). The content is solely the responsibility
of the authors and does not necessarily represent the official views of the
National Institutes of Health or PCORI.

## DATA AVAILABILITY STATEMENT

The simulation codes that generate the data supporting the findings of the
simulation study are openly available with the R  package riAFTBART. The NCDB
data used in the case study is publicly available upon approval of the NCDB
Participant User File application.

## References

  * 1 Ennis RD, Hu L, Ryemon SN, Lin J, Mazumdar M. Brachytherapy-based radiotherapy and radical prostatectomy are associated with similar survival in high-risk localized prostate cancer. Journal of Clinical Oncology 2018; 36(12): 1192–1198.
  * 2 Chen RC. Challenges of interpreting registry data in prostate cancer: interpreting retrospective results along with or in absence of clinical trial data. Journal of Clinical Oncology 2018; 36(12): 1181–1183.
  * 3 Hu L, Ji J, Li F. Estimating heterogeneous survival treatment effect in observational data using machine learning. Statistics in Medicine 2021; 40(21): 4691–4713.
  * 4 Zeng S, Li F, Hu L, Li F. Propensity Score Weighting Analysis of Survival Outcomes Using Pseudo-observations. Statistica Sinica 2022\.  In press.
  * 5 Feng P, Zhou XH, Zou QM, Fan MY, Li XS. Generalized propensity score for estimating the average treatment effect of multiple treatments. Statistics in Medicine 2012; 31(7): 681–697.
  * 6 McCaffrey DF, Griffin BA, Almirall D, Slaughter ME, Ramchand R, Burgette LF. A tutorial on propensity score estimation for multiple treatments using generalized boosted models. Statistics in Medicine 2013; 32(19): 3388–3414.
  * 7 Linden A, Uysal SD, Ryan A, Adams JL. Estimating causal effects for multivalued treatments: a comparison of approaches. Statistics in Medicine 2016; 35(4): 534–552.
  * 8 Hu L, Gu C, Lopez M, Ji J, Wisnivesky J. Estimation of causal effects of multiple treatments in observational studies with a binary outcome. Statistical Methods in Medical Research 2020; 29(11): 3218–3234.
  * 9 Hu L, 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.
  * 10 Yu Y, Zhang M, Shi X, Caram ME, Little RJ, Mukherjee B. A comparison of parametric propensity score-based methods for causal inference with multiple treatments and a binary outcome. Statistics in Medicine 2021; 40(7): 1653–1677.
  * 11 Von Elm E, Altman DG, Egger M, Pocock SJ, Gøtzsche PC, Vandenbroucke JP. The Strengthening the Reporting of Observational Studies in Epidemiology (STROBE) statement: guidelines for reporting observational studies. Bulletin of the World Health Organization 2007; 85(11): 867–872.
  * 12 Hu L, Zou J, Gu C, Ji J, Lopez M, Kale M. A flexible sensitivity analysis approach for unmeasured confounding with multiple treatments and a binary outcome with application to SEER-Medicare lung cancer data. The Annals of Applied Statistics 2022\.  In press.
  * 13 Chipman HA, George EI, McCulloch RE. BART: Bayesian additive regression trees. The Annals of Applied Statistics 2010; 4(1): 266–298.
  * 14 Hill JL. Bayesian nonparametric modeling for causal inference. Journal of Computational and Graphical Statistics 2011; 20(1): 217–240.
  * 15 Hu L, Lin J, 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.
  * 16 Hu L, Lin JY, Sigel K, Kale M. Estimating heterogeneous survival treatment effects of lung cancer screening approaches: A causal machine learning analysis. Annals of Epidemiology 2021; 62: 36–42.
  * 17 Chen PY, Tsiatis AA. Causal inference on the difference of the restricted mean lifetime between two groups. Biometrics 2001; 57(4): 1030–1038.
  * 18 Arpino B, 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 Hernán MA, Robins JM. Estimating causal effects from epidemiological data. Journal of Epidemiology & Community Health 2006; 60(7): 578–586.
  * 20 Hernán MA, Robins JM. Causal Inference: What If.  Boca Raton: Chapman & Hall/CRC .  2020\.
  * 21 Tan YV, Roy J. Bayesian additive regression trees and the General BART model. Statistics in Medicine 2019; 38(25): 5048–5069.
  * 22 Gelman A, Van Dyk DA, Huang Z, Boscardin JW. Using redundant parameterizations to fit hierarchical models. Journal of Computational and Graphical Statistics 2008; 17(1): 95–122.
  * 23 Henderson NC, Louis TA, Rosner GL, Varadhan R. Individualized treatment effects with censored data via fully nonparametric Bayesian accelerated failure time models. Biostatistics 2020; 21(1): 50–68.
  * 24 Roy J, Lum KJ, Daniels MJ. A Bayesian nonparametric approach to marginal structural models for point treatments and a continuous or survival outcome. Biostatistics 2017; 18(1): 32–47.
  * 25 Rosenbaum PR. Covariance adjustment in randomized experiments and observational studies. Statistical Science 2002; 17(3): 286–327.
  * 26 Kasza J, Wolfe R, Schuster T. Assessing the impact of unmeasured confounding for binary outcomes using confounding functions. International Journal of Epidemiology 2017; 46(4): 1303–1311.
  * 27 Robins JM. Association, causation, and marginal structural models. Synthese 1999: 151–179.
  * 28 VanderWeele TJ, Ding P. Sensitivity analysis in observational research: introducing the E-value. Annals of Internal Medicine 2017; 167(4): 268–274.
  * 29 Arpino B, Mealli F. The specification of the propensity score in multilevel observational studies. Computational Statistics & Data Analysis 2011; 55(4): 1770–1780.
  * 30 Li F, Zaslavsky AM, Landrum MB. Propensity score weighting with multilevel data. Statistics in Medicine 2013; 32(19): 3373–3387.
  * 31 Fuentes A, Lüdtke O, Robitzsch A. Causal Inference with Multilevel Data: A Comparison of Different Propensity Score Weighting Approaches. Multivariate Behavioral Research 2021\.  In press.
  * 32 Brumback BA, Hernán MA, Haneuse SJ, Robins JM. Sensitivity analyses for unmeasured confounding assuming a marginal structural model for repeated measures. Statistics in Medicine 2004; 23(5): 749–767.
  * 33 Rubin DB. Nested multiple imputation of NMES via partially incompatible MCMC. Statistica Neerlandica 2003; 57(1): 3–18.
  * 34 Zhou X, Reiter JP. A note on Bayesian inference after multiple imputation. The American Statistician 2010; 64(2): 159–163.
  * 35 Hogan JW, Daniels MJ, Hu L. A Bayesian perspective on assessing sensitivity to assumptions about unobserved data. In: Molenberghs G, Fitzmaurice G, Kenward MG, Tsiatis A, Verbeke G. , eds. Handbook of Missing Data MethodologyBoca Raton, FL: CRC Press. 2014 (pp. 405–434).
  * 36 Hastie TJ, Tibshirani RJ. Generalized additive models.  Boca Raton, FL: Chapman &Hall .  1990\.
  * 37 Royston P, 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.
  * 38 Cai J, Zeng D. Additive mixed effect model for clustered failure time data. Biometrics 2011; 67(4): 1340–1351.
  * 39 Laan V. dMJ, Polley EC, Hubbard AE. Super learner. Statistical Applications in Genetics and Molecular Biology 2007; 6(1).
  * 40 Jacobs BL, Lopa SH, Yabes JG, Nelson JB, Barnato AE, Degenholtz HB. Association of functional status and treatment choice among older men with prostate cancer in the Medicare Advantage population. Cancer 2016; 122(20): 3199–3206.
  * 41 Stommel M, Given BA, Given CW. Depression and functional status as predictors of death among cancer patients. Cancer 2002; 94(10): 2719–2727.
  * 42 Varenhorst E, Klaff R, Berglund A, Hedlund PO, Sandblom G, 5 SPCGSTN. Predictors of early androgen deprivation treatment failure in prostate cancer with bone metastases. Cancer Medicine 2016; 5(3): 407–414.
  * 43 Lehtonen M, Heiskanen L, Reinikainen P, Kellokumpu-Lehtinen PL. Both comorbidity and worse performance status are associated with poorer overall survival after external beam radiotherapy for prostate cancer. BMC Cancer 2020; 20(1): 1–8.
  * 44 Nag S, Beyer D, Friedland J, Grimm P, Nath R. American Brachytherapy Society (ABS) recommendations for transperineal permanent brachytherapy of prostate cancer. International Journal of Radiation Oncology, Biology, Physics 1999; 44(4): 789–799.
  * 45 Lee A, Schultz D, Renshaw A, Richie J, D’Amico A. Optimizing patient selection for prostate monotherapy. International Journal of Radiation Oncology, Biology, Physics 2000; 3(48): 306–307.
  * 46 Clay CA, Perera S, Wagner JM, Miller ME, Nelson JB, Greenspan SL. Physical function in men with prostate cancer on androgen deprivation therapy. Physical Therapy 2007; 87(10): 1325–1333.
  * 47 Lakens D. Calculating and reporting effect sizes to facilitate cumulative science: a practical primer for t-tests and ANOVAs. Frontiers in psychology 2013; 4: 863.
  * 48 Chan PH, Xu R, Chambers CD. A study of $R^{2}$ measure under the accelerated failure time models. Communications in Statistics-Simulation and Computation 2018; 47(2): 380–391.
  * 49 Hill J, Su YS. Assessing lack of common support in causal inference using Bayesian nonparametrics: Implications for evaluating the effect of breastfeeding on children’s cognitive outcomes. The Annals of Applied Statistics 2013; 7(3): 1386–1420.