Heteroscedasticity-Adjusted Ranking (HART)

• HART is a statistical framework for multiple testing that accounts for heteroscedasticity by using both summary statistics and their associated variances.
• It employs a conditional local false discovery rate to rank hypotheses, enabling more effective thresholding than standard z-score based methods.
• Empirical studies show that HART increases detection power in large, heterogeneous datasets by leveraging variance information ignored in conventional approaches.

HART, short for Heterocedasticity-Adjusted Ranking and Thresholding, is a two-step multiple testing framework for large-scale heterogeneous data that bypasses standardization and uses both a summary statistic and its unit-specific variance for ranking and thresholding. Its central claim is that standardization to $z$-scores or $p$-values can discard information carried by the variance structure of the alternative distribution, especially under heteroscedasticity. The methodology, theory, and empirical studies were developed by Fu, Gang, James, and Sun in “Heterocedasticity-Adjusted Ranking and Thresholding for Large-Scale Multiple Testing” (Fu et al., 2019).

1. Conceptual basis and statistical setting

HART is designed for settings in which, for $i=1,\dots,m$, one observes a summary statistic $X_i$ and a variance or standard deviation $s_i^2$ or $s_i$, with heteroscedasticity across units. The inferential task is

$$H_{0,i}:\mu_i=0 \qquad \text{versus} \qquad H_{1,i}:\mu_i\neq 0,$$

with signal indicator $\theta_i=1(\mu_i\neq 0)$ following a $\mathrm{Bernoulli}(\pi)$ model (Fu et al., 2019).

The canonical probabilistic model is Gaussian with heteroscedastic errors:

$$X_i \mid \mu_i,\sigma_i^2 \overset{ind}{\sim} N(\mu_i,\sigma_i^2),\qquad \mu_i \overset{iid}{\sim} (1-\pi)\delta_0+\pi g_\mu,\qquad \sigma_i^2 \overset{iid}{\sim} g_\sigma,$$

where $p$0 is a point mass at $p$1. In much of the analysis, $p$2 are treated as known or consistently estimated (Fu et al., 2019).

Conditioning on $p$3, the marginal mixture for $p$4 is

$p$5

with

$p$6

The motivating observation is that standardization to $p$7 maps all nulls to $p$8, but also transforms the alternative to

$p$9

This subsumes the variance structure of the alternative. Consequently, for the same $i=1,\dots,m$0-value, small-$i=1,\dots,m$1 units are more consistent with the alternative than large-$i=1,\dots,m$2 units. HART exploits this by conditioning on $i=1,\dots,m$3 in both ranking and thresholding (Fu et al., 2019).

2. Oracle ranking and thresholding rule

The core HART ranking statistic is the local false discovery rate conditional on $i=1,\dots,m$4:

$i=1,\dots,m$5

This is a Bayesian likelihood-ratio-based score; smaller $i=1,\dots,m$6 indicates stronger evidence against the null (Fu et al., 2019).

Let $i=1,\dots,m$7 denote the marginal false discovery rate level, $i=1,\dots,m$8, of the rule that rejects all $i=1,\dots,m$9 with $X_i$0. The oracle HART rule is

$X_i$1

This ranking-plus-threshold construction uses both $X_i$2 and $X_i$3 explicitly and is optimal within the class of rules based on $X_i$4 (Fu et al., 2019).

The contrast with standardized procedures is explicit. For standardized values $X_i$5, the classical local false discovery rate is

$X_i$6

and the oracle $X_i$7-value rule is $X_i$8. The oracle $X_i$9-value rule is $s_i^2$0. These procedures do not condition on $s_i^2$1 and therefore discard heteroscedastic information (Fu et al., 2019).

A key structural special case is the homoscedastic regime $s_i^2$2. Then $s_i^2$3 is a monotone transform of $s_i^2$4, and HART reduces to the standard local false discovery rate procedure (Fu et al., 2019).

3. Error control, power, and asymptotic theory

HART is developed for false discovery rate control while maximizing expected true positives. The paper defines

$s_i^2$5

where

$s_i^2$6

Power is measured by

$s_i^2$7

The oracle optimality result states that if $s_i^2$8 is the collection of all procedures based on $s_i^2$9 with $s_i$0, then

$s_i$1

In particular,

$s_i$2

The intuition given in the paper is that, conditional on $s_i$3, thresholding $s_i$4 is a monotone likelihood ratio decision under the mixture model, so its rejection region maximizes alternative mass at fixed expected false discoveries (Fu et al., 2019).

For the data-driven procedure, HART estimates both $s_i$5 and $s_i$6. Its asymptotic analysis assumes four conditions. First, bounded supports:

$s_i$7

Second, regularity of the kernel $s_i$8: positive, bounded, symmetric, with

$s_i$9

and $$H_{0,i}:\mu_i=0 \qquad \text{versus} \qquad H_{1,i}:\mu_i\neq 0,$$0 having bounded, continuous second derivative and being square-integrable. Third, bandwidth conditions:

$$H_{0,i}:\mu_i=0 \qquad \text{versus} \qquad H_{1,i}:\mu_i\neq 0,$$1

$$H_{0,i}:\mu_i=0 \qquad \text{versus} \qquad H_{1,i}:\mu_i\neq 0,$$2

Fourth,

$$H_{0,i}:\mu_i=0 \qquad \text{versus} \qquad H_{1,i}:\mu_i\neq 0,$$3

Under these conditions, the consistency proposition states

$$H_{0,i}:\mu_i=0 \qquad \text{versus} \qquad H_{1,i}:\mu_i\neq 0,$$4

The asymptotic validity and optimality theorem further gives

$$H_{0,i}:\mu_i=0 \qquad \text{versus} \qquad H_{1,i}:\mu_i\neq 0,$$5

and

$$H_{0,i}:\mu_i=0 \qquad \text{versus} \qquad H_{1,i}:\mu_i\neq 0,$$6

where $$H_{0,i}:\mu_i=0 \qquad \text{versus} \qquad H_{1,i}:\mu_i\neq 0,$$7 denotes the data-driven HART rule (Fu et al., 2019).

4. Estimation strategy and computational construction

HART estimates the conditional mixture density through a stabilized decomposition,

$$H_{0,i}:\mu_i=0 \qquad \text{versus} \qquad H_{1,i}:\mu_i\neq 0,$$8

where $$H_{0,i}:\mu_i=0 \qquad \text{versus} \qquad H_{1,i}:\mu_i\neq 0,$$9 is a weighted bivariate kernel estimator with variable bandwidth in $\theta_i=1(\mu_i\neq 0)$0 to account for heteroscedasticity:

$\theta_i=1(\mu_i\neq 0)$1

The weights $\theta_i=1(\mu_i\neq 0)$2 approximate $\theta_i=1(\mu_i\neq 0)$3 and are refined in two steps (Fu et al., 2019).

The initialization uses

$\theta_i=1(\mu_i\neq 0)$4

with

$\theta_i=1(\mu_i\neq 0)$5

The updated weights are

$\theta_i=1(\mu_i\neq 0)$6

The final score is

$\theta_i=1(\mu_i\neq 0)$7

Thresholding is empirical. After sorting

$\theta_i=1(\mu_i\neq 0)$8

HART selects

$\theta_i=1(\mu_i\neq 0)$9

and rejects $\mathrm{Bernoulli}(\pi)$0 (Fu et al., 2019).

The naive computational cost of evaluating $\mathrm{Bernoulli}(\pi)$1 for all $\mathrm{Bernoulli}(\pi)$2 is $\mathrm{Bernoulli}(\pi)$3. The paper notes that fast kernel summation, including KD-trees, FFT-based approximations, sparsity via compact kernels, batching, or GPUs, can reduce wall time. A jackknife variant, implemented by leaving one point out in the density estimate, is described as stabilizing estimates with negligible asymptotic cost (Fu et al., 2019).

The base model assumes $\mathrm{Bernoulli}(\pi)$4 independent of $\mathrm{Bernoulli}(\pi)$5, but the paper also gives a structured extension for dependence:

$\mathrm{Bernoulli}(\pi)$6

with

$\mathrm{Bernoulli}(\pi)$7

Misspecifying this dependence may reduce power, while the basic HART machinery still yields valid ranking under mild conditions (Fu et al., 2019).

5. Empirical behavior and applications

The paper’s finite-mixture illustration uses $\mathrm{Bernoulli}(\pi)$8, $\mathrm{Bernoulli}(\pi)$9, $$X_i \mid \mu_i,\sigma_i^2 \overset{ind}{\sim} N(\mu_i,\sigma_i^2),\qquad \mu_i \overset{iid}{\sim} (1-\pi)\delta_0+\pi g_\mu,\qquad \sigma_i^2 \overset{iid}{\sim} g_\sigma,$$0, and $$X_i \mid \mu_i,\sigma_i^2 \overset{ind}{\sim} N(\mu_i,\sigma_i^2),\qquad \mu_i \overset{iid}{\sim} (1-\pi)\delta_0+\pi g_\mu,\qquad \sigma_i^2 \overset{iid}{\sim} g_\sigma,$$1. The average power of three oracle procedures is reported as

$$X_i \mid \mu_i,\sigma_i^2 \overset{ind}{\sim} N(\mu_i,\sigma_i^2),\qquad \mu_i \overset{iid}{\sim} (1-\pi)\delta_0+\pi g_\mu,\qquad \sigma_i^2 \overset{iid}{\sim} g_\sigma,$$2

The full-data oracle therefore more than doubles the average power relative to the $$X_i \mid \mu_i,\sigma_i^2 \overset{ind}{\sim} N(\mu_i,\sigma_i^2),\qquad \mu_i \overset{iid}{\sim} (1-\pi)\delta_0+\pi g_\mu,\qquad \sigma_i^2 \overset{iid}{\sim} g_\sigma,$$3-value oracle (Fu et al., 2019).

In continuous-heteroscedastic simulations,

$$X_i \mid \mu_i,\sigma_i^2 \overset{ind}{\sim} N(\mu_i,\sigma_i^2),\qquad \mu_i \overset{iid}{\sim} (1-\pi)\delta_0+\pi g_\mu,\qquad \sigma_i^2 \overset{iid}{\sim} g_\sigma,$$4

with $$X_i \mid \mu_i,\sigma_i^2 \overset{ind}{\sim} N(\mu_i,\sigma_i^2),\qquad \mu_i \overset{iid}{\sim} (1-\pi)\delta_0+\pi g_\mu,\qquad \sigma_i^2 \overset{iid}{\sim} g_\sigma,$$5 and $$X_i \mid \mu_i,\sigma_i^2 \overset{ind}{\sim} N(\mu_i,\sigma_i^2),\qquad \mu_i \overset{iid}{\sim} (1-\pi)\delta_0+\pi g_\mu,\qquad \sigma_i^2 \overset{iid}{\sim} g_\sigma,$$6, all methods controlled FDR at nominal levels, while data-driven HART was slightly conservative at times. The reported power pattern is that data-driven HART uniformly exceeded BH and the $$X_i \mid \mu_i,\sigma_i^2 \overset{ind}{\sim} N(\mu_i,\sigma_i^2),\qquad \mu_i \overset{iid}{\sim} (1-\pi)\delta_0+\pi g_\mu,\qquad \sigma_i^2 \overset{iid}{\sim} g_\sigma,$$7-oracle proxy, often closing much of the gap to the full-data oracle; gains increased with stronger heterogeneity and sparser signals, with relative power improvements frequently in the 20–60% range versus BH and ZOR at matched FDR (Fu et al., 2019).

In the two-group heteroscedastic setting, with $$X_i \mid \mu_i,\sigma_i^2 \overset{ind}{\sim} N(\mu_i,\sigma_i^2),\qquad \mu_i \overset{iid}{\sim} (1-\pi)\delta_0+\pi g_\mu,\qquad \sigma_i^2 \overset{iid}{\sim} g_\sigma,$$8 equiprobable, $$X_i \mid \mu_i,\sigma_i^2 \overset{ind}{\sim} N(\mu_i,\sigma_i^2),\qquad \mu_i \overset{iid}{\sim} (1-\pi)\delta_0+\pi g_\mu,\qquad \sigma_i^2 \overset{iid}{\sim} g_\sigma,$$9, $p$00, $p$01, and $p$02, the data-driven procedure approximately matched the oracle in both FDR and power. HART applied differing $p$03-cutoffs per group, increasing the cutoff for the large-variance group, and the power gap relative to ZOR grew with $p$04 (Fu et al., 2019).

The principal empirical case study is a microarray analysis of myeloma. The dataset contained 12,625 genes; 36 myeloma patients had MRI-detected focal lesions and 137 had no lesions. Gene-wise contrasts $p$05 and standard errors $p$06 were computed. Because the theoretical $p$07 null was too narrow, BH and AZ used the empirical null $p$08. For stable HART density estimation, genes with $p$09 were analyzed, leaving 12,172 genes. At $p$10, BH discovered 8 genes (0.07%), AZ 25 (0.2%), and HART 122 (1%). The rejection boundaries for BH and AZ depended only on $p$11, whereas HART’s decisions varied with both $p$12 and $p$13, favoring small-$p$14 genes and avoiding large-$p$15 genes (Fu et al., 2019).

6. Relation to other procedures, practical regime, and limitations

HART is explicitly positioned against several established multiple-testing procedures. BH on standardized $p$16-values or $p$17-scores ignores structure in the alternative and often uses symmetric two-sided rejection. Weighted BH, IHW, and AdaPT use covariates to pre-order or weight $p$18-values, but naïvely ordering by $p$19 can be anti-informative because correctness depends on the relationship between $p$20 and signal prevalence or strength. HART instead performs a semiparametric empirical Bayes integration of $p$21 through conditional density modeling (Fu et al., 2019).

In the special case where $p$22 takes one of a few discrete values, HART reduces to conditional local false discovery rate analyzed separately per group, with group-specific thresholds. Relative to classical local false discovery rate procedures that estimate $p$23, HART estimates $p$24 and $p$25 conditionally, restoring the alternative variance structure lost after standardization (Fu et al., 2019).

The paper identifies the regime in which HART is most beneficial: clear heteroscedasticity, effect sizes whose discriminability depends on $p$26, and large $p$27, typically in the thousands or more, so that conditional density estimation is accurate. The recommended diagnostics are a histogram of $p$28, a scatter of $p$29 versus $p$30, empirical $p$31 plotted against $p$32 by $p$33 strata, and direct comparison of HART and $p$34-based rejection boundaries. For tuning, the paper recommends estimating $p$35 by methods robust to alternatives such as Jin–Cai, choosing bandwidths by rules of thumb on a subset enriched for signals, and considering jackknifed density estimation. For unknown $p$36, it recommends consistent variance estimators and, if needed, $p$37-approximations for $p$38 (Fu et al., 2019).

The limitations are equally explicit. HART requires large $p$39 for stable bivariate density estimation and is computationally heavier than BH or ZOR, although it may be faster than some covariate-adaptive procedures such as AdaPT in many settings. Its theoretical development assumes independence or weak dependence and known null scale per $p$40. The extension of empirical-null modeling under heteroscedasticity is nontrivial when the theoretical null is inadequate. Finally, misspecification of $p$41–$p$42 dependence can reduce power, motivating structured extensions with $p$43 and $p$44 (Fu et al., 2019).

HART’s defining methodological claim is therefore narrow but consequential: when unit-specific variances are themselves informative under the alternative, ranking by

$p$45

and thresholding to control FDR can strictly dominate procedures based only on standardized statistics. In the framework developed in (Fu et al., 2019), this gain is formalized by oracle optimality, asymptotic validity, and repeated empirical power improvements at fixed FDR.