## Control false discovery rate in multiple test ### install ``` cpan -i Statistics::Multtest ``` ### usage ```perl use Statistics::Multtest qw(bonferroni holm hommel hochberg BH BY qvalue); use Statistics::Multtest qw(:all); use strict; my $p; # p-values can be stored in an array by reference $p = [0.01, 0.02, 0.05,0.41,0.16,0.51]; # @$res has the same order as @$p my $res = BH($p); print join "\n", @$res; # p-values can also be stored in a hash by reference $p = {"a" => 0.01, "b" => 0.02, "c" => 0.05, "d" => 0.41, "e" => 0.16, "f" => 0.51 }; # $res is also a hash reference which is the same as $p $res = holm($p); foreach (sort {$res->{a} <=> $res->{$b}} keys %$res) { print "$_ => $res->{$_}\n"; } # since qvalue does not always run successfully, # it should be embeded in 'eval' $res = eval 'qvalue($p)'; if($@) { print $@; } else { print join "\n", @$res; } ``` ### description For statistical test, p-value is the probability of false positives. While there are many hypothesis for testing simultaneously, the probability of getting at least one false positive would be large. Therefore the origin p-values should be adjusted to decrease the false discovery rate. Seven procedures to controlling false positive rates is provided. The names of the methods are derived from `p.adjust` in `stat` package and `qvalue` in `qvalue` package (http://www.bioconductor.org/packages/release/bioc/html/qvalue.html) in R. Code is translated directly from R to Perl using `List::Vectorize` module. All seven subroutines receive one argument which can either be an array reference or a hash reference, and return the adjusted p-values in corresponding data structure. The order of items in the array does not change after the adjustment. ### subroutines - `bonferroni($pvalue)` Bonferroni single-step process. - `hommel($pvalue)` Hommel singlewise process. Hommel, G. (1988). A stagewise rejective multiple test procedure based on a modified Bonferroni test. Biometrika, 75, 383¨C386. - `holm($pvalue)` Holm step-down process. Holm, S. (1979). A simple sequentially rejective multiple test procedure. Scandinavian Journal of Statistics, 6, 65¨C70. - `hochberg($pvalue)` Hochberg step-up process. Hochberg, Y. (1988). A sharper Bonferroni procedure for multiple tests of significance. Biometrika, 75, 800¨C803. - `BH($pvalue)` Benjamini and Hochberg, controlling the FDR. Benjamini, Y., and Hochberg, Y. (1995). Controlling the false discovery rate: a practical and powerful approach to multiple testing. Journal of the Royal Statistical Society Series B, 57, 289¨C300. - `BY($pvalue)` Use Benjamini and Yekutieli. Benjamini, Y., and Yekutieli, D. (2001). The control of the false discovery rate in multiple testing under dependency. Annals of Statistics 29, 1165¨C1188. - `qvalue($pvalue, %setup)` Storey and Tibshirani. Storey JD and Tibshirani R. (2003) Statistical significance for genome-wide experiments. Proceedings of the National Academy of Sciences, 100: 9440-9445. The default method for estimating `pi0` in the origin `qvalue` package is to utilize cubic spline interpolation. However, there is no suitable perl module to do such work (external libraries should be installed if using `Math::GSL::Spline` and there seems to be some mistakes when I using `Math::Spline`). Therefore, in this module, we only provide 'bootstrap' method to estimate `pi0`, which is also the second `pi0` method in `qvalue` package. Some arguments which are the same in `qvalue` package can be set in `%setup` as follows. ```perl lambda => multiply(seq(0, 90, 5), # The value of the tuning parameter 0.01), # to estimate pi_0. Must be in [0,1). # It should be an array reference robust => 0, # An indicator of whether it is desired to make the estimate # more robust for small p-values and a direct finite sample # estimate of pFDR ``` For details, please see the Storey (2003) and the qvalue document in R. NOTE: The results of this subroutine are not always exactly consistent to the `qvalue` package due to the floating point number calculation. In some circumstance, the estimated `pi0 <= 0`, and the subroutine would throw an error. (try p-value list: `[0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1]`). So you should embeded this subroutine in 'eval' such as: ```perl my $qvalue; eval '$qvalue = qvalue($pvalue, %setup)'; if($@) { # do something } else { print join ", ", @$qvalue; } ```