# Information criteria for sparse methods in causal inference

Yoshiyuki Ninomiya
Department of Statistical Inference and Mathematics
The Institute of Statistical Mathematics 10-3 Midori-cho, Tachikawa-shi, Tokyo
190-8562, Japan. E-mail<EMAIL_ADDRESS>
###### Abstract

For propensity score analysis and sparse estimation, we develop an information
criterion for determining the regularization parameters needed in variable
selection. First, for Gaussian distribution-based causal inference models, we
extend Stein’s unbiased risk estimation theory, which leads to a generalized
Cp criterion that has almost no weakness in conventional sparse estimation,
and derive an inverse-probability-weighted sparse estimation version of the
criterion without resorting to asymptotics. Next, for general causal inference
models that are not necessarily Gaussian distribution-based, we extend the
asymptotic theory on LASSO for propensity score analysis, with the intention
of implementing doubly robust sparse estimation. From the asymptotic theory,
an AIC-type information criterion for inverse-probability-weighted sparse
estimation is given, and then a criterion with double robustness in itself is
derived for doubly robust sparse estimation. Numerical experiments compare the
proposed criterion with the existing criterion derived from a formal argument
and verify that the proposed criterion is superior in almost all cases, that
the difference is not negligible in many cases, and that the results of
variable selection differ significantly. Real data analysis confirms that the
difference between variable selection and estimation by these criteria is
actually large. Finally, generalizations to general sparse estimation using
group LASSO, elastic net, and non-convex regularization are made in order to
indicate that the proposed criterion is highly extensible.

Keywords: AIC; Cp; doubly robust; inverse-probability-weighted; LASSO; model
selection; propensity score analysis; statistical asymptotic theory; SURE
theory

##  1 Introduction

Sparse estimation is indispensable as a regression analysis method when there
are many candidate explanatory variables. Not only LASSO, proposed by
Tibshirani (1996), but also sparse estimation methods such as SCAD (Fan and Li
2001), elastic net (Zou and Hastie 2005) and MC+ (Zhang 2010), which were
developed later to improve LASSO, are already standard. The best subset
regression, in which all subsets of variables are candidates for models and
the optimal one is selected using an information criterion or a cross-
validation method, tends not to return appropriate results because the number
of candidates varies considerably for each size of variable set. This weakness
does not exist in sparse estimation, which performs variable selection
simultaneously with the estimation. On the other hand, an appropriate value of
the regularization parameter must be determined for appropriate variable
selection. Although the cross-validation method is often used for this
purpose, the generalized Cp criterion developed in Efron et al. (2004) and Zou
et al. (2007) has almost no weaknesses when it is used in LASSO. Since sparse
estimation uses an estimating function with non-differentiable points, the
conventional Cp criterion cannot be derived; however, by using Stein (1981)’s
unbiased risk estimation (SURE) theory, an unbiased estimator of the mean
squared error can be simply derived without resorting to asymptotics. Although
the generalized Cp criterion is obtained in Gaussian linear regression
settings, in generalized linear regression settings, one has only to rely on
asymptotics; indeed, Ninomiya and Kawano (2016) derived the AIC and
demonstrated its good performance numerically.

On the other hand, causal inference has long been a topic in biostatistics and
econometrics and has recently attracted attention in the machine-learning
community (see, for example, Peters et al. 2017, Hernán and Robins 2020).
Here, the semiparametric approach to propensity score analysis is a standard
theoretically guaranteed method that avoids the bias caused by the presence of
confounding. To illustrate the motivation for this paper, let us consider the
simplest model in causal inference, i.e.,

$\displaystyle
y_{i}=(1-t_{i})y_{i}^{[1]}+t_{i}y_{i}^{[2]}=(1-t_{i})\bm{x}_{i}^{{\rm
T}}\bm{\theta}^{[1]}+t_{i}\bm{x}_{i}^{{\rm T}}\bm{\theta}^{[2]}+\xi_{i},\qquad
i\in\\{1,2,\ldots,N\\},$ (1)

where the causal effect can be written as a linear sum of covariates. The
subscript $i$ denotes the $i$-th sample, $y_{i}\ (\in\mathbb{R})$ is an
outcome variable, $t_{i}\ (\in\\{0,1\\})$ is an assignment variable,
$\bm{x}_{i}\ (\in\mathbb{R}^{p})$ is an explanatory variable, and $\xi_{i}\
(\in\mathbb{R})$ is an error variable with expectation 0. In causal inference,
it is usually supposed that $y_{i}$ and $t_{i}$ are correlated, and an
estimation without considering the correlation will lead to severe bias. Let
us assume that a confounding variable resolving this problem is observed as
$\bm{z}_{i}\ (\in\mathbb{R}^{q})$. The vector $\bm{x}_{i}$ may contain some of
the $\bm{z}_{i}$ as elements, while the other elements are assumed to be
independent of $\bm{z}_{i}$ for simplicity. In addition, let us assume that
$\xi_{i}$ depends on $\bm{z}_{i}$ and its conditional distribution given
$\bm{z}_{i}$ is Gaussian. In this model, the causal effect of interest is
supposed to be $\bm{\theta}^{[2]}-\bm{\theta}^{[1]}$. Then, using the
propensity score $e(\bm{z}_{i})\equiv{\rm P}(t_{i}=1\mid\bm{z}_{i})$
introduced by Rosenbaum and Rubin (1983) and weighting $y_{i}$ by
$t_{i}/e(\bm{z}_{i})-(1-t_{i})/\\{1-e(\bm{z}_{i})\\}$, we can estimate
$\bm{\theta}^{[2]}-\bm{\theta}^{[1]}$ without bias by using the least squares
method. This is one of the inverse-probability-weighted estimation methods
(Rubin 1985, Robins et al. 1994). It also has the advantage of not requiring
an estimation of ${\rm E}(\xi_{i}\mid\bm{z}_{i})$, which is difficult to
model, and is what is relied upon in this study.

In this paper, sparse estimation is incorporated into propensity score
analysis, supposing that the dimension $p$ of the explanatory variables is
somewhat high. Specifically, for example, $\bm{\theta}$ which minimizes

$\displaystyle\sum_{i=1}^{N}\Bigg{[}\Bigg{\\{}\frac{t_{i}}{e(\bm{z}_{i})}-\frac{1-t_{i}}{1-e(\bm{z}_{i})}\Bigg{\\}}y_{i}-\bm{x}_{i}^{{\rm
T}}\bm{\theta}\Bigg{]}^{2}+\lambda\|\bm{\theta}\|_{1}$

is given as a LASSO-type estimator, where $\|\cdot\|_{1}$ denotes the
$\ell_{1}$ norm. This method is not new and has been part of proposals in,
e.g., Ning et al. (2020), Dukes and Vansteelandt (2021) and Zhao et al.
(2022). Then, for the essential problem of determining $\lambda$, we develop
an inverse-probability Cp criterion (IPCp) based on SURE theory as a method
with almost no weaknesses. Although the outcome variable $y_{i}$ does not
follow a Gaussian distribution, our development is made possible by
conditioning on the assignment variable $t_{i}$ earlier, a concept not usually
used in propensity score analysis.

When using doubly robust estimation (Scharfstein et al. 1999, Bang and Robins
2005), which has also been the focus of attention in propensity score
analysis, or when using general causal inference models that are not
necessarily Gaussian distribution-based, the information criterion must be
derived by relying on asymptotics. For simplicity of explanation, a Gaussian
distribution-based model in (1) is considered, while we treat a general causal
effect and provide it by estimating $\bm{\theta}^{[1]}$ and
$\bm{\theta}^{[2]}$ separately. In doubly robust estimation, by first assuming
some parametric function $e(\bm{z}_{i};\bm{\alpha})$ for the unknown
$e(\bm{z}_{i})$ and next estimating $\bm{\alpha}$ as $\hat{\bm{\alpha}}$
appropriately, we substitute $e(\bm{z}_{i};\hat{\bm{\alpha}})$ for
$e(\bm{z}_{i})$. That is, when incorporating LASSO, we have only to find
$(\bm{\theta}^{[1]},\bm{\theta}^{[2]})$ that minimizes

$\displaystyle\sum_{i=1}^{N}\Bigg{\\{}\frac{1-t_{i}}{1-e(\bm{z}_{i};\hat{\bm{\alpha}})}(y_{i}-\bm{x}_{i}^{{\rm
T}}\bm{\theta}^{[1]})^{2}+\frac{t_{i}}{e(\bm{z}_{i};\hat{\bm{\alpha}})}(y_{i}-\bm{x}_{i}^{{\rm
T}}\bm{\theta}^{[2]})^{2}\Bigg{\\}}+\lambda(\|\bm{\theta}^{[1]}\|_{1}+\|\bm{\theta}^{[2]}\|_{1}).$

Furthermore, considering that this model may be misspecified, we also assume
some parametric function for the conditional distribution of the outcome
variable given the confounding variable. Then, adding

$\displaystyle\sum_{i=1}^{N}\Bigg{[}\Bigg{\\{}1-\frac{1-t_{i}}{1-e(\bm{z}_{i};\hat{\bm{\alpha}})}\Bigg{\\}}{\rm
E}\\{(y_{i}^{[1]}-\bm{x}_{i}^{{\rm
T}}\bm{\theta}^{[1]})^{2}\mid\bm{z}_{i};\hat{\bm{\gamma}}\\}+\Bigg{\\{}1-\frac{t_{i}}{e(\bm{z}_{i};\hat{\bm{\alpha}})}\Bigg{\\}}{\rm
E}\\{(y_{i}^{[2]}-\bm{x}_{i}^{{\rm
T}}\bm{\theta}^{[2]})^{2}\mid\bm{z}_{i};\hat{\bm{\gamma}}\\}\Bigg{]}$

to the above term, we minimize it with respect to
$(\bm{\theta}^{[1]},\bm{\theta}^{[2]})$. Here, $\hat{\bm{\gamma}}$ is an
appropriate estimator of the parameters appearing in the parametric function,
and the expectation is taken with the conditional distribution. In so doing,
as long as the modeling of either $e(\bm{z}_{i};\bm{\alpha})$ or ${\rm
E}(\cdot\mid\bm{z}_{i};\bm{\gamma})$ is correct, the estimation of
$(\bm{\theta}^{[1]},\bm{\theta}^{[2]})$ becomes valid.

In order to obtain the asymptotic properties of the estimator in this setting,
it is necessary to generalize Knight and Fu (2000) and Ninomiya and Kawano
(2016). Then, as in Baba et al. (2017), we define a divergence suitable for
causal inference, i.e., the risk based on the loss function used in the
estimation, and derive an asymptotically unbiased estimator of it by using the
asymptotic properties. This means that the same procedure as for the
conventional Akaike information criterion (AIC; Akaike 1973) is used to derive
the information criterion. However, without special care, double robustness is
merely a property of the estimator. Therefore, as in Baba and Ninomiya (2021),
we derive an asymptotically unbiased estimator so that the information
criterion itself is doubly robust; we call it the doubly robust information
criterion (DRIC). DRIC becomes an asymptotically unbiased estimator of
divergence for causal inference as long as either the model of the assignment
variable or the model of the outcome variable is correctly specified.

Deriving an information criterion, a fundamental tool in statistical analysis,
has been a long-sought goal of propensity score analysis and sparse
estimation. Although fundamental, no good information criterion has been given
to date. The reason may be that certain information criteria that are
extensions of AIC, such as the regularization information criterion (RIC;
Shibata 1989), generalized information criterion (GIC; Konishi and Kitagawa
1996) and deviance information criterion (DIC; Spiegelhalter et al. 2002),
while definitely breakthroughs, do not give significantly different analytical
results from those of AIC. They are equivalent to AIC in special cases and are
significant in their expanded frameworks. Moreover, in propensity score
analysis, Platt et al. (2013) devised an information criterion that sets the
goodness-of-fit term to twice the inverse-probability-weighted log-likelihood
and the penalty term to twice the number of parameters, which was thought to
be reasonable for model selection. However, while the goodness-of-fit term is
appropriate, the penalty term is considerably underestimated in terms of bias
correction, and the penalty term proposed in this paper differs significantly
from it. Note that while the criterion proposed in Rolling and Yang (2014) is
also a strong choice in causal inference if we do not use the propensity score
as an inverse probability, this paper places more importance on the
semiparametric approach.

The remainder of this paper is organized as follows. Section 2 explains the
model and assumptions and introduces semiparametric estimation using
propensity scores and sparse estimation including SURE theory. In Section 3,
SURE theory is extended for propensity score analysis, and a criterion not
relying on asymptotics (IPCp) is developed. In Section 4, the asymptotic
theory of LASSO is extended for propensity score analysis, and for general
causal inference models that are not necessarily Gaussian distribution-based,
we develop an information criterion for inverse-probability-weighted sparse
estimation (IPIC). Subsequently, for doubly robust sparse estimation, we
develop DRIC in which the criterion itself is also doubly robust. In Section
5, numerical experiments are conducted to examine the approximation accuracy
of the proposed criterion to the divergence, and then, the proposed criterion
is compared with an existing criterion by measuring the mean square error and
the degree of variable selection. Section 6 describes a real data analysis
comparing the estimates given by the proposed and existing criteria in terms
of the magnitude of their differences. In Section 7, to explore the
possibilities for generalization of the proposed criterion, we first attempt
to derive IPCp for group LASSO and elastic net. Next, we attempt to derive
IPIC for general sparse estimation using nonconvex regularization. Finally, we
summarize our conclusions in Section 8.

##  2 Preparation

###  2.1 Model and assumption

First, we treat

$\displaystyle
y=\sum_{h=1}^{H}t^{[h]}y^{[h]}=\sum_{h=1}^{H}t^{[h]}\\{\bm{x}^{{\rm
T}}\bm{\theta}^{[h]}+\eta(\bm{z})+\varepsilon^{[h]}\\}$ (2)

as a generalized version of Rubin’s causal inference model. Here, there are
$H$ kinds of treatment, $t^{[h]}\ (\in\\{0,1\\})$ is an assignment variable
that is 1 when the $h$-th treatment is assigned ($\sum_{h=1}^{H}t^{[h]}=1$),
$y^{[h]}\ (\in\mathbb{R})$ is a potential outcome variable when the $h$-th
treatment is assigned, $\bm{x}\ (\in\mathbb{R}^{p})$ is an explanatory
variable, $\bm{\theta}^{[h]}\ (\in\mathbb{R}^{p})$ is a regression coefficient
parameter vector for $y^{[h]}$, $\bm{z}\ (\in\mathbb{R}^{q})$ is a confounding
variable affecting both $y^{[h]}$ and $t^{[h]}$,
$\eta:\mathbb{R}^{q}\to\mathbb{R}$ is an unknown nonparametric function, and
$\varepsilon^{[h]}\ (\in\mathbb{R})$ is an error variable distributed
according to ${\rm N}(0,\sigma^{2})$ independently of $t^{[h]}$ and $\bm{z}$
($h\in\\{1,2,\ldots,H\\}$). Note that $y$ on the left-hand side is an observed
outcome variable. We also assume that $\bm{x}$ may contain part of $\bm{z}$
and is otherwise independent of $\bm{z}$ for simplicity. Letting
$(c^{[1]},c^{[2]},\ldots,c^{[H]})$ be a contrast function satisfying
$\sum_{h=1}^{H}c^{[h]}=0$ and $\sum_{h=1}^{H}c^{[h]2}=1$, we target the causal
effect defined by $\bm{\theta}\equiv\sum_{h=1}^{H}c^{[h]}\bm{\theta}^{[h]}$
for estimation. The true value $\bm{\theta}^{*}$ of this is regarded as
$\bm{\theta}$, which is obtained by solving

$\displaystyle{\rm E}(\bm{x}\bm{x}^{{\rm
T}})\bm{\theta}-\sum_{h=1}^{H}c^{[h]}{\rm E}(\bm{x}y^{[h]})=\bm{0}_{p}.$

If $H=2$ and $(c^{[1]},c^{[2]})=(-1/\sqrt{2},1/\sqrt{2})$, then $\bm{\theta}$
is the familiar causal effect of $\bm{\theta}^{[2]}-\bm{\theta}^{[1]}$ divided
by $\sqrt{2}$.

Furthermore, we consider

$\displaystyle y=\sum_{h=1}^{H}t^{[h]}y^{[h]},\qquad y^{[h]}\sim
f(\cdot\mid\bm{x}^{[h]};\bm{\beta})$ (3)

as a generalized causal inference model of (2). Here, the notation is changed
from (2), where the explanatory variables are expressed for each $y^{[h]}$, as
in $\bm{x}^{[h]}$, and the parameters of the regression coefficients are
grouped for all groups, as in $\bm{\beta}$. For $f:\mathbb{R}\to\mathbb{R}$,
we treat a parametric function that is convex and differentiable with respect
to $\bm{\beta}$, for which the regularity conditions are appropriately
satisfied, such as for a generalized linear model with a natural link
function. Its true value $\bm{\beta}^{*}$ is expected to satisfy

$\displaystyle{\rm
E}\Bigg{\\{}\sum_{h=1}^{H}\frac{\partial}{\partial\bm{\beta}}\log
f(y^{[h]}\mid\bm{x}^{[h]};\bm{\beta}^{*})\Bigg{\\}}=\bm{0}_{p}.$ (4)

In these models, $H-1$ potential outcome variables $y^{[h]}$ with $t^{[h]}=0$
are regarded as missing. In general, since ${\rm E}(y^{[h]})\neq{\rm
E}(y^{[h]}\mid t^{[h]}=1)$, naively estimating $\bm{\theta}$ from the observed
values will result in a bias. Here, we suppose that a confounding variable
$\bm{z}$ affecting $y^{[h]}$ and $t^{[h]}$ is observed such that this bias can
be eliminated. In addition, we assume an ignorable treatment assignment
condition,

$\displaystyle\\{y^{[1]},y^{[2]},\ldots,y^{[H]}\\}\mathop{\perp\\!\\!\\!\\!\perp}\\{t^{[1]},t^{[2]},\ldots,t^{[H]}\\}\mid\bm{z},$

and a positive value condition ${\rm P}(t^{[h]}=1)>0$, which are what make the
elimination possible. There are $N$ samples subject to this model, and we will
subscript the variables in the $i$-th sample with $i$. We assume that the
samples are independent, i.e.,

$\displaystyle(y_{i}^{[1]},y_{i}^{[2]},\ldots,y_{i}^{[H]},t_{i}^{[1]},t_{i}^{[2]},\ldots,t_{i}^{[H]},\bm{x}_{i},\bm{z}_{i})\mathop{\perp\\!\\!\\!\\!\perp}(y_{i^{\prime}}^{[1]},y_{i^{\prime}}^{[2]},\ldots,y_{i^{\prime}}^{[H]},t_{i^{\prime}}^{[1]},t_{i^{\prime}}^{[2]},\ldots,t_{i^{\prime}}^{[H]},\bm{x}_{i^{\prime}},\bm{z}_{i^{\prime}})$

for $i\neq i^{\prime}\ (i,i^{\prime}\in\\{1,2,\ldots,N\\})$. It follows
automatically that $y_{i}\mathop{\perp\\!\\!\\!\\!\perp}y_{i^{\prime}}$.

###  2.2 Semiparametric estimation

If the relationship between the potential outcome variable $y^{[h]}$ and the
confounding variable $\bm{z}$ can be modeled correctly, we can make valid
estimates of the causal effect under the ignorable treatment assignment
condition; however, this modeling is generally difficult. In recent years, a
semiparametric approach using the propensity score $e^{[h]}(\bm{z})\equiv{\rm
P}(t^{[h]}=1\mid\bm{z})$ has been often used, which does not necessarily
require this correct modeling. This subsection describes two such estimation
methods without sparse estimation.

The first type is inverse-probability-weighted estimation (Robins et al.
1994). In this method, the observed values are multiplied by the inverse of
the propensity scores as weights to pseudo-recover the missing values; then,
an ordinary estimation is conducted. Specifically, for (2), by solving

$\displaystyle\frac{\partial}{\partial\bm{\theta}}\Bigg{[}\sum_{i=1}^{N}\Bigg{\\{}\sum_{h=1}^{H}\frac{c^{[h]}t_{i}^{[h]}y_{i}^{[h]}}{e^{[h]}(\bm{z}_{i})}-\bm{x}_{i}^{{\rm
T}}\bm{\theta}\Bigg{\\}}^{2}\Bigg{]}=\sum_{i=1}^{N}\bm{x}_{i}\Bigg{\\{}\bm{x}_{i}^{{\rm
T}}\bm{\theta}-\sum_{h=1}^{H}\frac{c^{[h]}t_{i}^{[h]}y_{i}^{[h]}}{e^{[h]}(\bm{z}_{i})}\Bigg{\\}}=\bm{0}_{p},$

we obtain the inverse-probability-weighted estimator,

$\displaystyle\hat{\bm{\theta}}=\Bigg{(}\sum_{i=1}^{N}\bm{x}_{i}\bm{x}_{i}^{{\rm
T}}\Bigg{)}^{-1}\sum_{i=1}^{N}\sum_{h=1}^{H}\frac{c^{[h]}t_{i}^{[h]}y_{i}^{[h]}\bm{x}_{i}}{e^{[h]}(\bm{z}_{i})}.$

Under the ignorable treatment assignment condition, if $\bm{z}_{i}$ is
conditioned, $t_{i}^{[h]}$ is independent of $y_{i}^{[h]}$ and has expectation
$e^{[h]}(\bm{z}_{i})$, and therefore it follows that $\hat{\bm{\theta}}$ is
unbiased for $\bm{\theta}^{*}$. On the other hand, for (3), the inverse-
probability-weighted estimator $\hat{\bm{\beta}}$ is obtained by solving

$\displaystyle\frac{\partial}{\partial\bm{\beta}}\Bigg{\\{}\frac{1}{N}\sum_{i=1}^{N}\sum_{h=1}^{H}\frac{t_{i}^{[h]}}{e^{[h]}(\bm{z}_{i})}\log
f(y_{i}^{[h]}\mid\bm{x}_{i}^{[h]};\bm{\beta})\Bigg{\\}}=\bm{0}_{p}.$ (5)

Under the ignorable treatment assignment condition, the left-hand side of (5)
with $\bm{\beta}=\bm{\beta}^{*}$ converges to the left-hand side of (4), and,
therefore, the inverse-probability-weighted estimator satisfies that
$\hat{\bm{\beta}}\stackrel{{\scriptstyle\rm p}}{{\to}}\bm{\beta}^{*}$. By
substituting $\hat{\bm{\beta}}$ for the left-hand side of (5) and using the
expansion, we obtain

$\displaystyle\hat{\bm{\beta}}-\bm{\beta}^{*}=\bm{J}(\bm{\beta}^{*})^{-1}\frac{1}{N}\sum_{i=1}^{N}\sum_{h=1}^{H}\frac{t_{i}^{[h]}}{e^{[h]}(\bm{z}_{i})}\frac{\partial}{\partial\bm{\beta}}\log
f(y_{i}^{[h]}\mid\bm{x}_{i}^{[h]};\bm{\beta}^{*})\\{1+{\rm o}_{\rm P}(1)\\},$

where

$\displaystyle\bm{J}(\bm{\beta}^{*})\equiv{\rm
E}\Bigg{\\{}-\sum_{h=1}^{H}\frac{\partial^{2}}{\partial\bm{\beta}\partial\bm{\beta}^{{\rm
T}}}\log f(y^{[h]}\mid\bm{x}^{[h]};\bm{\beta}^{*})\Bigg{\\}}.$ (6)

In the above, although the propensity score is assumed to be known, it is
often unknown in practice. In that case, a parametric function
$e^{[h]}(\bm{z};\bm{\alpha})$ with parameter $\bm{\alpha}\
(\in\mathbb{R}^{q})$ is used to model ${\rm P}(t^{[h]}=1\mid\bm{z})$ and some
estimator $\hat{\bm{\alpha}}$ is substituted into $\bm{\alpha}$. For
$\hat{\bm{\alpha}}$, for example, we only need to use the maximum likelihood
estimator after constructing the likelihood for $t^{[h]}$ from the multinomial
distribution of probability $e^{[h]}(\bm{z};\bm{\alpha})$. Although $y^{[h]}$
is correlated with $\bm{z}$ in the causal inference model, the inverse-
probability-weighted estimation does not directly use information on $\bm{z}$.
The doubly robust estimation (Scharfstein et al. 1999, Bang and Robins 2005)
is an improved version of the inverse-probability-weighted estimation.
Denoting the conditional distribution of $y^{[h]}$ given $\bm{z}$ as
$f^{[h]}(y^{[h]}\mid\bm{z};\bm{\gamma})$ with parameter $\bm{\gamma}\
(\in\mathbb{R}^{r})$, we will use the expectation of $\log
f(y^{[h]}\mid\bm{x}^{[h]};\bm{\beta})$ taken by
$f^{[h]}(y^{[h]}\mid\bm{z};\bm{\gamma})$. In practice, a consistent estimator
$\hat{\bm{\gamma}}$ is substituted for $\bm{\gamma}$. Specifically, the doubly
robust estimator $\hat{\bm{\beta}}$ is obtained by adding

$\displaystyle\frac{\partial}{\partial\bm{\beta}}\Bigg{[}\frac{1}{N}\sum_{i=1}^{N}\sum_{h=1}^{H}\Bigg{\\{}1-\frac{t_{i}^{[h]}}{e^{[h]}(\bm{z}_{i};\hat{\bm{\alpha}})}\Bigg{\\}}{\rm
E}\\{\log
f(y^{[h]}\mid\bm{x}^{[h]};\bm{\beta})\mid\bm{z};\hat{\bm{\gamma}}\\}\Bigg{]}$

to the left-hand side of (5) and solving for $\bm{\beta}$ such that the sum of
the sum equals $\bm{0}_{p}$, in the setting where the propensity score is
unknown. Not only does this estimator actually improve the inverse-
probability-weighted estimator; it is also asymptotically semiparametric
efficient (Robins and Rotnitzky 1995). In addition, if either the propensity
score or the conditional expectation is modeled correctly, it is consistent
and hence said to be doubly robust.

###  2.3 Sparse estimation

In this subsection, we treat a simple model,

$\displaystyle y_{i}^{[h]}=\bm{x}_{i}^{{\rm
T}}\bm{\theta}^{[h]}+\varepsilon_{i}^{[h]},\qquad\varepsilon_{i}^{[h]}\sim{\rm
N}(0,\sigma^{2}),$

in order to introduce Stein (1981)’s unbiased risk estimation (SURE) theory
($i\in\\{1,2,\ldots,N\\}$, $h\in\\{1,2,\ldots,H\\}$). The number of candidate
explanatory variables is also $p$, i.e.,
$\bm{x}_{i},\bm{\theta}^{[h]}\in\mathbb{R}^{p}$ in the above. Tibshirani
(1996) proposed LASSO as for a way of selecting explanatory variables and
estimating parameters at the same time. In this model, if we want to estimate
the causal effect $\bm{\theta}\equiv\sum_{h=1}^{H}c^{[h]}\bm{\theta}^{[h]}$,
LASSO gives

$\displaystyle\hat{\bm{\theta}}_{\lambda}\equiv\mathop{\rm
argmin}_{\bm{\theta}}\Bigg{\\{}\sum_{i=1}^{N}\Bigg{(}\sum_{h=1}^{H}c^{[h]}y_{i}^{[h]}-\bm{x}_{i}^{{\rm
T}}\bm{\theta}\Bigg{)}^{2}+\lambda\|\bm{\theta}\|_{1}\Bigg{\\}}.$

Here, $\|\cdot\|_{1}$ is the $\ell_{1}$ norm, i.e.
$\|\bm{\theta}\|_{1}=\sum_{j=1}^{p}|\theta_{j}|$. In addition, $\lambda$ is a
regularization parameter, and the larger it is, the stronger the effect of
$\|\bm{\theta}\|_{1}$ becomes and the more $j$’s there are such that
$\hat{\theta}_{\lambda,j}$ is strictly $0$, corresponding to a stronger
reduction in the variables. Therefore, determining the size of $\lambda$
implies model selection, which is an essential issue. The collection of $j$
such that $\hat{\theta}_{\lambda,j}$ is not equal to 0 is called the active
set, and the number of elements is denoted by $\hat{p}$.

Since $\|\bm{\theta}\|_{1}$ has non-differentiable points, the methods used to
derive the Cp criterion and AIC cannot be applied. For such a situation, Efron
et al. (2004) and Zou et al. (2007) give an unbiased estimator of the mean
squared error for the LASSO estimator minus a constant, which can be regarded
as a generalized Cp criterion, in an explicit form. They argue that a
reasonable $\lambda$ can be selected by minimizing it. As with the
conventional Cp criterion, the mean squared error consists of three terms:

$\displaystyle{\rm E}\Bigg{[}\sum_{i=1}^{N}\Bigg{\\{}\bm{x}_{i}^{{\rm
T}}\hat{\bm{\theta}}_{\lambda}-{\rm
E}\Bigg{(}\sum_{h=1}^{H}c^{[h]}y_{i}^{[h]}\Bigg{)}\Bigg{\\}}^{2}\Bigg{]}$
$\displaystyle={\rm
E}\Bigg{\\{}\sum_{i=1}^{N}\Bigg{(}\sum_{h=1}^{H}c^{[h]}y_{i}^{[h]}-\bm{x}_{i}^{{\rm
T}}\hat{\bm{\theta}}_{\lambda}\Bigg{)}^{2}\Bigg{\\}}-{\rm
E}\Bigg{[}\sum_{i=1}^{N}\Bigg{\\{}\sum_{h=1}^{H}c^{[h]}y_{i}^{[h]}-{\rm
E}\Bigg{(}\sum_{h=1}^{H}c^{[h]}y_{i}^{[h]}\Bigg{)}\Bigg{\\}}^{2}\Bigg{]}$
$\displaystyle\ \phantom{=}+{\rm
E}\Bigg{[}2\sum_{i=1}^{N}\Bigg{\\{}\sum_{h=1}^{H}c^{[h]}y_{i}^{[h]}-{\rm
E}\Bigg{(}\sum_{h=1}^{H}c^{[h]}y_{i}^{[h]}\Bigg{)}\Bigg{\\}}\Bigg{\\{}\bm{x}_{i}^{{\rm
T}}\hat{\bm{\theta}}_{\lambda}-{\rm
E}\Bigg{(}\sum_{h=1}^{H}c^{[h]}y_{i}^{[h]}\Bigg{)}\Bigg{\\}}\Bigg{]}.$

In the first term, the expectation is removed because its content does not
depend on the parameters. The second term can be ignored because it does not
depend on the model. The third term reduces to

$\displaystyle{\rm
E}\Bigg{[}2\sum_{i=1}^{N}\Bigg{(}\sum_{h=1}^{H}c^{[h]}\varepsilon_{i}^{[h]}\Bigg{)}\Bigg{\\{}\bm{x}_{i}^{{\rm
T}}\hat{\bm{\theta}}_{\lambda}-{\rm
E}\Bigg{(}\sum_{h=1}^{H}c^{[h]}y_{i}^{[h]}\Bigg{)}\Bigg{\\}}\Bigg{]}={\rm
E}\Bigg{(}2\sum_{i=1}^{N}\sum_{h=1}^{H}c^{[h]}\varepsilon_{i}^{[h]}\bm{x}_{i}^{{\rm
T}}\hat{\bm{\theta}}_{\lambda}\Bigg{)},$

and applying SURE theory yields

$\displaystyle{\rm
E}\Bigg{(}2\sigma^{2}\sum_{i=1}^{N}\sum_{h=1}^{H}c^{[h]}\bm{x}_{i}^{{\rm
T}}\frac{\partial\hat{\bm{\theta}}_{\lambda}}{\partial
y_{i}^{[h]}}\Bigg{)}={\rm
E}\Bigg{\\{}2\sigma^{2}\sum_{i=1}^{N}\sum_{h=1}^{H}c^{[h]}\bm{x}_{i}^{(2){\rm
T}}\Bigg{(}\sum_{i^{\prime}=1}^{N}\bm{x}_{i^{\prime}}^{(2)}\bm{x}_{i^{\prime}}^{(2){\rm
T}}\Bigg{)}^{-1}\bm{x}_{i}^{(2)}c^{[h]}\Bigg{\\}}={\rm
E}(2\hat{p}\sigma^{2}).$

This means that

$\displaystyle\sum_{i=1}^{N}\Bigg{(}\sum_{h=1}^{H}c^{[h]}y_{i}^{[h]}-\bm{x}_{i}^{{\rm
T}}\hat{\bm{\theta}}_{\lambda}\Bigg{)}^{2}+2\hat{p}\sigma^{2}$

is an unbiased estimator of the mean squared error minus a constant, and the
value of $\lambda$ that minimizes it should be selected as the optimal one.
Since $\hat{p}$ is the number of explanatory variables actually used in the
prediction, this generalized Cp criterion seems like the conventional Cp
criterion. However, the penalty term is likely to be larger because LASSO
considers all explanatory variables, while the penalty term is likely to be
smaller because LASSO is a shrinkage estimation, and the two just cancel and
$\hat{p}$ appears in the penalty (Lockhart et al. 2014). When $\sigma^{2}$ is
unknown, we only need to substitute the estimator $\hat{\sigma}^{2}$ for
$\sigma^{2}$ as in the conventional Cp criterion.

##  3 Information criterion based on SURE theory

For the model in (2), it is natural to consider the LASSO estimator,

$\displaystyle\hat{\bm{\theta}}_{\lambda}\equiv\mathop{\rm
argmin}_{\bm{\theta}}\Bigg{[}\sum_{i=1}^{N}\Bigg{\\{}\sum_{h=1}^{H}\frac{c^{[h]}t_{i}^{[h]}y_{i}^{[h]}}{e^{[h]}(\bm{z}_{i})}-\bm{x}_{i}^{{\rm
T}}\bm{\theta}\Bigg{\\}}^{2}+\lambda\|\bm{\theta}\|_{1}\Bigg{]},$

if the explanatory variables are of somewhat high dimensional. Here, the set
of the subscripts with zero components is denoted as
$\hat{\mathcal{J}}_{\lambda}^{(1)}$ and that with non-zero components is
denoted as $\hat{\mathcal{J}}_{\lambda}^{(2)}$. According to this notation,
$(\hat{\theta}_{\lambda,j})_{j\in\hat{\mathcal{J}}_{\lambda}^{(1)}}$ is
denoted as $\hat{\bm{\theta}}_{\lambda}^{(1)}$,
$(\hat{\theta}_{\lambda,j})_{j\in\hat{\mathcal{J}}_{\lambda}^{(2)}}$ is
denoted as $\hat{\bm{\theta}}_{\lambda}^{(2)}$,
$(x_{i,j})_{j\in\hat{\mathcal{J}}_{\lambda}^{(2)}}$ is denoted as
$\bm{x}_{i}^{(2)}$, and so on. Note that $\hat{\mathcal{J}}_{\lambda}^{(2)}$
is the active set. Of course, $\hat{\bm{\theta}}_{\lambda}^{(1)}$ is
$\bm{0}_{|\hat{\mathcal{J}}_{\lambda}^{(1)}|}$, and we see that
$\hat{\bm{\theta}}_{\lambda}^{(2)}$ satisfies

$\displaystyle\hat{\bm{\theta}}_{\lambda}^{(2)}=\Bigg{(}\sum_{i=1}^{N}\bm{x}_{i}^{(2)}\bm{x}_{i}^{(2){\rm
T}}\Bigg{)}^{-1}\Bigg{\\{}\sum_{i=1}^{N}\sum_{h=1}^{H}\frac{c^{[h]}t_{i}^{[h]}y_{i}^{[h]}\bm{x}_{i}^{(2)}}{e^{[h]}(\bm{z}_{i})}+\lambda\hskip
1.42262pt{\rm sgn}(\hat{\bm{\theta}}_{\lambda}^{(2)})\Bigg{\\}}$ (7)

by differentiation. Now let us derive a generalized Cp criterion to give an
appropriate $\lambda$. Existing SURE theory, which targets the mean squared
error, requires the data to follow a Gaussian distribution; this may seem
impossible in our setting where $y^{[h]}$ is non-Gaussian. However, by
conditioning the assignment variable $t^{[h]}$ earlier, which is not usually
done in propensity score analysis, SURE theory can be applied because of the
conditional Gaussianity.

Before defining the risk function for deriving the information criterion in
causal inference, we describe QICw proposed in Platt et al. (2013). This
criterion is for when there are missing data, such as potential outcome
variables, and twice the number of parameters is used as the second term. The
first term is defined as follows: First, the log-likelihood function of the
complete data without missing is differentiated, the terms in which the
outcome variable $y_{i}^{[h]}$ appears are multiplied by
$t_{i}^{[h]}/e^{[h]}(\bm{z}_{i})$, and the derivative is restored by
integration. Next, the estimator is substituted for the parameters and
doubled. Although it is not treated in the framework of sparse estimation, we
write

$\displaystyle{\rm
QICw}\equiv\frac{1}{\sigma^{2}}\sum_{i=1}^{N}\Bigg{\\{}\sum_{h=1}^{H}\frac{c^{[h]}t_{i}^{[h]}y_{i}^{[h]}}{e^{[h]}(\bm{z}_{i})}-\bm{x}_{i}^{(2){\rm
T}}\hat{\bm{\theta}}_{\lambda}^{(2)}\Bigg{\\}}^{2}+2|\hat{\mathcal{J}}_{\lambda}^{(2)}|$

according to this definition and use it for comparison. In this model, if the
inverse-probability $1/e^{[h]}(\bm{z}_{i})$ is removed by supposing no
confounding, we can confirm that QICw is an unbiased estimator of the mean
squared error. However, if confounding exists and the inverse-probability
$1/e^{[h]}(\bm{z}_{i})$ is added, the second term must be larger because of
the larger variation of the first term.

As a risk function for all groups combined, i.e., one based on the loss
function used in estimation, we consider the mean squared error,

$\displaystyle{\rm E}\Bigg{(}\sum_{i=1}^{N}\Bigg{[}\sum_{h=1}^{H}{\rm
E}\Bigg{\\{}\frac{c^{[h]}t_{i}^{[h]}y_{i}^{[h]}}{e^{[h]}(\bm{z}_{i})}\
\Bigg{|}\ t_{i}^{[h]},\bm{x}_{i},\bm{z}_{i}\Bigg{\\}}-\bm{x}_{i}^{{\rm
T}}\hat{\bm{\theta}}_{\lambda}\Bigg{]}^{2}\Bigg{)}.$ (8)

If $t_{i}^{[h]}$ is not given in this conditional expectation, the sum becomes
$\sum_{h=1}^{H}{\rm
E}(c^{[h]}y_{i}^{[h]}\mid\bm{x}_{i},\bm{z}_{i})=\sum_{h=1}^{H}c^{[h]}\\{\bm{x}_{i}^{{\rm
T}}\bm{\theta}^{[h]}+\eta(\bm{z}_{i})\\}=\bm{x}_{i}^{{\rm T}}\bm{\theta}$,
which may be easy to understand. Note, however, that we are considering a risk
function that is more tailored to actual data by giving $t_{i}^{[h]}$. By
inserting $\sum_{h=1}^{H}c^{[h]}t_{i}^{[h]}y_{i}^{[h]}/e^{[h]}(\bm{z}_{i})$
between the first and second terms, the risk function can be decomposed into
three terms:

$\displaystyle{\rm
E}\Bigg{[}\sum_{i=1}^{N}\Bigg{\\{}\sum_{h=1}^{H}\frac{c^{[h]}t_{i}^{[h]}y_{i}^{[h]}}{e^{[h]}(\bm{z}_{i})}-\bm{x}_{i}^{{\rm
T}}\hat{\bm{\theta}}_{\lambda}\Bigg{\\}}^{2}\Bigg{]}$ $\displaystyle-{\rm
E}\Bigg{(}\sum_{i=1}^{N}\Bigg{[}\sum_{h=1}^{H}\frac{c^{[h]}t_{i}^{[h]}y_{i}^{[h]}}{e^{[h]}(\bm{z}_{i})}-\sum_{h=1}^{H}{\rm
E}\Bigg{\\{}\frac{c^{[h]}t_{i}^{[h]}y_{i}^{[h]}}{e^{[h]}(\bm{z}_{i})}\
\Bigg{|}\ t_{i}^{[h]},\bm{x}_{i},\bm{z}_{i}\Bigg{\\}}\Bigg{]}^{2}\Bigg{)}$
$\displaystyle+2{\rm
E}\Bigg{(}\sum_{i=1}^{N}\Bigg{[}\sum_{h=1}^{H}\frac{c^{[h]}t_{i}^{[h]}y_{i}^{[h]}}{e^{[h]}(\bm{z}_{i})}-\sum_{h=1}^{H}{\rm
E}\Bigg{\\{}\frac{c^{[h]}t_{i}^{[h]}y_{i}^{[h]}}{e^{[h]}(\bm{z}_{i})}\
\Bigg{|}\ t_{i}^{[h]},\bm{x}_{i},\bm{z}_{i}\Bigg{\\}}\Bigg{]}\bm{x}_{i}^{{\rm
T}}\hat{\bm{\theta}}_{\lambda}\Bigg{)}.$

In the first term of this decomposition, the expectation is removed because
its content does not depend on the parameters, while the second term can be
ignored because it does not depend on the model. Since the expectation of the
third term can be expressed as

$\displaystyle{\rm
E}\Bigg{(}\sum_{i=1}^{N}\Bigg{[}\sum_{h=1}^{H}\frac{c^{[h]}t_{i}^{[h]}y_{i}^{[h]}}{e^{[h]}(\bm{z}_{i})}-\sum_{h=1}^{H}\frac{c^{[h]}t_{i}^{[h]}\\{\bm{x}_{i}^{{\rm
T}}\bm{\theta}^{[h]}+\eta(\bm{z}_{i})\\}}{e^{[h]}(\bm{z}_{i})}\Bigg{]}\bm{x}_{i}^{{\rm
T}}\hat{\bm{\theta}}_{\lambda}\Bigg{)}$ $\displaystyle={\rm
E}\Bigg{[}\sum_{i=1}^{N}\sum_{h=1}^{H}\frac{c^{[h]}t_{i}^{[h]}}{e^{[h]}(\bm{z}_{i})}\\{y_{i}^{[h]}-\bm{x}_{i}^{{\rm
T}}\bm{\theta}^{[h]}-\eta(\bm{z}_{i})\\}\bm{x}_{i}^{{\rm
T}}\hat{\bm{\theta}}_{\lambda}\Bigg{]},$

by conditioning on $(t_{i}^{[h]},\bm{x}_{i},\bm{z}_{i})$ earlier, we see that
SURE theory can be applied. This expectation is evaluated as

$\displaystyle{\rm E}\Bigg{[}\sum_{i=1}^{N}\sum_{h=1}^{H}{\rm
E}\Bigg{\\{}\frac{c^{[h]}t_{i}^{[h]}}{e^{[h]}(\bm{z}_{i})}\sigma^{2}\bm{x}_{i}^{{\rm
T}}\frac{\partial\hat{\bm{\theta}}_{\lambda}}{\partial y_{i}^{[h]}}\ \Bigg{|}\
t_{i}^{[h]},\bm{x}_{i},\bm{z}_{i}\Bigg{\\}}\Bigg{]}$ $\displaystyle={\rm
E}\Bigg{(}\sum_{i=1}^{N}\sum_{h=1}^{H}{\rm
E}\Bigg{[}\frac{c^{[h]}t_{i}^{[h]}}{e^{[h]}(\bm{z}_{i})}\sigma^{2}\bm{x}_{i}^{(2){\rm
T}}\Bigg{(}\sum_{i^{\prime}=1}^{N}\bm{x}_{i^{\prime}}^{(2)}\bm{x}_{i^{\prime}}^{(2){\rm
T}}\Bigg{)}^{-1}\frac{c^{[h]}t_{i}^{[h]}}{e^{[h]}(\bm{z}_{i})}\bm{x}_{i}^{(2)}\
\Bigg{|}\ t_{i}^{[h]},\bm{x}_{i},\bm{z}_{i}\Bigg{]}\Bigg{)}.$

If the expectation is taken only for $t_{i}^{[h]}$, then, given $\bm{z}_{i}$,
we see that $t_{i}^{[h]2}=t_{i}^{[h]}$ follows a Bernoulli distribution with
probability $e^{[h]}(\bm{z}_{i})$. This expectation can be further evaluated
as

$\displaystyle{\rm E}\Bigg{(}\sigma^{2}{\rm
tr}\Bigg{[}\sum_{i=1}^{N}\sum_{h=1}^{H}{\rm
E}\Bigg{\\{}\Bigg{(}\sum_{i^{\prime}=1}^{N}\bm{x}_{i^{\prime}}^{(2)}\bm{x}_{i^{\prime}}^{(2){\rm
T}}\Bigg{)}^{-1}\frac{c^{[h]2}}{e^{[h]}(\bm{z}_{i})}\bm{x}_{i}^{(2)}\bm{x}_{i}^{(2){\rm
T}}\ \Bigg{|}\ \bm{x}_{i},\bm{z}_{i}\Bigg{\\}}\Bigg{]}\Bigg{)}$
$\displaystyle={\rm E}\Bigg{(}\sigma^{2}{\rm
tr}\Bigg{[}\Bigg{(}\sum_{i=1}^{N}\bm{x}_{i}^{(2)}\bm{x}_{i}^{(2){\rm
T}}\Bigg{)}^{-1}\Bigg{\\{}\sum_{i=1}^{N}\sum_{h=1}^{H}\frac{c^{[h]2}}{e^{[h]}(\bm{z}_{i})}\bm{x}_{i}^{(2)}\bm{x}_{i}^{(2){\rm
T}}\Bigg{\\}}\Bigg{]}\Bigg{)}.$

This leads us to the following theorem.

######  Theorem 1.

In the model given by (2),

$\displaystyle{\rm IPCp}\equiv\ $
$\displaystyle\sum_{i=1}^{N}\Bigg{\\{}\sum_{h=1}^{H}\frac{c^{[h]}t_{i}^{[h]}y_{i}^{[h]}}{e^{[h]}(\bm{z}_{i})}-\bm{x}_{i}^{(2){\rm
T}}\hat{\bm{\theta}}_{\lambda}^{(2)}\Bigg{\\}}^{2}$
$\displaystyle+2\sigma^{2}{\rm
tr}\Bigg{[}\Bigg{(}\sum_{i=1}^{N}\bm{x}_{i}^{(2)}\bm{x}_{i}^{(2){\rm
T}}\Bigg{)}^{-1}\Bigg{\\{}\sum_{i=1}^{N}\sum_{h=1}^{H}\frac{c^{[h]2}}{e^{[h]}(\bm{z}_{i})}\bm{x}_{i}^{(2)}\bm{x}_{i}^{(2){\rm
T}}\Bigg{\\}}\Bigg{]}$

is an unbiased estimator of the mean squared error in (8).

This theorem does not rely on asymptotics, so the information criterion (IPCp)
is valid even when the dimension $p$ is large. If we set
$e^{[h]}(\bm{z}_{i})=1$ in the penalty term, the content of the trace becomes
the identity matrix and we obtain the conventional generalized Cp criterion,
or, in other words, we obtain QICw by dividing both sides by $\sigma^{2}$. In
reality, however, the penalty term becomes several times larger than that
because $e^{[h]}(\bm{z}_{i})$ is attached.

##  4 Information criterion based on asymptotics

###  4.1 AIC for inverse-probability-weighted LASSO

Similar to the AIC for LASSO in Ninomiya and Kawano (2016), for the
generalized causal inference model in (3), the weak limit $b^{\rm limit}$,
which is the source of the asymptotic bias, can be obtained from the
asymptotic properties of the semiparametric estimator and an information
criterion can be derived. Supposing that the dimension $p$ of the explanatory
variable $\bm{x}$ is still somewhat high and that the propensity score is
known. Later, we will consider the case where the propensity score is unknown.
The LASSO estimator is given by

$\displaystyle\hat{\bm{\beta}}_{\lambda}\equiv\mathop{\rm
argmin}_{\bm{\beta}}\Bigg{\\{}-\frac{1}{N}\sum_{i=1}^{N}\sum_{h=1}^{H}\frac{t_{i}^{[h]}}{e^{[h]}(\bm{z}_{i})}\log
f(y_{i}^{[h]}\mid\bm{x}_{i}^{[h]};\bm{\beta})+\lambda\|\bm{\beta}\|_{1}\Bigg{\\}}.$

Note that what we are going to do is to derive an information criterion by
relying on asymptotics and the order of the $\ell_{1}$ penalty is set to ${\rm
O}(1)$ in order to faithfully capture the behavior of
$\hat{\bm{\beta}}_{\lambda}$. To demonstrate the validity of a sparse
estimator such as an oracle property, its order is often set to ${\rm o}(1)$.
If we define

$\displaystyle\bm{\beta}_{\lambda}^{*}\equiv\mathop{\rm
argmin}_{\bm{\beta}}\Bigg{[}-{\rm E}\Bigg{\\{}\sum_{h=1}^{H}\log
f(y^{[h]}\mid\bm{x}^{[h]};\bm{\beta})\Bigg{\\}}+\lambda\|\bm{\beta}\|_{1}\Bigg{]},$
(9)

the convex lemma in Andersen and Gill (1982) and Pollard (1991) shows that
$\hat{\bm{\beta}}_{\lambda}$ converges to $\bm{\beta}_{\lambda}^{*}$ in
probability. In this $\bm{\beta}_{\lambda}^{*}$, the set of the subscripts
with $0$ components is denoted as $\mathcal{J}_{\lambda}^{*(1)}$ and that
without $0$ components is denoted as $\mathcal{J}_{\lambda}^{*(2)}$.
Furthermore, in this subsection,
$(\hat{\beta}_{\lambda,j})_{j\in\mathcal{J}_{\lambda}^{*(1)}}$ is written as
$\hat{\bm{\beta}}_{\lambda}^{(1)}$,
$(\hat{\beta}_{\lambda,j})_{j\in\mathcal{J}_{\lambda}^{*(2)}}$ as
$\hat{\bm{\beta}}_{\lambda}^{(2)}$,
$({\beta}_{\lambda,j}^{*})_{j\in\mathcal{J}_{\lambda}^{*(2)}}$ as
${\bm{\beta}}_{\lambda}^{*(2)}$,
$(J_{jj^{\prime}})_{j,j^{\prime}\in\mathcal{J}_{\lambda}^{*(2)}}$ as
$\bm{J}^{(22)}$, and so on. Since the objective function in (9) is convex and
differentiable, we have

$\displaystyle
j\in\mathcal{J}_{\lambda}^{*(1)}\quad\Leftrightarrow\quad\beta_{\lambda,j}^{*}=0\quad\Rightarrow\quad-\lambda<\frac{\partial}{\partial\beta_{j}}{\rm
E}\Bigg{\\{}\sum_{h=1}^{H}\log
f(y^{[h]}\mid\bm{x}^{[h]};\bm{\beta}_{\lambda}^{*})\Bigg{\\}}<\lambda$ (10)

and

$\displaystyle
j\in\mathcal{J}_{\lambda}^{*(2)}\quad\Leftrightarrow\quad\beta_{\lambda,j}^{*}\neq
0\quad\Rightarrow\quad\frac{\partial}{\partial\beta_{j}}{\rm
E}\Bigg{\\{}\sum_{h=1}^{H}\log
f(y^{[h]}\mid\bm{x}^{[h]};\bm{\beta}_{\lambda}^{*})\Bigg{\\}}=\lambda\hskip
1.42262pt{\rm sgn}(\beta_{\lambda,j}^{*})$ (11)

from the Karush-Kuhn-Tucker condition. Then, from the convex lemma in Hjort
and Pollard (2011), we obtain

$\displaystyle N\hat{\bm{\beta}}_{\lambda}^{(1)}={\rm o}_{\rm P}(1)$ (12)

and

$\displaystyle\sqrt{N}(\hat{\bm{\beta}}_{\lambda}^{(2)}-\bm{\beta}_{\lambda}^{*(2)})$
$\displaystyle=\bm{J}^{(22)}(\bm{\beta}_{\lambda}^{*})^{-1}\Bigg{\\{}\frac{1}{\sqrt{N}}\sum_{i=1}^{N}\sum_{h=1}^{H}\frac{\partial}{\partial\bm{\beta}^{(2)}}\frac{t_{i}^{[h]}}{e^{[h]}(\bm{z}_{i})}\log
f(y_{i}^{[h]}\mid\bm{x}_{i}^{[h]};\bm{\beta}_{\lambda}^{*})-\sqrt{N}\lambda\hskip
1.42262pt{\rm sgn}(\bm{\beta}_{\lambda}^{*(2)})\Bigg{\\}}+{\rm o}_{\rm P}(1)$
(13)

from the same argument as in Knight and Fu (2000).

Expressing QICw proposed in Platt et al. (2013) as in the previous section
gives

$\displaystyle{\rm
QICw}\equiv-2\sum_{i=1}^{N}\sum_{h=1}^{H}\frac{t_{i}^{[h]}}{e^{[h]}(\bm{z}_{i})}\log
f(y_{i}^{[h]}\mid\bm{x}_{i}^{[h]};\hat{\bm{\beta}}_{\lambda})+2|\hat{\mathcal{J}}_{\lambda}^{(2)}|,$

where $\hat{\mathcal{J}}_{\lambda}^{(2)}$ is the active set for
$\hat{\bm{\beta}}_{\lambda}$. With the first term in mind, let us consider the
risk,

$\displaystyle-2{\rm
E}\Bigg{\\{}\sum_{i=1}^{N}\sum_{h=1}^{H}\frac{t_{i}^{\dagger[h]}}{e^{[h]}(\bm{z}_{i}^{\dagger})}\log
f(y_{i}^{\dagger[h]}\mid\bm{x}_{i}^{\dagger[h]};\hat{\bm{\beta}}_{\lambda})\Bigg{\\}},$

based on the weighted log-likelihood loss, as in Baba and Ninomiya (2021),
where
$(t_{i}^{\dagger[h]},y_{i}^{\dagger[h]},\bm{x}_{i}^{\dagger},\allowbreak\bm{z}_{i}^{\dagger})$
is a copy of $(t_{i}^{[h]},y_{i}^{[h]},\bm{x}_{i},\bm{z}_{i})$. This risk is
based on the same loss function used for the parameter estimation and can be
regarded as the risk when the population is the sum of all groups from $h=1$
to $h=H$. As in the usual derivation of AIC, it can be expanded into three
terms:

$\displaystyle-2{\rm
E}\Bigg{\\{}\sum_{i=1}^{N}\sum_{h=1}^{H}\frac{t_{i}^{[h]}}{e^{[h]}(\bm{z}_{i})}\log
f(y_{i}^{[h]}\mid\bm{x}_{i}^{[h]};\hat{\bm{\beta}}_{\lambda})\Bigg{\\}}$
$\displaystyle+2{\rm
E}\Bigg{[}\sum_{i=1}^{N}\sum_{h=1}^{H}\frac{t_{i}^{[h]}}{e^{[h]}(\bm{z}_{i})}\\{\log
f(y_{i}^{[h]}\mid\bm{x}_{i}^{[h]};\hat{\bm{\beta}}_{\lambda})-\log
f(y_{i}^{[h]}\mid\bm{x}_{i}^{[h]};\bm{\beta}_{\lambda}^{*})\\}\Bigg{]}$
$\displaystyle+2{\rm
E}\Bigg{[}\sum_{i=1}^{N}\sum_{h=1}^{H}\frac{t_{i}^{\dagger[h]}}{e^{[h]}(\bm{z}_{i}^{\dagger})}\\{\log
f(y_{i}^{\dagger[h]}\mid\bm{x}_{i}^{\dagger[h]};\bm{\beta}_{\lambda}^{*})-\log
f(y_{i}^{\dagger[h]}\mid\bm{x}_{i}^{\dagger[h]};\hat{\bm{\beta}}_{\lambda})\\}\Bigg{]}.$

We remove the expectation of the first term and asymptotically evaluate the
second and third terms, as usual, regarding the sum as AIC. For the asymptotic
evaluation, in the expectation, we perform a Taylor expansion with respect to
$\hat{\bm{\beta}}_{\lambda}$ around $\bm{\beta}_{\lambda}^{*}$, ignore ${\rm
o}_{\rm P}(1)$ terms, and then instead of
$\hat{\bm{\beta}}_{\lambda}-\bm{\beta}_{\lambda}^{*}$, we substitute
asymptotically what we consider to be the main term, i.e., the one obtained
from (12) and (13). Letting $b^{\rm limit}$ be the weak limit of the
substitution, we use ${\rm E}(b^{\rm limit})$, which is regarded as an
asymptotic bias, for the correction. Specifically, if we use

$\displaystyle\sum_{i=1}^{N}\sum_{h=1}^{H}\frac{t_{i}^{[h]}}{e^{[h]}(\bm{z}_{i})}\\{\log
f(y_{i}^{[h]}\mid\bm{x}_{i}^{[h]};\hat{\bm{\beta}}_{\lambda})-\log
f(y_{i}^{[h]}\mid\bm{x}_{i}^{[h]};\bm{\beta}_{\lambda}^{*})\\}-(\hat{\bm{\beta}}_{\lambda}^{(2)}-\bm{\beta}_{\lambda}^{*(2)})^{{\rm
T}}N\lambda\hskip 1.42262pt{\rm sgn}(\bm{\beta}_{\lambda}^{*(2)})$
$\displaystyle=\sum_{h=1}^{H}(\hat{\bm{\beta}}_{\lambda}^{(2)}-\bm{\beta}_{\lambda}^{*(2)})^{{\rm
T}}\Bigg{\\{}\sum_{i=1}^{N}\frac{\partial}{\partial\bm{\beta}^{(2)}}\frac{t_{i}^{[h]}}{e^{[h]}(\bm{z}_{i})}\log
f(y_{i}^{[h]}\mid\bm{x}_{i}^{[h]};\bm{\beta}_{\lambda}^{*})-N\lambda\hskip
1.42262pt{\rm sgn}(\bm{\beta}_{\lambda}^{*(2)})\Bigg{\\}}$ $\displaystyle\
\phantom{=}+\sum_{h=1}^{H}(\hat{\bm{\beta}}_{\lambda}^{(2)}-\bm{\beta}_{\lambda}^{*(2)})^{{\rm
T}}\Bigg{\\{}\sum_{i=1}^{N}\frac{\partial^{2}}{\partial\bm{\beta}^{(2)}\partial\bm{\beta}^{(2){\rm
T}}}\frac{t_{i}^{[h]}}{e^{[h]}(\bm{z}_{i})}\log
f(y_{i}^{[h]}\mid\bm{x}_{i}^{[h]};\bm{\beta}_{\lambda}^{*})\Bigg{\\}}(\hat{\bm{\beta}}_{\lambda}^{(2)}-\bm{\beta}_{\lambda}^{*(2)})$

for the second term and

$\displaystyle\sum_{i=1}^{N}\sum_{h=1}^{H}\frac{t_{i}^{\dagger[h]}}{e^{[h]}(\bm{z}_{i}^{\dagger})}\\{\log
f(y_{i}^{\dagger[h]}\mid\bm{x}_{i}^{\dagger[h]};\hat{\bm{\beta}}_{\lambda})-\log
f(y_{i}^{\dagger[h]}\mid\bm{x}_{i}^{\dagger[h]};\bm{\beta}_{\lambda}^{*})\\}-(\hat{\bm{\beta}}_{\lambda}^{(2)}-\bm{\beta}_{\lambda}^{*(2)})^{{\rm
T}}N\lambda\hskip 1.42262pt{\rm sgn}(\bm{\beta}_{\lambda}^{*(2)})$
$\displaystyle=\sum_{h=1}^{H}(\hat{\bm{\beta}}_{\lambda}^{(2)}-\bm{\beta}_{\lambda}^{*(2)})^{{\rm
T}}\Bigg{\\{}\sum_{i=1}^{N}\frac{\partial}{\partial\bm{\beta}^{(2)}}\frac{t_{i}^{\dagger[h]}}{e^{[h]}(\bm{z}_{i}^{\dagger})}\log
f(y_{i}^{\dagger[h]}\mid\bm{x}_{i}^{\dagger[h]};\bm{\beta}_{\lambda}^{*})-N\lambda\hskip
1.42262pt{\rm sgn}(\bm{\beta}_{\lambda}^{*(2)})\Bigg{\\}}$ $\displaystyle\
\phantom{=}+\sum_{h=1}^{H}(\hat{\bm{\beta}}_{\lambda}^{(2)}-\bm{\beta}_{\lambda}^{*(2)})^{{\rm
T}}\Bigg{\\{}\sum_{i=1}^{N}\frac{\partial^{2}}{\partial\bm{\beta}^{(2)}\partial\bm{\beta}^{(2){\rm
T}}}\frac{t_{i}^{\dagger[h]}}{e^{[h]}(\bm{z}_{i}^{\dagger})}\log
f(y_{i}^{\dagger[h]}\mid\bm{x}_{i}^{\dagger[h]};\bm{\beta}_{\lambda}^{*})\Bigg{\\}}(\hat{\bm{\beta}}_{\lambda}^{(2)}-\bm{\beta}_{\lambda}^{*(2)})$

for the third term, the aggregate of the two removes the ${\rm o}_{\rm P}(1)$
terms. The aggregate is asymptotically equivalent to

$\displaystyle(\hat{\bm{\beta}}_{\lambda}^{(2)}-\bm{\beta}_{\lambda}^{*(2)})^{{\rm
T}}\Bigg{\\{}\sum_{i=1}^{N}\sum_{h=1}^{H}\frac{\partial}{\partial\bm{\beta}^{(2)}}\frac{t_{i}^{[h]}}{e^{[h]}(\bm{z}_{i})}\log
f(y_{i}^{[h]}\mid\bm{x}_{i}^{[h]};\bm{\beta}_{\lambda}^{*})-N\lambda\hskip
1.42262pt{\rm sgn}(\bm{\beta}_{\lambda}^{*(2)})\Bigg{\\}},$

and if $\hat{\bm{\beta}}_{\lambda}^{(2)}-\bm{\beta}_{\lambda}^{*(2)}$ is
replaced by $\sqrt{N}$ times (13), we obtain

$\displaystyle b^{\rm limit}=\ $
$\displaystyle\frac{1}{N}\Bigg{\\{}\sum_{i=1}^{N}\sum_{h=1}^{H}\frac{\partial}{\partial\bm{\beta}^{(2)}}\frac{t_{i}^{[h]}}{e^{[h]}(\bm{z}_{i})}\log
f(y_{i}^{[h]}\mid\bm{x}_{i}^{[h]};\bm{\beta}_{\lambda}^{*})-N\lambda\hskip
1.42262pt{\rm sgn}(\bm{\beta}_{\lambda}^{*(2)})\Bigg{\\}}$
$\displaystyle\bm{J}^{(22)}(\bm{\beta}_{\lambda}^{*})^{-1}\Bigg{\\{}\sum_{i^{\prime}=1}^{N}\sum_{h^{\prime}=1}^{H}\frac{\partial}{\partial\bm{\beta}^{(2)}}\frac{t_{i^{\prime}}^{[h^{\prime}]}}{e^{[h^{\prime}]}(\bm{z}_{i^{\prime}})}\log
f(y_{i^{\prime}}^{[h^{\prime}]}\mid\bm{x}_{i^{\prime}}^{[h^{\prime}]};\bm{\beta}_{\lambda}^{*})-N\lambda\hskip
1.42262pt{\rm sgn}(\bm{\beta}_{\lambda}^{*(2)})\Bigg{\\}}.$

From independence, we can divide $i\neq i^{\prime}$, ${\rm E}(b^{\rm limit})$
into the expectation about $i$ and the expectation about $i^{\prime}$ and see
that the expectations about both are $\bm{0}_{p}$ from (11). Therefore, it
follows that

$\displaystyle{\rm E}(b^{\rm limit})={\rm E}\Bigg{[}$
$\displaystyle\bm{J}^{(22)}(\bm{\beta}_{\lambda}^{*})^{-1}\frac{1}{N}\sum_{i=1}^{N}\Bigg{\\{}\sum_{h=1}^{H}\frac{\partial}{\partial\bm{\beta}^{(2)}}\frac{t_{i}^{[h]}}{e^{[h]}(\bm{z}_{i})}\log
f(y_{i}^{[h]}\mid\bm{x}_{i}^{[h]};\bm{\beta}_{\lambda}^{*})-\lambda\hskip
1.42262pt{\rm sgn}(\bm{\beta}_{\lambda}^{*(2)})\Bigg{\\}}$
$\displaystyle\Bigg{\\{}\sum_{h^{\prime}=1}^{H}\frac{\partial}{\partial\bm{\beta}^{(2)}}\frac{t_{i}^{[h^{\prime}]}}{e^{[h^{\prime}]}(\bm{z}_{i})}\log
f(y_{i}^{[h^{\prime}]}\mid\bm{x}_{i}^{[h^{\prime}]};\bm{\beta}_{\lambda}^{*})-\lambda\hskip
1.42262pt{\rm sgn}(\bm{\beta}_{\lambda}^{*(2)})\Bigg{\\}}\Bigg{]}.$ (14)

When $h\neq h^{\prime}$, this expectation is $0$ from
$t_{i}^{[h]}t_{i}^{[h^{\prime}]}=0$, and when $h=h^{\prime}$, that term
remains in the expectation since $t_{i}^{[h]2}=t_{i}^{[h]}$. If we write

$\displaystyle\bm{Q}_{\lambda}(\bm{\beta})\equiv\ $ $\displaystyle{\rm
E}\Bigg{\\{}\sum_{h=1}^{H}\frac{1}{e^{[h]}(\bm{z})}\frac{\partial}{\partial\bm{\beta}}\log
f(y^{[h]}\mid\bm{x}^{[h]};\bm{\beta})\frac{\partial}{\partial\bm{\beta}^{{\rm
T}}}\log f(y^{[h]}\mid\bm{x}^{[h]};\bm{\beta})\Bigg{\\}}$
$\displaystyle-\lambda{\rm E}\Bigg{\\{}\sum_{h=1}^{H}\hskip 1.42262pt{\rm
sgn}(\bm{\beta})\frac{\partial}{\partial\bm{\beta}^{{\rm T}}}\log
f(y^{[h]}\mid\bm{x}^{[h]};\bm{\beta})\Bigg{\\}},$ (15)

and further set
$\bm{Q}_{\lambda}^{(22)}\equiv(Q_{jj^{\prime}})_{j,j^{\prime}\in\mathcal{J}_{\lambda}^{*(2)}}$,
we obtain the following theorem.

######  Theorem 2.

Defining $\bm{J}(\bm{\beta})$ and $\bm{Q}_{\lambda}(\bm{\beta})$ by (6) and
(15), respectively, the asymptotic bias of the information criterion is given
by

$\displaystyle{\rm E}(b^{\rm limit})={\rm
tr}\\{\bm{J}^{(22)}(\bm{\beta}_{\lambda}^{*})^{-1}\bm{Q}_{\lambda}^{(22)}(\bm{\beta}_{\lambda}^{*})\\}$
(16)

for the inverse-probability-weighted LASSO estimation.

From this, we propose

$\displaystyle{\rm
IPIC}\equiv-2\sum_{i=1}^{N}\sum_{h=1}^{H}\frac{t_{i}^{[h]}}{e^{[h]}(\bm{z}_{i})}\log
f(y_{i}^{[h]}\mid\bm{x}_{i}^{[h]};\hat{\bm{\beta}}_{\lambda})+2{\rm
tr}\\{\hat{\bm{J}}^{(22)}(\hat{\bm{\beta}}_{\lambda})^{-1}\hat{\bm{Q}}_{\lambda}^{(22)}(\hat{\bm{\beta}}_{\lambda})\\}$

as the AIC for the inverse-probability-weighted LASSO estimation. Here,

$\displaystyle\hat{\bm{J}}(\bm{\beta})\equiv-\frac{1}{N}\sum_{i=1}^{N}\sum_{h=1}^{H}\frac{t_{i}^{[h]}}{e^{[h]}(\bm{z}_{i})}\frac{\partial^{2}}{\partial\bm{\beta}\partial\bm{\beta}^{{\rm
T}}}\log f(y_{i}^{[h]}\mid\bm{x}_{i}^{[h]};\bm{\beta})$ (17)

and

$\displaystyle\hat{\bm{Q}}_{\lambda}(\bm{\beta})\equiv\ $
$\displaystyle\frac{1}{N}\sum_{i=1}^{N}\sum_{h=1}^{H}\frac{t_{i}^{[h]}}{e^{[h]}(\bm{z}_{i})^{2}}\frac{\partial}{\partial\bm{\beta}}\log
f(y_{i}^{[h]}\mid\bm{x}_{i}^{[h]};\bm{\beta})\frac{\partial}{\partial\bm{\beta}^{{\rm
T}}}\log f(y_{i}^{[h]}\mid\bm{x}_{i}^{[h]};\bm{\beta})$
$\displaystyle-\frac{\lambda}{N}\sum_{i=1}^{N}\sum_{h=1}^{H}\hskip
1.42262pt{\rm
sgn}(\bm{\beta})\frac{t_{i}^{[h]}}{e^{[h]}(\bm{z}_{i})}\frac{\partial}{\partial\bm{\beta}^{{\rm
T}}}\log f(y_{i}^{[h]}\mid\bm{x}_{i}^{[h]};\bm{\beta}),$

and we set
$\hat{\bm{J}}^{(22)}\equiv(\hat{J}_{jj^{\prime}})_{j,j^{\prime}\in\hat{\mathcal{J}}_{\lambda}^{(2)}}$
and
$\hat{\bm{Q}}_{\lambda}^{(22)}\equiv(\hat{Q}_{\lambda,jj^{\prime}})_{j,j^{\prime}\in\hat{\mathcal{J}}_{\lambda}^{(2)}}$,
where $\hat{\mathcal{J}}_{\lambda}^{(2)}=\\{j:\hat{\beta}_{\lambda,j}\neq
0\\}$.

###  4.2 AIC for doubly robust LASSO

Next, let us consider the case where the propensity score is unknown.
Specifically, in its model $e^{[h]}(\bm{z}_{i};\bm{\alpha})\equiv{\rm
P}(t_{i}^{[h]}=1\mid\bm{z}_{i};\bm{\alpha})$, the parameter $\bm{\alpha}\
(\in\mathbb{R}^{q})$ is supposed to be unknown. In addition, supposing that
this modeling may be wrong, we use the doubly robust estimation. We write the
model of the conditional distribution of $y^{[h]}$ given $\bm{z}$ as
$f^{[h]}(y^{[h]}\mid\bm{z};\bm{\gamma})$ with parameter $\bm{\gamma}\
(\in\mathbb{R}^{r})$. To simplify the notation in the following, we will
replace $e^{[h]}(\bm{z}_{i};\bm{\alpha})$, $\log
f(y_{i}^{[h]}\mid\bm{x}_{i}^{[h]};\bm{\beta})$,
$f^{[h]}(y_{i}^{[h]}\mid\bm{z}_{i};\bm{\gamma})$ and $\int
f_{i}^{[h]}(\bm{\gamma})\ell_{i}^{[h]}(\bm{\beta}){\rm d}y_{i}^{[h]}$ by
$e_{i}^{[h]}(\bm{\alpha})$, $\ell_{i}^{[h]}(\bm{\beta})$,
$f_{i}^{[h]}(\bm{\gamma})$ and ${\rm
E}_{i}^{[h]}\\{\ell_{i}^{[h]}(\bm{\beta})\mid\bm{\gamma}\\}$, respectively.
Accordingly, the estimator given by

$\displaystyle\hat{\bm{\beta}}_{\lambda}\equiv\mathop{\rm
argmin}_{\bm{\beta}}\Bigg{(}-\frac{1}{N}\sum_{i=1}^{N}\sum_{h=1}^{H}\Bigg{[}\frac{t_{i}^{[h]}}{e_{i}^{[h]}(\hat{\bm{\alpha}})}\ell_{i}^{[h]}(\bm{\beta})+\Bigg{\\{}1-\frac{t_{i}^{[h]}}{e_{i}^{[h]}(\hat{\bm{\alpha}})}\Bigg{\\}}{\rm
E}_{i}^{[h]}\\{\ell_{i}^{[h]}(\bm{\beta})\mid\hat{\bm{\gamma}}\\}\Bigg{]}+\lambda\|\bm{\beta}\|_{1}\Bigg{)}$

is regarded as doubly robust. That is, if either $e_{i}^{[h]}(\bm{\alpha})$ or
$f_{i}^{[h]}(\bm{\gamma})$ is the correct model, since this objective function
converges uniformly, $\hat{\bm{\beta}}_{\lambda}$ converges in probability to
$\bm{\beta}_{\lambda}^{*}$.

Similar to (12) and (13) in the previous subsection, we have

$\displaystyle N\hat{\bm{\beta}}_{\lambda}^{(1)}={\rm o}_{\rm P}(1)$ (18)

and

$\displaystyle\sqrt{N}(\hat{\bm{\beta}}_{\lambda}^{(2)}-\bm{\beta}_{\lambda}^{*(2)})=\
$
$\displaystyle\bm{J}^{(22)}(\bm{\beta}_{\lambda}^{*})^{-1}\Bigg{(}\frac{1}{\sqrt{N}}\sum_{i=1}^{N}\sum_{h=1}^{H}\frac{\partial}{\partial\bm{\beta}^{(2)}}\Bigg{[}\frac{t_{i}^{[h]}}{e_{i}^{[h]}(\bm{\alpha}^{*})}\ell_{i}^{[h]}(\bm{\beta}_{\lambda}^{*})$
$\displaystyle+\Bigg{\\{}1-\frac{t_{i}^{[h]}}{e_{i}^{[h]}(\bm{\alpha}^{*})}\Bigg{\\}}{\rm
E}_{i}^{[h]}\\{\ell_{i}^{[h]}(\bm{\beta}_{\lambda}^{*})\mid\bm{\gamma}^{*}\\}\Bigg{]}-\sqrt{N}\lambda\hskip
1.42262pt{\rm sgn}(\bm{\beta}_{\lambda}^{*(2)})\Bigg{)}+{\rm o}_{\rm P}(1)$
(19)

by following the same argument as in Knight and Fu (2000). See the Appendix
for their derivations. Accordingly, the asymptotic bias in the risk can be
evaluated in the same way as in (14). Indeed, if we write

$\displaystyle\bm{R}_{\lambda}(\bm{\beta};\bm{\alpha},\bm{\gamma})\equiv{\rm
E}\Bigg{(}$
$\displaystyle\Bigg{\\{}\sum_{h=1}^{H}\frac{\partial}{\partial\bm{\beta}}\frac{t^{[h]}}{e^{[h]}(\bm{\alpha})}\ell^{[h]}(\bm{\beta})-\lambda\hskip
1.42262pt{\rm sgn}(\bm{\beta})\Bigg{\\}}$
$\displaystyle\sum_{h=1}^{H}\frac{\partial}{\partial\bm{\beta}^{{\rm
T}}}\Bigg{[}\frac{t^{[h]}}{e^{[h]}(\bm{\alpha})}\ell^{[h]}(\bm{\beta})+\Bigg{\\{}1-\frac{t^{[h]}}{e^{[h]}(\bm{\alpha})}\Bigg{\\}}{\rm
E}^{[h]}\\{\ell^{[h]}(\bm{\beta})\mid\bm{\gamma}\\}\Bigg{]}\Bigg{)},$ (20)

and set
$\bm{R}_{\lambda}^{(22)}\equiv(R_{\lambda,jj^{\prime}})_{j,j^{\prime}\in\mathcal{J}_{\lambda}^{*(2)}}$,
we obtain the following lemma.

######  Lamma 1.

Defining $\bm{J}(\bm{\beta})$ and
$\bm{R}_{\lambda}(\bm{\beta};\bm{\alpha},\bm{\gamma})$ by (6) and (20),
respectively, the asymptotic bias of the information criterion is given by

$\displaystyle{\rm E}(b^{\rm limit})={\rm
tr}\\{\bm{J}^{(22)}(\bm{\beta}_{\lambda}^{*})^{-1}\bm{R}_{\lambda}^{(22)}(\bm{\beta}_{\lambda}^{*};\bm{\alpha}^{*},\bm{\gamma}^{*})\\}$

for the doubly robust LASSO estimation.

From this, we propose

$\displaystyle-2\sum_{i=1}^{N}\sum_{h=1}^{H}\frac{t_{i}^{[h]}}{e_{i}^{[h]}(\hat{\bm{\alpha}})}\log
f(y_{i}^{[h]}\mid\bm{x}_{i}^{[h]};\hat{\bm{\beta}}_{\lambda})+2{\rm
tr}\\{\hat{\bm{J}}^{(22)}(\hat{\bm{\beta}}_{\lambda})^{-1}\hat{\bm{R}}_{\lambda}^{(22)}(\hat{\bm{\beta}}_{\lambda};\hat{\bm{\alpha}},\hat{\bm{\gamma}})\\}$

as the AIC for the doubly robust LASSO estimation, where
$\hat{\bm{J}}(\bm{\beta})$ is the matrix defined by (17),

$\displaystyle\hat{\bm{R}}_{\lambda}(\bm{\beta};\bm{\alpha},\bm{\gamma})\equiv\frac{1}{N}\sum_{i=1}^{N}$
$\displaystyle\Bigg{\\{}\sum_{h=1}^{H}\frac{t_{i}^{[h]}}{e_{i}^{[h]}(\bm{\alpha})}\frac{\partial}{\partial\bm{\beta}}\ell_{i}^{[h]}(\bm{\beta})-\lambda\hskip
1.42262pt{\rm sgn}(\bm{\beta})\Bigg{\\}}$
$\displaystyle\Bigg{(}\sum_{h=1}^{H}\frac{\partial}{\partial\bm{\beta}^{{\rm
T}}}\Bigg{[}\frac{t_{i}^{[h]}}{e_{i}^{[h]}(\bm{\alpha})}\ell_{i}^{[h]}(\bm{\beta})+\Bigg{\\{}1-\frac{t_{i}^{[h]}}{e_{i}^{[h]}(\bm{\alpha})}\Bigg{\\}}{\rm
E}_{i}^{[h]}\\{\ell_{i}^{[h]}(\bm{\beta})\mid\bm{\gamma}\\}\Bigg{]}\Bigg{)},$
(21)

and
$\hat{\bm{R}}_{\lambda}^{(22)}=(\hat{R}_{\lambda,jj^{\prime}})_{j,j^{\prime}\in\hat{\mathcal{J}}_{\lambda}^{(2)}}$.

###  4.3 DRIC for doubly robust LASSO

The doubly robust estimator is consistent even if the model for the assignment
variable or the model for the outcome variable is misspecified. Although the
previous subsection proposed an information criterion for the doubly robust
estimation, it was derived under the assumption that the two models are both
correct and thus is not valid when only one of them is correct. In this
subsection, we derive a so-called doubly robust information criterion with
validity when only one of them is correct, based on the idea in Baba and
Ninomiya (2021). Note that in the setting of doubly robust estimation, the
propensity score is unknown, i.e. $\bm{\alpha}$ is unknown.

First, in place of
$\sqrt{N}(\hat{\bm{\beta}}_{\lambda}^{(2)}-\bm{\beta}_{\lambda}^{*(2)})$, we
only need to use

$\displaystyle{\rm
E}\Bigg{(}-\sum_{h=1}^{H}\frac{\partial^{2}}{\partial\bm{\beta}^{(2)}\partial\bm{\beta}^{(2){\rm
T}}}\Bigg{[}\frac{t^{[h]}}{e^{[h]}(\bm{\alpha}^{*})}\ell^{[h]}(\bm{\beta}_{\lambda}^{*})+\Bigg{\\{}1-\frac{t^{[h]}}{e^{[h]}(\bm{\alpha}^{*})}\Bigg{\\}}{\rm
E}^{[h]}\\{\ell^{[h]}(\bm{\beta}_{\lambda}^{*})\mid\bm{\gamma}^{*}\\}\Bigg{]}\Bigg{)}^{-1}$
$\displaystyle\Bigg{\\{}\frac{1}{\sqrt{N}}\sum_{i=1}^{N}\sum_{h=1}^{H}\frac{\partial}{\partial\bm{\beta}^{(2)}}\Bigg{[}\frac{t_{i}^{[h]}}{e_{i}^{[h]}(\bm{\alpha}^{*})}\ell_{i}^{[h]}(\bm{\beta}_{\lambda}^{*})+\Bigg{\\{}1-\frac{t_{i}^{[h]}}{e_{i}^{[h]}(\bm{\alpha}^{*})}\Bigg{\\}}{\rm
E}_{i}^{[h]}\\{\ell_{i}^{[h]}(\bm{\beta}_{\lambda}^{*})\mid\bm{\gamma}^{*}\\}\Bigg{]}-\sqrt{N}\lambda\hskip
1.42262pt{\rm sgn}(\bm{\beta}_{\lambda}^{*(2)})$
$\displaystyle\phantom{\Bigg{\\{}}+{\rm
E}\Bigg{(}\sum_{h=1}^{H}\frac{\partial^{2}}{\partial\bm{\beta}^{(2)}\partial\bm{\alpha}^{{\rm
T}}}\frac{t^{[h]}}{e^{[h]}(\bm{\alpha}^{*})}[\ell^{[h]}(\bm{\beta}_{\lambda}^{*})-{\rm
E}^{[h]}\\{\ell^{[h]}(\bm{\beta}_{\lambda}^{*})\mid\bm{\gamma}^{*}\\}]\Bigg{)}\sqrt{N}(\hat{\bm{\alpha}}-\bm{\alpha}^{*})$
$\displaystyle\phantom{\Bigg{\\{}}+{\rm
E}\Bigg{[}\sum_{h=1}^{H}\Bigg{\\{}1-\frac{t^{[h]}}{e^{[h]}(\bm{\alpha}^{*})}\Bigg{\\}}\frac{\partial^{2}}{\partial\bm{\beta}^{(2)}\partial\bm{\gamma}^{{\rm
T}}}{\rm
E}^{[h]}\\{\ell^{[h]}(\bm{\beta}_{\lambda}^{*})\mid\bm{\gamma}^{*}\\}\Bigg{]}\sqrt{N}(\hat{\bm{\gamma}}-\bm{\gamma}^{*})\Bigg{\\}}.$
(22)

See the Appendix for the derivation. Now, using the fact that

$\displaystyle\sqrt{N}(\hat{\bm{\alpha}}-\bm{\alpha}^{*})={\rm
E}\Bigg{\\{}-\sum_{h=1}^{H}\frac{\partial^{2}}{\partial\bm{\alpha}\partial\bm{\alpha}^{{\rm
T}}}t^{[h]}\log
e^{[h]}(\bm{\alpha}^{*})\Bigg{\\}}^{-1}\frac{1}{\sqrt{N}}\sum_{i=1}^{N}\sum_{h=1}^{H}\frac{\partial}{\partial\bm{\alpha}}t_{i}^{[h]}\log
e_{i}^{[h]}(\bm{\alpha}^{*})$ (23)

and

$\displaystyle\sqrt{N}(\hat{\bm{\gamma}}-\bm{\gamma}^{*})={\rm
E}\Bigg{\\{}-\sum_{h=1}^{H}\frac{\partial^{2}}{\partial\bm{\gamma}\partial\bm{\gamma}^{{\rm
T}}}t^{[h]}\log
f^{[h]}(\bm{\gamma}^{*})\Bigg{\\}}^{-1}\frac{1}{\sqrt{N}}\sum_{i=1}^{N}\sum_{h=1}^{H}\frac{\partial}{\partial\bm{\gamma}}t_{i}^{[h]}\log
f_{i}^{[h]}(\bm{\gamma}^{*})$ (24)

from the conventional asymptotics, we can rewrite (22) as

$\displaystyle\frac{1}{\sqrt{N}}\sum_{i=1}^{N}\sum_{h=1}^{H}t_{i}^{[h]}\Bigg{\\{}\bm{C}_{1}^{[h](2)}(\bm{\beta}_{\lambda}^{*},\bm{\alpha}^{*},\bm{\gamma}^{*})^{{\rm
T}}\frac{\partial}{\partial\bm{\alpha}}\log
e_{i}^{[h]}(\bm{\alpha}^{*})+\bm{C}_{2}^{[h](2)}(\bm{\beta}_{\lambda}^{*},\bm{\alpha}^{*},\bm{\gamma}^{*})^{{\rm
T}}\frac{\partial}{\partial\bm{\gamma}}\log
f_{i}^{[h]}(\bm{\gamma}^{*})\Bigg{\\}},$

where

$\displaystyle\bm{C}_{1}^{[h]}(\bm{\beta},\bm{\alpha},\bm{\gamma})\equiv\ $
$\displaystyle{\rm
E}\Bigg{\\{}-\sum_{h^{\prime}=1}^{H}\frac{\partial^{2}}{\partial\bm{\alpha}\partial\bm{\alpha}^{{\rm
T}}}t_{i}^{[h^{\prime}]}\log e_{i}^{[h^{\prime}]}(\bm{\alpha})\Bigg{\\}}^{-1}$
$\displaystyle{\rm
E}\Bigg{(}\frac{t_{i}^{[h]}}{e_{i}^{[h]}(\bm{\alpha})}\frac{\partial}{\partial\bm{\alpha}}\log
e_{i}^{[h]}(\bm{\alpha})\Bigg{[}\frac{\partial}{\partial\bm{\beta}}\ell_{i}^{[h]}(\bm{\beta})-\frac{\partial}{\partial\bm{\beta}}{\rm
E}_{i}^{[h]}\\{\ell_{i}^{[h]}(\bm{\beta})\mid\bm{\gamma}\\}\Bigg{]}^{{\rm
T}}\Bigg{)}$

and

$\displaystyle\bm{C}_{2}^{[h]}(\bm{\beta},\bm{\alpha},\bm{\gamma})\equiv\ $
$\displaystyle{\rm
E}\Bigg{\\{}-\sum_{h^{\prime}=1}^{H}\frac{\partial^{2}}{\partial\bm{\gamma}\partial\bm{\gamma}^{{\rm
T}}}t_{i}^{[h^{\prime}]}\log
f_{i}^{[h^{\prime}]}(\bm{\gamma})\Bigg{\\}}^{-1}{\rm
E}\Bigg{[}\Bigg{\\{}1-\frac{t_{i}^{[h]}}{e_{i}^{[h]}(\bm{\alpha})}\Bigg{\\}}\frac{\partial^{2}}{\partial\bm{\gamma}\partial\bm{\beta}^{{\rm
T}}}{\rm E}^{[h]}\\{\ell_{i}^{[h]}(\bm{\beta})\mid\bm{\gamma}\\}\Bigg{]}.$

If we write

$\displaystyle\bm{K}(\bm{\beta};\bm{\alpha},\bm{\gamma})\equiv{\rm
E}\Bigg{(}-\sum_{h=1}^{H}\frac{\partial^{2}}{\partial\bm{\beta}\partial\bm{\beta}^{{\rm
T}}}\Bigg{[}\frac{t^{[h]}}{e^{[h]}(\bm{\alpha})}\ell^{[h]}(\bm{\beta})+\Bigg{\\{}1-\frac{t^{[h]}}{e^{[h]}(\bm{\alpha})}\Bigg{\\}}{\rm
E}^{[h]}\\{\ell^{[h]}(\bm{\beta})\mid\bm{\gamma}\\}\Bigg{]}\Bigg{)}$ (25)

and

$\displaystyle\bm{S}_{\lambda}(\bm{\beta};\bm{\alpha},\bm{\gamma})\equiv{\rm
E}\Bigg{[}\sum_{h=1}^{H}$
$\displaystyle\Bigg{\\{}\frac{\partial}{\partial\bm{\beta}}\frac{t^{[h]}}{e^{[h]}(\bm{\alpha})}\ell^{[h]}(\bm{\beta})-\lambda\hskip
1.42262pt{\rm sgn}(\bm{\beta})\Bigg{\\}}$
$\displaystyle\Bigg{\\{}\frac{\partial}{\partial\bm{\alpha}^{{\rm T}}}\log
e^{[h]}(\bm{\alpha})\bm{C}_{1}^{[h]}(\bm{\beta},\bm{\alpha},\bm{\gamma})+\frac{\partial}{\partial\bm{\gamma}^{{\rm
T}}}\log
f^{[h]}(\bm{\gamma})\bm{C}_{2}^{[h]}(\bm{\beta},\bm{\alpha},\bm{\gamma})\Bigg{\\}}\Bigg{]},$
(26)

and let
$\bm{K}^{(22)}\equiv(K_{jj^{\prime}})_{j,j^{\prime}\in\mathcal{J}_{\lambda}^{*(2)}}$
and
$\bm{S}_{\lambda}^{(22)}\equiv(S_{\lambda,jj^{\prime}})_{j,j^{\prime}\in\mathcal{J}_{\lambda}^{*(2)}}$,
we obtain the following theorem.

######  Theorem 3.

Defining $\bm{R}_{\lambda}(\bm{\beta}_{\lambda};\bm{\alpha},\bm{\gamma})$,
$\bm{K}(\bm{\beta};\bm{\alpha},\bm{\gamma})$ and
$\bm{S}_{\lambda}(\bm{\beta};\bm{\alpha},\bm{\gamma})$ by (20), (25) and (26),
respectively, the asymptotic bias of the information criterion is given by

$\displaystyle{\rm E}(b^{\rm limit})={\rm
tr}[\bm{K}^{(22)}(\bm{\beta}_{\lambda}^{*};\bm{\alpha}^{*},\bm{\gamma}^{*})^{-1}\\{\bm{R}_{\lambda}^{(22)}(\bm{\beta}_{\lambda}^{*};\bm{\alpha}^{*},\bm{\gamma}^{*})+\bm{S}_{\lambda}^{(22)}(\bm{\beta}_{\lambda}^{*};\bm{\alpha}^{*},\bm{\gamma}^{*})\\}]$
(27)

for the doubly robust LASSO estimation even if either of
$e^{[h]}(\bm{\alpha})$ or $f^{[h]}(\bm{\gamma})$ is misspecified as a model.

From this, we propose

$\displaystyle{\rm DRIC}\equiv$
$\displaystyle-2\sum_{i=1}^{N}\sum_{h=1}^{H}\frac{t_{i}^{[h]}}{e_{i}^{[h]}(\hat{\bm{\alpha}})}\log
f(y_{i}^{[h]}\mid\bm{x}_{i}^{[h]};\hat{\bm{\beta}}_{\lambda})$
$\displaystyle+2{\rm
tr}[\hat{\bm{K}}^{(22)}(\hat{\bm{\beta}}_{\lambda};\hat{\bm{\alpha}},\hat{\bm{\gamma}})^{-1}\\{\hat{\bm{R}}_{\lambda}^{(22)}(\hat{\bm{\beta}}_{\lambda};\hat{\bm{\alpha}},\hat{\bm{\gamma}})+\hat{\bm{S}}_{\lambda}^{(22)}(\hat{\bm{\beta}}_{\lambda};\hat{\bm{\alpha}},\hat{\bm{\gamma}})\\}]$

as the doubly robust information criterion for the doubly robust LASSO
estimation, where $\hat{\bm{R}}_{\lambda}(\bm{\beta};\bm{\alpha},\bm{\gamma})$
is the matrix defined by (21),

$\displaystyle\hat{\bm{K}}(\bm{\beta};\bm{\alpha},\bm{\gamma})\equiv-\frac{1}{N}\sum_{i=1}^{N}\sum_{h=1}^{H}\frac{\partial^{2}}{\partial\bm{\beta}\partial\bm{\beta}^{{\rm
T}}}\Bigg{[}\frac{t_{i}^{[h]}}{e_{i}^{[h]}(\bm{\alpha})}\ell_{i}^{[h]}(\bm{\beta})+\Bigg{\\{}1-\frac{t_{i}^{[h]}}{e_{i}^{[h]}(\bm{\alpha})}\Bigg{\\}}{\rm
E}_{i}^{[h]}\\{\ell_{i}^{[h]}(\bm{\beta})\mid\bm{\gamma}\\}\Bigg{]}$

and

$\displaystyle\hat{\bm{S}}_{\lambda}(\bm{\beta};\bm{\alpha},\bm{\gamma})\equiv\frac{1}{N}\sum_{i=1}^{N}\sum_{h=1}^{H}$
$\displaystyle\Bigg{\\{}\frac{\partial}{\partial\bm{\beta}}\frac{t_{i}^{[h]}}{e_{i}^{[h]}(\bm{\alpha})}\ell_{i}^{[h]}(\bm{\beta})-\lambda\hskip
1.42262pt{\rm sgn}(\bm{\beta})\Bigg{\\}}$
$\displaystyle\Bigg{\\{}\frac{\partial}{\partial\bm{\alpha}^{{\rm T}}}\log
e_{i}^{[h]}(\bm{\alpha})\bm{C}_{1}^{[h]}(\bm{\beta},\bm{\alpha},\bm{\gamma})+\frac{\partial}{\partial\bm{\gamma}^{{\rm
T}}}\log
f_{i}^{[h]}(\bm{\gamma})\bm{C}_{2}^{[h]}(\bm{\beta},\bm{\alpha},\bm{\gamma})\Bigg{\\}},$

as well as
$\hat{\bm{K}}^{(22)}\equiv(\hat{K}_{jj^{\prime}})_{j,j^{\prime}\in\hat{\mathcal{J}}_{\lambda}^{(2)}}$
and
$\hat{\bm{S}}^{(22)}_{\lambda}\equiv(\hat{S}_{\lambda,jj^{\prime}})_{j,j^{\prime}\in\hat{\mathcal{J}}_{\lambda}^{(2)}}$.

##  5 Numerical experiment

###  5.1 Case of continuous outcome

For the model given by (2), suppose there is a control group represented by
$h=1$ and a treatment group represented by $h=2$ ($H=2$). In addition, let us
assume that the explanatory variable $\bm{x}=(x_{1},x_{2},\ldots,x_{p})^{{\rm
T}}$ follows ${\rm N}(\bm{0}_{p},\bm{I}_{p})$ and that the confounding
variable $z$, which affects both $y^{[h]}$ and $t^{[h]}$, independently
follows ${\rm N}(0,1)$. For $t^{[h]}$ given $z$, we use a logit model with $z$
as the covariate, defined by

$\displaystyle{\rm E}(t^{[1]}\mid z)=\frac{\exp(z)}{1+\exp(z)},\qquad
t^{[2]}=1-t^{[1]}.$

On the other hand, for $y^{[h]}$ given $\bm{x}$ and $z$, we use a linear
regression model with them as covariates, defined by

$\displaystyle y^{[h]}\mid(\bm{x},z)\sim{\rm N}\\{\bm{x}^{{\rm
T}}\bm{\theta}^{[h]}+\eta(z),\sigma^{2}\\}.$

The true value used to generate the data is
$\bm{\theta}^{[1]}=(\theta_{1}^{*}\bm{1}_{p/4}^{{\rm
T}},\allowbreak\theta_{2}^{*}\bm{1}_{p/4}^{{\rm T}},\bm{0}_{p/2}^{{\rm T}})$,
$\bm{\theta}^{[2]}=-(\theta_{1}^{*}\bm{1}_{p/4}^{{\rm
T}},\allowbreak\theta_{2}^{*}\bm{1}_{p/4}^{{\rm T}},\bm{0}_{p/2}^{{\rm T}})$,
$\eta(z)=z$ and $\sigma^{2}=1$, where $\bm{1}_{p}$ is a $p$-dimensional $1$
vector. Let $0$, $0.2$ or $0.4$ as $\theta_{1}^{*}$ and $\theta_{2}^{*}$, let
$8$, $16$ or $32$ be the dimension of the explanatory variable $p$, and let
$40$, $80$ or $120$ be the sample size $N$. In addition, the causal effect of
the target is supposed to be $(\bm{\theta}^{[2]}-\bm{\theta}^{[1]})/\sqrt{2}$,
that is, $(c^{[1]},c^{[2]})=(-1/\sqrt{2},1/\sqrt{2})$.

Table 1 shows whether the penalty term of the information criterion IPCp given
by Theorem 1 approximates the third term, i.e., the bias, in the expansion of
its original risk in (8). The number of Monte Carlo iterations is set to
1,000. In all cases, we see that the bias can be evaluated with very high
accuracy. Since QICw evaluates these values at $2\times(p/8)$, i.e., the first
four rows are all evaluated at 2, the middle four rows are all evaluated at 4,
and the last four rows are all evaluated at 8, it is clear that the bias is
underestimated. The degree of underestimation is low for $H=2$, which is the
setting in this table, and becomes stronger as $H$ increases.

Table 1: Bias evaluation in causal inference models with Gaussian noise. The $j\ (\in\\{1,2,\ldots,8\\})$ column shows the bias when the number of elements in the active set is $\hat{p}=j$ minus the bias when $\hat{p}=j-1$. On the other hand, the True row gives the true value evaluated by the Monte Carlo method, and the IPCp row gives the average of estimators for the true value. $(p,\ N,\ \theta_{1}^{*},\ \theta_{2}^{*})$ |  | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8
---|---|---|---|---|---|---|---|---|---
(8, 40, 0.2, 0.2) | True | 3.36 | 3.94 | 3.90 | 4.18 | 4.42 | 5.32 | 3.90 | 5.71
IPCp | 2.98 | 4.02 | 3.91 | 4.25 | 4.43 | 5.22 | 4.30 | 5.50
(8, 40, 0.4, 0.0) | True | 3.33 | 4.03 | 3.84 | 4.32 | 4.36 | 5.14 | 4.02 | 5.69
IPCp | 2.97 | 4.03 | 3.88 | 4.13 | 4.53 | 5.18 | 4.26 | 5.63
(8, 120, 0.2, 0.2) | True | 3.44 | 3.68 | 3.69 | 4.11 | 4.97 | 5.01 | 4.51 | 5.97
IPCp | 2.90 | 3.95 | 3.88 | 4.07 | 4.92 | 4.86 | 4.68 | 5.41
(8, 120, 0.4, 0.0) | True | 3.39 | 3.82 | 4.24 | 4.08 | 4.99 | 4.65 | 4.04 | 6.39
IPCp | 2.98 | 4.08 | 3.78 | 4.39 | 4.59 | 4.72 | 4.14 | 5.90
(16, 40, 0.2, 0.2) | True | 5.63 | 5.70 | 6.89 | 7.99 | 8.12 | 10.44 | 13.29 | 10.27
IPCp | 5.57 | 6.72 | 6.36 | 8.21 | 7.74 | 10.99 | 13.04 | 10.67
(16, 40, 0.4, 0.0) | True | 5.59 | 5.80 | 6.51 | 7.98 | 7.92 | 11.14 | 12.22 | 10.90
IPCp | 5.80 | 6.19 | 6.54 | 8.29 | 7.94 | 11.25 | 11.93 | 11.32
(16, 120, 0.2, 0.2) | True | 5.67 | 5.44 | 7.20 | 7.35 | 9.50 | 10.83 | 11.60 | 11.68
IPCp | 5.36 | 5.77 | 7.58 | 6.97 | 9.21 | 10.92 | 11.46 | 11.95
(16, 120, 0.4, 0.0) | True | 6.52 | 6.18 | 5.86 | 8.33 | 7.05 | 10.91 | 13.89 | 10.66
IPCp | 6.00 | 6.55 | 6.15 | 8.11 | 7.25 | 9.90 | 14.52 | 10.70
(32, 40, 0.2, 0.2) | True | 9.61 | 10.60 | 11.46 | 11.52 | 19.35 | 19.21 | 30.89 | 24.54
IPCp | 9.81 | 10.86 | 11.58 | 11.97 | 19.87 | 18.34 | 30.74 | 25.60
(32, 40, 0.4, 0.0) | True | 10.14 | 9.87 | 13.56 | 12.43 | 13.39 | 25.05 | 28.32 | 24.05
IPCp | 10.48 | 9.99 | 14.04 | 12.28 | 14.46 | 23.09 | 29.52 | 24.82
(32, 120, 0.2, 0.2) | True | 9.15 | 8.31 | 11.38 | 13.33 | 16.92 | 16.90 | 32.18 | 32.35
IPCp | 9.00 | 8.50 | 11.07 | 13.19 | 17.45 | 15.96 | 31.34 | 32.05
(32, 120, 0.4, 0.0) | True | 10.06 | 8.27 | 12.96 | 10.89 | 18.11 | 19.38 | 29.57 | 31.70
IPCp | 10.13 | 8.64 | 12.21 | 10.95 | 16.90 | 19.09 | 29.82 | 30.72

Table 2 compares the proposed criterion IPCp with QICw. Since both criteria
basically try to estimate the causal effect $\bm{\theta}$ with good accuracy,
the mean squared error of $\hat{\bm{\theta}}$ is used as the main indicator of
the comparison. In addition, since the objective of sparse estimation is to
list the explanatory variables related to the outcome variable, it is also
important to narrow down the number of explanatory variables, with a small
mean square error as a prerequisite. An indicator of this is the number of
elements in the active set, $\hat{p}=|\hat{\mathcal{J}}_{\lambda}^{(2)}|$. As
reference indicators, the mean squared error of the $p/2$ nonzero components
in $\bm{\theta}$ and the number of elements in the active set and the mean
squared error of the $p/2$ zero components in $\bm{\theta}$ and the number of
elements in the active set are also included. In this table, the true values
of the regression parameters are set to
$\theta_{1}^{*}=\theta_{2}^{*}=\theta^{*}$, and the number of Monte Carlo
iterations is 200. It can be seen from the main indicators that IPCp
outperforms QICw in almost all settings, especially when there is the sample
size is moderate in terms of dimension. It can also be seen that IPCp narrows
down the explanatory variables considerably, even if it does not differ much
in the value of the mean squared error. Note that the setting dealt with in
this subsection is one in which $H=2$ and the difference between IPCp and QICw
is the smallest; however, the results are considerably different from each
other. The next subsection discusses how the superiority of the proposed
criterion becomes clearer as $H$ increases.

Table 2: Performance comparison of QICw and IPCp in causal inference models with Gaussian noise. The $\hat{p}_{1}$ and $\sqrt{{\rm MSE}_{1}}$ columns show the dimension of selection and the square root of the mean squared error ($\times 10$) for nonzero parameters, while the $\hat{p}_{2}$ and $\sqrt{{\rm MSE}_{2}}$ columns show those for zero parameters and the $\hat{p}$ and $\sqrt{{\rm MSE}}$ columns show those for all parameters. $(p,\ N,\ \theta^{*})$ |  | $\hat{p}_{1}$ | $\sqrt{{\rm MSE}_{1}}$ | $\hat{p}_{2}$ | $\sqrt{{\rm MSE}_{2}}$ | $\hat{p}$ | $\sqrt{{\rm MSE}}$
---|---|---|---|---|---|---|---
(8, 40, 0.2) | QICw | 3.1 [1.1] | 1.30 [0.56] | 3.0 [1.1] | 1.24 [0.61] | 6.1 [1.9] | 1.86 [0.67]
IPCp | 1.1 [1.3] | 0.97 [0.38] | 1.0 [1.3] | 0.52 [0.68] | 2.2 [2.5] | 1.21 [0.60]
(8, 80, 0.2) | QICw | 3.1 [1.1] | 0.89 [0.34] | 2.8 [1.2] | 0.83 [0.43] | 5.9 [2.0] | 1.26 [0.43]
IPCp | 1.0 [1.3] | 0.82 [0.19] | 0.7 [1.1] | 0.28 [0.41] | 1.7 [2.2] | 0.93 [0.29]
(8, 120, 0.2) | QICw | 3.1 [1.1] | 0.76 [0.26] | 2.7 [1.2] | 0.66 [0.34] | 5.8 [1.9] | 1.04 [0.32]
IPCp | 1.3 [1.4] | 0.74 [0.16] | 0.7 [1.1] | 0.23 [0.35] | 2.0 [2.2] | 0.84 [0.22]
(8, 40, 0.4) | QICw | 3.4 [0.8] | 1.39 [0.54] | 3.0 [1.1] | 1.29 [0.61] | 6.5 [1.6] | 1.97 [0.62]
IPCp | 1.9 [1.5] | 1.38 [0.37] | 1.3 [1.3] | 0.63 [0.69] | 3.2 [2.5] | 1.64 [0.47]
(8, 80, 0.4) | QICw | 3.7 [0.6] | 0.97 [0.37] | 3.0 [1.0] | 0.89 [0.42] | 6.7 [1.4] | 1.37 [0.42]
IPCp | 2.4 [1.5] | 1.13 [0.37] | 1.2 [1.2] | 0.44 [0.46] | 3.6 [2.5] | 1.30 [0.34]
(8, 120, 0.4) | QICw | 3.8 [0.5] | 0.77 [0.32] | 3.0 [1.1] | 0.73 [0.33] | 6.8 [1.3] | 1.11 [0.32]
IPCp | 2.9 [1.3] | 0.95 [0.38] | 1.3 [1.2] | 0.40 [0.39] | 4.2 [2.2] | 1.11 [0.35]
(16, 40, 0.2) | QICw | 7.0 [1.3] | 2.27 [0.84] | 7.0 [1.3] | 2.28 [0.95] | 14.0 [2.2] | 3.27 [1.11]
IPCp | 4.3 [2.8] | 1.71 [0.79] | 3.9 [2.7] | 1.41 [1.15] | 8.2 [5.3] | 2.32 [1.21]
(16, 80, 0.2) | QICw | 6.6 [1.4] | 1.36 [0.40] | 6.2 [1.6] | 1.28 [0.46] | 12.8 [2.6] | 1.90 [0.49]
IPCp | 2.5 [2.4] | 1.16 [0.27] | 1.9 [2.3] | 0.47 [0.56] | 4.4 [4.5] | 1.33 [0.42]
(16, 120, 0.2) | QICw | 6.8 [1.2] | 1.12 [0.28] | 6.2 [1.5] | 1.01 [0.36] | 13.0 [2.4] | 1.53 [0.35]
IPCp | 2.5 [2.5] | 1.07 [0.16] | 1.6 [2.1] | 0.34 [0.43] | 4.0 [4.3] | 1.19 [0.23]
(16, 40, 0.4) | QICw | 7.4 [1.0] | 2.38 [0.84] | 7.0 [1.2] | 2.37 [0.92] | 14.4 [1.8] | 3.41 [1.09]
IPCp | 5.3 [2.5] | 2.18 [0.69] | 4.4 [2.6] | 1.62 [1.16] | 9.8 [4.7] | 2.84 [1.06]
(16, 80, 0.4) | QICw | 7.5 [0.7] | 1.49 [0.41] | 6.5 [1.5] | 1.39 [0.46] | 14.1 [1.9] | 2.08 [0.47]
IPCp | 5.2 [2.5] | 1.65 [0.39] | 3.1 [2.4] | 0.75 [0.61] | 8.2 [4.5] | 1.91 [0.39]
(16, 120, 0.4) | QICw | 7.7 [0.7] | 1.22 [0.31] | 6.5 [1.5] | 1.08 [0.38] | 14.2 [1.8] | 1.66 [0.37]
IPCp | 6.0 [2.2] | 1.42 [0.40] | 3.0 [2.2] | 0.59 [0.42] | 9.0 [3.9] | 1.61 [0.33]
(32, 40, 0.2) | QICw | 15.6 [0.8] | 6.00 [2.32] | 15.6 [0.7] | 6.14 [2.27] | 31.2 [1.2] | 8.65 [3.07]
IPCp | 14.8 [1.9] | 5.65 [2.38] | 14.9 [1.9] | 5.74 [2.39] | 29.7 [3.6] | 8.12 [3.21]
(32, 80, 0.2) | QICw | 14.6 [1.6] | 2.32 [0.62] | 14.3 [1.7] | 2.26 [0.64] | 28.9 [2.9] | 3.26 [0.79]
IPCp | 9.5 [4.2] | 1.81 [0.48] | 8.3 [4.4] | 1.37 [0.75] | 17.8 [8.3] | 2.33 [0.72]
(32, 120, 0.2) | QICw | 14.4 [1.7] | 1.74 [0.37] | 13.7 [2.0] | 1.65 [0.45] | 28.1 [3.2] | 2.42 [0.48]
IPCp | 8.4 [4.5] | 1.49 [0.22] | 6.3 [4.4] | 0.85 [0.57] | 14.8 [8.5] | 1.78 [0.38]
(32, 40, 0.4) | QICw | 15.5 [0.8] | 6.16 [2.36] | 15.6 [0.6] | 6.29 [2.33] | 31.1 [1.2] | 8.87 [3.14]
IPCp | 15.0 [1.5] | 5.82 [2.35] | 14.9 [1.9] | 5.90 [2.40] | 29.9 [3.1] | 8.35 [3.20]
(32, 80, 0.4) | QICw | 15.2 [1.0] | 2.50 [0.61] | 14.4 [1.6] | 2.39 [0.66] | 29.6 [2.1] | 3.49 [0.78]
IPCp | 12.7 [3.1] | 2.35 [0.47] | 9.3 [4.1] | 1.66 [0.77] | 22.0 [6.6] | 2.95 [0.61]
(32, 120, 0.4) | QICw | 15.5 [0.8] | 1.90 [0.39] | 13.9 [2.1] | 1.76 [0.48] | 29.4 [2.5] | 2.61 [0.49]
IPCp | 13.5 [2.5] | 1.97 [0.38] | 8.8 [3.8] | 1.20 [0.54] | 22.3 [5.7] | 2.36 [0.41]

It should be emphasized that IPCp was derived without resorting to
asymptotics, so it should work well even in settings where the conventional
large-sample theory would not hold, for example, when $p=32$ and $N=40$. On
the other hand, it does not deal with doubly robust estimation. As shown in
Section 4, asymptotics can be used to handle the doubly robust estimation
without any difficulty, and with a certain amount of effort, it is possible to
construct DRIC with double robustness in the criterion itself. See the
Appendix for a numerical verification of this assertion. To show that the
information criterion can be derived in various settings, the setup is
slightly changed; i.e., each $\bm{\theta}^{[h]}$ is estimated instead of
$\bm{\theta}$ being estimated directly, while $\eta(z)$ is not required to be
common to $h$.

###  5.2 Case of discrete outcome

For the model in (3), we assume that the explanatory variable
$\bm{x}=(x_{1},x_{2},\allowbreak\ldots,x_{p})^{{\rm T}}$ follows ${\rm
N}(\bm{0}_{p},\bm{I}_{p})$ and that the confounding variable $z$ follows ${\rm
N}(0,1)$ independently of $\bm{x}$. Also, to deal with a setting that may
differ in aspect from the previous subsection, for $y^{[h]}\
(\in\\{0,1,\ldots,m\\})$ given $\bm{x}$, we suppose a logit model with $m$
trials defined by

$\displaystyle{\rm E}(y^{[h]}\mid\bm{x})=\frac{m\exp(\bm{x}^{{\rm
T}}\bm{\beta}^{[h]})}{1+\exp(\bm{x}^{{\rm T}}\bm{\beta}^{[h]})},$

where
$\bm{\beta}^{[h]}\equiv(\beta_{1}^{[h]},\beta_{2}^{[h]},\ldots,\beta_{p}^{[h]})^{{\rm
T}}$ is the parameter vector. Denoting the unit vector with $1$ in the $h$-th
element as $\bm{e}^{[h]}$ and letting
$\bm{x}^{[h]}\equiv\bm{x}\otimes\bm{e}^{[h]}$ and
$\bm{\beta}\equiv(\bm{\beta}^{[1]{\rm T}},\bm{\beta}^{[2]{\rm
T}},\ldots,\bm{\beta}^{[H]{\rm T}})^{{\rm T}}$, it holds that $\bm{x}^{{\rm
T}}\bm{\beta}^{[h]}=\bm{x}^{[h]{\rm T}}\bm{\beta}$, so this model can also be
expressed using (3). For $t^{[h]}$ and $y^{[h]}$ given $z$, we use the logit
models defined by

$\displaystyle{\rm E}(t^{[h]}\mid
z)=\frac{\exp(z\alpha^{[h]})}{\sum_{h^{\prime}=1}^{H}\exp(z\alpha^{[h^{\prime}]})},\qquad{\rm
E}(y^{[h]}\mid z)=\frac{m\exp(z\gamma)}{1+\exp(z\gamma)}$

with $\alpha^{[h]}$ and $\gamma$ as parameters ($h\in\\{1,2,\ldots,H\\}$). As
a natural progression from these models, letting
$\bm{z}\equiv(z,z^{2}/\sqrt{2}-1/\sqrt{2})^{{\rm T}}$,
$\bm{\alpha}\equiv(\alpha_{1}^{[h]},\alpha_{2}^{[h]})^{{\rm T}}$ and
$\bm{\gamma}\equiv(\gamma_{1},\gamma_{2})^{{\rm T}}$, we generate data from

$\displaystyle{\rm E}(t^{[h]}\mid\bm{z})=\frac{\exp(\bm{z}^{{\rm
T}}\bm{\alpha}^{[h]})}{\sum_{h^{\prime}=1}^{H}\exp(\bm{z}^{{\rm
T}}\bm{\alpha}^{[h^{\prime}]})},\qquad{\rm
E}(y^{[h]}\mid\bm{x},\bm{z})=\frac{m\exp(\bm{x}^{{\rm
T}}\bm{\beta}^{[h]}+\bm{z}^{{\rm T}}\bm{\gamma})}{1+\exp(\bm{x}^{{\rm
T}}\bm{\beta}^{[h]}+\bm{z}^{{\rm T}}\bm{\gamma})}.$

The true values used here are $\bm{\beta}^{[1]}=(2\beta^{*},0)^{{\rm T}}$,
$\bm{\beta}^{[2]}=(0,\beta^{*})^{{\rm T}}$,
$\bm{\beta}^{[3]}=(-\beta^{*},0)^{{\rm T}}$ and
$\bm{\beta}^{[4]}=(0,-2\beta^{*})^{{\rm T}}$ when $(p,H)=(2,4)$. Moreover,
$\bm{\beta}^{[5]}=(2\beta^{*},0)^{{\rm T}}$ and
$\bm{\beta}^{[6]}=(0,\beta^{*})^{{\rm T}}$ are added to these when
$(p,H)=(2,6)$, and $\bm{\beta}^{[7]}=(-\beta^{*},0)^{{\rm T}}$ and
$\bm{\beta}^{[8]}=(0,-2\beta^{*})^{{\rm T}}$ are further added when
$(p,H)=(2,8)$. The third component with the same value as the first component
is added in $\bm{\beta}^{[j]}$ when $(p,H)=(3,4)$, and the fourth component
with the same value of the second component is further added in
$\bm{\beta}^{[j]}$ when $(p,H)=(4,4)$ ($j\in\\{1,2,3,4\\}$). The true value of
$\bm{\alpha}$ is $(\alpha^{*},0)$ or $(\alpha^{*},\alpha^{*})$, and the true
value of $\bm{\gamma}$ is $(\gamma^{*},0)$ or $(\gamma^{*},\gamma^{*})$. When
the true value of $\bm{\alpha}$ is $(\alpha^{*},\alpha^{*})$, the model of the
assignment variable is misspecified, and when the true value of $\bm{\gamma}$
is $(\gamma^{*},\gamma^{*})$, the model of the outcome variable is
misspecified. In addition, $\alpha^{*}$ and $\beta^{*}$ are 0.1 or 0.2,
$\gamma^{*}$ is 0.2 or 0.4, the sample size $N$ is 200 or 400, and the number
of trials $m$ is 5 or 10.

Table 3 shows whether the asymptotic bias in (16) given by Theorem 2 or the
asymptotic bias in (27) given by Theorem 3 successfully approximates the bias,
i.e., the third term in the expansion of the risk in (8). The number of Monte
Carlo iterations is set to 200. The table shows that the bias can be evaluated
with a high degree of accuracy, except for the lowest setting. Since QICw
evaluates all these values at $pH/4=2$, the approximation accuracy of the
asymptotic bias is much better, at least in this comparison. In the lowest
setting, both the model of the assignment variable and the model of the
outcome variable are misspecified, in which case neither the doubly robust
estimation nor the bias evaluation in this paper is theoretically guaranteed
and, thus, the approximation accuracy is low.

Table 3: Bias evaluation in logit causal inference models. The upper entry of the parameters is $(p,H,N,m)$, the middle entry is $(\beta^{*},\alpha^{*},\gamma^{*})$, and the lower entry represents the misspecified model. The $j\ (\in\\{1,2,\ldots,8\\})$ column gives the bias when the number of elements in the active set is $\hat{p}=j$ minus the bias when $\hat{p}=j-1$. On the other hand, the True1 and True2 rows give the Monte Carlo evaluated true values for the inverse-probability-weighted and doubly robust estimation, while the IPIC and DRIC rows give the average of the unbiased estimators for True1 and True2 based on the respective criteria. parameters |  | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8
---|---|---|---|---|---|---|---|---|---
(2, 4, 200, 10) | True 1 | 6.42 | 6.81 | 10.78 | 7.67 | 10.74 | 3.59 | -14.20 | 38.25
(0.1, 0.1, 0.2) | IPIC | 5.46 | 7.22 | 7.48 | 7.30 | 7.31 | 5.02 | -11.01 | 27.75
none | True2 | 6.01 | 7.37 | 9.88 | 8.37 | 5.50 | 8.38 | -9.05 | 30.99
DRIC | 5.86 | 7.79 | 7.60 | 8.12 | 5.53 | 7.50 | -8.31 | 25.88
(2, 4, 200, 10) | True 1 | 6.20 | 11.34 | 0.90 | 7.77 | 23.10 | 2.49 | -10.46 | 28.55
(0.2, 0.1, 0.2) | IPIC | 5.45 | 8.16 | 6.98 | 5.56 | 7.46 | 5.42 | -8.20 | 21.50
none | True2 | 6.19 | 7.35 | 11.01 | 6.93 | 6.99 | 10.81 | -19.32 | 37.15
DRIC | 6.00 | 9.05 | 5.12 | 8.20 | 6.92 | 5.67 | -10.58 | 24.77
(2, 4, 200, 10) | True 1 | 7.92 | 7.34 | 7.72 | 8.34 | 6.75 | 8.31 | -12.22 | 35.29
(0.1, 0.2, 0.2) | IPIC | 5.83 | 7.60 | 7.76 | 6.47 | 7.08 | 6.71 | -9.35 | 27.82
none | True2 | 7.36 | 7.30 | 10.79 | 5.54 | 8.70 | 1.15 | -9.45 | 33.94
DRIC | 6.26 | 8.39 | 8.54 | 8.68 | 7.68 | 3.04 | -9.51 | 30.42
(2, 4, 200, 10) | True 1 | 9.39 | 8.60 | 9.36 | 13.52 | 11.70 | 5.07 | -6.24 | 39.54
(0.1, 0.2, 0.4) | IPIC | 5.23 | 7.10 | 6.99 | 7.89 | 7.90 | 3.37 | -6.07 | 24.18
none | True2 | 7.30 | 8.08 | 8.82 | 9.77 | 10.09 | 8.10 | -5.10 | 31.45
DRIC | 5.56 | 8.02 | 8.05 | 7.91 | 9.24 | 6.16 | -8.38 | 27.07
(2, 4, 200, 5) | True 1 | 6.42 | 7.01 | 7.54 | 7.14 | 9.14 | 7.22 | -8.20 | 29.80
(0.1, 0.2, 0.2) | IPIC | 5.27 | 7.86 | 6.38 | 6.88 | 8.49 | 5.17 | -7.16 | 23.50
none | True2 | 5.79 | 5.54 | 9.14 | 5.84 | 9.05 | 8.67 | -7.59 | 27.87
DRIC | 6.16 | 6.52 | 8.00 | 7.87 | 7.52 | 6.67 | -7.90 | 24.72
(2, 4, 400, 10) | True 1 | 7.89 | 6.37 | 12.88 | 12.37 | 9.04 | 2.58 | -14.65 | 42.42
(0.1, 0.1, 0.2) | IPIC | 5.34 | 7.19 | 7.36 | 8.58 | 6.82 | 3.89 | -12.17 | 29.97
none | True2 | 7.72 | 10.17 | 8.75 | 10.67 | 6.83 | 11.12 | -22.45 | 42.54
DRIC | 5.73 | 8.04 | 6.98 | 7.66 | 5.54 | 5.98 | -12.15 | 30.06
(2, 4, 200, 10) | True 1 | 7.11 | 8.79 | 11.73 | 10.45 | 9.41 | 5.99 | -8.60 | 31.44
(0.1, 0.1, 0.2) | IPIC | 6.79 | 7.72 | 7.48 | 9.21 | 7.88 | 5.84 | -11.36 | 28.85
potential | True2 | 4.98 | 7.16 | 11.57 | 6.69 | 10.00 | 4.40 | -5.24 | 25.95
DRIC | 5.91 | 8.10 | 7.99 | 6.75 | 7.59 | 4.92 | -6.35 | 25.53
(2, 4, 200, 10) | True 1 | 8.34 | 7.66 | 9.51 | 10.34 | 10.02 | 10.74 | -10.35 | 36.56
(0.1, 0.1, 0.2) | IPIC | 5.74 | 7.02 | 6.25 | 8.34 | 6.99 | 6.92 | -9.05 | 24.86
outcome | True2 | 8.74 | 8.62 | 8.36 | 7.97 | 12.71 | 8.13 | -9.09 | 36.53
DRIC | 6.17 | 7.41 | 7.52 | 6.93 | 8.60 | 3.76 | -6.02 | 24.88
(2, 4, 200, 10) | True 1 | 8.21 | 15.31 | 10.08 | 11.10 | 14.62 | 3.21 | -21.76 | 39.10
(0.1, 0.1, 0.2) | IPIC | 7.15 | 7.73 | 7.53 | 6.86 | 7.16 | 3.94 | -6.48 | 21.77
both | True2 | 6.96 | 5.74 | 10.93 | 7.67 | 20.86 | 4.73 | -14.00 | 25.52
DRIC | 5.85 | 8.36 | 7.00 | 5.55 | 8.78 | 4.20 | -7.09 | 21.67

Table 4 compares the proposed criteria IPIC and DRIC with the existing
criterion QICw. Since these criteria basically try to evaluate the causal
effect $\bm{\beta}$ with good accuracy, the mean squared error of
$\hat{\bm{\beta}}$ is used as the main indicator of the comparison. How
narrowly the explanatory variables are selected by considering the number of
elements in the active set, $\hat{p}=|\hat{\mathcal{J}}_{\lambda}^{(2)}|$, is
also an important indicator. As reference indicators, we include the mean
squared error and the number of active set elements for the $pH/2$ nonzero
components in $\bm{\beta}$, as well as the mean squared error and the number
of active set elements for the $pH/2$ zero components in $\bm{\beta}$. Looking
at the main indicators, we see that IPIC and DRIC outperform QICw in all
settings, although in some cases the difference is slight. Especially when
$H=6$ or $H=8$, since the difference between the penalty terms of the
propossed and existing criteria is larger, the difference in the mean squared
error is also larger. Even if the proposed criterion does not change the mean
squared error much, it narrows down the explanatory variables considerably and
changes the estimation results significantly. Among the proposed criteria,
DRIC is superior to IPIC. This is because the doubly robust estimation is not
only robust; it also has local asymptotic efficiency. Nevertheless, there is
no significant difference in the number of elements in the active set, and
even if model misspecification exists, the degree of improvement in DRIC is
not great. That is, the benefits of DRIC for IPIC may not be as large as those
of IPIC for QICw.

Table 4: Performance comparison of QICw, IPIC and DRIC in logit causal inference models. The upper entry of the parameters is $(p,H,N,m)$, the middle entry is $(\beta^{*},\alpha^{*},\gamma^{*})$, and the lower entry represents the misspecified model. The $\hat{p}_{1}$ and $\sqrt{{\rm MSE}_{1}}$ columns list the dimension of selection and the square root of the mean squared error ($\times 10$) for nonzero parameters, while the $\hat{p}_{2}$ and $\sqrt{{\rm MSE}_{2}}$ columns show those for zero parameters and the $\hat{p}$ and $\sqrt{{\rm MSE}}$ columns show those for all parameters. parameters |  | $\hat{p}_{1}$ | $\sqrt{{\rm MSE}_{1}}$ | $\hat{p}_{2}$ | $\sqrt{{\rm MSE}_{2}}$ | $\hat{p}$ | $\sqrt{{\rm MSE}}$
---|---|---|---|---|---|---|---
(2, 4, 200, 10) | QICw | 3.6 [0.6] | 1.84 [0.70] | 3.0 [1.0] | 1.53 [0.70] | 7.0 [1.3] | 2.48 [0.73]
(0.1, 0.1, 0.2) | IPIC | 2.4 [1.3] | 2.09 [0.73] | 1.3 [1.3] | 0.88 [0.90] | 4.0 [2.5] | 2.45 [0.70]
none | DRIC | 2.4 [1.3] | 2.05 [0.73] | 1.3 [1.3] | 0.83 [0.86] | 4.0 [2.4] | 2.39 [0.67]
(2, 4, 200, 10) | QICw | 4.0 [0.2] | 1.92 [0.65] | 3.0 [0.9] | 1.59 [0.71] | 7.4 [0.9] | 2.58 [0.70]
(0.2, 0.1, 0.2) | IPIC | 3.8 [0.4] | 2.05 [0.72] | 2.1 [1.2] | 1.32 [0.85] | 6.2 [1.5] | 2.58 [0.72]
none | DRIC | 3.8 [0.4] | 1.99 [0.70] | 2.0 [1.3] | 1.26 [0.84] | 6.1 [1.6] | 2.51 [0.69]
(2, 4, 200, 10) | QICw | 3.5 [0.7] | 1.79 [0.71] | 3.0 [1.0] | 1.68 [0.69] | 6.9 [1.4] | 2.55 [0.73]
(0.1, 0.2, 0.2) | IPIC | 2.6 [1.3] | 2.00 [0.71] | 1.5 [1.3] | 1.09 [0.90] | 4.3 [2.3] | 2.46 [0.68]
none | DRIC | 2.5 [1.3] | 2.01 [0.72] | 1.4 [1.3] | 0.96 [0.83] | 4.2 [2.4] | 2.39 [0.69]
(2, 4, 200, 10) | QICw | 3.6 [0.7] | 1.99 [0.75] | 3.0 [1.1] | 1.74 [0.86] | 7.0 [1.5] | 2.77 [0.81]
(0.1, 0.1, 0.4) | IPIC | 2.7 [1.2] | 2.11 [0.71] | 1.6 [1.2] | 1.20 [0.99] | 4.5 [2.2] | 2.62 [0.73]
none | DRIC | 2.6 [1.3] | 2.02 [0.69] | 1.4 [1.4] | 0.93 [0.92] | 4.2 [2.4] | 2.41 [0.67]
(2, 4, 400, 10) | QICw | 3.8 [0.5] | 1.32 [0.49] | 3.0 [1.0] | 1.09 [0.47] | 7.2 [1.2] | 1.78 [0.47]
(0.1, 0.1, 0.2) | IPIC | 3.4 [0.8] | 1.41 [0.53] | 1.9 [1.2] | 0.85 [0.55] | 5.6 [1.8] | 1.75 [0.47]
none | DRIC | 3.4 [0.7] | 1.39 [0.52] | 1.7 [1.2] | 0.79 [0.57] | 5.5 [1.8] | 1.71 [0.47]
(2, 4, 200, 5) | QICw | 3.3 [0.8] | 2.42 [0.93] | 2.8 [1.2] | 2.08 [1.15] | 6.3 [1.6] | 3.34 [1.11]
(0.1, 0.1, 0.2) | IPIC | 1.9 [1.4] | 2.59 [0.73] | 1.1 [1.3] | 1.11 [1.31] | 3.1 [2.5] | 3.06 [0.87]
none | DRIC | 1.8 [1.5] | 2.67 [0.73] | 1.0 [1.3] | 1.06 [1.27] | 3.0 [2.6] | 3.11 [0.84]
(2, 6, 200, 10) | QICw | 5.3 [0.9] | 2.77 [0.84] | 4.9 [1.1] | 2.66 [0.98] | 10.7 [1.7] | 3.94 [0.94]
(0.1, 0.1, 0.2) | IPIC | 3.1 [1.9] | 2.81 [0.63] | 2.0 [1.8] | 1.40 [1.30] | 5.3 [3.5] | 3.35 [0.82]
none | DRIC | 3.0 [2.0] | 2.79 [0.62] | 1.8 [1.8] | 1.31 [1.27] | 5.1 [3.6] | 3.31 [0.73]
(2, 8, 200, 10) | QICw | 7.3 [0.8] | 3.92 [1.11] | 6.9 [1.1] | 3.72 [1.20] | 14.7 [1.6] | 5.50 [1.23]
(0.1, 0.1, 0.2) | IPIC | 4.0 [2.5] | 3.81 [0.91] | 2.7 [2.4] | 1.92 [1.73] | 7.1 [4.6] | 4.53 [1.23]
none | DRIC | 3.4 [2.5] | 3.81 [0.82] | 2.2 [2.3] | 1.60 [1.65] | 6.0 [4.5] | 4.39 [1.07]
(3, 4, 200, 10) | QICw | 5.6 [0.6] | 2.35 [0.68] | 4.7 [1.2] | 1.97 [0.76] | 10.7 [1.6] | 3.15 [0.72]
(0.1, 0.1, 0.2) | IPIC | 4.6 [1.3] | 2.47 [0.70] | 2.5 [1.7] | 1.33 [0.94] | 7.5 [2.8] | 2.96 [0.70]
none | DRIC | 4.6 [1.3] | 2.44 [0.71] | 2.5 [1.7] | 1.25 [0.90] | 7.4 [2.6] | 2.89 [0.70]
(4, 4, 200, 10) | QICw | 7.4 [0.9] | 2.65 [0.72] | 6.2 [1.5] | 2.42 [0.81] | 14.0 [2.0] | 3.67 [0.77]
(0.1, 0.1, 0.2) | IPIC | 5.8 [2.1] | 2.87 [0.76] | 3.7 [2.4] | 1.69 [1.12] | 9.9 [4.1] | 3.51 [0.74]
none | DRIC | 5.5 [2.2] | 2.86 [0.78] | 3.0 [2.3] | 1.42 [1.05] | 8.9 [3.9] | 3.38 [0.73]
(2, 4, 200, 10) | QICw | 3.6 [0.6] | 1.88 [0.69] | 3.0 [1.0] | 1.72 [0.90] | 7.0 [1.4] | 2.65 [0.86]
(0.1, 0.1, 0.2) | IPIC | 2.6 [1.3] | 2.09 [0.72] | 1.6 [1.4] | 1.13 [1.09] | 4.4 [2.5] | 2.59 [0.81]
potential | DRIC | 2.5 [1.3] | 2.03 [0.69] | 1.4 [1.3] | 0.95 [0.91] | 4.1 [2.5] | 2.43 [0.65]
(2, 4, 200, 10) | QICw | 3.6 [0.6] | 2.02 [0.71] | 3.1 [1.0] | 1.79 [0.85] | 7.1 [1.3] | 2.82 [0.77]
(0.1, 0.1, 0.2) | IPIC | 2.6 [1.3] | 2.13 [0.70] | 1.7 [1.4] | 1.23 [1.01] | 4.6 [2.4] | 2.67 [0.66]
outcome | DRIC | 2.6 [1.3] | 2.16 [0.72] | 1.6 [1.4] | 1.18 [1.00] | 4.4 [2.5] | 2.67 [0.69]

##  6 Real data analysis

The LaLonde dataset (LaLonde 1986) is included in the package Matching in R.
The group that took part in the U.S. job training program in 1976 is defined
as $t=1$, and the group that did not take part in the program is defined as
$t=0$. We will estimate the average treatment effect as the difference in
annual income in 1978 between the two groups after the training. Covariates
include age (age), years of education (educ), black or not (black), Hispanic
or not (hisp), married or not (married), high school graduate or not (nodegr),
annual income in year 1974 (re74), annual income in year 1975 (re75), no
income in year 1974 or not (u74), and no income in year 1975 or not (u75). The
outcome variable is the annual income in year 1978 (re78). The sample size is
$N=445$.

For these data, we use the models in Sections 5.1 and 5.2. The optimal value
of the regularization parameter is given for each criterion, and the
regression coefficients of the selected variables are obtained. Table 5(a)
compares the QICw and IPCp results when the causal eﬀect is directly estimated
for the group with ${\rm black}=0$ and for the group with ${\rm hisp}=1$,
treating the outcome variable as continuous. Although the sample size is small
as $N=74$ or even $N=39$ because of the stratification, the results of
variable selection are reasonable. While the results of variable selection for
${\rm black}=0$ are not so different between the two criteria, those for ${\rm
hisp}=1$ are very different, suggesting that the sample size is not large
enough to select variables. Table 5(b) shows whether, after transforming the
outcome variable to binary, the full sample can still be used in an analysis.
The regression coefficients for the $t=0$ group and the $t=1$ group are
estimated separately for the causal effect estimation, and the results of
QICw, IPIC, and DRIC are compared. Depending on how the transformation to
binary is conducted, the results of variable selection and estimation can vary
greatly with these criteria. Since the true structure of the data is unknown,
it is not possible to assign superiority or inferiority at this point;
however, the difference in the analysis results between the proposed and
existing criteria is considerably large.

Table 5: Analyses of LaLonde data based on each criterion. Each variable is
pre-standardized. A variable with an estimate of 0 means that it was not
selected. In (a), the case using the sample with ${\rm black}=0$ is Case 1
($N=74$) and case using the sample with ${\rm hisp}=1$ is Case 2 ($N=39$). In
(b), the whole sample is used, and the case in which the outcome variable is
dichotomous between income and non-income is Case 1
($\sum_{i=1}^{445}t_{i}=308$, $\sum_{i=1}^{445}(1-t_{i})=137$), while the case
in which the outcome variable is dichotomous between above 1,500 and below
1,500 is Case 2 ($\sum_{i=1}^{445}t_{i}=270$,
$\sum_{i=1}^{445}(1-t_{i})=175$).

(a) Estimated regression coefficients for continuous outcome variables
($\times 10$).

|  | age | educ | hisp | marri | nodeg | re74 | re75 | u74 | u75
---|---|---|---|---|---|---|---|---|---|---
Case 1 | QICw | 3.97 | 1.08 | 0.00 | 0.00 | 4.42 | -2.06 | 0.00 | 0.34 | 1.27
IPCp | 3.61 | 0.00 | 0.00 | 0.00 | 3.30 | -1.83 | 0.00 | 0.16 | 1.10
Case 2 | QICw | 13.03 | 0.00 | —– | -2.77 | 0.00 | -5.04 | -4.83 | 0.00 | 0.00
IPCp | 4.19 | 0.00 | —– | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00

(b) Estimated regression coefficients for binary outcome variable ($\times
10$).

|  |  | age | educ | black | hisp | marri | nodeg | u74 | u75
---|---|---|---|---|---|---|---|---|---|---
| QICw | case | 0.00 | 0.00 | -2.21 | 0.47 | 0.98 | -0.23 | 0.00 | -1.08
| control | 0.00 | 0.00 | -2.27 | 0.00 | 0.00 | 0.00 | -1.56 | 0.00
Case 1 | IPIC | case | 0.00 | 0.00 | -2.22 | 0.56 | 1.07 | -0.33 | 0.21 | -1.30
control | 0.00 | 0.00 | -2.36 | 0.00 | 0.00 | 0.00 | -1.77 | 0.18
| DRIC | case | 0.00 | 0.00 | -1.81 | 0.02 | 0.35 | 0.00 | 0.00 | -0.51
| control | 0.00 | 0.00 | -1.55 | 0.00 | 0.00 | 0.00 | -0.85 | 0.00
| QICw | case | 0.00 | 0.00 | -2.51 | 0.00 | 1.18 | -0.56 | 0.82 | -0.57
| control | 0.00 | 0.00 | -3.27 | 0.00 | 0.00 | 0.00 | -1.11 | 0.00
Case 2 | IPIC | case | 0.00 | 0.00 | -1.97 | 0.00 | 0.59 | 0.00 | 0.00 | 0.00
control | 0.00 | 0.00 | -2.68 | 0.00 | 0.00 | 0.00 | -0.56 | 0.00
| DRIC | case | 0.00 | 0.00 | -1.48 | 0.00 | 0.12 | 0.00 | 0.00 | 0.00
| control | 0.00 | 0.00 | -2.19 | 0.00 | 0.00 | 0.00 | -0.08 | 0.00

##  7 Extension

###  7.1 Generalization of criterion based on SURE theory

Although the outcome variable is not basically supposed to follow a Gaussian
distribution in propensity score analysis, we can apply SURE theory by
conditioning on the assignment variable earlier. This application is valid for
sparse estimations not limited to LASSO. Here, we will introduce some of them.

First, let us replace the $\ell_{1}$ norm of the penalty term
$\|\bm{\theta}\|_{1}$ with another convex function that induces a different
kind of sparsity. Specifically, consider the group LASSO estimator given by

$\displaystyle(\hat{\bm{\theta}}_{\lambda}^{[1]},\hat{\bm{\theta}}_{\lambda}^{[2]},\ldots,\hat{\bm{\theta}}_{\lambda}^{[H]})\equiv\mathop{\rm
argmin}_{\bm{\theta}^{[1]},\bm{\theta}^{[2]},\ldots,\bm{\theta}^{[H]}}\Bigg{[}\sum_{i=1}^{N}\sum_{h=1}^{H}\Bigg{\\{}\frac{t_{i}^{[h]}y_{i}^{[h]}}{e^{[h]}(\bm{z}_{i})}-\bm{x}_{i}^{{\rm
T}}\bm{\theta}^{[h]}\Bigg{\\}}^{2}+\lambda\sum_{j=1}^{p}\Bigg{(}\sum_{h=1}^{H}\theta_{j}^{[h]2}\Bigg{)}^{1/2}\Bigg{]}$

in the model of (2) (Yuan and Lin 2006), as it can be expected to be useful.
In group LASSO, the parameters in the group tend to be zero at once, and the
fact that $\theta_{j}^{[1]},\theta_{j}^{[2]},\ldots,\theta_{j}^{[H]}$ are zero
means that the $j$-th explanatory variable $x_{j}$ does not affect the outcome
variable.

In this estimator, we denote the collection of subscripts $j$ such that
$(\hat{\theta}_{\lambda,j}^{[1]},\hat{\theta}_{\lambda,j}^{[2]},\ldots,\hat{\theta}_{\lambda,j}^{[H]})$
is $\bm{0}_{H}$ as $\hat{\mathcal{J}}_{\lambda}^{(1)}$ and the collection of
other subscripts as $\hat{\mathcal{J}}_{\lambda}^{(2)}$. According to this
notation,
$(\hat{\theta}_{\lambda,j}^{[h]})_{j\in\hat{\mathcal{J}}_{\lambda}^{(1)}}$,
$(\hat{\theta}_{\lambda,j}^{[h]})_{j\in\hat{\mathcal{J}}_{\lambda}^{(2)}}$,
$(x_{ij})_{j\in\hat{\mathcal{J}}_{\lambda}^{(2)}}$, etc. are written as
$\hat{\bm{\theta}}_{\lambda}^{[h](1)}$,
$\hat{\bm{\theta}}_{\lambda}^{[h](2)}$, $\bm{x}_{i}^{(2)}$, etc.,
respectively. Of course, $\hat{\bm{\theta}}_{\lambda}^{[h](1)}$ is
$\bm{0}_{|\hat{\mathcal{J}}_{\lambda}^{(1)}|}$, and differentiation with
respect to $\bm{\theta}^{[h](2)}$ gives

$\displaystyle\hat{\bm{\theta}}_{\lambda}^{[h](2)}=\Bigg{(}\sum_{i=1}^{N}\bm{x}_{i}^{(2)}\bm{x}_{i}^{(2){\rm
T}}\Bigg{)}^{-1}\Bigg{\\{}\sum_{i=1}^{N}\frac{t_{i}^{[h]}y_{i}^{[h]}\bm{x}_{i}^{(2)}}{e^{[h]}(\bm{z}_{i})}+\lambda\Bigg{(}\Bigg{(}\sum_{h^{\prime}=1}^{H}\hat{\theta}_{\lambda,j}^{[h^{\prime}]2}\Bigg{)}^{-1/2}\hat{\theta}_{\lambda,j}^{[h]}\Bigg{)}_{j\in\hat{\mathcal{J}}_{\lambda}^{(2)}}\Bigg{\\}}.$
(28)

As in Section 3, we consider a risk based on the loss function used in this
estimation and decompose it into three terms, removing the expectation in the
first term and ignoring the second term. The third term can be rewritten as

$\displaystyle{\rm
E}\Bigg{(}\sum_{i=1}^{N}\sum_{h=1}^{H}\Bigg{[}\frac{t_{i}^{[h]}y_{i}^{[h]}}{e^{[h]}(\bm{z}_{i})}-{\rm
E}\Bigg{\\{}\frac{t_{i}^{[h]}y_{i}^{[h]}}{e^{[h]}(\bm{z}_{i})}\ \Bigg{|}\
t_{i}^{[h]},\bm{x}_{i},\bm{z}_{i}\Bigg{\\}}\Bigg{]}\bm{x}_{i}^{{\rm
T}}\hat{\bm{\theta}}_{\lambda}^{[h]}\Bigg{)}$ $\displaystyle={\rm
E}\Bigg{(}\sum_{i=1}^{N}\sum_{h=1}^{H}\Bigg{[}\frac{t_{i}^{[h]}}{e^{[h]}(\bm{z}_{i})}\\{y_{i}^{[h]}-\bm{x}_{i}^{{\rm
T}}\bm{\theta}^{[h]}-\eta(\bm{z}_{i})\\}\bm{x}_{i}^{{\rm
T}}\hat{\bm{\theta}}_{\lambda}^{[h]}\Bigg{]}\Bigg{)},$ (29)

and by conditioning on $(t_{i}^{[h]},\bm{x}_{i},\bm{z}_{i})$ earlier, we use
SURE theory. Since it follows from (28) that

$\displaystyle\frac{\partial\hat{\bm{\theta}}_{\lambda}^{[h^{\prime}](2)}}{\partial
y_{i}^{[h]}}=\ $
$\displaystyle\Bigg{(}\sum_{i^{\prime}=1}^{N}\bm{x}_{i^{\prime}}^{(2)}\bm{x}_{i^{\prime}}^{(2){\rm
T}}\Bigg{)}^{-1}\Bigg{\\{}\frac{\delta_{\\{h=h^{\prime}\\}}t_{i}^{[h]}\bm{x}_{i}^{(2)}}{e^{[h]}(\bm{z}_{i})}$
$\displaystyle+\lambda\Bigg{(}\Bigg{(}\sum_{h^{\prime\prime}=1}^{H}\hat{\theta}_{\lambda,j}^{[h^{\prime\prime}]2}\Bigg{)}^{-1/2}\frac{\partial\hat{\theta}_{\lambda,j}^{[h^{\prime}]}}{\partial
y_{i}^{[h]}}-\Bigg{(}\sum_{h^{\prime\prime}=1}^{H}\hat{\theta}_{\lambda,j}^{[h^{\prime\prime}]2}\Bigg{)}^{-3/2}\sum_{h^{\prime\prime}=1}^{H}\Bigg{(}\hat{\theta}_{\lambda,j}^{[h^{\prime\prime}]}\frac{\partial\hat{\theta}_{\lambda,j}^{[h^{\prime\prime}]}}{\partial
y_{i}^{[h]}}\Bigg{)}\hat{\theta}_{\lambda,j}^{[h^{\prime}]}\Bigg{)}_{j\in\hat{\mathcal{J}}_{\lambda}^{(2)}}\Bigg{\\}},$

if we write

$\displaystyle\bm{R}_{h,h}^{(2)}\equiv\bm{I}_{|\hat{\mathcal{J}}_{\lambda}^{(2)}|}-\lambda\Bigg{(}\sum_{i^{\prime}=1}^{N}\bm{x}_{i^{\prime}}^{(2)}\bm{x}_{i^{\prime}}^{(2){\rm
T}}\Bigg{)}^{-1}{\rm
diag}\Bigg{\\{}\Bigg{(}\sum_{h^{\prime}=1}^{H}\hat{\theta}_{\lambda,j}^{[h^{\prime}]2}\Bigg{)}^{-1/2}-\Bigg{(}\sum_{h^{\prime}=1}^{H}\hat{\theta}_{\lambda,j}^{[h^{\prime}]2}\Bigg{)}^{-3/2}\hat{\theta}_{\lambda,j}^{[h]2}\Bigg{\\}}_{j\in\hat{\mathcal{J}}_{\lambda}^{(2)}}$

and

$\displaystyle\bm{R}_{h,h^{\prime}}^{(2)}\equiv\lambda\Bigg{(}\sum_{i^{\prime}=1}^{N}\bm{x}_{i^{\prime}}^{(2)}\bm{x}_{i^{\prime}}^{(2){\rm
T}}\Bigg{)}^{-1}{\rm
diag}\Bigg{\\{}\Bigg{(}\sum_{h^{\prime}=1}^{H}\hat{\theta}_{\lambda,j}^{[h^{\prime}]2}\Bigg{)}^{-3/2}\hat{\theta}_{\lambda,j}^{[h]}\hat{\theta}_{\lambda,j}^{[h^{\prime}]}\Bigg{\\}}_{j\in\hat{\mathcal{J}}_{\lambda}^{(2)}},\qquad
h\neq h^{\prime},$

and let
$\bm{R}^{(2)}\equiv(\bm{R}_{h,h^{\prime}}^{(2)})_{h,h^{\prime}\in\\{1,2,\ldots,H\\}}$
be the
$|\hat{\mathcal{J}}_{\lambda}^{(2)}|H\times|\hat{\mathcal{J}}_{\lambda}^{(2)}|H$
matrix, we can rewrite the derivative as

$\displaystyle\frac{\partial\hat{\bm{\theta}}_{\lambda}^{[h^{\prime}](2)}}{\partial
y_{i}^{[h]}}=(\bm{I}_{|\hat{\mathcal{J}}_{\lambda}^{(2)}|}\otimes\bm{e}_{h})^{{\rm
T}}\bm{R}^{(2)-1}\Bigg{[}\Bigg{\\{}\Bigg{(}\sum_{i^{\prime}=1}^{N}\bm{x}_{i^{\prime}}^{(2)}\bm{x}_{i^{\prime}}^{(2){\rm
T}}\Bigg{)}^{-1}\frac{t_{i}^{[h]}\bm{x}_{i}^{(2)}}{e^{[h]}(\bm{z}_{i})}\Bigg{\\}}\otimes\bm{e}_{h}\Bigg{]},$

where $\bm{e}_{h}$ is the unit vector with $1$ in the $h$-th component.
Therefore, (29) is evaluated to be

$\displaystyle{\rm E}\Bigg{\\{}\sum_{i=1}^{N}\sum_{h=1}^{H}{\rm
E}\Bigg{(}\frac{t_{i}^{[h]}}{e^{[h]}(\bm{z}_{i})}\bm{x}_{i}^{(2){\rm
T}}(\bm{I}_{|\hat{\mathcal{J}}_{\lambda}^{(2)}|}\otimes\bm{e}_{h})^{{\rm
T}}\bm{R}^{(2)-1}$ $\displaystyle\phantom{{\rm
E}\Bigg{[}\sum_{i=1}^{N}\sum_{h=1}^{H}{\rm
E}\Bigg{\\{}}\Bigg{[}\Bigg{\\{}\Bigg{(}\sum_{i^{\prime}=1}^{N}\bm{x}_{i^{\prime}}^{(2)}\bm{x}_{i^{\prime}}^{(2){\rm
T}}\Bigg{)}^{-1}\frac{t_{i}^{[h]}\bm{x}_{i}^{(2)}}{e^{[h]}(\bm{z}_{i})}\Bigg{\\}}\otimes\bm{e}_{h}\Bigg{]}\
\Bigg{|}\ t_{i}^{[h]},\bm{x}_{i},\bm{z}_{i}\Bigg{)}\Bigg{\\}}$
$\displaystyle={\rm E}\Bigg{[}{\rm
tr}\Bigg{\\{}\sum_{h=1}^{H}(\bm{I}_{|\hat{\mathcal{J}}_{\lambda}^{(2)}|}\otimes\bm{e}_{h})^{{\rm
T}}\bm{R}^{(2)-1}\Bigg{(}\Bigg{[}\Bigg{(}\sum_{i=1}^{N}\bm{x}_{i}^{(2)}\bm{x}_{i}^{(2){\rm
T}}\Bigg{)}^{-1}\Bigg{\\{}\sum_{i=1}^{N}\frac{\bm{x}_{i}^{(2)}\bm{x}_{i}^{(2){\rm
T}}}{e^{[h]}(\bm{z}_{i})}\Bigg{\\}}\Bigg{]}\otimes\bm{e}_{h}\Bigg{)}\Bigg{\\}}\Bigg{]}.$

Accordingly, we can derive

$\displaystyle{\rm IPCp}\equiv$
$\displaystyle\sum_{i=1}^{N}\sum_{h=1}^{H}\Bigg{\\{}\frac{t_{i}^{[h]}y_{i}^{[h]}}{e^{[h]}(\bm{z}_{i})}-\bm{x}_{i}^{(2){\rm
T}}\hat{\bm{\theta}}_{\lambda}^{[h](2)}\Bigg{\\}}^{2}$
$\displaystyle+2\sigma^{2}{\rm
tr}\Bigg{\\{}\sum_{h=1}^{H}(\bm{I}_{|\hat{\mathcal{J}}_{\lambda}^{(2)}|}\otimes\bm{e}_{h})^{{\rm
T}}\bm{R}^{(2)-1}\Bigg{(}\Bigg{[}\Bigg{(}\sum_{i=1}^{N}\bm{x}_{i}^{(2)}\bm{x}_{i}^{(2){\rm
T}}\Bigg{)}^{-1}\Bigg{\\{}\sum_{i=1}^{N}\frac{\bm{x}_{i}^{(2)}\bm{x}_{i}^{(2){\rm
T}}}{e^{[h]}(\bm{z}_{i})}\Bigg{\\}}\Bigg{]}\otimes\bm{e}_{h}\Bigg{)}\Bigg{\\}}$

as an unbiased estimator of the weighted mean squared error.

Next, let us consider the elastic net, where there are multiple tuning
parameters and the solution of the estimating equations is not necessarily
used directly as the estimator. In the elastic net, the setting in which
$\bm{\theta}$ is estimated directly, as in Section 3, is expected to be used
to the same extent as the setting in which each $\bm{\theta}^{[h]}$ is
estimated individually. We will treat the former setting, which has simpler
notation. That is, as $\hat{\bm{\theta}}_{\lambda_{1},\lambda_{2}}$, we use
the estimator that multiplies

$\displaystyle\mathop{\rm
argmin}_{\bm{\theta}}\Bigg{[}\sum_{i=1}^{N}\Bigg{\\{}\sum_{h=1}^{H}\frac{c^{[h]}t_{i}^{[h]}y_{i}^{[h]}}{e^{[h]}(\bm{z}_{i})}-\bm{x}_{i}^{{\rm
T}}\bm{\theta}\Bigg{\\}}^{2}+\lambda_{1}\|\bm{\theta}\|_{1}+\lambda_{2}\|\bm{\theta}\|_{2}^{2}\Bigg{]}$
(30)

by $(1+\lambda_{2})$ in (2). This estimator consists of
$\hat{\bm{\theta}}_{\lambda_{1},\lambda_{2}}^{(1)}=\bm{0}_{|\hat{\mathcal{J}}_{\lambda}^{(1)}|}$
and

$\displaystyle\hat{\bm{\theta}}_{\lambda_{1},\lambda_{2}}^{(2)}=(1+\lambda_{2})\Bigg{(}\sum_{i=1}^{N}\bm{x}_{i}^{(2)}\bm{x}_{i}^{(2){\rm
T}}+\lambda_{2}\bm{I}_{|\hat{\mathcal{J}}_{\lambda}^{(2)}|}\Bigg{)}^{-1}\Bigg{\\{}\sum_{i=1}^{N}\sum_{h=1}^{H}\frac{c^{[h]}t_{i}^{[h]}y_{i}^{[h]}\bm{x}_{i}^{(2)}}{e^{[h]}(\bm{z}_{i})}+\lambda_{1}\hskip
1.42262pt{\rm sgn}(\hat{\bm{\theta}}_{\lambda}^{(2)})\Bigg{\\}}.$

By generalizing Zou et al. (2007), it can be shown that SURE theory holds even
if there are two tuning parameters. See the Appendix for the details.
According to the theory, if we use

$\displaystyle\frac{\partial\hat{\bm{\theta}}_{\lambda_{1},\lambda_{2}}^{(2)}}{\partial
y_{i}^{[h]}}=\Bigg{(}\sum_{i=1}^{N}\bm{x}_{i}^{(2)}\bm{x}_{i}^{(2){\rm
T}}+\lambda_{2}\bm{I}_{|\hat{\mathcal{J}}_{\lambda}^{(2)}|}\Bigg{)}^{-1}\frac{(1+\lambda_{2})c^{[h]}t_{i}^{[h]}\bm{x}_{i}^{(2)}}{e^{[h]}(\bm{z}_{i})}$

in the third term of the decomposition of the risk, it is evaluated as

$\displaystyle{\rm E}\Bigg{(}{\rm tr}\Bigg{[}\sum_{i=1}^{N}\sum_{h=1}^{H}{\rm
E}\Bigg{\\{}\Bigg{(}\sum_{i^{\prime}=1}^{N}\bm{x}_{i^{\prime}}^{(2)}\bm{x}_{i^{\prime}}^{(2){\rm
T}}+\lambda_{2}\bm{I}_{|\hat{\mathcal{J}}_{\lambda}^{(2)}|}\Bigg{)}^{-1}\frac{(1+\lambda_{2})c^{[h]2}\bm{x}_{i}^{(2)}\bm{x}_{i}^{(2){\rm
T}}}{e^{[h]}(\bm{z}_{i})}\ \Bigg{|}\
\bm{x}_{i},\bm{z}_{i}\Bigg{\\}}\Bigg{]}\Bigg{)}$ $\displaystyle={\rm
E}\Bigg{(}{\rm
tr}\Bigg{[}\Bigg{(}\sum_{i=1}^{N}\bm{x}_{i}^{(2)}\bm{x}_{i}^{(2){\rm
T}}+\lambda_{2}\bm{I}_{|\hat{\mathcal{J}}_{\lambda}^{(2)}|}\Bigg{)}^{-1}\Bigg{\\{}\sum_{i=1}^{N}\sum_{h=1}^{H}\frac{(1+\lambda_{2})c^{[h]2}\bm{x}_{i}^{(2)}\bm{x}_{i}^{(2){\rm
T}}}{e^{[h]}(\bm{z}_{i})}\Bigg{\\}}\Bigg{]}\Bigg{)}.$

Consequently, we derive

$\displaystyle{\rm IPCp}\equiv\ $
$\displaystyle\sum_{i=1}^{N}\Bigg{\\{}\sum_{h=1}^{H}\frac{c^{[h]}t_{i}^{[h]}y_{i}^{[h]}}{e^{[h]}(\bm{z}_{i})}-\bm{x}_{i}^{(2){\rm
T}}\hat{\bm{\theta}}_{\lambda_{1},\lambda_{2}}^{(2)}\Bigg{\\}}^{2}$
$\displaystyle+2\sigma^{2}{\rm
tr}\Bigg{[}\Bigg{(}\sum_{i=1}^{N}\bm{x}_{i}^{(2)}\bm{x}_{i}^{(2){\rm
T}}+\lambda_{2}\bm{I}_{|\hat{\mathcal{J}}_{\lambda}^{(2)}|}\Bigg{)}^{-1}\Bigg{\\{}\sum_{i=1}^{N}\sum_{h=1}^{H}\frac{(1+\lambda_{2})c^{[h]2}\bm{x}_{i}^{(2)}\bm{x}_{i}^{(2){\rm
T}}}{e^{[h]}(\bm{z}_{i})}\Bigg{\\}}\Bigg{]}$

as an unbiased estimator of the weighted mean squared error.

###  7.2 Generalization of criterion based on asymptotics

The extensions described in the previous subsection assumed global convexity
of the objective function. On the other hand, sparse estimations using
nonconvex penalties such as SCAD (Fan and Li 2001) and MC+ (Zhang 2010) are in
use these days. In this subsection, we confirm that we can use asymptotics to
derive a valid information criterion for the causal inference version of
sparse estimation using nonconvex penalties under some mild conditions.

As in the setting of Section 4.1, we will treat a nonconvex penalty,
$\rho_{\lambda}(\bm{\beta})$, instead of $\lambda\|\bm{\beta}\|_{1}$. If
$\rho_{\lambda}(\bm{\beta})$ is monotonically increasing with respect to
$|\beta_{j}|$ ($j\in\\{1,2,\ldots,p\\}$), then $\hat{\bm{\beta}}_{\lambda}$
converges in probability to

$\displaystyle\bm{\beta}_{\lambda}^{*}\equiv\mathop{\rm
argmin}_{\bm{\beta}}\Bigg{[}-{\rm E}\Bigg{\\{}\sum_{h=1}^{H}\log
f(y_{i}^{[h]}\mid\bm{x}_{i}^{[h]};\bm{\beta})\Bigg{\\}}+\rho_{\lambda}(\bm{\beta})\Bigg{]}.$

Then, if
$\bm{J}(\bm{\beta})+\partial^{2}\rho_{\lambda}(\bm{\beta})/\partial\bm{\beta}\partial\bm{\beta}^{{\rm
T}}$ is positive and continuous at $\bm{\beta}=\bm{\beta}_{\lambda}^{*}$, we
obtain

$\displaystyle N\hat{\bm{\beta}}_{\lambda}^{(1)}={\rm o}_{\rm P}(1)$

and

$\displaystyle\sqrt{N}(\hat{\bm{\beta}}_{\lambda}^{(2)}-\bm{\beta}_{\lambda}^{*(2)})=$
$\displaystyle\Bigg{\\{}\bm{J}^{(22)}(\bm{\beta}_{\lambda}^{*})+\frac{\partial^{2}}{\partial\bm{\beta}^{(2)}\partial\bm{\beta}^{(2){\rm
T}}}\rho_{\lambda}(\bm{\beta}_{\lambda}^{*})\Bigg{\\}}^{-1}$
$\displaystyle\Bigg{\\{}\frac{1}{\sqrt{N}}\sum_{i=1}^{N}\sum_{h=1}^{H}\frac{\partial}{\partial\bm{\beta}^{(2)}}\frac{t_{i}^{[h]}}{e^{[h]}(\bm{z}_{i})}\log
f(y_{i}^{[h]}\mid\bm{x}_{i}^{[h]};\bm{\beta}_{\lambda}^{*})-\sqrt{N}\frac{\partial}{\partial\bm{\beta}^{(2)}}\rho_{\lambda}(\bm{\beta}_{\lambda}^{*})\Bigg{\\}}+{\rm
o}_{\rm P}(1)$

similarly to (12) and (13) by following the same argument as in Knight and Fu
(2000). Therefore, if $b^{\rm limit}$ is defined as it is in Section 4.1, it
reduces to

$\displaystyle\frac{1}{N}\frac{\partial}{\partial\bm{\beta}^{(2){\rm
T}}}\Bigg{\\{}\sum_{i=1}^{N}\sum_{h=1}^{H}\frac{t_{i}^{[h]}}{e^{[h]}(\bm{z}_{i})}\log
f(y_{i}^{[h]}\mid\bm{x}_{i}^{[h]};\bm{\beta}_{\lambda}^{*})-N\rho_{\lambda}(\bm{\beta}_{\lambda}^{*})\Bigg{\\}}$
$\displaystyle\Bigg{\\{}\bm{J}^{(22)}(\bm{\beta}_{\lambda}^{*})+\frac{\partial^{2}}{\partial\bm{\beta}^{(2)}\partial\bm{\beta}^{(2){\rm
T}}}\rho_{\lambda}(\bm{\beta}_{\lambda}^{*})\Bigg{\\}}^{-1}\frac{\partial}{\partial\bm{\beta}^{(2)}}\Bigg{\\{}\sum_{i^{\prime}=1}^{N}\sum_{h^{\prime}=1}^{H}\frac{t_{i^{\prime}}^{[h^{\prime}]}}{e^{[h^{\prime}]}(\bm{z}_{i^{\prime}})}\log
f(y_{i^{\prime}}^{[h^{\prime}]}\mid\bm{x}_{i^{\prime}}^{[h^{\prime}]};\bm{\beta}_{\lambda}^{*})-N\rho_{\lambda}(\bm{\beta}_{\lambda}^{*})\Bigg{\\}}.$

From independence, the expectation of this is $0$ when $i\neq i^{\prime}$,
$t_{i}^{[h]}t_{i}^{[h^{\prime}]}=0$ when $h\neq h^{\prime}$, and
$t_{i}^{[h]2}=t_{i}^{[h]}$ when $h=h^{\prime}$. Thus, ${\rm E}(b^{\rm limit})$
can be written as

$\displaystyle\Bigg{\\{}\bm{J}^{(22)}(\bm{\beta}_{\lambda}^{*})+\frac{\partial^{2}}{\partial\bm{\beta}^{(2)}\partial\bm{\beta}^{(2){\rm
T}}}\rho_{\lambda}(\bm{\beta}_{\lambda}^{*})\Bigg{\\}}^{-1}$
$\displaystyle\Bigg{[}\frac{\partial}{\partial\bm{\beta}^{(2)}}\rho_{\lambda}(\bm{\beta}_{\lambda}^{*})\frac{\partial}{\partial\bm{\beta}^{(2){\rm
T}}}\rho_{\lambda}(\bm{\beta}_{\lambda}^{*})+\sum_{h=1}^{H}{\rm
E}\Bigg{\\{}\frac{1}{e^{[h]}(\bm{z})}\frac{\partial}{\partial\bm{\beta}}\log
f(y^{[h]}\mid\bm{x}^{[h]};\bm{\beta}_{\lambda}^{*})\frac{\partial}{\partial\bm{\beta}^{{\rm
T}}}\log f(y^{[h]}\mid\bm{x}^{[h]};\bm{\beta}_{\lambda}^{*})\Bigg{\\}}$
$\displaystyle-\sum_{h=1}^{H}{\rm
E}\Bigg{\\{}\frac{\partial}{\partial\bm{\beta}^{(2)}}\log
f(y^{[h]}\mid\bm{x}^{[h]};\bm{\beta}_{\lambda}^{*})\frac{\partial}{\partial\bm{\beta}^{(2){\rm
T}}}\rho_{\lambda}(\bm{\beta}_{\lambda}^{*})+\frac{\partial}{\partial\bm{\beta}^{(2)}}\rho_{\lambda}(\bm{\beta}_{\lambda}^{*})\frac{\partial}{\partial\bm{\beta}^{(2){\rm
T}}}\log
f(y^{[h]}\mid\bm{x}^{[h]};\bm{\beta}_{\lambda}^{*})\Bigg{\\}}\Bigg{]}.$

Consequently, the information criterion for the inverse-probability-weighted
sparse estimation has a form similar to that of IPIC.

##  8 Conclusion

For a causal inference model in which the outcome variable $y$ follows a
Gaussian distribution conditional on the assignment variable $t$ and the
covariate $(\bm{x},\bm{z})$, we have shown that SURE theory can be extended by
conditioning on $(t,\bm{x},\bm{z})$ earlier, which is not usually considered
in propensity score analysis. On the basis of this extension, we have derived
a generalized Cp criterion for causal inference, IPCp, and we have shown that
$\sum_{i=1}^{N}\bm{x}_{i}\bm{x}_{i}^{{\rm T}}$ in the conventional generalized
Cp criterion becomes $\sum_{i=1}^{N}\sum_{h=1}^{H}\bm{x}_{i}\bm{x}_{i}^{{\rm
T}}/\\{He^{[h]}(\bm{z}_{i})\\}$ in IPCp, which leads to a substantially
different result. Next, we have extended the asymptotic theory of sparse
estimation in Knight and Fu (2000) for propensity score analysis for the case
where $y$ does not follow a Gaussian distribution even after conditioning and
for the case of doubly robust estimation. Then, we have derived an information
criterion (IPIC) for inverse-probability-weighted estimation and an
information criterion (DRIC) for doubly robust estimation that is doubly
robust in itself from a divergence customized for causal inference. When the
conditional distribution of $y$ is Gaussian and there is no model
misspecification, the penalty term in DRIC is equal to the penalty term in
IPCp, which is basically the main quantity. This means that the results of
analyses using the proposed criteria will differ significantly from those of
the existing criteria. Numerical experiments have confirmed that the accuracy
of IPCp is very high and that of DRIC is moderately high as an approximation
of risk. Furthermore, we have compared IPCp and DRIC with the existing
criterion and found that they reduced the mean square error in almost all
cases and that the difference is clear in many cases. As well, an analysis
conducted on a real data set have confirmed that the proposed criteria show
large differences from the existing criteria in terms of variable selection
and estimates, suggesting again the significance of using the proposed
criteria.

To explore possible generalizations of the proposed criteria, we have derived
IPCp and IPIC for general sparse estimation using group LASSO, elastic net,
and nonconvex regularization. These extensions proved feasible, so we plan to
conduct further customizations of the method in the future. For example, it is
important to make confounding and assignment time-dependent and to handle
dynamic treatment regimens (see e.g. Chakraborty and Moodie 2013, Tsiatis et
al. 2019). Also, in regard to dynamic treatment regimens, it is necessary to
find out what covariates should be modeled for the outcome variable, and since
the number of variables tends to be large, model selection is a critical
issue. In addition, the topic of estimating causal effects using propensity
score analysis and sparse estimation when the covariates are high dimensional
is currently a popular one in econometrics (e.g., Belloni et al. 2017,
Chernozhukov et al. 2018, Athey et al. 2018, Ning et al. 2020); model
selection would be necessary in this setting for efficiency.

## Acknowledgements

This research was supported by a JSPS Grant-in-Aid for Scientific Research
(16K00050).

## References

  * (1)
  * Akaike (1973) Akaike, H. (1973). Information theory and an extension of the maximum likelihood principle, In 2nd Intr. Symp. on Information Theory, Budapest, 1973, Akademiai Kiado.
  * Andersen and Gill (1982) Andersen, P. K. and Gill, R. D. (1982). Cox’s regression model for counting processes: a large sample study, The Annals of Statistics, 1100–1120.
  * Athey et al. (2018) Athey, S., Imbens, G., and Wager, S. (2018). Approximate residual balancing: debiased inference of average treatment effects in high dimensions, Journal of the Royal Statistical Society: Series B, 80, 597–623.
  * Baba et al. (2017) Baba, T., Kanemori, T., and Ninomiya, Y. (2017). A Cp criterion for semiparametric causal inference, Biometrika, 104, 845–861.
  * Baba and Ninomiya (2021) Baba, T. and Ninomiya, Y. (2021). Doubly robust criterion for causal inference, arXiv preprint, 2110.14525.
  * Bang and Robins (2005) Bang, H. and Robins, J. M. (2005). Doubly robust estimation in missing data and causal inference models, Biometrics, 61, 962–972.
  * Belloni et al. (2017) Belloni, A., Chernozhukov, V., Fernández-Val, I., and Hansen, C. (2017). Program evaluation and causal inference with high-dimensional data, Econometrica, 85, 233–298.
  * Chakraborty and Moodie (2013) Chakraborty, B. and Moodie, E. E. (2013). Statistical Methods for Dynamic Treatment Regimes: Springer.
  * Chernozhukov et al. (2018) Chernozhukov, V., Chetverikov, D., Demirer, M., Duflo, E., Hansen, C., Newey, W., and Robins, J. (2018). Double/debiased machine learning for treatment and structural parameters, The Econometrics Journal, 21.
  * Dukes and Vansteelandt (2021) Dukes, O. and Vansteelandt, S. (2021). Inference for treatment effect parameters in potentially misspecified high-dimensional models, Biometrika, 108, 321–334.
  * Efron et al. (2004) Efron, B., Hastie, T., Johnstone, I., and Tibshirani, R. (2004). Least angle regression, The Annals of Statistics, 32, 407–499.
  * Fan and Li (2001) Fan, J. and Li, R. (2001). Variable selection via nonconcave penalized likelihood and its oracle properties, Journal of the American Statistical Association, 96, 1348–1360.
  * Hernán and Robins (2020) Hernán, M. A. and Robins, J. M. (2020). Causal inference: what if: Boca Raton: Chapman & Hall/CRC.
  * Hjort and Pollard (2011) Hjort, N. L. and Pollard, D. (2011). Asymptotics for minimisers of convex processes, arXiv preprint, 1107.3806.
  * Knight and Fu (2000) Knight, K. and Fu, W. (2000). Asymptotics for lasso-type estimators, The Annals of Statistics, 28, 1356–1378.
  * Konishi and Kitagawa (1996) Konishi, S. and Kitagawa, G. (1996). Generalised information criteria in model selection, Biometrika, 83, 875–890.
  * LaLonde (1986) LaLonde, R. J. (1986). Evaluating the econometric evaluations of training programs with experimental data, The American Economic Review, 604–620.
  * Lockhart et al. (2014) Lockhart, R., Taylor, J., Tibshirani, R. J., and Tibshirani, R. (2014). A significance test for the lasso, Annals of Statistics, 42, 413.
  * Ning et al. (2020) Ning, Y., Sida, P., and Imai, K. (2020). Robust estimation of causal effects via a high-dimensional covariate balancing propensity score, Biometrika, 107, 533–554.
  * Ninomiya and Kawano (2016) Ninomiya, Y. and Kawano, S. (2016). AIC for the Lasso in generalized linear models, Electronic Journal of Statistics, 10, 2537–2560.
  * Peters et al. (2017) Peters, J., Janzing, D., and Schölkopf, B. (2017). Elements of causal inference: foundations and learning algorithms: The MIT Press.
  * Platt et al. (2013) Platt, R. W., Brookhart, M. A., Cole, S. R., Westreich, D., and Schisterman, E. F. (2013). An information criterion for marginal structural models, Statistics in Medicine, 32, 1383–1393.
  * Pollard (1991) Pollard, D. (1991). Asymptotics for least absolute deviation regression estimators, Econometric Theory, 7, 186–199.
  * 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.
  * Rolling and Yang (2014) Rolling, C. A. and Yang, Y. (2014). Model selection for estimating treatment effects, Journal of the Royal Statistical Society: Series B, 749–769.
  * 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.
  * Shibata (1989) Shibata, R. (1989). Statistical aspects of model selection, In From Data to Model: Springer, 215–240.
  * Spiegelhalter et al. (2002) Spiegelhalter, D. J., Best, N. G., Carlin, B. P., and van der Linde, A. (2002). Bayesian measures of model complexity and fit, Journal of the Royal Statistical Society: Series B, 64, 583–639.
  * Stein (1981) Stein, C. M. (1981). Estimation of the mean of a multivariate normal distribution, The Annals of Statistics, 1135–1151.
  * Tibshirani (1996) Tibshirani, R. (1996). Regression shrinkage and selection via the lasso, Journal of the Royal Statistical Society: Series B, 58, 267–288.
  * Tsiatis et al. (2019) Tsiatis, A. A., Davidian, M., Holloway, S. T., and Laber, E. B. (2019). Dynamic Treatment Regimes: Statistical Methods for Precision Medicine: CRC press.
  * Yuan and Lin (2006) Yuan, M. and Lin, Y. (2006). Model selection and estimation in regression with grouped variables, Journal of the Royal Statistical Society: Series B, 68, 49–67.
  * Zhang (2010) Zhang, C.-H. (2010). Nearly unbiased variable selection under minimax concave penalty, The Annals of Statistics, 38, 894–942.
  * Zhao et al. (2022) Zhao, Q., Small, D. S., and Ertefaie, A. (2022). Selective inference for effect modification via the lasso, Journal of the Royal Statistical Society: Series B, https://doi.org/10.1111/rssb.12483.
  * Zou et al. (2007) Zou, H., Hastie, T., and Tibshirani, R. (2007). On the “degrees of freedom” of the lasso, The Annals of Statistics, 35, 2173–2192.
  * 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.

##  Appendix A Appendix

###  A.1 Asymptotic expression of error for doubly robust estimator

Here, we derive the asymptotic expressions for the doubly robust estimators
used in Sections 4.2 and 4.3. The notation is the same as in those sections.
Letting $\bm{u}=(y^{[1]},y^{[2]},\ldots,y^{[H]},\allowbreak
t^{[1]},t^{[2]},\ldots,t^{[H]},\bm{x}^{[1]},\bm{x}^{[2]},\ldots,\bm{x}^{[H]},\bm{z})$,
and defining

$\displaystyle\bm{m}_{\lambda}(\bm{u};\bm{\beta},\bm{\alpha},\bm{\gamma})\equiv\sum_{h=1}^{H}\frac{\partial}{\partial\bm{\beta}}\Bigg{[}\frac{t^{[h]}}{e^{[h]}(\bm{\alpha})}\ell^{[h]}(\bm{\beta})+\Bigg{\\{}1-\frac{t^{[h]}}{e^{[h]}(\bm{\alpha})}\Bigg{\\}}{\rm
E}^{[h]}\\{\ell^{[h]}(\bm{\beta})\mid\bm{\gamma}\\}-\lambda\|\bm{\beta}\|_{1}\Bigg{]},$

the double-robust estimating equation is expressed as

$\displaystyle\sum_{i=1}^{N}\begin{pmatrix}\bm{m}_{\lambda}(\bm{u}_{i};\bm{\beta},\bm{\alpha},\bm{\gamma})\\\
\displaystyle\sum_{h=1}^{H}\frac{\partial}{\partial\bm{\alpha}}t_{i}^{[h]}\log
e_{i}^{[h]}(\bm{\alpha})\\\
\displaystyle\sum_{h=1}^{H}\frac{\partial}{\partial\bm{\gamma}}t_{i}^{[h]}\log
f_{i}^{[h]}(\bm{\gamma})\end{pmatrix}=\bm{0}_{p+q+r}.$

The limit of the solution
$(\hat{\bm{\beta}}_{\lambda},\hat{\bm{\alpha}},\hat{\bm{\gamma}})$ of this
estimating equation is denoted as
$(\bm{\beta}_{\lambda}^{*},\bm{\alpha}^{*},\bm{\gamma}^{*})$. Although
$\bm{\alpha}^{*}$ and $\bm{\gamma}^{*}$ are not necessarily true values, if
only one of them is wrong, then $\bm{\beta}_{\lambda}^{*}$ becomes the common
correct one given by (9). Through the conventional asymptotics, (23)is
obtained from lines $p+1$ to $p+q$ of this estimating equation and (24) from
lines $p+q+1$ to $p+q+r$ . Accordingly, we prepare a random function,

$\displaystyle\sum_{i=1}^{N}\sum_{h=1}^{H}\Bigg{[}$
$\displaystyle\frac{t_{i}^{[h]}}{e_{i}^{[h]}(\hat{\bm{\alpha}})}\ell_{i}^{[h]}\Bigg{(}\frac{\bm{u}^{(1)}}{N},\frac{\bm{u}^{(2)}}{\sqrt{N}}+\bm{\beta}_{\lambda}^{*(2)}\Bigg{)}-\frac{t_{i}^{[h]}}{e_{i}^{[h]}(\bm{\alpha}^{*})}\ell_{i}^{[h]}(\bm{\beta}_{\lambda}^{*(1)},\bm{\beta}_{\lambda}^{*(2)})$
$\displaystyle+\Bigg{\\{}1-\frac{t_{i}^{[h]}}{e_{i}^{[h]}(\hat{\bm{\alpha}})}\Bigg{\\}}{\rm
E}_{i}^{[h]}\Bigg{\\{}\ell_{i}^{[h]}\Bigg{(}\frac{\bm{u}^{(1)}}{N},\frac{\bm{u}^{(2)}}{\sqrt{N}}+\bm{\beta}_{\lambda}^{*(2)}\Bigg{)}\
\Bigg{|}\ \hat{\bm{\gamma}}\Bigg{\\}}$
$\displaystyle-\Bigg{\\{}1-\frac{t_{i}^{[h]}}{e_{i}^{[h]}(\bm{\alpha}^{*})}\Bigg{\\}}{\rm
E}_{i}^{[h]}\\{\ell_{i}^{[h]}(\bm{\beta}_{\lambda}^{*(1)},\bm{\beta}_{\lambda}^{*(2)})\mid\bm{\gamma}^{*}\\}$
$\displaystyle-\lambda\Bigg{\|}\frac{\bm{u}^{(1)}}{N}\Bigg{\|}_{1}-\lambda\Bigg{\|}\frac{\bm{u}^{(2)}}{\sqrt{N}}+\bm{\beta}_{\lambda}^{*(2)}\Bigg{\|}_{1}+\lambda\|\bm{\beta}_{\lambda}^{*(1)}\|_{1}+\lambda\|\bm{\beta}_{\lambda}^{*(2)}\|_{1}\Bigg{]}$
(31)

as in Knight and Fu (2000), where the asymptotic distribution of the LASSO
estimator is evaluated. Note that the minimum point of this random function
with respect to $(\bm{u}^{(1)},\bm{u}^{(2)})$ is
$(N\hat{\bm{\beta}}_{\lambda}^{(1)},\sqrt{N}(\hat{\bm{\beta}}_{\lambda}^{(2)}-\bm{\beta}_{\lambda}^{*(2)}))$.
By taking a Taylor expansion around
$(\bm{u}^{(1)},\allowbreak\bm{u}^{(2)})=(\bm{0},\bm{0})$, we can rewrite (31)
as

$\displaystyle\frac{1}{N}\sum_{i=1}^{N}\sum_{h=1}^{H}\frac{t_{i}^{[h]}}{e_{i}^{[h]}(\hat{\bm{\alpha}})}\frac{\partial}{\partial\bm{\beta}^{(1){\rm
T}}}\ell_{i}^{[h]}(\bm{\beta}_{\lambda}^{*})\bm{u}^{(1)}+\frac{1}{\sqrt{N}}\sum_{i=1}^{N}\sum_{h=1}^{H}\frac{t_{i}^{[h]}}{e_{i}^{[h]}(\hat{\bm{\alpha}})}\frac{\partial}{\partial\bm{\beta}^{(2){\rm
T}}}\ell_{i}^{[h]}(\bm{\beta}_{\lambda}^{*})\bm{u}^{(2)}$
$\displaystyle+\frac{1}{N}\sum_{i=1}^{N}\sum_{h=1}^{H}\Bigg{\\{}1-\frac{t_{i}^{[h]}}{e_{i}^{[h]}(\hat{\bm{\alpha}})}\Bigg{\\}}\frac{\partial}{\partial\bm{\beta}^{(1){\rm
T}}}{\rm
E}_{i}^{[h]}\\{\ell_{i}^{[h]}(\bm{\beta}_{\lambda}^{*})\mid\hat{\bm{\gamma}}\\}\bm{u}^{(1)}$
$\displaystyle+\frac{1}{\sqrt{N}}\sum_{i=1}^{N}\sum_{h=1}^{H}\Bigg{\\{}1-\frac{t_{i}^{[h]}}{e_{i}^{[h]}(\hat{\bm{\alpha}})}\Bigg{\\}}\frac{\partial}{\partial\bm{\beta}^{(2){\rm
T}}}{\rm
E}_{i}^{[h]}\\{\ell_{i}^{[h]}(\bm{\beta}_{\lambda}^{*})\mid\hat{\bm{\gamma}}\\}\bm{u}^{(2)}$
$\displaystyle+\frac{1}{N}\sum_{i=1}^{N}\sum_{h=1}^{H}\frac{t_{i}^{[h]}}{e_{i}^{[h]}(\hat{\bm{\alpha}})}\bm{u}^{(2){\rm
T}}\frac{\partial^{2}}{\partial\bm{\beta}^{(2)}\partial\bm{\beta}^{(2){\rm
T}}}\ell_{i}^{[h]}(\bm{\beta}_{\lambda}^{*})\bm{u}^{(2)}$
$\displaystyle+\frac{1}{N}\sum_{i=1}^{N}\sum_{h=1}^{H}\Bigg{\\{}1-\frac{t_{i}^{[h]}}{e_{i}^{[h]}(\hat{\bm{\alpha}})}\Bigg{\\}}\bm{u}^{(2){\rm
T}}\frac{\partial}{\partial\bm{\beta}^{(2)}\partial\bm{\beta}^{(2){\rm
T}}}{\rm
E}_{i}^{[h]}\\{\ell_{i}^{[h]}(\bm{\beta}_{\lambda}^{*})\mid\hat{\bm{\gamma}}\\}\bm{u}^{(2)}$
$\displaystyle-\lambda\|\bm{u}^{(1)}\|_{1}-\sqrt{N}\lambda\hskip 1.42262pt{\rm
sgn}(\bm{\beta}_{\lambda}^{*(2)})\bm{u}^{(2)}+{\rm o}_{\rm P}(1),$

and furthermore if $\hat{\bm{\alpha}}$ and $\hat{\bm{\gamma}}$ are expanded
around $\bm{\alpha}^{*}$ and $\bm{\gamma}^{*}$, we can express it as

$\displaystyle\frac{1}{N}\sum_{i=1}^{N}\sum_{h=1}^{H}\frac{t_{i}^{[h]}}{e_{i}^{[h]}(\bm{\alpha}^{*})}\frac{\partial}{\partial\bm{\beta}^{(1){\rm
T}}}\ell_{i}^{[h]}(\bm{\beta}_{\lambda}^{*})\bm{u}^{(1)}+\frac{1}{N}\sum_{i=1}^{N}\sum_{h=1}^{H}\Bigg{\\{}1-\frac{t_{i}^{[h]}}{e_{i}^{[h]}(\bm{\alpha}^{*})}\Bigg{\\}}\frac{\partial}{\partial\bm{\beta}^{(1){\rm
T}}}{\rm
E}_{i}^{[h]}\\{\ell_{i}^{[h]}(\bm{\beta}_{\lambda}^{*})\mid\bm{\gamma}^{*}\\}\bm{u}^{(1)}$
$\displaystyle+\frac{1}{\sqrt{N}}\sum_{i=1}^{N}\sum_{h=1}^{H}\Bigg{[}\frac{t_{i}^{[h]}}{e_{i}^{[h]}(\bm{\alpha}^{*})}\frac{\partial}{\partial\bm{\beta}^{(2){\rm
T}}}\ell_{i}^{[h]}(\bm{\beta}_{\lambda}^{*})\bm{u}^{(2)}+\Bigg{\\{}1-\frac{t_{i}^{[h]}}{e_{i}^{[h]}(\bm{\alpha}^{*})}\Bigg{\\}}\frac{\partial}{\partial\bm{\beta}^{(2){\rm
T}}}{\rm
E}_{i}^{[h]}\\{\ell_{i}^{[h]}(\bm{\beta}_{\lambda}^{*})\mid\bm{\gamma}^{*}\\}\bm{u}^{(2)}\Bigg{]}$
$\displaystyle+\frac{1}{\sqrt{N}}\sum_{i=1}^{N}\sum_{h=1}^{H}(\hat{\bm{\alpha}}-\bm{\alpha}^{*})^{{\rm
T}}\frac{\partial^{2}}{\partial\bm{\alpha}\partial\bm{\beta}^{(2){\rm
T}}}\frac{t_{i}^{[h]}\ell_{i}^{[h]}(\bm{\beta}_{\lambda}^{*})}{e_{i}^{[h]}(\bm{\alpha}^{*})}\bm{u}^{(2)}$
$\displaystyle+\frac{1}{\sqrt{N}}\sum_{i=1}^{N}\sum_{h=1}^{H}(\hat{\bm{\alpha}}-\bm{\alpha}^{*})^{{\rm
T}}\frac{\partial^{2}}{\partial\bm{\alpha}\partial\bm{\beta}^{(2){\rm
T}}}\Bigg{\\{}1-\frac{t_{i}^{[h]}}{e_{i}^{[h]}(\bm{\alpha}^{*})}\Bigg{\\}}{\rm
E}_{i}^{[h]}\\{\ell_{i}^{[h]}(\bm{\beta}_{\lambda}^{*})\mid\bm{\gamma}^{*}\\}\bm{u}^{(2)}$
$\displaystyle+\frac{1}{\sqrt{N}}\sum_{i=1}^{N}\sum_{h=1}^{H}\Bigg{\\{}1-\frac{t_{i}^{[h]}}{e_{i}^{[h]}(\bm{\alpha}^{*})}\Bigg{\\}}(\hat{\bm{\gamma}}-\bm{\gamma}^{*})^{{\rm
T}}\frac{\partial^{2}}{\partial\bm{\gamma}\partial\bm{\beta}^{(2){\rm T}}}{\rm
E}_{i}^{[h]}\\{\ell_{i}^{[h]}(\bm{\beta}_{\lambda}^{*})\mid\bm{\gamma}^{*}\\}\bm{u}^{(2)}$
$\displaystyle+\frac{1}{N}\sum_{i=1}^{N}\sum_{h=1}^{H}\bm{u}^{(2){\rm
T}}\frac{\partial^{2}}{\partial\bm{\beta}^{(2)}\partial\bm{\beta}^{(2){\rm
T}}}\Bigg{[}\frac{t_{i}^{[h]}\ell_{i}^{[h]}(\bm{\beta}_{\lambda}^{*})}{e_{i}^{[h]}(\bm{\alpha}^{*})}+\Bigg{\\{}1-\frac{t_{i}^{[h]}}{e_{i}^{[h]}(\bm{\alpha}^{*})}\Bigg{\\}}{\rm
E}_{i}^{[h]}\\{\ell_{i}^{[h]}(\bm{\beta}_{\lambda}^{*})\mid\bm{\gamma}^{*}\\}\Bigg{]}\bm{u}^{(2)}$
$\displaystyle-\lambda\|\bm{u}^{(1)}\|_{1}-\sqrt{N}\lambda\hskip 1.42262pt{\rm
sgn}(\bm{\beta}_{\lambda}^{*(2)})\bm{u}^{(2)}+{\rm o}_{\rm P}(1).$

Moreover, using the law of large numbers, we find that

$\displaystyle{\rm
E}\Bigg{(}\sum_{h=1}^{H}\Bigg{[}\frac{t^{[h]}}{e^{[h]}(\bm{\alpha}^{*})}\frac{\partial}{\partial\bm{\beta}^{(1){\rm
T}}}\ell^{[h]}(\bm{\beta}_{\lambda}^{*})+\Bigg{\\{}1-\frac{t^{[h]}}{e^{[h]}(\bm{\alpha}^{*})}\Bigg{\\}}\frac{\partial}{\partial\bm{\beta}^{(1){\rm
T}}}{\rm
E}^{[h]}\\{\ell^{[h]}(\bm{\beta}_{\lambda}^{*})\mid\bm{\gamma}^{*}\\}\Bigg{]}\Bigg{)}\bm{u}^{(1)}-\lambda\|\bm{u}^{(1)}\|_{1}$
$\displaystyle+\frac{1}{\sqrt{N}}\sum_{i=1}^{N}\sum_{h=1}^{H}\Bigg{[}\frac{t_{i}^{[h]}}{e_{i}^{[h]}(\bm{\alpha}^{*})}\frac{\partial}{\partial\bm{\beta}^{(2){\rm
T}}}\ell_{i}^{[h]}(\bm{\beta}_{\lambda}^{*})\bm{u}^{(2)}+\Bigg{\\{}1-\frac{t_{i}^{[h]}}{e_{i}^{[h]}(\bm{\alpha}^{*})}\Bigg{\\}}\frac{\partial}{\partial\bm{\beta}^{(2){\rm
T}}}{\rm
E}_{i}^{[h]}\\{\ell_{i}^{[h]}(\bm{\beta}_{\lambda}^{*})\mid\bm{\gamma}^{*}\\}\bm{u}^{(2)}\Bigg{]}$
$\displaystyle+\sqrt{N}(\hat{\bm{\alpha}}-\bm{\alpha}^{*})^{{\rm T}}{\rm
E}\Bigg{\\{}\sum_{h=1}^{H}\frac{\partial^{2}}{\partial\bm{\alpha}\partial\bm{\beta}^{(2){\rm
T}}}\frac{t^{[h]}\ell^{[h]}(\bm{\beta}_{\lambda}^{*})}{e^{[h]}(\bm{\alpha}^{*})}\Bigg{\\}}\bm{u}^{(2)}$
$\displaystyle-\sqrt{N}(\hat{\bm{\alpha}}-\bm{\alpha}^{*})^{{\rm T}}{\rm
E}\Bigg{\\{}\sum_{h=1}^{H}\frac{\partial^{2}}{\partial\bm{\alpha}\partial\bm{\beta}^{(2){\rm
T}}}\frac{t^{[h]}}{e^{[h]}(\bm{\alpha}^{*})}{\rm
E}^{[h]}\\{\ell^{[h]}(\bm{\beta}_{\lambda}^{*})\mid\bm{\gamma}^{*}\\}\Bigg{\\}}\bm{u}^{(2)}$
$\displaystyle+\sqrt{N}(\hat{\bm{\gamma}}-\bm{\gamma}^{*})^{{\rm T}}{\rm
E}\Bigg{[}\sum_{h=1}^{H}\Bigg{\\{}1-\frac{t^{[h]}}{e^{[h]}(\bm{\alpha}^{*})}\Bigg{\\}}\frac{\partial^{2}}{\partial\bm{\gamma}\partial\bm{\beta}^{(2){\rm
T}}}{\rm
E}^{[h]}\\{\ell^{[h]}(\bm{\beta}_{\lambda}^{*})\mid\bm{\gamma}^{*}\\}\Bigg{]}\bm{u}^{(2)}$
$\displaystyle+\bm{u}^{(2){\rm T}}{\rm
E}\Bigg{(}\sum_{h=1}^{H}\frac{\partial^{2}}{\partial\bm{\beta}^{(2)}\partial\bm{\beta}^{(2){\rm
T}}}\Bigg{[}\frac{t^{[h]}\ell^{[h]}(\bm{\beta}_{\lambda}^{*})}{e^{[h]}(\bm{\alpha}^{*})}+\Bigg{\\{}1-\frac{t^{[h]}}{e^{[h]}(\bm{\alpha}^{*})}\Bigg{\\}}{\rm
E}^{[h]}\\{\ell^{[h]}(\bm{\beta}_{\lambda}^{*})\mid\bm{\gamma}^{*}\\}\Bigg{]}\Bigg{)}\bm{u}^{(2)}$
$\displaystyle-\sqrt{N}\lambda\hskip 1.42262pt{\rm
sgn}(\bm{\beta}_{\lambda}^{*(2)})\bm{u}^{(2)}+{\rm o}_{\rm P}(1).$