# Sharp Inference on Selected Subgroups in Observational Studies

Xinzhou Guo1 Linqing Wei2 Chong Wu3 Jingshen Wang2 Correspondence:
<EMAIL_ADDRESS>The research of Jingshen Wang was supported in part
by the National Science Foundation (DMS 2015325).

(1Harvard-MIT Center for Regulatory Science
2Division of Biostatistics, UC Berkeley
3Department of Statistics, Florida State University
)

###### Abstract

In modern drug development, the broader availability of high-dimensional
observational data provides opportunities for scientist to explore subgroup
heterogeneity, especially when randomized clinical trials are unavailable due
to cost and ethical constraints. However, a common practice that naively
searches the subgroup with a high treatment level is often misleading due to
the “subgroup selection bias.” More importantly, the nature of high-
dimensional observational data has further exacerbated the challenge of
accurately estimating the subgroup treatment effects. To resolve these issues,
we provide new inferential tools based on resampling to assess the
replicability of post-hoc identified subgroups from observational studies.
Through careful theoretical justification and extensive simulations, we show
that our proposed approach delivers asymptotically sharp confidence intervals
and debiased estimates for the selected subgroup treatment effects in the
presence of high-dimensional covariates. We further demonstrate the merit of
the proposed methods by analyzing the UK Biobank data. The R package
“debiased.subgroup” implementing the proposed procedures is available on
GitHub.

Keywords: Bootstrap; Precision Medicine; Debiased Inference.

##  1 Introduction

###  1.1 Motivation and our contribution

Subgroup analysis, broadly speaking, aims to uncover and confirm treatment
effect heterogeneity within a population, and it has been frequently applied
to randomized clinical trials (RCTs) for drug development (Hébert et al.,
2002; Alosh et al., 2017). In recent years, the combination of increasing
availability of observational data and advancements in statistical methods and
computing capacity has stimulated researchers’ interest in using observational
data for subgroup analysis. This has spawned the hope for deeper findings and
more targeted recommendations or interventions in a wide variety of areas from
education, health care, to marketing. Although RCTs remain the gold standard
for assessing the efficacy of clinical interventions, observational data offer
broader opportunities to explore subgroup heterogeneity when we cannot afford
RCTs due to cost and ethical constraints.

While observational studies can be useful to identify subgroups with favorable
treatment efficacy or unfavorable adverse effects, we must account for the
impact of subgroup selection on any subsequent evaluation of the subgroup. We
have been frequently reminded by the failure of follow-up trials to confirm a
seemingly promising subgroup (see Kubota et al., 2014, for example), and by
the discussions in the domain science journals about the appeals and pitfalls
of subgroup analysis (Petticrew et al., 2012, e.g.). When subgroups are
identified post hoc from RCTs, it has been recognized in the literature (Cook
et al., 2014; Bornkamp et al., 2017; Guo and He, 2019) that the subgroup
selection bias, if not accounted for, is likely to lead to overtreatments and
false discoveries. When subgroups are identified from observational studies,
the problem becomes even more challenging as subgroup treatment effects need
to be estimated upon adjusting for possibly high-dimensional confounders. To
address these issues, in this paper, we provide new inferential tools to help
assess the replicability of post-hoc identified subgroups from observational
studies without having to resort to simultaneous inference methods that are
often too conservative to start with.

From a statistical methodological standpoint, our proposed approach delivers
asymptotically sharp confidence intervals as well as debiased estimates for
the selected subgroup treatment effects in the presence of high-dimensional
covariates. We break down our methodological contributions as follows:

First, we allow well-understood debiased regression parameter estimates to be
used as the building blocks for the subgroup treatment effects from
observational studies, thus enabling adjustments for high-dimensional
covariates in the study. In particular, we investigate how debiased Lasso and
repeated data splitting can be used in the subgroup analysis. We find that the
former works well at handling a large number of subgroup-specific parameters,
whereas the latter achieves better statistical efficiency at sparse models
especially when the subgroup assignments are highly correlated with other
covariates in the model.

Second, we propose new bootstrap-calibrated procedures for quantifying the
selected subgroup treatment effect (Section 3.1 and 3.2) and provide a
theoretical guarantee of its asymptotic validity in high dimensions (Section
4). Most importantly, the resulting confidence intervals are asymptotically
sharp in the sense that they approximate the nominal coverage probability
without overshooting it (Theorem 1 and Corollary 1). While the existing
statistical literature (Zhang and Cheng, 2017; Dezeure et al., 2017) has
argued for the use of multiplicity adjustment to ensure post-selection
inference validity in high dimensions, most recommended methods tend to
sacrifice power as they aim to protect the family-wise error rates over all
possible subgroups. The proposed method accounts automatically for subgroup
selection bias as well as increasingly complex correlation structures among
the subgroups, leading to targeted (instead of simultaneous) statistical
validity in subgroup analysis.

###  1.2 Related literature

Our paper is closely related to subgroup analysis. Here, the literature is
often divided into exploratory subgroup analysis that focuses on subgroup
identification (Shen and He, 2015; Su et al., 2009), and confirmatory subgroup
analysis that aims to validate subgroups identified from an earlier stage
(Jenkins et al., 2011; Friede et al., 2012). These classical methods are
usually made for randomized trials and are not applicable for observational
studies. While some recent literature on subgroup analysis (Fan et al., 2017;
Yang et al., 2020; Izem et al., 2020) accounts for the pre-treatment
information, their generalization to high-dimensional observation studies may
require substantial modification. Beyond subgroup analysis, some methods on
analyzing heterogeneous treatment effects (Imai et al., 2013; Wager and Athey,
2017) might be also applicable for subgroup analysis in observational studies.
Different from our goal of simultaneously identifying and estimating subgroup
treatment effect, their goal mainly focuses on subgroup identification or
making inference on a pre-defined subgroup.

Our methodology also contributes to high-dimensional inference and post-
selection inference literature, simply because accurate point estimation of
treatment effects is usually not available without regularization (e.g.
Tibshirani, 1996; Fan and Li, 2001) in high dimensions. Selective inference
(Lee et al., 2016; Tian et al., 2018) constructs exact confidence intervals
for the selected regression coefficients based on Lasso conditional on the
selected model, indicating that their framework tends to provide conservative
confidence intervals overall. While recent developments in debiased inference
(Zhang and Zhang, 2014; Van de Geer et al., 2014; Belloni et al., 2014)
removes the Lasso regularization bias by using an estimate of the inverse
population covariance matrix, they do not address the issue of subgroup
selection bias.

Our framework is connected to the multiple comparison literature. Although
several attempts have been made to address the multiple comparison issues in
subgroup analysis (Hall and Miller, 2010; Fuentes et al., 2018; Stallard et
al., 2008), those procedures are either conservative or poorly grounded. Guo
and He (2020) propose a asymptotically sharp procedure to adjust the subgroup
selection bias but it is not designed for high-dimensional observational
studies. In addition, most false discovery controls that are developed for
independent tests cannot be well justified for multiple analyses based on the
same data set. When they protect false positive rates, they tend to be
conservative too. In our proposal, we handle the dependence among tests by
using the bootstrap-based calibration.

##  2 Model setup and challenges

In classical subgroup analysis in RCTs, a subgroup is usually defined as a
subpopulation of patients decided by their baseline characteristics, and the
corresponding subgroup treatment effect is the effect of a certain medical
treatment for this given subpopulation. As observational studies often allow
us to investigate multiple treatment strategies at the same time, subgroups
can be more broadly defined. In this paper, other than simply viewing
subgroups in the classical RCT setting, we also consider subgroups that are
defined as subpopulations of patients that receive different treatments. We
refer the variables that capture the differential subgroup treatment effects
as treatment effect heterogeneity variables.

Suppose we have a random sample of $n$ i.i.d. observations
$\\{(Y_{i},Z_{i},X_{i})\\}_{i=1}^{n}$, where $Y_{i}$ is the response variable,
$Z_{i}$ is a $p_{1}$-dimensional vector of treatment effect heterogeneity
variables (See Remark 1 for its construction in different scenarios), and
$X_{i}$ is a $p_{2}$-dimensional vector of covariates for the $i$-th subject.
Without loss of generality, we assume that a larger value of $\beta_{i}$ means
a better treatment effect. We work under the high-dimensional setting where
$p_{1}+p_{2}\gg n$. We assume that

$Y_{i}=Z_{i}^{\prime}\beta+X_{i}^{\prime}\gamma+\varepsilon_{i},\hskip
28.45274pti=1,\dots,n,$ (1)

where $\beta=(\beta_{1},\dots,\beta_{p_{1}})^{\prime}$ captures the subgroup
treatment effects, $\gamma$ is a sparse vector of coefficients, and
$\varepsilon_{i}$’s are independent random errors with
$\mathbb{E}(\varepsilon_{i}|Z_{i},W_{i})=0$. Define
$[p_{1}]=\\{1,\dots,p_{1}\\}$. In practice, scientists often hope to identify
the subgroup with the highest treatment level. Our goal is to find a
consistent point estimate and a confidence interval with asymptotically sharp
coverage probability on:

(i) the selected subgroup treatment effect: $\beta_{\widehat{s}},\text{ where
}\widehat{s}=\arg\max_{j\in[p_{1}]}\widehat{\beta}_{j}$,
(ii) the maximal subgroup treatment effect:
$\beta_{\max}=\max_{j\in[p_{1}]}\beta_{j}$,
where $\widehat{\beta}_{i}$ is an estimate of $\beta_{i}$ (shall be specified
in Sections 3.1 and 3.2).

While one may debate which perspective is practically more relevant in
subgroup analysis, our proposed inferential procedure works for both
quantities. We will start with inference on $\beta_{\max}$ and show in Section
4 (Corollary 1) that the same procedure works for making inference on
$\beta_{\widehat{s}}$.

To fully realize the challenges of making inference on $\beta_{\max}$ in high
dimensions, we now provide an illustrative example. By showing this example,
our take-away point for practitioners is clear: Even if $\beta$ can be
estimated accurately, methods without appropriate adjustments of random
selection cannot produce valid inference on $\beta_{\max}$ (or
$\beta_{\widehat{s}}$).

######  Example 1 (Selection bias and regularization bias in estimating
$\beta_{\max}$).

In practice, $\beta_{\max}$ is frequently estimated in a two-step procedure:
One first obtains an estimate $\widehat{\beta}$ and then estimates
$\beta_{\max}$ by taking the maximum
$\max\\{\widehat{\beta}_{1},\ldots,\widehat{\beta}_{p_{1}}\\}$. Here, we use
two widely adopted procedures to estimate $\beta$ in high dimensions: (1)
Lasso, which refers to simply using the estimates
$\big{(}\widehat{\beta}_{\text{Lasso}}^{\prime},\widehat{\gamma}_{\text{Lasso}}^{\prime}\big{)}^{\prime}$
from the $\ell_{1}$-penalized regression program without any adjustments, and
(2) Refitted Lasso, which is obtained by refitting the linear model by the
ordinary least squares (OLS) based on the covariates in the support set of
$\big{(}\widehat{\beta}_{\text{Lasso}}^{\prime},\widehat{\gamma}_{\text{Lasso}}^{\prime}\big{)}^{\prime}$.
As a benchmark, we also report the performance of the oracle estimator
$\big{(}\widehat{\beta}_{\text{Oracle}}^{\prime},\widehat{\gamma}_{\text{Oracle}}^{\prime}\big{)}^{\prime}$
which pretend the true support set of $\gamma$ is known and is estimated by
OLS. We generate Monte Carlo samples following the setup in Model (1). We
generate $Z_{i1}=\mathbf{1}_{\\{W_{i1}>0.5\\}}$,
$Z_{i2}=\mathbf{1}_{\\{W_{i2}>0.5\\}}$, $X_{ij}=W_{ij}$ for $3\leq j\leq p$,
and $W_{i}\sim N(0,\Sigma)$ where $\Sigma=(\Sigma_{jk})_{j,k=1}^{p}$ and
$\Sigma_{jk}=0.5^{|j-k|}$ for $i=1,\ldots,n$. We set the sample size $n=100$
and the dimension $p=500$ and set the coefficients $\beta=(0.5,0.5)^{\prime}$
and $\gamma=(1,0,\ldots)\in\mathbb{R}^{p-2}$. In Figure 1, we report the
root-$n$ scaled bias along with its standard error bands based on 10,000 Monte
Carlo samples.

Figure 1:  Root-$n$ scaled bias along with its Monte Carlo standard error
bands for Example 1.

From the results in Figure 1, we observe that all three estimators are biased
towards estimating $\beta_{\max}$. Although the oracle estimator
$\widehat{\beta}_{\text{Oracle}}$ is an unbiased estimator of $\beta$, its
maximum is usually not centered around $\beta_{\max}$. In fact, following some
explicit evidence given in Nadarajah and Kotz (2008), this simple plug-in
estimate is usually biased upward. The magnitude of this bias also crucially
depends on the unknown parameters $\beta$ and the covariance structure of
$\widehat{\beta}_{\text{Oracle}}$. Therefore, any inference procedure based on
the naive plug-in estimates of $\beta_{\max}$ cannot be valid. This
discrepancy between the naive plug-in estimate and the true maximum effect is
referred to as the “subgroup selection bias” in the literature (Zöllner and
Pritchard, 2007; Cook et al., 2014; Bornkamp et al., 2017; Guo and He, 2019).
The Lasso and the Refitted estimators are typically biased for $\beta$ due to
regularization (see Wang et al., 2019, for detailed discussion), and they
cannot avoid the selection bias issue either. The inferential framework
proposed in the following sections simultaneously adjusts for the selection
bias and the penalization bias, and it produces a bias-reduced estimate as
well as an asymptotically sharp confidence interval of $\beta_{\max}$.

######  Remark 1 (Construction of $Z_{i}$).

When subgroups are defined similar to RCTs, $Z_{i}$ is composed of the
interaction terms between the treatment indicator variable and the subgroups
indicator variables. Under the Neyman-Rubin causal model (Neyman, 1923; Rubin,
1974), we show in Supplementary Materials (Section C) that such a construction
allows us to interpret $\beta$ as the treatment effects of the considered
subgroups. For now, we require that the considered subgroups do not overlap,
and we show in Section 5.2 how our framework can be naturally extended to
overlapping subgroups. When subgroups are defined as subpopulations of
patients with different treatments, $Z_{i}$ consists of $p_{1}$ indicator
variables each of which represents a different treatment, and $\beta$
represents the corresponding treatment effects. Here, we allow each individual
to receive multiple treatments. A similar model setup has been considered in
Imai et al. (2013); Wang and Ware (2013); Lipkovich et al. (2017).

Notation. Denote $p=p_{1}+p_{2}$ as the dimension of the vector
$(Z_{i}^{\prime},X_{i}^{\prime})^{\prime}$. Define the covariate matrix
$\bm{X}=(X_{1},\ldots,X_{n})^{\prime}\in\mathbb{R}^{n\times p_{2}}$ and the
subgroup design matrix
$\bm{Z}=(Z_{1},\ldots,Z_{n})^{\prime}\in\mathbb{R}^{n\times p_{1}}$. Denote
the maximal of any vector $a$ as $a_{\max}=\max\\{a_{j},j=1,\cdots p\\}$, and
denote a collection of integers from $1$ to $p$ as $[p]$. Suppose $M$ is a
subset of $\\{1,\cdots,p\\}$. Then for any $p$-dimensional vector $a$, $a_{M}$
is defined to be the sub-vector of $a$ indexed by $M$ and
$\bm{X}_{M}=\\{X_{\cdot j},j\in M\\}$, where $X_{\cdot j}$ is the $j$th column
of $\bm{X}$. Define $I_{p}$ to be the $p-$dimensional identity matrix. Define
$I_{M}$ and $I_{Z}$ to be index matrices, whose dimensions are context-
specific and satisfy
$(Z_{i}^{\prime},X_{i}^{\prime})I_{M}^{\prime}=(Z_{i}^{\prime},X_{i,M}^{\prime})$
and $I_{Z}(Z_{i}^{\prime},X_{i,M}^{\prime})^{\prime}=Z_{i}^{\prime}$,
respectively. Define the sample covariance matrix as
$\widehat{\Sigma}_{M}=\frac{1}{n}\sum_{i=1}^{n}(Z_{i}^{\prime},X_{i,M}^{\prime})^{\prime}(Z_{i}^{\prime},X_{i,M}^{\prime})$
and the projection matrix as
$\mathbf{P}_{M}=\bm{X}_{M}\big{(}\bm{X}_{M}^{\prime}\bm{X}_{M}\big{)}^{-1}\bm{X}_{M}^{\prime}$.
Lastly, $Q_{Y}(\tau)$ denotes the $\tau$th quantile of a random variable $Y$.

##  3 Methodology

In this section, we propose two bootstrap-based inferential frameworks for
$\beta_{\max}$. We allow two well-understood debiased estimators to be used as
building blocks: One built on the debiased Lasso (Zhang and Zhang, 2014; Van
de Geer et al., 2014; Belloni et al., 2014) for large $p_{1}$ (Section 3.1),
and the other one built on the repeated data splitting of Wang et al. (2019)
for fixed $p_{1}$ (Section 3.2). While both procedures have appealing
statistical properties, they have different strengths. On the one hand, the
procedure based on the debiased Lasso allows $p_{1}$ to increase with $n$ and
does not require the non-zero coefficients in $\gamma$ to be rather large.
Therefore, it equips researchers with a flexible choice of $Z_{i}$. On the
other hand, the procedure based on the repeated data splitting (R-Split) takes
advantages of a sparse model. When subgroups are highly correlated with the
covariates, it often provides a more efficient estimate of the target
parameter. As this section focuses on the methodology side of each method, we
delay the discussion on their comparison to Section 4.3.

###  3.1 Bootstrap assisted debiased Lasso adjustment for large $p_{1}$

When $p_{1}$ is a large and increases with the sample size $n$, we propose the
following procedure to simultaneous address the selection bias and the
penalization issues:

Step 1.

Construct the debiased estimator $\widehat{b}$ of $\beta$ through (2).

Step 2.

For $j\in[p_{1}]$, calculate the calibration term:

$\widehat{c}_{j}(r)=(1-n^{r-0.5})(\widehat{\beta}_{\max}-\widehat{\beta}_{j,\text{Lasso}}),$

where $r$ is a positive tuning parameter between 0 and 0.5, whose order is
motivated by our theoretical analysis (See Section 4). In practice, we choose
$r$ via cross-validation (see Section 5.1 for implementation).

Step 3.

For $b\leftarrow 1$ to $B$ do

  1. 1.

Generate a consistent bootstrap replicate of $\widehat{b}$, denoted as
$\widehat{b}^{*}$.

  2. 2.

Recalibrate bootstrap statistics via

$\displaystyle
T_{b}^{*}=\underset{j\in[p_{1}]}{\max}(\widehat{b}^{*}_{j}+\widehat{c}_{j}(r))-\widehat{\beta}_{\max}.$

Step 4.

The level-$\alpha$ one-sided confidence interval for $\beta_{\max}$ is
$[\widehat{b}_{\max}-{Q}_{T^{*}_{b}}(\alpha),+\infty)$, and a bias-reduced
estimate for $\beta_{\max}$ is
$\widehat{b}_{\max}-\frac{1}{B}\sum_{b=1}^{B}T_{b}^{*}$.

To cast some insight into the proposed framework, we comment on its
methodological details from three perspectives:

First, to address the penalization bias issue, we estimate $\beta$ via the de-
sparsified Lasso procedure (Zhang and Zhang, 2014; Van de Geer et al., 2014)

$\displaystyle\widehat{b}_{j}=\widehat{\beta}_{j,\text{Lasso}}+\widetilde{Z}_{\cdot
j}^{\prime}\big{(}Y-\bm{Z}\widehat{\beta}_{\text{Lasso}}-\bm{X}\widehat{\gamma}_{\text{Lasso}}\big{)},\quad
j=1,\ldots,p_{1},$ (2)

where $\widetilde{Z}_{\cdot j}$ can be viewed as a standardized residuals of
$Z_{\cdot j}$ after regressing it on $\bm{X}$, and it equals
$\widetilde{Z}_{\cdot j}=\frac{V_{j}}{V_{j}^{\prime}Z_{\cdot j}}$ with
$V_{j}=Z_{\cdot j}-(\bm{Z},\bm{X})_{-j}\widehat{\zeta}_{j}$ and

$\widehat{\zeta}_{j}=\mathrm{\arg\min}\Big{\\{}||Z_{\cdot
j}-(\bm{Z},\bm{X})_{-j}\zeta_{j}||_{2}^{2}/n+\lambda_{n_{j}}||\zeta_{j}||_{1}:\zeta_{j}\in\mathbb{R}^{p-1}\Big{\\}}.$
(3)

Second, in Step 3, our procedure relies on a valid bootstrap approximation of
the debiased estimate $\widehat{b}$. In high dimensions, we adopt the wild
bootstrap approach proposed in Dezeure et al. (2017):

$Y_{i}^{*}=Z_{i}^{\prime}\widehat{\beta}_{\text{Lasso}}+X_{i}^{\prime}\widehat{\gamma}_{\text{Lasso}}+\varepsilon^{*}_{i},\hskip
28.45274pti=1,\dots,n,$ (4)

where $\varepsilon_{i}^{*}=u_{i}\widehat{\varepsilon}_{i}$, and
$\widehat{\varepsilon}_{i}$ is the residual from the Lasso estimate for the
original sample and $u_{i}$ is i.i.d and independent of the data with
$E[u_{i}]=0$, $E[u_{i}^{2}]=1$ and $E[u_{i}^{4}]<\infty$.
$\widehat{b}^{*}=(\widehat{b}^{*}_{1},\dots,\widehat{b}^{*}_{p_{1}})$ is
defined as the debiased lasso estimate for $\beta$ based on the bootstrap
sample $\\{(Y_{i}^{*},Z_{i},X_{i})\\}_{i=1}^{n}$. We note that in (4), instead
of using the debiased Lasso estimate $\widehat{b}$ to generate bootstrap
sample, Dezeure et al. (2017) adopts $\widehat{\beta}_{\text{Lasso}}$ to
ensure the bootstrap sample is indeed generated from a sparse model. In
addition, this wild bootstrap procedure can incorporate heteroscedastic errors
without much increase of the computational cost as the design matrix remains
unchanged across different bootstrap samples.

Under certain regularity conditions, Dezeure et al. (2017) have shown that the
bootstrap procedure in (4) is consistent in the sense that conditional on the
data, the asymptotic distribution of
$\sqrt{n}(\widehat{b}^{*}-\widehat{\beta}_{\text{Lasso}})$ is the same as the
limiting distribution of the debiased Lasso estimator
$\sqrt{n}(\widehat{b}-\beta)$. However, similar to $\widehat{b}_{\max}$ not
centering around $\beta_{\max}$,
$\widehat{b}^{*}_{\max}=\underset{j\in[p_{1}]}{\max}\widehat{b}^{*}_{j}$ is
also not centered around
$\widehat{\beta}_{\max}=\underset{j\in[p_{1}]}{\max}\widehat{\beta}_{j,\text{Lasso}}$,
and usual high dimensional bootstrap procedures do not estimate the selection
bias in $\widehat{b}_{\max}$ correctly. As a result, any inference simply
based on $\widehat{b}^{*}_{\max}$ cannot be valid for making inference on
$\beta_{\max}$ and appropriate adjustments to the bootstrap procedure is
needed for a valid inference on $\beta_{\max}$.

Third, Step 3 constitutes the core of our bootstrap calibration procedure,
which addresses the issue of the selection bias. There, we propose to modify
$\widehat{b}^{*}_{\max}$ as the following,

$\widehat{b}^{*}_{\mathrm{modified};\max}=\underset{j\in[p_{1}]}{\max}(\widehat{b}^{*}_{j}+\widehat{c}_{j}(r)),$
(5)

where
$\widehat{c}_{j}(r)=(1-n^{r-0.5})(\widehat{\beta}_{\max}-\widehat{\beta}_{j,\text{Lasso}})$,
and $r$ is a positive tuning parameter whose theoretical order is
characterized in Assumptions Assumption 7 and Assumption 8. In practice, we
provide a simple cross validation procedure to adaptively choosing $r$ from
data (Section 5.1). In $\widehat{b}^{*}_{\mathrm{modified};\max}$, clearly we
make an adjustment to each debiased estimate $\widehat{b}^{*}_{j}$ by the
amount of $\widehat{c}_{j}(r)$, which measures the distance between
$\widehat{\beta}_{j,\text{Lasso}}$ and the observed maximal regression
coefficient $\widehat{\beta}_{\max}$ from Lasso. The amount of this adjustment
is large if $\widehat{\beta}_{j,\text{Lasso}}$ is small, and is small if
$\widehat{\beta}_{j,\text{Lasso}}$ is large. By adding the correction term
$\widehat{c}_{j}(r)$, we show in Section 4 that the asymptotic distributions
of $\sqrt{n}(\widehat{b}^{*}_{\mathrm{modified};\max}-\widehat{\beta}_{\max})$
and $\sqrt{n}(\widehat{b}_{\max}-\beta_{\max})$ are equivalent, meaning that
the bootstrap modification can successfully adjust for the regularization bias
and selection bias simultaneously.

###  3.2 Bootstrap assisted R-Split adjustment for fixed $p_{1}$

When $\beta$ is a low-dimensional parameter of interest, penalizing $\beta$
may not be necessary. In fact, in this case, inference on $\beta$ is
frequently carried out after a sufficiently small model is selected (Belloni
et al., 2013, 2014). Although any reasonable model selection procedure can be
adopted, we only penalize $\gamma$ since $\beta$ is the target parameter and
should always be kept in the selected model. In this case, given a properly
chosen data-dependent model $\widehat{M}$, the refitted OLS estimator
$\widehat{\beta}_{\text{OLS}}$

$\displaystyle(\widehat{\beta}_{\text{OLS}}^{\prime},\widehat{\gamma}^{\prime}_{\text{OLS}})^{\prime}=\underset{\beta\in\mathbb{R}^{p_{1}},\gamma\in\mathbb{R}^{|\widehat{M}|}}{\arg\min}\left\\{\frac{1}{n}\sum_{i=1}^{n}(Y_{i}-Z_{i}^{\prime}\beta-
X_{i,\widehat{M}}^{\prime}\gamma)^{2}\right\\},$

is a popular choice in practice. However, under the impact of the random model
$\widehat{M}$ entering the estimation process, $\widehat{\beta}_{\text{OLS}}$
is usually biased unless a perfect model is selected. As it is not the main
focus of this article, we only briefly discuss the reason of this issue in
this section and refer interested readers to Wang et al. (2019) for more
examples and simulation studies. Following a similar derivation in Wang et al.
(2019), we decompose $\widehat{\beta}_{\text{OLS}}$ as

$\displaystyle\sqrt{n}(\widehat{\beta}_{\text{OLS}}-\beta)=\underbrace{I_{Z}(\widehat{\Sigma}_{\widehat{M}})^{-1}\cdot\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\begin{pmatrix}Z_{i}\\\
X_{i,\widehat{M}}\end{pmatrix}\varepsilon_{i}}_{=:b_{n1}\text{ (Over-fitting
bias)}}+\underbrace{\vphantom{I_{Z}(\widehat{\Sigma}_{\widehat{M}})^{-1}\cdot\frac{1}{n}\sum_{i=1}^{n}\begin{pmatrix}rZ_{i}\\\
X_{i,\widehat{M}}\end{pmatrix}\varepsilon_{i}}(\bm{Z}^{\prime}(\bm{I}-\bm{P}_{\widehat{M}})\bm{Z}/n)^{-1}\bm{Z}^{\prime}(\bm{I}-\bm{P}_{\widehat{M}})\bm{X}\beta/\sqrt{n}}_{=:b_{n2}\text{
(Under-fitting bias)}}.$

Due to the correlation between $\varepsilon_{i}$ and a data dependent model
$\widehat{M}$, we typically have
$\mathbb{E}(\varepsilon_{i}|X_{i,\widehat{M}})\neq 0$, meaning that the first
term $b_{n1}$ is roughly the mean of $n$ random variables that do not have
mean zero. The second term $b_{n2}$ captures the impact of the unselected
variables in $\widehat{M}$ in estimating $\beta$, and it vanishes whenever the
selected model $\widehat{M}$ covers the support set of $\gamma$. To
simultaneously control the over- and under-fitting biases induced by
regularization procedures used in model selection, we adopt the repeated data
splitting approach (R-Split) proposed by Wang et al. (2019) to estimate
$\beta$. A detailed description of R-Split has been provided in Steps 1 and 2.
Under certain regularity conditions, R-Split estimator, denoted as
$\widetilde{b}$, removes the over- and under-fitting bias, and it converges to
a normal distribution centered around $\beta$ at a root-$n$ rate.

Note that the model selection procedure adopted in Step 1(b) can be any easily
accessible procedure, but ideally, we hope this selection procedure by-pass
the model selection mistakes in a high probability to avoid the risk of under-
fitting. Following the result given in Theorem 1 of Wang et al. (2019), the
smoothed estimator $\widetilde{b}$ satisfies the following linear expansion:

$\displaystyle\widetilde{b}-\beta=\Gamma_{n}\cdot\frac{1}{n}\sum_{i=1}^{n}(Z_{i}^{\prime},X_{i}^{\prime})^{\prime}\varepsilon_{i}+o_{p}(1/\sqrt{n}),$

where $\Gamma_{n}$ is a $p_{1}$ by $p$ dimensional matrix that is independent
of $\varepsilon$, and $\Gamma_{n}$ can be well approximated by its sample
analogy $\widetilde{\Gamma}_{n}$ defined in Step 2 in the previous section.
After generating the residuals
$(\varepsilon^{*}_{1},\ldots,\varepsilon^{*}_{n})$ following (4), such a
linear expansion allows us to quickly generate the bootstrap replicates of
$\widetilde{b}$ without implementing double bootstrap:

$\displaystyle\widetilde{b}^{*}=\widetilde{b}+\widetilde{\Gamma}_{n}\cdot\frac{1}{n}\sum_{i=1}^{n}\begin{pmatrix}Z_{i}\\\
X_{i}\end{pmatrix}\varepsilon^{*}_{i}.$ (6)

We provide the theoretical justification of this bootstrap procedure in the
Supplementary Material Section B.

Similar to (5) in the previous section, with the help of a valid bootstrap
procedure in replicating $\widetilde{b}$, we again modify the bootstrap
statistics by adding
$\widetilde{c}_{j}(r)=(1-n^{r-0.5})(\widetilde{b}_{\max}-\widetilde{b}_{j})$
to have a consistent approximation for the distribution of
$\widetilde{b}_{\max}$. A detailed description of the proposed procedure is
summarized below.

Step 1.

For $b\leftarrow 1$ to $B_{1}$ do

  1. 1.

Randomly split the data $\\{(Y_{i},X_{i},Z_{i})\\}_{i=1}^{n}$ into group
$T_{1}$ of size $n_{1}$ and group $T_{2}$ of size $n_{2}=n-n_{1}$, and let
$v_{bi}=\mathbf{1}_{(i\in T_{2})}$, for $i=1,\cdots,n$.

  2. 2.

Select a model $\widehat{{M}}_{b}$ to predict $Y$ based on $T_{1}$.

  3. 3.

Refit the model with the data in $T_{2}$ to get
$(\widetilde{b}_{b}^{\prime},\widetilde{\gamma}_{b}^{\prime})^{\prime}={\arg\min}\
\sum_{j\in
T_{2}}(Y_{j}-Z_{j}^{\prime}\beta-X^{\prime}_{j,\widehat{{M}}_{b}}{\gamma})^{2}.$

The R-Split estimate is obtained by averaging over $\widetilde{b}_{b}$:

$\displaystyle\widetilde{b}=\frac{1}{B_{1}}\sum_{b=1}^{B_{1}}\widetilde{b}_{b}.$

Step 2.

For $j\in[p_{1}]$, calculate:

$\widetilde{\Gamma}_{n}=\frac{1}{B_{1}}\sum_{b=1}^{B_{1}}I_{z}\left(\frac{1}{n_{1}}\sum_{i=1}^{n}v_{bi}\begin{pmatrix}Z_{i}\\\
X_{i,\widehat{M}}\end{pmatrix}(Z_{i}^{\prime},X^{\prime}_{i,\widehat{M}})^{\prime}\right)^{-1}I_{\widehat{M}},\quad\widetilde{c}_{j}(r)=(1-n^{r-0.5})(\widetilde{b}_{\max}-\widetilde{b}_{j}),$

where $r$ is a positive tuning parameter between 0 to 0.5.

Step 3.

For $b\leftarrow 1$ to $B_{2}$ do

  1. 1.

Generate bootstrap replicate $\widetilde{b}^{*}$ from (6).

  2. 2.

Recalibrate bootstrap statistics via

$\displaystyle
T_{b}^{*}=\underset{j\in[p_{1}]}{\max}(\widetilde{b}^{*}_{j}+\widetilde{c}_{j}(r))-\widetilde{b}_{\max}.$

Step 4.

The level-$\alpha$ one-sided confidence interval for $\beta_{\max}$ is
$[\widetilde{b}_{\max}-{Q}_{T^{*}_{b}}(\alpha),+\infty)$, and a bias-reduced
estimate for $\beta_{\max}$ is
$\widetilde{b}_{\max}-\frac{1}{B_{2}}\sum_{b=1}^{B_{2}}{T^{*}_{b}}$.

##  4 Theoretical investigation and a comparison

###  4.1 Theoretical investigation for bootstrap assisted debiased Lasso
adjustment

In this section, we discuss the theoretical properties of the bootstrap
assisted debiased Lasso adjustment in detail. Define the set of indexes with
the maximal regression coefficient as $H=\\{j\in[p_{1}]:\
\beta_{j}=\max_{k\in[p_{1}]}\beta_{k}\\}$. To establish the asymptotic
validity of the bootstrap assisted debiased Lasso procedure, under the fixed
design, we work under the following assumptions:

######  Assumption 1.

The Lasso estimates satisfy
$||\widehat{\beta}_{\textit{Lasso}}-\beta||_{1}+||\widehat{\gamma}_{\textit{Lasso}}-\gamma||_{1}=o_{p}(1/\sqrt{\log(p)\log(1+p_{1})}).$

######  Assumption 2.

The tuning parameter for Lasso satisfies $\lambda_{n}\asymp\sqrt{\log(p)/n}$.

######  Assumption 3.

The noise variables $\varepsilon_{1},\dots,\varepsilon_{n}$ are mutually
independent, and satisfy (i) $\mathbb{E}(\varepsilon_{i})=0$, (ii)
$\mathbb{E}\sum_{i=1}^{n}\varepsilon_{i}^{2}/n=\sigma_{\varepsilon}^{2}$ and
$\sigma_{i}^{2}=\mathbb{E}\varepsilon_{i}^{2}$ bounded away from zero, and
(iii) $\mathbb{E}|\varepsilon_{i}|^{4}$ is bounded.

######  Assumption 4.

The bootstrapped Lasso estimates satisfy
$||\widehat{\beta}^{*}_{\text{Lasso}}-\widehat{\beta}_{\text{Lasso}}||_{1}+||\widehat{\gamma}^{*}_{\text{Lasso}}-\widehat{\gamma}_{\text{Lasso}}||_{1}=o_{p}(1/\sqrt{\log(p)\log(1+p_{1})})$.

######  Assumption 5.

The treatment variables $Z_{i}$ and the covariates $X_{i}$ are bounded, for
$i=1,\cdots,n$.

######  Assumption 6.

For the de-sparsifed Lasso step, we require (i) $||V_{j}||_{2}^{2}/n$ is
bounded away from zero, (ii) $||V_{j}||_{4}^{4}=o(||V_{j}||_{2}^{4})$,
$j\in[p_{1}]$, and (iii) $\max_{j\in[p_{1}]}||V_{j}||_{\infty}$ is bounded
above, $\log(p_{1})=o(n^{1/7})$.

######  Assumption 7.

The tuning parameter and the Lasso estimate satisfy
$n^{r}||\widehat{\beta}_{\text{Lasso}}-\beta||_{\infty}=$
$o_{p}(1/\sqrt{\log(1+p_{1})})$.

######  Assumption 8.

The tuning parameter satisfies
$\liminf_{n\to\infty}\left(n^{r}(\max_{j\in[p_{1}]}\beta_{j}-\max_{j\notin
H}\beta_{j})-\log(p_{1})\right)=+\infty$.

Assumptions Assumption 1-Assumption 6 guarantee a valid (bootstrap) debiased
Lasso procedure. Dezeure et al. (2017) requires similar conditions and
identify sufficient conditions for Assumptions Assumption 1, Assumption 2,
Assumption 4 and Assumption 6. We refer interested readers to Dezeure et al.
(2017) for a comprehensive discussion of Assumptions Assumption 1-Assumption
6. Assumptions Assumption 7-Assumption 8 are required for a valid bias
correction procedure for $\widehat{b}_{\max}$. Intuitively, they require the
tuning parameter to be neither too large (Assumption Assumption 7) nor too
small (Assumption Assumption 8), so that the modified bootstrap can correct
the penalization bias and selection bias simultaneously. In fact, these two
assumptions are immediate results of Assumptions Assumption 1-Assumption 6
under certain conditions: if $n^{r}=O(\sqrt{\log(p)})$, then Assumption
Assumption 1 directly implies Assumption Assumption 7; if
$\max_{j\in[p_{1}]}\beta_{j}-\max_{j\notin H}\beta_{j}$ is a constant and
$r\geq 1/7$, then Assumption Assumption 6 implies Assumption Assumption 8.
Lastly, when $p_{1}$ is fixed and $r\in(0,0.5)$, Assumptions Assumption
7-Assumption 8 are automatically satisfied.

######  Theorem 1.

Under Assumptions Assumption 1-Assumption 8, the modified bootstrap maximal
treatment effect estimator $\widehat{b}^{*}_{\mathrm{modified};\max}$
satisfies:

$\sup_{c\in\mathbb{R}}|\mathbb{P}(\sqrt{n}(\widehat{b}_{\max}-\beta_{\max})\leq
c)-\mathbb{P}^{*}(\sqrt{n}(\widehat{b}^{*}_{\mathrm{modified};\max}-\widehat{\beta}_{\max})\leq
c)|=o_{p}(1).$

The proof of Theorem 1 is provided in Supplementary Material Section A.2.
Theorem 1 confirms that the proposed one-sided confidence interval is
asymptotically sharp in the sense that the proposed confidence interval
achieves the exact nominal level as the sample size goes to infinity under
mild conditions. In addition, any choice of the tuning parameter satisfying
Assumptions Assumption 7-Assumption 8 guarantees the asymptotic validity.

Besides estimating $\beta_{\max}$, the proposed confidence interval also
serves as an asymptotically sharp prediction interval for the selected
treatment effect, i.e. $\beta_{\widehat{s}}$. We formalize this notion in the
following corollary and its proof can be found in Supplementary Material
Section A.3.

######  Corollary 1 (Selected subgroup with the maximal treatment effect).

Under Assumptions Assumption 1-Assumption 8, we have

$\sup_{c\in\mathbb{R}}|\mathbb{P}(\sqrt{n}(\widehat{b}_{\max}-\beta_{\widehat{s}})\leq
c)-\mathbb{P}^{*}(\sqrt{n}(\widehat{b}^{*}_{\mathrm{modified};\max}-\widehat{\beta}_{\max})\leq
c)|=o_{p}(1).$

###  4.2 Theoretical investigation for bootstrap assisted R-Split adjustment

In this section, we provide the theoretical property of the bootstrap assisted
R-Split adjustment. As we work under the same assumption of the ones discussed
in Wang et al. (2019), we only state the theoretical conclusion to avoid
redundancy.

######  Theorem 2.

Under the assumptions given in Wang et al. (2019), when $p_{1}$ is a fixed
number, the modified bootstrap maximal treatment effect estimator
$\widetilde{b}^{*}_{\mathrm{modified};\max}=\underset{j\in[p_{1}]}{\max}(\widetilde{b}^{*}_{j}+\widetilde{c}_{j}(r))$
satisfies:

$\sup_{c\in\mathbb{R}}|\mathbb{P}(\sqrt{n}(\widetilde{b}_{\max}-\beta_{\max})\leq
c)-\mathbb{P}^{*}(\sqrt{n}(\widetilde{b}^{*}_{\mathrm{modified};\max}-\widetilde{b}_{\max})\leq
c)|=o_{p}(1),$

and

$\sup_{c\in\mathbb{R}}|\mathbb{P}(\sqrt{n}(\widetilde{b}_{\max}-\beta_{\widehat{s}})\leq
c)-\mathbb{P}^{*}(\sqrt{n}(\widetilde{b}^{*}_{\mathrm{modified};\max}-\widetilde{b}_{\max})\leq
c)|=o_{p}(1).$

###  4.3 A power analysis: comparison between debiased Lasso and R-Split
assisted bootstrap calibration

While both procedures provided in Sections 3.1 and 3.2 share similar
statistical guarantees, R-Split has some advantages when $p_{1}$ is small. As
discussed in Wang et al. (2019), when $Z_{i}$ and $X_{i}$ are correlated,
R-Split provides more efficient point estimates than the debiased Lasso method
at sparse models. In this case, we expect that the proposed method built on
R-Split provides shorter confidence intervals and hence is more powerful in
detecting subgroup treatment effect heterogeneity. Moreover, R-Split does not
penalize the target parameter $\beta$, meaning it is more transparent compared
to the debiased Lasso method.

To verify our heuristic claim, we provide a simple simulation study which is
similar to the two-stage model setup from the ones considered in Wang et al.
(2019) and Belloni et al. (2014). We generate random samples from the model
$Y_{i}=0.5+Z_{i}^{\prime}\beta+X_{i}^{\prime}\gamma+\varepsilon_{i},\quad
i=1,\dots,n$, where $\varepsilon_{i}$’s are i.i.d. standard normal random
variables. As for the covariates in the model, we consider two cases: (1) When
$Z_{i}$ is a vector of binary random variables, we generate $Z_{i}$ and
$X_{i}$ from $Z_{ij}=X_{i,2j-1}+X_{i,2j}+\nu_{ij},\quad j=1,\ldots,p_{1}$,
where $\\{\nu_{ij}\\}_{i=1}^{n}$ are i.i.d. Bernoulli random variables with
mean $p$, $\\{(X_{i,2j-1},X_{i,2j})^{\prime}\\}_{i=1}^{n}$ are i.i.d.
Bernoulli random vectors that satisfy
$X_{i,2j-1}+X_{i,2j}+\nu_{ij}\in\\{0,1\\}$ for all $j=1,\ldots,p_{1}$, and
$\\{(X_{i,p_{1}+1},X_{ip})^{\prime}\\}_{i=1}^{n}$ are i.i.d. random vectors
follows a multivariate normal distribution with covariance matrix
$I_{p_{2}-p_{1}}$; (2) When $Z_{i}=(Z_{i1},\ldots,Z_{ip_{1}})^{\prime}$ is a
vector of continuous random variables, we generate $Z_{i}$ and $X_{i}$ from
$Z_{ij}=0.5X_{i,2j+3}+\frac{0.5}{\sqrt{2}}X_{i,2j+4}+\nu_{ij},\quad
j=1,\ldots,p_{1}$, where $\nu_{i}\sim N(0,0.1\cdot I_{p_{1}})$ and $X_{i}\sim
N(0,I_{p_{2}})$. We fix the tuning parameter $r=0.1$ for simplicity. We
compare the power of the proposed bootstrap calibration methods based on
R-Split and the debiased Lasso. For R-Split, we choose the model size via
cross-validation with a minimum model size equals 3. For the debiased Lasso-
based method, we present its power curves with four different choice of
$\lambda$, which is the tuning parameter used to obtain
$\widehat{\beta}_{\text{Lasso}}$. Here, $\lambda.{\text{1se}}$ and
$\lambda.{\text{min}}$ are the default tuning parameters provided by the R
package glmnet, $\lambda_{0}$ is the tuning parameter adopted by the hdi
package. To fully illustrate the impact of tuning parameter in the debiased
Lasso assisted bootstrap calibration, we also include a larger tuning
parameter $1.1\cdot\lambda.{\text{1se}}$.

Figure 2: Power comparison for R-split and the debiased Lasso assisted
bootstrap calibrations. Our simulation results are evaluated through 1,000
Monte Carlo samples.

The results in Figure 2 indicate that both approaches control the type-I error
rate at the nominal level when $\beta_{\max}=0$. When the covariates $X_{i}$
are highly correlated with $Z_{i}$, R-Split assisted bootstrap calibration is
more powerful than the one based on the debiased Lasso. While the performance
of the debiased Lasso assisted bootstrap calibration is sensitive to the
choice of the tuning parameter $\lambda$, R-Split based approach is rather
robust. These findings are in line with our theoretical analysis.

##  5 Practical guide for the bootstrap calibration

###  5.1 Implementation detail and the choice of the tuning parameter $r$

We implement the debiased Lasso-based method in Section 3.1 using the R
package hdi. R-Split assisted bootstrap method in Section 3.2 is implemented
following the recommendation in Wang et al. (2019): we select the tuning
parameter (hence the model selection) in Lasso by cross-validation with a
constraint on the maximum and minimum model sizes. We generate $200$ bootstrap
replications for both methods and set $B_{1}=1000$ for R-Split following the
recommendation in Wang et al. (2019).

1: Suppose $R=\\{r_{1},\dots,r_{m}\\}$ is a set of candidate tuning parameters
with $r_{1}<\dots<r_{m}$ and $m$ is a finite integer. Randomly partition the
data into $v$ equal-sized subsamples;

2: for $l=1,\dots,m$ do

3: for $j=1,\dots,v$ do

4: Basic setup: use the $j$th subsample as the reference data and the rest as
the training data;

5: Bias-reduced estimator: use the training data to obtain the bias-reduced
estimator of the maximum treatment effect,
$\widehat{b}_{\max,\mathrm{reduced},j}(r_{l})$, with $r_{l}$ as the tuning
parameter;

6: for $i=1,\dots,k$ do

7: Calculations on the reference data: use the reference data to calculate the
debiased lasso estimate of the $i$-th covariate of interest,
$\widehat{b}_{i,j}$, and its standard error $\widehat{\sigma}_{i,j}$;

8: Evaluation of accuracy: calculate
$h_{i,j}(r_{l})=(\widehat{b}_{\max,\mathrm{reduced},j}(r_{l})-\widehat{b}_{i,j})^{2}-\widehat{\sigma}_{i,j}^{2}$;

9: end for

10: end for

11: end for

12: The tuning parameter is chosen to be
$\arg\min_{r_{l}}\\{\min_{i\in[k]}[\sum_{j=1}^{j=v}h_{i,j}(r_{l})/v]\\}$.

Algorithm 1 Cross-validated choice of the tuning parameter $r$

As for the tuning parameter $r$, we propose a data-adaptive cross-validated
algorithm to select $r$ (Algorithm 1). Here, the idea is to choose $r$ that
minimizes the mean square error between the proposed bias reduced estimate
$\widehat{b}_{\max,\mathrm{reduced}}$
($\widetilde{b}_{\max,\mathrm{reduced}}$) and $\beta_{\max}$. To make this
possible without knowing the true value of $\beta_{\max}$, we provide an
approximation of the mean square error that can be computed from the data
(Line 8). The justification of this cross validation method for fixed $p_{1}$
can be found in Guo and He (2020). When $p_{1}$ increases, we calibrate the
selected tuning parameter via a simple adjustment:
$r^{*}_{\text{cv}}=\frac{r_{\text{cv}}}{\sqrt{p_{1}/2}}$. Our simulation
results suggest this calibration yields confidence intervals with near optimal
coverage probabilities in finite samples. We leave theoretical justification
for such a choice to future research. In practice, we suggest implementing
Algorithm 1 via a three-fold cross-validation with a candidate set
$R=\\{1/3,1/6,\dots,1/30\\}$.

###  5.2 Overlapping subgroups

When a subgroup is defined as a subpopulation of patients decided by their
baseline characteristics, following the derivation in the Supplementary
Material, Model (1) requires each subject to fall into only one of the
candidate subgroups; in other words we require the considered subgroups to be
non-overlapped, so that the validity of the proposed bootstrap calibration
procedure is guaranteed. Overlapping subgroups obviously induce complex
correlation structure and pose difficulties in modelling, estimating, and
making inference on the (selected) subgroup treatment effects. In this
section, we provide a natural extension of our procedure to accommodate
overlapping subgroups. Our procedure entails following steps:

Step 1.

Separate the $K$ original (possibly overlapping) subgroups
$S_{1},\ldots,S_{K}$ into non-overlapping subgroups
$S_{1}^{*},\ldots,S_{p_{1}}^{*}$. Let $\theta_{j}$ denote the subgroup
treatment effects for the constructed non-overlapping subgroups for
$j=1,\dots,p_{1}$, and let $Z_{i}$ denote the interaction term between the
indicators of the constructed non-overlapping subgroups and the treatment
indicator variable;

Step 2.

Estimate the non-overlapping subgroup treatment effects either by the debiased
Lasso ($\widehat{\theta}$) or R-Split ($\widetilde{\theta}$), and generate
their corresponding bootstrap replications ($\widehat{\theta}^{*}$ and
$\widetilde{\theta}^{*}$) via Model (1). For the debiased lasso procedure, we
also calculate the Lasso estimator ($\widehat{\theta}_{\text{Lasso}}$);

Step 3.

Define a matrix $A\in\mathbb{R}^{K\times p_{1}}$ with

$\displaystyle A_{kj}=\mathbb{P}\big{(}\text{the individual }i\text{ belongs
to subgroup }S_{j}^{*}|\text{the individual }i\text{ belongs to subgroup
}S_{k}\big{)},$

for $k=1,\ldots,K$ and $j=1,\ldots,p_{1}$. In practice, $A$ is either known to
us given population demographic information or could be estimated externally
with some prior knowledge.

Step 4.

Estimate the original subgroup treatment effects $\beta$ by applying a linear
transformation: For the debiased Lasso procedure, we estimate the original
subgroup treatment effects via
$\widehat{\beta}_{\text{Lasso}}=A\widehat{\theta}_{\text{Lasso}}$,
$\widehat{b}=A\widehat{\theta}$, and construct their bootstrap replicates via
$\widehat{b}^{*}=A\widehat{\theta}^{*}$. Similarly, for R-Split, we obtain
$\widetilde{b}=A\widetilde{\theta}$ and
$\widetilde{b}^{*}=A\widetilde{\theta}^{*}$;

Step 5.

Proceed with the proposed bootstrap calibration procedure discussed in Section
3 based on these reconstructed subgroup treatment effects.

To further illustrate the proposed algorithm, we provide some explanation via
a simple toy example. Suppose we consider four candidate subgroups: (1) male
group $S_{1}$, (2) female group $S_{2}$, (3) young adult group $S_{3}$ (18 to
35 years old), and (4) senior group $S_{4}$ (older than 65 years old). Since
these four subgroups are obviously overlapped, the coefficient of $Z_{i}$ in
Model (1) cannot capture the subgroup treatment effect if we naively construct
$Z_{i}$ as the interaction term between the indicators of these four candidate
subgroups and the treatment indicator variable. Step 1 of the above algorithm
says we need to separate these overlapping subgroups into non-overlapping
subgroups: (1) young male adult group $S_{1}^{*}$, (2) senior male group
$S_{2}^{*}$, (3) young female adult group $S_{3}^{*}$, and (4) senior female
group $S_{4}^{*}$. Then in Step 2, we estimate the non-overlapping subgroup
treatment effects via Model (1). In Step 3 and 4, we reconstruct the original
subgroup treatment effects of interest by applying a simple linear
transformation $A\in\mathbb{R}^{4\times 4}$. There, the matrix ${A}$ encodes
the conditional proportion, and its component ${A}_{kj}$ represents the
proportion of individuals in $S_{k}$ who also belong to $S_{j}^{*}$. Take the
first male group for example: ${A}_{11}$ equals the proportion of young male
adults in the male group, ${A}_{12}$ equals the proportion of senior males in
the male group, ${A}_{13}$ and ${A}_{14}$ equal zero as there is no male in
the young female adult group or the senior female group. These proportions are
available to us if we have access to the population demographic information
(such as those from the Census Bureau). Given these reconstructed subgroup
treatment effect estimates via the linear transformation, we then proceed with
the proposed bootstrap calibration procedure as usual. We provide a detailed
implementation of the proposed algorithm in the Supplementary Material
(Section D.1).

Step 4 is the key element of the above algorithm, because it allows us to
estimate the original subgroup treatment effects of interest by applying a
simple linear transformation. As shown in the Supplementary Material (Section
D.3), under Assumption Assumption 9 which indicates the separation is
complete, we have $\beta=A\theta$, meaning that the original true subgroup
treatment effects of interest $\beta$ is connected with the non-overlapping
subgroup treatment effects $\theta$ via the linear transformation. We will
further argue in the Supplementary Material (Section D.2) that such complete
separation required in Assumption Assumption 9 always exists. As this section
focuses on the practical implementation side of our proposal, we leave the
theoretical justification of the proposed algorithm to the Supplementary
Material (Section D.4).

######  Assumption 9.

The separation of $S_{1},\ldots,S_{K}$ into non-overlapping subgroups
$S_{1}^{*},\ldots,S_{p_{1}}^{*}$ is complete in the sense that for any
$j=1,\ldots,p_{1}$ and any $k=1,\ldots,K$, $S^{*}_{j}\cap
S_{k}=S^{*}_{j}\text{ or }\emptyset$.

##  6 Simulation studies

In this section, we consider various simulation designs to demonstrate the
merit of our proposal. The main takeaway from the simulation study is as
follows: When $p_{1}$ is a large number, the debiased Lasso assisted bootstrap
calibration is more robust than the one based on R-Split. When many of the
covariates are correlated with $Z_{i}$ or $p_{1}$ is rather small, R-Split
assisted bootstrap calibration is more preferable to practitioners given its
high detection power.

We generate data from the model:

$\displaystyle
Y_{i}=0.5+Z_{i}^{\prime}\beta+X_{i}^{\prime}\gamma+\varepsilon_{i},\quad
i=1,\dots,n,$

where $n=600$ and $\varepsilon_{i}$’s are i.i.d. $N(0,1)$ random variables. We
consider two cases for $\beta$: (1) heterogeneous case with
$\beta=(0,\ldots,0,1)^{\prime}\in\mathbb{R}^{p_{1}}$, meaning that the
treatment effects differ across different subgroups; and (2) spurious
heterogeneous case with $\beta=(0,\ldots,0,0)^{\prime}\in\mathbb{R}^{p_{1}}$,
meaning that there is no significant subgroup in the population. We set
$\gamma=(1,1,1,1,0,\dots,0)\in\mathds{R}^{p_{2}}$, and we leave the non-sparse
case with $\gamma_{j}=1/j^{2},j=1,\dots,p_{2}$ to Supplementary Materials
(Section E). In all considered simulation designs, we set
$p_{1}\in\\{2,6,20\\}$ and $(n,p_{2})=(600,800)$.

As for the covariate design, we again consider two cases similar to Section
4.3. When $Z_{i}$ is a vector of binary random variables, we generate $Z_{i}$
and $X_{i}$ from

$\displaystyle
Z_{ij}\sim\text{Bernoulli}\Big{(}\frac{\exp(X_{i,2j-1}+X_{i,2j})}{1+\exp(X_{i,2j-1}+X_{i,2j})}\Big{)},\quad
j=1,\ldots,p_{1},$

where $X_{i}\sim N(0,\Sigma)$ with $\Sigma_{ij}=0.5^{|i-j|}$. When
$Z_{i}=(Z_{i1},\ldots,Z_{ip_{1}})^{\prime}$ is a vector of continuous random
variables, we generate $Z_{i}$ and $X_{i}$ from

$\displaystyle
Z_{ij}=0.5X_{i,2j+3}+\frac{0.5}{\sqrt{2}}X_{i,2j+4}+\nu_{ij},\quad
j=1,\ldots,p_{1},$

where $\nu_{i}\sim N(0,I_{p_{1}})$ and $X_{i}\sim N(0,I_{p_{2}})$.

We compare the finite sample performance of the proposed R-Split and the
debiased Lasso assisted bootstrap calibrations with two benchmark methods: (1)
a naive method with no adjustment, which uses the estimated maximal
coefficient along with its point estimates; and (2) the simultaneous method as
discussed in Dezeure et al. (2017) and Fuentes et al. (2018). For R-Split
estimator, following the recommendation in Wang et al. (2019), we choose the
model size via cross-validation with a minimum model size equals 5. As shown
in Section 4.3, the debiased Lasso based bootstrap calibration is sensitive to
the choice of $\lambda$. We report the results when it has the highest
detection power: we choose $\lambda=1.1\cdot\lambda_{\text{1se}}$ in the
continuous case and $\lambda=\lambda_{0}$ in the binary case. We report the
coverage probability, the $\sqrt{n}$ scaled Monte Carlo bias along with their
standard errors based on 500 Monte Carlo samples in Table 1.

From the results in Table 1, we observe that under spurious heterogeneity,
regardless of whether $Z_{i}$ are binary or continuous, the method with no
adjustment is clearly biased and under covered–especially when $p_{1}$ is
rather large. Because the simultaneous method provides conservative confidence
intervals, its testing power can be compromised in practice (also see real
data analysis in next section). In contrast, the proposed bootstrap-assisted
debiased Lasso and R-split methods perform well regardless of the underlying
data generating processes.

When comparing the proposed method based on R-Split with the one based on the
debiased Lasso, we observe that the procedure based on the debiased Lasso is
more robust than the one based on R-Split when $p_{1}$ is large. We conjecture
that the reason is twofold. First, a large $p_{1}$ increases the chance of
getting a singular refitting covariance matrix, as a result, R-Split can be
numerically unstable if some splits do not produce feasible refitted OLS
estimates. Second, the R-Split estimator may have an intractable asymptotic
distribution when $p_{1}$ increases with $n$. In contrast, the debiased Lasso
estimators are uniformly normally distributed for all components in $\beta$,
suggesting that the resulting confidence intervals and tests based on the
debiased Lasso are not sensitive to the change of $p_{1}$ as demonstrated in
Section 4.1.

Table 1: Simulation results for Section 6

|  | Debiased Lasso |  | Repeated Data Splitting |
---|---|---|---|---|---
|  | Boot-Calibrated | No adjustment | Simultaneous |  | Boot-Calibrated | No adjustment | Simultaneous |
| $Z_{i}$’s are binary random variables, and $\beta=(0,\dots,0,1)\in\mathbb{R}^{p_{1}}$ (heterogeneity) |
$p_{1}=2$ | Cover | 0.95(0.01) | 0.92(0.02) | 0.98(0.01) |  | 0.96(0.01) | 0.95(0.01) | 0.97(0.01) |
|  $\sqrt{n}$Bias | -0.09(0.00) | 0.22(0.00) | -1.61(0.00) |  | -0.08(0.00) | -0.09(0.00) | -1.02(0.00) |
$p_{1}=6$ | Cover | 0.96(0.01) | 0.94(0.02) | 0.99(0.00) |  | 0.95(0.01) | 0.94(0.01) | 0.99(0.00) |
|  $\sqrt{n}$Bias | -0.07(0.00) | 0.19(0.00) | -2.35(0.00) |  | -0.10(0.00) | -0.15(0.00) | -1.91(0.00) |
$p_{1}=20$ | Cover | 0.96(0.01) | 0.96(0.01) | 0.99(0.00) |  | 0.93(0.01) | 0.93(0.01) | 0.99(0.00) |
|  $\sqrt{n}$Bias | -0.04(0.00) | 0.13(0.00) | -3.16(0.00) |  | -0.12(0.00) | -0.20(0.00) | -2.69(0.00) |
|  $Z_{i}$’s are binary random variables, and $\beta=(0,\dots,0)\in\mathbb{R}^{p_{1}}$ (spurious heterogeneity) |
$p_{1}=2$ | Cover | 0.96(0.01) | 0.79(0.03) | 0.98(0.01) |  | 0.95(0.01) | 0.87(0.02) | 0.95(0.02) |
|  $\sqrt{n}$Bias | -0.02(0.00) | 1.18(0.00) | -0.15(0.00) |  | 0.04(0.00) | 0.91(0.00) | 0.10(0.00) |
$p_{1}=6$ | Cover | 0.96(0.01) | 0.69(0.03) | 0.97(0.01) |  | 0.93(0.01) | 0.77(0.03) | 0.97(0.01) |
|  $\sqrt{n}$Bias | -0.03(0.00) | 1.91(0.00) | -0.25(0.00) |  | 0.14(0.00) | 1.69(0.00) | 0.15(0.00) |
$p_{1}=20$ | Cover | 0.94(0.01) | 0.30(0.03) | 0.97(0.01) |  | 0.91(0.01) | 0.43(0.03) | 0.98(0.01) |
|  $\sqrt{n}$Bias | -0.06(0.00) | 2.72(0.00) | -0.30(0.00) |  | 0.17(0.00) | 2.52(0.00) | 0.18(0.00) |
|  $Z_{i}$’s are continuous random variables, and $\beta=(0,\dots,0,1)\in\mathbb{R}^{p_{1}}$ (heterogeneity) |
$p_{1}=2$ | Cover | 0.96(0.01) | 0.90(0.01) | 0.98(0.01) |  | 0.96(0.01) | 0.94(0.01) | 0.98(0.00) |
|  $\sqrt{n}$Bias | -0.07(0.00) | 0.12(0.00) | -0.43(0.00) |  | -0.06(0.00) | -0.09(0.00) | -0.28(0.00) |
$p_{1}=6$ | Cover | 0.96(0.01) | 0.91(0.01) | 0.99(0.00) |  | 0.96(0.01) | 0.93(0.01) | 0.99(0.00) |
|  $\sqrt{n}$Bias | -0.05(0.00) | 0.09(0.00) | -0.67(0.00) |  | -0.08(0.00) | -0.11(0.00) | -0.72(0.00) |
$p_{1}=20$ | Cover | 0.96(0.01) | 0.94(0.01) | 0.99(0.00) |  | 0.94(0.01) | 0.92(0.01) | 0.99(0.00) |
|  $\sqrt{n}$Bias | -0.02(0.00) | 0.07(0.00) | -0.96(0.00) |  | -0.09(0.00) | -0.12(0.00) | -1.03(0.00) |
| $Z_{i}$’s are continuous random variables, and $\beta=(0,\dots,0)\in\mathbb{R}^{p_{1}}$ (spurious heterogeneity) |
$p_{1}=2$ | Cover | 0.94(0.01) | 0.80(0.01) | 0.95(0.01) |  | 0.96(0.01) | 0.89(0.02) | 0.97(0.02) |
|  $\sqrt{n}$Bias | 0.06(0.00) | 1.06(0.00) | 0.02(0.00) |  | 0.01(0.00) | 0.24(0.00) | -0.16(0.00) |
$p_{1}=6$ | Cover | 0.95(0.01) | 0.84(0.02) | 0.96(0.01) |  | 0.95(0.01) | 0.77(0.03) | 0.97(0.01) |
|  $\sqrt{n}$Bias | 0.04(0.00) | 0.77(0.00) | 0.05(0.00) |  | 0.03(0.00) | 0.72(0.00) | -0.27(0.00) |
$p_{1}=20$ | Cover | 0.95(0.01) | 0.91(0.02) | 0.97(0.01) |  | 0.93(0.01) | 0.70(0.03) | 0.98(0.01) |
|  $\sqrt{n}$Bias | 0.02(0.00) | 0.16(0.00) | 0.08(0.00) |  | 0.08(0.00) | 1.10(0.00) | -0.38(0.00) |

  * •

Note: “Cover” is the empirical coverage of the 95% lower bound for
$\beta_{\max}$ and “ $\sqrt{n}$Bias ” captures the root-$n$ scaled Monte Carlo
bias for estimating $\beta_{\max}$.

##  7 Real data analysis

In precision medicine to treat high blood pressure, lifestyle change is
frequently recommended by doctors to prevent mild hypertension and to reduce
the dose levels of drugs needed to control hypertension (Whelton et al.,
2002). Such lifestyle factors include, but are not limited to, sodium intake,
smoking, dietary patterns, alcohol use, and physical inactivity. In addition,
blood pressure is also highly heritable, and to date, more than 600 risk loci
have been identified (Evangelou et al., 2018). So far, there has been no
unified answer as to (1) which genetic-risk subgroup has the most positive
response to a particular lifestyle change (e.g., smoking cessation status),
and (2) when subgroups are defined by different lifestyle factors, which
subgroup has the best blood pressure control. In response to these questions,
the identified risk loci need to be adjusted in the model-building process to
avoid potential confounding issues. Our framework, which incorporates
lifestyle factors as well as those high-dimensional risk loci into covariates,
provides natural solutions to the questions raised above.

To demonstrate the effectiveness of our proposed method on controlling the
selection and regularization bias and to cast some insights into the effects
of lifestyle change on lowering systolic blood pressure, we perform cross-
sectional observational studies from UK Biobank resources. The lifestyle
changes we consider include smoking cessation status, three physical activity
modifications, alcohol intake, and 12 diet changes.

The UK Biobank study is a long-term observational cohort study that recruited
about $500,000$ individuals aged between 40 and 69 in the United Kingdom. We
choose 285,582 unrelated white British individuals, for whom systolic blood
pressure, 17 lifestyle covariates, and genotype data are available. For each
individual, we calculate the genetic risk scores (GRS). GRS are calculated as
the weighted sum of the number of risk alleles, where the risk alleles and
their weights are defined by Evangelou et al. (2018). GRS profile individual-
level risks of disease at birth. We focus on individuals in either high (top
0.2% of GRS) or low (bottom 0.2% of GRS) genetic risk groups, resulting in 857
individuals in our final analyses. We choose 0.2% to make sure the difference
between the two subgroups is significant and meaningful: the mean systolic
blood pressure of the high-risk group is 144.3 mm Hg, which is 16.7 mm Hg
higher than that of the low-risk group. Because blood pressure is highly
heritable, 612 risk genetic variants are adjusted in our models. These genetic
variants are obtained from the GWAS Catalog (Buniello et al., 2019), with
minor allele frequency $>0.05$, and have passed standard quality control
steps. Several additional covariates, including lifestyle variables, BMI, age,
gender, top 40 principal components of genotype matrix, and genotype array,
are adjusted as well. The response is systolic blood pressure. Our empirical
comparisons are provided in Table 2. As the simultaneous method and the method
with no adjustment provide similar results, we only provide their results
based on the debiased Lasso adjustment.

For the two-subgroup comparison, we consider low- and high-genetic risk
subgroups and for whom the treatment (or lifestyle change) is smoking
cessation. While the result from the naive method indicates that blood
pressure control among high-genetic risk subgroup individuals is particularly
reactive to smoking, the proposed methods suggest such a discovery may not be
credible as the lower confidence bounds are both below zero. While both
smoking and blood pressure are major risk factors for cardiovascular disease
and some spurious relationships have been reported (Groppelli et al., 1992),
current belief is that smoking has no independent effect on blood pressure
control (Primatesta et al., 2001) and that smoking cessation does not lower
blood pressure (Williams et al., 2018). Our results not only align with the
current scientific belief, but also highlight the necessity of bias correction
to avoid spurious conclusions drawn from the naive method.

For the multiple-subgroup comparison, we compare the subgroup treatment
effects of 17 lifestyle changes with a focus on the high-generic risk
individuals. From Table 2, we observe that the proposed method and the naive
method reach the same conclusion: the identified best treatment–moderate
physical activity–shows a significant effect in controlling high blood
pressure. Moderate physical activity in our analysis is defined as two or more
days having 10 of minutes or more per week of moderate physical activities
like carrying light loads or cycling at normal pace.

In fact, this is a well-grounded conclusion (Williams et al., 2018). Current
guidelines recommend regular physical activity for preventing hypertension
(Williams et al., 2018) and several randomized controlled trials have further
confirmed the favorable effects of physical activity on reducing blood
pressure (Cornelissen and Smart, 2013). While the effect of physical activity
seems small,  in population level, an systolic blood pressure decrease of 1 mm
Hg decreases the risk of stroke by 5% (Lawes et al., 2004). In both cases,
because the simultaneous method is typically overly conservative and tends to
lose power due to widened confidence intervals, we see that it is unable to
identify any significant effect.

Table 2:  Results for real data analysis

| Debiased | R-Split | No adjustment | Simultaneous
---|---|---|---|---
| Two subgroup comparison
Lower bound | -2.29 | -1.20 | 1.27 | -3.78
Point estimate | -0.13 | 1.17 | 3.26 | -0.82
| Multiple treatment comparison
Lower bound | 0.15 | 0.03 | 1.96 | -0.38
Point estimate | 1.75 | 1.51 | 4.44 | 1.26

  * •

Note: The unit of the above results is mm Hg. “Lower bound” is the lower bound
of the 95%-coverage confidence interval.

##  8 Conclusion

In the presence of high-dimensional covariates from observational data, when a
seemingly promising subgroup has been identified from the data, a naive
estimation and inference for the identified subgroup that ignores the
selection process can lead to biased and overly optimistic conclusions. The
salient point of the present paper is that appropriate statistical analysis of
the post-hoc identified subgroup treatment effect must take the selection
process into account. We propose two bootstrap assisted debiasing procedures
to make valid inferences on the maximal or the selected subgroup treatment
effect. The two proposed methods are not only easy to implement but also
simultaneously control the regularization bias and the selection bias. The
resulting statistical inference is asymptotically sharp. They have different
strengths: the R-Split assisted bootstrap calibration has higher detection
power when the number of subgroups is small, while the debiased Lasso assisted
bootstrap calibration tends to be more robust when we want to investigate
multiple subgroups. Since the proposed approaches have their own strengths, it
is worth exploring the merits of these two approaches and further combine
them. We leave this to our future research.

#### Software and reproducibility

R code for the proposed procedures can be found in the package
“debiased.subgroup” that is publicly available at
https://github.com/WaverlyWei/debiased.subgroup. Simulation examples can be
reproduced by running examples in the R package.

#### Acknowledgement

We thank the individuals involved in the UK Biobank for their participation
and the research teams for their work on collecting, processing, and sharing
these datasets. This research has been conducted using the UK Biobank Resource
(application number 48240), subject to a data transfer agreement. We are also
grateful for the helpful discussions with Peng Ding, Avi Feller, Xuming He,
Xinwei Ma, Gongjun Xu, and the paticipants of the UC Berkeley Causal Inference
Reading group.

## References

  * Alosh et al. (2017) Alosh, M., Huque, M. F., Bretz, F., and D’Agostino Sr, R. B. (2017).  Tutorial on statistical considerations on subgroup analysis in confirmatory clinical trials.  Statistics in medicine, 36(8):1334–1360.
  * Belloni et al. (2013) Belloni, A., Chernozhukov, V., et al. (2013).  Least squares after model selection in high-dimensional sparse models.  Bernoulli, 19(2):521–547.
  * Belloni et al. (2014) Belloni, A., Chernozhukov, V., and Hansen, C. (2014).  Inference on treatment effects after selection among high-dimensional controls.  The Review of Economic Studies, 81(2):608–650.
  * Bornkamp et al. (2017) Bornkamp, B., Ohlssen, D., Magnusson, B. P., and Schmidli, H. (2017).  Model averaging for treatment effect estimation in subgroups.  Pharmaceutical statistics, 16(2):133–142.
  * Buniello et al. (2019) Buniello, A., MacArthur, J. A. L., Cerezo, M., Harris, L. W., Hayhurst, J., Malangone, C., McMahon, A., Morales, J., Mountjoy, E., Sollis, E., et al. (2019).  The nhgri-ebi gwas catalog of published genome-wide association studies, targeted arrays and summary statistics 2019.  Nucleic Acids Research, 47(D1):D1005–D1012.
  * Cook et al. (2014) Cook, D., Brown, D., Alexander, R., March, R., Morgan, P., Satterthwaite, G., and Pangalos, M. N. (2014).  Lessons learned from the fate of astrazeneca’s drug pipeline: a five-dimensional framework.  Nature reviews Drug discovery, 13(6):419–431.
  * Cornelissen and Smart (2013) Cornelissen, V. A. and Smart, N. A. (2013).  Exercise training for blood pressure: a systematic review and meta-analysis.  Journal of the American heart association, 2(1):e004473.
  * Dezeure et al. (2017) Dezeure, R., Bühlmann, P., and Zhang, C.-H. (2017).  High-dimensional simultaneous inference with the bootstrap.  TEST, 26(4):685–719.
  * Evangelou et al. (2018) Evangelou, E., Warren, H. R., Mosen-Ansorena, D., Mifsud, B., Pazoki, R., Gao, H., Ntritsos, G., Dimou, N., Cabrera, C. P., Karaman, I., et al. (2018).  Genetic analysis of over 1 million people identifies 535 new loci associated with blood pressure traits.  Nature Genetics, 50(10):1412–1425.
  * Fan et al. (2017) Fan, A., Song, R., and Lu, W. (2017).  Change-plane analysis for subgroup detection and sample size calculation.  Journal of the American Statistical Association, 112(518):769–778.
  * Fan and Li (2001) Fan, J. and Li, R. (2001).  Variable selection via nonconcave penalized likelihood and its oracle properties.  Journal of the American statistical Association, 96(456):1348–1360.
  * Friede et al. (2012) Friede, T., Parsons, N., and Stallard, N. (2012).  A conditional error function approach for subgroup selection in adaptive clinical trials.  Statistics in Medicine, 31(30):4309–4320.
  * Fuentes et al. (2018) Fuentes, C., Casella, G., Wells, M. T., et al. (2018).  Confidence intervals for the means of the selected populations.  Electronic Journal of Statistics, 12(1):58–79.
  * Groppelli et al. (1992) Groppelli, A., Giorgi, D., Omboni, S., Parati, G., and Mancia, G. (1992).  Persistent blood pressure increase induced by heavy smoking.  Journal of Hypertension, 10(5):495–499.
  * Guo and He (2019) Guo, X. and He, X. (2019).  Inference on the best selected subgroup with a case study of monet1 trial.  preprint.
  * Guo and He (2020) Guo, X. and He, X. (2020).  Inference on selected subgroups in clinical trials.  Journal of the American Statistical Association, (just-accepted):1–18.
  * Hall and Miller (2010) Hall, P. and Miller, H. (2010).  Bootstrap confidence intervals and hypothesis tests for extrema of parameters.  Biometrika, 97(4):881–892.
  * Hébert et al. (2002) Hébert, P. C., Cook, D. J., Wells, G., and Marshall, J. (2002).  The design of randomized clinical trials in critically ill patients.  Chest, 121(4):1290–1300.
  * Imai et al. (2013) Imai, K., Ratkovic, M., et al. (2013).  Estimating treatment effect heterogeneity in randomized program evaluation.  The Annals of Applied Statistics, 7(1):443–470.
  * Izem et al. (2020) Izem, R., Liao, J., Hu, M., Wei, Y., Akhtar, S., Wernecke, M., MaCurdy, T. E., Kelman, J., and Graham, D. J. (2020).  Comparison of propensity score methods for pre-specified subgroup analysis with survival data.  Journal of Biopharmaceutical Statistics, pages 1–18.
  * Jenkins et al. (2011) Jenkins, M., Stone, A., and Jennison, C. (2011).  An adaptive seamless phase ii/iii design for oncology trials with subpopulation selection using correlated survival endpoints.  Pharmaceutical statistics, 10(4):347–356.
  * Kubota et al. (2014) Kubota, K., Ichinose, Y., Scagliotti, G., Spigel, D., Kim, J., Shinkai, T., Takeda, K., Kim, S.-W., Hsia, T.-C., Li, R., et al. (2014).  Phase iii study (monet1) of motesanib plus carboplatin/paclitaxel in patients with advanced nonsquamous nonsmall-cell lung cancer (nsclc): Asian subgroup analysis.  Annals of oncology, 25(2):529–536.
  * Lawes et al. (2004) Lawes, C. M., Bennett, D. A., Feigin, V. L., and Rodgers, A. (2004).  Blood pressure and stroke: an overview of published reviews.  Stroke, 35(3):776–785.
  * Lee et al. (2016) Lee, J. D., Sun, D. L., Sun, Y., Taylor, J. E., et al. (2016).  Exact post-selection inference, with application to the lasso.  The Annals of Statistics, 44(3):907–927.
  * Lipkovich et al. (2017) Lipkovich, I., Dmitrienko, A., and B D’Agostino Sr, R. (2017).  Tutorial in biostatistics: data-driven subgroup identification and analysis in clinical trials.  Statistics in medicine, 36(1):136–196.
  * Nadarajah and Kotz (2008) Nadarajah, S. and Kotz, S. (2008).  Exact distribution of the max/min of two gaussian random variables.  IEEE Transactions on very large scale integration (VLSI) systems, 16(2):210–212.
  * Neyman (1923) Neyman, J. (1923).  On the application of probability theory to agricultural experiments. essay on principles. section 9.(tlanslated and edited by dm dabrowska and tp speed, statistical science (1990), 5, 465-480).  Annals of Agricultural Sciences, 10:1–51.
  * Petticrew et al. (2012) Petticrew, M., Tugwell, P., Kristjansson, E., Oliver, S., Ueffing, E., and Welch, V. (2012).  Damned if you do, damned if you don’t: subgroup analysis and equity.  J Epidemiol Community Health, 66(1):95–98.
  * Primatesta et al. (2001) Primatesta, P., Falaschetti, E., Gupta, S., Marmot, M. G., and Poulter, N. R. (2001).  Association between smoking and blood pressure: evidence from the health survey for england.  Hypertension, 37(2):187–193.
  * Rubin (1974) Rubin, D. B. (1974).  Estimating causal effects of treatments in randomized and nonrandomized studies.  Journal of educational Psychology, 66(5):688.
  * Shen and He (2015) Shen, J. and He, X. (2015).  Inference for subgroup analysis with a structured logistic-normal mixture model.  Journal of the American Statistical Association, 110(509):303–312.
  * Stallard et al. (2008) Stallard, N., Todd, S., and Whitehead, J. (2008).  Estimation following selection of the largest of two normal means.  Journal of Statistical Planning and Inference, 138(6):1629–1638.
  * Su et al. (2009) Su, X., Tsai, C.-L., Wang, H., Nickerson, D. M., and Li, B. (2009).  Subgroup analysis via recursive partitioning.  Journal of Machine Learning Research, 10(Feb):141–158.
  * Tian et al. (2018) Tian, X., Taylor, J., et al. (2018).  Selective inference with a randomized response.  The Annals of Statistics, 46(2):679–710.
  * Tibshirani (1996) Tibshirani, R. (1996).  Regression shrinkage and selection via the lasso.  Journal of the Royal Statistical Society. Series B (Methodological), pages 267–288.
  * Van de Geer et al. (2014) Van de Geer, S., Bühlmann, P., Ritov, Y., Dezeure, R., et al. (2014).  On asymptotically optimal confidence regions and tests for high-dimensional models.  The Annals of Statistics, 42(3):1166–1202.
  * Wager and Athey (2017) Wager, S. and Athey, S. (2017).  Estimation and inference of heterogeneous treatment effects using random forests.  Journal of the American Statistical Association, (just-accepted).
  * Wang et al. (2019) Wang, J., He, X., and Xu, G. (2019).  Debiased inference on treatment effect in a high dimensional model.  Journal of the American Statistical Association, (just-accepted):1–000.
  * Wang and Ware (2013) Wang, R. and Ware, J. H. (2013).  Detecting moderator effects using subgroup analyses.  Prevention Science, 14(2):111–120.
  * Whelton et al. (2002) Whelton, P. K., He, J., Appel, L. J., Cutler, J. A., Havas, S., Kotchen, T. A., Roccella, E. J., Stout, R., Vallbona, C., Winston, M. C., et al. (2002).  Primary prevention of hypertension: clinical and public health advisory from the national high blood pressure education program.  Jama, 288(15):1882–1888.
  * Williams et al. (2018) Williams, B., Mancia, G., Spiering, W., Agabiti Rosei, E., Azizi, M., Burnier, M., Clement, D. L., Coca, A., De Simone, G., Dominiczak, A., et al. (2018).  2018 esc/esh guidelines for the management of arterial hypertension: The task force for the management of arterial hypertension of the european society of cardiology (esc) and the european society of hypertension (esh).  European heart journal, 39(33):3021–3104.
  * Yang et al. (2020) Yang, S., Lorenzi, E., Papadogeorgou, G., Wojdyla, D. M., Li, F., and Thomas, L. E. (2020).  Propensity score weighting for causal subgroup analysis.  arXiv preprint arXiv:2010.02121.
  * Zhang and Zhang (2014) Zhang, C.-H. and Zhang, S. S. (2014).  Confidence intervals for low dimensional parameters in high dimensional linear models.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 76(1):217–242.
  * Zhang and Cheng (2017) Zhang, X. and Cheng, G. (2017).  Simultaneous inference for high-dimensional linear models.  Journal of the American Statistical Association, 112(518):757–768.
  * Zöllner and Pritchard (2007) Zöllner, S. and Pritchard, J. K. (2007).  Overcoming the winner’s curse: estimating penetrance parameters from case-control data.  The American Journal of Human Genetics, 80(4):605–615.