Patient stratification in multi-arm trials: a two-stage procedure with Bayesian profile regression

TL;DR

This paper introduces a novel two-stage Bayesian nonparametric method for patient stratification in multi-arm clinical trials, combining outcome prediction and subgroup clustering to improve personalized treatment analysis.

Contribution

It develops a new Bayesian framework that integrates outcome prediction with patient clustering and variable selection specifically for multi-arm trial settings.

Findings

Successfully identified five meaningful donor subgroups in a blood donation trial.

Demonstrated superior performance over existing methods in simulation studies.

Provides a flexible approach adaptable to various multi-arm clinical trial designs.

Abstract

Precision medicine is an emerging field that takes into account individual heterogeneity to inform better clinical practice. In clinical trials, the evaluation of treatment effect heterogeneity is an important component, and recently, many statistical methods have been proposed for stratifying patients into different subgroups based on such heterogeneity. However, the majority of existing methods developed for this purpose focus on the case with a dichotomous treatment and are not directly applicable to multi-arm trials. In this paper, we consider the problem of patient stratification in multi-arm trial settings and propose a two-stage procedure within the Bayesian nonparametric framework. Specifically, we first use Bayesian additive regression trees (BART) to predict potential outcomes (treatment responses) under different treatment options for each patient, and then we leverage Bayesian profile regression to cluster patients into subgroups according to their baseline characteristics and predicted potential outcomes. We further embed a variable selection procedure into our proposed framework to identify the patient characteristics that actively "drive" the clustering structure. We conduct simulation studies to examine the performance of our proposed method and demonstrate the method by applying it to a UK-based multi-arm blood donation trial, wherein our method uncovers five clinically meaningful donor subgroups.

Tables

Table 1. (a) Scenario 1

	cluster 1	cluster 2	cluster 3
$E\left(X_1\right)$	2	4	6
$E\left(X_2\right)$	4	6	1
$E\left(X_3\right)$	5	1	3
$E\left\{Y^{\ast }\left(1\right)\right\}$	4	4	4
$E\left\{Y^{\ast }\left(2\right)\right\}$	4	7	6
$E\left\{Y^{\ast }\left(3\right)\right\}$	4	9	5

Table 2. (a) Scenario 1

	cluster 1	cluster 2	cluster 3
$E\left(X_1\right)$	2	4	6
$E\left(X_2\right)$	4	6	1
$E\left(X_3\right)$	5	1	3
$E\left\{Y^{\ast }\left(1\right)\right\}$	4	4	4
$E\left\{Y^{\ast }\left(2\right)\right\}$	4	7	6
$E\left\{Y^{\ast }\left(3\right)\right\}$	4	9	5

Table 3. (b) Scenario 2

	cluster 1	cluster 2	cluster 3	cluster 4
$P\left(X_1=1\right)$	0.2	0.4	0.6	0.8
$P\left(X_2=1\right)$	0.4	0.4	0.6	0.6
$P\left(X_3=0\right)$	0.1	0.2	0.3	0.4
$P\left(X_3=1\right)$	0.15	0.3	0.15	0.3
$E\left(X_4\right)$	2	4	8	6
$E\left\{Y^{\ast }\left(1\right)\right\}$	2	2	2	3
$E\left\{Y^{\ast }\left(2\right)\right\}$	2	5	4	6
$E\left\{Y^{\ast }\left(3\right)\right\}$	2	4	5	8
$E\left\{Y^{\ast }\left(4\right)\right\}$	2	3	6	8

Equations

$$p({\mathbf{Y}}^{\ast },X|\mathbf{\pi },\mathbf{\Theta },\mathbf{\Phi })=\sum _{c=1}^{\infty }{\pi }_{c}f({\mathbf{Y}}^{\ast }|{\Theta }_c)f(X|{\Phi }_c),$$

$$\begin{array}{rl}{V}_c& \sim \text{Beta}(1,\alpha )\alpha >0,\\ {\pi }_1& =V_1,{\pi }_c=V_c\prod _{r=1}^{c-1}(1-V_r)\text{for }c\ge 2,\end{array}$$

$$f(X|{\Phi }_c)=f({\mathbf{X}}^{\text{cont}}|{\mathbf{\mu }}_{(X)c},{\Sigma }_{(X)c})f({\mathbf{X}}^{\text{disc}}|{\mathbf{\Psi }}_c),$$

$$\begin{split}f(\mathbf{X}^{\text{cont}}\big{|}\,\bm{\mu}_{(\operatorname*{\mathbf{X}})\,c},\Sigma_{(\operatorname*{\mathbf{X}})\,c})=\Big{\{}(2\pi)^{{p}_{1}}\big{|}\Sigma_{(\operatorname*{\mathbf{X}})\,c}\big{|}\Big{\}}^{-\frac{1}{2}}\exp{\bigg{\{}-\frac{1}{2}\big{(}\mathbf{X}^{\text{cont}}-\bm{\mu}_{(\operatorname*{\mathbf{X}})\,c}\big{)}^{\!\mathsf{T}}\Sigma^{-1}_{(\operatorname*{\mathbf{X}})\,c}\big{(}\mathbf{X}^{\text{cont}}-\bm{\mu}_{(\operatorname*{\mathbf{X}})\,c}\big{)}\bigg{\}}}.\end{split}$$

$$\begin{split}\bm{\mu}_{(\operatorname*{\mathbf{X}})\,c}\,|\,\Sigma_{(\operatorname*{\mathbf{X}})\,c}&\sim\mathcal{N}(\bm{\mu}_{(\mathbf{X})},\Sigma_{(\operatorname*{\mathbf{X}})\,c}/\kappa_{(\mathbf{X})}),\\ \Sigma_{(\operatorname*{\mathbf{X}})\,c}&\sim\mathcal{W}^{-1}(\Lambda_{(\mathbf{X})},\nu_{(\mathbf{X})}),\\ \big{(}\bm{\mu}_{(\operatorname*{\mathbf{X}})\,c},\Sigma_{(\operatorname*{\mathbf{X}})\,c}\big{)}&\sim{NIW}(\bm{\mu}_{(\mathbf{X})},\kappa_{(\mathbf{X})},\Lambda_{(\mathbf{X})},\nu_{(\mathbf{X})}),\end{split}$$

$$f({\mathbf{X}}^{\text{disc}}|{\mathbf{\Psi }}_c)=\prod _{j=p_1+1}^{p_1+p_2}{\psi }_{c,j,X_{j,}}$$

$${\Psi }_{c,j}=({\psi }_{c,j,1},{\psi }_{c,j,2},\dots ,{\psi }_{c,j,K_j})\sim \text{Dirichlet}({\mathbf{a}}_j),$$

$$f(\operatorname*{\mathbf{Y^{*}}}\big{|}\,\bm{\mu}_{(\mathbf{Y^{*}})\,c},\Sigma_{(\mathbf{Y^{*}})\,c})=\Big{\{}(2\pi)^{K}\big{|}\Sigma_{(\mathbf{Y^{*}})\,c}\big{|}\Big{\}}^{-\frac{1}{2}}\exp{\bigg{\{}-\frac{1}{2}\big{(}\operatorname*{\mathbf{Y^{*}}}-\bm{\mu}_{(\mathbf{Y^{*}})\,c}\big{)}^{\!\mathsf{T}}\Sigma^{-1}_{(\operatorname*{\mathbf{Y^{*}}})\,c}\big{(}\operatorname*{\mathbf{Y^{*}}}-\bm{\mu}_{(\mathbf{Y^{*}})\,c}\big{)}\bigg{\}}}.$$

$$\begin{split}\bm{\mu}_{(\mathbf{Y^{*}})\,c}\,|\,\Sigma_{(\mathbf{Y^{*}})\,c}&\sim\mathcal{N}(\bm{\mu}_{(\operatorname*{\mathbf{Y^{*}}})},\Sigma_{(\mathbf{Y^{*}})\,c}/\kappa_{(\operatorname*{\mathbf{Y^{*}}})}),\\ \Sigma_{(\mathbf{Y^{*}})\,c}&\sim\mathcal{W}^{-1}(\Lambda_{(\operatorname*{\mathbf{Y^{*}}})},\nu_{(\operatorname*{\mathbf{Y^{*}}})}),\\ \big{(}\bm{\mu}_{(\mathbf{Y^{*}})\,c},\Sigma_{(\mathbf{Y^{*}})\,c}\big{)}&\sim{NIW}(\bm{\mu}_{(\operatorname*{\mathbf{Y^{*}}})},\kappa_{(\operatorname*{\mathbf{Y^{*}}})},\Lambda_{(\operatorname*{\mathbf{Y^{*}}})},\nu_{(\operatorname*{\mathbf{Y^{*}}})}),\end{split}$$

$$\begin{array}{r}{\mu }_{c,j}^{\ast }={\gamma }_{c,j}{\mu }_{c,j}+(1-{\gamma }_{c,j}){\bar{x}}_j=({\mu }_{c,j}{)}^{{\gamma }_{c,j}}\times {({\bar{x}}_j)}^{(1-{\gamma }_{c,j})},\end{array}$$

$$\begin{array}{r}f({\mathbf{X}}^{\text{disc}}|{\mathbf{\Psi }}_c,{\mathbf{\gamma }}_c^{\text{disc}})=\prod _{j=p_1+1}^{p_1+p_2}{\gamma }_{c,j}{\psi }_{c,j,X_j}+(1-{\gamma }_{c,j}){\psi }_{0,j,X_j}=\prod _{j=p_1+1}^{p_1+p_2}({\psi }_{c,j,X_j}{)}^{{\gamma }_{c,j}}\times ({\psi }_{0,j,X_j}{)}^{(1-{\gamma }_{c,j})}.\end{array}$$

Full text

1 Introduction

Precision medicine is a paradigm that leverages patient heterogeneity to better inform clinical practice through data-driven, evidence-based, and scientifically-rigorous approaches, and it has gained much popularity recently partly due to the wealth of biological, clinical and epidemiological data now available. In particular, patient stratification based on patients’ treatment response profiles and baseline characteristics can greatly facilitate the understanding and characterization of the underlying patient heterogeneity, and thus can be of great interest and importance for a variety of medical applications (Sies et al., 2019).

Our research was motivated by a UK-based, multi-arm trial called INTERVAL (Moore et al., 2014), which was the first randomized trial to investigate the effect of different inter-donation intervals on blood supply and donor health. In INTERVAL, male donors were randomly assigned to 12-week (standard), 10-week, or 8-week inter-donation intervals, and female donors to 16-week (standard), 14-week, or 12-week inter-donation intervals. The outcomes of this trial include the amount of blood collected, and the number of low hemoglobin deferrals (temporary suspension of donors from giving blood) during the trial period. INTERVAL trial participants were well-characterized at baseline, providing an opportunity to investigate donor stratification and characterize different types of donors. This can be important for achieving more effective targeted recruitment of donors and increased efficiency of blood collection in the future. For example, it would be useful to stratify donors into those who have the capacity to donate more often than the current clinical practice (often referred to as “super donors” who have high donation capacity and minimal number of deferrals) and those who tend to be deferred more frequently than the average donors (often referred to as “brittle donors” who struggle more than the average donors to replace iron stores after blood donation).

A number of data-driven methods for investigating patient heterogeneity with clinical trial data have been proposed (for a review, see Sies et al. (2019)). However, to our knowledge, most of these methods were developed for two-arm trials, and they have not been extended to the situation with multi-level treatments, despite the popularity of multi-arm trials in clinical practice (Baron et al., 2013).

In this paper, we propose a method for leveraging treatment effect heterogeneity to identify patient subgroups in multi-arm trial settings from a clustering perspective. As a remark, we will use the terms “subgroup” and “cluster” interchangeably throughout this paper. In practice, the functional form of the relationship between treatment responses and covariates is typically unknown a priori . To this end, we leverage Bayesian nonparametric approaches. In particular, our proposed method is a two-stage procedure. In the first stage, we use Bayesian additive regression trees (BART) to predict each patient’s potential outcomes (treatment responses) under different treatment options. In the second stage, we employ a Bayesian mixture model called “profile regression” (Molitor et al., 2010) to cluster patients into subgroups according to their covariate profiles and predicted potential outcomes obtained in the first stage. We also extend the second stage clustering model to incorporate the variable selection feature, via which covariates that actively “drive” the clustering structure can be identified (distinguished from covariates that have similar profiles across all clusters).

Our method has several desirable properties. First, unlike conventional unsupervised clustering methods, our proposed approach allows the outcome variable to inform cluster membership, which ensures that resulting clusters are associated with the target outcome of clinical interest and are clinically meaningful. Second, our clustering model is based on the Dirichlet process prior, in which case the number of clusters can be directly inferred from the data, thus sidestepping the difficulty of pre-specifying the number of clusters in classical clustering approaches. Third, the proposed approach can handle correlated patient characteristics (models joint effects and the unit of inference is covariate profile) and correlations do not undermine model performance and interpretability. Finally, our method takes into account the uncertainty associated with the number of clusters and cluster assignments, and the uncertainty of the “representative” clustering can be quantified via model-averaging approaches (Molitor et al., 2010).

This paper is structured as follows. In Section 2, we introduce the statistical framework. The numerical performance of the method is evaluated by simulation studies in Section 3. An application of the proposed method to the INTERVAL data is presented in Section 4. We conclude the paper with a discussion in Section 5.

2 Methodology

We consider a clinical trial with $K$ treatments (randomized groups) and $n$ subjects. We let $A\in\mathcal{A}=\{1,\ldots,K\}$ denote the treatment assignment, $\operatorname*{\mathbf{X}}\in\mathcal{X}$ denote the covariate data, and $Y\in\mathbb{R}$ denote the observed outcome of interest. The potential outcome under treatment $A=a$ is denoted by $Y^{*}(a)$. Our proposed two-stage patient stratification method is detailed below in Sections 2.1 and 2.2. A graphical representation of the proposed approach is presented in Web Appendix A (Web Figure 1).

2.1 First stage: BART - predict potential outcomes

In the first stage, we fit a flexible regression model for $Y$ given $A$ and $\operatorname*{\mathbf{X}}$ in order to get the potential outcome estimates for subjects. We denote the predicted potential outcomes under treatments $1,\ldots,K$ by $\widehat{Y}^{*}(1),\ldots,\widehat{Y}^{*}(K)$, respectively. In theory, $\widehat{Y}^{*}(1),\ldots,\widehat{Y}^{*}(K)$ can be obtained from the observed data using any supervised machine learning or regression algorithms (Künzel et al., 2019). We use BART in our implementation (Hill, 2011).

We note that $\widehat{Y}^{*}(1),\ldots,\widehat{Y}^{*}(K)$ can also be obtained by fitting a separate model for $Y|\operatorname*{\mathbf{X}}$ under each treatment option (“separate-learner”) as opposed to fitting one single model for $Y|A,\operatorname*{\mathbf{X}}$ using the data from all subjects (“single-learner”) (Künzel et al., 2019). These two ways of predicting potential outcomes yield similar results in our numerical studies, and we only present the results based on the “single-learner” throughout this paper.

2.2 Second stage: Bayesian profile regression - clustering

Molitor et al. (2010) proposed a Bayesian mixture model called “profile regression”, which links the outcome vector to a set of possibly correlated covariates nonparametrically through cluster membership. Profile regression clusters subjects into subgroups by leveraging the Dirichlet process mixture model (DPMM). In profile regression, the covariates and the outcome are modeled jointly, and thus both influence cluster allocations. This is appealing, since covariate data typically give rise to many different clustering structures, and the inclusion of the target outcome information in clustering can be particularly helpful for ensuring that the clustering result has clinical utility with regard to our aim (Bair, 2013).

Our second stage clustering model follows the extension of profile regression to the case with a multivariate normal outcome (Rouanet et al., 2020). We use $\operatorname*{\mathbf{Y^{*}}}=(\widehat{Y}^{*}(1),\ldots,\widehat{Y}^{*}(K))$ to denote the multivariate normal response variable in the profile regression model. The likelihood function is given by

[TABLE]

where $c$ is the index of mixture components, $\bm\pi=(\pi_{1},\pi_{2},\ldots)$ is a vector of mixture weights, $\bm{\Theta}=(\Theta_{1},\Theta_{2},\ldots)$ denotes cluster-specific parameters in the density function for $\operatorname*{\mathbf{Y^{*}}}$, and $\bm{\Phi}=(\Phi_{1},\Phi_{2},\ldots)$ denotes cluster-specific parameters in the density function for $\operatorname*{\mathbf{X}}$. We note that by construction, DPMM allows infinite components. Therefore, the summation in (1) goes from $1$ to infinity.

2.2.1 Specification of model components

In the following, we describe the model for each component of (1), but we will not provide further details on the Bayesian computation aspect in this paper. We refer the interested readers to Liverani et al. (2015) and Rouanet et al. (2020) for computational details.

Mixture weights $\pi_{c}$

The mixture weights ($\pi_{c}$) are modeled according to the stick-breaking construction of the Dirichlet process (DP) as follows (Sethuraman, 1994):

[TABLE]

where $V_{1},V_{2},\ldots$ are independent random variables. Under this representation, $\alpha$ is the concentration parameter of DP, which reflects the dispersion level and controls the number of non-empty clusters implicitly (Teh, 2010; Hastie et al., 2015; Frühwirth-Schnatter and Malsiner-Walli, 2019). We follow Rouanet et al. (2020) and adopt a Gamma prior for $\alpha$.

The model for covariates $\operatorname*{\mathbf{X}}$

The profile regression model can handle both continuous and discrete covariates. In the following, we consider the case where $\operatorname*{\mathbf{X}}$ consists of $p_{1}$ continuous covariates and $p_{2}$ discrete covariates. In order to describe the complete covariate model, we assume that $p_{1},p_{2}\geq 1$ for now. The situations with $p_{1}=0$ or $p_{2}=0$ will be discussed later. We let $\mathbf{X}^{\text{cont}}$ and $\mathbf{X}^{\text{disc}}$ denote the subset of continuous and discrete covariates in $\operatorname*{\mathbf{X}}$, respectively. Without loss of generality, and for notational convenience, we assume that the first $p_{1}$ covariates in $\operatorname*{\mathbf{X}}$ are continuous, i.e. $X_{j}$ (the $j^{\text{\,th}}$ covariate in $\operatorname*{\mathbf{X}}$) is continuous for $j={1,\ldots,p_{1}}$, and $X_{j}$ is discrete for $j={p_{1}+1,\ldots,p_{1}+p_{2}}$.

We assume that $\mathbf{X}^{\text{cont}}$ and $\mathbf{X}^{\text{disc}}$ are independent conditional on the cluster assignments, and then the “density” for covariates $\operatorname*{\mathbf{X}}$ can be written as:

[TABLE]

where $\Phi_{c}=(\bm{\mu}_{(\operatorname*{\mathbf{X}})\,c},\Sigma_{(\operatorname*{\mathbf{X}})\,c},\bm{\Psi}_{c})$ represents the cluster-specific parameter set for covariates. In particular, $\bm{\mu}_{(\operatorname*{\mathbf{X}})\,c}$ and $\Sigma_{(\operatorname*{\mathbf{X}})\,c}$ are the mean vector and the covariance matrix for continuous covariates in cluster $c$, and $\bm{\Psi}_{c}$ denotes the parameter for discrete covariates in cluster $c$. Note that we add “$(\operatorname*{\mathbf{X}})$” in the subscripts of the parameters for continuous covariates (i.e. $\bm{\mu}_{(\operatorname*{\mathbf{X}})\,c}$ and $\Sigma_{(\operatorname*{\mathbf{X}})\,c}$) in order to distinguish them from the parameters for the outcome model. If all the covariates in $\operatorname*{\mathbf{X}}$ are continuous (i.e. $p_{2}$=0), the right-hand side of Equation (2) will be replaced by $f(\operatorname*{\mathbf{X}}|\bm{\mu}_{(\operatorname*{\mathbf{X}})\,c},\Sigma_{(\operatorname*{\mathbf{X}})\,c})$. On the other hand, if all the covariates in $\operatorname*{\mathbf{X}}$ are discrete (i.e. $p_{1}$=0), the right-hand side of Equation (2) will be replaced by $f(\operatorname*{\mathbf{X}}|\bm{\Psi}_{c})$.

For $\mathbf{X}^{\text{cont}}$, we assume the following probability density function:

[TABLE]

The conjugate normal-inverse-Wishart (NIW) prior is used for inference, i.e.

[TABLE]

where $\mathcal{N}$ denotes the multivariate normal distribution, $\mathcal{W}^{-1}$ denotes the inverse-Wishart distribution, ${NIW}$ denotes the normal-inverse-Wishart distribution, and $\bm{\mu}_{(\mathbf{X})}$, $\kappa_{(\mathbf{X})}$, $\nu_{(\mathbf{X})}$, and $\Lambda_{(\mathbf{X})}$ (matrix) are hyperparameters.

For discrete covariates, we assume that they are locally independent (independent conditional on the cluster assignments). Then the probability mass function of $\mathbf{X}^{\text{disc}}$ is given by

[TABLE]

where $\psi_{c,j,k}$ denotes the probability that covariate $j$ takes the value $k$ in cluster $c$, $j=p_{1}+1,\ldots,p_{1}+p_{2}$, $k=1,\ldots,K_{j}$ ($K_{j}$ denotes the number of categories for covariate $j$). Following Liverani et al. (2015), we adopt the conjugate Dirichlet prior, i.e.

[TABLE]

where $\bm{a}_{j}=(a_{j,1},a_{j,2},\ldots,a_{j,K_{j}})$, $j=p_{1}+1,\ldots,p_{1}+p_{2}$.

The model for the outcome $\operatorname*{\mathbf{Y^{*}}}$

We let $\Theta_{c}=(\bm{\mu}_{(\mathbf{Y^{*}})\,c},\Sigma_{(\mathbf{Y^{*}})\,c})$ denote the cluster-specific parameters for the multivariate normal outcome $\operatorname*{\mathbf{Y^{*}}}$ and we consider a multivariate normal (MVN) model for $\operatorname*{\mathbf{Y^{*}}}$, i.e.

[TABLE]

Similar to $\mathbf{X}^{\text{cont}}$, we specify the conjugate normal-inverse-Wishart prior for inference. Specifically,

[TABLE]

where $\bm{\mu}_{(\operatorname*{\mathbf{Y^{*}}})}$, $\kappa_{(\operatorname*{\mathbf{Y^{*}}})}$, $\nu_{(\operatorname*{\mathbf{Y^{*}}})}$ and $\Lambda_{(\operatorname*{\mathbf{Y^{*}}})}$ are hyperparameters. We note that in the profile regression model, we do not impose any specific covariance structure on $\operatorname*{\mathbf{Y^{*}}}$, and $\Sigma_{(\mathbf{Y^{*}})\,c}$ is estimated in a data-driven way.

2.2.2 Post-processing of the clustering output

Identify the “representative” clustering

We note that the proposed Bayesian mixture modeling framework (stochastic) takes into account the uncertainty associated with the number of clusters and cluster assignments, and the clustering output can vary across iterations of the Markov chain Monte Carlo (MCMC) sampler. To facilitate the interpretation of the clustering results, we “summarize” the clustering output across iterations and identify a “representative” clustering structure. In each MCMC iteration, we can construct an $n\times n$ score matrix, where the element $(i_{1},i_{2})$ equals to 1 if individuals $i_{1}$ and $i_{2}$ are allocated to the same cluster in this iteration, and equals to 0 if $i_{1}$ and $i_{2}$ are assigned to different clusters. We can then average the score matrices over all MCMC iterations to get a posterior similarity matrix, $\mathbf{S}$, which records the probability that two individuals are assigned to the same cluster (i.e. $\mathbf{S}$ records the posterior co-clustering probabilities of all pairs of individuals). The “representative” clustering can be identified as the partition of the data that best represents $\mathbf{S}$. Specifically, we use the partitioning around medoids (PAM) algorithm (Kaufman and Rousseeuw, 1990): PAM is directly applied to the posterior dissimilarity matrix $1-\mathbf{S}$, which allocates individuals to clusters in a way consistent with $\mathbf{S}$. For each fixed number of clusters up to a prespecified maximum, we select the best PAM partition that minimizes the sum of dissimilarities between the center of each cluster and all other members of the same cluster. The final “representative” clustering is then chosen by maximizing the average silhouette width (Rousseeuw, 1987) across these best PAM partitions (Molitor et al., 2010).

Quantify the uncertainty associated with the “representative” clustering

We can evaluate how confident we are about the “representative” clustering by examining whether or not across different iterations the second-stage clustering model consistently clusters individuals in a way similar to the “representative” clustering, and we would expect the credible intervals associated with cluster parameter estimates to be narrower for more consistent clustering (stronger clustering signal). To this end, the model-averaging approach discussed in Molitor et al. (2010) can be employed under our proposed framework.

2.2.3 Variable selection in profile regression

Variable selection methods can be embedded into profile regression in order to identify the covariates that contribute significantly to the formation of clusters. In our implementation, we follow the variable selection approach taken by Liverani et al. (2015).

Continuous covariates

We let $\bm{\gamma}^{\text{cont}}_{c}=(\gamma_{c,1},\gamma_{c,2},\ldots,\gamma_{c,\,p_{1}})$, where $\gamma_{c,j}$ is a binary random variable that determines whether or not covariate $j$, $j=1,\ldots,p_{1}$, is important for allocating subjects to cluster $c$ ($\gamma_{c,j}=1$ if the answer is yes and 0 otherwise). Let $\bar{x}_{j}$ denote the average value of covariate $j$ (sample average), for $j=1,\ldots,p_{1}$, and we define $\bm{\mu}^{*}_{c}=(\mu^{*}_{c,1},\mu^{*}_{c,2},\ldots,\mu^{*}_{c,\,p_{1}})$, where

[TABLE]

where $\mu_{c,j}$ is the $j^{\text{\,th}}$ element of $\bm{\mu}_{(\operatorname*{\mathbf{X}})\,c}$, $j=1,\ldots,p_{1}$. We then replace $\bm{\mu}_{(\operatorname*{\mathbf{X}})\,c}$ in (3) with $\bm{\mu}^{*}_{c}$. We assume that $\gamma_{c,j}\sim\text{Bernoulli}(\rho_{j})$, $j=1,\ldots,p_{1}$. A sparsity inducing prior is used for $\rho_{j}$ (Papathomas et al., 2012; Liverani et al., 2015).

Discrete covariates

We can perform variable selection on discrete covariates in a similar manner. We let $\bm{\gamma}^{\text{disc}}_{c}=(\gamma_{c,\,p_{1}+1},\gamma_{c,\,p_{1}+2},\ldots,\gamma_{c,\,p_{1}+p_{2}})$, where $\gamma_{c,j}$ denotes a binary variable indicating whether or not covariate $j$, $j=p_{1}+1,\ldots,p_{1}+p_{2}$, is important for assigning subjects to cluster $c$. Let $\psi_{c,j,k}$ denote the probability that covariate $j$ takes the value $k$ in cluster $c$, and $\psi_{0,j,k}$ be the observed proportion of covariate $j$ taking the value $k$, $j=p_{1}+1,\ldots,p_{1}+p_{2}$, $k=1,\ldots,K_{j}$. To incorporate the variable selection feature, the discrete covariate model (4) is modified accordingly to

[TABLE]

Similar to the continuous covariate case, we assume that $\gamma_{c,j}\sim\text{Bernoulli}(\rho_{j})$, and each $\rho_{j}$ is assigned a sparsity inducing prior, $j=p_{1}+1,\ldots,p_{1}+p_{2}$.

3 Simulation studies

In this section, we evaluate the performance of our proposed method via simulation studies.

3.1 Simulation design

We consider two simulation scenarios. The first scenario corresponds to the case where all covariates are continuous, and we evaluate the performance of the proposed method under different levels of correlations between covariates. In the second scenario, covariates include a mix of continuous and discrete ones and cluster sizes are unequal. We also extend this scenario to evaluate the performance of the variable selection procedure embedded in our method. We examine two sample sizes: $n=450$ and $n=900$, and we repeat the simulation 100 times for each scenario. Details on the simulation design are given in the following.

3.1.1 Scenario 1

In this setting, we consider a three-arm trial ($K=3$) where treatment $A$ is sampled from {1, 2, 3} with equal probabilities. $\operatorname*{\mathbf{X}}$ consists of 3 normally-distributed continuous covariates: $X_{1}$, $X_{2}$, and $X_{3}$. There are 3 underlying clusters of equal size under this scenario. The mean values of covariates and expected potential outcomes under each treatment option for 3 clusters are summarized in Table 1 (a). We also plot the expected treatment response profiles of subjects in each cluster in Figure 1 (a) for better visualization. The observed outcome $Y$ is simulated from a normal distribution with mean $\sum_{a=1}^{3}E\{Y^{*}(a)\}I(A=a)$ and variance $\sigma^{2}_{Y}$. We consider different noise levels in covariates ($\sigma_{\operatorname*{\mathbf{X}}}=0.5,1$) and noise levels in the outcome ($\sigma_{Y}=0.5,1$), which reflect different degrees of cluster separability. In addition, we vary the correlation between $X_{1}$ and $X_{2}$ ($\rho_{\operatorname*{\mathbf{X}}}=0$, 0.5, 0.8 conditional on the cluster assignment; other pairwise correlations are set as 0 conditional on the cluster assignment) to examine the influence of the degree of between-covariate correlations on our proposed method’s performance.

3.1.2 Scenario 2

In the second setting, we consider a four-arm trial ($K=4$) where treatment $A$ is sampled uniformly from {1, 2, 3, 4}. There exist 4 clusters, and the sizes of clusters 1-4 are ${n}/{9}$, ${2n}/{9}$, ${2n}/{9}$, and ${4n}/{9}$, respectively.

We consider the situation where there are four covariates, all of which inform the clustering structure. In particular, $X_{1}$ and $X_{2}$ are binary variables that take on the value of 0 or 1, $X_{3}$ is a categorical variable taking on values 0, 1, or 2, and $X_{4}$ is a normally-distributed continuous variable. In Table 1 (b), we summarize the covariate profiles and expected treatment response profiles for 4 clusters. The expected treatment response profiles (by cluster) are also plotted in Figure 1 (b). The observed outcome $Y$ follows a normal distribution with mean $\sum_{a=1}^{4}E\{Y^{*}(a)\}I(A=a)$ and variance $\sigma^{2}_{Y}$. Four different combinations of noise levels in $X_{4}$ and noise levels in the outcome $Y$ are considered, namely, $(\sigma_{Y},\sigma_{\operatorname*{\mathbf{X}}})=(0.2,0.2)$, $(0.5,0.5)$, $(1,0.5)$, $(1,1)$. We also consider a modification of scenario 2 to evaluate the performance of the proposed variable selection method. Details are provided in Web Appendix C.

3.2 Evaluation metrics

The performance of our proposed method in simulation studies is evaluated based on the following metrics:

Estimated number of clusters, $\widehat{N}_{\text{cluster}}\,$. 2. 2.

Adjusted Rand index (ARI), which measures the similarity between the computed clustering structure and the ground truth (i.e. ideal clustering, we will refer to it as “class structure” to distinguish it from the computed clustering structure) (Hubert and Arabie, 1985). ARI equals to 1 if two partitions of the data are identical and a larger value of ARI indicates a higher level of agreement between two partitions. 3. 3.

Conditional entropy-based metrics (Rosenberg and Hirschberg, 2007):

1. (a)

Homogeneity, which measures whether all members of a given computed cluster are from the same class (ground truth). Homogeneity score is larger (more desirable) when each cluster contains elements of fewer classes. 2. (b)

Completeness, which measures whether all members of a given class (ground truth) are in the same computed cluster. Completeness score is larger (more desirable) when the members of a given class are allocated to fewer clusters.

Both homogeneity and completeness are bounded between 0 and 1.

Details on how to compute ARI, homogeneity, and completeness are provided in Web Appendix B.

3.3 Simulation results

In each simulation replicate, BART is run for 6000 iterations to predict the potential outcomes and profile regression is run for 2000 iterations to cluster subjects into subgroups. For both BART and profile regression, the first 1000 iterations are discarded as burn-in. We use the default values specified in the BART package and the PReMiuMar package for all hyperparameters. Investigation of posterior samples suggests no evidence against convergence in our simulation studies.

3.3.1 Scenario 1

Simulation results for scenario 1 are presented in Table 2.

The results are fairly robust across different levels of correlations ($\rho_{\operatorname*{\mathbf{X}}}$) among covariates. This is a highly desirable property of our proposed approach compared to many other standard approaches, whose performance can be very sensitive to the degree of correlations among covariates due to the multicollinearity problem (Molitor et al., 2010).

The $\widehat{N}_{\text{cluster}}$ column shows that our proposed method over-estimates the true number of clusters (${N}_{\text{cluster}}=3$ in this scenario), especially when the sample size is large ($n=900$). This is not surprising given that the profile regression adopts a Dirichlet process mixture model (DPMM): it has been demonstrated that when the true number of clusters is finite and small, the posterior inference on the number of clusters by using the DPMM may be inconsistent, and DPMM tends to over-estimate the true number of clusters and produce some small extraneous clusters around the true components (Onogi et al., 2011; Miller and Harrison, 2013; Yang et al., 2019; Lu et al., 2018). This phenomenon is referred to as “over-clustering” by Lu et al. (2018), which might be due to the sensitivity of DPMM to even minor deviations that exist among clusters. As has been discussed in Section 2.2.1, the choice of the concentration parameter $\alpha$ implicitly affects the expected number of (non-empty) clusters. Results presented in Table 2 are obtained by assuming a Gamma(2,1) prior for $\alpha$ (i.e. the default in the PReMiuMar package; $E(\alpha)=2$). We also examine the empirical results obtained under other commonly-used priors for $\alpha$, for example, Gamma(2,4) ($E(\alpha)=0.5$) suggested by Escobar and West (1995), and the prior that is matched to the sparse finite mixture model (Frühwirth-Schnatter and Malsiner-Walli, 2019). The estimated number of clusters, $\widehat{N}_{\text{cluster}}$, are very similar across different prior specifications for $\alpha$ (not presented in the paper), even though based on theoretical considerations, the results might differ (Frühwirth-Schnatter and Malsiner-Walli, 2019).

Not surprisingly, homogeneity is higher (better) than completeness in most cases, given that $N_{\text{cluster}}$ is typically over-estimated by our method. In this sense, it is more likely that subjects who are clustered together are from the same class (the ground truth in simulation studies), and thus leading to “more homogeneous” clusters. On the other hand, small extraneous clusters (centered around true clusters) that are produced by DPMM may explain the “less complete” clustering results.

As expected, higher noise levels in covariates ($\sigma_{\operatorname*{\mathbf{X}}}$) and outcomes ($\sigma_{Y}$) result in worse clustering performance (i.e. lower ARI, completeness, and homogeneity, and a larger upward bias in the estimation of ${N}_{\text{cluster}}$), given that we would expect the underlying clustering structure to be less clear (lower cluster separability) with larger values of $\sigma_{\operatorname*{\mathbf{X}}}$ and $\sigma_{Y}$.

A comparison of the results obtained when $n=450$ with those obtained when $n=900$ suggests that larger sample sizes do not seem to improve the performance of our method in this setting. When the noises in covariates and the outcome are large ($\sigma_{Y}=\sigma_{\operatorname*{\mathbf{X}}}=1$), ARI and completeness even get worse as the sample size increases. One possible reason for this observation is that the over-estimation of $N_{\text{cluster}}$ by DPMM is more pronounced when $n=900$, and this has a subsequent (negative) effect on ARI and completeness.

In addition to measuring the overall clustering accuracy of our proposed method based on the clustering performance metrics described in Section 3.2, we also examine how well our method performs in terms of recovering the true underlying cluster-specific parameters (e.g. means). To this end, in each simulation replicate, we first obtain the posterior summaries of mean parameters for $X_{1}$, $X_{2}$, $X_{3}$, $Y^{*}(1)$, $Y^{*}(2)$, and $Y^{*}(3)$ in each resulting cluster, and these cluster-specific parameters are then re-weighted by the cluster sizes. We note that the re-weighting step is important: posterior summary statistics for larger clusters should be assigned more weight since larger clusters carry more information. The densities of the re-weighted results over 100 simulation replicates in low ($\sigma_{Y}=\sigma_{\operatorname*{\mathbf{X}}}=0.5$) and high ($\sigma_{Y}=\sigma_{\operatorname*{\mathbf{X}}}=1$) noise settings are plotted in Figure 2 ((a) for covariates $X_{1}$, $X_{2}$, and $X_{3}$ and (b) for potential outcomes under each treatment option). These plots correspond to the case with $n=900$ and $\rho_{\operatorname*{\mathbf{X}}}=0.5$. Density plots corresponding to other $n$ and $\rho_{\operatorname*{\mathbf{X}}}$ values look similar (do not alter our conclusions) and are thus omitted.

Figure 2 implies that the true cluster-specific means (dashed lines) for all covariates and potential outcomes can be recovered by our proposed method, even though the estimated number of clusters is greater than the truth. When the noise level is high ($\sigma_{Y}=\sigma_{\operatorname*{\mathbf{X}}}=1$), the densities are flatter (nosier) compared to the case with low noise levels ($\sigma_{Y}=\sigma_{\operatorname*{\mathbf{X}}}=0.5$), but our method still manages to recover the truth.

3.3.2 Scenario 2

Simulation results under scenario 2 are summarized in Table 3.

As in the first scenario, the homogeneity score is almost always higher than the completeness score due to “over-clustering”. When the noise level is low ($\sigma_{Y}=\sigma_{\operatorname*{\mathbf{X}}}=0.2$), the clusters are well-separated and our proposed method performs almost perfectly in terms of ARI, completeness and homogeneity, despite the over-estimation of $N_{\text{cluster}}$ remaining a problem (for this scenario, $N_{\text{cluster}}=4$). As the noise level increases, the clustering performance gets worse, and it seems to be more sensitive to the change in $\sigma_{\operatorname*{\mathbf{X}}}$ than to the change in $\sigma_{Y}$ (when we compare the results corresponding to $(\sigma_{Y},\sigma_{\operatorname*{\mathbf{X}}})=(1,0.5)$ with those corresponding to $(\sigma_{Y},\sigma_{\operatorname*{\mathbf{X}}})=(1,1)$). We also observe that $\widehat{N}_{\text{cluster}}$ first increases and then decreases as the noise gets larger, possibly because our proposed algorithm merges some clusters when the true underlying clustering structure has very low separability (i.e. very large noise).

4 Application to the INTERVAL trial

We apply our proposed patient stratification approach to the data from the INTERVAL trial. The purpose of this analysis is to uncover donor subgroups with different baseline characteristics and potential “treatment” (inter-donation interval in the blood donation context) response profiles. Specifically, we focus our analysis on a “much-in-demand but vulnerable” donor population of 884 female donors who were younger than 40 (more at risk of iron deficiency after donating blood) and had O negative blood type (the “universal” blood group that can be transfused to any patient in need and used in medical emergencies). The three randomized groups for female donors are 16-week, 14-week, and 12-week inter-donation intervals. The target outcome of our interest is a utility score (denoted by $U$) that accounts for the trade-off between the total units of blood collected per donor over the 2-year trial period (the benefit outcome, denoted by $G$) and the number of low hemoglobin (Hb) deferrals per donor during the same period (the risk outcome, denoted by $R$), i.e. $U=G-b\times R$, where $b$ is the trade-off parameter reflecting the equivalent benefit loss for one unit increase in the risk. We examine the case with $b=3$ to reflect the potential costs of reduced efficiency of blood collection and reduced donor retention due to low Hb deferrals. Seven baseline donor characteristics are considered, including the Short Form Health Survey version 2 (SF-36v2) physical component score (PCS), mental component score (MCS), ferritin level, red blood cell count (RBC), mean corpuscular volume (MCV), mean corpuscular hemoglobin (MCH), and body mass index (BMI).

In our analysis, we use the default setups that are specified in the BART package and the PReMiuMar package for priors and hyperparameters. We run BART for 6000 MCMC iterations with an initial burn-in of 1000 iterations to predict potential outcomes, and run profile regression for 40000 iterations with a burn-in of 10000 iterations in the clustering step. We do not find strong evidence against convergence. More details on convergence diagnostics are given in Web Appendix D.

Unlike in the simulation studies where we can evaluate the performance of our method by calculating external metrics such as ARI, homogeneity, and completeness, in the real data application, the ground truth is not known and thus these external validation metrics cannot be used. Instead, we focus on the interpretation of the “representative” clustering of 884 female donors in order to assess whether or not our proposed approach can give insights into donor heterogeneity and stratify donors into clinically meaningful subgroups.

Inspection of the raw output from the profile regression model (before applying the post-processing method described in Section 2.2.2) suggests that at each MCMC iteration, there are either 6 or 7 resulting non-empty clusters, with 5 moderately-sized clusters and 1 or 2 very small clusters. In particular, for all MCMC iterations, one of the resulting clusters contains only 1 donor whose “red blood cell-related” blood measurements are fairly extreme (lowest RBC, highest MCH and $3^{\text{rd}}$ highest MCV among all 884 female donors). The other small cluster (if it exists, i.e. when the total number of non-empty clusters is 7) includes donors with extremely low PCS (the size of this cluster varies slightly across different iterations, but it is always less than 8).

We apply the post-processing method discussed in Section 2.2.2 to identify the “representative” clustering based on posterior samples. The heatmap of the posterior similarity matrix for 884 female donors is presented in Web Appendix E (Web Figure 4). The “representative” clustering consists of 5 donor subgroups, and the sizes of these subgroups are 171, 126, 101, 93, and 393, respectively. This indicates that the 1 or 2 very small cluster(s) in the raw output are merged into larger clusters in this post-processing step.

Figure 3 shows the covariate profiles (posterior distributions of the mean parameters for donors’ baseline characteristics) and the potential outcome profiles (posterior distributions of the mean parameters for potential outcomes under 16-, 14-, and 12-week inter-donation intervals) corresponding to each of the five donor subgroups. These results reveal considerable evidence for the presence of heterogeneity within the donor population under investigation.

Cluster 4 (with 126 donors) represents the potentially “super donor” subgroup who can give blood more frequently than the current clinical practice and for these donors, more frequent donations lead to larger utilities (Figure 3 (b)). Comparing across subgroups, this group of donors has higher donation capacity since the utility scores are on average higher, especially under the 14-week and 12-week inter-donation intervals. These “super donors” are characterized by having high levels of ferritin, RBC, BMI, MCS, MCV and MCH and low levels of PCS. In particular, the most distinguishing feature that differentiates this subgroup from the remaining subgroups is the ferritin level (significantly higher in this subgroup). This is consistent with a priori expectation that donors with higher ferritin levels are in general more capable of donating blood more often (especially for female donors).

Cluster 2 is the largest subgroup with 393 donors. Figure 3 (b) suggests that donors in this subgroup have the potential to donate blood every 12 weeks, even though the gain in the utility score by switching from the 16-week inter-donation interval to the 12-week inter-donation interval is only moderate and not as significant compared to cluster 4. Consistent with cluster 4, donors in cluster 2 on average have high MCS, MCV, and MCH. However, their ferritin levels, RBC and BMI are lower than those of donors in cluster 4 (Figure 3 (a)).

Cluster 1 (with 101 donors) represents a potentially “brittle donor” subgroup. Utility score gets smaller as the inter-donation interval gets shorter (Figure 3 (b)). The low ferritin levels, RBC and MCS of these donors (Figure 3 (a)) may potentially explain why they are vulnerable. As an aside, the cluster with only one donor in the raw output before post-processing is merged into this “brittle donor” subgroup. This is sensible since blood-based measurements of this donor suggest that she may be “vulnerable”.

Cluster 5 (with 171 donors) also represents a subgroup of donors with low donation capacity. The utility scores of donors in this subgroup are similar under three inter-donation interval options and they are lower than the utility scores of donors in clusters 2, 3, and 4 (Figure 3 (b)). Donors in cluster 5 (on average) have a high MCS similar to that of donors in cluster 4. However, low values of donation capacity-related characteristics such as ferritin levels, MCV, MCH, and BMI (Figure 3 (a)) may lead to higher than average deferral rates and lower than average utility scores in this subgroup.

The utility scores of donors in cluster 3 (with 93 donors) are almost identical under 16-week and 14-week inter-donation intervals but decrease significantly when the donation frequency is every 12 weeks (Figure 3 (b)). The ferritin levels, MCV, MCH, and BMI are relatively high (Figure 3 (a)) for these donors, and the drop in the utility score when the inter-donation interval is 12-week may be attributed to the low MCS.

To summarize, our analysis of the INTERVAL data using the proposed approach reveals clinically meaningful donor subgroups within the “much-in-demand but vulnerable” donor population under investigation.

5 Discussion

The uncovering of subgroups from a heterogeneous population plays an important role in precision medicine applications. In this paper, we present a two-stage patient stratification approach that leverages Bayesian nonparametric techniques. Our proposed method captures the heterogeneity in the underlying population and clusters the population into subgroups of subjects who share similar covariate profiles and display similar treatment responses. Specifically, in the first stage, we predict the potential outcome under each treatment arm for each patient, and in the second stage, we apply profile regression to link the multivariate potential outcome vector to a set of covariates (can be continuous, discrete, or a mix of continuous and discrete ones) through cluster membership (Molitor et al., 2010). Based on the posterior samples, the resulting subgroups can be characterized in terms of covariate and treatment response profiles.

Our method offers several advantages. Firstly, while most existing methods for subgroup identification only cover two-treatment cases, our proposed approach is applicable to multi-arm trials. Secondly, the use of the Dirichlet process prior allows the number of clusters to be estimated from the data, thus bypassing the need for pre-specifying it. Thirdly, our method can properly handle correlated covariates (avoiding well-known problems caused by multicollinearity), which are common in clinical studies. Fourthly, a variable selection procedure is embedded into our model for identifying important covariates that actively “drive” the clustering structure (i.e. contribute significantly to the cluster patterns). Lastly, the proposed approach is built under the Bayesian framework and takes into account model uncertainties (Molitor et al., 2010).

The application of our proposed method to a subset of the INTERVAL data (a “much-in-demand but vulnerable” donor population) identifies 5 clinically meaningful donor subgroups with different donation capacities and covariate (donors’ baseline characteristics) profiles. These results provide insight into the underlying donor heterogeneity by highlighting the differences between donors in terms of both baseline characteristics and potential response (to three different inter-donation intervals) patterns, and can be leveraged to inform and guide targeted donor recruitment and donor management strategy. For example, donors who are identified as belonging to the “super donor” subgroup may be asked to give blood more frequently if there is a blood shortage or if their blood group is rare or universal. In contrast, donors who belong to the “brittle donor” subgroup will be allowed longer time between donations to ensure donor health and safety (BTRU, 2019).

The use of the Dirichlet process prior in the second-stage profile regression model allows the number of clusters to be discovered in a data-driven way. However, as demonstrated by our simulation studies, Dirichlet process mixture model (DPMM) tends to over-estimate the number of clusters (produce some superfluous small-sized clusters). Indeed, the inconsistent inference (over-estimation) on the number of components by DPMM is a well-known problem (referred to as “over-clustering” in Lu et al. (2018)) that has been empirically observed and reported in the literature before, and it appears that in most cases, the extra clusters are centered around the true components and only include a very small number of individuals (Onogi et al., 2011; Miller and Harrison, 2013). To our knowledge, how to correct for such inconsistency remains an open question in the field of Bayesian mixture modeling (Miller and Harrison, 2013; Yang et al., 2019). In general, we think that a slight over-estimation of the true number of clusters is not an issue of major concern in the context of patient stratification, since in this case, the primary interest typically lies in the characterization and the interpretation of clusters rather than the inference on the exact number of underlying clusters. As has been noted by Onogi et al. (2011), even though the extra clusters produced by DPMM are considered as redundant and interpreted as over-estimation in simulation studies, they may provide useful information in real data applications since they reflect some relatively subtle heterogeneity that might be clinically interesting. In practice, in order to achieve better interpretability of the clustering results, the model-based “representative” clustering should be coupled with practitioners’ subject-matter knowledge when determining the optimal number of clusters. If the inspection of the output indicates that the small clusters produced by DPMM are not of much clinical interest because they do not reflect a general pattern, the other larger and more representative clusters will be given more emphasis in the interpretation. In addition, if the profiles of some small clusters are very similar to those of some much larger clusters in terms of clinical meaningfulness, we can merge them a posteriori. From our point of view, this would be preferable to a method of low granularity, in which case some useful information on the population heterogeneity may be overlooked.

The work in this paper has raised new research questions that are worthy of further investigation. Although our method is primarily developed for the setting with a univariate and continuous target outcome, it can be extended to the multi-outcome setting.

Since we use the posterior mean of the potential outcome predictions (from BART) rather than the full posterior samples as the outcome in the profile regression model, the uncertainties associated with the predicted potential outcomes (the first stage) are not carried forward to the second stage and thus not reflected in the final output. We may employ the idea of Markov melding (a generic Bayesian computational method for evidence synthesis) to allow uncertainty propagation (Goudie et al., 2019).

Numerical experiments suggest that our proposed patient stratification method does not scale well to large datasets, mainly due to the computational inefficiency of the current implementation of profile regression with a multivariate normal outcome in the PReMiuMar package. MCMC sampling can be prohibitively slow, and a variational inference algorithm for Dirichlet process mixture models may be considered to speed up the computation (Blei and Jordan, 2006).

In the profile regression model, both covariates ($\operatorname*{\mathbf{X}}$) and outcome ($\operatorname*{\mathbf{Y^{*}}}$) inform the clustering structure. If the dimension of $\operatorname*{\mathbf{X}}$ is much higher than that of $\operatorname*{\mathbf{Y^{*}}}$, the contribution of $\operatorname*{\mathbf{Y^{*}}}$ to the likelihood is likely to be overwhelmed by that of $\operatorname*{\mathbf{X}}$ (i.e. covariates might dominate the likelihood and the relative contribution of the outcome may be undermined). Consequently, the impact of the response data on the cluster allocation will be small and the resulting clusters will be formed mainly based on the similarity in the covariate space. Bigelow and Dunson (2009) argued that this might be a desirable property for some epidemiological studies. However, in some other cases where we expect the outcome to play a more important role in the clustering, we might need to upweight the outcome likelihood. In practice, the weight may be subjective and may depend heavily on the contexts and research aims.

Acknowledgments

This work was supported by the UK Medical Research Council programme MC_UU_00002/2 and the Cambridge International Scholarship. Participants in the INTERVAL trial were recruited with the active collaboration of NHS Blood and Transplant England (www.nhsbt.nhs.uk), which has supported field work and other elements of the trial. The academic coordinating centre for INTERVAL was supported by core funding from: NIHR Blood and Transplant Research Unit in Donor Health and Genomics (NIHR BTRU-2014-10024), UK Medical Research Council (MR/L003120/1), British Heart Foundation (SP/09/002; RG/13/13/30194; RG/18/13/33946) and the NIHR [Cambridge Biomedical Research Centre at the Cambridge University Hospitals NHS Foundation Trust]. A complete list of the investigators and contributors to the INTERVAL trial is provided in Di Angelantonio et al. (2017). The academic coordinating centre would like to thank blood donor centre staffs and blood donors for participating in the INTERVAL trial. This work was also supported by Health Data Research UK, which is funded by the UK Medical Research Council, Engineering and Physical Sciences Research Council, Economic and Social Research Council, Department of Health and Social Care (England), Chief Scientist Office of the Scottish Government Health and Social Care Directorates, Health and Social Care Research and Development Division (Welsh Government), Public Health Agency (Northern Ireland), British Heart Foundation and Wellcome. The views expressed in this paper are those of the authors and not necessarily those of the NHS, the NIHR or the Department of Health and Social Care.

Supporting Information

Additional supporting information may be found online in the Supporting Information section at the end of the article. The R code for implementing the proposed patient stratification method is available at https://github.com/yx299/stratification.