R: Unexpected Results from P.Adjust (Fdr)

R: unexpected results from p.adjust (FDR)

The method is working fine. You just don't understand what the algorithm is doing fully apparently. I would suggest reading up a little bit about the procedure (Wikipedia Link).

n <- 40
p <- runif(n, 0, 1)
adjusted.p <- p.adjust(p, "fdr") #same as p.adjust(p, "BH")

## Let's recreate.  Going back to the paper we
## can construct a q-value by sorting the p-values
## and then applying the transformation:
## q_i = p_i*n/i
## Where p_i is the ith smallest p-value

## Sort the p-values and keep track of the order
id <- order(p)
tmp <- p[id]
(q <- tmp*n/(1:n))
# [1] 2.0322183 1.0993436 1.2357457 1.1113084 0.9282536 0.7877181 0.7348436
# [8] 0.7204077 0.6587102 0.6289974 0.9205222 0.8827835 0.8600234 0.8680704
#[15] 1.1532693 1.1055951 1.0451330 1.1314681 1.1167659 1.1453142 1.1012847
#[22] 1.0783747 1.0672471 1.0500311 1.0389407 1.0089661 1.0117881 0.9917817
#[29] 0.9892826 0.9668636 0.9844052 0.9792733 1.0000320 0.9967073 0.9707149
#[36] 0.9747723 0.9968153 0.9813367 1.0058445 0.9942353

## However there is one more portion to the adjustment
## We don't want a p-value of .02 getting
## a larger q-value than a p-value of .05
## so we make sure that if a smaller q-value
## shows up in the list we set all of
## the q-values corresponding to smaller p-values
## to that corresponding q-value.

## We can do this by reversing the list
## Keeping a running tally of the minimum value
## and then re-reversing
new.q <- rev(cummin(rev(q)))

## Put it back in the original order
new.q[order(id)]
# [1] 0.6289974 0.9668636 0.6289974 0.6289974 0.6289974 0.9707149 0.9707149
# [8] 0.8600234 0.9668636 0.6289974 0.9707149 0.9668636 0.9668636 0.9813367
#[15] 0.9668636 0.9668636 0.9668636 0.9707149 0.9668636 0.9668636 0.6289974
#[22] 0.9942353 0.8600234 0.6289974 0.9668636 0.6289974 0.6289974 0.8680704
#[29] 0.9668636 0.9668636 0.9668636 0.9942353 0.9707149 0.9813367 0.9668636
#[36] 0.8600234 0.6289974 0.9668636 0.9668636 0.9747723

## This should match the adjustment
adjusted.p
# [1] 0.6289974 0.9668636 0.6289974 0.6289974 0.6289974 0.9707149 0.9707149
# [8] 0.8600234 0.9668636 0.6289974 0.9707149 0.9668636 0.9668636 0.9813367
#[15] 0.9668636 0.9668636 0.9668636 0.9707149 0.9668636 0.9668636 0.6289974
#[22] 0.9942353 0.8600234 0.6289974 0.9668636 0.6289974 0.6289974 0.8680704
#[29] 0.9668636 0.9668636 0.9668636 0.9942353 0.9707149 0.9813367 0.9668636
#[36] 0.8600234 0.6289974 0.9668636 0.9668636 0.9747723

BH correction procedure giving monotonous/repeatitive adjusted p-values

This is not a "problem" per se. Please have a look at these two posts for details on BH correction at
CrossValidated and
StackOverflow.

The data you have (assuming that it is properly preprocessed / normalized) suggests that there is no significant difference based on sex. There might be other confounding effects that you need to take into account though (batch, age, disease status, etc...).

How Does R Calculate the False Discovery Rate

The reason is that the FDR calculation ensures that FDR never increases as the p-value decreases. That's because you can always choose to set a higher threshold for your rejection rule if that higher threshold will get you a lower FDR.

In your case, your second hypothesis had a p-value of 0.0006 and an FDR of 0.010971585, but the third hypothesis had a larger p-value and a smaller FDR. If you can achieve an FDR of 0.009120228 by setting your p-value threshold to 0.0019, there is never a reason to set a lower threshold just to get a higher FDR.

You can see this in the code by typing p.adjust:

...
}, BH = {
    i <- lp:1L
    o <- order(p, decreasing = TRUE)
    ro <- order(o)
    pmin(1, cummin(n/i * p[o]))[ro]

The cummin function takes the cumulative minimum of the vector, going backwards in the order of p.

You can see this in the Benjamini-Hochberg paper you link to, including in the definition of the procedure on page 293, which states (emphasis mine):

let k be the largest i for which P(i) <= i / m q*;

then reject all H_(i) i = 1, 2, ..., k

How to find orphan plugins in eclipse RCPs?

Here's a starting point (this applies for Eclipse 3.4 and later, when the P2 repository was introduced, earlier versions store their configuration differently. IIRC you could see all the plugins and features in platform.xml).

Create a new plugin project (New->Other->Plug-in Development->Plug-in Project) with the "Hello World" template then drop this code into the run method of the SampleAction.

Run the plugin as a test Eclipse Application and select Sample Menu->Sample Action, the plugins that don't belong to a feature will be output to the parent workspace's console. When I ran this there were quite a few more than I expected, I've had a few looks through and can't spot the logic error.

Edit, found logic error, was using the wrong array index used in innermost loop. Still not quite right though.

Edit 2. (Facepalm moment) Found the problem. Be sure to run the target workspace with all workspace and enabled target plugins enabled, or it will skew your results, obviously. If you install the plugin and dress it up a little bit you won't have this problem.

//get all the plugins that belong to features
IBundleGroupProvider[] providers = Platform.getBundleGroupProviders();

Map<Long, IBundleGroup> bundlesMap = new HashMap<Long, IBundleGroup>();

if (providers != null) {
    for (int i = 0; i < providers.length; i++) {
        IBundleGroup[] bundleGroups = providers[i].getBundleGroups();

        System.out.println("Bundle groups:");
        for (int j = 0; j < bundleGroups.length; j++) {
            Bundle[] bundles = bundleGroups[j] == null ? new Bundle[0] : bundleGroups[j]
                    .getBundles();
            System.out.println(bundleGroups[j].getIdentifier());
            for (int k = 0; k < bundles.length; k++) {
                bundlesMap.put(bundles[k].getBundleId(), bundleGroups[j]);
            }                
        }
    }
}

BundleContext bundleContext = Activator.getDefault().getBundle().getBundleContext();

if(bundleContext instanceof BundleContextImpl) {            
    Bundle[] bundles = ((BundleContextImpl)bundleContext).getBundles();

    System.out.println("Orphan Bundles:");
    for (int i = 0; i < bundles.length; i++) {
        if(!bundlesMap.containsKey(bundles[i].getBundleId())) {
            System.out.println(bundles[i].getSymbolicName());
        }
    }            
}

---

Related Topics