# Selective Inference in Propensity Score Analysis

Yoshiyuki Ninomiya
Department of Statistical Inference and Mathematics
The Institute of Statistical Mathematics  Yuta Umezu
School of Information and Data Science, Nagasaki University  Ichiro Takeuchi
Department of Computer Science, Nagoya Institute of Technology
RIKEN Center for Advanced Intelligence Project

###### Abstract

Selective inference (post-selection inference) is a methodology that has
attracted much attention in recent years in the fields of statistics and
machine learning. Naive inference based on data that are also used for model
selection tends to show an overestimation, and so the selective inference
conditions the event that the model was selected. In this paper, we develop
selective inference in propensity score analysis with a semiparametric
approach, which has become a standard tool in causal inference. Specifically,
for the most basic causal inference model in which the causal effect can be
written as a linear sum of confounding variables, we conduct Lasso-type
variable selection by adding an $\ell_{1}$ penalty term to the loss function
that gives a semiparametric estimator. Confidence intervals are then given for
the coefficients of the selected confounding variables, conditional on the
event of variable selection, with asymptotic guarantees. An important property
of this method is that it does not require modeling of nonparametric
regression functions for the outcome variables, as is usually the case with
semiparametric propensity score analysis.

Keywords: causal inference; confidence interval; inverse-probability-weighted
estimation; Lasso; pivot statistic; post-selection inference; variable
selection

##  1 Introduction

It has often been pointed out that generating a hypothesis of interest from
data and then using the same data to construct a statistical test or a
confidence interval for that hypothesis is inappropriate. Breiman (2001)
called it a “quiet scandal of statistics.” The so-called “file drawer problem”
can be regarded as this kind of problem. A typical example in regression is
naive inference on the coefficients of the selected variables after variable
selection. To deal with this variable selection problem, “simultaneous
inference” was proposed by Berk et al. (2013). Simultaneous inference gives a
reasonable inference no matter what model is selected, but as one can imagine
from the methods of multiple testing, it is rather conservative. In these
situations, a “selective inference” approach, which is based on the work of
Lee and Taylor (2014), Lee et al. (2016), and Tibshirani et al. (2016), has
recently gained attention in the field of statistics and machine learning.
Selective inference conditions on that the model has been selected. At about
the same time as these papers were published, Taylor and Tibshirani (2015)
explained the great potential of selective inference. It has been extensively
developed since then.

On the other hand, causal inference has been a hot topic in biostatistics and
econometrics and has recently been of importance in machine learning (see,
e.g., Peters et al. 2017; Hernán and Robins 2020). In particular, propensity
score analysis with a semiparametric approach has become a standard tool, as
it is a theoretically guaranteed method that avoids the bias introduced by the
existence of confounding. To illustrate the motivation for this paper, let us
examine one of the simplest models in causal inference, where the causal
effect can be written as a linear sum of confounding variables, as follows:

$\displaystyle
Y_{i}=T_{i}\bm{X}_{i}^{\prime}\bm{\beta}^{(1)}+(1-T_{i})\bm{X}_{i}^{\prime}\bm{\beta}^{(0)}+\xi_{i}\qquad(1\leq
i\leq n).$

The subscript $i\ (\in\\{1,\ldots,n\\})$ represents the $i$-th sample, $Y_{i}\
(\in\mathbb{R})$ is the outcome variable, $T_{i}\ (\in\\{0,1\\})$ is the
assignment variable, and $\bm{X}_{i}=(X_{1i},\ldots,X_{pi})^{\prime}$ is the
confounding variable vector for $Y_{i}$ and $T_{i}$. In addition, $\xi_{i}\
(\in\mathbb{R})$ is an unobserved random variable that depends on
$\bm{X}_{i}$, and the conditional distribution given $\bm{X}_{i}$ is supposed
to be normal. We assume that
$T_{i}\mathop{\perp\\!\\!\\!\\!\perp}Y_{i}\mid\bm{X}_{i}$, which is called the
ignorable treatment assignment condition. In this model, the parameter of the
causal effect is $\bm{\beta}^{(1)}-\bm{\beta}^{(0)}$, and the basic
semiparametric method for consistent estimation is the inverse-probability-
weighted estimation using the propensity score $e(\bm{X}_{i})\equiv{\rm
P}(T_{i}=1\mid\bm{X}_{i})$ (Rubin 1985; Robins et al. 1994). It has the
advantage of not requiring an estimation of ${\rm E}(\xi_{i}\mid\bm{X}_{i})$,
which is difficult to model, and we rely upon it in this paper.

Now, supposing that the dimension $p$ of the confounding variable vector is
high, let us try to estimate the causal effect by selecting variables to
improve the estimation accuracy. Specifically, letting $\|\cdot\|_{1}$ be the
$\ell_{1}$ norm, we consider

$\displaystyle\sum_{i=1}^{n}\left[\left\\{\frac{T_{i}}{e(\bm{X}_{i})}-\frac{1-T_{i}}{1-e(\bm{X}_{i})}\right\\}Y_{i}-\bm{X}_{i}^{\prime}\bm{\beta}\right]^{2}+\lambda\|\bm{\beta}\|_{1}$

and give its minimizer with respect to $\bm{\beta}$ as a Lasso-type estimator
(Tibshirani 1996). As mentioned at the beginning, this type of variable
selection affects the subsequent inference of the causal effect, for example,
in constructing a confidence interval for $\bm{\beta}^{(1)}-\bm{\beta}^{(0)}$.
Therefore, we will attempt to combine causal inference with selective
inference. Specifically, on the basis of the method of Lee et al. (2016), we
can give a confidence interval with a confidence coefficient of $\alpha$
conditional on the event of variable selection for each component of
$\bm{\beta}^{(1)}-\bm{\beta}^{(0)}$ with the guarantee of asymptotics. The
same problem has been addressed by Zhao et al. (2017), who give a
nonparametric estimation of ${\rm E}(\xi_{i}\mid\bm{X}_{i})$ without using
propensity scores. One of the reasons why propensity score analysis has become
a hot topic is that it does not necessarily require this nonparametric
estimation; we put much value on it in this paper.

As a preparation, Section 2 introduces propensity score analysis in causal
inference models and selective inference with a focus on Lee et al. (2016).
Section 3 deals with selective inference in propensity score analysis and
builds a foundation for selective inference by conditioning on the assignment
variables before the outcome variables, an operation that is not usually
considered in propensity score analysis. Then, asymptotically guaranteed
confidence intervals are given by using higher-order asymptotic theory.
Sections 4 and 5 describe numerical experiments and real data analyses
comparing the proposed method with a naive method that ignores the influence
of model selection. We find that there is a significant difference between the
two methods and that the naive method lacks validity to a large extent. In
Section 6, we extend the proposed method to handle especially developed
versions of propensity score analysis and conclude in Section 7.

##  2 Preparation

###  2.1 Causal inference model

Let us assume that there are $H$ types of treatments, and let ${Y}^{(h)}\
(\in\mathbb{R})$ be the potential outcome variable when the $h$-th treatment
is received. In addition, let $T^{(h)}$ be the assignment variable that is $1$
when $Y^{(h)}$ is observed and $0$ when it is not ($h\in\\{1,\ldots,H\\},\
\sum_{h=1}^{H}T^{(h)}=1$). Moreover, let $\bm{X}\ (\in\mathbb{R}^{p})$ be the
confounding variable vector, i.e., $p$ kinds of confounding variables are
observed; we consider the model,

$\displaystyle{Y}\equiv\sum_{h=1}^{H}T^{(h)}{Y}^{(h)}=\sum_{h=1}^{H}T^{(h)}\\{{\mu}^{(h)}(\bm{X})+{f}(\bm{X})+{\epsilon}^{(h)}\\}.$
(1)

Here, $f:\mathbb{R}^{p}\to\mathbb{R}$ is a nonlinear function and
${\epsilon}^{(h)}$ is the error variable that follows a normal distribution
${\rm N}({0},{\sigma^{2}})$ independently of $\bm{X}$. Note that ${Y}$ is the
observed outcome variable. Letting $(c^{(1)},\ldots,c^{(H)})^{\prime}\
(\in\mathbb{R}^{H})$ be a contrast satisfying $\sum_{h=1}^{H}c^{(h)}=0$, we
suppose that the object of interest is the conditional average treatment
effect $\sum_{h=1}^{H}c^{(h)}{\mu}^{(h)}(\bm{x})$ for $\bm{X}=\bm{x}$. If
$H=2$, and $h=1$ and $h=2$ represent the control and treatment groups,
respectively, and if $\mu^{(1)}(\bm{x})$ and $\mu^{(2)}(\bm{x})$ do not depend
on $\bm{x}$, then this model is the simplest Rubin’s causal inference model,
and the causal effect $\mu^{(2)}-\mu^{(1)}$ defined by
$(c^{(1)},c^{(2)})^{\prime}=(-1,1)^{\prime}$ becomes the simple one estimated
in Rosenbaum and Rubin (1983) or Rubin (1985).

Suppose there are $n$ samples that follow this model and the variables in the
$i$-th sample are identified by the subscript i. Letting
$\tilde{\bm{Y}}=({Y}_{1},\ldots,{Y}_{n})^{\prime}$, $\tilde{\bm{T}}^{(h)}={\rm
diag}(T_{1}^{(h)},\ldots,T_{n}^{(h)})$,
$\tilde{\bm{Y}}^{(h)}=({Y}_{1}^{(h)},\ldots,{Y}_{n}^{(h)})^{\prime}$,
$\tilde{\bm{X}}=(\bm{X}_{1},\ldots,\bm{X}_{n})^{\prime}$,
$\tilde{\bm{\mu}}^{(h)}(\tilde{\bm{X}})=({\mu}^{(h)}(\bm{X}_{1}),\ldots,{\mu}^{(h)}(\bm{X}_{n}))^{\prime}$,
$\tilde{\bm{f}}(\tilde{\bm{X}})=({f}(\bm{X}_{1}),\allowbreak\ldots,{f}(\bm{X}_{n}))^{\prime}$
and
$\tilde{\bm{\epsilon}}^{(h)}=({\epsilon}_{1}^{(h)},\ldots,{\epsilon}_{n}^{(h)})^{\prime}$,
we express the model as

$\displaystyle\tilde{\bm{Y}}=\sum_{h=1}^{H}\tilde{\bm{T}}^{(h)}\tilde{\bm{Y}}^{(h)}=\sum_{h=1}^{H}\tilde{\bm{T}}^{(h)}\\{\tilde{\bm{\mu}}^{(h)}(\tilde{\bm{X}})+\tilde{\bm{f}}(\tilde{\bm{X}})+\tilde{\bm{\epsilon}}^{(h)}\\}.$
(2)

We will assume the usual conditions for this setting. First, we assume the
weakly ignorable treatment assignment condition (Imbens 2000),

$\displaystyle{Y}_{i}^{(h)}\mathop{\perp\\!\\!\\!\\!\perp}T_{i}^{(h)}\mid\bm{X}_{i}\qquad(i\in\\{1,\ldots,n\\},\
h\in\\{1,\ldots,H\\}).$ (3)

Note that ${Y}_{i}^{(h)}$ in this condition can be changed to
${\epsilon}_{i}^{(h)}$. Second, we assume the positivity condition $0<{\rm
P}(T_{i}^{(h)}=1\mid\bm{X}_{i})<1$ ($i\in\\{1,\ldots,n\\}$). Third, we assume
that the samples are independent, that is,

$\displaystyle(T_{i}^{(h)},\bm{X}_{i},{\epsilon}_{i}^{(h)})\mathop{\perp\\!\\!\\!\\!\perp}(T_{l}^{(h)},\bm{X}_{l},{\epsilon}_{l}^{(h)})\qquad(i\neq
l,\ i,l\in\\{1,\ldots,n\\},\ h\in\\{1,\ldots,H\\}).$

As a natural consequence, this means that
${Y}_{i}\mathop{\perp\\!\\!\\!\\!\perp}{Y}_{l}$ ($i\neq l,\
i,l\in\\{1,\ldots,n\\}$).

In this paper, we will suppose that the dimension $p$ of the confounding
variable vector is high and try to estimate the conditional average treatment
effect with higher accuracy by simultaneously conducting variable selection
and estimation. Specifically, we will use a Lasso-type method (Tibshirani
1996) in Section 3. Note that a model selection method for causal inference
was recently proposed by Baba et al. (2017), but it deals with the selection
of marginal structures and is not used here.

In (2), the $H-1$ potential outcome variables ${Y}_{i}^{(h)}$ with
$T_{i}^{(h)}=0$ are considered to be missing. Since ${\rm
E}({Y}_{i}^{(h)})\neq{\rm E}({Y}_{i}^{(h)}\mid T_{i}^{(h)}=1)$ in general,
simply estimating the causal effect $\sum_{h=1}^{H}c^{(h)}{\mu}^{(h)}(\bm{x})$
by using the least-squares method based on minimization of
$\|\sum_{h=1}^{H}c^{(h)}\tilde{\bm{T}}^{(h)}\tilde{\bm{Y}}-\sum_{h=1}^{H}c^{(h)}\tilde{\bm{\mu}}^{(h)}(\tilde{\bm{X}})\|_{2}^{2}$
would result in a large bias. Here, $\|\cdot\|_{2}$ is the $\ell_{2}$ norm. As
described in Section 1, this paper deals with inverse-probability-weighted
estimation using the propensity score $e^{(h)}(\bm{X}_{i})\equiv{\rm
P}(T_{i}^{(h)}=1\mid\bm{X}_{i})$ to avoid the bias.

In this estimation method, the missing values are pseudo-recovered by
multiplying the observed values by the inverse of the propensity score as
weights; then, the usual estimation is conducted. Specifically, letting
$\tilde{\bm{W}}^{(h)}(\tilde{\bm{T}}^{(h)},\tilde{\bm{X}})\equiv{\rm
diag}\\{T^{(h)}_{1}/e^{(h)}(\bm{X}_{1}),\ldots,T^{(h)}_{n}/e^{(h)}(\bm{X}_{n})\\}$
and $\tilde{\bm{T}}=(\tilde{\bm{T}}^{(1)},\ldots,\tilde{\bm{T}}^{(H)})$ and
defining
$\tilde{\bm{W}}(\tilde{\bm{T}},\tilde{\bm{X}})\equiv\sum_{h=1}^{H}c^{(h)}\tilde{\bm{W}}^{(h)}(\tilde{\bm{T}}^{(h)},\tilde{\bm{X}})$
as the weight matrix, the inverse-probability-weighted estimator can be found
by minimizing the squared loss function,

$\displaystyle\bigg{\|}\tilde{\bm{W}}(\tilde{\bm{T}},\tilde{\bm{X}})\tilde{\bm{Y}}-\sum_{h=1}^{H}c^{(h)}\tilde{\bm{\mu}}^{(h)}(\tilde{\bm{X}})\bigg{\|}_{2}^{2}.$
(4)

When the propensity scores are unknown, we assume some parametric function for
$e^{(h)}(\bm{X}_{i})$, maximize the likelihood function
$\prod_{i=1}^{n}\prod_{h=1}^{H}e^{(h)}(\bm{X}_{i})^{T_{i}^{(h)}}$ to obtain
the estimator $\hat{e}^{(h)}(\bm{X}_{i})$, and use it instead. The inverse-
probability-weighted estimator is consistent under the weakly ignorable
treatment assignment condition in (3).

###  2.2 Selective inference

In Section 2.1, we suggested the use of variable selection. After the variable
selection, in order to measure the extent to which the selected variables have
an impact on the causal effects, tests are performed or confidence intervals
are constructed. However, the problem is that the $p$-values used in this
process are no longer reliable, because the selected variables are likely to
be significant, or in other words, the model using the selected variables
likely overfits the data.

To explain the inference after variable selection in detail, we will omit the
superscript (h) as $H=1$ in (2); i.e., $\tilde{\bm{T}}$ is an $n\times n$
identity matrix $\bm{I}_{n}$. Let $\tilde{\bm{X}}$ be a non-random matrix
$\tilde{\bm{x}}$, and let $\tilde{\bm{f}}(\tilde{\bm{x}})$ be an
$n$-dimensional zero vector $\bm{0}_{n}$. Moreover, we will consider a linear
function of $\tilde{\bm{x}}$ as a model for
$\tilde{\bm{\mu}}(\tilde{\bm{x}})$. Then, for the subset
$\bm{M}\subset\\{1,\ldots,p\\}$, we define an estimand as

$\displaystyle\bm{\beta}^{\ddagger\bm{M}}\equiv\mathop{\rm
argmin}_{\bm{b}^{\ddagger\bm{M}}}{\rm
E}(\|\tilde{\bm{Y}}-\tilde{\bm{x}}_{\bm{M}}\bm{b}^{\ddagger\bm{M}}\|_{2}^{2})=(\tilde{\bm{x}}_{\bm{M}}^{\prime}\tilde{\bm{x}}_{\bm{M}})^{-1}\tilde{\bm{x}}_{\bm{M}}^{\prime}\tilde{\bm{\mu}}(\tilde{\bm{x}}),$

which minimizes the expected squared error for the linear sum of the
confounding variables belonging to $\bm{M}$, where
$\tilde{\bm{x}}_{\bm{M}}=(\tilde{x}_{ij})_{i\in\\{1,\ldots,n\\},j\in\bm{M}}$.
If the models $\bm{M}$ and $\bm{M}^{*}$ are different, then
$\beta_{j}^{\ddagger\bm{M}}$ and $\beta_{j}^{\ddagger\bm{M}^{*}}$ are
generally different, i.e., the target of inference is different for each model
selected. Consequently, we can say that inference after selection has some
ambiguity.

In terms of confidence interval construction, the simultaneous inference
mentioned in Section 1 is a method of creating an interval in which all the
regression coefficients regardless of which model is selected are included
with probability $1-\alpha$ or greater. When $\hat{\bm{M}}$ is the selected
model and the regression coefficients in that model are denoted as
$\beta_{j}^{\ddagger\hat{\bm{M}}}\ (j\in\hat{\bm{M}})$, this method finds
$C_{j}^{\ddagger\hat{\bm{M}}}\ (j\in\hat{\bm{M}})$ such that

$\displaystyle{\rm P}(\forall j\in\hat{\bm{M}},\
\beta_{j}^{\ddagger\hat{\bm{M}}}\in C_{j}^{\ddagger\hat{\bm{M}}})\geq
1-\alpha.$

Since the interval is guaranteed for all $j$, a problem arises in that the
interval becomes considerably wide when the size $|\hat{\bm{M}}|$ is large.

On the other hand, selective inference, in terms of confidence interval
construction, is a method of creating an interval in which the regression
coefficients in the selected model are included with probability $1-\alpha$ or
greater under the condition that the model was selected. The method is to find
$C_{j}^{\ddagger\bm{M}}\ (j\in\bm{M})$ such that

$\displaystyle{\rm P}(\beta_{j}^{\ddagger\bm{M}}\in
C_{j}^{\ddagger\bm{M}}\mid\hat{\bm{M}}=\bm{M})\geq 1-\alpha.$

Accordingly, the false coverage rate is controlled to be

$\displaystyle{\rm
E}\bigg{(}\frac{|\\{j\in\hat{\bm{M}}:\beta_{j}^{\ddagger\hat{\bm{M}}}\notin
C_{j}^{\ddagger\hat{\bm{M}}}\\}|}{|\hat{\bm{M}}|};\
|\hat{\bm{M}}|>0\bigg{)}\leq\alpha.$ (5)

In the cornerstone papers, beautiful inferences were presented for marginal
screening and Lasso model selection in the setting of normal distributions.
Then, their validity was theoretically demonstrated in more generalized
settings. For example, Fithian et al. (2014) showed the optimality of the
method in terms of uniformly most powerful unbiasedness in the setting of
exponential families of distributions, while Tian and Taylor (2017), Tian and
Taylor (2018), and Tibshirani et al. (2018) showed that the distributions of
the pivot statistics for getting $p$-values or confidence intervals converge
uniformly to the actual distribution when a normal distribution is used, even
when the noise has a non-normal distribution. Furthermore, Liu et al. (2018)
proposed to increase the power of the original method by reducing the over-
conditioning. In this paper, we develop selective inference for causal
inference on the basis of Lee et al. (2016).

We will explain the method of Lee et al. (2016). Let us denote the usual Lasso
estimator by
$\hat{\bm{\beta}}^{\ddagger}=(\hat{\beta}_{1}^{\ddagger},\ldots,\hat{\beta}_{p}^{\ddagger})^{\prime}$,
the collection of non-zero estimators by
$\hat{\bm{\beta}}^{\ddagger\hat{\bm{M}}}$, and its sign by
$\hat{\bm{s}}^{\ddagger\hat{\bm{M}}}={\rm
sign}(\hat{\bm{\beta}}^{\ddagger\hat{\bm{M}}})$. In this case, from the
Karuch-Kuhn-Tucker conditions, for any $\bm{s}\in\\{-1,1\\}^{|\bm{M}|}$, there
exists an $n\times n$ matrix $\bm{A}(\bm{M},\bm{s})$ and an $n$-dimensional
vector $\bm{b}(\bm{M},\bm{s})$ such that

$\displaystyle\\{\hat{\bm{M}}^{\ddagger}=\bm{M},\
\hat{\bm{s}}^{\ddagger\hat{\bm{M}}}=\bm{s}\\}=\\{\bm{A}(\bm{M},\bm{s})\tilde{\bm{Y}}\leq\bm{b}(\bm{M},\bm{s})\\}.$

Now, using an appropriate unit vector $\bm{e}_{j}\
(\in\mathbb{R}^{|\bm{M}|})$, we define
$\tilde{\bm{\eta}}_{j}^{\ddagger}\equiv\tilde{\bm{x}}_{\bm{M}}(\tilde{\bm{x}}_{\bm{M}}^{\prime}\tilde{\bm{x}}_{\bm{M}})^{-1}\bm{e}_{j}$.
Since the target parameter can be written as
$\beta_{j}^{\ddagger\bm{M}}=\tilde{\bm{\eta}}_{j}^{\ddagger}{}^{\prime}\tilde{\bm{\mu}}(\tilde{\bm{x}})$,
we use the statistic
$\tilde{\bm{\eta}}_{j}^{\ddagger}{}^{\prime}\tilde{\bm{Y}}$ to create its
confidence interval. Then, by using the strong properties of a normal
distribution, we obtain

$\displaystyle
F_{\beta_{j}^{\ddagger\bm{M}},\sigma^{2}\tilde{\bm{\eta}}_{j}^{\ddagger}{}^{\prime}\tilde{\bm{\eta}}_{j}^{\ddagger}}^{[{\mathcal{V}}_{\bm{s},j}^{-}(\tilde{\bm{Z}}),{\mathcal{V}}_{\bm{s},j}^{+}(\tilde{\bm{Z}})]}(\tilde{\bm{\eta}}_{j}^{\ddagger}{}^{\prime}\tilde{\bm{Y}})\mid\\{\bm{A}(\bm{M},\bm{s})\tilde{\bm{Y}}\leq\bm{b}(\bm{M},\bm{s})\\}\sim{\rm
Unif}(0,1).$

Here, $F_{\mu,\sigma^{2}}^{[a,b]}(\cdot)$ denotes the cumulative distribution
function of ${\rm TN}(\mu,\sigma^{2},a,b)$, which is ${\rm N}(\mu,\sigma^{2})$
truncated into the interval $[a,b]$,
$\tilde{\bm{Z}}=\\{\bm{I}_{n}-\tilde{\bm{\eta}}_{j}^{\ddagger}(\tilde{\bm{\eta}}_{j}^{\ddagger}{}^{\prime}\tilde{\bm{\eta}}_{j}^{\ddagger})^{-1}\tilde{\bm{\eta}}_{j}^{\ddagger}{}^{\prime}\\}\tilde{\bm{Y}}$,

${\mathcal{V}}_{\bm{s},j}^{-}(\tilde{\bm{Z}})=\max_{k:(\bm{A}(\bm{M},\bm{s})\tilde{\bm{\eta}}_{j}^{\ddagger}(\tilde{\bm{\eta}}_{j}^{\ddagger}{}^{\prime}\tilde{\bm{\eta}}_{j}^{\ddagger})^{-1})_{k}<0}\\{b(\bm{M},\bm{s})_{k}-(\bm{A}(\bm{M},\bm{s})\tilde{\bm{Z}})_{k}\\}/(\bm{A}(\bm{M},\bm{s})\tilde{\bm{\eta}}_{j}^{\ddagger}(\tilde{\bm{\eta}}_{j}^{\ddagger}{}^{\prime}\tilde{\bm{\eta}}_{j}^{\ddagger})^{-1})_{k},$
${\mathcal{V}}_{\bm{s},j}^{+}(\tilde{\bm{Z}})=\min_{k:(\bm{A}(\bm{M},\bm{s})\tilde{\bm{\eta}}_{j}^{\ddagger}(\tilde{\bm{\eta}}_{j}^{\ddagger}{}^{\prime}\tilde{\bm{\eta}}_{j}^{\ddagger})^{-1})_{k}>0}\\{b(\bm{M},\bm{s})_{k}-(\bm{A}\allowbreak(\bm{M},\bm{s})\tilde{\bm{Z}})_{k}\\}/(\bm{A}(\bm{M},\bm{s})\tilde{\bm{\eta}}_{j}^{\ddagger}(\tilde{\bm{\eta}}_{j}^{\ddagger}{}^{\prime}\tilde{\bm{\eta}}_{j}^{\ddagger})^{-1})_{k},$

and ${\rm Unif}(0,1)$ denotes a continuous uniform distribution in $[0,1]$. We
have obtained a pivot statistic that is unique to a normal distribution, and
from it we can construct a statistical test or a confidence interval for
$\beta_{j}^{\ddagger\bm{M}}$, with the significance level or coverage
probability controlled in a conditional sense.

##  3 Main result

Consider selective inference for the causal inference model in (2). This
problem has been addressed by Zhao et al. (2017), who do not use propensity
scores, but perform nonparametric estimation of
$\tilde{\bm{f}}(\tilde{\bm{X}})$. One of the reasons why causal inference has
become a hot topic is the rapid development of propensity score analysis,
which does not necessarily require this nonparametric estimation, and so, in
this section, we will avoid the nonparametric estimation as well. This
avoidance would require an essential generalization of selective inference,
unlike Zhao et al. (2017).

Now, for $\tilde{\bm{\mu}}^{(h)}(\tilde{\bm{X}})$ in (2), let us choose a
linear function of $\tilde{\bm{X}}$ as a model. For the subset
$\bm{M}\subset\\{1,\ldots,p\\}$, because the causal effect is
$\sum_{h=1}^{H}c^{(h)}\mu^{(h)}(\bm{X})$, we define an estimand as

$\displaystyle\bm{\beta}^{\bm{M}}\equiv\mathop{\rm
argmin}_{\bm{b}^{\bm{M}}}{\rm
E}\bigg{(}\bigg{\|}\sum_{h=1}^{H}c^{(h)}\tilde{\bm{Y}}^{(h)}-\tilde{\bm{X}}_{\bm{M}}\bm{b}^{\bm{M}}\bigg{\|}_{2}^{2}\
\bigg{|}\
\tilde{\bm{X}}\bigg{)}=(\tilde{\bm{X}}_{\bm{M}}^{\prime}\tilde{\bm{X}}_{\bm{M}})^{-1}\tilde{\bm{X}}_{\bm{M}}^{\prime}\sum_{h=1}^{H}c^{(h)}\tilde{\bm{\mu}}^{(h)}(\tilde{\bm{X}}),$

which minimizes the expected squared error for the linear sum of the
confounding variables belonging to $\bm{M}$, where
$\tilde{\bm{X}}_{\bm{M}}=(\tilde{X}_{ij})_{i\in\\{1,\ldots,n\\},j\in\bm{M}}$.
Note that, in this model, the causal effect is captured by
$\bm{x}_{\bm{M}}^{\prime}\bm{\beta}^{\bm{M}}$. Then, denoting the selected
model as $\hat{\bm{M}}$, we try to find $C_{j}^{\bm{M}}\ (j\in\bm{M})$, which
gives the conditional coverage,

$\displaystyle{\rm P}(\beta_{j}^{\bm{M}}\in
C_{j}^{\bm{M}}\mid\hat{\bm{M}}=\bm{M})\geq 1-\alpha$ (6)

for the regression coefficient $\beta_{j}^{\bm{M}}\ (j\in\bm{M})$.

This allows us to construct a statistical test or a confidence interval for
$\beta_{j}^{\bm{M}}$. Using an appropriate unit vector $\bm{e}_{j}\
(\in\mathbb{R}^{|\bm{M}|})$, we have
$\beta_{j}^{\bm{M}}=\bm{e}_{j}^{\prime}(\tilde{\bm{X}}_{\bm{M}}^{\prime}\tilde{\bm{X}}_{\bm{M}})^{-1}\tilde{\bm{X}}_{\bm{M}}^{\prime}\sum_{h=1}^{H}c^{(h)}\tilde{\bm{\mu}}^{(h)}(\tilde{\bm{X}})$.
Then, defining
$\tilde{\bm{\eta}}_{j}\equiv\tilde{\bm{X}}_{\bm{M}}(\tilde{\bm{X}}_{\bm{M}}^{\prime}\tilde{\bm{X}}_{\bm{M}})^{-1}\bm{e}_{j}$,
later we consider the distribution of
$\tilde{\bm{\eta}}_{j}^{\prime}\tilde{\bm{W}}(\tilde{\bm{T}},\tilde{\bm{X}})\tilde{\bm{Y}}$
conditional on $\hat{\bm{M}}=\bm{M}$. Since it is not possible to perfectly
capture the true model by $\hat{\bm{M}}$ in finite data, and the true
structure may not be included in the model,
$\tilde{\bm{X}}_{\bm{M}}\bm{\beta}^{\bm{M}}$ is not itself the conditional
average treatment effect
$\sum_{h=1}^{H}c^{(h)}\tilde{\bm{\mu}}^{(h)}(\tilde{\bm{X}})$. However, it is
a problem with selective inference itself, and it is not addressed in this
paper.

As a model selection method, by letting
$\bm{\beta}=\sum_{h=1}^{H}c^{(h)}\bm{\beta}^{(h)}$ after using
$\tilde{\bm{\mu}}^{(h)}(\tilde{\bm{X}})=\tilde{\bm{X}}\bm{\beta}^{(h)}$ in
(4), we propose a Lasso-type one that estimates $\bm{\beta}$ by

$\displaystyle\hat{\bm{\beta}}=\mathop{\rm
argmin}_{\bm{\beta}}\bigg{\\{}\frac{1}{2}\|\tilde{\bm{W}}(\tilde{\bm{T}},\tilde{\bm{X}})\tilde{\bm{Y}}-\tilde{\bm{X}}\bm{\beta}\|_{2}^{2}+\lambda\|\bm{\beta}\|_{1}\bigg{\\}}.$
(7)

Note that $\tilde{\bm{Y}}$ has a nonlinear component
$\tilde{\bm{f}}(\tilde{\bm{X}})$, however, this nonlinear component is
canceled out in $\tilde{\bm{W}}(\tilde{\bm{T}},\tilde{\bm{X}})\tilde{\bm{Y}}$
and there is no problem in using Lasso for linear models. Since Lasso gives a
sparse solution, the selected model can be expressed as

$\displaystyle\hat{\bm{M}}=\\{j:\ \hat{\beta}_{j}\neq 0\\}.$

Here, we write
$\hat{\bm{\beta}}^{\hat{\bm{M}}}=(\hat{\beta}_{j})_{j\in\hat{\bm{M}}}$ for the
collection of non-zero estimators and $\hat{\bm{s}}^{\hat{\bm{M}}}={\rm
sign}(\hat{\bm{\beta}}^{\hat{\bm{M}}})$ for their sign.

The fact that $\hat{\bm{\beta}}$ is the solution of (7) is equivalent to the
existence of $\hat{\bm{s}}=(\hat{s}_{1},\ldots,\hat{s}_{p})^{\prime}$
satisfying the Karush-Kuhn-Tucker conditions,

$\displaystyle\hat{\beta}_{j}\neq 0\quad\Rightarrow\quad\hat{s}_{j}={\rm
sign}(\hat{\beta}_{j}),$
$\displaystyle\hat{\beta}_{j}=0\quad\Rightarrow\quad\hat{s}_{j}\in[-1,1]$

and

$\displaystyle\tilde{\bm{X}}^{\prime}\\{\tilde{\bm{X}}\hat{\bm{\beta}}-\tilde{\bm{W}}(\tilde{\bm{T}},\tilde{\bm{X}})\tilde{\bm{Y}}\\}+\lambda\hat{\bm{s}}=\bm{0}.$

For the model $\bm{M}$, if the candidate for the sign of its regression
coefficients is $\bm{s}\ (\in\\{-1,1\\}^{|\bm{M}|})$, then we can use this
equivalence conditions to show

$\displaystyle\\{\hat{\bm{M}}=\bm{M},\
\hat{\bm{s}}^{\hat{\bm{M}}}=\bm{s}\\}=\\{\bm{A}(\bm{M},\bm{s})\tilde{\bm{W}}(\tilde{\bm{T}},\tilde{\bm{X}})\tilde{\bm{Y}}\leq\bm{b}(\bm{M},\bm{s})\\}.$
(8)

Here, we use the definitions,

$\displaystyle\bm{A}(\bm{M},\bm{s})=\frac{1}{\lambda}\left(\begin{array}[]{c}\tilde{\bm{X}}_{\bm{M}^{\rm
c}}^{\prime}\\{\bm{I}_{n}-\tilde{\bm{X}}_{\bm{M}}(\tilde{\bm{X}}_{\bm{M}}^{\prime}\tilde{\bm{X}}_{\bm{M}})^{-1}\tilde{\bm{X}}_{\bm{M}}^{\prime}\\}\\\
-\tilde{\bm{X}}_{\bm{M}^{\rm
c}}^{\prime}\\{\bm{I}_{n}-\tilde{\bm{X}}_{\bm{M}}(\tilde{\bm{X}}_{\bm{M}}^{\prime}\tilde{\bm{X}}_{\bm{M}})^{-1}\tilde{\bm{X}}_{\bm{M}}^{\prime}\\}\\\
-{\rm
diag}(\bm{s})(\tilde{\bm{X}}_{\bm{M}}^{\prime}\tilde{\bm{X}}_{\bm{M}})^{-1}\tilde{\bm{X}}_{\bm{M}}^{\prime}\end{array}\right)$

and

$\displaystyle\bm{b}(\bm{M},\bm{s})=\left(\begin{array}[]{c}\bm{1}_{|\bm{M}^{\rm
c}|}-\tilde{\bm{X}}_{\bm{M}^{\rm
c}}^{\prime}\tilde{\bm{X}}_{\bm{M}}(\tilde{\bm{X}}_{\bm{M}}^{\prime}\tilde{\bm{X}}_{\bm{M}})^{-1}\bm{s}\\\
\bm{1}_{|\bm{M}^{\rm c}|}+\tilde{\bm{X}}_{\bm{M}^{\rm
c}}^{\prime}\tilde{\bm{X}}_{\bm{M}}(\tilde{\bm{X}}_{\bm{M}}^{\prime}\tilde{\bm{X}}_{\bm{M}})^{-1}\bm{s}\\\
-{\rm
diag}(\bm{s})(\tilde{\bm{X}}_{\bm{M}}^{\prime}\tilde{\bm{X}}_{\bm{M}})^{-1}\bm{s}\end{array}\right),$

where $\bm{M}^{\rm c}=\\{1,\ldots,p\\}\setminus\bm{M}$,
$\tilde{\bm{\bm{X}}}_{\bm{M}^{\rm
c}}=(\tilde{X}_{ij})_{i\in\\{1,\ldots,n\\},j\in\bm{M}^{\rm c}}$ and
$\bm{1}_{|\bm{M}^{\rm c}|}$ is a $|\bm{M}^{\rm c}|$-dimensional one-vector.
Its derivation method is the same as the one used by Lee et al. (2016). Hence,
we also have

$\displaystyle\\{\hat{\bm{M}}=\bm{M}\\}=\bigcup_{\bm{s}\in\\{-1,1\\}^{\bm{M}}}\\{\bm{A}(\bm{M},\bm{s})\tilde{\bm{W}}(\tilde{\bm{T}},\tilde{\bm{X}})\tilde{\bm{Y}}\leq\bm{b}(\bm{M},\bm{s})\\}.$
(9)

Although we mentioned that we would consider the conditional coverage in (6),
conditioning also on the sign $\hat{\bm{s}}^{\hat{\bm{M}}}$ leads to an easier
problem; hence, we will consider

$\displaystyle{\rm P}(\beta_{j}^{\bm{M}}\in
C_{j}^{\bm{M}}\mid\hat{\bm{M}}=\bm{M},\ \hat{\bm{s}}^{\hat{\bm{M}}}=\bm{s})$
$\displaystyle={\rm P}\\{\beta_{j}^{\bm{M}}\in
C_{j}^{\bm{M}}\mid\bm{A}(\bm{M},\bm{s})\tilde{\bm{W}}(\tilde{\bm{T}},\tilde{\bm{X}})\tilde{\bm{Y}}\leq\bm{b}(\bm{M},\bm{s})\\}\geq
1-\alpha$

as in Lee et al. (2016). This interval $C_{j}^{\bm{M}}$ also satisfies (6),
but if one wants tighter coverage, one can condition on (9). To find
$C_{j}^{\bm{M}}$, we consider the distribution of
$\tilde{\bm{\eta}}_{j}^{\prime}\tilde{\bm{W}}(\tilde{\bm{T}},\tilde{\bm{X}})\tilde{\bm{Y}}$
conditional on (8). Let
$\tilde{\bm{D}}(\tilde{\bm{T}},\tilde{\bm{X}})\equiv\tilde{\bm{W}}(\tilde{\bm{T}},\tilde{\bm{X}})^{2}\tilde{\bm{\eta}}_{j}\\{\tilde{\bm{\eta}}_{j}^{\prime}\tilde{\bm{W}}(\tilde{\bm{T}},\tilde{\bm{X}})^{2}\tilde{\bm{\eta}}_{j}\\}^{-1}$.
Using the same derivation method as in Lee et al. (2016), we can rewrite the
conditioned polyhedral region by using
$\tilde{\bm{\eta}}_{j}^{\prime}\tilde{\bm{W}}(\tilde{\bm{T}},\tilde{\bm{X}})\tilde{\bm{Y}}$
and

$\displaystyle\tilde{\bm{Z}}\equiv\\{\bm{I}_{n}-\tilde{\bm{D}}(\tilde{\bm{T}},\tilde{\bm{X}})\tilde{\bm{\eta}}_{j}^{\prime}\\}\tilde{\bm{W}}(\tilde{\bm{T}},\tilde{\bm{X}})\tilde{\bm{Y}},$

from which we get

$\displaystyle\\{\bm{A}(\bm{M},\bm{s})\tilde{\bm{W}}(\tilde{\bm{T}},\tilde{\bm{X}})\tilde{\bm{Y}}\leq\bm{b}(\bm{M},\bm{s})\\}$
$\displaystyle=\\{{\mathcal{V}}_{\bm{s},j}^{-}(\tilde{\bm{Z}},\tilde{\bm{T}},\tilde{\bm{X}})\leq\tilde{\bm{\eta}}_{j}^{\prime}\tilde{\bm{W}}(\tilde{\bm{T}},\tilde{\bm{X}})\tilde{\bm{Y}}\leq{\mathcal{V}}_{\bm{s},j}^{+}(\tilde{\bm{Z}},\tilde{\bm{T}},\tilde{\bm{X}}),\
{\mathcal{V}}_{\bm{s},j}^{0}(\tilde{\bm{Z}},\tilde{\bm{T}},\tilde{\bm{X}})\geq
0\\},$

where

$\displaystyle{\mathcal{V}}_{\bm{s},j}^{-}(\tilde{\bm{Z}},\tilde{\bm{T}},\tilde{\bm{X}})\equiv\max_{k:(\bm{A}(\bm{M},\bm{s})\tilde{\bm{D}}(\tilde{\bm{T}},\tilde{\bm{X}}))_{k}<0}\frac{\bm{b}(\bm{M},\bm{s})_{k}-(\bm{A}(\bm{M},\bm{s})\tilde{\bm{Z}})_{k}}{(\bm{A}(\bm{M},\bm{s})\tilde{\bm{D}}(\tilde{\bm{T}},\tilde{\bm{X}}))_{k}},$
(10)
$\displaystyle{\mathcal{V}}_{\bm{s},j}^{+}(\tilde{\bm{Z}},\tilde{\bm{T}},\tilde{\bm{X}})\equiv\min_{k:(\bm{A}(\bm{M},\bm{s})\tilde{\bm{D}}(\tilde{\bm{T}},\tilde{\bm{X}}))_{k}>0}\frac{\bm{b}(\bm{M},\bm{s})_{k}-(\bm{A}(\bm{M},\bm{s})\tilde{\bm{Z}})_{k}}{(\bm{A}(\bm{M},\bm{s})\tilde{\bm{D}}(\tilde{\bm{T}},\tilde{\bm{X}}))_{k}}$
(11)

and

$\displaystyle{\mathcal{V}}_{\bm{s},j}^{0}(\tilde{\bm{Z}},\tilde{\bm{T}},\tilde{\bm{X}})\equiv\min_{k:(\bm{A}(\bm{M},\bm{s})\tilde{\bm{D}}(\tilde{\bm{T}},\tilde{\bm{X}}))_{k}=0}\\{\bm{b}(\bm{M},\bm{s})_{k}-(\bm{A}(\bm{M},\bm{s})\tilde{\bm{Z}})_{k}\\}.$

Given this setup, we want to construct a statistical test or a confidence
interval for the target parameter
$\beta_{j}^{\bm{M}}=\bm{e}_{j}^{\prime}(\tilde{\bm{X}}_{\bm{M}}^{\prime}\tilde{\bm{X}}_{\bm{M}})^{-1}\tilde{\bm{X}}_{\bm{M}}^{\prime}\sum_{h=1}^{H}c^{(h)}\allowbreak\tilde{\bm{\mu}}^{(h)}(\tilde{\bm{X}})$.
However, unlike in Section 2.2, $\tilde{\bm{Y}}$ is not normally distributed
and $\tilde{\bm{f}}(\tilde{\bm{X}})$ exists. Also, $\tilde{\bm{X}}$ is random,
although it is a trivial difference. Therefore, we consider using the
realizations $\tilde{\bm{t}}$ and $\tilde{\bm{x}}$ to further condition on
that $\tilde{\bm{T}}=\tilde{\bm{t}}$ and $\tilde{\bm{X}}=\tilde{\bm{x}}$. This
is a difficult operation to understand because, $\tilde{\bm{T}}$ is not
conditioned first in the usual propensity score analysis, but it is
essentially necessary to extract the properties of a normal distribution. It
is also difficult to intuitively understand how such conditioning guarantees
valid inferences, but as we will see later, it is easy to show that it does.
In general, under arbitrary conditioning, inference about the target parameter
becomes impossible, but in this setting, under the weakly ignorable treatment
assignment condition, inference based on asymptotics is possible. Also, as in
the usual propensity score analysis, it is independent of
$\tilde{\bm{f}}(\tilde{\bm{X}})$.

Under this new conditioning,
$\tilde{\bm{\eta}}_{j}^{\prime}\tilde{\bm{W}}(\tilde{\bm{T}},\tilde{\bm{X}})\tilde{\bm{Y}}$
and $\tilde{\bm{Z}}$ become independent; thus, letting $\tilde{\bm{z}}$ be the
realization of $\tilde{\bm{Z}}$, we have

$\displaystyle[\tilde{\bm{\eta}}_{j}^{\prime}\tilde{\bm{W}}(\tilde{\bm{T}},\tilde{\bm{X}})\tilde{\bm{Y}}\mid\bm{A}(\bm{M},\bm{s})\tilde{\bm{W}}(\tilde{\bm{t}},\tilde{\bm{x}})\tilde{\bm{Y}}\leq\bm{b}(\bm{M},\bm{s}),\
\tilde{\bm{Z}}=\tilde{\bm{z}},\ \tilde{\bm{T}}=\tilde{\bm{t}},\
\tilde{\bm{X}}=\tilde{\bm{x}}]$ $\displaystyle\stackrel{{\scriptstyle\rm
d}}{{=}}[\tilde{\bm{\eta}}_{j}^{\prime}\tilde{\bm{W}}(\tilde{\bm{t}},\tilde{\bm{x}})\tilde{\bm{Y}}\mid{\mathcal{V}}_{\bm{s},j}^{-}(\tilde{\bm{z}},\tilde{\bm{t}},\tilde{\bm{x}})\leq\tilde{\bm{\eta}}_{j}^{\prime}\tilde{\bm{W}}(\tilde{\bm{t}},\tilde{\bm{x}})\tilde{\bm{Y}}\leq{\mathcal{V}}_{\bm{s},j}^{+}(\tilde{\bm{z}},\tilde{\bm{t}},\tilde{\bm{x}}),\
\tilde{\bm{T}}=\tilde{\bm{t}},\ \tilde{\bm{X}}=\tilde{\bm{x}}].$ (12)

In addition, because the ignorable treatment assignment condition is
satisfied, $\tilde{\bm{\epsilon}}^{(h)}$ follows ${\rm
N}(\bm{0}_{n},\allowbreak\sigma^{2}\bm{I}_{n})$ under this conditioning.
Hence, by letting

$\displaystyle\kappa_{j}^{\bm{M}}(\tilde{\bm{t}},\tilde{\bm{x}})\equiv\bm{e}_{j}^{\prime}(\tilde{\bm{x}}_{\bm{M}}^{\prime}\tilde{\bm{x}}_{\bm{M}})^{-1}\tilde{\bm{x}}_{\bm{M}}^{\prime}\sum_{h=1}^{H}c^{(h)}\tilde{\bm{W}}^{(h)}(\tilde{\bm{t}}^{(h)},\tilde{\bm{x}})\\{\tilde{\bm{\mu}}^{(h)}(\tilde{\bm{x}})+\tilde{\bm{f}}(\tilde{\bm{x}})\\}$
(13)

and

$\displaystyle\zeta_{j}^{\bm{M}}(\tilde{\bm{t}},\tilde{\bm{x}})\equiv\sigma^{2}\bm{e}_{j}^{\prime}(\tilde{\bm{x}}_{\bm{M}}^{\prime}\tilde{\bm{x}}_{\bm{M}})^{-1}\tilde{\bm{x}}_{\bm{M}}^{\prime}\tilde{\bm{W}}(\tilde{\bm{t}},\tilde{\bm{x}})^{2}\tilde{\bm{x}}_{\bm{M}}(\tilde{\bm{x}}_{\bm{M}}^{\prime}\tilde{\bm{x}}_{\bm{M}})^{-1}\bm{e}_{j},$
(14)

it can be seen that the conditional distribution is a truncated normal
distribution,

$\displaystyle{\rm
TN}(\kappa_{j}^{\bm{M}}(\tilde{\bm{t}},\tilde{\bm{x}}),\zeta_{j}^{\bm{M}}(\tilde{\bm{t}},\tilde{\bm{x}}),{\mathcal{V}}_{\bm{s},j}^{-}(\tilde{\bm{z}},\tilde{\bm{t}},\tilde{\bm{x}}),{\mathcal{V}}_{\bm{s},j}^{+}(\tilde{\bm{z}},\tilde{\bm{t}},\tilde{\bm{x}})).$

Thus, letting
$F^{\tilde{\bm{Z}},\tilde{\bm{T}},\tilde{\bm{X}}}(\tilde{\bm{\eta}}_{j}^{\prime}\tilde{\bm{W}}(\tilde{\bm{T}},\tilde{\bm{X}})\tilde{\bm{Y}})\equiv
F_{\kappa_{j}^{\bm{M}}(\tilde{\bm{T}},\tilde{\bm{X}}),\zeta_{j}^{\bm{M}}(\tilde{\bm{T}},\tilde{\bm{X}})}^{[{\mathcal{V}}_{\bm{s},j}^{-}(\tilde{\bm{Z}},\tilde{\bm{T}},\tilde{\bm{X}}),{\mathcal{V}}_{\bm{s},j}^{+}(\tilde{\bm{Z}},\tilde{\bm{T}},\tilde{\bm{X}})]}\allowbreak(\tilde{\bm{\eta}}_{j}^{\prime}\tilde{\bm{W}}(\tilde{\bm{T}},\tilde{\bm{X}})\tilde{\bm{Y}})$,
we obtain

$\displaystyle[F^{\tilde{\bm{Z}},\tilde{\bm{T}},\tilde{\bm{X}}}(\tilde{\bm{\eta}}_{j}^{\prime}\tilde{\bm{W}}(\tilde{\bm{T}},\tilde{\bm{X}})\tilde{\bm{Y}})\mid\bm{A}(\bm{M},\bm{s})\tilde{\bm{W}}(\tilde{\bm{T}},\tilde{\bm{X}})\tilde{\bm{Y}}\leq\bm{b}(\bm{M},\bm{s}),\
\tilde{\bm{Z}}=\tilde{\bm{z}},\ \tilde{\bm{T}}=\tilde{\bm{t}},\
\tilde{\bm{X}}=\tilde{\bm{x}}]$ $\displaystyle\stackrel{{\scriptstyle\rm
d}}{{=}}[F^{\tilde{\bm{z}},\tilde{\bm{t}},\tilde{\bm{x}}}(\tilde{\bm{\eta}}_{j}^{\prime}\tilde{\bm{W}}(\tilde{\bm{t}},\tilde{\bm{x}})\tilde{\bm{Y}})\mid\bm{A}(\bm{M},\bm{s})\tilde{\bm{W}}(\tilde{\bm{t}},\tilde{\bm{x}})\tilde{\bm{Y}}\leq\bm{b}(\bm{M},\bm{s}),\
\tilde{\bm{Z}}=\tilde{\bm{z}},\ \tilde{\bm{T}}=\tilde{\bm{t}},\
\tilde{\bm{X}}=\tilde{\bm{x}}]$ $\displaystyle\sim{\rm Unif}(0,1).$

Denoting the conditional density function of $V$ by $p_{V}$, we have shown
that

$p_{F^{\tilde{\bm{Z}},\tilde{\bm{T}},\tilde{\bm{X}}}(\tilde{\bm{\eta}}_{j}^{\prime}\tilde{\bm{W}}(\tilde{\bm{T}},\tilde{\bm{X}})\tilde{\bm{Y}})\mid\tilde{\bm{Z}},\tilde{\bm{T}},\tilde{\bm{X}}}(v\mid\tilde{\bm{z}},\tilde{\bm{t}},\tilde{\bm{x}})=1_{[0,1]}(v),$

from which it can be seen that

$\displaystyle
p_{F^{\tilde{\bm{Z}},\tilde{\bm{T}},\tilde{\bm{X}}}(\tilde{\bm{\eta}}_{j}^{\prime}\tilde{\bm{W}}(\tilde{\bm{T}},\tilde{\bm{X}})\tilde{\bm{Y}})}(v)$
$\displaystyle=\int\sum_{\tilde{\bm{t}}}\left\\{\int
p_{F^{\tilde{\bm{Z}},\tilde{\bm{T}},\tilde{\bm{X}}}(\tilde{\bm{\eta}}_{j}^{\prime}\tilde{\bm{W}}(\tilde{\bm{T}},\tilde{\bm{X}})\tilde{\bm{Y}})\mid\tilde{\bm{Z}},\tilde{\bm{T}},\tilde{\bm{X}}}\left(v\mid\tilde{\bm{z}},\tilde{\bm{t}},\tilde{\bm{x}}\right)p_{\tilde{\bm{Z}}\mid\tilde{\bm{T}},\tilde{\bm{X}}}\left(\tilde{\bm{z}}\mid\tilde{\bm{t}},\tilde{\bm{x}}\right){\rm
d}\tilde{\bm{z}}\right\\}p_{\tilde{\bm{T}},\tilde{\bm{X}}}\left(\tilde{\bm{t}},\tilde{\bm{x}}\right){\rm
d}\tilde{\bm{x}}$ $\displaystyle=1_{[0,1]}(v).$

This leads us to the following theorem.

######  Theorem 1.

Let $F_{\mu,\sigma^{2}}^{[a,b]}(\cdot)$ be the cumulative distribution
functions of ${\rm TN}(\mu,\sigma^{2},a,b)$, which is ${\rm
N}(\mu,\sigma^{2})$ truncated into the interval $[a,b]$, and let
${\mathcal{V}}_{\bm{s},j}^{-}$, ${\mathcal{V}}_{\bm{s},j}^{+}$,
$\kappa_{j}^{\bm{M}}$ and $\zeta_{j}^{\bm{M}}$ be the functions defined in
(10), (11), (13) and (14), respectively. Then,

$\displaystyle[F_{\kappa_{j}^{\bm{M}}(\tilde{\bm{T}},\tilde{\bm{X}}),\zeta_{j}^{\bm{M}}(\tilde{\bm{T}},\tilde{\bm{X}})}^{[{\mathcal{V}}_{\bm{s},j}^{-}(\tilde{\bm{Z}},\tilde{\bm{T}},\tilde{\bm{X}}),{\mathcal{V}}_{\bm{s},j}^{+}(\tilde{\bm{Z}},\tilde{\bm{T}},\tilde{\bm{X}})]}(\tilde{\bm{\eta}}_{j}^{\prime}\tilde{\bm{W}}(\tilde{\bm{T}},\tilde{\bm{X}})\tilde{\bm{Y}})\mid\bm{A}(\bm{M},\bm{s})\tilde{\bm{W}}(\tilde{\bm{T}},\tilde{\bm{X}})\tilde{\bm{Y}}\leq\bm{b}(\bm{M},\bm{s})]\sim{\rm
Unif}(0,1)$

holds.

The key point here is that, like the conditioning of $\tilde{\bm{Z}}$, the
conditioning of $\tilde{\bm{T}}$ and $\tilde{\bm{X}}$ does not interfere with
the inference.

This theorem cannot be used for the inference about $\beta_{j}^{\bm{M}}$
without further consideration, because unlike the selective inference in Lee
et al. (2016), $\beta_{j}^{\bm{M}}$ does not appear explicitly in this pivot
statistic. Therefore, we consider higher-order asymptotic theory and extract
$\beta_{j}^{\bm{M}}$ from
$\kappa_{j}^{\bm{M}}(\tilde{\bm{T}},\tilde{\bm{X}})$, which appears as the
mean parameter of the truncated normal distribution. Then, we construct an
inference guaranteed asymptotically by evaluating the asymptotic behavior of
$\kappa_{j}^{\bm{M}}(\tilde{\bm{T}},\tilde{\bm{X}})$. Note that this operation
is peculiar to the selective inference for propensity score analysis.

In (13) of $(\tilde{\bm{t}},\tilde{\bm{x}})$ as
$(\tilde{\bm{T}},\tilde{\bm{X}})$, we decompose the parts other than
$\bm{e}_{j}^{\prime}(\tilde{\bm{X}}_{\bm{M}}^{\prime}\tilde{\bm{X}}_{\bm{M}})^{-1}$
and evaluate their expectations. Then, we obtain

$\displaystyle{\rm
E}\bigg{[}\bm{X}_{\bm{M},i}\sum_{h=1}^{H}c^{(h)}\frac{T_{i}^{(h)}}{e^{(h)}(\bm{X}_{i})}\\{\mu^{(h)}(\bm{X}_{i})+f(\bm{X}_{i})\\}\
\bigg{|}\ \tilde{\bm{X}}\bigg{]}$
$\displaystyle=\bm{X}_{\bm{M},i}\sum_{h=1}^{H}c^{(h)}\\{\mu^{(h)}(\bm{X}_{i})+f(\bm{X}_{i})\\}$
$\displaystyle=\bm{X}_{\bm{M},i}\sum_{h=1}^{H}c^{(h)}\mu^{(h)}(\bm{X}_{i}),$

since the sum of contrasts is 0, where
$\bm{X}_{\bm{M},i}=(X_{ij})_{j\in\bm{M}}$.

Thus, we have

$\displaystyle\kappa_{j}^{\bm{M}}(\tilde{\bm{T}},\tilde{\bm{X}})-\beta_{j}^{\bm{M}}$
$\displaystyle=\bm{e}_{j}^{\prime}(\tilde{\bm{X}}_{\bm{M}}^{\prime}\tilde{\bm{X}}_{\bm{M}})^{-1}\tilde{\bm{X}}_{\bm{M}}^{\prime}\sum_{h=1}^{H}c^{(h)}\\{\tilde{\bm{W}}^{(h)}(\tilde{\bm{T}}^{(h)},\tilde{\bm{X}})-\bm{I}_{n}\\}\\{\tilde{\bm{\mu}}^{(h)}(\tilde{\bm{X}})+\tilde{\bm{f}}(\tilde{\bm{X}})\\}={\rm
o}_{\rm P}(1)$

from the law of large numbers. We prepare an estimator $Y_{i}^{(h)\dagger}$
for $\mu^{(h)}(\bm{x}_{i})+f(\bm{x}_{i})$, and let
$\tilde{\bm{Y}}^{(h)\dagger}=(Y_{1}^{(h)\dagger},\ldots,Y_{n}^{(h)\dagger})^{\prime}$
and
$\tilde{\bm{Y}}^{\dagger}=(\tilde{\bm{Y}}^{(1)\dagger}{}^{\prime},\ldots,\tilde{\bm{Y}}^{(H)\dagger}{}^{\prime})^{\prime}$.
Accordingly, we define

$\displaystyle\tau_{j}^{\bm{M}}(\tilde{\bm{Y}}^{\dagger},\tilde{\bm{T}},\tilde{\bm{X}})\equiv\bm{e}_{j}^{\prime}(\tilde{\bm{X}}_{\bm{M}}^{\prime}\tilde{\bm{X}}_{\bm{M}})^{-1}\tilde{\bm{X}}_{\bm{M}}^{\prime}\sum_{h=1}^{H}c^{(h)}\\{\tilde{\bm{W}}^{(h)}(\tilde{\bm{T}}^{(h)},\tilde{\bm{X}})-\bm{I}_{n}\\}\tilde{\bm{Y}}^{(h)\dagger}.$
(15)

This leads us to the following lemma.

######  Lemma 1.

Let $F_{\mu,\sigma^{2}}^{[a,b]}(\cdot)$ be the cumulative distribution
functions of ${\rm TN}(\mu,\sigma^{2},a,b)$, which is ${\rm
N}(\mu,\sigma^{2})$ truncated into the interval $[a,b]$, and let
${\mathcal{V}}_{\bm{s},j}^{-}$, ${\mathcal{V}}_{\bm{s},j}^{+}$,
$\tau_{j}^{\bm{M}}$ and $\zeta_{j}^{\bm{M}}$ be the functions defined in (10),
(11), (15) and (14), respectively. Then, if
$Y_{i}^{(h)\dagger}-\\{\mu^{(h)}(\bm{X}_{i})+f(\bm{X}_{i})\\}={\rm o}_{\rm
P}(1)$ under the condition,

$\displaystyle[F_{\beta_{j}^{\bm{M}}+\tau_{j}^{\bm{M}}(\tilde{\bm{Y}}^{\dagger},\tilde{\bm{T}},\tilde{\bm{X}}),\zeta_{j}^{\bm{M}}(\tilde{\bm{T}},\tilde{\bm{X}})}^{[{\mathcal{V}}_{\bm{s},j}^{-}(\tilde{\bm{Z}},\tilde{\bm{T}},\tilde{\bm{X}}),{\mathcal{V}}_{\bm{s},j}^{+}(\tilde{\bm{Z}},\tilde{\bm{T}},\tilde{\bm{X}})]}(\tilde{\bm{\eta}}_{j}^{\prime}\tilde{\bm{W}}(\tilde{\bm{T}},\tilde{\bm{X}})\tilde{\bm{Y}})\mid\bm{A}(\bm{M},\bm{s})\tilde{\bm{W}}(\tilde{\bm{T}},\tilde{\bm{X}})\tilde{\bm{Y}}\leq\bm{b}(\bm{M},\bm{s})]\stackrel{{\scriptstyle\rm
d}}{{\to}}{\rm Unif}(0,1)$

holds as $n\to\infty$.

Note that although this lemma makes use of asymptotics unlike Theorem 1, it is
essentially different from those of Tian and Taylor (2017) and Tibshirani et
al. (2018). If we try to use the results of these papers simply because the
error in the context of causal inference, $f(\bm{X})+\epsilon^{(h)}$, has a
non-Gaussian distribution, the estimation will have a large bias due to the
existence of confounding. On the other hand, in this paper, we assume a
Gaussian distribution for $\epsilon^{(h)}$, and this lemma does not mean that
our result is valid even if we make it non-Gaussian.

Similarly to Lee et al. (2016), if we condition on (9), we only have to use
the normal distribution ${\rm N}(\mu,\sigma^{2})$ truncated into the union set
$S$ of intervals. Writing its cumulative distribution function as
$F^{S}_{\mu,\sigma^{2}}(\cdot)$, we get

$\displaystyle\bigg{[}F_{\beta_{j}^{\bm{M}}+\tau_{j}^{\bm{M}}(\tilde{\bm{Y}}^{\dagger},\tilde{\bm{T}},\tilde{\bm{X}}),\zeta_{j}^{\bm{M}}(\tilde{\bm{T}},\tilde{\bm{X}})}^{\bigcup_{\bm{s}}[{\mathcal{V}}_{\bm{s},j}^{-}(\tilde{\bm{Z}},\tilde{\bm{T}},\tilde{\bm{X}}),{\mathcal{V}}_{\bm{s},j}^{+}(\tilde{\bm{Z}},\tilde{\bm{T}},\tilde{\bm{X}})]}(\tilde{\bm{\eta}}_{j}^{\prime}\tilde{\bm{W}}(\tilde{\bm{T}},\tilde{\bm{X}})\tilde{\bm{Y}})\
\bigg{|}\
\bigcup_{\bm{s}}\\{\bm{A}(\bm{M},\bm{s})\tilde{\bm{W}}(\tilde{\bm{T}},\tilde{\bm{X}})\tilde{\bm{Y}}\leq\bm{b}(\bm{M},\bm{s})\\}\bigg{]}$
$\displaystyle\stackrel{{\scriptstyle\rm d}}{{\to}}{\rm Unif}(0,1).$

In addition, because
$F_{\beta_{j}^{\bm{M}}+\tau_{j}^{\bm{M}}(\tilde{\bm{Y}}^{\dagger},\tilde{\bm{T}},\tilde{\bm{X}}),\zeta_{j}^{\bm{M}}(\tilde{\bm{T}},\tilde{\bm{X}})}^{\bigcup_{\bm{s}}[{\mathcal{V}}_{\bm{s},j}^{-}(\tilde{\bm{Z}},\tilde{\bm{T}},\tilde{\bm{X}}),{\mathcal{V}}_{\bm{s},j}^{+}(\tilde{\bm{Z}},\tilde{\bm{T}},\tilde{\bm{X}})]}(\tilde{\bm{\eta}}_{j}^{\prime}\tilde{\bm{W}}(\tilde{\bm{T}},\tilde{\bm{X}})\tilde{\bm{Y}})$
is a monotonically decreasing function with respect to $\beta_{j}^{\bm{M}}$,
we can see that if we set $L$ and $U$ to satisfy

$\displaystyle
F_{L+\tau_{j}^{\bm{M}}(\tilde{\bm{Y}}^{\dagger},\tilde{\bm{T}},\tilde{\bm{X}}),\zeta_{j}^{\bm{M}}(\tilde{\bm{T}},\tilde{\bm{X}})}^{\bigcup_{\bm{s}}[{\mathcal{V}}_{\bm{s},j}^{-}(\tilde{\bm{Z}},\tilde{\bm{T}},\tilde{\bm{X}}),{\mathcal{V}}_{\bm{s},j}^{+}(\tilde{\bm{Z}},\tilde{\bm{T}},\tilde{\bm{X}})]}(\tilde{\bm{\eta}}_{j}^{\prime}\tilde{\bm{W}}(\tilde{\bm{T}},\tilde{\bm{X}})\tilde{\bm{Y}})=1-\frac{\alpha}{2}$

and

$\displaystyle
F_{U+\tau_{j}^{\bm{M}}(\tilde{\bm{Y}}^{\dagger},\tilde{\bm{T}},\tilde{\bm{X}}),\zeta_{j}^{\bm{M}}(\tilde{\bm{T}},\tilde{\bm{X}})}^{\bigcup_{\bm{s}}[{\mathcal{V}}_{\bm{s},j}^{-}(\tilde{\bm{Z}},\tilde{\bm{T}},\tilde{\bm{X}}),{\mathcal{V}}_{\bm{s},j}^{+}(\tilde{\bm{Z}},\tilde{\bm{T}},\tilde{\bm{X}})]}(\tilde{\bm{\eta}}_{j}^{\prime}\tilde{\bm{W}}(\tilde{\bm{T}},\tilde{\bm{X}})\tilde{\bm{Y}})=\frac{\alpha}{2},$

we get

$\displaystyle{\rm P}(\beta_{j}^{\bm{M}}\in[L,U]\mid\hat{\bm{M}}=\bm{M})\to
1-\alpha,$

which gives the asymptotically guaranteed conditional coverage.

Next, we develop a theory using a concrete and easily conceivable estimator
for $\mu^{(h)}(\bm{X}_{i})+f(\bm{X}_{i})$. Suppose we have a sequence of real
numbers $\\{\delta_{n}\\}$ that converges to 0, and let
${\mathcal{N}}_{i}^{(h)}\equiv\\{l\neq
i:\|\bm{X}_{l}-\bm{X}_{i}\|_{2}\leq\delta_{n},\ T_{l}^{(h)}=1\\}$. We take
$\delta_{n}$ such that ${\rm O}(1)\neq{\mathcal{N}}_{i}^{(h)}={\rm o}(n)$. For
example, if the set of possible values of $\bm{X}_{i}$ is a closed set of
bounded open sets in $\mathbb{R}^{p}$, then we only have to take $\delta_{n}$
such that ${\rm O}(n^{-1})\neq\delta_{n}={\rm o}(1)$. Letting
$\hat{\bm{\beta}}^{(h)}\equiv\mathop{\rm
argmin}_{\bm{\beta}}\|\tilde{\bm{T}}^{(h)}(\tilde{\bm{Y}}-\tilde{\bm{X}}\bm{\beta})\|_{2}^{2}$
and
$\tilde{\bm{Y}}^{(h)*}\equiv(\bm{X}_{1}^{\prime}\hat{\bm{\beta}}^{(h)}+\sum_{l\in{\mathcal{N}}_{1}^{(h)}}(Y_{l}-\bm{X}_{l}^{\prime}\hat{\bm{\beta}}^{(h)})/|{\mathcal{N}}_{1}^{(h)}|,\ldots,\bm{X}_{n}^{\prime}\hat{\bm{\beta}}^{(h)}+\sum_{l\in{\mathcal{N}}_{n}^{(h)}}(Y_{l}-\bm{X}_{l}^{\prime}\hat{\bm{\beta}}^{(h)})/|{\mathcal{N}}_{n}^{(h)}|)^{\prime}$,
if $f(\cdot)$ is Lipschitz continuous, then it is obvious that
$\kappa_{j}^{\bm{M}}(\tilde{\bm{T}},\tilde{\bm{X}})-\\{\beta_{j}^{\bm{M}}+\tau_{j}^{\bm{M}}(\tilde{\bm{Y}}^{*},\tilde{\bm{T}},\tilde{\bm{X}})\\}={\rm
o}_{\rm P}(n^{-1/2})$, where $\tau_{j}^{\bm{M}}$ is the function defined in
(15) and
$\tilde{\bm{Y}}^{*}=(\tilde{\bm{Y}}^{(1)*}{}^{\prime},\ldots,\tilde{\bm{Y}}^{(H)*}{}^{\prime})^{\prime}$.
Thus, we obtain

$\displaystyle{\rm
E}\\{\tilde{\bm{\eta}}_{j}^{\prime}\tilde{\bm{W}}(\tilde{\bm{T}},\tilde{\bm{X}})\tilde{\bm{Y}}-\tau_{j}^{\bm{M}}(\tilde{\bm{Y}}^{*},\tilde{\bm{T}},\tilde{\bm{X}})\mid\tilde{\bm{T}},\tilde{\bm{X}}\\}=\beta_{j}^{\bm{M}}+{\rm
o}_{\rm P}(n^{-1/2}).$ (16)

For the evaluation of the variance,

$\displaystyle{\rm
V}\\{\tilde{\bm{\eta}}_{j}^{\prime}\tilde{\bm{W}}(\tilde{\bm{T}},\tilde{\bm{X}})\tilde{\bm{Y}}-\tau_{j}^{\bm{M}}(\tilde{\bm{Y}}^{*},\tilde{\bm{T}},\tilde{\bm{X}})\mid\tilde{\bm{T}},\tilde{\bm{X}}\\}$
$\displaystyle={\rm
E}([\\{\tilde{\bm{\eta}}_{j}^{\prime}\tilde{\bm{W}}(\tilde{\bm{T}},\tilde{\bm{X}})\tilde{\bm{Y}}-\kappa_{j}^{\bm{M}}(\tilde{\bm{T}},\tilde{\bm{X}})\\}-\tau_{j}^{\bm{M}}(\tilde{\bm{\epsilon}}^{*},\tilde{\bm{T}},\tilde{\bm{X}})]^{2}\mid\tilde{\bm{T}},\tilde{\bm{X}}),$

the conditional expectation of the square of the first term is
$\zeta_{j}^{\bm{M}}(\tilde{\bm{T}},\tilde{\bm{X}})$, the conditional
expectation of the square of the second term is

$\displaystyle\tilde{\bm{\eta}}_{j}^{\prime}\sum_{h=1}^{H}\sum_{k=1}^{H}c^{(h)}c^{(k)}{\rm
E}[\\{\tilde{\bm{W}}^{(h)}(\tilde{\bm{T}}^{(h)},\tilde{\bm{X}})-\bm{I}_{n}\\}\tilde{\bm{\epsilon}}^{(h)*}\tilde{\bm{\epsilon}}^{(k)*}{}^{\prime}\\{\tilde{\bm{W}}^{(k)}(\tilde{\bm{T}}^{(k)},\tilde{\bm{X}})-\bm{I}_{n}\\}\mid\tilde{\bm{T}},\tilde{\bm{X}}]\tilde{\bm{\eta}}_{j},$
(17)

and the conditional expectation of the product of the first and second terms
is

$\displaystyle\tilde{\bm{\eta}}_{j}^{\prime}\sum_{h=1}^{H}\sum_{k=1}^{H}c^{(h)}c^{(k)}{\rm
E}[\tilde{\bm{W}}^{(h)}(\tilde{\bm{T}}^{(h)},\tilde{\bm{X}})\tilde{\bm{\epsilon}}^{(h)}\tilde{\bm{\epsilon}}^{(k)*}{}^{\prime}\\{\tilde{\bm{W}}^{(k)}(\tilde{\bm{T}}^{(k)},\tilde{\bm{X}})-\bm{I}_{n}\\}\mid\tilde{\bm{T}},\tilde{\bm{X}}]\tilde{\bm{\eta}}_{j},$
(18)

where $\tilde{\bm{\epsilon}}^{*}=\tilde{\bm{Y}}^{*}-{\rm
E}(\tilde{\bm{Y}}^{*}\mid\tilde{\bm{T}},\tilde{\bm{X}})$.

In (17), if $h\neq k$, then $\tilde{\bm{\epsilon}}^{(h)*}$ and
$\tilde{\bm{\epsilon}}^{(k)*}$ are independent, and if $h=k$, then
$\tilde{\bm{\epsilon}}^{(h)*}\tilde{\bm{\epsilon}}^{(k)*}{}^{\prime}$ has at
most $\sum_{i=1}^{n}|{\mathcal{N}}_{i}^{(h)}|^{2}$ non-zero components. Given
the expectation with respect to $\tilde{\bm{T}}$ of the conditional
expectation, the off-diagonal components of this matrix are 0 due to the
weakly ignorable treatment assignment condition and the independence of
$\\{T_{i}^{(h)}\\}$; moreover, the diagonal components are
$(\sigma^{2}/|{\mathcal{N}}_{i}^{(h)}|)\\{1-e^{(h)}(\bm{X}_{i})\\}/e^{(h)}(\bm{X}_{i})$
from the assignment condition. Hence, we can see that (17) becomes

$\displaystyle\sigma^{2}\sum_{h=1}^{H}c^{(h)2}\tilde{\bm{\eta}}_{j}^{\prime}\tilde{\bm{K}}^{(h)}(\tilde{\bm{X}})\tilde{\bm{\eta}}_{j}+{\rm
o}_{\rm P}(n^{-3/2}),$ (19)

where $\tilde{\bm{K}}^{(h)}(\tilde{\bm{X}})\equiv{\rm
diag}[\\{1/e^{(h)}(\bm{X}_{1})-1\\}/|{\mathcal{N}}_{1}^{(h)}|,\ldots,\\{1/e^{(h)}(\bm{X}_{n})-1\\}/|{\mathcal{N}}_{n}^{(h)}|]$.
In (18), if $h\neq k$, then $\tilde{\bm{\epsilon}}^{(h)}$ and
$\tilde{\bm{\epsilon}}^{(k)*}$ are independent, and if $h=k$, then
$\tilde{\bm{\epsilon}}^{(h)}\tilde{\bm{\epsilon}}^{(k)*}{}^{\prime}$ has at
most $\sum_{i=1}^{n}|{\mathcal{N}}_{i}^{(h)}|$ non-zero components. Given the
expectation with respect to $\tilde{\bm{T}}$ of the conditional expectation,
the off-diagonal components of this matrix are 0 due to the weakly ignorable
treatment assignment condition and the independence of $\\{T_{i}^{(h)}\\}$.
Moreover, the diagonal components are also 0 because of the assignment
condition and the independence of $\epsilon_{i}^{(h)}$ and
$\epsilon_{i}^{(h)*}$. Hence, we can see that (18) becomes ${\rm o}_{\rm
P}(n^{-3/2})$.

Thus, by defining

$\displaystyle\rho_{j}^{\bm{M}}(\tilde{\bm{T}},\tilde{\bm{X}})\equiv\zeta_{j}^{\bm{M}}(\tilde{\bm{T}},\tilde{\bm{X}})+\sigma^{2}\sum_{h=1}^{H}c^{(h)2}\tilde{\bm{\eta}}_{j}^{\prime}\tilde{\bm{K}}^{(h)}(\tilde{\bm{X}})\tilde{\bm{\eta}}_{j},$
(20)

it can be seen that

$\displaystyle{\rm
V}\\{\tilde{\bm{\eta}}_{j}^{\prime}\tilde{\bm{W}}(\tilde{\bm{T}},\tilde{\bm{X}})\tilde{\bm{Y}}-\tau_{j}^{\bm{M}}(\tilde{\bm{Y}}^{*},\tilde{\bm{T}},\tilde{\bm{X}})\mid\tilde{\bm{T}},\tilde{\bm{X}}\\}=\rho_{j}^{\bm{M}}(\tilde{\bm{T}},\tilde{\bm{X}})+{\rm
o}_{\rm P}(n^{-3/2}).$ (21)

From (16) and (21), we obtain the following theorem, which is an asymptotic
version of (5).

######  Theorem 2.

Let $F_{\mu,\sigma^{2}}^{[a,b]}(\cdot)$ be the cumulative distribution
functions of ${\rm TN}(\mu,\sigma^{2},a,b)$, which is ${\rm
N}(\mu,\sigma^{2})$ truncated into the interval $[a,b]$, and let
${\mathcal{V}}_{\bm{s},j}^{-}$, ${\mathcal{V}}_{\bm{s},j}^{+}$,
$\tau_{j}^{\bm{M}}$ and $\rho_{j}^{\bm{M}}$ be the functions defined in (10),
(11), (15) and (20), respectively. Then, if we define $L_{j}^{\hat{\bm{M}}}$
and $U_{j}^{\hat{\bm{M}}}$ such that

$\displaystyle
F_{L_{j}^{\hat{\bm{M}}}+\tau_{j}^{\hat{\bm{M}}}(\tilde{\bm{Y}}^{*},\tilde{\bm{T}},\tilde{\bm{X}}),\rho_{j}^{\hat{\bm{M}}}(\tilde{\bm{T}},\tilde{\bm{X}})}^{\bigcup_{\bm{s}}[{\mathcal{V}}_{\bm{s},j}^{-}(\tilde{\bm{Z}},\tilde{\bm{T}},\tilde{\bm{X}}),{\mathcal{V}}_{\bm{s},j}^{+}(\tilde{\bm{Z}},\tilde{\bm{T}},\tilde{\bm{X}})]}(\tilde{\bm{\eta}}_{j}^{\prime}\tilde{\bm{W}}(\tilde{\bm{T}},\tilde{\bm{X}})\tilde{\bm{Y}})=1-\frac{\alpha}{2}$

and

$\displaystyle
F_{U_{j}^{\hat{\bm{M}}}+\tau_{j}^{\hat{\bm{M}}}(\tilde{\bm{Y}}^{*},\tilde{\bm{T}},\tilde{\bm{X}}),\rho_{j}^{\hat{\bm{M}}}(\tilde{\bm{T}},\tilde{\bm{X}})}^{\bigcup_{\bm{s}}[{\mathcal{V}}_{\bm{s},j}^{-}(\tilde{\bm{Z}},\tilde{\bm{T}},\tilde{\bm{X}}),{\mathcal{V}}_{\bm{s},j}^{+}(\tilde{\bm{Z}},\tilde{\bm{T}},\tilde{\bm{X}})]}(\tilde{\bm{\eta}}_{j}^{\prime}\tilde{\bm{W}}(\tilde{\bm{T}},\tilde{\bm{X}})\tilde{\bm{Y}})=\frac{\alpha}{2},$

then

$\displaystyle{\rm
E}\bigg{(}\frac{|\\{j\in\hat{\bm{M}}:\beta_{j}^{\hat{\bm{M}}}\notin[L_{j}^{\hat{\bm{M}}},U_{j}^{\hat{\bm{M}}}]\\}|}{|\hat{\bm{M}}|};\
|\hat{\bm{M}}|>0\bigg{)}\leq\alpha$

holds as $n\to\infty$.

The key point here is that the pivot statistic does not depend on the
nonparametric function $\tilde{\bm{f}}(\tilde{\bm{X}})$.

######  Remark.

Regarding $\tilde{\bm{Y}}^{*}$, for example, when $H=2$ and
$(c^{(1)},c^{(2)})=(-1,1)$, if we define
$\tilde{\bm{Y}}^{(1)*}\equiv(\bm{X}_{1}^{\prime}\hat{\bm{\beta}}^{(1)}+\\{\sum_{l\in{\mathcal{N}}_{1}^{(1)}}(Y_{l}-\bm{X}_{l}^{\prime}\hat{\bm{\beta}}^{(1)})+\sum_{l\in{\mathcal{N}}_{1}^{(2)}}(Y_{l}-\bm{X}_{\hat{\bm{M}},l}^{\prime}\hat{\bm{\beta}}^{\hat{\bm{M}}}-\bm{X}_{l}^{\prime}\hat{\bm{\beta}}^{(1)})\\}/(|{\mathcal{N}}_{1}^{(1)}|+|{\mathcal{N}}_{1}^{(2)}|),\ldots,\bm{X}_{n}^{\prime}\hat{\bm{\beta}}^{(1)}+\\{\sum_{l\in{\mathcal{N}}_{n}^{(1)}}(Y_{l}-\bm{X}_{l}^{\prime}\hat{\bm{\beta}}^{(1)})+\sum_{l\in{\mathcal{N}}_{n}^{(2)}}(Y_{l}-\bm{X}_{\hat{\bm{M}},l}^{\prime}\hat{\bm{\beta}}^{\hat{\bm{M}}}-\bm{X}_{l}^{\prime}\hat{\bm{\beta}}^{(1)})\\}/(|{\mathcal{N}}_{n}^{(1)}|+|{\mathcal{N}}_{n}^{(2)}|))^{\prime}$
and
$\tilde{\bm{Y}}^{(2)*}\equiv(\bm{X}_{1}^{\prime}\hat{\bm{\beta}}^{(2)}+\\{\sum_{l\in{\mathcal{N}}_{1}^{(1)}}(Y_{l}+\bm{X}_{\hat{\bm{M}},l}^{\prime}\hat{\bm{\beta}}^{\hat{\bm{M}}}-\bm{X}_{l}^{\prime}\hat{\bm{\beta}}^{(2)})+\sum_{l\in{\mathcal{N}}_{1}^{(2)}}(Y_{l}-\bm{X}_{l}^{\prime}\hat{\bm{\beta}}^{(2)})\\}/(|{\mathcal{N}}_{1}^{(1)}|+|{\mathcal{N}}_{1}^{(2)}|),\ldots,\bm{X}_{n}^{\prime}\hat{\bm{\beta}}^{(2)}+\\{\sum_{l\in{\mathcal{N}}_{n}^{(1)}}(Y_{l}+\bm{X}_{\hat{\bm{M}},l}^{\prime}\hat{\bm{\beta}}^{\hat{\bm{M}}}-\bm{X}_{l}^{\prime}\hat{\bm{\beta}}^{(2)})+\sum_{l\in{\mathcal{N}}_{n}^{(2)}}(Y_{l}-\bm{X}_{l}^{\prime}\hat{\bm{\beta}}^{(2)})\\}/(|{\mathcal{N}}_{n}^{(1)}|+|{\mathcal{N}}_{n}^{(2)}|))^{\prime}$,
the accuracy tends to increase. Note that in the definition of
$\tilde{\bm{K}}^{(h)}(\tilde{\bm{X}})$, $|{\mathcal{N}}_{i}^{(h)}|$ becomes
$|{\mathcal{N}}_{i}^{(1)}|+|{\mathcal{N}}_{i}^{(2)}|$.

##  4 Numerical experiment

We numerically verified the usefulness of the proposed method. Letting $H=2$,
for the confounding variable $\bm{X}\ (\in\mathbb{R}^{p})$, the outcome
variable $Y\ (\in\mathbb{R})$ was observed through the model,

$\displaystyle Y=T\mu(\bm{X})+f(\bm{X})+\epsilon.$

That is, the treatment group was represented by $T=1$ and the control group by
$T=0$. The error variable $\epsilon$ was assumed to follow a normal
distribution ${\rm N}(0,\sigma^{2})$ independently of $\bm{X}$, and the
variance was assumed to be known. This model is a special case of (1), where
the causal effect is $\mu(\bm{x})$.

The setup of the numerical experiment is as follows. In each run, we choose
one of the following

(F1)

$f(\bm{x})=0$

(F2)

$f(\bm{x})=3x_{1}+x_{2}+x_{3}+x_{4}+x_{5}-3.5$

(F3)

$f(\bm{x})=x_{1}+0.5x_{2}+0.5x_{3}+\pi\sin(\pi x_{4})/32+\pi\sin(\pi
x_{5})/32-1.125$

as an infinite-dimensional nuisance parameter $f(\cdot)$. In addition, we
chose one of the following

(E1)

$e(\bm{x})=0.5$

(E2)

$e(\bm{x})=1/\\{1+\exp(-x_{1}-0.5x_{2}-0.5x_{3}-0.5x_{4}-0.5x_{5}+1.5)\\}$

as the propensity score $e(\cdot)$ and one of the following

(M1)

$\mu(\bm{x})=0$

(M2)

$\mu(\bm{x})=3x_{1}+x_{2}+x_{3}+x_{4}+x_{5}-3.5$

as the causal effect $\mu(\cdot)$. Note that (E1) indicates a completely
random assignment. In the case of (M1), all variables are inactive for
$\mu(\bm{x})$, while in the case of (M2), the variables
$\\{x_{1},x_{2},x_{3},x_{4},x_{5}\\}$ are active for $\mu(\bm{x})$. Each
component of the confounding variable vector independently follows a
continuous uniform distribution ${\rm U}(0,1)$ or Bernoulli distribution ${\rm
Ber}(1/2)$, and the data are generated from the above model with a sample size
of $n=1000$. Here, we set $p$ to $5$ or $25$ and the error variance
$\sigma^{2}$ to $0.25^{2}$.

We compared the proposed method (SI) with a method that ignores the influence
of model selection (Naive). Model selection in Naive is conducted by using
Lasso and then confidence intervals are constructed using a normal
distribution instead of a truncated normal distribution. Because Naive ignores
the fact that the selected variables are likely to be significant, it is
supposed that it will have a large number of false positives and the validity
of its inference will be impaired.

We evaluated the size $|\hat{M}|$ of the model selected by Lasso and the false
coverage rate (FCR) by conducting 1,000 Monte Carlo simulations. We also
evaluated the probability that $\beta_{1},\ldots,\beta_{5}$ is significant,
i.e., the probability that $\beta_{1}=0,\ldots,\beta_{5}=0$ is not covered in
the confidence interval when $p=5$, and numbers of true positive (TP) and
false positive (FP) when $p=25$. In all cases, the significance level was set
to 0.05; i.e., the confidence coefficient was set to 0.95. Therefore, the
closer the FCR is to 0.05, the more valid the method is. Note that with
$\hat{M}$ as the selected model and $C_{j}^{\hat{\bm{M}}}$ as the confidence
interval for $\beta_{j}^{\hat{\bm{M}}}$, TP is the number of intervals that
did not cover zero when $\beta_{j}^{\hat{\bm{M}}}$ is non-zero,
$|\\{j\in\hat{\bm{M}}:\beta_{j}^{\hat{\bm{M}}}\neq 0,\ 0\notin
C_{j}^{\hat{\bm{M}}}\\}|$, and FP is the number intervals that did not cover
zero when $\beta_{j}^{\hat{\bm{M}}}$ is zero,
$|\\{j\in\hat{\bm{M}}:\beta_{j}^{\hat{\bm{M}}}=0,\ 0\notin
C_{j}^{\hat{\bm{M}}}\\}|$.

###  4.1 Case of $p=5$

First, let us consider (M1), where all variables are inactive. The confounding
variables followed a Bernoulli distribution, and the sequence $\delta_{n}$
used to define the neighborhood ${\cal N}_{i}^{(h)}$ was set to $0$. The
tuning parameter of Lasso was $\lambda=\sigma n^{-1/2}(\log p)^{1/2}$ when the
infinite-dimensional nuisance parameter was (F1), and $\lambda=2\sigma
n^{-1/2}(\log p)^{1/2}$ when it was (F2) or (F3). Table 1 shows $|\hat{M}|$,
the non-coverage probability, and the FCR. Since $\beta_{j}=0$ is correct for
all variables, the closer the non-coverage probability is to 0.05, the better.
Although the non-coverage probability of SI seems to be somewhat large when
the propensity score and the nuisance parameter are (E2) and (F2), it is
approximately 0.05 in most other cases. On the other hand, if we look at
Naive, we can see that when the nuisance parameter is (F2) or (F3), 0 is not
covered for all the selected variables, and the validity of the inference is
impaired. Note that when the nuisance parameter is (F1), the non-coverage
probability is greater than 0.05, although this value is not as extreme as
when the nuisance parameter is (F2) or (F3). The FCR shows a similar trend,
where SI keeps the FCR at almost 0.05, but Naive results in a very large
value.

Table 1: Case in which $p=5$, $X_{ij}\sim{\rm Ber}(1/2)$ and (M1) is correct. The figures in parentheses are standard deviations. In all settings, SI controls the FCR appropriately, but the inference of Naive is invalid.  $e(\bm{x})$ | $f(\bm{x})$ | $|\hat{M}|$ |  | $\beta_{1}$ | $\beta_{2}$ | $\beta_{3}$ | $\beta_{4}$ | $\beta_{5}$ | FCR
---|---|---|---|---|---|---|---|---|---
|  |  | SI | 0.039 | 0.031 | 0.026 | 0.045 | 0.047 | 0.035
| (F1) | 2.593 (1.164) | (0.194) | (0.175) | (0.159) | (0.208) | (0.212) | (0.124)
| Naive | 0.080 | 0.085 | 0.070 | 0.089 | 0.088 | 0.079
|  |  | (0.272) | (0.279) | (0.255) | (0.285) | (0.284) | (0.178)
|  |  | SI | 0.049 | 0.055 | 0.038 | 0.050 | 0.059 | 0.051
(E1) | (F2) | 3.641 (1.012) | (0.215) | (0.227) | (0.190) | (0.218) | (0.235) | (0.123)
Naive | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 0.999
|  |  | (0.000) | (0.000) | (0.000) | (0.000) | (0.000) | (0.032)
|  |  | SI | 0.040 | 0.051 | 0.042 | 0.043 | 0.048 | 0.039
| (F3) | 1.675 (1.141) | (0.196) | (0.219) | (0.201) | (0.204) | (0.213) | (0.156)
| Naive | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 0.842
|  |  | (0.000) | (0.000) | (0.000) | (0.000) | (0.000) | (0.365)
|  |  | SI | 0.043 | 0.049 | 0.040 | 0.050 | 0.044 | 0.042
| (F1) | 2.734 (1.154) | (0.202) | (0.216) | (0.196) | (0.217) | (0.204) | (0.134)
| Naive | 0.081 | 0.076 | 0.080 | 0.064 | 0.070 | 0.068
|  |  | (0.274) | (0.265) | (0.272) | (0.246) | (0.255) | (0.157)
|  |  | SI | 0.065 | 0.067 | 0.063 | 0.050 | 0.079 | 0.062
(E2) | (F2) | 3.744 (1.004) | (0.246) | (0.250) | (0.244) | (0.218) | (0.270) | (0.153)
Naive | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 0.999
|  |  | (0.000) | (0.000) | (0.000) | (0.000) | (0.000) | (0.032)
|  |  | SI | 0.045 | 0.066 | 0.056 | 0.053 | 0.050 | 0.042
| (F3) | 1.916 (1.191) | (0.207) | (0.248) | (0.230) | (0.224) | (0.219) | (0.145)
| Naive | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 0.886
|  |  | (0.000) | (0.000) | (0.000) | (0.000) | (0.000) | (0.318)

Table 2 shows results for when (M2) is the causal effect. Here, the
confounding variables followed a continuous uniform distribution,
$\delta_{n}=(p/6)^{1/2}$ and $\lambda=5\sigma n^{-1/2}(\log p)^{1/2}$. In this
case, all variables are active, so the closer the non-coverage probability is
to 1, the better. The non-coverage probability of both methods is close to 1,
but Naive does not keep the FCR at 0.05. On the other hand, SI keeps it at
almost 0.05.

Table 2: Case in which $p=5$, $X_{ij}\sim{\rm U}(0,1)$ and (M2) is correct. The figures in parentheses are standard deviations. In all settings, SI removes a high percentage of 0’s from the confidence intervals while appropriately controlling the FCR.  $e(\bm{x})$ | $f(\bm{x})$ | $|\hat{M}|$ |  | $\beta_{1}$ | $\beta_{2}$ | $\beta_{3}$ | $\beta_{4}$ | $\beta_{5}$ | FCR
---|---|---|---|---|---|---|---|---|---
|  |  | SI | 1.000 | 0.998 | 0.998 | 0.999 | 0.999 | 0.054
| (F1) | 4.987 (0.113) | (0.000) | (0.045) | (0.045) | (0.032) | (0.032) | (0.099)
| Naive | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 0.417
|  |  | (0.000) | (0.000) | (0.000) | (0.000) | (0.000) | (0.233)
|  |  | SI | 1.000 | 0.994 | 0.995 | 0.992 | 0.990 | 0.067
(E1) | (F2) | 4.460 (0.754) | (0.000) | (0.076) | (0.068) | (0.089) | (0.101) | (0.127)
Naive | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 0.737
|  |  | (0.000) | (0.000) | (0.000) | (0.000) | (0.000) | (0.215)
|  |  | SI | 1.000 | 0.998 | 0.991 | 0.999 | 0.997 | 0.064
| (F3) | 4.888 (0.334) | (0.000) | (0.045) | (0.096) | (0.032) | (0.055) | (0.114)
| Naive | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 0.591
|  |  | (0.000) | (0.000) | (0.000) | (0.000) | (0.000) | (0.230)
|  |  | SI | 1.000 | 1.000 | 0.998 | 0.999 | 0.999 | 0.054
| (F1) | 4.980 (0.140) | (0.000) | (0.000) | (0.045) | (0.032) | (0.032) | (0.101)
| Naive | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 0.458
|  |  | (0.000) | (0.000) | (0.000) | (0.000) | (0.000) | (0.227)
|  |  | SI | 1.000 | 0.992 | 0.987 | 0.993 | 0.991 | 0.066
(E2) | (F2) | 4.445 (0.778) | (0.000) | (0.089) | (0.112) | (0.083) | (0.097) | (0.125)
Naive | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 0.749
|  |  | (0.000) | (0.000) | (0.000) | (0.000) | (0.000) | (0.221)
|  |  | SI | 1.000 | 0.998 | 0.998 | 0.997 | 0.996 | 0.057
| (F3) | 4.864 (0.384) | (0.000) | (0.045) | (0.046) | (0.056) | (0.064) | (0.108)
| Naive | 1.000 | 1.000 | 1.000 | 1.000 | 1.000 | 0.614
|  |  | (0.000) | (0.000) | (0.000) | (0.000) | (0.000) | (0.228)

Figure 1 compares the confidence intervals obtained by SI and Naive. In all
panels, the nuisance parameter is (F3). Figures 1LABEL:sub@f:1a and
1LABEL:sub@f:1b respectively use (E1) and (E2) as the propensity scores in the
setting of Table 1. On the other hand, Figures 1LABEL:sub@f:1c and
1LABEL:sub@f:1d respectively use (E1) and (E2) as the propensity scores in the
setting of Table 2. Note that the horizontal axis represents the variables
selected by Lasso, and in Figure 1LABEL:sub@f:1a, the variables
$x_{2},x_{3},x_{4},x_{5}$ are selected, while in the remaining figures, all
variables are selected. Since (M1) is correct in Figures 1LABEL:sub@f:1a and
1LABEL:sub@f:1b, the confidence intervals should contain 0. In both panels, it
can be seen that most of Naive’s results do not contain 0, but those of SI do
contain 0. In Figure 1LABEL:sub@f:1a, the two confidence intervals for SI do
not include the estimates, which may seem somewhat strange, but a similar
phenomenon has been reported in Liu et al. (2018). Figure 1LABEL:sub@f:1c
shows that the confidence intervals for SI appropriately include the true
values of the parameters except for $x_{2}$, but the confidence intervals for
Naive do not include the true values for any variable.

(a) Case of (E1) and (F3) in the setup of Table 1

(b) Case of (E2) and (F3) in the setup of Table 1

(c) Case of (E1) and (F3) in the setup of Table 2

(d) Case of (E2) and (F3) in the setup of Table 2

Figure 1: Comparison of confidence intervals for SI and Naive. The black (red)
dots are the estimates for SI (Naive), and the solid lines are the confidence
intervals with a confidence coefficient of 0.95.

###  4.2 Case of $p=25$

Table 3 shows the results for when (M2) is the causal effect; i.e., the first
5 variables are active and the remaining 20 variables are inactive. In all
settings, the tuning parameter of Lasso was set to $\lambda=2\sigma
n^{-1/2}(\log p)^{1/2}$. The sequence $\delta_{n}$ used to define the
neighborhood ${\cal N}_{i}^{(h)}$ was set to $(p/2)^{1/2}$ for $X_{ij}\sim{\rm
Ber}(1/2)$ and $(p/6)^{1/2}$ for $X_{ij}\sim{\rm U}(0,1)$. As in the case of
$p=5$, we can see that SI appropriately controls the FCR, keeping it at almost
0.05 in all settings. In addition, since TP is close to 5, the confidence
intervals for the active variables are appropriately away from 0, and FP is
sufficiently small. On the other hand, the FCR for Naive is well above 0.05,
suggesting that reasonable confidence intervals were not obtained. The values
of TP, FP, and $|\hat{M}|$ indicate that almost all of the confidence
intervals for the variables selected by Lasso are away from 0.

Table 3: Cases in which $p=25$ and (M2) is correct. The figures in parentheses are standard deviations. In all settings, SI controls the FCR appropriately, but the inference of Naive is invalid.  |  | $X_{ij}\sim{\rm Ber}(1/2)$ | $X_{ij}\sim{\rm U}(0,1)$ |
---|---|---|---|---
$e(\bm{x})$ | $f(\bm{x})$ |  | $|\hat{M}|$ | TP | FP | FCR | $|\hat{M}|$ | TP | FP | FCR
|  | SI |  | 4.991 | 0.323 | 0.049 |  | 5.000 | 0.002 | 0.051
| (F1) | 11.775 | (0.114) | (0.592) | (0.068) | 5.143 | (0.000) | (0.045) | (0.098)
| Naive | (2.154) | 5.000 | 6.775 | 0.830 | (0.380) | 5.000 | 0.143 | 0.436
|  |  | (0.000) | (2.154) | (0.106) |  | (0.000) | (0.380) | (0.234)
|  | SI |  | 4.972 | 0.879 | 0.058 |  | 4.827 | 0.392 | 0.058
(E1) | (F2) | 19.804 | (0.213) | (1.213) | (0.074) | 11.399 | (0.423) | (0.647) | (0.079)
Naive | (1.979) | 4.982 | 14.822 | 0.965 | (2.023) | 4.852 | 6.547 | 0.893
|  |  | (0.133) | (1.979) | (0.040) |  | (0.385) | (1.997) | (0.093)
|  | SI |  | 4.990 | 0.573 | 0.053 |  | 4.992 | 0.076 | 0.053
| (F3) | 15.960 | (0.148) | (0.847) | (0.065) | 6.676 | (0.126) | (0.273) | (0.089)
| Naive | (2.198) | 5.000 | 10.960 | 0.921 | (1.230) | 4.998 | 1.678 | 0.690
|  |  | (0.000) | (2.198) | (0.066) |  | (0.045) | (1.230) | (0.187)
|  | SI |  | 4.979 | 0.405 | 0.048 |  | 4.994 | 0.013 | 0.053
| (F1) | 13.841 | (0.211) | (0.729) | (0.080) | 5.287 | (0.089) | (0.113) | (0.098)
| Naive | (2.322) | 5.000 | 8.841 | 0.874 | (0.539) | 5.000 | 0.287 | 0.491
|  |  | (0.000) | (2.322) | (0.089) |  | (0.000) | (0.539) | (0.229)
|  | SI |  | 4.938 | 0.834 | 0.053 |  | 4.780 | 0.391 | 0.054
(E2) | (F2) | 20.380 | (0.280) | (0.935) | (0.056) | 11.844 | (0.498) | (0.673) | (0.077)
Naive | (1.906) | 4.954 | 15.425 | 0.968 | (2.121) | 4.818 | 7.026 | 0.904
|  |  | (0.214) | (1.887) | (0.038) |  | (0.435) | (2.085) | (0.084)
|  | SI |  | 4.963 | 0.647 | 0.052 |  | 4.974 | 0.101 | 0.054
| (F3) | 17.238 | (0.309) | (1.028) | (0.075) | 7.050 | (0.222) | (0.359) | (0.099)
| Naive | (2.244) | 5.000 | 12.238 | 0.933 | (1.361) | 4.997 | 2.053 | 0.724
|  |  | (0.000) | (2.244) | (0.059) |  | (0.055) | (1.359) | (0.174)

Figure 2 compares the confidence intervals for SI and Naive with a confidence
coefficient of 0.95 where the nuisance parameters and propensity scores are
(F3) and (E2). As explained in Table 3, Naive has extremely short confidence
intervals due to loss of the validity of the inference, and the confidence
interval does not include 0 even when the true value is 0. On the other hand,
the confidence intervals for SI contain a high percentage of true values.

(a) Case of $X_{ij}\sim{\rm U}(0,1)$ (b) Case of $X_{ij}\sim{\rm Ber}(1/2)$

Figure 2: Comparison of confidence intervals for SI and Naive for cases using
(F3), (E2), and (M2) as models. The black (red) dots are the estimates for SI
(Naive), and the solid lines are the confidence intervals with a confidence
coefficient of 0.95.

##  5 Data analysis

In practice, the error variance $\sigma^{2}$ and propensity score $e(\bm{x})$
must be estimated in some way. As described in Sections 6.1 and 6.3,
appropriate estimators can be substituted for them. Thereafter, the error
variance estimate $\hat{\sigma}^{2}$ can be calculated using the method in
Section 6.1 as $\delta_{n}=(p/2)^{1/2}$, and the propensity score estimate
$\hat{e}(\bm{x})$ can be calculated from a logistic regression model. Using
the estimates, the tuning parameter of Lasso was set to
$\lambda=\hat{\sigma}n^{-1/2}(\log p)^{1/2}$. The following subsections
describe the application of the proposed method to two real datasets.

###  5.1 Application to the lalonde dataset

The lalonde dataset is treated in LaLonde (1986) and is included in the R
package Matching. We set the group that took the U.S. job training program in
1976 as $t=1$ and the group that did not take the program as $t=0$. We
predicted the difference in annual income in 1978 with training as a causal
effect. The confounding variables were age (age), years of education (educ),
black or not (black), Hispanic or not (hisp), married or not (married), high
school graduate or higher (nodegr), income in 1974 (re74), income in 1975
(re75), zero income in 1974 or not (u74), zero income in 1975 or not (u75),
i.e., $p=10$, and the outcome variable was income in 1978 (re78). The sample
size was $n=445$.

When Lasso was used, 5 variables were selected: age, educ, re74, re75 and u74.
For these coefficients, we compared the confidence intervals with a confidence
coefficient of 0.95 given by the proposed method (SI) and a method ignoring
the influence of model selection (Naive). As shown in Table 4, for all
variables, the confidence interval for Naive is shorter than that for SI. This
is because, as mentioned in the previous section, Naive ignores that the
selected variable is likely to be significant. In both methods, the confidence
intervals for educ and u74 do not contain 0, indicating that they are
significant, while the confidence intervals for age and re75 contain 0,
indicating that they are not significant. On the other hand, SI and Naive give
different results for re74. Since u74 is a variable that is 0 if there is
income in 1974 and 1 if there is not, if u74 is 1, the corresponding re74 is
0. Therefore, the SI result can be interpreted as that the causal effect is
affected by whether there was income or not in 1974, rather than by the size
of the income in 1974 itself.

Table 4:  Confidence intervals with a confidence coefficient of 0.95 for the 5 variables obtained by applying Lasso to the lalonde dataset.  | SI | Naive
---|---|---
Selected variable | Estimates | lower | upper | Estimates | lower | upper
age | $49.353$ | $-164.745$ | $165.970$ | $51.345$ | $-66.302$ | $168.991$
educ | $721.188$ | $228.418$ | $1196.124$ | $862.309$ | $390.500$ | $1334.117$
re74 | $0.240$ | $-0.173$ | $0.479$ | $0.236$ | $0.019$ | $0.454$
re75 | $0.275$ | $-0.341$ | $0.645$ | $0.293$ | $-0.040$ | $0.625$
u74 | $6459.088$ | $3709.181$ | $8956.314$ | $6480.235$ | $3997.591$ | $8962.880$

###  5.2 Application to E-commerce dataset

The MineThatData dataset is published by Hillstrom (2008) and contains 12
variables such as purchase amount and district classification code for 64,000
customers. The objective is to predict the purchase price of a product from
the type of email delivered (male, female, or unsent). The sample size of
64,000 was too large computational load for our resources, so we set the
assignment variable to that of emails delivered to males ($t=1$) and unsent
emails ($t=0$). Moreover, we set the sample size to $n=6,382$ and the district
code (zip_code) to Rural. The confounding variables were the number of months
since the last purchase (recency), the amount of purchases made in the last
year (history), whether the user purchased men’s products in the last year
(mens), whether the user purchased women’s products in the last year (womens),
whether the user became a new user in the last 12 months (newbie), whether the
user received the email within two weeks of the delivery (newbie), whether the
user visited the site within two weeks after the mail was delivered (visit),
and whether the user made a purchase within two weeks after the mail was
delivered (conversion), i.e., $p=7$. The outcome variable was the purchase
amount when the user made a purchase (spend).

After the preprocessing described above, the error variance and propensity
scores were estimated, and Lasso selected four variables: recency, history,
newbie and conversion. Table 5 compares the confidence intervals with a
confidence coefficient of 0.95 for these coefficients. In the case of Naive,
all four variables are significant, while in the case of SI, only history and
conversion are significant. A natural interpretation of the SI result is that
it is important for a customer to make a large number of purchases within two
weeks of receiving an email in order to predict the purchase amount. In the
case of Naive, history is negative and significant, indicating that a large
purchase amount over the preceding year has a negative effect on the
prediction of the purchase amount, which is rather unnatural. On the other
hand, in SI, although its estimate is negative, its confidence interval is
positive; i.e., it is a natural interval estimation.

Table 5:  Confidence intervals with a confidence coefficient of 0.95 for the 4 variables obtained by applying Lasso to the E-commerce dataset.  | SI | Naive
---|---|---
Selected variable | Estimates | lower | upper | Estimates | lower | upper
recency | $-0.094$ | $-0.059$ | $0.225$ | $-0.143$ | $-0.274$ | $-0.011$
history | $-0.002$ | $0.000$ | $0.004$ | $-0.003$ | $-0.005$ | $-0.001$
newbie | $0.704$ | $-0.640$ | $1.631$ | $1.419$ | $0.494$ | $2.343$
conversion | $13.378$ | $8.661$ | $17.857$ | $20.191$ | $15.730$ | $24.651$

##  6 Extension

###  6.1 Case of unknown variance

Up to this point, we have treated the variance $\sigma^{2}$ of
$\epsilon^{(h)}$ in (1) as known. However, it is usually unknown in practical
applications. As such, it would seem that we only have to substitute a
consistent estimator, but since the model contains $f(\cdot)$ and we are
trying to develop a method that avoids its identification, its determination
is more difficult than usual. Here, let us choose an appropriate sequence of
real numbers $\\{\delta_{n}\\}$ that converges to $0$, define
${\mathcal{N}}_{i}^{\dagger}=\\{l\neq i:\
\|\bm{X}_{l}-\bm{X}_{i}\|_{2}<\delta_{n},\ \bm{T}_{l}=\bm{T}_{i}\\}$ using
$\bm{T}_{i}=(T_{i}^{(1)},\ldots,T_{i}^{(H)})$, and use

$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\frac{|{\mathcal{N}}_{i}^{\dagger}|}{1+|{\mathcal{N}}_{i}^{\dagger}|}\bigg{(}Y_{i}-\frac{1}{|{\mathcal{N}}_{i}^{\dagger}|}\sum_{l\in{\mathcal{N}}_{i}^{\dagger}}Y_{l}\bigg{)}^{2}.$

If the set of possible values of $\bm{X}_{i}$ is finite, or if the set is
compact and $f(\cdot)$ is continuous, this estimator is obviously a consistent
estimator of $\sigma^{2}$.

###  6.2 Generalization of error structure, model selection method, and causal
effect

In Section 3, we developed selective inference for basic causal inference
models. There is no difficulty in extending it to the case of certain error
structures, model selection methods, and causal effects. Here, we assume that,
in (1), the variance of $\bm{\epsilon}^{(h)}\ (\in\mathbb{R}^{r})$ is
$\bm{\Sigma}$ in the setup of multivariate regression, conduct a model
selection by using an elastic net (Zou and Hastie 2005) as is done in Lee et
al. (2016), and consider the causal effect ${\rm
E}(\sum_{h=1}^{H}c^{(h)}Y^{(h)}\mid T^{(1)}=1)$ in the group receiving the
first treatment as an estimand. Note that this causal effect is the average
treatment effect on the treated or the average treatment effect on the
untreated when $H=2$ (Imbens 2004). Below, we assume the strongly ignorable
treatment assignment condition that assures a valid estimation. The variance
of $\tilde{\bm{\epsilon}}^{(h)}\ (\in\mathbb{R}^{nr})$ in (2) is
$\tilde{\bm{\Sigma}}\otimes\bm{I}_{n}$, which we will denote as
$\tilde{\bm{\Sigma}}$. Let $\lambda_{1}$ and $\lambda_{2}$ be the tuning
parameters. We will use

$\displaystyle\hat{\bm{\beta}}=\mathop{\rm
argmin}_{\bm{\beta}}\bigg{[}\frac{1}{2}\\{\tilde{\bm{W}}(\tilde{\bm{T}},\tilde{\bm{X}})\tilde{\bm{Y}}-\tilde{\bm{X}}\bm{\beta}\\}^{\prime}\tilde{\bm{\Sigma}}^{-1}\\{\tilde{\bm{W}}(\tilde{\bm{T}},\tilde{\bm{X}})\tilde{\bm{Y}}-\tilde{\bm{X}}\bm{\beta}\\}+\lambda_{1}\|\bm{\beta}\|_{1}+\frac{\lambda_{2}}{2}\|\bm{\beta}\|_{2}^{2}\bigg{]}$

instead of (7). Here, $\tilde{\bm{W}}(\tilde{\bm{T}},\tilde{\bm{X}})$ is
$\sum_{h=1}^{H}c^{(h)}\tilde{\bm{W}}^{(h)}(\tilde{\bm{T}}^{(h)},\tilde{\bm{X}})$,
and $\tilde{\bm{W}}^{(h)}(\tilde{\bm{T}}^{(h)},\tilde{\bm{X}})\equiv{\rm
diag}[T_{1}^{(h)}\allowbreak e^{(1)}(\bm{X}_{1})/\\{e^{(h)}(\bm{X}_{1}){\rm
E}(T_{1}^{(1)})\\},\ldots,T_{n}^{(h)}e^{(1)}(\bm{X}_{n})/\\{e^{(h)}(\bm{X}_{n}){\rm
E}(T_{n}^{(1)})\\}]$ is the weight matrix for estimating the causal effect
described above.

In this setup, the coefficients in (9) are given by

$\displaystyle\bm{A}(\bm{M},\bm{s})=\frac{1}{\lambda_{1}}\left(\begin{array}[]{c}\tilde{\bm{X}}_{\bm{M}^{\rm
c}}^{\prime}\tilde{\bm{\Sigma}}^{-1}\\{\bm{I}_{n}-\tilde{\bm{X}}_{\bm{M}}(\tilde{\bm{X}}_{\bm{M}}^{\prime}\tilde{\bm{\Sigma}}^{-1}\tilde{\bm{X}}_{\bm{M}}+\lambda_{2}\bm{I}_{\bm{M}})^{-1}\tilde{\bm{X}}_{\bm{M}}^{\prime}\tilde{\bm{\Sigma}}^{-1}\\}\\\
-\tilde{\bm{X}}_{\bm{M}^{\rm
c}}^{\prime}\tilde{\bm{\Sigma}}^{-1}\\{\bm{I}_{n}-\tilde{\bm{X}}_{\bm{M}}(\tilde{\bm{X}}_{\bm{M}}^{\prime}\tilde{\bm{\Sigma}}^{-1}\tilde{\bm{X}}_{\bm{M}}+\lambda_{2}\bm{I}_{\bm{M}})^{-1}\tilde{\bm{X}}_{\bm{M}}^{\prime}\tilde{\bm{\Sigma}}^{-1}\\}\\\
-{\rm
diag}(\bm{s})(\tilde{\bm{X}}_{\bm{M}}^{\prime}\tilde{\bm{\Sigma}}^{-1}\tilde{\bm{X}}_{\bm{M}}+\lambda_{2}\bm{I}_{\bm{M}})^{-1}\tilde{\bm{X}}_{\bm{M}}^{\prime}\tilde{\bm{\Sigma}}^{-1}\end{array}\right)$

and

$\displaystyle\bm{b}(\bm{M},\bm{s})=\left(\begin{array}[]{c}\bm{1}_{|\bm{M}^{\rm
c}|}-\tilde{\bm{X}}_{\bm{M}^{\rm
c}}^{\prime}\tilde{\bm{\Sigma}}^{-1}\tilde{\bm{X}}_{\bm{M}}(\tilde{\bm{X}}_{\bm{M}}^{\prime}\tilde{\bm{\Sigma}}^{-1}\tilde{\bm{X}}_{\bm{M}}+\lambda_{2}\bm{I}_{\bm{M}})^{-1}\bm{s}\\\
\bm{1}_{|\bm{M}^{\rm c}|}+\tilde{\bm{X}}_{\bm{M}^{\rm
c}}^{\prime}\tilde{\bm{\Sigma}}^{-1}\tilde{\bm{X}}_{\bm{M}}(\tilde{\bm{X}}_{\bm{M}}^{\prime}\tilde{\bm{\Sigma}}^{-1}\tilde{\bm{X}}_{\bm{M}}+\lambda_{2}\bm{I}_{\bm{M}})^{-1}\bm{s}\\\
-{\rm
diag}(\bm{s})(\tilde{\bm{X}}_{\bm{M}}^{\prime}\tilde{\bm{\Sigma}}^{-1}\tilde{\bm{X}}_{\bm{M}}+\lambda_{2}\bm{I}_{\bm{M}})^{-1}\bm{s}\end{array}\right).$

Then, by letting $\tilde{\bm{E}}(\tilde{\bm{X}})\equiv{\rm
diag}\\{e^{(1)}(\bm{X}_{1}),\ldots,e^{(1)}(\bm{X}_{n})\\}$,
$\tilde{\bm{K}}^{(h)}(\tilde{\bm{X}})\equiv{\rm
diag}[e^{(1)}(\bm{X}_{1})\allowbreak\\{1/e^{(h)}(\bm{X}_{1})-1\\}/|{\mathcal{N}}_{1}^{(h)}|,\ldots,e^{(1)}(\bm{X}_{n})\\{1/e^{(h)}(\bm{X}_{n})-1\\}/|{\mathcal{N}}_{n}^{(h)}|]$,

$\displaystyle\tilde{\bm{\eta}}_{j}=\tilde{\bm{\Sigma}}^{-1}\tilde{\bm{X}}_{\bm{M}}(\tilde{\bm{X}}_{\bm{M}}^{\prime}\tilde{\bm{\Sigma}}^{-1}\tilde{\bm{X}}_{\bm{M}})^{-1}\bm{e}_{j},$
$\displaystyle\kappa_{j}^{\bm{M}}(\tilde{\bm{T}},\tilde{\bm{X}})=\bm{e}_{j}^{\prime}(\tilde{\bm{X}}_{\bm{M}}^{\prime}\tilde{\bm{\Sigma}}^{-1}\tilde{\bm{X}}_{\bm{M}})^{-1}\tilde{\bm{X}}_{\bm{M}}^{\prime}\tilde{\bm{\Sigma}}^{-1}\sum_{h=1}^{H}c^{(h)}\tilde{\bm{W}}^{(h)}(\tilde{\bm{T}}^{(h)},\tilde{\bm{X}})\\{\tilde{\bm{\mu}}^{(h)}(\tilde{\bm{X}})+\tilde{\bm{f}}(\tilde{\bm{X}})\\},$
$\displaystyle\zeta_{j}^{\bm{M}}(\tilde{\bm{T}},\tilde{\bm{X}})=\bm{e}_{j}^{\prime}(\tilde{\bm{X}}_{\bm{M}}^{\prime}\tilde{\bm{\Sigma}}^{-1}\tilde{\bm{X}}_{\bm{M}})^{-1}\tilde{\bm{X}}_{\bm{M}}^{\prime}\tilde{\bm{\Sigma}}^{-1}\tilde{\bm{W}}(\tilde{\bm{T}},\tilde{\bm{X}})\tilde{\bm{\Sigma}}\tilde{\bm{W}}(\tilde{\bm{T}},\tilde{\bm{X}})\tilde{\bm{\Sigma}}^{-1}\tilde{\bm{X}}_{\bm{M}}(\tilde{\bm{X}}_{\bm{M}}^{\prime}\tilde{\bm{\Sigma}}^{-1}\tilde{\bm{X}}_{\bm{M}})^{-1}\bm{e}_{j},$
$\displaystyle\tau_{j}^{\bm{M}}(\tilde{\bm{Y}}^{\dagger},\tilde{\bm{T}},\tilde{\bm{X}})=\bm{e}_{j}^{\prime}(\tilde{\bm{X}}_{\bm{M}}^{\prime}\tilde{\bm{\Sigma}}^{-1}\tilde{\bm{X}}_{\bm{M}})^{-1}\tilde{\bm{X}}_{\bm{M}}^{\prime}\tilde{\bm{\Sigma}}^{-1}\sum_{h=1}^{H}c^{(h)}\\{\tilde{\bm{W}}^{(h)}(\tilde{\bm{T}}^{(h)},\tilde{\bm{X}})-\tilde{\bm{E}}(\tilde{\bm{X}})\otimes\bm{I}_{r}\\}\tilde{\bm{Y}}^{(h)\dagger}$

and

$\displaystyle\rho_{j}^{\bm{M}}(\tilde{\bm{T}},\tilde{\bm{X}})=\zeta_{j}^{\bm{M}}(\tilde{\bm{T}},\tilde{\bm{X}})+\tilde{\bm{\eta}}_{j}^{\prime}\sum_{h=1}^{H}c^{(h)2}\\{\tilde{\bm{K}}^{(h)}(\tilde{\bm{X}})\otimes\bm{\Sigma}\\}\tilde{\bm{\eta}}_{j},$

it can be shown that Theorem 1, Lemma 1 and Theorem 2 hold.

###  6.3 Case of unknown propensity score

In Section 3, we treated the propensity scores as known, but in reality, they
are often unknown and are usually estimated by performing some kind of
discriminant analysis. We can formally use the results of Section 3 as if the
estimated propensity scores were known, but for the sake of a more rigorous
methodology, we will take the effect of the estimation into account. We will
regress $\tilde{\bm{x}}$ on $\tilde{\bm{t}}$ with some function and write
$\hat{e}^{(h)}(\bm{X}_{i})$ as the estimator of the propensity score
$e^{(h)}(\bm{X}_{i})={\rm P}(T_{i}^{(h)}=1\mid\bm{X}_{i})$. Here, we assume
that there is no misspecification in the regression model of the propensity
scores and that $\hat{e}^{(h)}(\bm{X}_{i})-e^{(h)}(\bm{X}_{i})={\rm O}_{\rm
P}(n^{-1/2})$. Note that the results in this section do not depend on the
representation of this regression model.

In $\tilde{\bm{\bm{W}}}(\tilde{\bm{\bm{T}}},\tilde{\bm{\bm{X}}})$ which
appears in Section 3, we replace $e^{(h)}(\bm{X}_{i})$ with
$\hat{e}^{(h)}(\bm{X}_{i})$ and write it as
$\hat{\tilde{\bm{W}}}(\tilde{\bm{T}},\tilde{\bm{X}})$. Note that under the
conditioning of (8), $\hat{e}^{(h)}(\bm{X}_{i})$ is a random variable, but
under the conditioning of (12), which is unique to this paper,
$\hat{e}^{(h)}(\bm{X}_{i})$ is a non-random variable. This does not require
using the asymptotic variance of the estimator when the propensity score is
unknown, which has an interesting paradox of being smaller than when it is
known (Hirano et al. 2003). Then, Theorem 1 and Lemma 1 still hold. The
different consideration is needed in the derivation of Theorem 2, where we
evaluate the variance of
$\tilde{\bm{\eta}}_{j}^{\prime}\hat{\tilde{\bm{W}}}(\tilde{\bm{T}},\tilde{\bm{X}})\tilde{\bm{Y}}-\tau_{j}^{\bm{M}}(\tilde{\bm{Y}}^{*},\tilde{\bm{T}},\tilde{\bm{X}})$.
Because of the fluctuations in $\hat{e}^{(h)}(\bm{X}_{i})$, the evaluation of
the matrix in (17) is modified to include ${\rm
E}(\tilde{\bm{\epsilon}}^{(h)*}\tilde{\bm{\epsilon}}^{(k)*}{}^{\prime})\times{\rm
O}_{\rm P}(n^{-1})$. As a result, the conditional expectation of the second
terms still has the form (19). As well, the evaluation of the matrix in (18)
is modified to include ${\rm
E}(\tilde{\bm{\epsilon}}^{(h)}\tilde{\bm{\epsilon}}^{(k)*}{}^{\prime})\times{\rm
O}_{\rm P}(n^{-1/2})$. As a result, ${\rm o}_{\rm P}(n^{-3/2})$ changes to
${\rm O}_{\rm P}(n^{-3/2})$ in the evaluation of (21). Despite this change,
Theorem 2 still holds.

###  6.4 Case of doubly robust estimation

Doubly robust estimation is often used as a potentially more efficient
alternative to inverse-probability-weighted estimation (Robins and Rotnitzky
1995). In this estimation, we directly model the regression of the outcome
variable on the confounding variable, i.e., $\mu^{(h)}(\bm{X})+f(\bm{X})$, and
we write its regression function as $g^{(h)}(\bm{X})$. If either the modeling
or the propensity score modeling is correct, this method can give a consistent
estimator of the causal effect, which is why it is called doubly robust
(Scharfstein et al. 1999). Note that unlike in the previous section, in doubly
robust estimation, it is not necessarily true that
$\hat{e}^{(h)}(\bm{X}_{i})-e^{(h)}(\bm{X}_{i})={\rm O}_{\rm P}(n^{-1/2})$.

For simplicity, we will use
$\\{\bm{A}(\bm{M},\bm{s})\hat{\tilde{\bm{W}}}(\tilde{\bm{T}},\tilde{\bm{X}})\tilde{\bm{Y}}\leq\bm{b}(\bm{M},\bm{s})\\}$
for the conditioned event. Note that, although we will not consider it here,
there is no difficulty in using double robust estimating equation for this
model selection when the estimator of the function $g^{(h)}(\bm{X})$, written
as $\hat{g}^{(h)}(\bm{X})$, is a linear function of $\bm{Y}$. The doubly
robust estimator of $\beta_{j}$ is given by

$\displaystyle\bm{e}_{j}^{\prime}(\tilde{\bm{X}}_{\bm{M}}^{\prime}\tilde{\bm{X}}_{\bm{M}})^{-1}\tilde{\bm{X}}_{\bm{M}}^{\prime}\sum_{h=1}^{H}c^{(h)}[\hat{\tilde{\bm{W}}}^{(h)}(\tilde{\bm{T}}^{(h)},\tilde{\bm{X}})\tilde{\bm{Y}}^{(h)}+\\{\bm{I}_{n}-\hat{\tilde{\bm{W}}}^{(h)}(\tilde{\bm{T}}^{(h)},\tilde{\bm{X}})\\}\hat{\tilde{\bm{g}}}^{(h)}(\tilde{\bm{X}})].$
(22)

Here,
$\hat{\tilde{\bm{g}}}^{(h)}(\tilde{\bm{X}})=(\hat{g}^{(h)}(\bm{X}_{1}),\ldots,\hat{g}^{(h)}(\bm{X}_{n}))^{\prime}$,
and we assume that when this modeling is correct, the differences between it
and $\tilde{\bm{\mu}}^{(h)}(\tilde{\bm{X}})+\tilde{\bm{f}}(\tilde{\bm{X}})$
and between it and ${\rm
E}\\{\hat{\tilde{\bm{g}}}^{(h)}(\tilde{\bm{X}})\mid\tilde{\bm{X}}\\}$ are each
${\rm O}_{\rm P}(n^{-1/2})$. Then, by defining

$\displaystyle\tau_{j}^{\bm{M}}(\tilde{\bm{Y}}^{\dagger},\tilde{\bm{T}},\tilde{\bm{X}})\equiv\bm{e}_{j}^{\prime}(\tilde{\bm{X}}_{\bm{M}}^{\prime}\tilde{\bm{X}}_{\bm{M}})^{-1}\tilde{\bm{X}}_{\bm{M}}^{\prime}\sum_{h=1}^{H}c^{(h)}\\{\hat{\tilde{\bm{W}}}^{(h)}(\tilde{\bm{T}}^{(h)},\tilde{\bm{X}})-\bm{I}_{n}\\}[\tilde{\bm{Y}}^{(h)\dagger}-{\rm
E}\\{\hat{\tilde{\bm{g}}}^{(h)}(\tilde{\bm{X}})\mid\tilde{\bm{X}}\\}],$

and replacing the estimator
$\tilde{\bm{\eta}}_{j}^{\prime}\tilde{\bm{W}}(\tilde{\bm{T}},\tilde{\bm{X}})\tilde{\bm{Y}}$
with (22), we can see that Lemma 1 holds.

Next, as $\tilde{\bm{Y}}^{(h)\dagger}$ above, we consider using its natural
estimator like $\tilde{\bm{Y}}^{(h)*}$ in Section 3. However,
$\tilde{\bm{Y}}^{(h)*}$ converges slowly, so if the modeling of the propensity
score is incorrect and the expectation of
$\hat{\tilde{\bm{W}}}^{(h)}(\tilde{\bm{T}}^{(h)},\tilde{\bm{X}})-\bm{I}_{n}$
does not converge to the zero matrix, we cannot expect quick convergence in
evaluations like those of (16) and (21). As a naive solution, we propose to
use a goodness-of-fit statistic for $\hat{e}^{(h)}(\bm{X})$ and substitute
$\tilde{\bm{Y}}^{(h)*}-\hat{\tilde{\bm{g}}}^{(h)}(\tilde{\bm{X}})$ or
$\bm{0}_{n}$ for $\tilde{\bm{Y}}^{(h)\dagger}-{\rm
E}\\{\hat{\tilde{\bm{g}}}^{(h)}(\tilde{\bm{X}})\mid\tilde{\bm{X}}\\}$ by case
separation. For example, suppose that the possible values of
$\hat{e}^{(h)}(\bm{X})$ estimated using a $q$-dimensional discriminant model
are discrete, and the estimate is $\hat{e}^{(h)}_{r}$ when $\bm{X}\in{\cal
X}_{r}\ (\subset\mathbb{R}^{p})$, where $\bigcup_{r=1}^{R}{\cal
X}_{r}=\mathbb{R}^{p}$. Then, by letting
$\hat{p}^{(h)}_{r}\equiv\sum_{i:\bm{X}_{i}\in{\cal
X}_{r}}T_{i}^{(h)}/|\\{i:\bm{X}_{i}\in{\cal X}_{r}\\}|$, we only have to
consider
$d(\tilde{\bm{T}},\tilde{\bm{X}})\equiv-2\sum_{r=1}^{R}\sum_{i:\bm{X}_{i}\in{\cal
X}_{r}}\sum_{h=1}^{H}T_{i}^{(h)}\log(\hat{p}^{(h)}_{r}/\hat{e}^{(h)}_{r})+(LH-q)\log
n$, which is equivalent to a difference between Bayesian information criteria
(BIC), as the goodness-of-fit statistic. If this goodness-of-fit statistic is
positive, the propensity score model is regarded as valid and we use
$\tilde{\bm{Y}}^{(h)*}-\hat{\tilde{\bm{g}}}^{(h)}(\tilde{\bm{X}})$. If it is
negative, the propensity score model is regarded as not valid and we use
$\bm{0}_{n}$. Then, from the consistency of BIC, we can ensure the convergence
as in Section 3. Specifically, we define

$\displaystyle\tau_{j}^{\bm{M}}(\tilde{\bm{Y}}^{*},\tilde{\bm{T}},\tilde{\bm{X}})$
$\displaystyle\equiv\bm{e}_{j}^{\prime}(\tilde{\bm{X}}_{\bm{M}}^{\prime}\tilde{\bm{X}}_{\bm{M}})^{-1}\tilde{\bm{X}}_{\bm{M}}^{\prime}I\\{d(\tilde{\bm{T}},\tilde{\bm{X}})>0\\}\sum_{h=1}^{H}c^{(h)}\\{\hat{\tilde{\bm{W}}}^{(h)}(\tilde{\bm{T}}^{(h)},\tilde{\bm{X}})-\bm{I}_{n}\\}\\{\tilde{\bm{Y}}^{(h)*}-\hat{\tilde{\bm{g}}}^{(h)}(\tilde{\bm{X}})\\}.$

In line with this, we will let $\tilde{\bm{K}}^{(h)}(\tilde{\bm{X}})\equiv{\rm
diag}[\\{1/e^{(h)}(\bm{X}_{1})-1\\}/|{\mathcal{N}}_{1}^{(h)}|,\ldots,\\{1/e^{(h)}(\bm{X}_{n})-1\\}/|{\mathcal{N}}_{n}^{(h)}|]$.
Moreover, we define

$\displaystyle\rho_{j}^{\bm{M}}(\tilde{\bm{T}},\tilde{\bm{X}})\equiv\zeta_{j}^{\bm{M}}(\tilde{\bm{T}},\tilde{\bm{X}})+\sigma^{2}I\\{d(\tilde{\bm{T}},\tilde{\bm{X}})>0\\}\sum_{h=1}^{H}c^{(h)2}\tilde{\bm{\eta}}_{j}^{\prime}\tilde{\bm{K}}^{(h)}(\tilde{\bm{X}})\tilde{\bm{\eta}}_{j},$

where $I\\{\cdot\\}$ is the indicator function. Using these definitions, it
can be seen that Theorem 2 holds.

###  6.5 Case of time-varying confounding

Dealing with time-varying confounding is a fundamental topic in causal
inference (Bang and Robins 2005; Daniel et al. 2013). Here, suppose that there
are $m$ time points, indexed by $j\ (\in\\{1,\ldots,m\\})$. Let $T_{j}$ be the
assignment variable representing the binary time-varying treatment and
$\bm{X}_{j}$ be the time-varying confounding variable at time $j$. The
collections $\bm{T}=(T_{1},\ldots,T_{m})^{\prime}$ and
$\bm{X}=(\bm{X}_{1}^{\prime},\ldots,\bm{X}_{m}^{\prime})^{\prime}$ are called
the treatment history and the confounding history, respectively. Let us denote
the realization of $\bm{T}$ as $\bm{t}=(t_{1},\ldots,t_{m})^{\prime}$ and the
potential outcome variable for $\bm{t}$ as $Y^{(\bm{t})}$. We assume the
following structure,

$\displaystyle
Y=\sum_{\bm{t}\in\\{0,1\\}^{m}}I(\bm{T}=\bm{t})Y^{(\bm{t})}=\sum_{\bm{t}\in\\{0,1\\}^{m}}I(\bm{T}=\bm{t})\\{\mu^{(\bm{t})}+f(\bm{X})+\epsilon\\}$

as in Daniel et al. (2013). Then, for this $\mu^{(\bm{t})}$, we consider the
model $(\bm{t}-\bm{1}_{m}/2)^{\prime}\bm{\beta}$. Note that $\beta_{j}$ is a
parameter that represents the treatment effect at time $j$. As in the
preceding sections, the variables for the $i$-th sample are denoted by the
subscript i and their collection is denoted by the tilde symbol $\tilde{\ }$.
In this setting, the ignorable treatment assignment condition defined by (3)
does not give a consistent estimation, so

$\displaystyle
Y_{i}^{(\bm{t})}\mathop{\perp\\!\\!\\!\\!\perp}T_{j,i}\mid\bm{X}_{1,i},\ldots,\bm{X}_{j,i},T_{1,i},\ldots,T_{j-1,i}\qquad(\bm{t}\in\\{0,1\\}^{m};\
i=1,\ldots,n;\ j=1,\ldots,m)$

is assumed instead. Here, when $j=1$, the assignment variable is assumed to be
unconditioned.

The propensity score in this setting is $e^{(\bm{t})}(\bm{X})={\rm
P}(T_{1}=t_{1}\mid\bm{X}_{1})\prod_{j=2}^{m}{\rm
P}(T_{j}=t_{j}\mid\bm{X}_{1},\ldots,\bm{X}_{j},T_{1}=t_{1},\ldots,T_{j-1}=t_{j-1})$,
and the weight matrix is
$\tilde{\bm{W}}^{(\bm{t})}(\tilde{\bm{T}},\tilde{\bm{X}})={\rm
diag}\\{I(\bm{T}_{1}=\bm{t})/e^{(\bm{t})}(\bm{X}_{1}),\ldots,I(\bm{T}_{n}\allowbreak=\bm{t})/e^{(\bm{t})}(\bm{X}_{n})\\}$.
We conduct model selection by minimizing
$\sum_{\bm{t}\in\\{0,1\\}^{m}}\allowbreak\|\tilde{\bm{W}}^{(\bm{t})}(\tilde{\bm{T}},\tilde{\bm{X}})\tilde{\bm{Y}}-\bm{1}_{n}(\bm{t}-\bm{1}_{m}/2)^{\prime}\bm{\beta}\|_{2}^{2}/2+\lambda\|\bm{\beta}\|_{1}$.
Once a model $\bm{M}$ is selected, we use the inverse-probability-weighted
estimator,

$\displaystyle\frac{1}{n2^{|\bm{M}|-2}}\bm{e}_{j}^{\prime}\sum_{\bm{t}\in\\{0,1\\}^{m}}\bigg{(}\bm{t}^{\bm{M}}-\frac{\bm{1}_{|\bm{M}|}}{2}\bigg{)}\bm{1}_{n}^{\prime}\tilde{\bm{W}}^{(\bm{t})}(\tilde{\bm{T}},\tilde{\bm{X}})\tilde{\bm{Y}},$

as a consistent estimator of $\beta_{j}^{\bm{M}}$. Although the model
selection and estimation have different expressions from those in Section 3,
it makes no difference that the squared losses are based on a linear function
of $\tilde{\bm{Y}}$ minus a linear function of $\bm{\beta}$. Therefore, for
example, an event in which model $\bm{M}$ with sign $\bm{s}$ is selected is
expressed as
$\bm{A}(\bm{M},\bm{s},\tilde{\bm{T}},\tilde{\bm{X}})\tilde{\bm{Y}}\leq\bm{b}(\bm{M},\bm{s})$,
and ${\mathcal{V}}_{\bm{s},j}^{-}$ and ${\mathcal{V}}_{\bm{s},j}^{+}$ are
easily determined, i.e., the argument can be developed in the same way as in
Section 3. Specifically, by letting
$\tilde{\bm{K}}^{(\bm{t})}(\tilde{\bm{X}})\equiv{\rm
diag}[\\{1/e^{(\bm{t})}(\bm{X}_{1})-1\\}/|{\cal
N}_{1}^{(\bm{t})}|,\ldots,\\{1/e^{(\bm{t})}(\bm{X}_{n})-1\\}/|{\cal
N}_{n}^{(\bm{t})}|]$,

$\displaystyle\kappa_{j}^{\bm{M}}(\tilde{\bm{T}},\tilde{\bm{X}})=\frac{1}{n2^{|\bm{M}|-2}}\bm{e}_{j}^{\prime}\sum_{\bm{t}\in\\{0,1\\}^{m}}\bigg{(}\bm{t}^{\bm{M}}-\frac{\bm{1}_{|\bm{M}|}}{2}\bigg{)}\bm{1}_{n}^{\prime}\tilde{\bm{W}}^{(\bm{t})}(\tilde{\bm{T}},\tilde{\bm{X}})\\{\tilde{\bm{\mu}}^{(\bm{t})}+\tilde{\bm{f}}(\bm{X})\\},$
$\displaystyle\zeta_{j}^{\bm{M}}(\tilde{\bm{T}},\tilde{\bm{X}})=\frac{\sigma^{2}}{n^{2}2^{2|\bm{M}|-4}}\sum_{\bm{t}\in\\{0,1\\}^{m}}\bm{e}_{j}^{\prime}\bigg{(}\bm{t}^{\bm{M}}-\frac{\bm{1}_{|\bm{M}|}}{2}\bigg{)}\bm{1}_{n}^{\prime}\tilde{\bm{W}}^{(\bm{t})}(\tilde{\bm{T}},\tilde{\bm{X}})^{2}\bm{1}_{n}\bigg{(}\bm{t}^{\bm{M}}-\frac{\bm{1}_{|\bm{M}|}}{2}\bigg{)}^{\prime}\bm{e}_{j},$
$\displaystyle\tau_{j}^{\bm{M}}(\tilde{\bm{Y}}^{\dagger},\tilde{\bm{T}},\tilde{\bm{X}})=\frac{1}{n2^{|\bm{M}|-2}}\bm{e}_{j}^{\prime}\sum_{\bm{t}\in\\{0,1\\}^{m}}\bigg{(}\bm{t}^{\bm{M}}-\frac{\bm{1}_{|\bm{M}|}}{2}\bigg{)}\bm{1}_{n}^{\prime}\\{\tilde{\bm{W}}^{(\bm{t})}(\tilde{\bm{T}},\tilde{\bm{X}})-\bm{I}_{n}\\}\tilde{\bm{Y}}^{(\bm{t})\dagger}$

and

$\displaystyle\rho_{j}^{\bm{M}}(\tilde{\bm{T}},\tilde{\bm{X}})=\zeta_{j}^{\bm{M}}(\tilde{\bm{T}},\tilde{\bm{X}})+\frac{\sigma^{2}}{n^{2}2^{2|\bm{M}|-4}}\sum_{\bm{t}\in\\{0,1\\}^{m}}\bm{e}_{j}^{\prime}\bigg{(}\bm{t}^{\bm{M}}-\frac{\bm{1}_{|\bm{M}|}}{2}\bigg{)}\bm{1}_{n}^{\prime}\tilde{\bm{K}}^{(\bm{t})}(\tilde{\bm{X}})\bm{1}_{n}\bigg{(}\bm{t}^{\bm{M}}-\frac{\bm{1}_{|\bm{M}|}}{2}\bigg{)}^{\prime}\bm{e}_{j},$

it can be seen that Theorem 1, Lemma 1 and Theorem 2 hold.

##  7 Conclusion

We developed selective inference in a basic causal inference model where the
least-squares method is inappropriate for when inverse-probability-weighted
estimation is used as the standard for propensity score analysis. Conditioning
on the assignment variables first is usually not done in propensity score
analysis, but doing so, we were able to construct a pivot statistic that
follows a uniform distribution. Moreover, by dealing with higher-order
asymptotics, we gave asymptotically guaranteed confidence intervals in the
context of selective inference. Numerical experiments showed that a method
that ignores the quiet scandal of statistics results in significant deviations
from the preset coverage of the confidence intervals, whereas our method
maintains the coverage. By comparing the two methods on benchmark data, we
showed that there is a significant difference between them.

To carry on the beautiful methodology described by Lee et al. (2016), we
considered in this paper model selection with basic sparse estimation in a
setting that follows a normal distribution if there is no assignment. However,
it is also important to consider model selection with an information criterion
in a setting that includes a non-normal distribution. The information
criterion for this case was recently developed by Baba et al. (2017), and it
is expected that there will be no difficulty in implementing it in the method
of Charkhi and Claeskens (2018). On the other hand, in order to treat the
error distribution non-parametrically and to construct an elaborate theory
based on the uniform convergence of pivot statistics, it is necessary to
develop the works of Tian and Taylor (2017) and Tibshirani et al. (2018).

The selective inference is developing rapidly, and there are ongoing
discussions on increasing the power of test. In fact, Kivaranovic and Leeb
(2021) has shown that the standard method of selective inference can give too
wide confidence intervals due to ${\rm
E}(U_{j}^{\hat{\bm{M}}}-L_{j}^{\hat{\bm{M}}})=\infty$. Then, for example, Tian
and Taylor (2018) proposed a modern version of sample splitting, specifically
adding an appropriate randomization to the data, which is further discussed in
Rasines and Young (2021). Depending on how the randomization is realized,
these methods may yield different results, and so they cannot be accepted in
all the fields of application, however, a considerable improvement in the
power has been observed. In addition, Liu et al. (2018) pointed out that it is
over-conditioning to condition on $\hat{\bm{M}}=\bm{M}$ when inferring each
$\beta_{j}$, and increased the power by conditioning on a wider region. Since
the selective inference in propensity score analysis involves many non-trivial
arguments in constructing valid inferences themselves, and since the above-
mentioned improvements are expected to be realized by simply combining them,
we have not dealt with them in this paper. However, they are significant and
will be an important subject in the future.

## Acknowledgement

Yoshiyuki Ninomiya was supported by JSPS KAKENHI (16K00050), Yuta Umezu was
supported by JSPS KAKENHI (18K18010), and Ichiro Takeuchi was supported by
JSPS KAKENHI (17H00758, 20H00601), JST CREST (JPMJCR1502) and RIKEN Center for
Advanced Intelligence Project.

## References

  * (1)
  * Baba et al. (2017) Baba, T., Kanemori, T., and Ninomiya, Y. (2017). A Cp criterion for semiparametric causal inference, Biometrika, 104, 845–861.
  * Bang and Robins (2005) Bang, H. and Robins, J. M. (2005). Doubly robust estimation in missing data and causal inference models, Biometrics, 61, 962–973.
  * Berk et al. (2013) Berk, R., Brown, L., Buja, A., Zhang, K., and Zhao, L. (2013). Valid post-selection inference, The Annals of Statistics, 41, 802–837.
  * Breiman (2001) Breiman, L. (2001). Statistical modeling: The two cultures, Statistical Science, 16, 199–231.
  * Charkhi and Claeskens (2018) Charkhi, A. and Claeskens, G. (2018). Asymptotic post-selection inference for the Akaike information criterion, Biometrika, 105, 645–664.
  * Daniel et al. (2013) Daniel, R. M., Cousens, S. N., De Stavola, B. L., Kenward, M. G., and Sterne, J. A. C. (2013). Methods for dealing with time-dependent confounding, Statistics in Medicine, 32, 1584–1618.
  * Fithian et al. (2014) Fithian, W., Sun, D., and Taylor, J. (2014). Optimal inference after model selection, arXiv preprint:1410.2597.
  * Hernán and Robins (2020) Hernán, M. A. and Robins, J. M. (2020). Causal inference: what if: Boca Raton: Chapman & Hall/CRC.
  * Hillstrom (2008) Hillstrom, K. (2008). The MineThatData e-mail analytics and data mining challenge, MineThatData blog, http://blog.minethatdata.com/2008/03/minethatdata-e-mail-analytics-and-data.html.
  * Hirano et al. (2003) Hirano, K., Imbens, G. W., and Ridder, G. (2003). Efficient estimation of average treatment effects using the estimated propensity score, Econometrica, 71, 1161–1189.
  * Imbens (2000) Imbens, G. W. (2000). The role of the propensity score in estimating dose-response functions, Biometrika, 87, 706–710.
  * Imbens (2004) (2004). Nonparametric estimation of average treatment effects under exogeneity: A review, Review of Economics and Statistics, 86, 4–29.
  * Kivaranovic and Leeb (2021) Kivaranovic, D. and Leeb, H. (2021). On the length of post-model-selection confidence intervals conditional on polyhedral constraints, Journal of the American Statistical Association, 116, 845–857.
  * LaLonde (1986) LaLonde, R. J. (1986). Evaluating the econometric evaluations of training programs with experimental data, The American Economic Review, 604–620.
  * 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, The Annals of Statistics, 44, 907–927.
  * Lee and Taylor (2014) Lee, J. D. and Taylor, J. E. (2014). Exact post model selection inference for marginal screening, In Advances in Neural Information Processing Systems, 27, 136–144.
  * Liu et al. (2018) Liu, K., Markovic, J., and Tibshirani, R. (2018). More powerful post-selection inference, with application to the lasso, arXiv preprint:1801.09037.
  * Peters et al. (2017) Peters, J., Janzing, D., and Schölkopf, B. (2017). Elements of causal inference: foundations and learning algorithms: The MIT Press.
  * Rasines and Young (2021) Rasines, D. G. and Young, G. A. (2021). Splitting strategies for post-selection inference, arXiv preprint arXiv:2102.02159.
  * Robins et al. (1994) Robins, J. M., Rotnitzky, A., and Zhao, L. P. (1994). Estimation of regression coefficients when some regressors are not always observed, Journal of the American Statistical Association, 89, 846–866.
  * Robins and Rotnitzky (1995) Robins, J. M. and Rotnitzky, A. (1995). Semiparametric efficiency in multivariate regression models with missing data, Journal of the American Statistical Association, 90, 122–129.
  * Rosenbaum and Rubin (1983) Rosenbaum, P. R. and Rubin, D. B. (1983). The central role of the propensity score in observational studies for causal effects, Biometrika, 70, 41–55.
  * Rubin (1985) Rubin, D. B. (1985). The use of propensity scores in applied Bayesian inference, Bayesian Statistics, 2, 463–472.
  * Scharfstein et al. (1999) Scharfstein, D. O., Rotnitzky, A., and Robins, J. M. (1999). Adjusting for nonignorable drop-out using semiparametric nonresponse models, Journal of the American Statistical Association, 94, 1096–1120.
  * Taylor and Tibshirani (2015) Taylor, J. and Tibshirani, R. J. (2015). Statistical learning and selective inference, Proceedings of the National Academy of Sciences, 112, 7629–7634.
  * Tian and Taylor (2017) Tian, X. and Taylor, J. (2017). Asymptotics of selective inference, Scandinavian Journal of Statistics, 44, 480–499.
  * Tian and Taylor (2018) (2018). Selective inference with a randomized response, The Annals of Statistics, 46, 679–710.
  * Tibshirani (1996) Tibshirani, R. (1996). Regression shrinkage and selection via the lasso, Journal of the Royal Statistical Society: Series B, 58, 267–288.
  * Tibshirani et al. (2016) Tibshirani, R. J., Taylor, J., Lockhart, R., and Tibshirani, R. (2016). Exact post-selection inference for sequential regression procedures, Journal of the American Statistical Association, 111, 600–620.
  * Tibshirani et al. (2018) Tibshirani, R. J., Rinaldo, A., Tibshirani, R., and Wasserman, L. (2018). Uniform asymptotic inference and the bootstrap after model selection, Annals of Statistics, 46, 1255–1287.
  * Zhao et al. (2017) Zhao, Q., Small, D. S., and Ertefaie, A. (2017). Selective inference for effect modification via the lasso, arXiv preprint:1705.08020.
  * Zou and Hastie (2005) Zou, H. and Hastie, T. (2005). Regularization and variable selection via the elastic net, Journal of the Royal Statistical Society: Series B, 67, 301–320.