Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
47,503 characters · 15 sections · 39 citation commands
Average Adjusted Association: Efficient Estimation with High Dimensional Confounders
\doublespacing
There are several statistical measures of association, among which the (log) odds ratio is one of the most popular ones. The odds ratio has been frequently used in medicine, biostatistics and epidemiology because it is simple, it has a natural interpretation in the standard logistic regression model, and it is invariant under various sampling designs that include case-control studies. See breslow1976regression,breslow1978there,breslow1980statistical for early work and Bland1468,OR-JAMA for how to use the odds ratio in practice.
To describe the object of interest and the background more concretely, let $Y$ denote a binary outcome, $T$ a binary exposure and $X$ a vector of measured confounders/covariates. The conditional odds ratio between $Y$ and $T$ given $X=x$, which we denote by $\mathrm{OR}(x)$, provides a complete picture of the adjusted association between $Y$ and $T$ for different subpopulations defined by different values of $x$. Conditioning on $X=x$ matters for two reasons: first, it is more plausible to make a causal interpretation out of $\mathrm{OR}(x)$ than the unconditional odds ratio between $Y$ and $T$ Holland:Rubin,greenland1999confounding; second, the association can be heterogeneous over different subpopulations so that $\mathrm{OR}(x)$ is a complicated function of $x$ in general.
There are several approaches to modeling $\mathrm{OR}(\cdot)$. It is most common to parametrize the function $\mathrm{OR}(\cdot)$ typically via a logistic regression model, but it may suffer from misspecification. Therefore, effort has been made to mitigate the problem from the perspective of doubly robust estimation yun2007semiparametric,tchetgen2010doubly,tchetgen2013closed. On the contrary, some authors have emphasized heterogeneity, advocating nonparametric estimation of the function $\mathrm{OR}(\cdot)$ chen2011nonparametric,hui2013nonparametric. However, fully nonparametric approaches generally suffer from the curse of dimensionality and they are often not a practical option in finite samples. In fact, all the numerical experiments in chen2011nonparametric and hui2013nonparametric are limited with the scalar $X$.
In this paper, we consider the case where the dimension of $X$ is high. We introduce a new summary measure of association, which we call the average adjusted association, by taking the average of conditional log odds ratios over covariates. That is, we think of $\theta_0 := \mathbb{E}\{ \log\mathrm{OR}(X) \}$ as a summary measure of adjusted association {in a heterogeneous population}. The issues of what to condition on for $X$ and how to summarize $\mathrm{OR}(\cdot)$ have been recognized in the literature muller2014estimating,greenland1999confounding. However, to the best of our knowledge, {the average adjusted association has not been formally studied in the literature,} let alone how to deal with high dimensional $X$. {Indeed, many recently published papers in both natural and social sciences have used the odds ratio to report their findings but, to the best of our knowledge, they are all based on strong parametric assumptions that restrict heterogeneous association: a logistic regression model is the most popular, where heterogeneity in association needs to be pre-specified. For example, recent studies on determinants of callbacks to job applications farber2016determinants, association of Alzheimer's dementia with genotypes reiman2020exceptionally, and police use of force with respect to race hoekstra2022does are all based on logistic regression and pre-specified odds ratios. It is worth noting that odds ratios are natural in these examples because callbacks, dementia, and use of force are all rare events that occur with small probabilities.}
{As we depart from parametric models, the odds ratio can be a complicated function of observed characteristics.} It is a natural and common practice to summarize a heterogeneous quantity by taking an average. For instance, when the treatment effect is heterogeneous, we frequently focus on the average treatment effect HIR:2003,ray2019debiased,Shi:Blei:Veitch {as a summary parameter}. Despite the popularity and usefulness of the average treatment effect and related quantities, the literature on such a statistic based on the odds ratios is surprisingly sparse. {The only exception we are aware of is the Mantel-Haenszel approach that averages the (log) odds ratio across finite strata.} This paper fills this important gap.
Treating the function $\log\{\mathrm{OR}(\cdot)\}$ nonparametrically, we derive two equivalent forms of the efficient influence function for $\theta_0$, which we then explore to construct efficient double/debiased machine learning (DML) estimators of $\theta_0$. We take this approach because the class of DML estimators are generally more advocated than conventional plug-in efficient estimators particularly when $X$ is high-dimensional chernozhukov2018double, lewis2021double. The efficient influence function can be expressed in terms of either prospective probabilities (i.e., the probabilities of the outcome conditional on the exposure and confounders) or retrospective ones (i.e., the probabilities of the exposure given the outcome and confounders). Therefore, we have two types of DML estimators: they are asymptotically equivalent and efficient, but they are not the same in finite samples. Our work is the first to derive the forms of the efficient influence function and to propose suitable DML estimators of $\theta_0$. We also provide easy-to-follow computational and inferential algorithms for implementation. {Given the basic nature of this research, we do not see potential negative societal impacts of our work, although we should mention that it generally requires much caution to draw causal inference from statistical association.}
\paragraph{Related Literature}
There are a couple of recent papers on automatic DML estimators (see, e.g. pmlr-v162-chernozhukov22a and CNS:2022). However, applying automatic orthogonalization to our setting is not necessarily better because (i) automation may induce additional approximation errors and (ii) the explicit formula obtained in this paper can provide useful insight into the estimation problem.
General theory for semiparametric efficiency is well-developed in the literature Ai:Chen:2012,ACHL:2014. However, if we used the general framework, we would need to verify all the regularity conditions and our proof would not be self-contained. This type of verification is not needed for our setting because we can directly calculate the efficient influence function for $\theta_0$ with a more preliminary proof technique.
\paragraph{Replication Files}
The replication files for all the numerical results are available at \url{https://github.com/sokbae/replication-JunLee-2023-AISTATS}.
The odds ratio can be expressed by using either prospective or retrospective probabilities: i.e., for all $x$ in the support $\mathcal{X}$ of $X$,
where the second equality is known as the invariance property of the odds ratio, which can be verified by the Bayes rule cornfield1951method. The two expressions of $\mathrm{OR}$ can be used to develop two different machine learning estimators.
Since the function $\mathrm{OR}$ is an infinite-dimensional object, it is generally difficult to estimate with high precision or even to communicate estimation results in a fully nonparametric manner, unless the dimension of $x$ is limited to 1 or 2. In this context, we propose, as a scalar summary measure of association, to take the expectation of $\log\mathrm{OR}(\cdot)$ using the probability distribution of $X$. Specifically, we define
which can be understood as the Average Adjusted Association (AAA) for the entire population. We have taken the logarithm before taking expectation because $\log\mathrm{OR}(x)$ corresponds to a coefficient in the traditional logistic model; for instance, if $\mathbb{P}(Y=1 \mid T=t,X=x) = G(\alpha_0 + \alpha_1 t + \alpha_2^{\mathpalette\raiseT{\intercal}} x)$, where $G(s) = \exp(s)/\{1+\exp(s) \}$, then $\log\mathrm{OR}(x) = \alpha_1$ for all $x\in \mathcal{X}$. But the essence of our results does not rely on this structure; all of our results can be straightforwardly modified to the case of aggregating without the logarithm. Since $\mathrm{OR}(\cdot)$ is a nonlinear function in general, $\theta_0$ differs from the log odds evaluated at the average value of $X$.
We now characterize the efficient influence function for estimating $\theta_0$ and present two equivalent expressions. We first state {an assumption} that will be used throughout the paper. Recall that $Y$ and $T$ are binary variables, and let $\mathcal{X}$ be the support of $X$.
{(ref) is to ensure that all the four joint outcomes of $(Y,T)$ occur with positive probability conditional on any value of $X=x\in\mathcal{X}$. Therefore, conditioning on any outcome of $(Y,T)$ does not exclude any value of $X$ from the support: i.e., (ref) implies that the joint support of $(Y,T,X)$ is given by $\{0,1\}\times\{0,1\}\times \mathcal{X}$. This requirement may not be trivial in some applications, unless we restrict out attention to a certain subpopulation. For example, if $Y$ represents prostate cancer, then it is reasonable to focus on the subpopulation of men. We are implicit about this type of (extra) conditioning throughout the analysis.}
{Also, under (ref), both of the prospective and retrospective representations of $\mathrm{OR}$ are well-defined for any $x\in \mathcal{X}$. Further, the existence of such an $\epsilon>0$ in (ref) guarantees that $x \mapsto \mathrm{OR}(x)$ is uniformly bounded from below by zero and from above by infinity.}
{Below we discuss the efficient influence function for $\theta_0$ when (ref) is the only restriction imposed on the distribution of $(Y,T,X)$.}
For $y, t \in \{0,1\}$, define
Further, define
The efficient influence function for $\theta_0$ is given in the following theorem.
(ref) includes equality between $F_p(Y,T,X)$ and $F_r(Y,T,X)$; $F_p$ is based on the prospective expression of $\mathrm{OR}$, whereas $F_r$ uses the retrospective one. Theorem (ref) is proved in two steps: (i) by direct calculation, as in hahn1998role, $\theta_0$ is pathwise differentiable along regular parametric submodels in the sense of newey1990semiparametric,newey1994asymptotic; (ii) the pathwise derivative is an element of the tangent space, from which we can obtain the semiparametric efficiency bound $V_{\textrm{eff}}$ for $\theta_0$:
The bound $V_{\textrm{eff}}$ can be achieved by double/debiased machine learning (DML) estimators based on the efficient influence function, i.e., the representation in either (ref) or (ref). The DML approach has an advantage that it is robust to local perturbation on the unknown functions that need to be estimated in the first step, which is known as the Neyman orthogonality property chernozhukov2018double.
In the remaining part of this section we formally show that the moment condition based on either (ref) or (ref) is indeed robust to local perturbation on the nonparametric components. For this purpose, note first that each of $F_p$ and $F_r$ depends on three nonparametric elements, i.e.,
respectively. Let $\mathcal{G}$ be a space of (measurable) functions on $\mathcal{X}$ such that $g\in \mathcal{G}$ satisfies $0<\inf g(x) \leq \sup g(x) < 1$: see (ref). Let $\tilde F_p(\cdot)[Y,T,X]$ and $\tilde F_r(\cdot)[Y,T,X]$ denote the functionals defined on $\mathcal{G}^3$ such that $\tilde F_p(\eta_{p0})[Y,T,X] = F_p(Y,T,X)$ and $\tilde F_r(\eta_{r0})[Y,T,X] = F_r(Y,T,X)$. Then, Neyman orthogonality is concerned about the Gateaux derivatives of $\tilde F_p$ and $\tilde F_r$ at $\eta_{p0}$ and $\eta_{r0}$, respectively.
(ref) says that both $F_p(Y,T,X)$ and $F_r(Y,T,X)$ provide Neyman orthogonal moments conditional on $X$. One implication of the local robustness property is that the first step nonparametric estimation of $\eta_{p0}$ (or $\eta_{r0}$) in estimating $\theta_0$ will not have any first-order consequence, i.e., the limiting distribution would be the same as if $\eta_{p0}$ (respectively, $\eta_{r0}$) were known. In other words, all the adjustment terms that are needed to address the effect of the first step estimation are already reflected in $F_p$ (respectively, in $F_r$).
We now describe a couple of double/debiased machine learning (DML) estimators of $\theta_0$, for which we use the functions $F_p$ and $F_r$ as estimating equations. We assume that a random sample $\{ (Y_i, T_i, X_i^{\mathpalette\raiseT{\intercal}})^{\mathpalette\raiseT{\intercal}}:\ i=1,2,\dots, n\}$ is available, where $X_i$ is allowed to be high dimensional.
Let $K \geq 2$ be some fixed integer (say, 5, 10 or 20). For simplicity, assume that $n$ is divisible by $K$. Let $\{ I_k : k=1,\ldots,K \}$ denote a $K$-fold partition of $\{1,\ldots,n\}$ such that $| I_k | = n/K$ for each $k$. Suppose that one estimates all the conditional probabilities appearing in (ref) or (ref), depending on which one to use for estimation, via machine-learning estimators by using observations that belong to $I_k^c = \{1,\ldots,n\} \setminus I_k$ for each $k$. Using data in $I_k^c$ to estimate conditional probabilities evaluated at the points in $I_k$ is reminiscent of the traditional leave-one-out method.
In using the prospective formula in (ref), we start with
where $\widehat{\mathbb{P}}_{\mathrm{ML},k}$ denotes a machine-learning estimator of a probability model using observations that belong to $I_k^c$. We then define the prospective DML estimator $\widehat{\theta}_{p} $ of $\theta_0$ by
where
The estimator $\hat \theta_p$ is asymptotically normal and efficient as we will show in (ref). We have summarized the estimation procedure in (ref).
The retrospective DML estimator $\widehat{\theta}_{r} $ of $\theta_0$ is defined analogously. That is, we start with
and we define
where
The algorithm for the retrospective DML estimator can be stated easily by making simple modifications in (ref). We omit details for brevity.
Before we finish this subsection, we remark that the consistency of $\widehat{\theta}_{p}$ and $\widehat{\theta}_{r}$ does not require that $\widehat w_{p,k}$ and $\widehat w_{r,k}$ consistently estimate $w_p$ and $w_r$; indeed, $F_p(Y,T,X)$ and $F_r(Y,T,X)$ have mean zero even if $w_p$ and $w_r$ deviate from the truth.
Let $\| \cdot \|_{P, 2}$ denote the $L_2(P)$-norm, where $P$ is the probability distribution of $(Y,T,X)$: i.e., $\| a \|_{P, 2} := \max_{1 \leq \ell \leq d} \left\{ \mathbb{E} [ a_\ell^2 (Y,T,X) ] \right\}^{1/2}$ for a $d$-dimensional vector-valued function $a := (a_1,\ldots,a_d)$. For each $k$, let
denote the vector of machine learning estimators of $\eta_{p0}$ defined in (ref), using observations belonging to $I_k^c$. The first step estimators $\widehat\eta_{n,p,k}$ are inputs to the second step in (ref). Likewise, let $\widehat\eta_{n,r,k}$ be the vector of machine learning estimators of $\theta_{r0}$ defined in (ref), which will be used for the retrospective DML estimator of $\theta_0$.
(ref) is a high-level assumption that may not be trivial if $X$ is high dimensional. For instance, it may fail even when all the conditional probabilities are logistically specified and they are estimated by the method of maximum likelihood if $X$ is high dimensional sur2019modern,zhao2022asymptotic. However, (ref) is known to be attainable for a variety of machine learning methods. The primitive conditions for $\ell_1$-penalized logit estimators are worked out by vandegeer2008 and Belloni:2016:JBES among others. A maximum likelihood approach with some adjustment for the dimensionality and signal strength of $X$ as in e.g., yadlowsky2021sloe is another possibility although we do not pursue the latter in this paper.
An application of Theorems 3.1 and 3.2 of chernozhukov2018double gives the following result that formally justifies the estimation and inference methods proposed in (ref).
(ref) establishes that both the prospective and retrospective DML estimators are asymptotically normal and efficient. However, asymptotic equivalence does not imply that it is irrelevant which one to use between $\widehat{\theta}_{p}$ and $\widehat{\theta}_{r}$ in finite samples. In fact, the first steps of the two estimators involve different nonparametric regression functions, and therefore they generally lead to different estimates. Comparing the estimates and their standard errors can be a useful diagnostic check in practice; we recommend reporting inference results based on both estimators.
Since we do not use any parametric specification to identify and estimate $\theta_0$, potential misspecification is not a concern, at least asymptotically{: it is a concern only to the extent that our choices of nonparametric estimators must satisfy (ref).} However, nonparametric approaches do rely on several input parameters, which are important for the performance of the estimators in finite samples. In this regard we show that there is a nonparametric version of double robustness.
For every $x\in \mathcal{X}$, we implicitly define $\varphi_{p0}(x), \varphi_{r0}(x)$, and $\vartheta_0(x)$ by the following equations: for $t,y\in\{0,1\}$,
For instance, (ref) with $t=0$ defines $\varphi_{p0}(x)$. Here, we have four conditional probabilities to define three objects. This is because the four conditional probabilities are restricted by the invariance property of the odds ratio. Indeed, simple algebra shows that $\vartheta_0(x) = \log \mathrm{OR}(x)$, on which no restrictions have been imposed.
Let $\mathcal{H}$ be the class of functions on $\mathcal{X}$, which contains $\varphi_{p0}, \varphi_{r0}$, and $\vartheta_0$. We then consider the function $m:\mathcal{H}^3\times \{0,1\}^2 \times \mathcal{X} \rightarrow \mathbb{R}$ defined by
where $\Lambda_0(\varphi, x) = \exp\{\varphi(x)\} / \bigl[ 1+\exp\{\varphi(x)\} \bigr]$.
(ref) is a nonparametric extension of the double robustness idea of tchetgen2010doubly and tchetgen2013closed: i.e., either $\mathbb{P}(Y=1\mid T=0,X=x ) = \varphi_{p0}(x)$ or $\mathbb{P}(T=1\mid Y=0,X=x ) = \varphi_{r0}(x)$ (but not both) needs to be correctly specified to estimate $\vartheta_0(x) = \log\mathrm{OR}(x)$ consistently. Further, the second assertion shows that the function $m$ can be used to identify $\vartheta_0$ once $\varphi_{p0}$ and $\varphi_{r0}$ are given by their definition. It can be shown that the efficiency bound for estimating $\theta_0 = \mathbb{E}\{ \vartheta_0(X) \}$ by using the conditional moment equation in (ref) is given by
which coincides with $F_p(Y,T,X) = F_r(Y,T,X)$. Therefore, we do not lose anything in terms of semiparametric efficiency in using the moment function $m$ to estimate $\theta_0$.
Therefore, it is possible to have an extra DML-based doubly robust algorithm for efficient estimation of $\theta_0$. However, we do not pursue this possibility in the current paper for a couple of reasons. As we described in (ref), the DML approach requires splitting the sample into multiple subsamples, but it often causes computational issues in practice when multiple tiers of nonparametric estimation are involved as in the current setup. Also, since our approach does not require any parametric specification at all, double robustness seems to have limited merits; there is no misspecification in the limit in the approach described in (ref).
One can rely on (ref) to conduct statistical inference on $\theta_0$. For instance, using the prospective estimator $\hat \theta_p$, a (symmetric two-sided) $95\%$ confidence interval for $\theta_0$ can be obtained in the usual manner, i.e., $\hat\theta_p \pm 1.96\cdot \hat\sigma_p/\sqrt{n}$. Here, we emphasize that $\theta_0$ is understood as a simple association parameter to which we do not give any causal interpretation at this stage. However, $\theta_0$ can be used for causal inference under a few extra assumptions.
In order to discuss causal inference, let $Y(t)$ denote the potential outcome when the treatment is exogenously fixed at $t$; i.e., the observed outcome $Y$ is equal to $Y(T) = Y(1)T + Y(0)(1-T)$. Then, $\mathrm{OR}(x)$ is related with the following causal parameters:
That is, $\vartheta_\mathrm{OR}(x)$ represents the adjusted causal odds ratio, while $\vartheta_\mathrm{RR}(x)$ is the adjusted causal relative risk parameter. Below we summarize some facts that are known in the literature; see, e.g., Holland:Rubin and jun2021causal for more detail.
Therefore, the average parameter $\theta_0$ and (ref) can be used for causal inference on the entire population. For example, if the researcher is willing to assume that the treatment was randomly assigned conditional on $X$, then the usual symmetric confidence interval such as $\hat\theta_r\pm 1.96 \cdot\hat\sigma_r$ is a confidence interval for $\mathbb{E}\{ \vartheta_\mathrm{OR}(X)\}$. There is no general relationship between $\mathbb{E}\{ \vartheta_\mathrm{OR}(X)\}$ and $\mathbb{E}\{ \vartheta_\mathrm{RR}(X)\}$, but if $\mathbb{P}(Y=1)$ is close to zero (known as the rare-disease assumption), then the two causal parameters are known to be close to each other as long as there is no confounding given $X$. So, the symmetric confidence interval of $\theta_0$ can also be understood as an approximate confidence interval of $\mathbb{E}\{ \vartheta_\mathrm{RR}(X) \}$ in this scenario.
The no-confounding assumption or the rare-disease assumption can be unrealistic in some applications. For instance, many treatments of interest in social sciences such as education choices are deliberate decisions, and in such cases the no-confounding assumption is unrealistic. Even so, if one is willing to assume that education is potentially beneficial but it never hurts and that those who deliberately chose to take higher education is generally no less likely to “succeed” than those who did not, then $\theta_0$ can be understood as a sharp upper bound on $\mathbb{E}\{\vartheta_\mathrm{RR}(X)\}$. Therefore, a one-sided confidence interval such as $[1,\ \hat\theta_r + 1.64\cdot \hat \sigma_r]$ can be reported if the causal parameter $\mathbb{E}\{ \vartheta_\mathrm{RR}(X)\}$ is of interest.
In practice there may be a subpopulation of particular interest, in which case averaging over the entire population may not provide the most relevant summary statistic. For example, consider the population of patients with a certain type of cancer. Suppose that $Y$ and $T$, respectively, indicate five-year survival (say, $Y=1$ for survival and $Y=0$ for death) and a certain type of treatment (say, $T=1$ for treatment and $T=0$ for no treatment). Here, it may be relevant to summarize the association between $Y$ and $T$ for those who received the treatment, i.e., $\theta_T(1) := \mathbb{E}\{\log\mathrm{OR}(X) \mid T=1 \}$. Alternatively, the association between $Y$ and $T$ for those who survived the cancer, which can be captured by $\theta_Y(1) := \mathbb{E}\{ \log\mathrm{OR}(X) \mid Y=1 \}$, can be an interesting quantity to look at. Of course, if the average adjusted association between $T$ and $Y$ is homogeneous across $X$, then there will be no difference among $\theta_0, \theta_T(1)$, and $\theta_Y(1)$. However, in general, they are all distinct and can be substantially different, depending on the degree of heterogeneity.
Also, $\theta_T(1)$ and $\theta_Y(1)$ can be of interest if our access to a random sample is limited. For example, $\theta_T(1)$ is point identifiable even when we only have a treatment-based sample, where a half of the sample comes from the patients who received the treatment and the other half is from those who did not. Similarly, an outcome-based sample (e.g., case-control studies) is sufficient to identify $\theta_Y(1)$. The statistical analysis of $\theta_T(1)$ and $\theta_Y(1)$ is similar to that of $\theta_0$ and we do not repeat it here.
In addition, there has been increasing interest in estimating individual level treatment effects pmlr-v70-shalit17a. {Averaging over the entire population may have a risk of over-simplification: for example, if the association is positive for some values of $X$, and it is negative for other values of $X$, then the overall average association may be close to zero.} With this motivation in mind, suppose that there is a vector of low-dimensional confounders (say, $Z$) such that we are interested in estimating $z \mapsto \mathbb{E}\{\log\mathrm{OR}(X) \mid Z=z \}$. Then it can be estimated by projecting $\log\mathrm{OR}(X)$ on a low-dimensional space of $Z$ as in Ogburn:2015, Lee:Okui:Whang, and Semenov:Chernozhukov:20. However, it is a topic of future research to develop this idea formally.
In this section, we provide numerical results to illustrate the usefulness of our approach.
We start with a real-data example. Table (ref) summarizes data from American Community Survey (ACS) 2018, cross-tabulating the likelihood of top income by educational attainment. The sample is restricted to white males residing in California with at least a bachelor's degree. It is extracted from IPUMS USA IPUMS. The ACS is an ongoing annual survey by the US Census Bureau that provides key information about the US population.
The binary outcome variable `Top income' ($Y$) is defined to be one if a respondent's annual total pre-tax wage and salary income is top-coded. In ACS 2018, the threshold income for top-coding is different across states. In our sample extract, the top-coded income bracket has median income \$565,000 and the next highest income that is not top-coded is \$327,000. The binary exposure variable ($T$) is defined to be one if a respondent has a master's degree, a professional degree, or a doctoral degree.
To adjust for individual differences, we include age and industry code as covariates ($X$). In particular, cubic B-splines of age with 17 inner knots as well as 254 industry dummies are included in this specification, which can be viewed as a high-dimensional setting.
Specifically, we implement $\ell_1$-penalized logistic estimation with glmnet package in R glmnet to estimate $\mathbb{P}(Y=1|T=t,X=x), t=0,1$ and $\mathbb{P}(T=1|X=x)$ for the prospective model (respectively, $\mathbb{P}(T=1|Y=y,X=x), y=0,1$ and $\mathbb{P}(Y=1|X=x)$ for the retrospective model) with 10-fold cross-fitting. The underlying assumption here is that the B-spline terms plus the industry dummies are rich enough to approximate $\mathbb{P}(Y=1|T=t,X=x)$ as well as $\mathbb{P}(T=1|X=x)$ for the prospective model (respectively, $\mathbb{P}(T=1|Y=y,X=x)$ as well as $\mathbb{P}(Y=1|X=x)$ for the retrospective model). The penalization tuning parameter is chosen by cross-validation (that is, lambda.min in the glmnet package). Here, we focus on $\ell_1$-penalized estimators among other possible machine learning estimators because the primitive conditions for Assumption (ref) are well established for $\ell_1$-penalized logit estimators vandegeer2008,Belloni:2016:JBES, as we mentioned in Section (ref).
(ref) reports estimation results. Looking at Panel A, the prospective estimate of $\theta_0$ is 0.72, which is almost the same as the retrospective estimate of 0.71. In Panel B, we present point estimates of $\exp(\theta_0)$ and its confidence intervals using asymptotic normality obtained in (ref).
The estimates of $\exp(\theta_0)$ are comparable to the usual odds ratio in terms of its scale; therefore, they can be interpreted similarly. Furthermore, as can be seen from (ref), top income is a rare event (that is, the sample proportion of $\mathbb{P}(Y=1)$ is approximately $0.05$). When the outcome of interest is rare, an average of the conditional odds ratio approximates the average of conditional relative risk, namely the average of the ratio between $\mathbb{P}(Y=1|T=1,X)/\mathbb{P}(Y=1|T=0,X)$. Hence, obtaining a higher-level degree is associated with doubling the chance of earning very high incomes. The 95% confidence interval for the prospective estimate is $[1.61,2.63]$ (respectively, [1.67,2.49] for the retrospective model). As discussed in Section (ref), $\theta_0$ can be interpreted as the upper bound on the causal parameter $\mathbb{E}\{ \log \vartheta_\mathrm{RR}(X) \}$ if one assumes the MTR/MTS assumptions here.
Recall that the covariates consists of cubic B-splines of age as well as industry dummies, resulting in 274 regressors. As $n = 17,816$, one may simply try to estimate a parametric logistic regression model with the same set of regressors. However, it turns out that this flexible parametric approach suffers from a couple of numerical issues: (i) a very small number of estimated coefficients are NA due to multicollinearity; (ii) some of predicted probabilities are numerically zero. As a result, the conditional odds ratios are not defined for all values of the regressors. To resolve these problems, we make some ad hoc adjustments: (i) we ignore the problematic regressors by setting their coefficients to be zero; (ii) we set a lower bound on the fitted probabilities. Specifically, any fitted probability less than 1e-6 is set to be 1e-6. Then, we estimate $\theta_0$ by simply plugging the fitted probabilities into the formula of $\theta_0 = \mathbb{E}\{ \log\mathrm{OR}(X) \}$. The parametric plug-in estimates with the ad hoc adjustments turn out to be 0.38 (prospective estimate) and 0.88 (retrospective estimate). The large difference between the prospective and the retrospective estimates indicates that there is an anomaly in the plug-in estimates. In addition, we also consider plug-in estimation of $\theta_0$ using $\ell_1$-penalized logistic estimation with the same specifications and tuning parameters as in DML estimation. Hence, in this case, the plug-in estimator is different from the DML estimator in that (i) it uses a different estimating equation and (ii) it does not use cross-fitting. The resulting plug-in estimates are 0.78 (prospective estimate) and 0.68 (retrospective estimate). They look more similar to the DML estimators; however, there is no theoretically proven result regarding how to conduct inference with the $\ell_1$-penalized plug-in estimators.
We turn to a Monte Carlo experiment to make a more systematic comparison between the plug-in and DML estimators. We generate observations in the following way: (i) the covariates are randomly drawn from the empirical distribution of the ACS sample; (ii) the binary exposure variable is generated from a logit model with $\mathbb{P}(T=1|X) = G( \alpha_0 + \alpha_1 Age + \alpha_2 Age^2)$, where $G(\cdot)$ is the logit link function and the parameters $(\alpha_0, \alpha_1, \alpha_2)$ are chosen by fitting the logit model with the ACS sample; (iii) the binary outcome variable is drawn from a logit model with $\mathbb{P}(Y=1|T,X) = G( \beta_0 + \beta_1 T + \beta_2 Age + \beta_3 Age^2)$, where the parameters $(\beta_0, \beta_1, \beta_2, \beta_3)$ are again chosen by fitting a logit model with the ACS sample. In this experimental design, the true model is such that $\theta_0 = \beta_1$ and the industry effects are null. However, we fit exactly the same specifications as in the previous real-data example to examine the differences between the plug-in and DML estimators. The only change made here is that 5-fold cross-validation is adopted for DML estimation to speed up Monte Carlo simulations. The sample size is $n=5,000$ and the number of Monte Carlo replications is 500.
Table (ref) summarizes the results of the experiments. The prospective DML estimator has a much smaller mean bias than the prospective plug-in estimator without increasing the standard deviation. Further, its size (the probability of excluding the true value of $\theta_0$ in the confidence interval) is close to the 10% nominal level. The retrospective DML estimator does not perform as well as the prospective DML estimator. This is due to the fact that experimental data are generated from a prospective logit model. Overall, the results of the experiment verify that the DML estimators are superior to the plug-in estimators when the underlying machine learning estimators are the $\ell_1$-penalized logistic regression estimators.
Our proposed DML estimators offer a novel way of estimating the summary measure of association, namely the AAA functional $\theta_0$. In particular, we provide a method for statistical inference on $\theta_0$ based on asymptotic normality of our efficient DML estimators.
This paper has focused on binary outcome and exposure. However, it is possible to define the AAA functional beyond the current setup. For example, following tchetgen2010doubly, we can define the conditional odds ratio function as \[ \mathrm{OR}(x):=\frac{f(y \mid t, x)}{f(y \mid t_0, x)} \frac{f\left(y_0 \mid t_0, x\right)}{f\left(y_0 \mid t, x\right)}, \] where $Y$ and $T$ can take either discrete values, continuous values, or a mixture of both; $\left(y_0, t_0\right)$ is a user specified point in the sample space; and $f(y \mid t, x)$ is the conditional density of $Y$ given $T=t$ and $X=x$ with respect to a dominating measure $\mu$. It is a topic of future research to develop this idea formally.
We would like to thank a meta-reviewer and four anonymous reviewers for helpful comments.