# Assessing Heterogeneous Risk of Type II Diabetes Associated with Statin
Usage: Evidence from Electronic Health Record Data

Xinzhou Guo1 These authors contributed equally to this work and are
alphabetically ordered. Waverly Wei211footnotemark: 1 Molei Liu3 Tianxi Cai3
Chong Wu4 Jingshen Wang2 Correspondence<EMAIL_ADDRESS>
(1Department of Mathematics, Hong Kong University of Science and Technology
2Division of Biostatistics, UC Berkeley
3Department of Biostatistics, Harvard University
4Department of Statistics, Florida State University
)

###### Abstract

There have been increased concerns that the use of statins, one of the most
commonly prescribed drugs for treating coronary artery disease, is potentially
associated with the increased risk of new-onset type II diabetes (T2D).
However, because existing clinical studies with limited sample sizes often
suffer from selection bias issues, there is no robust evidence supporting as
to whether and what kind of populations are indeed vulnerable for developing
T2D after taking statins. In this case study, building on the biobank and
electronic health record data in the Partner Health System, we introduce a new
data analysis pipeline from a biological perspective and a novel statistical
methodology that address the limitations in existing studies to: (i)
systematically examine heterogeneous treatment effects of stain use on T2D
risk, (ii) uncover which patient subgroup is most vulnerable to T2D after
taking statins, and (iii) assess the replicability and statistical
significance of the most vulnerable subgroup via bootstrap calibration. Our
proposed bootstrap calibration approach delivers asymptotically sharp
confidence intervals and debiased estimates for the treatment effect of the
most vulnerable subgroup in the presence of possibly high-dimensional
covariates. By implementing our proposed approach, we find that females with
high T2D genetic risk at baseline are indeed at high risk of developing T2D
due to statin use, which provides evidences to support future clinical
decisions with respect to statin use.
Keywords: Bootstrap; Causal Inference; Debiased Inference; Subgroup Analysis;
Precision Medicine

##  1 Introduction

###  1.1 Motivation and objectives

Coronary artery disease (CAD), a disease affecting the function of heart, is
the leading cause of deaths worldwide (Skourtis et al., 2020). Over the past
decades, efforts have been made in developing effective and safe drugs in
preventing and treating CAD (Povsic et al., 2017). Among those novel agents,
statins are perhaps the most commonly prescribed drugs due to their clear
benefits in reducing the level of low–density lipoprotein (LDL) and
subsequently lowering CAD risks through 3-hydroxy-3-methylglutaryl-coenzyme A
reductase (HMGCR) inhibition (Nissen et al., 2005). Despite their clear
benefits in reducing CAD risks, the use of statins is potentially associated
with the increased risk of new-onset type II diabetes (T2D) (Waters et al.,
2013; Macedo et al., 2014; Mansi et al., 2015).

Although many studies have been conducted to investigate the potential side
effects of statins in developing T2D, to date, there is still no robust
evidence as to whether and on what kind of population statin use increases the
risk of T2D. Take some frequently cited studies as examples. Rajpathak et al.
(2009) find through meta-analysis that there is a small increase in T2D risk
111Relative risk (RR): 1.13, 95% CI [1.03, 1.23] associated with the use of
statins, but this association is no longer significant after including the
results from the WOSCOPS trial (Packard et al., 1998)–the first study
investigating the association between T2D and statins. Recent studies also
suggest that the effect of statins on T2D risk might be heterogeneous across
different sub-populations and be more pronounced in certain subgroups defined
by gender and baseline T2D risk (Mora et al., 2010; Goodarzi et al., 2013).
Nevertheless, existing studies may lead to untrustworthy findings in
subgroups, because their statistical analyses either are conducted under
randomized controlled trials with limited sample sizes whose results might not
be generalizable beyond study population, or do not adjust for multiple
comparisons issue when several subgroups are under consideration (Waters et
al., 2013).

In this case study, leveraging Partner Health System (PHS) biobank and
electronic health record (EHR) data, we conduct subgroup analyses from a novel
biological perspective aiming for the following three objectives: (i)
designing a rigorous causal framework that systematically examines the
(heterogeneous) causal effects of stain use on T2D risk in observational data,
(ii) uncovering which patient subgroup is most vulnerable for developing T2D
after taking statins, and (iii) assessing the replicability and statistical
significance of the most vulnerable subgroup (identified from the data) via a
bootstrap calibration procedure.

###  1.2 Overview of research methods and findings

From a study design perspective, to establish valid causal conclusions, our
research methods not only circumvent some common issues in randomized
controlled trials (RCTs), but also alleviate the concern of unmeasured
confounding bias issues in observation studies.

On the one hand, while existing studies often investigate the adverse effects
of stains on T2D risk in RCTs with limited sample sizes, our study design
leverages the large PHS biobank with linked EHR data, providing robust
evidence for assessing the heterogeneous adverse effect of statin use from
observational data. Concretely, we extract and link genotype information and
diagnostics from consented subjects in the Partner Health System(PHS) biobank
and EHR data respectively, which leads to an EHR virtual cohort of $17,023$
subjects, a much larger cohort than those from usual clinical trials, and
$337$ features; see Section 2.3 for detailed descriptions. Compared to study
cohorts enrolled in RCTs, our study cohort can be a more representative sample
of the general population.

On the other hand, because causal conclusions derived from observational
studies can be susceptible to unmeasured confounding and reverse causation
bias issues, to alleviate these concerns, our study design adopts a randomly
inherited single nucleotide polymorphism (SNP), rs12916-T, as a surrogate
treatment variable of statin use. rs12916-T is a reliable surrogate treatment
variable because it resides in the HMGCR gene encoding the drug target of
statins, and has been recently used as an unbiased, unconfounded proxy for
pharmacological action on the target of statins (i.e., HMG-CoA reductase
inhibition) (Swerdlow et al., 2015; Consortium et al., 2012). Furthermore,
adopting a randomly inherited genetic variant at conception as a surrogate for
statin usage in EHRs allows us to establish a clear temporal precedence
between the treatment and T2D onset (which is a prerequisite to concluding
causality, see Holland (1986) for example), avoiding the potential issue of
reverse causation. Lastly, because the surrogate treatment variable is
naturally inherited, variables observed after birth are independent of
rs12916-T and can at most be mediators that belong to a different causal
pathway. We shall discuss the reasoning of using rs12916-T in more detail in
Section 2.1.

Leveraging the above study design, we further conduct subgroup analyses to
assess if subjects carrying the rs12916-T allele (i.e. taking statins) have
heterogeneous risks at developing T2D, and, in particular, to what extent the
most vulnerable subgroup suffers from statin use. Inspired by the study
designs adopted in Mora et al. (2010) and Wang and Ware (2013), we divide the
study cohort into six subgroups based on gender and baseline T2D genetic risk
profiles (measured by the number of T2D risk alleles of variants rs35011184-A
and rs1800961-T each individual carries). While numerous methods have been
proposed for identifying subgroups (Lipkovich et al., 2011; Ma and Huang,
2017) and testing subgroup homogeneity (Shen and He, 2015; Fan et al., 2017),
in our case study, we focus on post-hoc inference on the most vulnerable
subgroup.

The limitation of post-hoc inference on the most vulnerable subgroup is well
recognized. Due to subgroup selection bias induced by multiple comparisons,
such post-hoc inference often leads to false positive results (Thomas and
Bornkamp, 2017). Although several attempts have been made to address the
subgroup selection bias issue, existing procedures are either poorly grounded
(Stallard et al., 2008; Rosenkranz, 2016) or tend to be conservative (Hall and
Miller, 2010; Fuentes et al., 2018), as the latter is typically built on
simultaneous inference which aims to control the family-wise error rate for
all candidate subgroups. These conservative simultaneous inference procedures
are usually undesirable in subgroup analysis, because they often yield false
negative discoveries and have inadequate power to confirm the most vulnerable
subgroup (Magnusson and Turnbull, 2013; Burke et al., 2015). In our context,
because subgroup analyses need to be conducted in observational studies (Yang
et al., 2020; Lu et al., 2018), the problem becomes even more challenging as
subgroup treatment effects need to take possibly high-dimensional confounders
into account. To address the above-mentioned issues, we provide new
inferential tools to help assess the efficacy of the post-hoc identified
subgroups from observational studies without having to resort to simultaneous
inference methods, which are often too conservative to start with.

By applying the proposed method to our study cohort, we find that, although
the overall adverse effect of statin use on developing T2D is not significant,
the female subgroup with more than two T2D risk alleles is identified as the
most vulnerable subgroup from the data. The statistical significance of such a
finding is also confirmed by the proposed bootstrap calibration method. In
sum, our case study not only provides new evidence supporting the
heterogeneous adverse effect of statin use from a biological perspective, but
also suggests that more caution should be taken when statins are prescribed,
especially for females who are already at higher risk of developing T2D. The
specific actions may include preventive treatments for diabetes and
recommendations on lifestyle changes.

The remainder of the manuscript is organized as follows. Section 2 describes
our study design, data source, and model setup. We then discuss our proposed
methodology in Section 3, followed by its theoretical investigations and
empirical studies in Section 4. In Section 5, we discuss our data analysis in
depth and verify our causal evidence with sensitivity analyses. Finally, we
summarize our findings and discuss the potential impacts of this case study in
Section 6.

##  2 Data description and model setup

###  2.1 Study design

We adopt the genetic variant (rs12916-T) as the surrogate treatment variable
of statin use. When the treatment indicator variable $\mbox{t}=1$, this means
that “the subject carries the variant rs12916-T.” When $\mbox{t}=0$, this
means that “the subject does not carry the variant rs12916-T.” Carrying
rs12916-T is a proxy for pharmacological action of statin use, because
rs12916-T, which resides in the HMGCR gene encoding the drug target of
statins, has been recently adopted as an unbiased, unconfounded proxy for
pharmacological action on the target of statins (i.e., HMG-CoA reductase
inhibition) (Swerdlow et al., 2015; Consortium et al., 2012). In particular,
Würtz et al. (2016) shows that the metabolic changes (e.g. decreased LDL
cholesterol level) associated with statin use have been found to resemble the
association between rs12916-T and the metabolic changes with $R^{2}=0.94$.
Therefore, we believe that rs12916-T is a credible surrogate treatment
variable of statin use.

Furthermore, following the causal diagram of our study design in Figure 1, we
argue that adopting rs12916-T as the surrogate treatment variable of statin
use invests our framework with two key benefits.

First, because our study design ensures the temporal precedence between
inheriting the variant rs12916-T (surrogate of statin use) and T2D onset, this
lifts the concern of reverse causation in conducting causal inference from
observational data. In particular, when studying the causal effect of statin
use on T2D risk from observational data, one normally assumes that statin use
causes the change in T2D risk (Pan et al., 2020; Swerdlow, 2016). This
suggests that only data collected from subjects who have recorded statin use
status before T2D onsets can be used for credible casual analyses.
Unfortunately, in our EHR data, it is impossible to establish the temporal
precedence between statin use and T2D onset. Using a genetic variant
(rs12916-T) as the treatment variable circumvents the above-mentioned issue.
Because genetic variants are randomly inherited at conception, our study
design guarantees that the cause (carrying rs12916-T [proxy for
pharmacological action of statin use] or not) must occur before T2D onset.
This clear temporal precedence makes the established causal relationships more
plausible (see Figure 1).

Figure 1:  The causal diagram under our study design.

Second, adopting rs12916-T as a surrogate treatment variable attenuates the
concern of unmeasured confounding biases. Given rs12916-T is naturally
inherited, variables observed after birth can at most be mediators that belong
to a different causal pathway. Since we aim to estimate the direct causal
effect between the treatment and the outcome, conditional on such (unmeasured)
mediators is not necessary. Thus, unmeasured confounding bias issues are
alleviated under our study design. Yet we still include additional potential
confounders that are obtained before birth, such as genetic variants
associated with T2D and potentially associated with rs12916-T, in order to
further robustify our causal conclusion and improve statistical estimation
efficiency of the underlying causal effect.

Note that while rs12916-T is sometimes used as an instrumental variable in
Mendelian Randomization analyses (Würtz et al., 2016), we use rs12916-T as a
surrogate treatment variable, instead of an instrumental variable. Our study
design is thus different from Mendelian randomization (MR) (Kang et al., 2016;
Windmeijer et al., 2019). Furthermore, because MR often assumes that the
causal direction between two traits is known as a priori (Xue and Pan, 2020)
and our data do not provide the temporal order between stain use and T2D
onset, the causal direction for MR analyses can not be specified and MR is not
suitable in our case study.

The outcome of interest is T2D status. We defer the technical details on
defining T2D status from EHR data to Section 2.3. Considering that the
existing literature (Mora et al., 2010; Waters et al., 2013) has suggested
that the treatment effect of statin use on T2D risk could be heterogeneous, we
conduct subgroup analyses to investigate the heterogeneous treatment effect.
Inspired by the study designs in Mora et al. (2010) and Waters et al. (2013)
where the authors divide patient population based on gender in the former and
the number of T2D baseline risk factors in the latter, we divide our PHS study
cohorts into six subgroups based on gender and baseline T2D genetic risk
profiles to conduct subgroup analyses. The T2D genetic risk is measured by the
number of copies of T2D risk allele of variants rs35011184-A and rs1800961-T
each subject has. The more T2D risk alleles the subject has, the higher
baseline T2D risk the subject bears (Lango et al., 2008). We define “low-risk”
as the total number of alleles $=0$, “mid-risk” as the total number of alleles
$=1$, and “high-risk” as the total number of alleles $\geq 2$. The six
subgroups are thereby divided as (1) high-risk female; (2) mid-risk female;
(3) low-risk female; (4) high-risk male; (5) mid-risk male, and (6) low-risk
male.

Here, we consider pre-defined subgroups instead of post-hoc identified
subgroups. Post-hoc identified subgroups are often adopted when there is no
prior information about the segregation of study population (Lipkovich et al.,
2011; Ma and Huang, 2017). In our setting, because previous studies (Mora et
al., 2010; Waters et al., 2013) suggest that T2D risk might be heterogeneous
across gender and baseline T2D genetic profiles, predefined subgroups
considered above are more suitable for the present case study.

###  2.2 Model setup

We work with the following sparse logistic regression model:

$\text{logit}\left\\{{\mathbb{P}}(\mbox{y}^{\intercal}=1\mid\mbox{z},\mbox{x})\right\\}=\mbox{z}^{\intercal}\beta+\mbox{x}^{\intercal}\gamma,\quad\left\|\gamma\right\|_{0}\ll
p.$ (1)

Here, y is the observed binary outcome representing the T2D status.
$\mbox{z}\in\mathbb{R}^{p_{1}\times n}$ includes variables representing
interactions between the treatment variable and all the six subgroup indicator
variables. $\mbox{x}\in\mathbb{R}^{p_{2}\times n}$ contains $336$ covariates
and an intercept (hence $p_{2}=337)$. The $336$ covariates contain $5$
subgroup indicator variables (the sixth subgroup, low-risk male, indicator
variable is dropped to avoid collinearity) and 331 potential confounders
(including race and age as baseline characteristics, and $329$ SNPs associated
with T2D related factors accounting for potential confounding issues). Note
that we do not include the treatment variable as a covariate because including
it causes collinearity issues. All observed covariates are obtained from
Partner Health System biobank.

Following the above setup, $\beta\in\mathbb{R}^{p_{1}}$ represents subgroup
treatment effects on the scale of log odds ratio (OR) (hence $p_{1}=6$). More
concretely, under Model (1),
$\beta=\big{(}\log\alpha_{1},\ldots,\log\alpha_{6}\big{)}$ with
$\log\alpha_{j}$ representing the log odds ratio of subgroup $j$, for
$j=1,\ldots,6$. Following the Neyman-Rubin causal model and our current study
design, we provide rigorous causal identification results to justify why the
model parametrization in (1) enables us to estimate the heterogeneous
treatment effects in the pre-defined subgroups. This theoretical justification
is provided in Supplementary Materials Section E. We further assume that
$\gamma\in\mathbb{R}^{p_{2}}$ is a sparse vector with the support set
$\texttt{M}_{0}$. The sparsity assumption not only provides a parsimonious
explanation of the data but also carries our prior belief that not every
genetic variant is predictive of the outcome.

###  2.3 Data description and exploration

Following our study design described in the previous section, we have
extracted information of $n=17,023$ subjects each with $p=337$ features
through linked diagnostics from the PHS EHR data and genotype information from
the PHS biobank data; see Table 1 for data summary.

Variable | Frequency (percent)
---|---
Sex |
Female | 7,592 (45)
Male | 9,431 (55)
Age (years) |
$<$40 | 2,333 (14)
40–50 | 1,733 (10)
50–60 | 2,824 (17)
60–70 | 3,971 (23)
70–80 | 4,042 (24)
$\geq$80 | 2,120 (12)
Race |
European | 15,048 (88)
African American | 1,004 (6)
Other/Unknown | 971 (6)
Ethnicity |
Hispanic or Latino | 698 (4)
Other/Unknown | 16,325 (96)
Number of rs12916-T allele |
2 | 6,245 (37)
1 | 8,079 (47)
0 | 2,699 (16)
With T2D |
Yes | 2,565 (15)
No | 14,458 (85)
Table 1: Demographics of 17,023 PHS subjects considered in our study.

Recall that the covariates contain age, race, subgroup indicators, and genetic
variants associated with T2D related factors (including LDL, high density
lipoprotein and obesity). As for the definition of the outcome, since the
diagnostic billing code for T2D has limited specificity in classifying the
true T2D status, we define the T2D status based on a previously validated
multimodal automated phenotyping (MAP) algorithm(Liao et al., 2019). The area
under the ROC curve (AUC) of MAP’s risk prediction score for classifying true
T2D status is $0.99$, and the specificity and sensitivity of its classifier
are $0.97$ and $0.92$ respectively. These suggest that the MAP classifier of
T2D can be reliably used to define the T2D outcome. Among our study cohort,
MAP classifies $2,565$ subjects having T2D.

To explore the association between statin use and T2D risk, we report
preliminary data exploration results from a “full” logistic regression model
for y against the treatment t and the covariates x in Table 2. There, although
a modest association is found between carrying the rs12916-T variant and T2D
status in the overall study cohort, unlike the results in Swerdlow et al.
(2015), this association is not statistically significant. Moreover, our
analysis reports only $16$ regression coefficients having $p$-values $<0.05$,
suggesting that the logistic regression coefficient vector is likely to be
sparse. Because the overall treatment effect is marginal (estimated marginal
treatment effect equals 0.04 with p-value 0.35), this motivates us to conduct
subgroup analyses to further investigate the heterogeneous adverse effect of
statin use on T2D risk.

Treatment effect | Est | SE |  $Z$-value |  $p$-value | # of significant coefficients
---|---|---|---|---|---
Full logistic regression | $0.04$ | $0.04$ | $0.93$ | $0.35$ | 16
Table 2: The estimated treatment effect (Est), standard error (SE), $Z$\- and
(two-sided) $p$-values of statins’ usage on the overall PHS study cohort. We
fit the “full” logistic regression model with $337$ features (age, race and
genetic information). “significant coefficients” are the estimated regression
coefficients with $p$-values $<0.05$.

###  2.4 Challenges in statistical inference

Because our goal is to uncover the patient subgroup most vulnerable for
developing T2D and assess the statistical replicability of our findings, our
methodological development hence centers around delivering valid inference
(accurate point estimate and valid confidence interval) on the effect size of
maximal regression coefficient
$\beta_{\max}=\underset{j\in\\{1,\ldots,p_{1}\\}}{\max}\beta_{j}$ in Model
(1). Making inference on this parameter allows us to inform public health
researchers and regulatory agencies whether extra caution should be taken for
the most vulnerable patient subgroup when prescribing statins.

In the presence of high-dimensional covariates, finding an accurate point
estimate and conducting inference on $\beta_{\max}$ can be a challenging task,
due to the presence of the regularization and selection biases. Regularization
bias occurs whenever penalization approaches are adopted to select a smaller
working model to enhance the estimation efficiency of $\beta$ (Wang et al.,
2019; Hong et al., 2018). Selection bias occurs whenever we use a simple
sample-analogue
$\widehat{\beta}_{\max}=\underset{j\in\\{1,\ldots,p_{1}\\}}{\max}\widehat{\beta}_{j}$
to estimate the true maximum effect $\beta_{\max}$, because
$\widehat{\beta}_{\max}$ tends to overestimate $\beta_{\max}$. Such an
overestimation phenomenon is well-recognized in post-hoc subgroup analyses
(see Zöllner and Pritchard, 2007; Cook et al., 2014, for example). While
several approaches have been proposed to address the regularization bias issue
(Zhang and Zhang, 2014; Li, 2020), and provide valid inference on a single
regression coefficient, these methods cannot account for the selection bias.
As for the selection bias, existing methods are mostly made for low-
dimensional data and are not directly applicable to observational data with
high dimensional covariates in this case study (Bornkamp et al., 2017; Guo and
He, 2020). While Guo et al. (2021) proposes bootstrap-based approaches to
simultaneously address the regularization and selection bias issues in high
dimensional linear models, we broaden its validity by providing a bootstrap
procedure that is asymptotically valid for high dimensional logistic
regression estimators with rigorous statistical guarantees. To the best of our
knowledge, debiasing procedures that simultaneously remove the regularization
bias and selection bias in high dimensional non-linear models have been
lacking. Technical discussions on these bias issues are deferred to the
Supplementary Materials Section A, and we demonstrate the selection and
regularization biases in estimating $\beta_{\max}$ within a simple simulated
example.

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

We use two widely adopted procedures to estimate $\beta$: (1) Lasso for
generalized linear models (GLM) (Park and Hastie, 2007), which estimates
$\beta$ with
$\big{(}\widehat{\beta}_{\texttt{GLasso}}^{\intercal},\widehat{\gamma}_{\texttt{GLasso}}^{\intercal}\big{)}^{\intercal}$
obtained from the $\ell_{1}$-penalized logistic regression program without any
adjustments, and (2) Refitted GLM Lasso, which estimates $\beta$ by refitting
the logistic regression model based on the covariates in the support set of
$\big{(}\widehat{\beta}_{\texttt{GLasso}}^{\intercal},\widehat{\gamma}_{\texttt{GLasso}}^{\intercal}\big{)}^{\intercal}$.
As a benchmark, we also report the performance of the oracle estimator
$\big{(}\widehat{\beta}_{\texttt{Oracle}}^{\intercal},\widehat{\gamma}_{\texttt{Oracle}}^{\intercal}\big{)}^{\intercal}$
which pretends the true support set of $\gamma$ is known and is estimated by
refitting the logistic regression model with the true support set.
$\beta_{\max}$ is then estimated in a two-step procedure: One first obtains an
estimate $\widehat{\beta}$ and then estimates $\beta_{\max}$ by taking the
maximum, that is
$\max\\{\widehat{\beta}_{1},\ldots,\widehat{\beta}_{p_{1}}\\}$. To mimic the
causal relationship in this case study, we generate Monte Carlo samples with
$\mbox{t}_{i}\sim\text{Bernoulli}(0.5)$ independent of the covariate
$\mbox{w}_{i}\sim N(0,\Sigma)$, where $\Sigma=(\Sigma_{jk})_{j,k=1}^{p-6}$ and
$\Sigma_{jk}=0.5^{|j-k|}$ for $i=1,\ldots,n$. We then generate
$\mbox{x}_{ij}=\mathds{1}(\mbox{w}_{ij}>0)$ for $1\leq j\leq p-6$, and
$\mbox{z}_{il}=t_{i}x_{il}$, $l=1,\ldots,6$. $\mbox{y}_{i}$ is generated
following Model (1). We set the sample size $n=1,000$ and the dimension
$p=200$ and set the coefficients $\beta=(0.5,0.5,0,0,0,0)^{\intercal}$ and
$\gamma=(1,1,0,\ldots)\in\mathbb{R}^{p-6}$. In Figure 2, we report the
root-$n$ scaled bias based on $500$ Monte Carlo samples.

From the results in Figure 2, we observe that all three estimators are biased.
Although $\widehat{\beta}_{\text{Oracle}}$ is a consistent estimator of
$\beta$, its maximum is not centered around $\beta_{\max}$. The residual bias
in the oracle estimator $\widehat{\beta}_{\text{Oracle}}$ is caused by the
selection bias issue. In fact, following some explicit evidence given in
Nadarajah and Kotz (2008), $\widehat{\beta}_{\max}$ is usually biased upward
for estimating $\beta_{\max}$. The magnitude of this bias crucially depends on
the unknown parameters $\beta$ and the covariance structure of
$\widehat{\beta}_{\text{Oracle}}$. Hence, any inference procedure without
adjusting for selection bias cannot be valid. On the top of the selection bias
issues, the GLM Lasso and the refitted estimators suffer from the
regularization bias and hence are also not correctly centered around
$\beta_{\max}$, unless in some special cases where the regularization bias and
the selection bias cancel out.

Figure 2:  Root-$n$ scaled bias for Example 1. The maximum Monte Carlo
standard error of the above results is $0.31$.

To simultaneously adjust for the selection bias and the regularization bias
without knowing the underlying true parameters, in the following section, we
propose an inferential framework that produces a bias-reduced estimate as well
as a valid confidence interval of $\beta_{\max}$.

##  3 Methodology

###  3.1 Bootstrap calibrated debiased estimates

We start with describing an estimation strategy of $\beta$ that resolves the
regularization bias induced by model selection and helps addressing our
research objectives (ii) and (iii) discussed in Section 1.1. Regularization
bias arises when the selected model is either over-fitted or under-fitted; see
detailed discussion provided in the Supplementary Materials Section A. While
the risk of under-fitting can be mitigated by aiming for a larger model for
parameter estimation, we resolve the issue of over-fitting by sample
splitting. Sample splitting divides a sample into two parts: The first part of
the sample is used for model selection and the remaining part is used for
estimation based on the selected model. When $\gamma$ is sparse and a larger
model is selected on the first half of the sample, we expect refitted GLM
estimator on the second part of the sample to be free of significant bias.
Nevertheless, sample splitting provides debiased estimator of $\beta$ at a
cost of increased variability, because only a part of the sample is used for
estimation. To minimize this efficiency loss due to sample splitting, we
consider the method of repeated sample splitting (R-Split) that averages
different estimates of $\beta$ across different splits. Our strategy, in a
spirit similar to bagging and ensemble algorithms in machine learning, helps
to stabilize and improve the accuracy of the estimated $\beta$ in a subsample.

Step 1

(Repeated sample splitting that accounts for the regularization bias) For
$b\leftarrow 1$ to $B_{1}$: (1) Randomly split the sample
$\\{(\mbox{y}_{i},\mbox{x}_{i},\mbox{z}_{i})\\}_{i=1}^{n}$ into two
subsamples: a subsample $T_{1}$ of size $n_{1}$ and a subsample $T_{2}$ of
size $n_{2}=n-n_{1}$; (2) select a model $\widehat{\texttt{M}}_{b}$ to predict
y based on $T_{1}$; (3) refit the selected model with the data in $T_{2}$ to
estimate $\beta_{b}$ and $\gamma_{b}$ via logistic regression:

$\displaystyle(\widehat{\beta}_{b}^{\intercal},\widehat{\gamma}_{b}^{\intercal})^{\intercal}=\arg\min\left\\{\sum_{l\in
T_{2}}\Big{(}\mbox{y}_{l}\cdot(\mbox{z}_{l}^{\intercal}\beta+\mbox{x}_{l,\widehat{\texttt{M}}_{b}}^{\intercal}\gamma)-\log\big{(}1+\exp(\mbox{z}_{l}^{\intercal}\beta+\mbox{x}_{l,\widehat{\texttt{M}}_{b}}^{\intercal}\gamma)\big{)}\Big{)}\right\\}.$

(4) obtain the R-Split estimate:
$\widetilde{\beta}=\frac{1}{B_{1}}\sum_{b=1}^{B_{1}}\widehat{\beta}_{b}$.

In this step, any reasonable model selection procedures may be used and the
choice of model size is subjective, but the selected model needs to be large
enough for the under-fitting bias to be negligible. In our simulation and real
data analyses, we use GLM Lasso for model selection, and choose the model size
from cross-validation (see Section 3.2 for detailed description). The choice
of splits $B_{1}$ needs to be sufficiently large so that the R-Split estimator
$\widetilde{\beta}$ has a tractable asymptotic distribution. Under appropriate
regularity conditions, we show that $\widetilde{\beta}$ converges to a normal
distribution centered around $\beta$ at a root-$n$ rate (statistical
justification is provided in the Supplementary Materials Section C.2).

As $\widetilde{\beta}$ provides an accurate estimate of $\beta$, we use
$\widetilde{\beta}$ to address our research objectives (ii) and (iii). In
particular, the subgroup with the largest coefficient,
$\arg\max_{j\in[p_{1}]}\widetilde{\beta}_{j}$, is most vulnerable for
developing T2D after taking statins. However, due to the selection bias,
simply relying on $\widetilde{\beta}$ will not lead to valid inference on
$\beta_{\max}$, and we need a second step to address objective (iii). Built
upon an accurate estimate of $\beta$, we store an inverse Hessian matrix for
the later bootstrap calibration to adjust for selection bias:

$\widetilde{\Gamma}_{n}=\frac{1}{B_{1}}\sum_{b=1}^{B_{1}}\mbox{I}_{z}\left(\frac{1}{n_{1}}\sum_{i\in
T_{2,b}}f_{i,b}\begin{pmatrix}\mbox{z}_{i}\\\
\mbox{x}_{i,\widehat{\texttt{M}}_{b}}\end{pmatrix}(\mbox{z}_{i}^{\intercal},\mbox{x}^{\intercal}_{i,\widehat{\texttt{M}}_{b}})^{\intercal}\right)^{-1}\mbox{I}_{\widehat{\texttt{M}}_{b}},$

where
$f_{i,b}=\text{expit}^{\prime}(\mbox{z}_{i}^{\intercal}\widehat{\beta}_{b}+\mbox{x}^{\intercal}_{i,\widehat{\texttt{M}}_{b}}\widehat{\gamma}_{b})$.
Its benefits will be apparent in the following step:

Step 2

(Calibrated bootstrap that accounts for the selection bias) For $b\leftarrow
1$ to $B_{2}$: generate bootstrap replicate $\widetilde{\beta}^{*}$ from:

$\displaystyle\widetilde{\beta}^{*}=\widetilde{\beta}+\widetilde{\Gamma}_{n}\cdot\frac{1}{n}\sum_{i=1}^{n}\begin{pmatrix}\mbox{z}_{i}\\\
\mbox{x}_{i}\end{pmatrix}\nu^{*}_{i},$ (2)

where $\nu_{i}^{*}=u_{i}\widehat{\nu}_{i}$ is the permuted GLasso residual,
$\widehat{\nu}_{i}=\mbox{y}_{i}-\text{expit}(\mbox{z}_{i}^{\intercal}\widehat{\beta}_{\texttt{GLasso}}+\mbox{x}_{i}^{\intercal}\widehat{\gamma}_{\texttt{GLasso}})$.
Then recalibrate bootstrap statistics via

$\displaystyle
T_{b}^{*}=\underset{j\in[p_{1}]}{\max}(\widetilde{\beta}^{*}_{j}+\widetilde{c}_{j}(r))-\widetilde{\beta}_{\max},\quad\widetilde{c}_{j}(r)=(1-n^{r-0.5})(\widetilde{\beta}_{\max}-\widetilde{\beta}_{j}),\text{
where }r\in(0,0.5).$

In this step, rather than adopting the simple bootstrap statistics
$\max_{j\in[p_{1}]}\widetilde{\beta}^{*}_{j}-\widetilde{\beta}_{\max}$ to make
inference on $\beta_{\max}$, we make an adjustment to each coordinate of
$\widetilde{\beta}^{*}$ by the amount $\widetilde{c}_{j}(r)$. This is because
just as $\widetilde{\beta}_{\max}$ is a biased estimator of $\beta_{\max}$,
the simple bootstrap statistics $\max_{j\in[p_{1}]}\widetilde{\beta}^{*}_{j}$
is also not centered at $\widetilde{\beta}_{\max}$. The amount of adjustment
$\widetilde{c}_{j}(r)$ is large when $\widetilde{\beta}_{j}$ is small, and is
small when $\widetilde{\beta}_{j}$ is large. By adding the correction term
$\widetilde{c}_{j}(r)$, under certain regularity conditions, the distributions
of
$\sqrt{n}(\widetilde{\beta}^{*}_{\text{modified;max}}-\widetilde{\beta}_{\max})$
and $\sqrt{n}(\widetilde{\beta}_{\max}-\beta_{\max})$ are asymptotically
equivalent, implying that our proposed method adjusts for the selection bias
and the regularization bias simultaneously, where
$\widetilde{\beta}^{*}_{\mathrm{modified};\max}=\underset{j\in[p_{1}]}{\max}(\widetilde{\beta}^{*}_{j}+\widetilde{c}_{j}(r))$.
We relegate the theoretical details of this bootstrap calibration procedure in
the Supplementary Materials Section C. Note that $r\in(0,0.5)$ is a positive
tuning parameter (see Section 3.2 for its data adaptive choice).

At this point, we note that our procedure adopts wild bootstrap to construct
bootstrapped statistics of the R-Split estimate $\widetilde{\beta}$. The wild
bootstrap procedure adopted here is not only computationally efficient in high
dimensions, as the Hessian matrix remains unchanged across different bootstrap
samples, but also provably consistent in our problem setup. Furthermore,
Dezeure et al. (2017) shows that wild bootstrap is more versatile than other
residual bootstrap methods because it correctly captures the asymptotic
variance in the presence of homoscedastic or heteroscedastic errors. With the
help of a valid bootstrap calibration procedure in replicating
$\widetilde{\beta}_{\max}$, we are now ready to propose our final step that
constructs confidence intervals and debiased estimate for $\beta_{\max}$:

Step 3

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

###  3.2 Practical implementation

In Step 1 (1), we recommend a split ratio of $0.6:0.4$ for $n_{1}:n_{2}$
because a larger sample size for subsample $T_{1}$ improves model selection
accuracy. In Step 1 (2), the model selection procedure can be any easily
accessible procedure. In our case, since the real data have binary outcomes,
we adopt GLM lasso for model selection with R package glmnet (Park and Hastie,
2007). The model size is selected via cross-validation with a constraint on
the maximal and minimal model sizes. We recommend to set $B_{1}=500$ for the
number of repeated splits and $B_{2}=1,000$ for the number of bootstrap
replications.

As for the tuning parameter $r$, we propose a data-adaptive cross-validated
algorithm to select $r$ as the following (Guo and He, 2020):

Step 1

Denote $R=\\{r_{1},\dots,r_{m}\\}$ as a set of candidate tuning parameters.
Randomly split the sample into $v$ equal-sized subsamples.

Step 2

For l $\leftarrow$ 1 to m:

For j $\leftarrow$ 1 to $v$:

(a)

Use subsample $j$ as reference data and the rest as training data. Obtain
$\widetilde{\beta}_{\max,\mathrm{reduced},j}(r_{l})$ on the training data,
where $r_{l}$ is the tuning parameter.

(b)

For i $\leftarrow$ 1 to k: Obtain R-Split estimate of
$\widetilde{\beta}_{i,j}$ and its standard error  $\widetilde{\sigma}_{i,j}$
on the reference data; evaluate the accuracy

$\displaystyle
h_{i,j}(r_{l})=(\widetilde{\beta}_{\max,\mathrm{reduced},j}(r_{l})-\widetilde{\beta}_{i,j})^{2}-\widetilde{\sigma}_{i,j}^{2}.$

Step 3

Select the tuning parameter via
$\arg\min_{r_{l}}\\{\min_{i\in[k]}[\sum_{j=1}^{j=v}h_{i,j}(r_{l})/v]\\}.$

Intuitively, we would like to choose $r$ that minimizes the mean squared error
between the proposed bias reduced estimate
$\widetilde{\beta}_{\max,\mathrm{reduced}}$ and $\beta_{\max}$. To make this
operational without knowing $\beta_{\max}$, we provide an approximation of the
mean squared error that can be computed from the data via cross-validation in
Step 2 (b). The justification of this cross validation method for fixed
$p_{1}$ can be found in Guo and He (2020). In our empirical work, we implement
the above tuning selection method via three-fold cross-validation with a
candidate set $R=\\{1/3,1/6,\dots,1/30\\}$.

##  4 Theoretical and empirical justification

In this section, we provide theoretical justifications of the proposed
bootstrap-assisted R-Split estimator along with a simple power analysis, where
we demonstrate that our approach not only has rigorous theoretical guarantee
but also shows high statistical detection power. We then examine the
performance of the proposed method through simulation studies.

###  4.1 Theoretical investigation and a power analysis

The following theorem confirms that the asymptotic distribution of
$\sqrt{n}(\widetilde{\beta}_{\text{modified;max}}-\widetilde{\beta}_{\max})$
converges to $\sqrt{n}(\widetilde{\beta}_{\max}-\beta_{\max})$. This suggests
that the proposed confidence interval constructed in Step 3 of Section 3.1 is
“asymptotically sharp,” meaning that it achieves the exact nominal level as
the sample size goes to infinity. This distinguishes the proposed procedure
from other conservative methods made for subgroup analysis. The proof of
Theorem 1 is provided in the Supplementary Materials Section C.3. To simplify
presentation, we relegate regularity assumptions to the Supplementary
Materials Section C.1.

######  Theorem 1.

Under Assumptions 1-9 given in the Supplementary Materials Section C.1, when
$p_{1}$ is a fixed number, the modified bootstrap maximum treatment effect
estimator,
$\widetilde{\beta}^{*}_{\mathrm{modified};\max}=\underset{j\in[p_{1}]}{\max}(\widetilde{\beta}^{*}_{j}+\widetilde{c}_{j}(r))$,
satisfies:

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

Figure 3:  Power comparison for bootstrap-assisted R-Split, bootstrap-assisted
logistic regression, R-Split with simultaneous confidence interval, and the
desparsified Lasso with simultaneous confidence interval evaluated over 500
Monte Carlo samples.

The above theoretical result has two direct implications. On the one hand, as
the proposed bootstrap calibration strategy successfully replicates the
distribution of $\sqrt{n}(\widetilde{\beta}_{\max}-\beta_{\max})$, our bias-
reduced estimator discussed in the Step 3 of Section 3.1 simultaneously
removes the regularization bias and the selection bias in
$\widetilde{\beta}_{\max}$. On the other hand, although simultaneous inference
also delivers valid inference on $\beta_{\max}$ with strict Type-I error rate
control, our proposal delivers valid inference on $\beta_{\max}$ without
sacrificing the statistical power. This property is more desirable in our
problem setup as we aim to look for the subgroup with the most severe side
effect of statin use while simultaneous methods often lead to overly
conservative conclusions for this purpose.

To further demonstrate the merit of constructing an asymptotically sharp
confidence interval for $\beta_{\max}$ and the benefit of conducting variable
selection in finite samples, we compare statistical power for testing the null
hypothesis $H_{0}:\beta_{\max}=0$ for four procedures: (1) the proposed
bootstrap-assisted R-Split, (2) R-Split with simultaneous confidence
intervals, (3) the proposed bootstrap-assisted logistic regression, and (4)
the desparsified Lasso estimator discussed in Zhang and Zhang (2014) with
simultaneous confidence interval (Dezeure et al., 2017; Fuentes et al., 2018).
We follow the same simulation setup as in Example 1. The tuning parameter is
fixed at $r=0.15$ for simplicity. For R-Split, we choose the model size via
cross-validation (see Section 3.2) with a minimal model size equals 3 and a
maximal model size equals 10.

From Figure 3, we observe that all considered approaches control the Type-I
error rate at the nominal level when $\beta_{\max}=0$. The bootstrap-assisted
R-Split has the highest detection power over a range of $\beta_{\max}$ among
all considered procedures. The bootstrap-assisted logistic regression has the
lowest detection power, which demonstrates the necessity of conducting
variable selection to screen out irrelevant predictors. As we have expected,
both the R-Split method with simultaneous confidence interval and the
desparsified Lasso with simultaneous confidence interval do not retain
sufficient statistical power to detect subgroup treatment effect
heterogeneity.

###  4.2 Simulation studies

In this section, we consider various simulation designs to demonstrate the
merit of our proposal. There are three main takeaways from this simulation.
First, our proposed bootstrap calibration procedure provides confidence
intervals with nominal coverage probabilities of $\beta_{\max}$ in finite
samples. Second, R-Split based methods provide more accurate point estimates
and shorter confidence intervals than the logistic regression based approaches
without variable selection. Third, the bootstrap-assisted methods have higher
statistical efficiency (shorter confidence intervals) compared to the
simultaneous methods.

We generate Monte Carlo samples from the following model:

$\displaystyle\text{logit}\left\\{\mathbb{P}(\mbox{y}_{i}=1\mid\mbox{z}_{i},\mbox{x}_{i})\right\\}=\mbox{z}_{i}^{\intercal}\beta+\mbox{x}_{i}^{\intercal}\gamma,\quad
i=1,\ldots,n,$

with $n=2,000$. We consider two cases for $\beta$: (1) heterogeneous case with
$\beta=(0,\ldots,0,1)^{\intercal}\in\mathbb{R}^{p_{1}}$, meaning that there
exists subgroup treatment effect heterogeneity and only one subgroup singles
out; and (2) spurious heterogeneous case with
$\beta=(0,\ldots,0,0)\in\mathbb{R}^{p_{1}}$, meaning that there is no subgroup
with significant treatment effect in the population. We set
$\gamma=(1,1,1,1,0,\dots,0)\in\mathds{R}^{p_{2}}$. In all considered
simulation designs, we set $p_{1}\in\\{4,10\\}$. We consider the case with
$(n,p_{2})=(2,000,~{}150)$ for logistic regression, R-Split, and the
desparsified Lasso (Zhang and Zhang, 2014), and consider the case with
$(n,p_{2})=(2,000,~{}500)$ for R-Split and the desparsified Lasso, since
logistic regression tends to provide inconsistent estimates in moderately high
dimensions (Sur and Candès, 2019).

As for the covariates design, we generate $\mbox{z}_{i}$ and $\mbox{x}_{i}$
from

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

where $\mbox{x}_{i}\sim N(0,\Sigma)$ with $\Sigma_{ij}=0.5^{|i-j|}$. We
compare the finite sample performance of the proposed bootstrap-assisted
R-Split and the bootstrap-assisted logistic regression with two benchmark
methods: (1) a naive method with no bootstrap calibration, which directly uses
the estimated maximum coefficient to estimate $\beta_{\max}$ and (2) the
simultaneous method as discussed in Dezeure et al. (2017) and Fuentes et al.
(2018). For the desparsified Lasso (Zhang and Zhang, 2014), we only consider
the above-mentioned two benchmark methods: the naive method and the
simultaneous method (without bootstrap calibration). For the R-Split method,
we choose the model size via cross-validated GLM Lasso with a minimal model
size equals 3. We report the coverage probability, the $\sqrt{n}$ scaled
confidence interval length and the $\sqrt{n}$ scaled Monte Carlo bias along
with their standard errors based on 1,000 Monte Carlo samples in Table 3 and
Table 4.

Comparing the bootstrap-assisted methods with the naive methods, we observe
that the bootstrap-assisted methods attain nominal-level coverage
probabilities, while the naive methods are biased and under-covered. This
comparison verifies the theoretical results in Section 4.1 that the proposed
bootstrap calibration successfully adjusts for the selection bias.

Comparing the bootstrap-assisted methods with the simultaneous methods, we
find that although simultaneous methods have higher coverage probabilities,
the confidence intervals are rather long, implying that simultaneous methods
are overly conservative. While our proposed inferential framework reaches the
nominal-level coverage probabilities and has shorter confidence intervals
leading to asymptotically sharp inference.

The comparison between the bootstrap-assisted R-Split with the bootstrap-
assisted logistic regression shows that the latter has larger biases and lower
coverage probabilities. The bootstrap-assisted logistic regression has
undesirable performance because logistic regression yields biased estimates in
moderately high dimensions (Sur and Candès, 2019). This comparison reveals the
benefit of conducting variable selection when $\gamma$ is sparse and the
dimension of covariates is large, and it confirms that R-Split alleviates the
regularization bias issue. Comparing R-Split with the desparsified Lasso, in
line with our earlier conjecture, we observe that the desparsified Lasso
approach has wider confidence intervals than those obtained by R-Split and
tends to provide conservative inference.

This simulation study verifies that our proposed inferential framework not
only achieves nominal coverage probabilities, but also mitigates the
regularization and selection biases. Thus, the proposed inferential framework
is sensible to consider for our case study.

Table 3:  Simulation results (heterogeneous case)

$\beta=(0,\dots,0,1)\in\mathbb{R}^{p_{1}}$ (heterogeneity) Logistic Regression
($\bm{p_{2}=150}$) Boot-Calibrated No adjustment Simultaneous $p_{1}=4$ Cover
0.95(0.01) 0.89(0.01) 0.99(0.01) $\sqrt{n}$Length 9.57(0.05) 8.83(0.03)
14.6(0.04) $\sqrt{n}$Bias -2.51(2.70) 4.91(4.60) — $p_{1}=10$ Cover 0.93(0.01)
0.86(0.01) 0.99(0.01) $\sqrt{n}$Length 10.4(0.04) 9.18(0.04) 16.7(0.03)
$\sqrt{n}$Bias -4.07(3.89) 5.17(4.67) — Repeated Sample Splitting
($\bm{p_{2}=150}$) Boot-Calibrated No adjustment Simultaneous $p_{1}=4$ Cover
0.96(0.01) 0.94(0.01) 0.99(0.00) $\sqrt{n}$ Length 3.56(0.07) 2.17(0.06)
5.14(0.07) $\sqrt{n}$Bias 0.11(0.26) 0.14(0.25) — $p_{1}=10$ Cover 0.95(0.02)
0.92(0.01) 0.99(0.01) $\sqrt{n}$Length 3.62(0.07) 2.57(0.05) 6.61(0.05)
$\sqrt{n}$Bias 0.25(0.39) 0.32(0.30) — Desparsified Lasso ($\bm{p_{2}=150}$)
Boot-Calibrated No adjustment Simultaneous $p_{1}=4$ Cover — 0.92(0.01)
0.99(0.00) $\sqrt{n}$Length — 2.13(0.06) 6.51(0.05) $\sqrt{n}$Bias —
0.29(0.22) — $p_{1}=10$ Cover — 0.93(0.01) 0.99(0.01) $\sqrt{n}$Length —
2.10(0.07) 6.98(0.07) $\sqrt{n}$Bias — 0.27(0.17) — Repeated Sample Splitting
($\bm{p_{2}=500}$) Boot-Calibrated No adjustment Simultaneous $p_{1}=4$ Cover
0.95(0.02) 0.92(0.03) 0.99(0.00) $\sqrt{n}$Length 4.44(0.06) 2.22(0.06)
6.08(0.05) $\sqrt{n}$Bias -0.68(0.80) 1.22(1.18) — $p_{1}=10$ Cover 0.93(0.02)
0.88(0.03) 0.98(0.01) $\sqrt{n}$Length 5.11(0.04) 2.95(0.05) 6.77(0.05)
$\sqrt{n}$Bias -0.90(0.85) 1.36(1.20) — Desparsified Lasso ($\bm{p_{2}=500}$)
Boot-Calibrated No adjustment Simultaneous $p_{1}=4$ Cover — 0.90(0.01)
0.99(0.00) $\sqrt{n}$Length — 2.19(0.05) 7.48(0.08) $\sqrt{n}$Bias —
1.29(1.13) — $p_{1}=10$ Cover — 0.91(0.01) 0.99(0.01) $\sqrt{n}$Length —
2.15(0.06) 7.60(0.08) $\sqrt{n}$Bias — 1.25(1.17) —

  * •

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

Table 4:  Simulation results (spurious heterogeneous case)

$\beta=(0,\dots,0,0)\in\mathbb{R}^{p_{1}}$ (spurious heterogeneity) Logistic
Regression ($\bm{p_{2}=150}$) Boot-Calibrated No adjustment Simultaneous
$p_{1}=4$ Cover 0.94(0.02) 0.87(0.03) 0.99(0.01) $\sqrt{n}$Length 7.90(0.05)
6.21(0.05) 11.1(0.04) $\sqrt{n}$Bias 3.10(3.44) 5.05(4.46) — $p_{1}=10$ Cover
0.91(0.01) 0.83(0.02) 0.98(0.01) $\sqrt{n}$Length 8.38(0.06) 7.25(0.05)
12.4(0.06) $\sqrt{n}$Bias 5.30(4.88) 7.32(6.58) — Repeated Sample Splitting
($\bm{p_{2}=150}$) Boot-Calibrated No adjustment Simultaneous $p_{1}=4$ Cover
0.95(0.02) 0.93(0.02) 0.98(0.02) $\sqrt{n}$ Length 1.87(0.04) 1.03(0.06)
5.08(0.04) $\sqrt{n}$Bias 0.15(0.24) 0.31(0.39) — $p_{1}=10$ Cover 0.95(0.01)
0.91(0.02) 0.96(0.01) $\sqrt{n}$Length 2.02(0.06) 1.47(0.06) 6.46(0.04)
$\sqrt{n}$Bias 0.29(0.40) 0.98(0.90) — Desparsified Lasso ($\bm{p_{2}=150}$)
Boot-Calibrated No adjustment Simultaneous $p_{1}=4$ Cover — 0.92(0.01)
0.99(0.01) $\sqrt{n}$ Length — 1.01(0.07) 5.52(0.05) $\sqrt{n}$Bias —
1.23(0.99) — $p_{1}=10$ Cover — 0.93(0.01) 0.99(0.01) $\sqrt{n}$Length —
1.39(0.06) 6.20(0.07) $\sqrt{n}$Bias — 0.97(0.85) — Repeated Sample Splitting
($\bm{p_{2}=500}$) Boot-Calibrated No adjustment Simultaneous $p_{1}=4$ Cover
0.95(0.02) 0.91(0.02) 0.98(0.01) $\sqrt{n}$Length 3.77(0.05) 3.18(0.04)
5.90(0.04) $\sqrt{n}$Bias 0.62(0.72) 1.58(1.41) — $p_{1}=10$ Cover 0.92(0.02)
0.85(0.01) 0.95(0.01) $\sqrt{n}$Length 3.54(0.06) 2.72(0.06) 6.52(0.05)
$\sqrt{n}$Bias 1.53(1.39) 2.82(1.97) — Desparsified Lasso ($\bm{p_{2}=500}$)
Boot-Calibrated No adjustment Simultaneous $p_{1}=4$ Cover — 0.89(0.01)
0.99(0.01) $\sqrt{n}$ Length — 3.10(0.05) 6.88(0.08) $\sqrt{n}$Bias —
2.30(1.90) — $p_{1}=10$ Cover — 0.90(0.01) 0.99(0.01) $\sqrt{n}$Length —
2.68(0.05) 7.63(0.08) $\sqrt{n}$Bias — 2.08(1.96) —

  * •

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

##  5 Case study

###  5.1 Case study results

In this section, we investigate the adverse effect of statin use in our pre-
specified six subgroups divided by gender and T2D genetic risk using the data
introduced in Section 2.3. We compare the results from three methods: (1)
repeated sample splitting (R-Split) without bootstrap calibration, (2) R-Split
based on the simultaneous method discussed in Dezeure et al. (2017), and (3)
the proposed bootstrap-assisted R-Split. We summarize our real data analyses
results in Table 5, in which we have reported the estimated subgroup treatment
effects from R-Split along with their $p$-values and two-sided confidence
intervals, adjusted $p$-values to account for the multiple comparisons issue
with simultaneous method and Bonferroni correction, and bootstrap calibrated
$p-$values for the subgroup with the largest treatment effect. The results
with one-sided confidence lower bounds are summarized in Supplementary
Materials Section F.

Method Subgroup (prevalence; # of case) Est (95% CI) $p$-value Bonf $p$-value
R-Split High-risk female $(0.14,100)$ $0.41~{}(0.04,~{}0.78)$ $0.030$ $0.180$
(without bootstrap calibration) Mid-risk female $(0.12,396)$
$0.10~{}(-0.03,~{}0.24)$ $0.132$ $0.792$ Low-risk female $(0.11,630)$
$-0.00~{}(-0.10,~{}0.09)$ $0.990$ $1$ High-risk male $(0.24,139)$
$-0.07~{}(-0.38,~{}0.25)$ $0.658$ $1$ Mid-risk male $(0.21,561)$
$0.02~{}(-0.07,~{}0.11)$ $0.673$ $1$ Low-risk male $(0.17,739)$
$-0.03~{}(-0.16,~{}0.10)$ $0.651$ $1$ Overall $0.07~{}(-0.16,~{}0.39)$ $0.545$
– Simultaneous High-risk female $(0.14,100)$ – $0.256$ – Bootstrap-assisted
R-Split High-risk female $(0.14,100)$ $0.35~{}(0.02,~{}0.70)$ $0.037$ –

Table 5:  Estimated treatment effects (Est) on the PHS cohort in six subgroups
divided by gender and T2D genetic risk, together with two-sided $95\%$
confidence intervals (CI), corresponding two-sided $p$-values and the
Bonferroni $p$-values in the last column. We also present the prevalence of
T2D in each subgroup.

From Table 5, the results of the R-Split estimator without bootstrap
calibration not only indicate that the treatment effect of statins tends to
vary across different subgroups, but also suggest that the high-genetic-risk
female subgroup is the most vulnerable group for developing T2D with estimated
log-odds ratio $0.41$, 95% two-sided confidence interval $0.04-0.78$ (${\rm
OR}=1.04-2.18$) with $p$-value $0.030$. For males with various genetic risk
levels and females with lower T2D genetic risk, the adverse effects of statin
use are not significant based on R-Split without bootstrap calibration. The
treatment effect in the overall study cohort is slightly positive but is not
significant, which is in-line with our expectation from the preliminary
analysis in Section 2.3.

Although the estimates and confidence intervals from the R-Split without
bootstrap calibration suggest that taking statins causes the increased risk of
developing T2D for the most vulnerable subgroup, the statistical significance
of this finding is unclear since R-Split is implemented without bootstrap
calibration and can not address the multiple comparisons issue. After
accounting for the multiple comparisons issue through conservative procedures
including the simultaneous method or Bonferroni correction, the $p$-values for
the female high-risk group are no longer significant, seemingly suggesting
that our data do not provide enough evidence to claim the existence of the
adverse effect of statin use in the female high-risk subgroup. This might be
due to the fact that both the simultaneous method and Bonferroni correction
are rather conservative and tend to provide false negative discoveries.
Fortunately, our proposed bootstrap assisted R-Split procedure directly
conducts inference on the most vulnerable group, and our results suggest that
among high-genetic-risk female patients, the odds of developing T2D after
taking statins are 1.42 times the odds of developing T2D for the patients
without taking statins ($p-$value $0.037$ for two-sided test).

Our findings are in-line with reported results in existing clinical studies.
For example, Mora et al. (2010) suggest that statin use incurs a larger T2D
risk increment on females than on males, and Waters et al. (2013) suggest that
statins only significantly increase the risk of T2D on those with at least
three out of four common T2D risk factors at baseline.222The risk factors used
by Waters et al. (2013) include high fasting blood glucose, history of
hypertension, high body mass index, and high fasting triglycerides. Compared
with the existing studies, our findings provide more robust evidence with the
new data analysis pipeline built under the causal inference framework. Our
data analysis pipeline addresses several limitations of existing studies; in
particular, limited sample size and multiple comparisons issue. Moreover,
compared to existing studies, our findings provide a more biologically driven
depiction of statins’ heterogeneous adverse effect, which can further support
effective and precise clinical decisions and actions concerning the
prescription of statins. Our study further demonstrates that in practice, the
genetic profiles could assist T2D prevention of statin receivers to improve
the quality of clinical practices.

###  5.2 Sensitivity analysis

To evaluate the validity of causal conclusions derived from our real data
analyses, we conduct sensitivity analyses with the E-value method. The E-value
method computes the minimal strength of an unmeasured confounder to explain
away the estimated causal effect (VanderWeele and Ding, 2017). Practitioners
could then evaluate if there exists such an unmeasured confounder with the
strength quantified by the E-value. A large E-value implies that the
unmeasured confounder needs to have a strong association with the outcome and
the treatment in order to explain away the causal evidence. The E-values for
our estimated subgroup causal effects are summarized in Table 6. Table 6 shows
that the E-value in the high-risk female group is 2.38, which implies that an
unmeasured confounder that is associated with both the treatment and the
outcome 2.38 times stronger than the measured confounders and could explain
away the estimated causal effect. In other words, this E-value indicates that
the association of the outcome and treatment with the unmeasured confounders
should be 2.38-fold larger than their association with the observed
confounders. According to a meta-study on E-value applications, most computed
E-values from existing literature are below $2.0$ (Ioannidis et al., 2019). In
sum, the results from Table 6 imply that the causal evidence collected from
our data is reasonably robust.

Method Subgroup (prevalence; # of case) E-value R-Split High-risk female
$(0.14,100)$ $2.38$ (without bootstrap calibration) Mid-risk female
$(0.12,396)$ $1.45$ Low-risk female $(0.11,630)$ $1.00$ High-risk male
$(0.24,139)$ $1.23$ Mid-risk male $(0.21,561)$ $1.11$ Low-risk male
$(0.17,739)$ $1.14$ Overall $1.23$ Bootstrap-assisted R-Split High-risk female
$(0.14,100)$ $2.19$

Table 6:  Sensitivity analysis of our causal evidence measured by the E-value.

##  6 Discussion

In this case study, we investigate the heterogeneous T2D risk associated with
statin use. To overcome the limitations of existing studies and generate more
trustworthy evidence, we introduce a rigorous study design under the causal
inference framework and based on the EHR and biobank data from the Partner
Health System. Built on this study design, we find that although the adverse
effect of statin use for developing T2D is only marginal for the overall study
cohort, taking statins significantly increases the risk of developing T2D for
female patients with high genetic predisposition to T2D. While the objective
of this case study is to make inference on the most vulnerable subgroup, a
natural question to ask is whether statin use will significantly increase the
T2D risk for the other subgroups. The point estimates of statins’ side effects
on T2D risk suggest that minimal risk is associated with male subjects and
female subjects with low genetic predisposition to T2D. However, rigorous
statistical inference for these subgroups adjusting for multiple comparisons
warrants future research.

This case study is made for subgroup analysis with an exploratory purpose, as
our analysis is built upon observational data. Although confirmatory clinical
trials are more desirable for confirming the heterogeneous T2D risk associated
with statin use, compared with existing studies, our findings can be helpful
for researchers to design further confirmatory clinical trials with a focus on
the high-genetic-risk female subpopulation.

This case study considers pre-defined candidate subgroups. While predefined
subgroups are more suitable in our case study (as discussed in Section 2),
extending the proposed method to data-adaptively identified subgroups warrant
future research. Data-adaptive subgroup identification approaches include
varying coefficient model (Chen and He, 2018), regression tree method
(Lipkovich et al., 2011) and fused lasso (Ma and Huang, 2017). When working
with data-adaptively identified subgroups, one needs to not only adjust for
the regularization and selection biases, but also account for the randomness
induced by subgroup identification. We leave the extension of the proposed
method to future research as the primary objective of this manuscript is to
investigate the heterogeneous causal effect of statin use on T2D risk.

#### 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.

## References

  * 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.
  * Burke et al. (2015) Burke, J. F., Sussman, J. B., Kent, D. M., and Hayward, R. A. (2015).  Three simple rules to ensure reasonably credible subgroup analyses.  Bmj, 351.
  * Chen and He (2018) Chen, X. and He, Y. (2018).  Inference of high-dimensional linear models with time-varying coefficients.  Statistica Sinica, pages 255–276.
  * Consortium et al. (2012) Consortium, I.-. R. M. R. A. I. M. et al. (2012).  The interleukin-6 receptor as a target for prevention of coronary heart disease: a mendelian randomisation analysis.  The Lancet, 379(9822):1214–1224.
  * 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.
  * 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.
  * 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.
  * Fuentes et al. (2018) Fuentes, C., Casella, G., and Wells, M. T. (2018).  Confidence intervals for the means of the selected populations.  Electronic Journal of Statistics, 12(1):58–79.
  * Goodarzi et al. (2013) Goodarzi, M. O., Li, X., Krauss, R. M., Rotter, J. I., and Chen, Y.-D. I. (2013).  Relationship of sex to diabetes risk in statin trials.  Diabetes Care, 36(7):e100–e101.
  * 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.
  * Guo et al. (2021) Guo, X., Wei, L., Wu, C., and Wang, J. (2021).  Sharp inference on selected subgroups in observational studies.  arXiv preprint arXiv:2102.11338.
  * 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.
  * Holland (1986) Holland, P. W. (1986).  Statistics and causal inference.  Journal of the American statistical Association, 81(396):945–960.
  * Hong et al. (2018) Hong, L., Kuffner, T. A., and Martin, R. (2018).  On overfitting and post-selection uncertainty assessments.  Biometrika, 105(1):221–224.
  * Ioannidis et al. (2019) Ioannidis, J. P., Tan, Y. J., and Blum, M. R. (2019).  Limitations and misinterpretations of e-values for sensitivity analyses of observational studies.  Annals of internal medicine, 170(2):108–111.
  * Kang et al. (2016) Kang, H., Zhang, A., Cai, T. T., and Small, D. S. (2016).  Instrumental variables estimation with some invalid instruments and its application to mendelian randomization.  Journal of the American statistical Association, 111(513):132–144.
  * Lango et al. (2008) Lango, H., Palmer, C. N., Morris, A. D., Zeggini, E., Hattersley, A. T., McCarthy, M. I., Frayling, T. M., and Weedon, M. N. (2008).  Assessing the combined impact of 18 common genetic variants of modest effect sizes on type 2 diabetes risk.  Diabetes, 57(11):3129–3135.
  * Li (2020) Li, S. (2020).  Debiasing the debiased lasso with bootstrap.  Electronic Journal of Statistics, 14(1):2298–2337.
  * Liao et al. (2019) Liao, K. P., Sun, J., Cai, T. A., Link, N., Hong, C., Huang, J., Huffman, J. E., Gronsbell, J., Zhang, Y., Ho, Y.-L., et al. (2019).  High-throughput multimodal automated phenotyping (map) with application to phewas.  Journal of the American Medical Informatics Association, 26(11):1255–1262.
  * Lipkovich et al. (2011) Lipkovich, I., Dmitrienko, A., Denne, J., and Enas, G. (2011).  Subgroup identification based on differential effect search—a recursive partitioning method for establishing response to treatment in patient subpopulations.  Statistics in Medicine, 30(21):2601–2621.
  * Lu et al. (2018) Lu, M., Sadiq, S., Feaster, D. J., and Ishwaran, H. (2018).  Estimating individual treatment effect in observational data using random forest methods.  Journal of Computational and Graphical Statistics, 27(1):209–219.
  * Ma and Huang (2017) Ma, S. and Huang, J. (2017).  A concave pairwise fusion approach to subgroup analysis.  Journal of the American Statistical Association, 112(517):410–423.
  * Macedo et al. (2014) Macedo, A. F., Douglas, I., Smeeth, L., Forbes, H., and Ebrahim, S. (2014).  Statins and the risk of type 2 diabetes mellitus: cohort study using the uk clinical practice pesearch datalink.  BMC cardiovascular disorders, 14(1):1–12.
  * Magnusson and Turnbull (2013) Magnusson, B. P. and Turnbull, B. W. (2013).  Group sequential enrichment design incorporating subgroup selection.  Statistics in medicine, 32(16):2695–2714.
  * Mansi et al. (2015) Mansi, I., Frei, C. R., Wang, C.-P., and Mortensen, E. M. (2015).  Statins and new-onset diabetes mellitus and diabetic complications: a retrospective cohort study of us healthy adults.  Journal of general internal medicine, 30(11):1599–1610.
  * Mora et al. (2010) Mora, S., Glynn, R. J., Hsia, J., MacFadyen, J. G., Genest, J., and Ridker, P. M. (2010).  Statins for the primary prevention of cardiovascular events in women with elevated high-sensitivity c-reactive protein or dyslipidemia: results from the justification for the use of statins in prevention: An intervention trial evaluating rosuvastatin (jupiter) and meta-analysis of women from primary prevention trials.  Circulation, 121(9):1069–1077.
  * 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.
  * Nissen et al. (2005) Nissen, S. E., Tuzcu, E. M., Schoenhagen, P., Crowe, T., Sasiela, W. J., Tsai, J., Orazem, J., Magorien, R. D., O’Shaughnessy, C., and Ganz, P. (2005).  Statin therapy, ldl cholesterol, c-reactive protein, and coronary artery disease.  New England Journal of Medicine, 352(1):29–38.
  * Packard et al. (1998) Packard, C., Shepherd, J., Cobbe, S., Ford, I., Isles, C., McKillop, J., Macfarlane, P., Lorimer, A., and Norrie, J. (1998).  Influence of pravastatin and plasma lipids on clinical events in the west of scotland coronary prevention study (woscops).  Circulation, 97(15):1440–1445.
  * Pan et al. (2020) Pan, W., Sun, W., Yang, S., Zhuang, H., Jiang, H., Ju, H., Wang, D., and Han, Y. (2020).  Ldl-c plays a causal role on t2dm: a mendelian randomization analysis.  Aging (Albany NY), 12(3):2584.
  * Park and Hastie (2007) Park, M. Y. and Hastie, T. (2007).  $l_{1}$ regularization path algorithm for generalized linear models.  Journal of the Royal Statistical Society. Series B (Statistical Methodology), 69(4):659–677.
  * Povsic et al. (2017) Povsic, T. J., Scott, R., Mahaffey, K. W., Blaustein, R., Edelberg, J. M., Lefkowitz, M. P., Solomon, S. D., Fox, J. C., Healy, K. E., Khakoo, A. Y., et al. (2017).  Navigating the future of cardiovascular drug development—leveraging novel approaches to drive innovation and drug discovery: summary of findings from the novel cardiovascular therapeutics conference.  Cardiovascular drugs and therapy, 31(4):445–458.
  * Rajpathak et al. (2009) Rajpathak, S. N., Kumbhani, D. J., Crandall, J., Barzilai, N., Alderman, M., and Ridker, P. M. (2009).  Statin therapy and risk of developing type 2 diabetes: a meta-analysis.  Diabetes care, 32(10):1924–1929.
  * Rosenkranz (2016) Rosenkranz, G. K. (2016).  Exploratory subgroup analysis in clinical trials by model selection.  Biometrical Journal, 58(5):1217–1228.
  * 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.
  * Skourtis et al. (2020) Skourtis, D., Stavroulaki, D., Athanasiou, V., Fragouli, P. G., and Iatrou, H. (2020).  Nanostructured polymeric, liposomal and other materials to control the drug delivery for cardiovascular diseases.  Pharmaceutics, 12(12):1160.
  * 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.
  * Sur and Candès (2019) Sur, P. and Candès, E. (2019).  A modern maximum-likelihood theory for high-dimensional logistic regression.  Proceedings of the National Academy of Sciences of the United States of America, 116(29):14516–14525.
  * Swerdlow (2016) Swerdlow, D. I. (2016).  Mendelian randomization and type 2 diabetes.  Cardiovascular Drugs and Therapy, 30(1):51–57.
  * Swerdlow et al. (2015) Swerdlow, D. I., Preiss, D., Kuchenbaecker, K. B., Holmes, M. V., Engmann, J. E., Shah, T., Sofat, R., Stender, S., Johnson, P. C., Scott, R. A., et al. (2015).  Hmg-coenzyme a reductase inhibition, type 2 diabetes, and bodyweight: evidence from genetic analysis and randomised trials.  The Lancet, 385(9965):351–361.
  * Thomas and Bornkamp (2017) Thomas, M. and Bornkamp, B. (2017).  Comparing approaches to treatment effect estimation for subgroups in clinical trials.  Statistics in Biopharmaceutical Research, 9(2):160–171.
  * VanderWeele and Ding (2017) VanderWeele, T. J. and Ding, P. (2017).  Sensitivity analysis in observational research: introducing the e-value.  Annals of internal medicine, 167(4):268–274.
  * 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.
  * Waters et al. (2013) Waters, D. D., Ho, J. E., Boekholdt, S. M., DeMicco, D. A., Kastelein, J. J., Messig, M., Breazna, A., and Pedersen, T. R. (2013).  Cardiovascular event reduction versus new-onset diabetes during atorvastatin therapy: effect of baseline risk factors for diabetes.  Journal of the American College of Cardiology, 61(2):148–152.
  * Windmeijer et al. (2019) Windmeijer, F., Farbmacher, H., Davies, N., and Davey Smith, G. (2019).  On the use of the lasso for instrumental variables estimation with some invalid instruments.  Journal of the American Statistical Association, 114(527):1339–1350.
  * Würtz et al. (2016) Würtz, P., Wang, Q., Soininen, P., Kangas, A. J., Fatemifar, G., Tynkkynen, T., Tiainen, M., Perola, M., Tillin, T., Hughes, A. D., et al. (2016).  Metabolomic profiling of statin use and genetic inhibition of hmg-coa reductase.  Journal of the American college of cardiology, 67(10):1200–1210.
  * Xue and Pan (2020) Xue, H. and Pan, W. (2020).  Inferring causal direction between two traits in the presence of horizontal pleiotropy with gwas summary data.  PLoS genetics, 16(11):e1009105.
  * 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.
  * 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.