# Splitting strategies for post-selection inference

Daniel G. Rasines111Corresponding author: daniel.garcia-
<EMAIL_ADDRESS>G. Alastair Young

###### Abstract

We consider the problem of providing valid inference for a selected parameter
in a sparse regression setting. It is well known that classical regression
tools can be unreliable in this context due to the bias generated in the
selection step. Many approaches have been proposed in recent years to ensure
inferential validity. Here, we consider a simple alternative to data splitting
based on randomising the response vector, which allows for higher selection
and inferential power than the former and is applicable with an arbitrary
selection rule. We provide a theoretical and empirical comparison of both
methods and extend the randomisation approach to non-normal settings. Our
investigations show that the gain in power can be substantial.

Keywords: Data splitting; Randomisation; Post-selection inference; Regression;
Variable selection.

##  1 Introduction

Suppose we have data $Y\sim N(\mu,\sigma^{2}I_{n})$, where
$\mu\in\mathbb{R}^{n}$, $\sigma^{2}>0$, and $I_{n}$ is the $n\times n$
identity matrix. We assume that the components of $\mu$ are modelled as a
function of $p$ covariates, $\mu_{i}=g(x_{i1},\ldots,x_{ip})$ for some unknown
$g$, and denote by $X=(x_{ij})\in\mathbb{R}^{n\times p}$ the known, fixed
design matrix. In many situations, it is suspected that only a few covariates
are truly active, and a preliminary variable-selection step is performed to
identify these. Having screened a set of potentially relevant variables, we
may want to provide inference for the regression coefficients of the best
linear approximation of $\mu$ in the selected model, or some other parameter
depending on the output of the selection step.

If the same data that was used for selection is also used to provide inference
for the selected parameter, standard inferential procedures are unreliable,
typically leading to overoptimistic results; see e.g. Hong et al. (2018). Data
splitting techniques, whereby a portion of the data is reserved for
uncertainty quantification, offer a simple yet effective way to circumvent
this problem. Unfortunately, data splitting often leads to procedures with
little power both for identifying the active covariates and for providing
inference for the selected parameters. In this paper, we study an alternative
to data splitting, motivated by the work by Tian and Taylor (2018), which
provides a more efficient way of splitting the sample information, resulting
in more powerful procedures, and which is easy to apply with a general
variable-selection method. Our main objective here is to provide a theoretical
and empirical comparison between the two splitting strategies.

Data splitting is a popular tool in prediction problems, where the hold-out
observations are used to assess the accuracy of a predictive model. When the
goal is inference rather than prediction, a frequent criticism of data
splitting is that different splits can produce different selected models, and
therefore two people analysing the same data may end up answering different
questions. While this is a valid concern, we stress that selection based on
the full data is not free from some level of arbitrariness, as the selection
process always involves subjective decisions, including the choice of the
selection rule itself, and in many cases a random input independent of the
data; see e.g. Wasserman and Roeder (2009); Meinshausen and Bühlmann (2010);
and Candès et al. (2018). Nevertheless, it is important to keep the effect of
the random components as low as possible, as failure to do so results in high
uncertainty about the relevance of the selected variables.

The splitting technique considered in this paper can be viewed as a variant of
data splitting which produces datasets which are more similar to the full
sample than those resulting from the latter method, and is therefore
potentially less affected by randomness. Furthermore, it entails no extra
computational cost with respect to data splitting. The method operates by
applying the variable selection algorithm to a randomised version of the data,
and then basing inference on the conditional distribution of the data given
its randomised form, thereby avoiding any selection bias. The general idea of
basing selection on an artificial perturbation of the data in this context was
proposed by Tian and Taylor (2018) as a way of deriving uniformly consistent
and powerful inferential procedures, and has become a popular device in the
literature. In the original paper, the authors compare the resulting
inferential power of randomisation and data splitting in circumstances where
it is possible to provide inference conditionally on the selection event,
showing the superiority of the former. Here we compare the methods in
circumstances where inference discards all the information of the selection
split, and is therefore unaffected by the complexity of the selection rule.
This is explained in more detail in §3, but the key point is that methods
discussed here may be applied in circumstances where the conditional sampling
distribution of the data given selection is intractable or unknown.

A large number of methods have been proposed in recent years to deal with
selection bias. They can be broadly divided into two categories: those which
assume that the selection algorithm is of a specific form; and those which
provide guarantees for an arbitrary selection rule. An important class of
methods in the first group is formed by conditional procedures, as considered
above, which are constructed by analysing the conditional distribution of the
data given the specified selection event, when this can be found explicitly.
This line of work was started by Lockhart et al. (2014) and has subsequently
been extended in multiple works such as Lee and Taylor (2014); Loftus and
Taylor (2014); Lee et al. (2016); Fithian et al. (2017); Tibshirani et al.
(2018); and Panigrahi et al. (2020). The second group of methods includes the
POSI approach of Berk et al. (2013) and extensions of it (Bachoc et al., 2017,
2020), which achieve uniformly valid inference by maximising over all possible
model selection procedures, and are very conservative as a result, as well as
the data splitting approach of Rinaldo et al. (2019), which provides model-
free procedures with asymptotically valid guarantees in a random-design
setting. Cox (1975) analysed data splitting in a simple inferential problem
involving many normal means and found it to be competitive against a natural
alternative. Methods based on data splitting have also been considered in more
recent works such as Rubin et al. (2006), Wasserman and Roeder (2009),
Ignatiadis et al. (2016), and DiCiccio et al. (2020). Fithian et al. (2017)
observed that inference after data splitting based only on the hold-out
observations is inadmissible, being always dominated by data-carving rules,
which consider the sampling distribution of the full data, conditional on the
selection event. Such data-carving may, however, be very complicated or
intractable to implement in many situations due to the complexity of the
conditional distribution.

##  2 Post-selection inference

Suppose that, for a given data vector $Y$, a variable-selection algorithm
selects a subset $s\subseteq\\{1,\ldots,p\\}$ of the covariates. In general,
determining an appropriate inferential objective post-selection is not
straightforward. Under a linearity assumption, $\mu=X\beta$ for some
$\beta=(\beta_{1},\ldots,\beta_{p})^{T}\in\mathbb{R}^{p}$, a natural
possibility is to provide inference for the components of $\beta$ associated
with the selected variables, $\\{\beta_{i}\colon i\in s\\}$. However, this is
a difficult problem when $p>n$ as the model is not identifiable. A popular
alternative target of inference is the projection parameter, proposed
originally by Berk et al. (2013). The projection parameter is the regression
parameter of $\mu$ projected onto the subspace spanned by the selected columns
of $X$:

$\beta_{s}(X)=\arg\min_{z\in\mathbb{R}^{|s|}}E\left(\|Y-X(s)z\|^{2}\right)=\\{X(s)^{T}X(s)\\}^{-1}X(s)^{T}\mu,$
(1)

where $X(s)$ is the submatrix of $X$ that contains the selected columns and
$|s|$ denotes the number of selected covariates. In other words, it is the
best linear predictor of $\mu$ in the selected model. When the model happens
to be linear in the selected covariates, $\mu=X(s)\beta_{s}$ for some
$\beta_{s}\in\mathbb{R}^{|s|}$, $\beta_{s}(X)$ is simply $\beta_{s}$.
Otherwise, the interpretation of the projection parameter is less transparent:
the $j$-th component of $\beta_{s}(X)$ may be viewed as the average change of
the response when the $j$-th selected covariate increases by one unit,
approximated in the selected model. An alternative interpretation is given in
§6.5 (see also §3.2 of Berk et al. (2013)). When the full model is not linear,
one can also consider the projection parameter using the full design matrix,
$\beta^{F}(X)=(X^{T}X)^{-1}X^{T}\mu$ (assuming $X$ has full column rank). Some
authors refer to this parameter as the full target, and to the previous one as
the partial target. In this case, post-selection inference may be provided for
the components of $\beta^{F}(X)$ associated with the selected variables,
$\\{\beta^{F}(X)_{i}\colon i\in s\\}$. Another option, presumably more common
in practice, is to proceed under the assumption that $\mu=X(s)\beta_{s}$ and
carry out inference on $\beta_{s}$. In the random-design case, Rinaldo et al.
(2019) develop inferential methods for other choices of the selected parameter
which depend on the distribution of the covariates.

Let us denote a generic selected parameter, possibly depending on the design
matrix as well as on the selected set, by $h_{s}(X)$. Adopting the conditional
approach (Fithian et al., 2017), we deem an inferential statement about
$h_{s}(X)$ valid if its error guarantees hold under the conditional
distribution of the data given the event that $h_{s}(X)$ was selected. For
example, a $1-\alpha$ confidence set $T$ for $h_{s}(X)$ is valid if it
satisfies

$\operatorname*{\mathrm{pr}}\left(h_{s}(X)\in T\mid S=s\right)\geq 1-\alpha,$
(2)

where $S$ is the random set of selected covariates. For simplicity, we will
assume that the choice of the interest parameters depends only on $s$, so that
inference on a given $h_{s}(X)$ is required if and only if $S=s$.

For some popular variable-selection algorithms, such as the lasso or stepwise
procedures with fixed tuning parameters, the conditioning event $\\{S=s\\}$
can be studied analytically (often, it can be written as a union of affine
sets; see e.g. Lee and Taylor (2014); Loftus and Taylor (2014); and Lee et al.
(2016)). In most cases, however, this event is too complicated to be explored
analytically. Furthermore, even when they can be implemented, conditional
methods tend to be very conservative (Kivaranovic and Leeb, 2020). In such
cases, data splitting offers an analytically simple and computationally light
solution to the inference problem.

##  3 Splitting methods

The most common form of data splitting is simple data splitting. Here, a
fraction $f=n_{1}/n$, $1\leq n_{1}<n$, is specified, and a set of indices $r$
is chosen uniformly at random from the subsets of $\\{1,\ldots,n\\}$ of size
$n_{1}$. Then, the sets of observations $(Y_{i}\colon i\in r)$ and
$(Y_{i}\colon i\in r^{c})$ are respectively used for selection and for
inference, where $r^{c}$ denotes the complement of $r$. Since all the
observations are independent, the conditional distribution of the inference
set given the output of the selection step is the same as the unconditional
one, so classical procedures can be used to provide valid inference for a
selected parameter. More elaborated splitting rules can be found in the
prediction literature; see e.g. Reitermanová (2010). These rules allocate the
samples to the selection and inferential sets according the observed values of
the covariates, usually trying to divide them as evenly as possible, to ensure
that the analyst has access to similar regions of the design space in both
stages. For example, when the observations come from different groups,
stratified data splitting rules allocate the samples so that both the
inference and the selection sets contain observations from all the observed
groups. Quite generally, then, a data splitting strategy can be formalised as
a random variable $R$, possibly depending on $X$, taking values in the power
set of $\\{1,\ldots,n\\}$.

In the current context, data splitting techniques can be viewed as particular
cases of a more flexible method. Suppose that $W$ is a random quantity,
possibly depending on $X$, and that in the selection step we only allow
ourselves to observe the value of a function $U\equiv u(Y,W)$. Since selection
depends on the data only through $U$, inference based on the conditional
distribution of the data given its observed value, $Y\mid U$, is free of
selection bias, and does not require knowledge of the selection mechanism.
Note that, by contrast, a conditional approach would base inference on
$Y\mid\\{S=s\\}$. Abstractly, we may think of the data under this conditional
distribution as the portion reserved for inference. Data splitting techniques
can be written this way by letting $W=R$ and $U$ be the set of observations
allocated to selection. A natural question to examine, then, is whether there
exist data decompositions of this form that are better than data splitting.

Tian and Taylor (2018) proposed randomisation schemes of the form $U=Y+\gamma
W$, where $W$ is $n$-dimensional artificial noise and $\gamma>0$ is a tuning
parameter which controls the amount of information reserved for inference.
Small values of $\gamma$ assign most of the sample information for selection,
while large values allocate most of it for inference. One clear advantage of
this approach over data splitting is that it gives access to all the observed
values of the covariates both at the selection and at the inferential stages,
while in data splitting we only have access to a subset of them at each stage.
To ensure high inferential power, Tian and Taylor (2018) recommend that the
distribution of $W$ has tails at least as heavy as the normal distribution. If
the observation variance $\sigma^{2}$ is known, the choice $W\sim
N(0_{n},\sigma^{2}I_{n})$, where $0_{n}$ is an $n$-dimensional vector of
zeroes, allows for a remarkably simple analysis. With this choice, $U\sim
N(\mu,(1+\gamma^{2})\sigma^{2}I_{n})$ is of the same parametric form as the
data. Furthermore, basing inference for the selected parameter on $Y\mid U$
amounts to basing it on the marginal distribution of $V=Y-\gamma^{-1}W\sim
N(\mu,(1+\gamma^{-2})\sigma^{2}I_{n})$. This follows because $U$ and $V$ are
independent, as they are uncorrelated and normal, and jointly sufficient.
Henceforth, we refer to this split as the $(U,V)$ decomposition. Tian (2020)
used such a decomposition to derive an unbiased estimator of the prediction
risk in a related context.

If the error variance is unknown but can be estimated with reasonable
precision, an approximate decomposition can be achieved by plugging in the
variance estimate in the variance of $W$. In section 5 we consider the
asymptotic validity of this approach. In the linear case, if $p$ is small
relative to $n$, $\sigma^{2}$ can be estimated in the classical way. Otherwise
we have to resort to high-dimensional alternatives. In our simulation studies
we used the estimator implemented in the selectiveInference R package of
Tibshirani et al. (2019), which estimates $\sigma^{2}$ using the residual sum
of squares from a lasso fit, with the penalty parameter tuned by cross-
validation. The good performance of this estimator in sparse models was
demonstrated in Reid et al. (2016). Many other methods are available; e.g. Fan
et al. (2012) and Bayati et al. (2013).

##  4 Theoretical comparison

In the low-dimensional linear case, $\mu=X\beta$, the optimality of a design
$X$ is commonly measured through summary statistics of the inverse Fisher
information, $\sigma^{2}(X^{T}X)^{-1}$. Popular criteria in the literature
include the D-optimality criterion, which minimises the determinant of the
inverse information, the A-optimality criterion, which minimises the trace,
and a few others. When a data split $r$ is used to partition the sample, the
inverse information matrices of the selection and inference sets are
$\sigma^{2}(X_{r}^{T}X_{r})^{-1}$ and
$\sigma^{2}(X_{r^{c}}^{T}X_{r^{c}})^{-1}$, where $X_{r}$ and $X_{r^{c}}$ are
the design matrices that contain only the observations in splits $r$ and
$r^{c}$, respectively. Correspondingly, an additive randomisation scheme
$U=Y+\gamma W$, $W\sim N(0_{n},\sigma^{2}I_{n})$, produces splits of the
inverse Fisher information $(1+\gamma^{2})\sigma^{2}(X^{T}X)^{-1}$ and
$(1+\gamma^{-2})\sigma^{2}(X^{T}X)^{-1}$.

Proposition 1, which we state for a generic parametric model, shows that the
latter split is more efficient than a class of data splitting strategies when
the quality of the inverse Fisher information is measured in a particular way.
For this analysis, we let $Y_{1},\ldots,Y_{n}$ be independent (not necessarily
identically distributed) observations from a generic regular parametric model
with $p$-dimensional parameter $\beta$, and write $I_{i}(\beta)$ for the
Fisher information provided by the $i$-th observation, and
$\mathcal{I}_{Y}(\beta)=\sum_{i=1}^{n}I_{i}(\beta)$ for the full-data Fisher
information. A generic data split $r$ distributes the information between the
selection and inferential tasks as

$\mathcal{I}_{Y}(\beta)=\sum_{i\in r}I_{i}(\beta)+\sum_{i\in
r^{c}}I_{i}(\beta)\equiv\mathcal{I}_{r}(\beta)+\mathcal{I}_{r^{c}}(\beta).$
(3)

The result concerns the behaviour of $\mathcal{I}_{R}(\beta)$ and
$\mathcal{I}_{R^{c}}(\beta)$ under repeated sampling from a random data
splitting strategy $R$, such as simple data splitting.

In what follows, all expectations are taken with respect to the distribution
of $R$, and, for two positive definite matrices $A$ and $B$, we use the
notation $A<B$ ($\leq$) to indicate that $A-B$ is negative (semi)definite.
Furthermore, we say that a real-valued matrix function $\varphi$ is convex if
it satisfies $\varphi\\{tA+(1-t)B\\}\leq t\varphi(A)+(1-t)\varphi(B)$ for all
matrices $A$ and $B$ in its domain and $t\in(0,1)$, and we say that it is
strictly increasing if $A<B$ implies $\varphi(A)<\varphi(B)$.

######  Proposition 1.

Let $R$ be any splitting strategy satisfying $\operatorname*{\mathrm{pr}}(i\in
R)=f$ for some $0<f<1$ and all $i=1,\ldots,n$, and let $\varphi$ be a real-
valued function defined on the set of $p\times p$ positive definite matrices
which is convex and strictly increasing. Furthermore, assume that
$\mathcal{I}_{r}(\beta)$ and $\mathcal{I}_{r^{c}}(\beta)$ are invertible for
all $r$ in the support of $R$, and that
$\mathcal{I}_{r_{1}}(\beta)\neq\mathcal{I}_{r_{2}}(\beta)$ for some $r_{1}$
and $r_{2}$. Then,

$\varphi\left\\{\frac{1}{f}\mathcal{I}_{Y}(\beta)^{-1}\right\\}<E\left[\varphi\left\\{\mathcal{I}_{R}(\beta)^{-1}\right\\}\right],\quad\varphi\left\\{\frac{1}{1-f}\mathcal{I}_{Y}(\beta)^{-1}\right\\}<E\left[\varphi\left\\{\mathcal{I}_{R^{c}}(\beta)^{-1}\right\\}\right].$
(4)

Examples of functions $\varphi$ that satisfy the conditions of Proposition 1
include: $\varphi(A)=\text{tr}(A)$, related to the average variance of the
estimates of the regression coefficients; $\varphi(A)=v^{T}Av$ for any vector
$v\in\mathbb{R}^{p}$, related to the estimation variance of a specific
combination of parameters $v^{T}\beta$; $\varphi(A)=\max\\{\text{diag}(A)\\}$,
related to the maximum estimation variance of the regression coefficients; and
$\varphi(A)=\lambda_{\text{max}}(A)$, where $\lambda_{\text{max}}(A)$ is the
maximum eigenvalue of $A$, related to the variance estimation of the first
principal component.

Proposition 1 says that, if we measure the quality of a design by such a
$\varphi$, a splitting scheme $u(Y,W)$ which splits the Fisher information as
$\mathcal{I}_{Y}(\beta)=f\mathcal{I}_{Y}(\beta)+(1-f)\mathcal{I}_{Y}(\beta)$
is more efficient, on average, than any data splitting strategy for which
every sample has the same probability of being allocated to either set.

In additive models, $Y=\mu+\varepsilon$, splitting strategies satisfying the
former condition are trivially given by rules of the form $U=Y+\gamma W$,
where $W$ has $E(W)=0_{n}$ and
$\operatorname*{\mathrm{var}}(W)=\operatorname*{\mathrm{var}}(\varepsilon)$,
and $\gamma$ satisfies $(1+\gamma^{2})^{-1}=f$. On the other hand, the set of
data splitting strategies for which $\operatorname*{\mathrm{pr}}(i\in R)$ is
constant include: simple data splitting with splitting fraction $f=n_{1}/n$;
the class of stratified sampling rules where the observations are divided into
$k$ equal-size sets, and the selection set is formed by choosing a fixed
number of observations uniformly at random from each set; and any splitting
rule with
$\operatorname*{\mathrm{pr}}(R=A)=\operatorname*{\mathrm{pr}}(R=A^{c})=1/2$
for some $A\subseteq\\{1,\ldots,n\\}$. For example, we might use a rule that
divides the samples into two sets as evenly as possible, and then label them
as selection and inference sets at random.

##  5 Generalised $(U,V)$ decomposition

In additive models, efficient splits are easily achievable via additive
perturbation of the data. However, when the observation error and the
randomisation noise are not normal, the distributions of $U$ and $Y|U$ are
generally not available in closed form, complicating the selection and
inferential analyses. Furthermore, when the observation variance has to be
estimated or the normality assumption is mildly violated, basing inference on
the marginal distribution of $V$, as described above, is not formally
justified. In this section we provide a set of conditions under which the
$(U,V)$ decomposition is asymptotically valid in more general settings.

Suppose that the model is $Y=\mu+\varepsilon$,
$\varepsilon=(\varepsilon_{1},\ldots,\varepsilon_{n})^{T}$, where the errors
are a random sample from an unknown distribution, not necessarily normal. All
the selected parameters considered in §2 can be written as $\eta^{T}\mu$ for
some contrast $\eta\in\mathbb{R}^{n}$. If the errors are $N(0,\sigma^{2})$ and
we know $\sigma^{2}$, exact post-selection inference for $\eta^{T}\mu$ after
selection based on $U=Y+\gamma W$ is provided via $\eta^{T}V\sim
N(\eta^{T}\mu,(1+\gamma^{-2})\sigma^{2}\|\eta\|^{2})$, where $\|\cdot\|$
denotes the Euclidean norm. Theorem 1 gives conditions under which this
approach is asymptotically valid when the errors are not normal and the
observation variance is unknown but can be estimated consistently.

The result considers asymptotic regimes with a fixed number of covariates and
a growing number of observations. Unless otherwise stated, all the elements
involved in the analysis depend on the sample size $n$. This dependence is not
made explicit to avoid saturating the notation. Also, for a vector or matrix
$A$, $\max(A)$ denotes the maximum absolute entry of $A$, and, for a square
matrix $A$, $\lambda_{\min}(A)$ and $\lambda_{\max}(A)$ denote the smallest
and largest eigenvalues of $A$.

######  Theorem 1.

Let $p$ be fixed, $n\to\infty$, and consider a sequence of selection sets
$s_{1},s_{2},\ldots\subseteq\\{1,\ldots,p\\}$ such that the sequence selection
probabilities, $\\{\operatorname*{\mathrm{pr}}(S=s_{i})\colon i\geq 1\\}$, is
bounded away from zero. Assume that all the selection events $\\{S=s_{i}\\}$
can be written in the form $\\{X^{T}U\in E\\}$, where $U=Y+\gamma W$ for a
fixed $\gamma\in(0,1)$, and $E$ is a finite union of convex sets. Furthermore,
assume that $X^{T}X$ is invertible for all $n$,
$\lambda_{\max}\\{(X^{T}X)^{-1}\\}=O(n^{-1})$, $\max(X)=o(n^{1/6})$, and that
the contrasts $\eta$ satisfy $\max(\eta)\|\eta\|^{-1}=O(n^{-1/2})$. If the
observation errors $\varepsilon_{1},\ldots,\varepsilon_{n}$ are independent
and identically distributed with mean zero, variance $\sigma^{2}$, and
$E(|\varepsilon_{1}|^{3})<\infty$, and $W=\hat{\sigma}Z$, where $\hat{\sigma}$
is a weakly consistent estimator of $\sigma$ and $Z\sim N(0_{n},I_{n})$ is
independent of the data, then

$(1+\gamma^{-2})^{-1/2}\sigma^{-1}\|\eta\|^{-1}(\eta^{T}V-\eta^{T}\mu)|\\{S=s_{n}\\}\xrightarrow{d}N(0,1),\quad
n\to\infty,$ (5)

where $V=Y-\gamma^{-1}W$.

The requirement that selection depends on the data only through $X^{T}U$
allows us to reduce the dimension of the conditioning event from $n$ to $p$,
which stays bounded. It is satisfied by all penalised-likelihood methods,
which estimate the active set of covariates by the support of

$\hat{\beta}=\arg\min_{\beta\in\mathbb{R}^{p}}\|U-X\beta\|^{2}+l(\beta)$ (6)

for an appropriate loss function $l$. This follows because the first term can
be written, up to an additive constant independent of $\beta$, as
$\|X(X^{T}X)^{-1}X^{T}U-X\beta\|^{2}$.

Regarding the convexity assumption, suppose that, for a given $s$, the
selection event in terms of the full vector $U$ is
$E^{\prime}\subseteq\mathbb{R}^{n}$. Under the previous condition, the
selection event in terms of $X^{T}U$ is $E=X^{T}E^{\prime}=\\{X^{T}u\colon
u\in E^{\prime}\\}$. Therefore, if $E^{\prime}=E^{\prime}_{1}\cup\ldots\cup
E^{\prime}_{m}$, with $E^{\prime}_{1},\ldots,E^{\prime}_{m}$ convex,
$E=X^{T}E^{\prime}_{1}\cup\ldots\cup X^{T}E^{\prime}_{m}$ is also a union of
convex sets, so the requirement can be imposed on $E^{\prime}$ directly.
Verifying that $E^{\prime}$ is a union of convex sets can still be difficult
problem in general, but it has been shown to be true for some widely used
selection rules such as forward stepwise regression, the lasso, and least
angle regression, all with fixed tuning parameters (Tibshirani et al., 2018).
We note, however, that this requirement is mainly convenient from a
mathematical perspective and can probably be weakened.

The asymptotic condition on $\eta$ is very natural. It ensures that the
asymptotic support of $\eta/\|\eta\|$ is unbounded. When inference is sought
for a selected coefficient or for a projection parameter, the contrasts are of
the form $X(s)\\{X(s)^{T}X(s)\\}^{-1}e_{j}$ for some
$s\subseteq\\{1,\ldots,p\\}$ and $1\leq j\leq p$, where
$e_{j}\in\mathbb{R}^{p}$ is the $j$-th vector of the canonical basis. These
vectors have

$\displaystyle\|\eta\|^{2}$ $\displaystyle=$ $\displaystyle
e_{j}^{T}\\{X(s)^{T}X(s)\\}^{-1}e_{j}$ (7) $\displaystyle\geq$
$\displaystyle\lambda_{\min}\left[\\{X(s)^{T}X(s)\\}^{-1}\right]$ (8)
$\displaystyle=$ $\displaystyle\lambda_{\max}\\{X(s)^{T}X(s)\\}^{-1}$ (9)
$\displaystyle\geq$ $\displaystyle\lambda_{\max}\\{X^{T}X\\}^{-1}$ (10)

by the min-max theorem and the fact that
$\lambda_{\max}\\{X(s)^{T}X(s)\\}\leq\lambda_{\max}\\{X^{T}X\\}$. Furthermore,
by similar arguments and the Cauchy–Schwarz inequality, it follows that

$\displaystyle\max(\eta)$ $\displaystyle=$ $\displaystyle\max_{1\leq i\leq
n}e_{i}^{T}X(s)\\{X(s)^{T}X(s)\\}^{-1}e_{j}$ (11) $\displaystyle\leq$
$\displaystyle\max_{1\leq i\leq
n}\|e_{i}^{T}X(s)\|\|\\{X(s)^{T}X(s)\\}^{-1}e_{j}\|$ (12) $\displaystyle\leq$
$\displaystyle
p^{1/2}\max(X)\left[e_{j}^{T}\\{X(s)^{T}X(s)\\}^{-2}e_{j}\right]^{1/2}$ (13)
$\displaystyle\leq$ $\displaystyle
p^{1/2}\max(X)\lambda_{\max}\left[\\{X(s)^{T}X(s)\\}^{-1}\right]$ (14)
$\displaystyle=$ $\displaystyle\max(X)O(n^{-1})$ (15)

under the assumptions of the theorem. If we further assume that the entries of
$X$ are bounded as $n$ grows, $\max(X)=O(1)$, putting all together it follows
that the asymptotic condition on the contrast is satisfied if
$\lambda_{\max}(X^{T}X)=O(n)$.

##  6 Simulation study

###  6.1 Setting

We compared data splitting with the $(U,V)$ decomposition in the context of
the normal linear model. The data generating process was
$Y=X\beta+N(0_{n},I_{n})$, where $\beta\in\mathbb{R}^{p}$ is an unknown sparse
vector of coefficients. At each replication of the simulations, the rows of
the design matrix were generated as independent samples from the distribution
$N(0_{p},\Gamma)$, where $\Gamma\in\mathbb{R}^{p\times p}$ is a Toeplitz
matrix with $(i,j)$ entry $\rho^{|i-j|}$ for some $\rho\geq 0$. The
observation variance, fixed at $1$, was assumed unknown, and was estimated in
the classical way when $p<n/4$, and using the high-dimensional alternative
otherwise, as per §3. We compared the splitting strategies in terms of
selection power, selection stability, and inferential power in low and high-
dimensional settings. All the simulation results can be found in the Appendix.

The data splitting strategy considered was the DUPLEX (Snee, 1977). In the
even case ($f=1/2$), the DUPLEX proceeds as follows. First, it finds the two
covariate observations that are farthest apart (in Euclidean distance) and
assigns them to the selection set. Then, it finds the next two observations
that are farthest apart and assigns them to the inference set. Finally, it
allocates the remaining samples one at a time, rotating between the selection
and inference sets, by selecting the remaining sample that is farthest apart
from the rest of the observations in the given set, starting with the
selection set. If the selection and inferential sets have unequal sizes
($f\neq 1/2$), the algorithm is applied until the smallest set is filled, and
all the remaining observations are assigned to the other one. When the
covariate observations are not uniformly distributed, this algorithm ensures
that neither the selection nor the inferential splits are imbalanced. We note
that this algorithm does not satisfy the condition of Proposition 1, but it is
nevertheless dominated by the $(U,V)$ decomposition, as we shall see.

Following the discussion of section 4, for a given fraction $f$, we compared
the performance of the data splitting algorithm with splitting fraction $f$
with the performance of the randomised procedure with tuning parameter
$\gamma=(f^{-1}-1)^{1/2}$. Specifically, for a given $f$ and a given dataset
$(Y,X)$, we considered two procedures. The first bases selection on
$(Y^{(1)},X^{(1)})$ and inference on $(Y^{(2)},X^{(2)})$, where
$Y^{(1)}=(Y_{i}\colon i\in r)^{T}$, $X^{(1)}=(x_{ij}\colon i\in
r,j=1,\ldots,p)$, $Y^{(2)}=(Y_{i}\colon i\in r^{c})^{T}$, and
$X^{(2)}=(x_{ij}\colon i\in r^{c},j=1,\ldots,p)$. The second one bases
selection on $(U,X)=(Y+\gamma W,X)$ and inference on
$(V,X)=(Y-\gamma^{-1}W,X)$, with $W=\hat{\sigma}Z$, where $Z\sim
N(0_{n},I_{n})$ is artificially generated noise independent of the data and
$\hat{\sigma}^{2}$ is the estimate of model variance.

Regarding selection, we considered two recently proposed algorithms: the
fixed-$X$ knockoff algorithm of Barber and Candès (2015) and the stability
selector paired with the lasso of Meinshausen and Bühlmann (2010), as
implemented in the R packages knockoff (Barber et al., 2020) and stabs (Hofner
and Hothorn, 2017). These are selection rules for which conditional inference
is analytically intractable and for which computational approaches are very
demanding (Markovic et al., 2019), so they constitute an example where
splitting approaches would likely be preferred in practice. These algorithms
aim at identifying the set of active coefficients while keeping the number of
false discoveries under control. The knockoff provides a guarantee on the
false discovery rate, while stability selection controls the expected number
of false discoveries. In all the simulations we set the false discovery rate
of the knockoff algorithm at $0.3$, and the expected number of false
discoveries and cutoff threshold of the stability algorithm at $3$ at $0.7$
(see the package documentation for details). Since the fixed-$X$ knockoff is
only applicable in settings with less covariates than observations, it was
only considered in the lower dimensional cases.

###  6.2 Selection power

In this simulation, we generated $10^{3}$ triplets $(\beta,Y,X)$ independently
for each combination of the following parameters: $n=200$; $\rho=0,0.5$;
$f=1/2,3/4$ (correspondingly $\gamma=1,3^{-1/2}$); $p=30,50$ for the knockoff
algorithm, and $p=200,1000$ for the stability selection algorithm. For each
combination of the parameters and each repetition, the true $\beta$ was
generated by sampling $10$ non-zero positions uniformly at random and filling
them with independent random variables distributed uniformly in the set
$\\{-1,-0.9,-0.8,\ldots,-0.1,0.1,\ldots,0.9,1\\}$. Then, for each $\beta$, a
pair $(Y,X)$ was generated, and the corresponding selection algorithm was
applied to $(Y^{(1)},X^{(1)})$ in the case of data splitting, and to $(U,X)$
in the case of the randomised procedure.

The selection ability of the splitting methods was compared according to two
criteria: true positive rate and power. The true positive rate is the average
number of correct discoveries divided by the total number of active
covariates, and the power is defined, for each possible value of
$|\beta_{i}|=0.1,\ldots,1$, as the average number of times a coefficient with
absolute value $|\beta_{i}|$ is selected, averaged over all the generated
coefficients of all $\beta$’s. We also computed, for comparison, the results
corresponding to applying the selection algorithms to the full dataset,
$(Y,X)$.

Table 1 shows the observed true positive rates of data-splitting and
randomisation divided by the observed true positive rates of full-data
selection, and Figs. 1 and 2 show the empirical power functions of the three
methods. The results clearly favour the randomised split procedure over data
splitting, particularly in the cases with larger values of $p$. Quite
remarkably, despite the fact that the choices of $f$ and $\gamma$ were
balanced (in the sense described in §3), selection based on a randomised split
yielded a performance which was in some cases closer to full-data selection
than to data-splitting selection.

###  6.3 Selection stability

In addition to enjoying high power, it is important that the selection method
does not depend strongly on ancillary components of the analysis: in the case
of data splitting, on the choice of the selection and inference sets, and, in
the case of randomisation, on the observed value of the artificial noise $W$.
To measure the stability of the selection strategies with respect to these
elements, a simulation was conducted in which the selection algorithms were
applied to multiple splits of the same dataset. Since the DUPLEX is a
deterministic allocation rule, in this section we considered a simple data
splitting scheme instead.

We set $\beta=(1,0.9,\ldots,0.1,0,\ldots,0)^{T}$, $\rho=0.5$, $f=1/2,3/4$,
$p=50$ for the knockoff, and $p=400$ for stability selection. For each
$(f,p)$, we generated $100$ pairs $(Y,X)$ and, for each pair, we sampled $50$
selection sets uniformly at random and $50$ realisations of the randomisation
noise $W$, and applied the selection algorithms to each data split and
perturbed instance of the data. For each $(Y,X)$, we recorded the average
number of times each active covariate was selected across the different
splits, an estimate of $\operatorname*{\mathrm{pr}}(i\in S|Y,X)$,
$i=1,\ldots,10$. The empirical mean and standard deviations of the $100$
estimated averages can be found in Table 2.

Randomisation was more stable than data splitting, because the corresponding
selection probabilities of the active covariates were more concentrated around
higher values. For example, for the knockoff with $f=1/2$, for most values of
$(Y,X)$ we had more than a $90\%$ probability of selecting $\beta_{3}$ under
repeated sampling from $W$, as opposed to the $\sim 77\%$ of data splitting.

###  6.4 Inference for selected coefficients

In this simulation, we consider the problem of constructing confidence
intervals for the selected coefficients $\\{\beta_{i}\colon i\in s\\}$.
Firstly, we show that a face-value approach, which reports the standard
confidence intervals ignoring selection, undercovers coefficients with small
effect size $|\beta_{i}|$, while the splitting techniques provide valid
intervals. This gives an illustration of the need to take selection into
account in the inferential stage. Secondly, we compare the lengths of the
intervals derived from the randomised procedure with respect to those obtained
by data splitting. To avoid complications related to confidence interval
construction in non-identifiable settings, we shall only consider here low-
dimensional cases.

Assume that $X$ has full rank and let
$\hat{\beta}_{i}=e_{i}^{T}(X^{T}X)^{-1}X^{T}Y$ be the ordinary least squares
estimator of $\beta_{i}$. In a classical, non-selective setting, where the
inferential goals are determined prior to the data collection, marginal
inference for $\beta_{i}$ is based on the pivot
$T_{i}=(\hat{\beta}_{i}-\beta_{i})/\hat{\sigma}$, which follows a scaled $t$
distribution with $n-p$ degrees of freedom, where $\hat{\sigma}^{2}$ is the
classical variance estimator. If, however, inference on a given $\beta_{i}$ is
only provided for some data samples, $T_{i}$ is no longer pivotal, and the
resulting confidence intervals can be miscalibrated. To exemplify this, we run
the following simulation. We fixed $n=200$, $p=30$, true parameter
$\beta=(1,-1,0.5,-0.5,0.2,-0.2,0,0,\ldots,0)^{T}$, $\rho=0,0.5$, and
$f=1/2,3/4$, and for each combination of $(\rho,f)$, we generated $5\times
10^{3}$ pairs $(Y,X)$. We then applied the selection algorithm to each dataset
as in the previous sections and constructed classical equal-tailed confidence
intervals based on $T_{i}$ for each selected coefficient, with a nominal
coverage of $90\%$.

Tables 3 and 4 contain the observed coverages of these intervals in the rows
indicated by FV (face-value), averaged across coefficients with equal absolute
effect. Clearly, the face-value intervals are unreliable when $|\beta_{i}|=0$
or $0.2$, with actual coverages significantly lower than the nominal one in
most cases. Note that here the face-value approach is being implemented after
selection based on a fraction of the sample information, so the selection
effect is not as large as if all the data had been used for selection.

The rows labelled with HD contain the coverages of the intervals constructed
using only the hold-out split of the data. In the case of data splitting, the
intervals were derived from the $t$-pivots of the hold-out observations,
$\hat{\beta}^{\text{DS}}_{i}=e_{i}^{T}(X^{(2)T}X^{(2)})^{-1}X^{(2)T}Y^{(2)}$.
With the current simulation parameters, the number of remaining observations
is always larger than the number of covariates, so $X^{(2)}$ has full rank
almost surely. For the randomised procedure, we used the approximation
$\hat{\beta}^{\text{R}}_{i}=e_{i}^{T}(X^{T}X)^{-1}X^{T}V\dot{\sim}N(\beta_{i},(1+\gamma^{-2})\sigma^{2}e_{i}^{T}(X^{T}X)^{-1}e_{i})$.
The intervals were thus constructed by studentisation of
$\hat{\beta}^{\text{R}}_{i}$. Despite the distributional approximation, we see
that the coverages of the resulting confidence intervals were very close to
the nominal one.

For the hold-out methods we also recorded the average length of the intervals,
which can be found in Tables 5 and 6. The intervals provided by the randomised
procedure were always shorter, on average, than those provided by the data
splitting procedure. In the cases where the amount of information reserved for
inference was small ($f=3/4$), the matrices $X^{(2)}$ were of dimension
$50\times 30$, and the resulting intervals were very wide as a consequence.
The randomised procedure performed significantly better in these cases, giving
intervals roughly $75\%$ shorter than data splitting. The maximum observed
standard deviations of the figures shown in the tables were, respectively,
$0.9,1.9,0.004$, and $0.007$.

###  6.5 Inference for projection parameters

We now consider settings with a number of covariates exceeding the sample
size. In these situations the full model is not identifiable, and the methods
used in the previous section are not applicable. Instead, we consider the
problem of constructing confidence intervals for the coefficients of a
projection parameter. As noted by Berk et al. (2013), the coefficients of a
projection parameter for a selected set $s$ can be interpreted as a measure of
the correlation between the selected coefficients and the mean $\mu$,
accounting for all the other selected covariates. This follows because
$\beta_{s}(X)_{i}$ can be written as the projection of $\mu$ onto the residual
vector of the $i$-th selected covariate regressed against the other selected
covariates. Indeed, if $X(s)_{i}$ is the $i$-th column of $X(s)$ and
$X(s)_{-i}$ is $X(s)$ with the $i$-th column removed, a simple calculation
shows that

$\beta_{s}(X)_{i}=\frac{1}{\|r_{i}\|^{2}}r_{i}^{T}\mu,\quad
r_{i}=\left\\{I_{n}-X(s)_{-i}\left(X(s)_{-i}^{T}X(s)_{-i}\right)^{-1}X(s)_{-i}^{T}\right\\}X(s)_{i}.$
(16)

Thus, $\beta_{s}(X)_{i}$ provides a measure of the relevance of the $i$-th
covariate relative to the selected set, with a null value indicating that the
variable is redundant and can be safely removed from $s$.

If $s$ had been fixed in advance and we knew the value of $\sigma^{2}$,
inference on the coefficients of $\beta_{s}(X)$ would be based on the
components of $\hat{\beta}_{s}(X)=\\{X(s)^{T}X(s)\\}^{-1}X(s)^{T}Y\sim
N(\beta_{s}(X),\sigma^{2}\\{X(s)^{T}X(s)\\}^{-1})$. So, for a nominal coverage
of $1-\alpha$, the confidence intervals would be given by
$[\hat{\beta}_{s}(X)_{i}\mp
q_{1-\alpha/2}\sigma[e_{i}^{T}\\{X(s)^{T}X(s)\\}^{-1}e_{i}]^{1/2}]$, where
$q_{1-\alpha/2}$ is the $1-\alpha/2$ quantile of a standard normal
distribution. When $\sigma^{2}$ is unknown, we can plug-in the high-
dimensional estimate used for the $(U,V)$ decomposition,
$\hat{\sigma}^{2}_{\text{HD}}$, say, and report $[\hat{\beta}_{s}(X)_{i}\mp
q_{1-\alpha/2}\hat{\sigma}_{\text{HD}}[e_{i}^{T}\\{X(s)^{T}X(s)\\}^{-1}e_{i}]^{1/2}]$.

When $s$ is data-dependent, hold-out inference can be provided similarly, with
the difference that, for data splitting, the inferential target is not the
full-projection parameter $\beta_{s}(X)$, but the projection parameter based
on the hold-out observations. Inference is thus provided for
$\beta_{s}(X^{(2)})=\\{X^{(2)}(s)^{T}X^{(2)}(s)\\}^{-1}X^{(2)}(s)^{T}\mu_{2}$,
where $\mu_{2}=E(Y^{(2)})$, and is based on
$\hat{\beta}_{s}^{\text{DS}}(X^{(2)})=\\{X^{(2)}(s)^{T}X^{(2)}(s)\\}^{-1}X^{(2)}(s)^{T}Y^{(2)}\sim
N(\beta_{s}(X^{(2)}),\sigma^{2}\\{X^{(2)}(s)^{T}X^{(2)}(s)\\}^{-1})$. In the
case of randomisation, the target parameter is still $\beta_{s}(X)$, and, as
before, inference is based on the normal approximation to the distribution of

$\hat{\beta}_{s}^{\text{R}}(X)=\\{X(s)^{T}X(s)\\}^{-1}X(s)^{T}V\dot{\sim}N(\beta_{s}(X),(1+\gamma^{-2})\sigma^{2}\\{X(s)^{T}X(s)\\}^{-1}).$
(17)

In both cases, $\sigma^{2}$ was approximated by
$\hat{\sigma}^{2}_{\text{HD}}$.

The simulation parameters were set as in §6.4 except for the number of
covariates, which was $p=400$. Here the knockoff was not considered, as it is
not applicable with $p>n$. The results can be found in Tables 7 and 8. Note
that in this case the columns indicate the absolute value of the full-model
coefficients $\beta_{i}$, not the coefficients of the projection parameters,
which in general depend on $s$. The coverage results were similar to the lower
dimensional case, with the difference that the effect of selection was more
pronounced here. In one of the cases considered, the coefficients of the
projection parameters associated with the null covariates in the full model
were almost guaranteed to be missed by the face-value intervals. The hold-out
intervals, on the other hand, remained well-calibrated across the different
coefficients. Regarding interval lengths, in this case we computed the average
of all the intervals produced for a given selection set $s$, and averaged the
results over selection sets of equal size in order to establish a more
equitable comparison. We see that the average lengths were very similar in
this case. Recall, however, that the selection power of the latter method is
substantially higher in high-dimensional settings, so in conjunction
randomisation dominates data splitting. The maximum standard deviation of the
coverage figures were $1.4,4.5,1.2$ and $0.4$, for $|\beta_{i}|=0,0.2,0.5$,
and $1$, respectively. The maximum standard deviation of the length figures
was $0.04$. The entries with a dash indicate that no selection set of the
corresponding size was selected in the simulation.

##  7 Discussion

Randomisation offers an alternative to data splitting which divides the sample
information more efficiently. We have compared randomisation and data
splitting with respect to their selection stability, and their selection and
inferential power, showing that the dominance of the former can be substantial
in settings with a limited amount of information. As an example, consider a
linear model with $n=2p+2$ and splitting fraction $f=1/2$. Data splitting
would produce an inferential set with more variables than observations, while
randomisation would avoid this at the expense of larger per-observation error.
One limitation of randomisation, however, is that implementation requires an
effective estimate of the model variance.

Here we have only focused on the use of randomisation for inference after
variable selection. Naturally, this idea has potential applicability as a
power enhancement device in other problems where data splitting is normally
used. Despite its shortcomings, data splitting remains a widely used technique
in modern statistical practice thanks to its ability to handle complex
selection issues in a straightforward way, and randomisation opens a way to
mitigate some of its drawbacks.

## Appendix

### Simulation results

Table 1: True positive rate of the selection algorithms applied after data
splitting (DS) and randomisation (R), normalised by the true positive rate of
selection applied to the full dataset

Knockoff | Split
---|---
$f$ | $\rho$ | $p$ | DS | R
1/2 | 0 | 30 | 0.903 | 0.942
1/2 | 0 | 50 | 0.738 | 0.923
1/2 | 0.5 | 30 | 0.820 | 0.914
1/2 | 0.5 | 50 | 0.648 | 0.890
3/4 | 0 | 30 | 0.972 | 0.978
3/4 | 0 | 50 | 0.935 | 0.971
3/4 | 0.5 | 30 | 0.940 | 0.964
3/4 | 0.5 | 50 | 0.890 | 0.970
Stability | Split
---|---
$f$ | $\rho$ | $p$ | DS | R
1/2 | 0 | 200 | 0.694 | 0.867
1/2 | 0 | 1000 | 0.486 | 0.821
1/2 | 0.5 | 200 | 0.691 | 0.858
1/2 | 0.5 | 1000 | 0.478 | 0.820
3/4 | 0 | 200 | 0.895 | 0.948
3/4 | 0 | 1000 | 0.826 | 0.933
3/4 | 0.5 | 200 | 0.886 | 0.948
3/4 | 0.5 | 1000 | 0.826 | 0.934

Table 2: Means and standard deviations of the estimated selection
probabilities for fixed values of $(Y,X)$

|  | $\beta_{i}$
---|---|---
$f$ | Split | 1 | 0.9 | 0.8 | 0.7 | 0.6 | 0.5 | 0.4 | 0.3 | 0.2 | 0.1
|  | Knockoff
1/2 | DS | 0.81 | 0.76 | 0.77 | 0.74 | 0.71 | 0.67 | 0.62 | 0.53 | 0.38 | 0.15
(0.08) | (0.09) | (0.09) | (0.09) | (0.09) | (0.11) | (0.11) | (0.16) | (0.16) | (0.13)
R | 0.98 | 0.92 | 0.98 | 0.95 | 0.95 | 0.92 | 0.89 | 0.77 | 0.57 | 0.25
(0.06) | (0.18) | (0.06) | (0.10) | (0.11) | (0.13) | (0.13) | (0.23) | (0.26) | (0.22)
3/4 | DS | 0.98 | 0.95 | 0.97 | 0.97 | 0.96 | 0.88 | 0.89 | 0.80 | 0.62 | 0.25
(0.14) | (0.22) | (0.17) | (0.17) | (0.20) | (0.33) | (0.31) | (0.40) | (0.49) | (0.44)
R | 0.99 | 0.93 | 0.99 | 0.97 | 0.98 | 0.96 | 0.96 | 0.86 | 0.69 | 0.34
(0.04) | (0.19) | (0.04) | (0.10) | (0.08) | (0.10) | (0.11) | (0.22) | (0.29) | (0.32)
|  | Stability
1/2 | DS | 1.00 | 1.00 | 1.00 | 1.00 | 0.99 | 0.95 | 0.87 | 0.58 | 0.15 | 0.01
(0.00) | (0.00) | (0.00) | (0.01) | (0.02) | (0.10) | (0.17) | (0.28) | (0.18) | (0.02)
R | 1.00 | 1.00 | 1.00 | 1.00 | 1.00 | 0.98 | 0.95 | 0.76 | 0.30 | 0.04
(0.00) | (0.00) | (0.00) | (0.00) | (0.01) | (0.05) | (0.10) | (0.25) | (0.27) | (0.07)
3/4 | DS | 1.00 | 1.00 | 1.00 | 1.00 | 1.00 | 1.00 | 0.98 | 0.88 | 0.31 | 0.08
(0.00) | (0.00) | (0.00) | (0.00) | (0.00) | (0.00) | (0.14) | (0.33) | (0.47) | (0.27)
R | 1.00 | 1.00 | 1.00 | 1.00 | 1.00 | 1.00 | 1.00 | 0.90 | 0.44 | 0.06
(0.00) | (0.00) | (0.00) | (0.00) | (0.00) | (0.01) | (0.02) | (0.20) | (0.37) | (0.14)

Figure 1: Power of the knockoff applied to the complete dataset (black), a
data split (red), and a randomised version of the data (blue). The titles
indicate $(f,\rho,p)$. Figure 2: Power of stability selection applied to the
complete dataset (black), a data split (red), and a randomised version of the
data (blue). The titles indicate $(f,\rho,p)$. Table 3: Coverages of
confidence intervals for the selected coefficients; knockoff selection

|  |  |  | $|\beta_{i}|$
---|---|---|---|---
$f$ | $\rho$ | Split | Method | 0 | 0.2 | 0.5 | 1
1/2 | 0 | DS | FV | 68.7 | 90.1 | 90.8 | 89.8
R | FV | 67.3 | 89.9 | 90.6 | 89.8
DS | HD | 89.8 | 91.3 | 90.0 | 90.0
R | HD | 89.6 | 89.9 | 89.9 | 90.3
1/2 | 0.5 | DS | FV | 73.9 | 85.2 | 91.0 | 89.9
R | FV | 72.0 | 84.0 | 91.2 | 89.7
DS | HD | 90.0 | 90.8 | 89.9 | 90.2
R | HD | 89.8 | 89.2 | 90.0 | 90.4
3/4 | 0 | DS | FV | 59.2 | 91.3 | 90.5 | 89.8
R | FV | 57.4 | 91.8 | 90.4 | 89.8
DS | HD | 90.3 | 89.7 | 89.7 | 89.6
R | HD | 89.8 | 90.2 | 90.0 | 89.6
3/4 | 0.5 | DS | FV | 66.4 | 85.8 | 91.5 | 90.0
R | FV | 65.0 | 85.2 | 91.4 | 89.9
DS | HD | 90.6 | 90.1 | 90.4 | 89.8
R | HD | 89.2 | 90.4 | 89.7 | 89.6

Table 4: Coverages of confidence intervals for the selected coefficients;
stability selection

|  |  |  | $|\beta_{i}|$
---|---|---|---|---
$f$ | $\rho$ | Split | Method | 0 | 0.2 | 0.5 | 1
1/2 | 0 | DS | FV | 39.1 | 80.2 | 91.1 | 89.8
R | FV | 34.9 | 82.6 | 91.1 | 89.8
DS | HD | 90.3 | 90.6 | 90.0 | 90.0
R | HD | 90.1 | 89.4 | 89.8 | 89.7
1/2 | 0.5 | DS | FV | 61.7 | 81.6 | 89.5 | 90.0
R | FV | 55.8 | 80.1 | 89.5 | 89.9
DS | HD | 90.6 | 92.5 | 89.9 | 90.1
R | HD | 90.8 | 89.6 | 90.5 | 89.9
3/4 | 0 | DS | FV | 21.9 | 83.7 | 90.7 | 89.8
R | FV | 14.2 | 85.0 | 90.6 | 89.8
DS | HD | 90.8 | 90.2 | 89.6 | 89.6
R | HD | 88.5 | 90.4 | 89.9 | 90.3
3/4 | 0.5 | DS | FV | 46.9 | 79.9 | 90.3 | 89.9
R | FV | 40.0 | 80.8 | 90.6 | 89.9
DS | HD | 88.9 | 89.7 | 90.3 | 89.8
R | HD | 90.8 | 91.0 | 89.6 | 90.3

Table 5: Average length of confidence intervals for the selected coefficients;
knockoff selection

|  |  |  | $|\beta_{i}|$
---|---|---|---|---
$f$ | $\rho$ | Split | Method | 0 | 0.2 | 0.5 | 1
1/2 | 0 | DS | HD | 0.391 | 0.390 | 0.391 | 0.391
R | HD | 0.356 | 0.355 | 0.358 | 0.358
1/2 | 0.5 | DS | HD | 0.508 | 0.507 | 0.510 | 0.482
R | HD | 0.457 | 0.455 | 0.459 | 0.436
3/4 | 0 | DS | HD | 0.653 | 0.651 | 0.650 | 0.651
R | HD | 0.503 | 0.502 | 0.506 | 0.506
3/4 | 0.5 | DS | HD | 0.899 | 0.904 | 0.903 | 0.846
R | HD | 0.644 | 0.645 | 0.650 | 0.617

Table 6: Average length of confidence intervals for the selected coefficients;
stability selection

|  |  |  | $|\beta_{i}|$
---|---|---|---|---
$f$ | $\rho$ | Split | Method | 0 | 0.2 | 0.5 | 1
1/2 | 0 | DS | HD | 0.389 | 0.391 | 0.391 | 0.392
R | HD | 0.360 | 0.354 | 0.358 | 0.358
1/2 | 0.5 | DS | HD | 0.506 | 0.506 | 0.509 | 0.483
R | HD | 0.459 | 0.458 | 0.457 | 0.438
3/4 | 0 | DS | HD | 0.657 | 0.648 | 0.651 | 0.651
R | HD | 0.507 | 0.500 | 0.507 | 0.507
3/4 | 0.5 | DS | HD | 0.890 | 0.904 | 0.902 | 0.847
R | HD | 0.648 | 0.648 | 0.646 | 0.620

Table 7: Coverages of confidence intervals for the coefficients of the
projection parameters; stability selection

|  |  |  | $|\beta_{i}|$
---|---|---|---|---
$f$ | $\rho$ | Split | Method | 0 | 0.2 | 0.5 | 1
1/2 | 0 | DS | FV | 38.6 | 70.9 | 91.2 | 89.4
R | FV | 25.0 | 71.6 | 91.1 | 89.2
DS | HD | 89.7 | 86.1 | 89.4 | 90.0
R | HD | 87.9 | 89.5 | 89.8 | 89.9
1/2 | 0.5 | DS | FV | 41.7 | 75.0 | 83.9 | 90.6
R | FV | 32.5 | 75.3 | 87.5 | 90.8
DS | HD | 91.2 | 88.0 | 88.0 | 90.4
R | HD | 88.4 | 88.9 | 89.2 | 90.2
3/4 | 0 | DS | FV | 11.7 | 71.2 | 90.9 | 89.1
R | FV | 05.5 | 72.8 | 90.3 | 89.1
DS | HD | 88.2 | 89.8 | 90.0 | 90.0
R | HD | 89.4 | 90.3 | 90.0 | 89.9
3/4 | 0.5 | DS | FV | 18.3 | 71.3 | 88.5 | 90.5
R | FV | 12.9 | 70.1 | 89.9 | 90.1
DS | HD | 90.0 | 92.3 | 91.4 | 90.5
R | HD | 89.0 | 91.3 | 90.4 | 90.3

Table 8: Average length of confidence intervals for the coefficients of the
projection parameters; stability selection

|  |  |  | $|s|$
---|---|---|---|---
$f$ | $\rho$ | Split | Method | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8
1/2 | 0 | DS | HD | - | 0.34 | 0.34 | 0.34 | 0.33 | 0.33 | 0.34 | -
R | HD | - | 0.35 | 0.34 | 0.33 | 0.33 | 0.33 | 0.33 | 0.30
1/2 | 0.5 | DS | HD | 0.34 | 0.38 | 0.37 | 0.37 | 0.36 | - | - | -
R | HD | 0.35 | 0.39 | 0.37 | 0.36 | 0.36 | 0.35 | 0.34 | -
3/4 | 0 | DS | HD | - | 0.48 | 0.48 | 0.48 | 0.48 | 0.47 | 0.47 | 0.48
R | HD | - | 0.52 | 0.49 | 0.47 | 0.47 | 0.46 | 0.46 | 0.45
3/4 | 0.5 | DS | HD | 0.48 | 0.55 | 0.54 | 0.52 | 0.52 | 0.53 | - | -
R | HD | 0.49 | 0.56 | 0.54 | 0.52 | 0.51 | 0.48 | 0.49 | -

### Proof of Proposition 1

###### Proof.

Let $\mathcal{R}$ be the support of $R$, write
$p_{r}=\operatorname*{\mathrm{pr}}(R=r)$ for $r\in R$, and let $1(\cdot)$
denote the indicator function. By the harmonic-arithmetic inequality for
Hermitian matrices (Bhagwat and Subramanian, 1978), it follows that

$\displaystyle\frac{1}{f}\mathcal{I}_{Y}(\beta)^{-1}$ $\displaystyle=$
$\displaystyle\left\\{\sum_{i=1}^{n}\operatorname*{\mathrm{pr}}(i\in
R)I_{i}(\beta)\right\\}^{-1}$ (18) $\displaystyle=$
$\displaystyle\left\\{\sum_{i=1}^{n}\sum_{r\in\mathcal{R}}p_{r}1(i\in
r)I_{i}(\beta)\right\\}^{-1}$ (19) $\displaystyle=$
$\displaystyle\left\\{\sum_{r\in\mathcal{R}}p_{r}\sum_{i=1}^{n}1(i\in
r)I_{i}(\beta)\right\\}^{-1}$ (20) $\displaystyle=$
$\displaystyle\left\\{\sum_{r\in\mathcal{R}}p_{r}\mathcal{I}_{r}(\beta)\right\\}^{-1}$
(21) $\displaystyle<$
$\displaystyle\sum_{r\in\mathcal{R}}p_{r}\mathcal{I}_{r}(\beta)^{-1}$ (22)
$\displaystyle=$ $\displaystyle
E\left\\{\mathcal{I}_{R}(\beta)^{-1}\right\\}.$ (23)

The inequality is strict because the $\mathcal{I}_{r}(\beta)$’s are not all
equal. By the monotonicity of $\varphi$, this implies that

$\varphi\left\\{\frac{1}{f}\mathcal{I}_{Y}(\beta)^{-1}\right\\}<\varphi\left[E\left\\{\mathcal{I}_{R}(\beta)^{-1}\right\\}\right].$
(24)

Furthermore, since $\varphi$ is convex, Jensen’s inequality gives

$\varphi\left[E\left\\{\mathcal{I}_{R}(\beta)^{-1}\right\\}\right]\leq
E\left[\varphi\left\\{\mathcal{I}_{R}(\beta)^{-1}\right\\}\right],$ (25)

which shows the first inequality. The second one can be shown with a analogous
argument. ∎

### Proof of Theorem 1

###### Proof.

Assume without loss of generality that the maximum absolute entry of $\eta$ is
$1$ for all $n$. Otherwise we can always divide it by its maximum absolute
entry. Furthermore, assume for the moment that $\hat{\sigma}=\sigma$ and that
$E$ is convex for all $n$.

In the first part of the proof we show that, under these conditions, the
standardised joint distribution of $X^{T}U$ and $\eta^{T}V$ is asymptotically
normal. Define the $2n$ vector and $2n\times(p+1)$ matrix

$T=\begin{pmatrix}\varepsilon\\\ W\end{pmatrix},\quad
A=\begin{pmatrix}X&\eta\\\ \gamma X&-\gamma^{-1}\eta\end{pmatrix}.$ (26)

Note that $A$ has full column rank and that

$A^{T}A=\begin{pmatrix}(1+\gamma^{2})X^{T}X&0_{p}\\\
0_{p}^{T}&(1+\gamma^{-2})\|\eta\|^{2}\end{pmatrix},\quad
A^{T}T=\begin{pmatrix}X^{T}U\\\
\eta^{T}V\end{pmatrix}-\begin{pmatrix}X^{T}\mu\\\ \eta^{T}\mu\end{pmatrix}.$
(27)

Define also the $(p+1)$-dimensional vectors
$P_{i}=\sigma^{-1}(A^{T}A)^{-1/2}A^{T}e_{i}\varepsilon_{i}$ for
$i=1,\ldots,n$, and $P_{i}=\sigma^{-1}(A^{T}A)^{-1/2}A^{T}e_{i}W_{i}$ for
$i=n+1,\ldots,2n$, where $e_{i}$ is the $i$-th canonical vector of
$\mathbb{R}^{2n}$. Clearly, $E(P_{i})=0_{p+1}$ for all $i=1,\ldots,2n$.
Furthermore,

$\sum_{i=1}^{2n}P_{i}=\sigma^{-1}(A^{T}A)^{-1/2}A^{T}T\equiv
Q,\quad\sum_{i=1}^{2n}\operatorname*{\mathrm{var}}(P_{i})=I_{p+1}.$ (28)

Consider events of the form $\\{X^{T}U\in E,\eta^{T}V\leq c\\}$ for
$c\in\mathbb{R}$. These events can be equivalently written as $\\{Q\in
B_{c}\\}$, where

$B_{c}=\left\\{x\in\mathbb{R}^{p+1}\colon\sigma(A^{T}A)^{1/2}x+\begin{pmatrix}X^{T}\mu\\\
\eta^{T}\mu\end{pmatrix}\in E\times(-\infty,c]\right\\}.$ (29)

Clearly, since $E$ is convex, $B_{c}$ is also convex. Thus, by a version of
the multivariate Berry–Esseen theorem for convex sets (Raič, 2019, Theorem
1.1), it follows that

$\displaystyle|\operatorname*{\mathrm{pr}}\\{Q\in
B_{c}\\}-\operatorname*{\mathrm{pr}}\\{N(0_{p+1},I_{p+1})\in B_{c}\\}|$ (30)
$\displaystyle\leq$
$\displaystyle\frac{\alpha}{\sigma^{3}}\left\\{E(|\varepsilon_{1}|^{3})\sum_{i=1}^{n}(e_{i}^{T}A(A^{T}A)^{-1}A^{T}e_{i})^{3/2}+E(|W_{1}|^{3})\sum_{i=n+1}^{2n}e_{i}^{T}A(A^{T}A)^{-1}A^{T}e_{i})^{3/2}\right\\}$
(31) $\displaystyle=$
$\displaystyle\tilde{\alpha}\sum_{i=1}^{2n}(e_{i}^{T}A(A^{T}A)^{-1}A^{T}e_{i})^{3/2},$
(32)

for some $\alpha$ and $\tilde{\alpha}$ are constants independent of $n$. By
the min-max theorem, we can bound the terms of the last sum as follows:

$e_{i}^{T}A(A^{T}A)^{-1}A^{T}e_{i}\leq\|e_{i}^{T}A\|^{2}\lambda_{\max}\left\\{\left(A^{T}A\right)^{-1}\right\\}\leq(p+1)\max(A)^{2}\lambda_{\max}\left\\{\left(A^{T}A\right)^{-1}\right\\}.$
(33)

Hence,

$\sum_{i=1}^{2n}(e_{i}^{T}A(A^{T}A)^{-1}A^{T}e_{i})^{3/2}\leq
2n(p+1)^{3/2}\max(A)^{3}\lambda_{\max}\left\\{\left(A^{T}A\right)^{-1}\right\\}^{3/2}=o(1)$
(34)

by the various asymptotic conditions on $X$ and $\eta$. Recall that we are
assuming that the entries of $\eta$ are bounded, and note that the eigenvalues
of $A^{T}A$ are those of $X^{T}X$ multiplied by $1+\gamma^{2}$, together with
$(1+\gamma^{-2})\|\eta\|^{2}$. The bound given by the theorem is uniform over
convex sets, so we have shown that

$|\operatorname*{\mathrm{pr}}\\{Q\in
B_{c}\\}-\operatorname*{\mathrm{pr}}\\{N(0_{p+1},I_{p+1})\in B_{c}\\}|=o(1)$
(35)

uniformly in $c$. Now, using that $A^{T}A$ is block-diagonal, the asymptotic
probability can be written as

$\displaystyle\operatorname*{\mathrm{pr}}\\{N(0_{p+1},I_{p+1})\in B_{c}\\}$
(36) $\displaystyle=$
$\displaystyle\operatorname*{\mathrm{pr}}\left\\{N\left(\begin{pmatrix}X^{T}\mu\\\
\eta^{T}\mu\end{pmatrix},\sigma^{2}A^{T}A\right)\in
E\times(-\infty,c]\right\\}$ (37) $\displaystyle=$
$\displaystyle\operatorname*{\mathrm{pr}}\left\\{N(X^{T}\mu,(1+\gamma^{2})\sigma^{2}X^{T}X)\in
E\right\\}\operatorname*{\mathrm{pr}}\left\\{N(\eta^{T}\mu,(1+\gamma^{-2})\sigma^{2}\|\eta\|^{2})\leq
c\right\\}.$ (38)

To conclude this part, let $c=kx+\eta^{T}\mu$ for an arbitrary $x$ independent
of $n$ ,where $k=(1+\gamma^{-2})^{1/2}\sigma\|\eta\|$. We have that

$\operatorname*{\mathrm{pr}}\\{k^{-1}(\eta^{T}V-\eta^{T}\mu)\leq x\mid
X^{T}U\in E\\}=\frac{\operatorname*{\mathrm{pr}}(X^{T}U\in E,\eta^{T}V\leq
c)}{\operatorname*{\mathrm{pr}}(X^{T}U\in E,\eta^{T}V<\infty)}.$ (39)

Using the normal approximation to the numerator and denominator, and the fact
that the denominator is bounded away from zero, it follows that

$\operatorname*{\mathrm{pr}}\\{k^{-1}(\eta^{T}V-\eta^{T}\mu)\leq x\mid
X^{T}U\in E\\}\rightarrow\operatorname*{\mathrm{pr}}\\{N(0,1)\leq x\\}.$ (40)

This establishes convergence in distribution under a known variance assumption
and a convex selection set. Extending the result to the conditions of the
theorem is straightforward. First, if $\hat{\sigma}\neq\sigma$, write
$V=\tilde{V}+\gamma^{-1}(\sigma-\hat{\sigma})Z$, with
$\tilde{V}=Y-\gamma^{-1}\sigma Z$. Then,

$k^{-1}(\eta^{T}V-\eta^{T}\mu)=k^{-1}(\eta^{T}\tilde{V}-\eta^{T}\mu)+k^{-1}\gamma^{-1}(\sigma-\hat{\sigma})\eta^{T}Z.$
(41)

We have shown that, conditionally on selection, the first term converges in
distribution to a standard normal. The second one is bounded by a constant
times $(\sigma-\hat{\sigma})\|\eta\|^{-1}\eta^{T}Z$. Note that,
unconditionally on selection, $\sigma-\hat{\sigma}=o_{p}(1)$, by assumption,
and $\|\eta\|^{-1}\eta^{T}Z=O_{p}(1)$, as it has zero mean and unit variance.
Hence, the last term is $o_{p}(1)$ unconditionally on selection. However,
since the selection probability is bounded away from zero, it follows
trivially the last term converges to zero also conditionally on selection, as

$\operatorname*{\mathrm{pr}}\left(|(\sigma-\hat{\sigma})\|\eta\|^{-1}\eta^{T}Z|>\delta\mid
S=s_{n}\right)\leq\frac{\operatorname*{\mathrm{pr}}\left(|(\sigma-\hat{\sigma})\|\eta\|^{-1}\eta^{T}Z|>\delta\right)}{\operatorname*{\mathrm{pr}}\left(S=s_{n}\right)}\rightarrow
0$ (42)

as $n\to\infty$ for any positive $\delta$. The result follows from Slutsky’s
theorem.

Finally, suppose that $E=E_{1}\cup E_{2}$, with $E_{1}$ and $E_{2}$ convex. We
can write

$\displaystyle\operatorname*{\mathrm{pr}}\\{k^{-1}(\eta^{T}V-\eta^{T}\mu)\leq
x\mid X^{T}U\in E\\}$ (43) $\displaystyle=$
$\displaystyle\frac{\operatorname*{\mathrm{pr}}(X^{T}U\in E_{1},\eta^{T}V\leq
c)+\operatorname*{\mathrm{pr}}(X^{T}U\in E_{2},\eta^{T}V\leq
c)-2\operatorname*{\mathrm{pr}}(X^{T}U\in E_{1}\cap E_{2},\eta^{T}V\leq
c)}{\operatorname*{\mathrm{pr}}(X^{T}U\in
E_{1})+\operatorname*{\mathrm{pr}}(X^{T}U\in
E_{2})-2\operatorname*{\mathrm{pr}}(X^{T}U\in E_{1}\cap E_{2})}$ (44)

and use the normal approximation as before, since all the events involved are
convex. When the union involves more than two sets we can proceed recursively.
∎

## References

  * (1)
  * Bachoc et al. (2017) Bachoc, F., Leeb, H. and Pötscher, B. M. (2017), Valid confidence intervals for post-model-selection predictors.  arXiv:1412.4605v3.
  * Bachoc et al. (2020) Bachoc, F., Preinerstorfer, D. and Steinberger, L. (2020), ‘Uniformly valid confidence intervals post-model-selection’, Ann. Stat. 48(1), 440–463.
  * Barber and Candès (2015) Barber, R. F. and Candès, E. (2015), ‘Controlling the false discovery rate via knockoffs’, Ann. Stat. 43(5), 2055–2085.
  * Barber et al. (2020) Barber, R. F., Candès, E., Janson, L., Patterson, E. and Sesia, M. (2020), The knockoff filter for controlled variable selection.  R package version 0.3.3.
https://CRAN.R-project.org/package=knockoff

  * Bayati et al. (2013) Bayati, M., Erdogdu, M. A. and Montanari, A. (2013), Estimating LASSO risk and noise level’, in C. J. C. Burges, L. Bottou, M. Welling, Z. Ghahramani and K. Q. Weinberger, eds, ‘Advances in Neural Information Processing Systems 26 (NIPS 2013)’.
  * Berk et al. (2013) Berk, R., Brown, L., Buja, A., Zhang, K. and Zhao, L. (2013), ‘Valid post-selection inference’, Ann. Stat. 41(2), 802–837.
  * Bhagwat and Subramanian (1978) Bhagwat, K. V. and Subramanian, R. (1978), ‘Inequalities between means of positive operators’, Math. Proc. Camb. Philos. Soc. 83(3), 393–401.
  * Candès et al. (2018) Candès, E., Fan, Y., Janson, L. and Lv, J. (2018), ‘Panning for gold: ‘model‐X’ knockoffs for high dimensional controlled variable selection’, J. R. Statist. Soc. B 80(3), 551–577.
  * Cox (1975) Cox, D. R. (1975), ‘A note on data-splitting for the evaluation of significance levels’, Biometrika 62(2), 441–444.
  * DiCiccio et al. (2020) DiCiccio, C. J., DiCiccio, T. J. and Romano, J. P. (2020), ‘Exact tests via multiple data splitting’, Stat. Probab. Lett. 166.
  * Fan et al. (2012) Fan, J., , Guo, S. and Hao, N. (2012), ‘Variance estimation using refitted cross-validation in ultrahigh dimensional regression’, J. R. Statist. Soc. B 74(1), 37–65.
  * Fithian et al. (2017) Fithian, W., Sun, D. L. and Taylor, J. E. (2017), ‘Optimal inference after model selection’.  arXiv:1410.2597v4.
  * Hofner and Hothorn (2017) Hofner, B. and Hothorn, T. (2017), Stability selection with error control.  R package version 0.6-3.
https://CRAN.R-project.org/package=stabs

  * 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.
  * Ignatiadis et al. (2016) Ignatiadis, N., Klaus, B., Zaugg, J. and Huber, W. (2016), ‘Data-driven hypothesis weighting increases detection power in genome-scale multiple testing’, Nat. Methods 13(7), 577–580.
  * Kivaranovic and Leeb (2020) Kivaranovic, D. and Leeb, H. (2020), ‘On the length of post-model-selection confidence intervals conditional on polyhedral constraints’, JASA 00(0), 1–13.
  * Lee et al. (2016) Lee, J. D., Sun, D. L., Sun, Y. and Taylor, J. E. (2016), ‘Exact post-selection inference, with application to the lasso’, Ann. Stat. 44(3), 907–927.
  * Lee and Taylor (2014) Lee, J. D. and Taylor, J. E. (2014), ‘Exact post model selection inference for marginal screening’, Adv. Neural Inf. Process Syst. 1, 136–144.
  * Lockhart et al. (2014) Lockhart, R., Taylor, J. E., Tibshirani, R. and Tibshirani, R. (2014), ‘A significance test for the lasso’, Ann. Stat. 42(2), 413–468.
  * Loftus and Taylor (2014) Loftus, J. R. and Taylor, J. E. (2014), ‘A significance test for forward stepwise model selection’.  arXiv:1405.3920v1.
  * Markovic et al. (2019) Markovic, J., Taylor, J. and Taylor, J. (2019), Inference after black box selection.  arXiv:1901.09973v1.
  * Meinshausen and Bühlmann (2010) Meinshausen, N. and Bühlmann, P. (2010), ‘Stability selection’, J. R. Statist. Soc. B 72(4), 417–473.
  * Panigrahi et al. (2020) Panigrahi, S., Taylor, J. and Weinstein, A. (2020), ‘Integrative methods for post-selection inference under convex constraints’.  arXiv:1605.08824v7.
  * Raič (2019) Raič, M. (2019), ‘A multivariate berry–esseen theorem with explicit constants’, Bernoulli 25(4A), 2824–2853.
  * Reid et al. (2016) Reid, S., Tibshirani, R. and Friedman, J. (2016), ‘A study of error variance estimation in lasso regression’, Stat. Sin. 26, 35–67.
  * Reitermanová (2010) Reitermanová, Z. (2010), ‘Data splitting’, WDS’10 Proceedings of Contributed Papers 1, 31–36.
  * Rinaldo et al. (2019) Rinaldo, A., Wasserman, L. and G’Sell, M. (2019), ‘Bootstrapping and sample splitting for high-dimensional, assumption-lean inference’, Ann. Stat. 47(6), 3438–3469.
  * Rubin et al. (2006) Rubin, D., Dudoit, S. and van der Laan, M. (2006), ‘A method to increase the power of multiple testing procedures though sample splitting’, Stat. Appl. Genet. Mol. Biol. 5(19).
  * Snee (1977) Snee, R. D. (1977), ‘Validation of regression models: methods and examples’, Technometrics 19, 415–428.
  * Tian (2020) Tian, X. (2020), ‘Prediction error after model search’, Ann. Stat. 48(2), 763–784.
  * Tian and Taylor (2018) Tian, X. and Taylor, J. E. (2018), ‘Selective inference with a randomized response’, Ann. Stat. 46(2), 679–710.
  * Tibshirani et al. (2018) Tibshirani, R., Rinaldo, A., Tibshirani, R. and Wasserman, L. (2018), ‘Uniform asymptotic inference and the bootstrap after model selection’, Ann. Stat. 46(3), 1255–1287.
  * Tibshirani et al. (2019) Tibshirani, R., Tibshirani, R., Taylor, J., Loftus, J., Reid, S. and Markovic, J. (2019), selectiveInference: Tools for post-selection inference.  R package version 1.2.5.
https://CRAN.R-project.org/package=selectiveInference

  * Wasserman and Roeder (2009) Wasserman, L. and Roeder, K. (2009), ‘High-dimensional variable selection’, Ann. Stat. 37(5A), 2178–2201.


# Precise Measurement of the Spin Parameter of the Stellar-Mass Black Hole M33
X-7

Jifeng Liu11affiliation: Harvard-Smithsonian Center for Astrophysics,60 Garden
Street, Cambridge, MA 02138 , Jeffrey E. McClintock11affiliation: Harvard-
Smithsonian Center for Astrophysics,60 Garden Street, Cambridge, MA 02138 ,
Ramesh Narayan11affiliation: Harvard-Smithsonian Center for Astrophysics,60
Garden Street, Cambridge, MA 02138 , Shane W. Davis22affiliation: Institute
for Advanced Study, Einstein Drive, Princeton, NJ 08540 , and Jerome A.
Orosz33affiliation: Department of Astronomy, San Diego State University, 5500
Campanile Drive, San Diego, CA 92182

# Erratum: “Precise Measurement of the Spin Parameter of the Stellar-mass
Black Hole M33 X-7” (ApJL, 679, 37L [2008])

Jifeng Liu11affiliation: Harvard-Smithsonian Center for Astrophysics,60 Garden
Street, Cambridge, MA 02138 , Jeffrey E. McClintock11affiliation: Harvard-
Smithsonian Center for Astrophysics,60 Garden Street, Cambridge, MA 02138 ,
Ramesh Narayan11affiliation: Harvard-Smithsonian Center for Astrophysics,60
Garden Street, Cambridge, MA 02138 , Shane W. Davis22affiliation: Institute
for Advanced Study, Einstein Drive, Princeton, NJ 08540 , and Jerome A.
Orosz33affiliation: Department of Astronomy, San Diego State University, 5500
Campanile Drive, San Diego, CA 92182

###### Abstract

In prior work, Chandra and Gemini-North observations of the eclipsing X-ray
binary M33 X-7 have yielded measurements of the mass of its black hole primary
and the system’s orbital inclination angle of unprecedented accuracy.
Likewise, the distance to the binary is known to a few percent. In an analysis
based on these precise results, fifteen Chandra and XMM-Newton X-ray spectra,
and our fully relativistic accretion disk model, we find that the
dimensionless spin parameter of the black hole primary is $a_{*}=0.77\pm
0.05$. The quoted 1-$\sigma$ error includes all sources of observational
uncertainty. Four Chandra spectra of the highest quality, which were obtained
over a span of several years, all lead to the same estimate of spin to within
statistical errors (2%), and this estimate is confirmed by 11 spectra of lower
quality. There are two remaining uncertainties: (1) the validity of the
relativistic model used to analyze the observations, which is being addressed
in ongoing theoretical work; and (2) our assumption that the black hole spin
is approximately aligned with the angular momentum vector of the binary, which
can be addressed by a future X-ray polarimetry mission.

Galaxies: individual(M33) — X-rays: binaries — black hole physics — binaries:
individual (M33 X-7)

††slugcomment: Published in ApJ Letters: 679, 37-40 (2008)

##  1 INTRODUCTION

M33 X-7 is the first stellar-mass black hole to be discovered that is eclipsed
by its companion [Pietsch et al., 2006]. The X-ray eclipse and the precisely
known distance of this system, $D=840\pm 20$ kpc, underpin the most accurate
dynamical model that has been achieved for any of the 21 known black hole
binaries (Orosz et al. 2007, hereafter O07). The two dynamical parameters of
interest in this Letter are the black hole mass $M=15.65\pm 1.45M_{\odot}$ and
the orbital inclination angle $i=74.6^{\circ}\pm 1.0^{\circ}$ (O07).

Our group has published spin estimates for three stellar-mass black holes
using the X-ray continuum fitting method: GRO J1655-40, $a_{*}=0.65-0.75$; 4U
1543-47, $a_{*}=0.75-0.85$; and GRS 1915+105, $a_{*}=0.98-1.0$ (Shafee et al.
2006, hereafter S06; McClintock et al. 2006, hereafter M06). For LMC X-3,
Davis et al. [2006] find $a_{*}<0.26$. Meanwhile, the Fe line method has been
used to obtain two additional estimates of black hole spin (Brenneman and
Reynolds 2006; Miller et al. 2008). The dimensionless spin parameter
$a_{*}\equiv a/M=cJ/GM^{2}$, where $M$ and $J$ are the mass and angular
momentum of the black hole; $-1\leq a_{*}\leq 1$ [Shapiro and Teukolsky,
1986].

The continuum-fitting method, which was pioneered by [Zhang et al., 1997]
(also see Gierliński et al. 2001), is based on the existence of an innermost
stable circular orbit (ISCO) for a particle orbiting a black hole, inside
which the particle suddenly plunges into the hole. In the continuum-fitting
method, one identifies the inner edge of the black hole’s accretion disk with
the ISCO and estimates the radius $R_{\rm ISCO}$ of this orbit by fitting the
X-ray continuum spectrum. Since the dimensionless radius $r_{\rm isco}\equiv
R_{\rm ISCO}/(GM/c^{2})$ is solely a monotonic function of the black hole spin
parameter [Shapiro and Teukolsky, 1986], knowing its value allows one to
immediately infer the black hole spin parameter $a_{*}$.

Our estimates of spin are based on our fully relativistic accretion disk model
[Li et al., 2005] and an advanced treatment of spectral hardening [Davis et
al., 2005]. We consider only rigorously-selected thermal-state data [Remillard
and McClintock, 2006], which are largely free of the effects of
Comptonization. Furthermore, we only accept data for which the bolometric disk
luminosity is moderate, $L/L_{\rm Edd}<0.3$, in order to ensure that the
standard geometrically-thin thermal disk model is applicable (Shafee et al.
2008; M06).

For the continuum-fitting method to succeed, it is essential to have accurate
values of the black hole mass, orbital inclination and distance (M06),
quantities that are known precisely in the case of M33 X-7. Other virtues of
M33 X-7 for the determination of spin, are the abundance of Chandra and XMM
data, the remarkably thermal and featureless spectrum of the X-ray source, and
its moderate luminosity (§3).

##  2 DATA SELECTION AND REDUCTION

There have been 17 Chandra ACIS observations and 12 XMM-Newton EPIC
observations of M33 X-7. We analyzed Chandra observations (downloaded from the
Chandra Data Archive) with CIAO 3.4 and extracted the spectra from source
ellipses enclosing 95% of the source photons as reported by wavdetect. The XMM
observations were downloaded from the HEASARC archive and analyzed with SAS
7.0.0 [Gabriel et al., 2004]. Using the standard procedures and excluding
intervals of high background, we extracted separate spectra from the PN and
two MOS chips using radii of 400 pixels (i.e., $20^{\prime\prime}$) and fitted
them independently.

The resultant count rate data were folded on the M33 X-7 X-ray eclipse
ephemeris [Pietsch et al., 2006]: HJD $(2453639.119\pm 0.005)\pm N(3.453014\pm
0.000020)$. The folded light curve for all the Chandra data are shown in
Figure $2a$ in O07. For the purpose of measuring the spin of M33 X-7, we
excluded (1) spectra obtained in the phase range –0.3 to 0.2, i.e., during the
eclipse or the pre-eclipse period of erratic X-ray variability (which is
presumably caused by the accretion stream; O07); (2) four ACIS spectra (ObsIDs
6384, 6385, 7170, 7171) that are severely affected by pile-up; and (3) all
spectra that contain less than 1000 counts. The 15 spectra so selected are
listed in Table 1 and comprise eight Chandra ACIS spectra and seven EPIC
PN/MOS spectra from five XMM observations. We refer throughout to the four
ACIS spectra with $\buildrel{\scriptstyle>}\over{\scriptstyle\sim}$ 5,000
counts as the “gold” spectra and to the rest as the “silver” spectra.

##  3 ANALYSIS AND RESULTS

The procedures used here are precisely the same as those that are described
fully in M06. Briefly, the relativistic accretion disk model kerrbb2 has just
two fit parameters, namely the black hole spin $a_{*}$ and the mass accretion
rate $\dot{M}$ (or equivalently, $a_{*}$ and the Eddington-scaled bolometric
luminosity, $l\equiv L_{\rm bol}(a_{*},\dot{M})/L_{\rm Edd}$; M06). In the
case of M33 X-7, we also fit for a third parameter, $N_{\rm H}$, the hydrogen
column density (phabs in XSPEC).

The spectral hardening factor $f\equiv T_{\rm col}/T_{\rm eff}$ was computed
as a function of $l$ for the appropriate metallicity of M33 X-7
($Z=0.1Z_{\odot}$; O07) using the model of Davis and Hubeny 2006 (bhspec in
XSPEC). These values of $f$ are contained in a pair of lookup tables, which
correspond to two representative values of the viscosity parameter
($\alpha=0.01,0.1$; M06) for a wide range of the spin parameter (e.g.,
$0<a_{*}<0.99$). We find that our results are quite insensitive to the choice
of $\alpha$ or an increase in metallicity. We have also experimented with
varying the input parameters $M$, $D$, $i$ and $N_{\rm H}$, and we find that
the values of $f$ are scarcely affected.

All of the spectra were well-fitted using a simple absorbed kerrbb2 model
[i.e., phabs(kerrbb2) in XSPEC]. Notably, neither Fe line/edge components nor
an additional nonthermal component was required, as they were in our earlier
work (R06, M06). We fitted each spectrum for $a_{*}$, the mass accretion rate
$\dot{M}$, and the neutral hydrogen column density $N_{\rm H}$ with the input
parameters fixed at their baseline values (see §1; O07). The normalization was
fixed at unity (as appropriate when $M$, $i$ and $D$ are held fixed). We
included the effects of limb darkening (lflag = 1) and returning radiation
effects (rflag = 1), and we set the torque at the inner boundary of the
accretion disk to zero ($\eta=0$).

The fits obtained for all 15 spectra are quite acceptable with
$\chi^{2}_{\nu}<1.2$; results for fits over the energy range 0.3–8 keV are
summarized in Table 1. An inspection of the fitting residuals for all the
spectra show them to be free of any systematic effects, as illustrated in
Figure 1. Given the modest luminosities, $0.07<l<0.11$ (Table 1), which are
well below our selection limit of $l=0.3$, the accretion disk in M33 X-7 is
quite thin, $H/R\leq 0.04$ (see Fig. 17 in Shafee et al. 2008), and our
assumption of zero torque at the inner boundary is likely to be valid.

Figure 2 shows plots of $a_{*}$ for all 15 observations, which are ordered by
the number of counts detected. Each of the four panels corresponds to a
different choice for the energy interval used in fitting the data (e.g., 0.3–8
keV, 0.5–8 keV, etc.); a comparison of the results in the four panels shows
that this choice is quite unimportant. The four gold spectra with $\lower
2.0pt\hbox{$\buildrel{\scriptstyle>}\over{\scriptstyle\sim}$}\ 5,000$ counts
each (solid symbols) yield spin estimates that agree with their mean value
(indicated by the dotted lines) typically to within their $\approx 2$%
statistical uncertainties. The stability of these four gold spectra is
especially remarkable given that three of the observations were separated by
3-month intervals in 2005–2006, and one of them was obtained five years
earlier in 2000 (Table 1). The dispersion for the 11 silver spectra that have
$\lower 2.0pt\hbox{$\buildrel{\scriptstyle<}\over{\scriptstyle\sim}$}\ 3,000$
counts (open symbols) is much larger. However, in each panel, the mean of
these 11 spin values agrees with the mean determined using the gold spectra to
within $\approx 1$%. As concluded in the caption of Figure 2, our adopted
average spin for the four gold spectra is $\bar{a}_{*}=0.77$ with a standard
deviation of $\Delta a_{*}=0.02$.

In order to determine the error in $a_{*}$ due to the combined uncertainties
in $M$, $i$ and $D$ (§1), we performed Monte Carlo simulations assuming that
the uncertainties in these parameters are normally and independently
distributed. The results for 3,000 simulation runs are plotted in Figure 3.
The histogram of $\bar{a}_{*}-\bar{a}_{*0}$ shows that the $1\sigma$ error in
the spin due to the combined uncertainties of the three input parameters is
about $\Delta a_{*}=0.05$. The error is dominated by the uncertainty in $M$;
the uncertainties in $i$ and $D$ are relatively unimportant. This error is
based on a readily available table that was computed for solar metallicity.
Despite this limitation, we believe that our error estimate is accurate
because the effects of going from $Z=0.1Z_{\odot}$ to $Z=Z_{\odot}$ are very
small at the luminosities in question, $l\approx 0.1$.

##  4 DISCUSSION

The largest error in our spin estimate arises from the uncertainties in the
validity of the disk model we employ. For example, the spin depends on
accurate model determinations of the hardening factor $f$; this problem is
quite tractable and vigorous theoretical efforts are underway (Davis et al.
2005, 2006, Blaes et al. 2006). Possibly more problematic is our assumption
that the viscous torque vanishes at the ISCO and that there is no significant
emission from the gas inside the ISCO.

Hydrodynamic models of the accretion disk indicate that the viscous torque at
the ISCO as well as emission from inside the ISCO should both be negligible
for the geometrically thin disks and low luminosities ($l\leq 0.3$) that we
restrict ourselves to (Afshordi and Paczyński 2003; S08). The emission from
inside the ISCO causes rather modest errors in spin estimates; in the case of
M33 X-7, the estimated error is $\Delta a_{*}\leq 0.01$ since $l\leq 0.1$ and
hence $R/H\leq 0.04$ (M06; S08). On the other hand, MHD simulations of
accretion flows around black holes [Hawley and Krolik, 2002, Beckwith et al.,
2008] find a large torque at the ISCO and substantial dissipation inside the
ISCO. We note, however, that these simulations carried out so far are for
geometrically thick systems, with $H/R\sim 0.2$; these flows are nearly an
order of magnitude thicker than the disk in M33 X-7. In the hydrodynamic
models of Shafee et al. [2008] the stress at the ISCO increases rapidly with
increasing disk thickness, so it is conceivable that there is no serious
disagreement between the hydrodynamic and MHD results. Numerical MHD
simulations of truly thin disks are necessary to resolve this issue. We note
that a recent MHD simulation of a geometrically thin accretion disk for a
pseudo-Newtonian potential does show a dramatic drop in the mid-plane density
and vertical column density over a narrow range of radii close to the ISCO
[Reynolds and Fabian, 2008].

Although there is theoretical uncertainty about conditions near the ISCO,
there is a long history of evidence suggesting that fitting the X-ray
continuum is a promising approach to measuring black hole spin. This history
begins in the mid-1980s with the simple non-relativistic multicolor disk model
[Mitsuda et al., 1984], which returns the color temperature $T_{\rm in}$ at
the inner-disk radius $R_{\rm in}$. Tanaka and Lewin [1995] summarize examples
of the steady decay (by factors of 10–100) of the thermal flux of transient
sources during which $R_{\rm in}$ remains quite constant (see their Fig.
3.14). More recently, this evidence for a constant inner radius in the thermal
state has been presented for a number of sources in several papers via plots
showing that the bolometric luminosity of the thermal component is
approximately proportional to $T^{4}$ (McClintock et al. 2007, and references
therein). Obviously, this non-relativistic analysis cannot provide a secure
value for the radius of the ISCO nor even establish that this stable radius is
the ISCO. Nevertheless, the presence of a fixed radius indicates that the
continuum-fitting method is a well-founded approach to measuring black hole
spin.

It is reasonable to assume that the inner X-ray-emitting portion of the disk
is aligned with the spin axis of the black hole by the Bardeen-Petterson
effect [lod05]. Throughout, in making use of the orbital inclination angle, we
have assumed that the black hole spin is aligned with the angular momentum
vector of the binary system. As Figure 3 indicates, if any misalignment is
$\lower 2.0pt\hbox{$\buildrel{\scriptstyle<}\over{\scriptstyle\sim}$}\
3^{\circ}$, then it will contribute an error in $a_{*}$ that is no larger than
our total observational error of $\Delta a_{*}=0.05$. There is no evidence for
significant misalignments despite the often-cited examples of GRO J1655-40 and
SAX J1819.3-2525 (see §2.2 in Narayan and McClintock 2005; but see Maccarone
2002). The clear-cut way to assess the degree of alignment is via X-ray
polarimetric observations of black hole systems in the thermal state (Li et
al. 2008, in preparation).

What is the origin of the spin of M33 X-7? Was the black hole born with its
present spin, or was it torqued up gradually via the accretion flow supplied
by its companion? In order to achieve a spin of $a_{*}=0.77$ via disk
accretion, an initially non-spinning black hole must accrete $4.9M_{\odot}$
from its donor [King and Kolb, 1999] in becoming the $M=15.65M_{\odot}$ that
we observe today (O07). However, to transfer this much mass even in the case
of Eddington-limited accretion ($\dot{M}_{\rm Edd}\equiv L_{\rm
Edd}/c^{2}\approx 4\times 10^{-8}M_{\odot}/{\rm yr}$) requires $\sim 120$
million years, whereas the age of the system is only 2–3 million years (O07).
Thus, it appears that the spin of M33 X-7 must be natal, which is the same
conclusion that has been reached for two other stellar black holes (S06, M06;
but see Bethe et al. 2003 on the possibility of hypercritical accretion)

M33 X-7’s secure dynamical data and distance, the X-ray source’s clean
thermal-state spectrum and moderate luminosity, and an abundance of Chandra
and XMM data have provided arguably the most secure estimate of black hole
spin that has been achieved to date: $a_{*}=0.77\pm 0.05$, where the error
estimate includes all sources of observational error. Since an astrophysical
black hole can be described by just the two parameters that specify its mass
and spin [Shapiro and Teukolsky, 1986], we now have a complete description of
an asteroid-size object that is situated at a distance of about one Mpc.

JFL and SWD acknowledge support from NASA through the Chandra Fellowship
Program, grants PF6-70043 and PF6-70045. JEM acknowledges support from NASA
grant AR8-9006X. We thank Rebecca Shafee for technical advice and Jack Steiner
for critical comments on the manuscript. In the paper “Precise Measurement of
the Spin Parameter of the Stellar-Mass Black Hole M33 X-7” by Jifeng Liu,
Jeffrey E. McClintock, Ramesh Narayan, Shane W. Davis, and Jerome A. Orosz
(ApJ, 679, L37 [2008]), the reported value of the black-hole spin parameter
$a_{*}=0.77\pm 0.05$ is in error. The correct value is larger by 0.068 and is
$a_{*}=0.84\pm 0.05$. The error is the result of a bug in the XSPEC accretion-
disk model
kerrbb.111http://heasarc.gsfc.nasa.gov/docs/xanadu/xspec/issues/archive/issues.12.5.0an.html
(patch 12.5.0a). Prior to 1 December 2008, the model’s two parameter flags
that switch limb darkening and self-irradiation of the disk on/off were
reversed (e.g., “par8” incorrectly controlled limb darkening rather than self-
irradiation). In computing tables of the spectral hardening factor $f$, we use
both kerrbb and the disk atmosphere model bhspec [McClintock et al., 2006].
Because the latter model does not include the effect of self-irradiation, we
switch this feature off in kerrbb when computing the $f$-tables. In this
instance, because of the bug we switched off limb darkening instead of self-
irradiation, which corrupted our results. Meanwhile, our earlier spin results
for GRS 1915+105 (McClintock et al. 2006) and for 4U 1543–47 and GRO J1655–40
[Shafee et al., 2006] are unaffected by the bug. The figures and tabular data
in the paper are essentially unaffected, apart from the increase in $a_{*}$
and corresponding decreases in $f$ and the Eddington-scaled luminosity $l$
(8.7% and 4.5%, respectively, for the four gold spectra). The higher spin
increases somewhat our estimate of how much mass ($4.9~{}M_{\odot}$) and time
($\sim 120$ million years) would be required to spin up an initially
nonspinning black hole to the present spin of M33 X-7. In order to achieve
$a_{*}=0.84$, the black hole must accrete $5.7~{}M_{\odot}$, which would
require $\sim 140$ million years. Because the age of the binary system is only
2–3 million years this change does not at all affect our conclusion that the
spin of M33 X-7 is natal.

## References

  * Afshordi and Paczyński [2003] N. Afshordi and B. Paczyński.  Geometrically Thin Disk Accreting into a Black Hole.  _ApJ_ , 592:354–367, July 2003.  10.1086/375559.
  * Beckwith et al. [2008] K. Beckwith, J. F. Hawley, and J. H. Krolik.  Where is the radiation edge in magnetized black hole accretion discs?  _MNRAS_ , 390:21–38, October 2008.  10.1111/j.1365-2966.2008.13710.x.
  * Bethe et al. [2003] H. A. Bethe, G. E. Brown, and C.-H. Lee.  _Formation and evolution of black holes in the Galaxy : selected papers with commentary_.  2003\.
  * Blaes et al. [2006] O. M. Blaes, S. W. Davis, S. Hirose, J. H. Krolik, and J. M. Stone.  Magnetic Pressure Support and Accretion Disk Spectra.  _ApJ_ , 645:1402–1407, July 2006.  10.1086/503741.
  * Brenneman and Reynolds [2006] L. W. Brenneman and C. S. Reynolds.  Constraining Black Hole Spin via X-Ray Spectroscopy.  _ApJ_ , 652:1028–1043, December 2006.  10.1086/508146.
  * Davis and Hubeny [2006] S. W. Davis and I. Hubeny.  A Grid of Relativistic, Non-LTE Accretion Disk Models for Spectral Fitting of Black Hole Binaries.  _ApJS_ , 164:530–535, June 2006.  10.1086/503549.
  * Davis et al. [2005] S. W. Davis, O. M. Blaes, I. Hubeny, and N. J. Turner.  Relativistic Accretion Disk Models of High-State Black Hole X-Ray Binary Spectra.  _ApJ_ , 621:372–387, March 2005.  10.1086/427278.
  * Davis et al. [2006] S. W. Davis, C. Done, and O. M. Blaes.  Testing Accretion Disk Theory in Black Hole X-Ray Binaries.  _ApJ_ , 647:525–538, August 2006.  10.1086/505386.
  * Gabriel et al. [2004] C. Gabriel, M. Denby, D. J. Fyfe, J. Hoar, A. Ibarra, E. Ojero, J. Osborne, R. D. Saxton, U. Lammers, and G. Vacanti.  The XMM-Newton SAS - Distributed Development and Maintenance of a Large Science Analysis System: A Critical Analysis.  In F. Ochsenbein, M. G. Allen, & D. Egret, editor, _Astronomical Data Analysis Software and Systems (ADASS) XIII_ , volume 314 of _Astronomical Society of the Pacific Conference Series_ , pages 759–+, July 2004.
  * Gierliński et al. [2001] M. Gierliński, A. Maciołek-Niedźwiecki, and K. Ebisawa.  Application of a relativistic accretion disc model to X-ray spectra of LMC X-1 and GRO J1655-40.  _MNRAS_ , 325:1253–1265, August 2001.  10.1046/j.1365-8711.2001.04540.x.
  * Hawley and Krolik [2002] J. F. Hawley and J. H. Krolik.  High-Resolution Simulations of the Plunging Region in a Pseudo-Newtonian Potential: Dependence on Numerical Resolution and Field Topology.  _ApJ_ , 566:164–180, February 2002.  10.1086/338059.
  * King and Kolb [1999] A. R. King and U. Kolb.  The evolution of black hole mass and angular momentum.  _MNRAS_ , 305:654–660, May 1999.  10.1046/j.1365-8711.1999.02482.x.
  * Li et al. [2005] L.-X. Li, E. R. Zimmerman, R. Narayan, and J. E. McClintock.  Multitemperature Blackbody Spectrum of a Thin Accretion Disk around a Kerr Black Hole: Model Computations and Comparison with Observations.  _ApJS_ , 157:335–370, April 2005.  10.1086/428089.
  * Maccarone [2002] T. J. Maccarone.  On the misalignment of jets in microquasars.  _MNRAS_ , 336:1371–1376, November 2002.  10.1046/j.1365-8711.2002.05876.x.
  * McClintock et al. [2006] J. E. McClintock, R. Shafee, R. Narayan, R. A. Remillard, S. W. Davis, and L.-X. Li.  The Spin of the Near-Extreme Kerr Black Hole GRS 1915+105.  _ApJ_ , 652:518–539, November 2006.  10.1086/508457.
  * McClintock et al. [2007] J. E. McClintock, R. Narayan, and R. Shafee.  Estimating the Spins of Stellar-Mass Black Holes.  _ArXiv e-prints_ , July 2007.
  * Miller et al. [2008] J. M. Miller, C. S. Reynolds, A. C. Fabian, E. M. Cackett, G. Miniutti, J. Raymond, D. Steeghs, R. Reis, and J. Homan.  Initial Measurements of Black Hole Spin in GX 339-4 from Suzaku Spectroscopy.  _ApJ_ , 679:L113–L116, June 2008.  10.1086/589446.
  * Mitsuda et al. [1984] K. Mitsuda, H. Inoue, K. Koyama, K. Makishima, M. Matsuoka, Y. Ogawara, K. Suzuki, Y. Tanaka, N. Shibazaki, and T. Hirano.  Energy spectra of low-mass binary X-ray sources observed from TENMA.  _PASJ_ , 36:741–759, 1984.
  * Narayan and McClintock [2005] R. Narayan and J. E. McClintock.  Inclination Effects and Beaming in Black Hole X-Ray Binaries.  _ApJ_ , 623:1017–1025, April 2005.  10.1086/428709.
  * Orosz et al. [2007] J. A. Orosz, J. E. McClintock, R. Narayan, C. D. Bailyn, J. D. Hartman, L. Macri, J. Liu, W. Pietsch, R. A. Remillard, A. Shporer, and T. Mazeh.  A 15.65-solar-mass black hole in an eclipsing binary in the nearby spiral galaxy M 33.  _Nature_ , 449:872–875, October 2007.  10.1038/nature06218.
  * Pietsch et al. [2006] W. Pietsch, F. Haberl, M. Sasaki, T. J. Gaetz, P. P. Plucinsky, P. Ghavamian, K. S. Long, and T. G. Pannuti.  M33 X-7: ChASeM33 Reveals the First Eclipsing Black Hole X-Ray Binary.  _ApJ_ , 646:420–428, July 2006.  10.1086/504704.
  * Remillard and McClintock [2006] R. A. Remillard and J. E. McClintock.  X-Ray Properties of Black-Hole Binaries.  _ARA &A_, 44:49–92, September 2006.  10.1146/annurev.astro.44.051905.092532.
  * Reynolds and Fabian [2008] C. S. Reynolds and A. C. Fabian.  Broad Iron-K$\alpha$ Emission Lines as a Diagnostic of Black Hole Spin.  _ApJ_ , 675:1048–1056, March 2008.  10.1086/527344.
  * Shafee et al. [2006] R. Shafee, J. E. McClintock, R. Narayan, S. W. Davis, L.-X. Li, and R. A. Remillard.  Estimating the Spin of Stellar-Mass Black Holes by Spectral Fitting of the X-Ray Continuum.  _ApJ_ , 636:L113–L116, January 2006.  10.1086/498938.
  * Shafee et al. [2008] R. Shafee, R. Narayan, and J. E. McClintock.  Viscous Torque and Dissipation in the Inner Regions of a Thin Accretion Disk: Implications for Measuring Black Hole Spin.  _ApJ_ , 676:549–561, March 2008.  10.1086/527346.
  * Shapiro and Teukolsky [1986] S. L. Shapiro and S. A. Teukolsky.  _Black Holes, White Dwarfs and Neutron Stars: The Physics of Compact Objects_.  June 1986.
  * Tanaka and Lewin [1995] Y. Tanaka and W. H. G. Lewin.  Black hole binaries.  In W. H. G. Lewin, J. van Paradijs, & E. P. J. van den Heuvel, editor, _X-ray binaries, p. 126 - 174_ , pages 126–174, 1995.
  * Zhang et al. [1997] S. N. Zhang, W. Cui, and W. Chen.  Black Hole Spin in X-Ray Binaries: Observational Consequences.  _ApJ_ , 482:L155+, June 1997.  10.1086/310705.

Table 1: Kerrbb2 fit results for M33 X-7 in 0.3-8 keVaaThe columns are (1) ID number; (2) date of observation; (3) exposure time in ksec; (4) no. of counts; (5) spin parameter; (6) mass accretion rate in $10^{18}$ g s-1; (7) hydrogen column density in $10^{20}$ cm-2; (8) spectral hardening factor; (9) Eddington-scaled bolometric luminosity; and (10) reduced $\chi^{2}$ per dof.  spectrum | obs-date | Texp | counts | $a_{*}$ | $\dot{M}$ | $n_{\rm H}$ | $f_{\rm col}$ | $\lg l$ |  $\chi^{2}_{\nu}$/dof
---|---|---|---|---|---|---|---|---|---
acis6376 | 2006-03-03 | 93.1 | 9748 | 0.751 $\pm$ 0.026 | 1.88 $\pm$ 0.12 | 11.1 $\pm$ 1.1 | 1.78 | -1.01 | 1.07/180
acis6387 | 2006-06-26 | 77.3 | 7271 | 0.782 $\pm$ 0.019 | 1.64 $\pm$ 0.10 | 11.4 $\pm$ 1.2 | 1.76 | -1.05 | 0.93/157
acis6382 | 2005-11-23 | 72.3 | 6515 | 0.772 $\pm$ 0.030 | 1.72 $\pm$ 0.14 | 9.9 $\pm$ 1.4 | 1.78 | -1.04 | 1.17/152
acis1730 | 2000-07-12 | 49.5 | 4855 | 0.800 $\pm$ 0.026 | 1.37 $\pm$ 0.11 | 6.1 $\pm$ 1.3 | 1.77 | -1.12 | 1.15/126
acis7344 | 2006-07-01 | 21.5 | 1711 | 0.873 $\pm$ 0.031 | 1.05 $\pm$ 0.14 | 8.6 $\pm$ 2.8 | 1.75 | -1.16 | 0.69/55
acis6386 | 2005-10-31 | 14.9 | 1491 | 0.786 $\pm$ 0.041 | 1.55 $\pm$ 0.21 | 12.8 $\pm$ 3.2 | 1.77 | -1.07 | 0.95/49
acis7197 | 2005-11-03 | 12.7 | 1117 | 0.892 $\pm$ 0.043 | 0.97 $\pm$ 0.20 | 9.1 $\pm$ 4.1 | 1.75 | -1.17 | 0.99/37
acis7208 | 2005-11-21 | 11.5 | 1014 | 0.678 $\pm$ 0.110 | 1.73 $\pm$ 0.39 | 13.9 $\pm$ 4.6 | 1.78 | -1.10 | 0.81/33
PN0102642301 | 2002-01-27 | 10.0 | 2724 | 0.832 $\pm$ 0.031 | 1.34 $\pm$ 0.14 | 7.3 $\pm$ 1.0 | 1.75 | -1.09 | 0.84/103
PN0102641201 | 2000-08-02 | 10.3 | 1836 | 0.618 $\pm$ 0.056 | 2.51 $\pm$ 0.29 | 12.0 $\pm$ 1.5 | 1.78 | -0.97 | 0.87/69
PN0141980801 | 2003-02-12 | 8.4 | 1596 | 0.636 $\pm$ 0.074 | 2.50 $\pm$ 0.36 | 10.7 $\pm$ 1.5 | 1.78 | -0.96 | 0.90/60
PN0141980601 | 2003-01-23 | 11.6 | 1545 | 0.656 $\pm$ 0.077 | 2.21 $\pm$ 0.32 | 12.5 $\pm$ 1.6 | 1.77 | -1.00 | 1.00/59
M10102642301 | 2002-01-27 | 12.3 | 1199 | 0.841 $\pm$ 0.042 | 1.30 $\pm$ 0.20 | 6.7 $\pm$ 2.1 | 1.75 | -1.10 | 0.89/40
PN0102640401 | 2000-08-02 | 9.2 | 1136 | 0.838 $\pm$ 0.039 | 1.42 $\pm$ 0.20 | 6.8 $\pm$ 1.8 | 1.75 | -1.06 | 1.15/44
M20102642301 | 2002-01-27 | 12.3 | 1135 | 0.839 $\pm$ 0.043 | 1.25 $\pm$ 0.20 | 7.2 $\pm$ 2.2 | 1.75 | -1.12 | 0.83/39

Figure 1: X-ray spectrum of M33 X-7. (upper panel) This spectrum (ObsID 6376)
is representative of the four gold spectra (see text). The histogram shows a
model that has been fitted to the spectrum (0.3–8.0 keV), which is comprised
of only the thermal disk component (kerrbb2) and a low-energy absorption model
(phabs). The fit parameters are summarized in Table 1. (lower panel) The fit
is good ($\chi_{\nu}^{2}/dof=1.074/180$) and the fit residuals show no
systematic structure; in particular, there is no evidence for an Fe-line,
absorption edges, or a nonthermal power-law/Comptonization component of
emission at higher energies. Figure 2: Spin results for all 15 spectra ordered
by total 0.3–8 keV counts. (a) Results based on spectral fits over the energy
interval 0.3–8 keV. Filled circles are for the gold Chandra spectra, and the
other plotting symbols are for the 11 silver spectra, which include both
Chandra ACIS spectra (open circles) and XMM-Newton EPIC spectra (crosses). The
indicated uncertainties are at the 90% level of confidence. The dotted line
indicates the average spin for the 4 gold spectra, $\bar{a}_{*}=0.776\pm
0.018$. The average for the 11 silver spectra is almost identical, although
the dispersion is much greater, $\bar{a}_{*}=0.772\pm 0.098$. (b–d) Same as
panel a except the fit interval is as indicated rather than 0.3–8 keV. Figure
3: Effect on the spin parameter $a_{*}$ of varying the input parameters $M$,
$i$ and $D$. (a) Spin versus mass $M$ for 3,000 sets of parameters drawn at
random. The black filled circle indicate our final adopted estimate of the
spin ($\bar{a}_{*0}=0.77\pm 0.02$) for $M$, $i$ and $D$ at their baseline
values. (b) Spin versus distance $D$. (c) Spin versus inclination angle $i$.
(d) Histogram of spin displacements for 3,000 parameters sets. The vertical
solid line indicates the average spin ($\bar{a}_{*}=0.77$). The two dotted
lines enclose 68.3% of the spin values centered on the solid line; the half-
separation, $\Delta(\bar{a}_{*}-\bar{a}_{*0})=0.053$, represents the $1\sigma$
error in the average spin of $\bar{a}_{*0}=0.77$.