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.
46,371 characters · 15 sections · 27 citation commands
The role of the geometric mean in case-control studies
Outcome-dependent sampling schemes like case-control sampling\footnote{In case-control sampling, `cases' are observations where the outcome occurred and `controls' are observations where the outcome did not occur. This is also known as case-referent sampling, case-noncase sampling, or choice-based sampling.} have long been used in medical, epidemiological, and survey research when the outcome is rare or the data are expensive to collect. For example, pancreatic cancer is a particularly deadly but rare cancer whose risk factors have been analyzed under case-control designs hassan2007risk, lucenteforte2012alcohol. Case-control methods are also often used in criminal justice settings to understand risk factors for crime and to audit for racial bias in police brutality loftin1988analysis, campbell2003risk, bogstrand2014drugs, wheeler2017factors, ridgeway2020role.
A case-control design samples from two distributions: the population where the outcome occurred (‘cases’) and the population where the outcome did not occur (‘controls’) cornfield1951method, miettinen1976estimability, breslow1996statistics. An alternative outcome-dependent sampling scheme, known as case-population sampling, samples from the general population instead of the control distribution lancaster1996case, jun2020causal. The broader class of biased sampling designs has a vast literature including choice-based sampling manski1977estimation, heckman2009note, stratified sampling imbens1996efficient, and exposure-based sampling such as the matched cohort design kennedy2015semiparametric.
Under outcome-dependent sampling, common effect measures such as the average risk difference and the average risk ratio are not identified, but the conditional odds ratio is identified. The conditional odds ratio is often estimated under constant effect and/or parametric assumptions. Common in practice, logistic regression estimates the conditional odds ratio under the strong assumption that the parametric form is correctly specified prentice1979logistic. A semiparametric approach relaxes these assumptions by placing restrictions directly on the odds ratio function, leaving other parts of the data-generating process unspecified chen2007semiparametric, tchetgen2010doubly. To avoid parametric assumptions altogether, a targeted maximum likelihood estimator (TMLE) with double robustness properties has been proposed for the population odds ratio, though only when the outcome rate is known van2008estimation.
The odds ratio can approximate or bound other For example, the odds ratio approximates the risk ratio under a rare disease assumption cornfield1951method, greenland1982need. A recent work uses the odds ratio to upper bound the risk ratio under monotonicity assumptions on the treatment response and treatment selection jun2020causal. To account for biased sampling, they target estimation of the outcome-conditional log odds ratio, and they propose a retrospective sieve logistic estimator to do so. The odds ratio can also be used to define a measure of interaction effect when there is more than one exposure of interest vanderweele2011weighting.
Our work considers odds ratios as the effect of interest, in a setting where the outcome rate is unknown. We compare several odds ratio-based measures of effect, discussing their usefulness (Section (ref)), identifiability (Section (ref)), and collapsibility (Section (ref)). Of note is our novel observation that the geometric mean odds ratio is collapsible, unlike the standard (arithmetic aggregated) odds ratio. We offer a new definition of collapsibility that makes explicit the choice of aggregation. We provide the partial identification of the geometric mean odds ratio under outcome-dependent sampling. We propose a doubly robust-style estimation approach for the partially identified geometric odds ratio and describe conditions under which the estimator is $\sqrt{n}$-consistent and asymptotically normal in Section (ref). Lastly we describe how to do inference for the geometric odds ratio over a range of possible values for the unknown outcome rate in Section (ref).
Define $Z = (X, A, Y)$ where $X \in \mathbb{R}^d$ are covariates, $A \in \{0,1\}$ is a binary treatment or exposure, and $Y\in \{0,1\}$ a binary outcome. Let $\mathcal{P}$ denote the space of probability distributions.
Our population of interest has distribution $P$. We observe data $(Z_1, Z_2,... Z_n)$ from a possibly biased distribution $Q$, where the bias results from outcome-dependent sampling. Let $\rho=P(Y=1)$ and $\omega=Q(Y=1)$ denote the probability that $Y=1$ under $P$ and $Q$ respectively. The quantity $\omega$ is fixed and known--or if not known, can be easily estimated from the data. The quantity $\rho$ is fixed but generally unknown. We make virtually no assumptions about knowledge of $\rho$, requiring only that it lie in a user-specified range $[\underline{\rho},~\overline{\rho}]$. We anticipate in many settings the user to generally have some knowledge about $\rho$ to inform the range, but if not, the user can specify $\underline{\rho} = \epsilon ,\overline{\rho} =1-\epsilon$.
While case-control data consists of samples from two different distributions, it can equivalently be viewed as independent and identically distributed draws from a modified distribution, known as Bernoulli sampling breslow2000semi. Under Bernoulli sampling, we first sample $Y \sim \mathrm{Bernoulli}(\omega)$, and then we sample $(X,A) \sim P(X, A \mid Y)$. We can relate $P$ and $Q$ as
We use potential (counterfactual) outcomes to define several measures of effect. Let $Y^1$ denote the potential outcome that would have been observed under treatment and $Y^0$ the potential outcome that would have been observed under no treatment. Expectations are taken with respect to the target distribution $P$ unless otherwise noted via subscript.
The population odds ratio describes the ratio of the odds of the outcome in a world where everyone is treated to the odds of the outcome in a world where no one is treated. In some settings, such as when the outcome is very rare, we may be interested in interventions targeting certain strata of the population. Then our target of interest is the conditional odds ratio, as given in the following definition.
Like many nonparametric function estimation problems, the conditional odds ratio is generally hard to estimate flexibly at fast rates. Logistic regression (without interactions) is probably the most commonly used method to estimate the conditional odds ratio in practice, but this approach will be biased when the conditional odds ratio varies across strata or when the linearity or logistic link function assumptions do not hold. On the other hand, nonparametric plugin methods suffer the curse of dimensionality. We may hope to avoid the curse of dimensionality by targeting the marginal effect since aggregated quantities can typically be estimated at faster rates even in nonparametric models. When we have randomly sampled iid observations, we can compare the importance of characterizing conditional effects against the reasonableness of the assumptions required, and we can then choose to target the marginal or conditional effect accordingly. When our data is from an outcome-dependent sample, however, this usual trade-off no longer applies because the marginal effect is not identified, as we show in the next section.
When the covariates are continuous or high-dimensional, we may desire an aggregation of the conditional odds ratio to describe a summary measure. The most common aggregation uses the arithmetic mean.
A complicating issue is the population odds ratio could be smaller (or larger) than all conditional odds ratios. In other words, the odds ratio is not collapsible greenland1986identifiability, miettinen1981confounding, robinson1991some.
\paragraph{Miscellaneous notation} $L \lesssim R$ indicates that $L \leq C \cdot R$ for some universal constant $C$. Define the squared $L_2(P)$ norm of a function $f$ as ${\left\lVertf\right\rVert^2 := \int (f(x))^2 dP(x)}$. $\mathbb{I}\{\}$ denotes the indicator function. For samples $Z_1, Z_2,.... Z_n \sim Q$, we denote sample averages $\frac{1}{n}\sum_{i=1}^n f(Z_i)$ as $Q_n(f(Z))$.
To identify our targets in terms of observable data, we make the following assumptions:
Assumption (ref) requires that there is no interference between units; e.g., the potential outcome for a unit does not depend on the treatment assignment of other units. Assumption (ref) requires that the treatment assignment is as good as random conditional on measured covariates. Assumption (ref) requires that all covariate strata have some non-zero probability of receiving both treatment decisions.\footnote{For estimation we will require a stronger boundedness assumption on overlap.} In settings where it is unreasonable to make these assumptions, sensitivity analysis can be conducted in order to assess whether the results are sensitive to violations of the assumptions rosenbaum2010design, luedtke2015statistics, Robins(00).
If our data consists of random samples from the target distribution, then we can identify the population odds ratio, conditional odds ratio, arithmetic odds ratio, and geometric odds ratio respectively as
The conditional, marginal, and aggregated effects are all identified under random sampling. The next section will consider outcome-dependent sampling, where the conditional effect remains identified, but the marginal and aggregated effects are partially identified.
When our samples $(Z_1, Z_2, ... Z_n)$ are drawn from the biased distribution $Q$, we can still identify the conditional odds ratio due to the symmetry of the odds ratio cornfield1951method:
We cannot point identify the arithmetic odds ratio and geometric odds ratio under outcome-dependent sampling because we sample from $P(X \mid Y)$, not $P(X)$. Without prior knowledge about the outcome rate $\rho$, we cannot estimate $P(X)$. We can partially identify the aggregation measures as a function of the unknown parameter $\rho$.
We first define additional notation for regression functions:
$\mu_a(x)$ is point identified under outcome-dependent sampling. We can partially identify $\nu_a(x)$ as a function of the unknown $\rho$ by applying Bayes' Rules to obtain
The population odds ratio is partially identified as
The partial identification of the arithmetic odds ratio is
The geometric odds ratio is partially identified as
Or alternatively,
This identification indicates that $\mathbb{E}[\mathrm{logit}(\mu_a(X)) \mid Y = y]$ is a key object for the geometric odds ratio, and we will see later that this function plays a central role in the efficiency theory and estimation. We use the following to denote this function:
$\psi_{a,y}$ is identified under outcome-dependent sampling. To recap, Eqs. (ref)-(ref) identify the marginal and aggregated odds ratios up to the unknown constant $\rho$.
A desirable quality for a measure of effect is that the marginal effect describes the effect for a representative unit. The property of collapsibility (Def. (ref)) formalizes this quality whittemore1978collapsibility, greenland1999confounding. In this section, we take a detour from the biased sampling design to discuss collapsibility in detail. For simplicity our example will use random sampling, but the ideas are generally applicable.
Collapsibility is often discussed with respect to the arithmetic mean, under which collapsibility requires that we can specify weights for conditional effects such that the marginal effects equals their weighted average hernan2021causal. For instance, the risk ratio, $RR = \mathbb{E}[Y^1]/\mathbb{E}[Y^0]$, is collapsible with respect to the arithmetic mean with weights $\frac{P(X)}{\mathbb{E}[Y^0]}\mathbb{E}[Y^0 \mid X]$. The odds ratio, however, is not collapsible with respect to the arithmetic mean greenland1999confounding. The population odds ratio is generally not equal to the conditional odds ratio, even if the conditional odds ratio is a constant.
\paragraph{Example:} $X = \mathbb{I}\{\mathrm{Female}\}$ with $P(X =1) = 0.5$. The risks under treatment and no treatment for women and men are given in the below table along with the corresponding conditional risk ratios, $RR(X)$, and conditional odds ratios, $OR(X)$. The marginal risk ratio with respect to arithmetic aggregation is $\frac{140}{1053} \approx 0.133$. Averaging the conditional risk ratios with weights $\frac{P(X)}{\mathbb{E}[Y^0]}\mathbb{E}[Y^0 \mid X]$ yields the marginal risk ratio under arithmetic aggregation.
\\ The conditional odds ratio for women equals the conditional odds ratio for men. However the marginal odds ratio under arithmetic aggregation is $\approx \frac{3}{2} *\frac{1}{45}$. It is not possible to find a weighted average of the conditionals that equals the marginal odds ratio.
While the odds ratio is not collapsible under the arithmetic mean, it is collapsible under the geometric mean. To demonstrate this, we introduce additional notation. Let ${f(a, b) \colon \mathbb{R}^2 \mapsto \mathbb{R}}$ denote an effect contrast and let $g_{w(x)}(P) \colon \mathcal{P} \mapsto \mathbb{R}$ denote a statistical functional that aggregates $X \sim P$ with weighting function $w(x)$. For example, letting $p(x)$ denote the density of random variable $X \sim P$, we describe the average risk difference (commonly referred to as average treatment effect) by specifying $g_{p(x)}(P) = \int x p(x) dx$ and $f(a,b) = a - b$. We consider aggregations that can be written as a Fréchet mean--that is, there is an associated distance function $d$ such that $$g_{w(x)}(P) = \displaystyle \operatorname*{arg\,min}_{z \in \mathcal{X}} \int_{\mathcal{X}} w(x) d^2(z, x) dx.$$ For ease of notation, we will write $g(X)$ to indicate $g(P)$ for the distribution $P$ over $X$.
The left hand side describes the marginal effect--that is, the contrast of the aggregations of $\mu_a(x)$ for $a \in \{0,1\}$. The right hand side describes a weighted aggregation of the conditional contrasts.
Returning to our example, the average risk difference is collapsible with respect to the arithmetic mean using as weights the density of $x$. We briefly remark on the weights $p(x)$ and $w(x)$. While $p(x) = w(x)$ for the average risk difference, this need not be the case. The risk ratio has contrast $f(a,b) = \frac{a}{b}$ and is collapsible under aggregation $g_{w(x)}(P) = \int x w(x) dx$ with $w(x) = \frac{p(x)\mathbb{E}[Y^0 \mid X = x]}{\mathbb{E}[Y^0]}$.
As far as we are aware, this expanded definition of collapsibility that explicitly incorporates the aggregation method has not appeared in the literature before. With this machinery in place, we make the novel observation that the odds ratio is collapsible under geometric aggregation.
If the conditional OR is a constant $c$, then the geometric odds ratio also equals $c$. Recall that this was not necessarily the case for the arithmetic odds ratio, which could take on a value other than $c$. Since the geometric mean exhibits the desirable property of collapsibility, the remainder of this paper will consider estimation of the geometric odds ratio. We focus on the outcome-dependent sampling design, deferring results for the random sampling design to the Appendix. To motivate our estimation approach, we start by studying the efficiency theory.
First, we restate and define additional notation for the nuisance functions:
Recall from Section (ref) (Eq.s (ref)-(ref)) that we can write our target geometric OR as a function of $\rho$ and $\psi_{a,y}$ for $(a,y) \in \{0,1\}^2$:
We will first provide a von Mises-type expansion for $\psi_{a,y}$ using as an example $a=0$ and $y=1$ (see Appendix (ref) for the general result) and subsequently provide the efficiency theory for $\gamma$. Functioning as a distributional analog to the Taylor expansion for real-valued functions, the von Mises-type expansion of a target parameter describes two key elements for efficiency theory: the influence function and a remainder term. Influence functions enable us to construct estimators with desirable properties, such as second-order bias, which can achieve fast convergence rates even in nonparametric settings. The remainder term plays an important role in characterizing the error of such estimators (see Section (ref)). In a fully nonparametric model, the singular influence function is called the efficient influence function because it characterizes the efficiency bound in a local asymptotic minimax sense. The efficient influence function is therefore instructive for constructing optimal estimators. We direct the interested reader to bickel1993efficient, tsiatis2006semiparametric, kennedy2022semiparametric, hines2022demystifying for more information on influence functions. We first define notation to refer to the nuisance functions on the distribution $\bar{Q}$: let $\bar{\eta}(x) := \bar{Q}(Y=1 \mid X=x)$ and similarly for $\bar{\pi}_0(x)$ and $\bar{\mu}_0(x)$.
For quick reference we will restate the influence function of $\gamma$ using notation that defines $\eta_{y}(x) := Q(Y=y \mid x)$ and $\omega_y := Q(Y=y)$:
where $\varphi_{a,y}(z;Q) = \mathrm{IF}(\psi_{a,y}(z;Q))$
The influence function for $\gamma$ indicates that in addition to requiring that the propensity scores be bounded away from zero and one, we will additionally require that the conditional variances $(1-\mu_1(x))\mu_1(x)$ and $(1-\mu_0(x)) \mu_0(x)$ be bounded away from zero. This is notably different from the usual risk difference (ATE) setting, where we want the conditional variances to be small to improve efficiency.
The efficiency bound describes the local asymptotic minimax lower bound on the mean squared error for any estimator of the target parameter, analogous to the Cramer-Rao bound for parametric settings. This bound provides a benchmark against which we can compare estimators. Additionally, the efficiency bound illuminates which components affect the difficulty of the estimation problem. For additional details we refer the reader to bickel1993efficient, van2003unified, tsiatis2006semiparametric, kennedy2022semiparametric. Before stating our efficiency bound, we introduce notation that will simplify the result. We define the distribution-corrected log conditional odds ratio:
and the distribution-corrected aggregate log odds ratio:
The proof is given in Appendix (ref).
The two terms that involve $\Psi^{*}(Y)$ result from having to estimate $\omega$; if $\omega$ is known by sampling design, then these terms drop from the bound. The coefficient $\big(\frac{\rho -\omega}{\omega(1-\omega)} \big)\eta(X) +\frac{1-\rho}{1-\omega} $ results from sampling bias. When $\rho=\omega$ this coefficient equals 1.\footnote{Compare to the bound under random sampling given in Appendix (ref)}
Theorem (ref) indicates that the difficulty of our estimation problem depends on the following factors:
The second line of the efficiency bound shows the efficiency bound decreases with the regression function variances, which is notably different from the average risk difference efficiency bound $\sigma^2_{ARD}$ hahn1998role:
We propose a doubly robust style estimator for $\gamma$ that relies on estimation of $\log(\gamma)$. Before providing the error analysis of our proposed estimator, we briefly remark on the role of sample splitting. Estimating our nuisance functions on a separate sample, which we denote by $\hat Q$, that is independent of the sample denoted by $Q_n$, enables us to avoid overfitting without having to rely on empirical process conditions. With iid data, we can obtain these independent samples simply by randomly partitioning the data into two or more folds. More generally, one can use cross-fitting, a procedure which swaps the samples and averages the results, to regain sample efficiency robins2008higher, zheng2010asymptotic, chernozhukov2018generic. For simplicity we present our analysis under single sample splitting. We note that the outcome rate $\omega$ can be estimated on the full data sample.
Theorem (ref) demonstrates that our proposed estimator has second-order errors in the nuisance estimation errors, yielding “doubly-fast” rates. That is, we obtain a faster rate for our estimator even when estimating the nuisance function at slower rates. For example, to obtain $n^{-1/2}$ rates for our estimator, it is sufficient to estimate the nuisance functions at $n^{-1/4}$, allowing us to use flexible machine learning methods to nonparametrically estimate the nuisance functions under smoothness or sparsity assumptions. Since our error involves squared terms, this is not the usual double-robustness property that guarantees fast rates when either of the propensity or regression function is estimated at fast rates.
Our proposed estimator for $\gamma(\rho)$ is
where $\hat{\psi}_{a,y}$ is defined in Eq. (ref).
This section discusses how to do inference when estimating $\gamma(\rho)$ over a user-specified range of values $[\underline{\rho},~\overline{\rho}]$ for the unknown outcome rate $\rho$. We note that $\gamma(\rho)$ is monotonic in $\rho$, so our bound on $\gamma(\rho)$ has as endpoints $\gamma(\underline{\rho})$ and $\gamma(\overline{\rho})$. First we discuss how to obtain a confidence interval on the endpoints, using $\gamma(\overline{\rho})$ as an example. Then we show how these imply a confidence interval on the bound for $\gamma(\rho)$. Guided by our theoretical results in the previous sections, our proposed approach is influence-function based. Alternatively, one could use the corrected confidence interval approach in imbens2004confidence to give a confidence interval for $\gamma(\rho)$.
Based on the efficient influence function of $\gamma(\rho)$ (Eq. (ref)), we create a pseudo-outcome $\zeta(Z_i; \hat \gamma(\overline{\rho}))$ for each observation $Z_i$ as
where $\hat \gamma(\overline{\rho})$ and $\hat \psi_{a,y}$ for $(a,y) \in \{0,1\}^2$ are estimates using the doubly robust approach in Sec. (ref) and where, for $a \in \{0,1\}$,
\\
We can then make use of the asymptotic normality results from the previous section. If the conditions in Corollary (ref) are met, then Corollary (ref) and Slutsky's theorem give that a $100(1-\alpha)\%$ asymptotic confidence interval for $\gamma(\overline{\rho})$ is
where $z_\beta$ is the standard normal quantile of $\beta$ and the empirical variance is over $Z$, holding as fixed the estimated nuisance functions. We can repeat this process on the same sample to obtain the confidence interval for $\gamma(\underline{\rho})$.
To obtain the confidence interval for our bound on $\rho$, we define
gives a $100(1-\alpha)\%$ asymptotic confidence interval for the bound on $\gamma(\rho)$ for $p \in [\underline{\rho}, \overline{\rho}]$.
The geometric mean odds ratio has the desirable property of collapsibility, unlike the more commonly used arithmetic mean. Under outcome-dependent sampling, the geometric odds ratio is not point identified, but we can estimate the geometric odds ratio as a function of the unknown outcome rate. We detail the efficiency theory for the geometric odds ratio, describe a doubly robust estimation procedure that is $\sqrt{n}$-consistent and asymptotically normal under mild conditions, and propose an inference procedure to construct confidence intervals for the geometric odds ratio over a range of possible values for the unknown outcome rate $\rho$.
Coston gratefully acknowledges financial support support from the National Science Foundation Graduate Research Fellowship Program under Grant No. DGE1745016. Any opinions, findings, and conclusions or recommendations expressed in this material are solely those of the authors.