Filter one-sided tests before FDR

Question

In this blog post, the author shows one example of "weird" p-value histogram ("Scenario C"):

One provided explanation is that a one-sided test was run, and the tests where the effect was in the opposite direction gave p-values close to 1. His recommendation is to filter out the tests based on direction of effect before computing the FDR.

However, I am aware that pre-filtering before FDR computation can be problematic if the filtering is not done the right way, as discussed e.g. here, and it feels like filtering on effect direction is "cheating", as we look at the results.

Now, if I run a quick simulation in R (shown below), filtering leads to p-values distributed from 0 to 0.5, so the assumption that the p-values are uniform between 0 and 1 is no longer valid. So, is it correct to filter based on effect direction before FDR computation?

Further, looking empirically at the FDR and power in my simulations, it looks like I keep a correct FDR control if I specify the total number of tests in p.adjust(), while still gaining a bit of power. Is this the correct way to do it?

Simulation

set.seed(1)

## Generate data ----
real_means <- c(
  rep(0, 1000), # no effect
  runif(200, -1, 0), # decreased
  runif(200, 0, 1) # increased
)
labels <- c(
  rep("none", 1000),
  rep("decreased", 200),
  rep("increased", 200)
)

samples <- lapply(real_means, \(mu) rnorm(30, mu, sd = 1))



## Compute p-values and fdr ----
p_values <- sapply(samples, \(x) t.test(x, alternative = "greater")$p.value )
p_values <- setNames(p_values, labels)

hist(p_values, breaks = 70)

fdr <- p.adjust(p_values, method = "BH")


## check results ----
false_positives <- sum( fdr[names(fdr) != "increased"] < 0.05 )
predicted_positives <- sum( fdr < 0.05 )
true_positives <- sum( fdr[names(fdr) == "increased"] < 0.05 )
positives <- sum( names(fdr) == "increased" )

# false discovery rate
false_positives / predicted_positives
#> [1] 0.01234568

# true positive rate
true_positives / positives
#> [1] 0.4



## With filtering ----

direction_increasing <- sapply(samples, \(x) mean(x) >= 0 )


filtered_p_values <- p_values[ direction_increasing ]

hist(filtered_p_values)

filt_fdr <- p.adjust(filtered_p_values, method = "BH" )


false_positives <- sum( filt_fdr[names(filt_fdr) != "increased"] < 0.05 )
predicted_positives <- sum( filt_fdr < 0.05 )
true_positives <- sum( filt_fdr[names(filt_fdr) == "increased"] < 0.05 )
positives <- sum( names(filt_fdr) == "increased" )

# false discovery rate
false_positives / predicted_positives
#> [1] 0.04301075

# true positive rate
true_positives / positives
#> [1] 0.4863388



## With filtering, but specifying total number of tests ----

filt_fdr_with_n <- p.adjust(filtered_p_values, method = "BH", n = length(p_values) )


false_positives <- sum( filt_fdr_with_n[names(filt_fdr_with_n) != "increased"] < 0.05 )
predicted_positives <- sum( filt_fdr_with_n < 0.05 )
true_positives <- sum( filt_fdr_with_n[names(filt_fdr_with_n) == "increased"] < 0.05 )
positives <- sum( names(filt_fdr_with_n) == "increased" )

# false discovery rate
false_positives / predicted_positives
#> [1] 0.01234568

# true positive rate
true_positives / positives
#> [1] 0.4371585

^{Created on 2024-07-19 with reprex v2.1.0}

Answer 1

Thanks for sharing the post. Your simulation is correct.

If some test objects have a positive mean and some have a negative mean. We can do hypothesis tests on $H_1:\mu<0$, $H_1:\mu=0$ and $H_1:\mu>0$. Two one-sided tests will generate the Scenario C plot and look like the mirrored version of each other. The two-sided test will generate the Scenario A plot. If this is the case, we may believe the situation is the same as the simulation.

Then we can apply FDR control on two-side test results. After that, objects with positve mean can be selected. There is no need for filtering before applying FDR control.

Both FDR and TPR improved slightly in this set.seed(1) case.

### Do two-side test
p_values_twoside <- sapply(samples, \(x) t.test(x)$p.value )
p_values_twoside <- setNames(p_values_twoside, labels)
## FDR control
fdr_twoside <- p.adjust(p_values_twoside, method = "BH")
means <- sapply(samples, mean)
positives_twoside <- fdr_twoside<0.05 & means>0

false_positives <- sum(positives_twoside[names(positives_twoside) != "increased"])
predicted_positives <- sum(positives_twoside)
true_positives <- sum(positives_twoside[names(positives_twoside) == "increased"])
positives <- sum( names(positives_twoside) == "increased" )

# false discovery rate
false_positives / predicted_positives
#> [1] 0.01219512

# true positive rate
true_positives / positives
#> [1] 0.405

Your method that filters the objects and specifies the total number of tests is also correct. You can prove it mathematically. I just don't like the filtering method. There is a loss of information, although the post's author claims there is not.