CODEWITHSUNDEEP
rstatistics
KnightofniDK
KnightofniDK

Reputation: 169

Same code yields different results (p.adjust)

I want to perform p adjustment for Benjamini and Hochberg’s FDR using R. However when I use the build-in function the result is not what I would expect (almost exclusively the same values). However, when I extract the function (with minor modifications) and have it as its own it gives a different result. I am curious as to why that is.

Full list of p-values can be found at http://pastebin.com/w2WRJD2s


Code:

BHAdjust = function(p) {
  lp = n = length(p)
  i = lp:1L
  o = order(p, decreasing = TRUE)
  ro = order(o)
  return(pmin(1, cummin(n/i * p[o]))[ro])
}

adjustedValues = data.frame(p,
                            p.adjust(p, "BH"),
                            BHAdjust(p))
head(adjustedValues, 25)


Output:

> head(adjustedValues, 25)
            p p.adjust.p...BH.. BHAdjust.p.
1          NA                NA          NA
2  0.74611808         0.9935715  0.86179530
3  0.67478588         0.9935715  0.82143053
4  0.02903144         0.9935715  0.05407105
5  0.91896732         0.9935715  0.95715994
6  0.99468266         0.9975246  0.99601960
7  0.84588305         0.9935715  0.91863392
8  0.45523400         0.9935715  0.65221025
9  0.02105568         0.9935715  0.03951256
10 0.97237337         0.9935715  0.98574978
11 0.75585835         0.9935715  0.86649987
12 0.98392280         0.9935715  0.98906222
13 0.51172583         0.9935715  0.69568566
14 0.68731809         0.9935715  0.82806060
15 0.65301206         0.9935715  0.80637908
16 0.59183087         0.9935715  0.75312367
17 0.34104743         0.9935715  0.52065642
18 0.34438992         0.9935715  0.52468404
19         NA                NA          NA
20 0.07132209         0.9935715  0.13023274
21 0.59003164         0.9935715  0.75283716
22 0.52619364         0.9935715  0.70506163
23 0.92279478         0.9935715  0.95883139
24         NA                NA          NA
25         NA                NA          NA

Upvotes: 0

Views: 272

Answers (2)

Roland
Roland

Reputation: 132874

Here is a barebones version that handles NAs:

BHAdjust <- function (p) { 
  p0 <- p
  if (all(nna <- !is.na(p))) 
    nna <- TRUE
  p <- p[nna]
  #note that n = length(p) in p.adjust gets evaluated when N is called the first time
  n <- lp <- length(p) 
  i <- lp:1L
  o <- order(p, decreasing = TRUE)
  ro <- order(o)
  p0[nna] <- pmin(1, cummin(n/i * p[o]))[ro]

  p0
}

all.equal(
  BHAdjust(adjustedValues[,1]),
  p.adjust(adjustedValues[,1], method="BH"))
#[1] TRUE

Upvotes: 1

pushStack
pushStack

Reputation: 4333

Try to use the last extra parameter :

p.adjust(p, method = p.adjust.methods, n = length(p))

Library Adjust

Upvotes: 0

Related Questions