3
$\begingroup$

I'm trying to demonstrate BH method for controlling false discovery rate in R.

The problem is that I get quite unexpected results (so I assume that something is going wrong with this)

set.seed(4)
x = matrix(rnorm(100*5),nrow=100)
set.seed(5)
y = matrix(rnorm(100*5),nrow=100)
set.seed(6)
y[81:100,]= rnorm(100, mean=3)

The p-values are produces from above distributions. (So the first 80 p-values are coming from true null distribution and the rest 20 follow the alternative distribution.)

p = sapply(1:100, function(i) t.test(x[i,],y[i,])$p.val)
sort.p<-sort(p)

Then the BH procedure is defined as follows (p_i<= q*i/n) (n=100, q=0.1) :

eq<-function(){(q*(1:n))/n}

rej.pvals <- sort.p[sort.p<=eq()]     #Find all the p-values which are smaller than q*i/n

After that I found the largest i which satisfy that equation. And in this example that is p_24. So I have 24 rejected p-values. The problem is that among these 24 rejected p-values, I have 4 true null hypothesis (=false discoveries). So the FDR is 4/24=0.167 which is greater than the selected q=0.1..

> sum(match(rej.p,p)<=80)
[1] 4

Can some help me to undestand what I'm missing there? Thank you in advance

$\endgroup$
1
$\begingroup$

What you seem to be computing is not the FDR, but the FDP - the false discovery proportion.

The FDR $Q_e$ is the expected value of the FDP, i.e., the proportion of false rejections $v$ to all rejections $r$,

$$Q_e=E\left(\frac{V}{R}\right)=E\left(\frac{V}{V+S}\right)=:E\left(Q\right).$$

The Benjamini-Hochberg procedure promises to control the FDR, i.e., that the expected proportion of false rejections is less than $q$, but, unfortunately, not necessarily the proportion itself.

$\endgroup$
2
  • $\begingroup$ Hi, Thanks for your answer. But I'm still little uncertain, since in this case I know the count of true null hypothesis, so in that case would it be that FDR=FDP? $\endgroup$
    – Laura
    Jun 14 '18 at 7:37
  • $\begingroup$ No. In this case, in which you generate the data, you are in the position to verify that the r.v. (the FDP) takes a larger value than its expectation, the FDR: $0.167>0.1$. In practice you would not know that as, as you point out, you would not know how many nulls are false. That is basically the same situation in a standard hypothesis test, in which you do not know if a rejection is the right decision or a type-I error. $\endgroup$ Jun 14 '18 at 8:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.