FDR (mafdr function) odd behavior

I noticed this odd result when using matlab's FDR function - thought maybe someone might want to comment. Maybe I am using the function badly - but I guess I might not be the only one if so - thanks.

clc;clear;
rng(0);
fp = [0,0];
N_sims = 1000;
N_samp = 20;
adj_p2 = nan(N_samp,1);adj_p1 = nan(N_samp,1);
for ix = 1:N_sims
    p = rand(N_samp,1);
    [adj_p1,Q] = mafdr(p);
%     [~, ~, ~, adj_p2] = fdr_bh(p);
    fp = fp + double([adj_p1,adj_p2]<0.05);
end
sum(fp)/N_sims

result: 3.3 (and near 0.05 as expected with version from https://www.mathworks.com/matlabcentral/fileexchange/27418-fdr-bh)

I found this which was relevant (although it's for 2012, while I am using 2018a): https://www.mathworks.com/matlabcentral/answers/44824-why-am-i-getting-a-warning-in-mafdr-when-using-my-list-of-p-values

There is a suggestion there to set lambda = 0.15. This does seem to work, infact, pretty much any value of lambda that is settable seems to give more reasonable answers than when it is not set.

It seems a bit dangerous that by default mafdr has this odd behavior if lambda is not set.

Answer 1

[Hua Xie]

For future reference, the number of p-values here is too small to use Storey's Q-value method. If you change N_samp to say 500, the results would be much more reasonable. The default FDR implementation in mafdr was originally developed for genomic data, as it is mentioned in https://www.ncbi.nlm.nih.gov/pmc/articles/PMC170937/, "Because we are considering many features (i.e., m is very large), it can be shown that ..."

For small sets of p-value (m<1,000) such as the one you are testing, it is safer to use Benjamini and Hochberg algorithm using FDR = mafdr(P,'BHFDR',1), which probably should be the default setting for mafdr.

For detailed discussion on this issue, read https://support.bioconductor.org/p/56959/