Effect Inference from Two-Group Data with Sampling Bias

TL;DR

This paper introduces a new inference method that accurately compares two populations despite sampling biases, maintaining low false positive rates where standard methods fail.

Contribution

The authors develop a bias-resilient inference technique that controls false positives under moderate sampling biases, improving reliability over traditional methods.

Findings

Method performs well on synthetic data

Effective on real biomarker datasets

Reduces false positives under bias

Abstract

In many applications, different populations are compared using data that are sampled in a biased manner. Under sampling biases, standard methods that estimate the difference between the population means yield unreliable inferences. Here we develop an inference method that is resilient to sampling biases and is able to control the false positive errors under moderate bias levels in contrast to the standard approach. We demonstrate the method using synthetic and real biomarker data.

Equations

$$b=E\left[\hat{\delta }\right]-\delta ,$$

$$\mathbf{y}=\left[\begin{array}{c}{\mathbf{y}}_0\\ {\mathbf{y}}_1\end{array}\right]\sim \mathcal{N}\left(\left[\begin{array}{cc}1& 0\\ 1& 1\end{array}\right]\left[\begin{array}{c}\mu \\ \delta \end{array}\right],\left[\begin{array}{cc}{v}_0\mathbf{I}& 0\\ 0& v_1\mathbf{I}\end{array}\right]\right)$$

$$\delta \sim \mathcal{N}(0,\lambda ),$$

$$\Pr\Big{\{}\>\delta\in C_{\alpha}(\mathbf{y})\>\Big{\}}\;\geq\;1-\alpha.$$

$$\mathbf{\theta }=\mathrm{col}\{\mu ,\lambda ,v_0,v_1\}.$$

$$\begin{split}\widehat{\delta}(\mathbf{y})&=\operatorname{E}_{\theta}[\delta|\mathbf{y}]\>\Big{|}_{\theta=\widehat{\theta}}\\ &=\frac{\rho n_{1}}{\rho n_{1}+1}(\overline{y}_{1}-\mu)\>\Big{|}_{\theta=\widehat{\theta}},\end{split}$$

$$c_{\theta }^2=\frac{\rho v_1}{\rho n_1+1}+\frac{{\rho }^2{n}_1^2}{(\rho n_1+1{)}^2}{\left(\frac{n_0}{v_0}+\frac{n_1}{v_1}\frac{1}{\rho n_1+1}\right)}^{-1}.$$

$$C_{\alpha}(\mathbf{y})=\big{\{}\delta^{\prime}:|\delta^{\prime}-\widehat{\delta}(\mathbf{y})|<\alpha^{-1/2}c_{\theta}\big{\}}.$$

$$p_{\theta }(\mathbf{y})=\int p_{\theta }(\mathbf{y}|\delta )p_{\theta }(\delta )d\delta$$

$${\mathrm{Cov}}_{\theta }[\mathbf{y}]=\mathrm{diag}(v_0\mathbf{I},\lambda 11^{\top }+v_1\mathbf{I}),$$

$$\hat{\mathbf{\theta }}=\underset{\mathbf{\theta }}{\mathrm{argmax}}{p}_{\theta }(\mathbf{y}),$$

$$\begin{array}{rl}\alpha & ={\mathbf{y}}_1^{\top }{\mathbf{y}}_1-\mu n_1(2{\overline{y}}_1-\mu )\\ \beta & =n_1^2({\overline{y}}_1-\mu {)}^2\\ \gamma & =\alpha n_1-\beta .\end{array}$$

$${\hat{v}}_0=\frac{1}{n_0}{\mathbf{y}}_0^{\top }{\mathbf{y}}_0-\mu (2{\overline{y}}_0-\mu ),$$

$$\begin{array}{r}{\hat{v}}_1=\frac{1}{n_1}\frac{\alpha +\rho \gamma }{1+\rho n_1},\end{array}$$

$$\begin{array}{r}\hat{\rho }=\{\begin{array}{ll}\frac{\beta -\alpha }{\gamma },& \beta -\alpha \ge 0.\\ 0,& \text{otherwise.}\end{array}\end{array}$$

$$\begin{array}{rl}f(\mu )& =n_0\ln {\hat{v}}_0+n_1\ln (\alpha +\hat{\rho }\gamma )-(n_1-1)\ln (1+\hat{\rho }{n}_1)\end{array}$$

$$\begin{array}{rl}E\left[|\delta -\hat{\delta }{|}^2\right]& =E_y\left[E_{\delta |y}\left[|\delta -\overline{\delta }+\overline{\delta }-\hat{\delta }{|}^2\right]\right]\\ & =E_y\left[\mathrm{Var}[\delta |\mathbf{y}]+|\overline{\delta }-\hat{\delta }{|}^2\right]\\ & =\frac{\lambda v_1}{\lambda n_1+v_1}+E_y\left[|\overline{\delta }-\hat{\delta }{|}^2\right].\end{array}$$

$$\begin{array}{r}\mathbf{\varphi }\triangleq {\partial }_{\theta }\ln p_{\theta }(\mathbf{y})\text{and}\mathbf{J}=E_y[\mathbf{\varphi }{\mathbf{\varphi }}^{\top }].\end{array}$$

$$\begin{array}{rl}\mathbf{J}& =\left[\begin{array}{cc}{J}_{1,1}& 0\\ 0& \ast \end{array}\right],\end{array}$$

$$\begin{array}{rl}{J}_{1,1}& =1^{\top }\left[\begin{array}{cc}{v}_0^{-1}\mathbf{I}& 0\\ 0& {\mathbf{V}}_1^{-1}\end{array}\right]1\\ & =v_0^{-1}{1}^{\top }1+1^{\top }{\mathbf{V}}_1^{-1}1=\frac{n_0}{v_0}+\frac{n_1}{\lambda n_1+v_1}\end{array}$$

$$\begin{array}{rl}0& \le E_y\left[|(\overline{\delta }-\hat{\delta })-{\mathbf{g}}^{\top }{\mathbf{J}}^{-1}\mathbf{\varphi }{|}^2\right]\\ & =E_y\left[|\overline{\delta }-\hat{\delta }{|}^2\right]-{\mathbf{g}}^{\top }{\mathbf{J}}^{-1}\mathbf{g}.\end{array}$$

$$\begin{array}{rl}\mathbf{g}& =\int [{\partial }_{\theta }\ln p_{\theta }](\overline{\delta }-\hat{\delta })p_{\theta }d\mathbf{y}\\ & =\int {\partial }_{\theta }[p_{\theta }(\overline{\delta }-\hat{\delta })]-p_{\theta }[{\partial }_{\theta }(\overline{\delta }-\hat{\delta })]d\mathbf{y}\\ & ={\partial }_{\theta }[\text{bias}(\mathbf{\theta })]-E_y[{\partial }_{\theta }(\overline{\delta }-\hat{\delta })]\\ & =-E_y\left[{\partial }_{\theta }\left(\frac{\lambda }{\lambda n_1+v_1}{1}^{\top }({\mathbf{y}}_1-\mu 1)\right)\right]\\ & =-\left[\begin{array}{c}\frac{\lambda }{\lambda n_1+v_1}{1}^{\top }1+0\\ {\partial }_{\lambda }\frac{\lambda }{\lambda n_1+v_1}{1}^{\top }{E}_y[{\mathbf{y}}_1-\mu 1]\\ 0\\ {\partial }_{v_1}\frac{\lambda }{\lambda n_1+v_1}{1}^{\top }{E}_y[{\mathbf{y}}_1-\mu 1]\end{array}\right]=-\left[\begin{array}{c}\frac{\lambda n_1}{\lambda n_1+v_1}\\ 0\\ 0\\ 0\end{array}\right],\end{array}$$

$$\begin{array}{r}E\left[|\delta -\hat{\delta }{|}^2\right]\ge \frac{\lambda v_1}{\lambda n_1+v_1}+{\left(\frac{\lambda n_1}{\lambda n_1+v_1}\right)}^2{J}_{1,1}^{-1}.\end{array}$$

$$\begin{array}{r}\delta \notin C_{\alpha }(\mathbf{y})\iff \alpha c_{\theta }^{-2}|\delta -\hat{\delta }(\mathbf{y}){|}^2\ge 1.\end{array}$$

$$\begin{split}\Pr\big{\{}\delta\not\in C_{\alpha}(\mathbf{y})\big{\}}&=\int_{\delta\not\in C_{\alpha}(\mathbf{y})}p(\mathbf{y},\delta)d\delta d\mathbf{y}\\ &\leq\int_{\delta\not\in C_{\alpha}(\mathbf{y})}\alpha c^{-2}_{\theta}|\delta-\widehat{\delta}(\mathbf{y})|^{2}p(\mathbf{y},\delta)d\delta d\mathbf{y}\\ &\leq\alpha c^{-2}_{\theta}\operatorname{E}\left[|\delta-\widehat{\delta}(\mathbf{y})|^{2}\right]=\alpha\frac{\text{MSE}}{c^{2}_{\theta}}.\end{split}$$

$$\begin{array}{rl}f(\mathbf{\theta })& =\underset{f_0(\mu ,v_0)}{\underbrace{n_0\ln v_0+\frac{1}{v_0}\Vert {\mathbf{y}}_0-\mu 1{\Vert }^2}}+\underset{f_1(\mu ,\lambda ,v_1)}{\underbrace{\ln |{\mathbf{V}}_1|+\Vert {\mathbf{y}}_1-\mu 1{\Vert }_{{\mathbf{V}}_1^{-1}}^2}}.\end{array}$$

$${\hat{v}}_0=\Vert {\mathbf{y}}_0-\mu 1{\Vert }^2/n_0$$

$$\begin{array}{rl}{f}_0(\mu ,{\hat{v}}_0)& =n_0\ln {\hat{v}}_0+n_0\end{array}$$

$$\begin{array}{rl}{f}_1(\mu ,\rho ,v)& =\ln (1+\rho n)+\ln v^n\\ & +\frac{1}{v}\left(\Vert {\mathbf{y}}_1-\mu 1{\Vert }^2-\frac{\rho |1^{\top }({\mathbf{y}}_1-\mu 1){|}^2}{1+\rho n}\right)\end{array}$$

$$\begin{array}{rl}{f}_1(\mu ,\rho ,\hat{v})& =\ln \frac{(\alpha +\rho \gamma {)}^n}{(1+\rho n{)}^{n-1}}+n.\end{array}$$

Full text

I Introduction

In many applications of statistical inference, the aim is to compare data from different populations. Specifically, given $n_{0}$ and $n_{1}$ samples from two groups, collected in vectors $\mathbf{y}_{0}$ and $\mathbf{y}_{1}$, the target quantity is often the difference between their means, denoted $\delta$, which we call the effect. For instance, in randomized trials and A/B testing, the data are outcomes from two populations and $\delta$ is the average causal effect of assigning subjects to a test group ‘$1$’ as compared to a control group ‘[math]’. [1, 2] The standard approach is to use the difference between sample averages in each group, viz. $\widehat{\delta}=\overline{y}_{1}-\overline{y}_{0}$, where $\overline{y}_{i}=\mathbf{1}^{\top}\mathbf{y}_{i}/n_{i}$. Confidence intervals for $\widehat{\delta}$ can be obtained using Welch’s method, which employs an approximating t-distribution [3, 4, 5]. Inferring $\delta\neq 0$ is equivalent to detecting that the means of two distributions differ, which is a classical problem in statistical signal processing [6, 7].

Ideally, the samples from both groups are representative of their target populations. Then the bias of the estimator,

[TABLE]

is zero. However, in nonideal conditions with finite samples this is not the case, e.g., when some units of the intended populations are less likely to be included than others. Under such conditions, $b$ decreases with sample sizes $n_{0}$ and $n_{1}$ but will nevertheless be nonzero. Sampling biases increase the risk of inferring spurious effects when using standard inference methods.

In this paper, we develop an inference method that is resilient to sampling biases. In contrast to the standard approach, the proposed method reduces the risk of reporting spurious effect estimates and is capable of controlling the false positive errors under moderate biases. The method relies on an effect estimator using a fully automatic and data-adaptive regularization. We demonstrate its performance on both synthetic and real data.

Remark 1*.*

Code for the method can be found at https://github.com/dzachariah/two-groups-data

II Problem formulation

We model the dataset as

[TABLE]

The model based on the Gaussian distribution yields the least favourable distribution for estimating the unknown effect $\delta$ [8]. We model the effect as a random variable, where different ranges of values of $\delta$ have different probabilities. To achieve resilliance to sampling biases, we adopt a conservative approach in which nonexistant or negligible effects are considered to be more probable. Specifically, we employ the following model:

[TABLE]

where $\lambda$ is an unknown parameter.

Our aim is to derive a confidence interval $C_{\alpha}(\mathbf{y})$ that contains the unknown $\delta$ with a coverage probability of at least $1-\alpha$. That is,

[TABLE]

The confidence interval is to be centered on an estimator $\widehat{\delta}(\mathbf{y})$ and should be resilient to sampling biases. That is, even if $b\neq 0$ the interval must not indicate nonzero effects with a probability greater than $\alpha$. Fig. 1 illustrates the ability of the method proposed below to ensure (3) under a range of biases, provided $b$ does not greatly exceed the dispersion of sample averages, i.e., $\sqrt{v_{i}/n_{i}}$.

We will derive a confidence interval using model (1) and (2), with nuisance parameters

[TABLE]

III Proposed method

Let $\operatorname{E}_{\theta}[\delta|\mathbf{y}]$ be the conditional mean of the effect given the data. Using an estimate $\widehat{\boldsymbol{\theta}}$ of the nuisance parameters, we propose the following effect estimator

[TABLE]

where we introduce the variable $\rho\equiv\lambda/v_{1}$ that can be interpreted as a signal-to-noise ratio, see [9] for a derivation.

Result 1** (Cramér-Rao bound).**

When the systematic error of $\widehat{\delta}(\mathbf{y})$ is invariant with respect to $\boldsymbol{\theta}$, then the mean-squared error over all possible effects and data has a Cramér-Rao bound $\operatorname{E}\left[|\delta-\widehat{\delta}(\mathbf{y})|^{2}\right]\geq c^{2}_{\theta},$ where

[TABLE]

Proof.

See Appendix -A. ∎

Result 2** (Confidence interval).**

Let

[TABLE]

When using an efficient estimator that attains the bound (5), the interval in (6) satifies the specified coverage probability (3).

Proof.

See Appendix -B. ∎

Evaluating $\widehat{\delta}(\mathbf{y})$ and $C_{\alpha}(\mathbf{y})$ requires estimates of the nuisance parameters $\boldsymbol{\theta}$. Here we adopt the maximum likelihood approach and estimate $\boldsymbol{\theta}$ using the marginalized data distribution,

[TABLE]

It can be shown that (7) is a Gaussian distribution [9] with mean $\operatorname{E}_{\theta}[\mathbf{y}]=\mathbf{1}$ and covariance

[TABLE]

The estimated parameters are given by

[TABLE]

which can be shown to yield an asymptotically efficient estimator (4) [10, corr. 9].

Interestingly, the problem (8) can be solved by a one-dimensional numerical search. Begin by defining the variables

[TABLE]

Note that $\gamma\geq 0$. Then the following result holds.

Result 3** (Nuisance parameter estimates).**

The estimated variances are given by

[TABLE]

[TABLE]

which are ensured to be nonnegative, and $\widehat{\lambda}=\widehat{\rho}\widehat{v}_{1}$, where

[TABLE]

All variables in (9)-(11) are functions of the mean $\mu$, whose estimate $\widehat{\mu}$ is obtained by minimizing the one-dimensional function

[TABLE]

Proof.

See Appendix -C. ∎

By plugging in $\widehat{\mu}$, $\widehat{\rho}$, $\widehat{v}_{0}$ and $\widehat{v}_{1}$ into (4) and (6), we obtain estimates $\widehat{\delta}(\mathbf{y})$ and $C_{\alpha}(\mathbf{y})$, respectively. We note that the overall mean $\mu$ is fitted to the data in a nonstandard manner using (12), which yields a fully automatic and data-adaptive regularization of the effect estimator (4). If the minimizing $\widehat{\mu}$ is such that $\beta<\alpha$, then the estimated signal-to-noise ratio is $\widehat{\rho}=0$. In this case, the method indicates that the data is not sufficiently informative to discriminate any systematic difference from noise. Consequently, $\widehat{\delta}(\mathbf{y})$ collapses to zero and $C_{\alpha}(\mathbf{y})=\emptyset$, indicating a case in which the effect cannot be reliably inferred.

IV Experimental results

We demonstrate the proposed inference method using both synthetic and real data.

IV-A Synthetic data

We generate two-group data using the model (1) and add a negative bias $b$ to the test group, using the setup parameters described in Fig. 1. The adaptive regularization of $\widehat{\delta}$ is illustrated in Fig. 2: when the unknown effect is nonexistent, $\delta=0$, the estimates are concentrated at zero, despite the bias $b$. As $\delta$ exceeds the dispersion of the sample averages, however, the regularized and standard estimators become nearly identical.

We report a significant effect estimate when a nonempty interval $C_{\alpha}(\mathbf{y})$ excludes the zero effect. Fig. 3 illustrates the ability of the proposed method to control the false positive error probability as $n_{0}$ increases, in contrast to the standard method. This is achieved while incurring a loss of statistical power that vanishes as the number of samples increases.

IV-B Prostate cancer data

We now consider real data from $n_{0}=50$ healthy individuals and $n_{1}=52$ individuals with prostate cancer [11, 12]. The data contains 6033 different biomarker responses. The inferred effects are shown in Fig. 4. For 6 markers, the effects were found to be significant at the $\alpha=0.05$ level. By contrast, the standard approach using Welch’s t-intervals yields 478 genes, but the inferences are less reliable under sampling biases.

V Conclusions

We developed a method for inferring effects in two-group data that, unlike the standard approach, is resilient to sampling biases. The method is able to control the false positive errors under moderate bias levels and its performance was demonstrated using both synthetic and real biomarker data.

-A The derivation of the Cramér-Rao bound

The mean-square error can be decomposed as

[TABLE]

where $\overline{\delta}$ is the conditional mean. Next, define the score function and the information matrix,

[TABLE]

Since the marginal pdf is Gaussian, we can compute $\mathbf{J}$ using Slepian-Bangs formula [13]. It has a block diagonal form

[TABLE]

where

[TABLE]

and $\mathbf{V}_{1}=\lambda\mathbf{1}\mathbf{1}^{\top}+v_{1}\mathbf{I}$.

Let $\mathbf{g}\triangleq\operatorname{E}_{y}[\boldsymbol{\phi}(\overline{\delta}-\widehat{\delta})]$ denote the correlation between the score function and estimation error. Then we have the general bound

[TABLE]

In our case, we obtain

[TABLE]

where the fourth line follows under the constant bias assumption. Inserting this expression for $\mathbf{g}$ in (17) yields

[TABLE]

This completes the proof.

-B The derivation of the confidence interval

We have that

[TABLE]

Let $p(\mathbf{y},\delta)=p_{\theta}(\mathbf{y}|\delta)p_{\theta}(\delta)$, then

[TABLE]

Thus \Pr\big{\{}\delta\in C_{\alpha}(\mathbf{y})\big{\}}\geq 1-\alpha when the estimator is efficient.

-C The derivation of the concentrated cost

Problem (8) can be formulated equivalently as the minimization of:

[TABLE]

The minimizer

[TABLE]

is inserted back to yield a concentrated cost function

[TABLE]

Next, using the Sherman-Morrison and matrix determinant lemmas we can reparametrize $f_{1}$ as

[TABLE]

where we dropped the subindices for notational convenience.

Using the identities $\alpha=\|\mathbf{y}_{1}-\mu\mathbf{1}\|^{2}$, $\beta=|\mathbf{1}^{\top}(\mathbf{y}_{1}-\mu\mathbf{1})|^{2}$ and $\gamma=\alpha n-\beta$, the minimizing $v$ of (23) is found as (10). Inserting the variance estimate back, yields a concentrated cost function

[TABLE]

To find the minimizing $\rho\geq 0$, we first consider the stationary point of

[TABLE]

Taking the derivative with respect to $\rho$, yields the following condition for a stationary point:

[TABLE]

or equivalently $\gamma(1+\rho n)-(n-1)(\alpha+\rho\gamma)=0$. Solving for $\rho\geq 0$, we obtain the estimate (11).

By evaluating the second derivative at this point, we verify that it is a minimum. Inserting (11) back into (24) and combining with (22), we can write (20) in the concentrated form (12) after omitting irrelevant constants.