P-value correction: Difference between revisions
Thundergnat (talk | contribs) m (→{{header|Perl 6}}: use a named variable for clarity) |
Thundergnat (talk | contribs) m (→{{header|Perl 6}}: whitespace regularization & minor simplifications) |
||
Line 1,449: | Line 1,449: | ||
sub check ( @p ) { die 'p-values must be in range 0.0 to 1.0' if @p.min < 0 or 1 < @p.max; @p } |
sub check ( @p ) { die 'p-values must be in range 0.0 to 1.0' if @p.min < 0 or 1 < @p.max; @p } |
||
multi ratchet ('up', @p){ my $m; @p[$_] min= $m, $m = @p[$_] for ^@p; @p } |
multi ratchet ( 'up', @p ) { my $m; @p[$_] min= $m, $m = @p[$_] for ^@p; @p } |
||
multi ratchet ('dn', @p){ my $m; @p[$_] max= $m, $m = @p[$_] for ^@p .reverse; @p } |
multi ratchet ( 'dn', @p ) { my $m; @p[$_] max= $m, $m = @p[$_] for ^@p .reverse; @p } |
||
sub schwartzian ( @p, &transform, :$ratchet ) { |
sub schwartzian ( @p, &transform, :$ratchet ) { |
||
my @pa = @p.map( {[$_, $++]} ).sort( -*.[0] ). |
my @pa = @p.map( {[$_, $++]} ).sort( -*.[0] ).map( {[transform(.[0]), .[1]]} ); |
||
.map( {[ transform(.[0]), .[1] ]} ).values; |
|||
@pa[*;0] = ratchet($ratchet, @pa»[0]); |
@pa[*;0] = ratchet($ratchet, @pa»[0]); |
||
@pa.sort( *.[1] )»[0] |
@pa.sort( *.[1] )»[0] |
||
Line 1,519: | Line 1,518: | ||
; |
; |
||
for < Benjamini-Hochberg Benjamini-Yekutieli Bonferroni Hochberg Holm Hommel Šidák> |
for < Benjamini-Hochberg Benjamini-Yekutieli Bonferroni Hochberg Holm Hommel Šidák > |
||
{ |
{ |
||
say adjusted @p-values, $_ |
say adjusted @p-values, $_ |
Revision as of 14:49, 11 November 2017
Given a list of p-Values, adjust the p-values for multiple comparisons. This is done in order to control the false positive, or Type 1 error rate. This is also known as the "FDR" After adjustment, the p-values will be significantly higher.
Task Description
Given one list of p-values, return the p-values correcting for multiple comparisons
p = {4.533744e-01, 7.296024e-01, 9.936026e-02, 9.079658e-02, 1.801962e-01, 8.752257e-01, 2.922222e-01, 9.115421e-01, 4.355806e-01, 5.324867e-01, 4.926798e-01, 5.802978e-01, 3.485442e-01, 7.883130e-01, 2.729308e-01, 8.502518e-01, 4.268138e-01, 6.442008e-01, 3.030266e-01, 5.001555e-02, 3.194810e-01, 7.892933e-01, 9.991834e-01, 1.745691e-01, 9.037516e-01, 1.198578e-01, 3.966083e-01, 1.403837e-02, 7.328671e-01, 6.793476e-02, 4.040730e-03, 3.033349e-04, 1.125147e-02, 2.375072e-02, 5.818542e-04, 3.075482e-04, 8.251272e-03, 1.356534e-03, 1.360696e-02, 3.764588e-04, 1.801145e-05, 2.504456e-07, 3.310253e-02, 9.427839e-03, 8.791153e-04, 2.177831e-04, 9.693054e-04, 6.610250e-05, 2.900813e-02, 5.735490e-03}
There are numerous implementations of how to do this, namely Benjamini-Hochberg, Benjamini-Yekutieli, Holm, Hochberg, Hommel, Bonferroni. Each of which has its own advantages and disadvantages.
C
This work is a translation of the R source code. In order to confirm that the new function is working correctly, each value is compared to R's output and a cumulative absolute error is returned.
The C function p_adjust
is designed to work as similarly to the R function p.adjust
as possible, and is able to do any one of the methods.
This program, for example, fdr.c, can be compiled by
gcc -o fdr fdr.c -Wall -pedantic -std=c11 -lm -O4
or
clang -o fdr fdr.c -Wall -pedantic -std=c11 -lm -O4
.
Link with -lm
<lang C>#include <stdio.h>//printf
- include <stdlib.h>//qsort
- include <math.h>//fabs
- include <stdbool.h>//bool data type
- include <strings.h>//strcasecmp
unsigned int *restrict seq_len(const size_t START, const size_t END) { //named after R function of same name, but simpler function size_t start = START; size_t end = END; if (START == END) { unsigned int *restrict sequence = malloc( (end+1) * sizeof(unsigned int)); if (sequence == NULL) { printf("malloc failed at %s line %u\n", __FILE__, __LINE__); perror(""); exit(EXIT_FAILURE); } for (size_t i = 0; i < end; i++) { sequence[i] = i+1; } return sequence; } if (START > END) { end = START; start = END; } const size_t LENGTH = end - start ; unsigned int *restrict sequence = malloc( (1+LENGTH) * sizeof(unsigned int)); if (sequence == NULL) { printf("malloc failed at %s line %u\n", __FILE__, __LINE__); perror(""); exit(EXIT_FAILURE); } if (START < END) { for (size_t index = 0; index <= LENGTH; index++) { sequence[index] = start + index; } } else { for (size_t index = 0; index <= LENGTH; index++) { sequence[index] = end - index; } } return sequence; }
//modified from https://phoxis.org/2012/07/12/get-sorted-index-orderting-of-an-array/
double *restrict base_arr = NULL;
static int compar_increase (const void *restrict a, const void *restrict b) { int aa = *((int *restrict ) a), bb = *((int *restrict) b); if (base_arr[aa] < base_arr[bb]) { return 1; } else if (base_arr[aa] == base_arr[bb]) { return 0; } else { return -1; } }
static int compar_decrease (const void *restrict a, const void *restrict b) { int aa = *((int *restrict ) a), bb = *((int *restrict) b); if (base_arr[aa] < base_arr[bb]) { return -1; } else if (base_arr[aa] == base_arr[bb]) { return 0; } else { return 1; } }
unsigned int *restrict order (const double *restrict ARRAY, const unsigned int SIZE, const bool DECREASING) { //this has the same name as the same R function unsigned int *restrict idx = malloc(SIZE * sizeof(unsigned int)); if (idx == NULL) { printf("failed to malloc at %s line %u.\n", __FILE__, __LINE__); perror(""); exit(EXIT_FAILURE); } base_arr = malloc(sizeof(double) * SIZE); if (base_arr == NULL) { printf("failed to malloc at %s line %u.\n", __FILE__, __LINE__); perror(""); exit(EXIT_FAILURE); } for (unsigned int i = 0; i < SIZE; i++) { base_arr[i] = ARRAY[i]; idx[i] = i; } if (DECREASING == false) { qsort(idx, SIZE, sizeof(unsigned int), compar_decrease); } else if (DECREASING == true) { qsort(idx, SIZE, sizeof(unsigned int), compar_increase); } free(base_arr); base_arr = NULL; return idx; }
double *restrict cummin(const double *restrict ARRAY, const unsigned int NO_OF_ARRAY_ELEMENTS) { //this takes the same name of the R function which it copies //this requires a free() afterward where it is used if (NO_OF_ARRAY_ELEMENTS < 1) { puts("cummin function requires at least one element.\n"); printf("Failed at %s line %u\n", __FILE__, __LINE__); exit(EXIT_FAILURE); } double *restrict output = malloc(sizeof(double) * NO_OF_ARRAY_ELEMENTS); if (output == NULL) { printf("failed to malloc at %s line %u.\n", __FILE__, __LINE__); perror(""); exit(EXIT_FAILURE); } double cumulative_min = ARRAY[0]; for (unsigned int i = 0; i < NO_OF_ARRAY_ELEMENTS; i++) { if (ARRAY[i] < cumulative_min) { cumulative_min = ARRAY[i]; } output[i] = cumulative_min; } return output; }
double *restrict cummax(const double *restrict ARRAY, const unsigned int NO_OF_ARRAY_ELEMENTS) { //this takes the same name of the R function which it copies //this requires a free() afterward where it is used if (NO_OF_ARRAY_ELEMENTS < 1) { puts("function requires at least one element.\n"); printf("Failed at %s line %u\n", __FILE__, __LINE__); exit(EXIT_FAILURE); } double *restrict output = malloc(sizeof(double) * NO_OF_ARRAY_ELEMENTS); if (output == NULL) { printf("failed to malloc at %s line %u.\n", __FILE__, __LINE__); perror(""); exit(EXIT_FAILURE); } double cumulative_max = ARRAY[0]; for (size_t i = 0; i < NO_OF_ARRAY_ELEMENTS; i++) { if (ARRAY[i] > cumulative_max) { cumulative_max = ARRAY[i]; } output[i] = cumulative_max; } return output; }
double *restrict pminx(const double *restrict ARRAY, const size_t NO_OF_ARRAY_ELEMENTS, const double X) { //named after the R function pmin if (NO_OF_ARRAY_ELEMENTS < 1) { puts("pmin requires at least one element.\n"); printf("Failed at %s line %u\n", __FILE__, __LINE__); exit(EXIT_FAILURE); } double *restrict pmin_array = malloc(sizeof(double) * NO_OF_ARRAY_ELEMENTS); if (pmin_array == NULL) { printf("failed to malloc at %s line %u.\n", __FILE__, __LINE__); perror(""); exit(EXIT_FAILURE); } for (unsigned int index = 0; index < NO_OF_ARRAY_ELEMENTS; index++) { if (ARRAY[index] < X) { pmin_array[index] = ARRAY[index]; } else { pmin_array[index] = X; } } return pmin_array; }
void double_say (const double *restrict ARRAY, const size_t NO_OF_ARRAY_ELEMENTS) { printf("[1] %e", ARRAY[0]); for (size_t i = 1; i < NO_OF_ARRAY_ELEMENTS; i++) { printf(" %.10f", ARRAY[i]); if (((i+1) % 5) == 0) { printf("\n[%zu]", i+1); } } puts("\n"); }
/*void uint_say (const unsigned int *restrict ARRAY, const size_t NO_OF_ARRAY_ELEMENTS) { //for debugging printf("%u", ARRAY[0]); for (size_t i = 1; i < NO_OF_ARRAY_ELEMENTS; i++) { printf(",%u", ARRAY[i]); } puts("\n"); }*/
double *restrict uint2double (const unsigned int *restrict ARRAY, const size_t NO_OF_ARRAY_ELEMENTS) { double *restrict doubleArray = malloc(sizeof(double) * NO_OF_ARRAY_ELEMENTS); if (doubleArray == NULL) { printf("Failure to malloc at %s line %u.\n", __FILE__, __LINE__); perror(""); exit(EXIT_FAILURE); } for (size_t index = 0; index < NO_OF_ARRAY_ELEMENTS; index++) { doubleArray[index] = (double)ARRAY[index]; } return doubleArray; }
double min2 (const double N1, const double N2) { if (N1 < N2) { return N1; } else { return N2; } }
double *restrict p_adjust (const double *restrict PVALUES, const size_t NO_OF_ARRAY_ELEMENTS, const char *restrict STRING) { //this function is a translation of R's p.adjust "BH" method // i is always i[index] = NO_OF_ARRAY_ELEMENTS - index - 1 if (NO_OF_ARRAY_ELEMENTS < 1) { puts("p_adjust requires at least one element.\n"); printf("Failed at %s line %u\n", __FILE__, __LINE__); exit(EXIT_FAILURE); } short int TYPE = -1; if (strcasecmp(STRING, "BH") == 0) { TYPE = 0; } else if (strcasecmp(STRING, "fdr") == 0) { TYPE = 0; } else if (strcasecmp(STRING, "by") == 0) { TYPE = 1; } else if (strcasecmp(STRING, "Bonferroni") == 0) { TYPE = 2; } else if (strcasecmp(STRING, "hochberg") == 0) { TYPE = 3; } else if (strcasecmp(STRING, "holm") == 0) { TYPE = 4; } else if (strcasecmp(STRING, "hommel") == 0) { TYPE = 5; } else { printf("%s doesn't match any accepted FDR methods.\n", STRING); printf("Failed at %s line %u\n", __FILE__, __LINE__); exit(EXIT_FAILURE); } //--------------------------------------------------------------------------- //--------------------------------------------------------------------------- if (TYPE == 2) {//Bonferroni method double *restrict bonferroni = malloc(sizeof(double) * NO_OF_ARRAY_ELEMENTS); if (bonferroni == NULL) { printf("failed to malloc at %s line %u.\n", __FILE__, __LINE__); perror(""); exit(EXIT_FAILURE); } for (size_t index = 0; index < NO_OF_ARRAY_ELEMENTS; index++) { const double BONFERRONI = PVALUES[index] * NO_OF_ARRAY_ELEMENTS; if (BONFERRONI >= 1.0) { bonferroni[index] = 1.0; } else if ((0.0 <= BONFERRONI) && (BONFERRONI < 1.0)) { bonferroni[index] = BONFERRONI; } else { printf("%g is outside of the interval I planned.\n", BONFERRONI); printf("Failure at %s line %u\n", __FILE__, __LINE__); exit(EXIT_FAILURE); } } return bonferroni; //--------------------------------------------------------------------------- //--------------------------------------------------------------------------- } else if (TYPE == 4) {//Holm method /*these values are computed separately from BH, BY, and Hochberg because they are computed differently*/ unsigned int *restrict o = order(PVALUES, NO_OF_ARRAY_ELEMENTS, false); //sorted in reverse of methods 0-3 double *restrict o2double = uint2double(o, NO_OF_ARRAY_ELEMENTS); double *restrict cummax_input = malloc(sizeof(double) * NO_OF_ARRAY_ELEMENTS); for (size_t index = 0; index < NO_OF_ARRAY_ELEMENTS; index++) { cummax_input[index] = (NO_OF_ARRAY_ELEMENTS - index ) * PVALUES[o[index]]; // printf("cummax_input[%zu] = %e\n", index, cummax_input[index]); } free(o); o = NULL; unsigned int *restrict ro = order(o2double, NO_OF_ARRAY_ELEMENTS, false); free(o2double); o2double = NULL;
double *restrict cummax_output = cummax(cummax_input, NO_OF_ARRAY_ELEMENTS); free(cummax_input); cummax_input = NULL;
double *restrict pmin = pminx(cummax_output, NO_OF_ARRAY_ELEMENTS, 1); free(cummax_output); cummax_output = NULL; double *restrict qvalues = malloc(sizeof(double) * NO_OF_ARRAY_ELEMENTS); for (size_t index = 0; index < NO_OF_ARRAY_ELEMENTS; index++) { qvalues[index] = pmin[ro[index]]; } free(pmin); pmin = NULL; free(ro); ro = NULL; return qvalues; //--------------------------------------------------------------------------- //--------------------------------------------------------------------------- } else if (TYPE == 5) {//Hommel method //i <- seq_len(n) //o <- order(p) unsigned int *restrict o = order(PVALUES, NO_OF_ARRAY_ELEMENTS, false);//false is R's default //p <- p[o] double *restrict p = malloc(sizeof(double) * NO_OF_ARRAY_ELEMENTS); if (p == NULL) { printf("failed to malloc at %s line %u.\n", __FILE__, __LINE__); perror(""); exit(EXIT_FAILURE); } for (size_t index = 0; index < NO_OF_ARRAY_ELEMENTS; index++) { p[index] = PVALUES[o[index]]; } //ro <- order(o) double *restrict o2double = uint2double(o, NO_OF_ARRAY_ELEMENTS); free(o); o = NULL; unsigned int *restrict ro = order(o2double, NO_OF_ARRAY_ELEMENTS, false); free(o2double); o2double = NULL; // puts("ro"); //q <- pa <- rep.int(min(n * p/i), n) double *restrict q = malloc(sizeof(double) * NO_OF_ARRAY_ELEMENTS); if (q == NULL) { printf("failed to malloc at %s line %u.\n", __FILE__, __LINE__); perror(""); exit(EXIT_FAILURE); } double *restrict pa = malloc(sizeof(double) * NO_OF_ARRAY_ELEMENTS); if (pa == NULL) { printf("failed to malloc at %s line %u.\n", __FILE__, __LINE__); perror(""); exit(EXIT_FAILURE); } double min = (double)NO_OF_ARRAY_ELEMENTS * p[0]; for (size_t index = 1; index < NO_OF_ARRAY_ELEMENTS; index++) { const double TEMP = (double)NO_OF_ARRAY_ELEMENTS * p[index] / (1+index); if (TEMP < min) { min = TEMP; } } for (size_t index = 0; index < NO_OF_ARRAY_ELEMENTS; index++) { pa[index] = min; q[index] = min; } // puts("q & pa"); // double_say(q, NO_OF_ARRAY_ELEMENTS); /*for (j in (n - 1):2) {
ij <- seq_len(n - j + 1) i2 <- (n - j + 2):n q1 <- min(j * p[i2]/(2:j)) q[ij] <- pmin(j * p[ij], q1) q[i2] <- q[n - j + 1] pa <- pmax(pa, q) }
- /
for (size_t j = (NO_OF_ARRAY_ELEMENTS-1); j >= 2; j--) { // printf("j = %zu\n", j); unsigned int *restrict ij = seq_len(1,NO_OF_ARRAY_ELEMENTS - j + 1); for (size_t i = 0; i < NO_OF_ARRAY_ELEMENTS - j + 1; i++) { ij[i]--;//R's indices are 1-based, C's are 0-based } const size_t I2_LENGTH = j - 1; unsigned int *restrict i2 = malloc(I2_LENGTH * sizeof(unsigned int)); for (size_t i = 0; i < I2_LENGTH; i++) { i2[i] = NO_OF_ARRAY_ELEMENTS-j+2+i-1; //R's indices are 1-based, C's are 0-based, I added the -1 }
double q1 = j * p[i2[0]] / 2.0; for (size_t i = 1; i < I2_LENGTH; i++) {//loop through 2:j const double TEMP_Q1 = (double)j * p[i2[i]] / (2 + i); if (TEMP_Q1 < q1) { q1 = TEMP_Q1; } }
for (size_t i = 0; i < (NO_OF_ARRAY_ELEMENTS - j + 1); i++) {//q[ij] <- pmin(j * p[ij], q1) q[ij[i]] = min2( j*p[ij[i]], q1); } free(ij); ij = NULL;
for (size_t i = 0; i < I2_LENGTH; i++) {//q[i2] <- q[n - j + 1] q[i2[i]] = q[NO_OF_ARRAY_ELEMENTS - j];//subtract 1 because of starting index difference } free(i2); i2 = NULL;
for (size_t i = 0; i < NO_OF_ARRAY_ELEMENTS; i++) {//pa <- pmax(pa, q) if (pa[i] < q[i]) { pa[i] = q[i]; } } // printf("j = %zu, pa = \n", j); // double_say(pa, N); }//end j loop free(p); p = NULL; for (size_t index = 0; index < NO_OF_ARRAY_ELEMENTS; index++) { q[index] = pa[ro[index]];//Hommel q-values } //now free memory free(ro); ro = NULL; free(pa); pa = NULL; return q; } //The methods are similarly computed and thus can be combined for clarity unsigned int *restrict o = order(PVALUES, NO_OF_ARRAY_ELEMENTS, true); if (o == NULL) { printf("failed to malloc at %s line %u.\n", __FILE__, __LINE__); perror(""); exit(EXIT_FAILURE); } double *restrict o_double = uint2double(o, NO_OF_ARRAY_ELEMENTS); for (unsigned int index = 0; index < NO_OF_ARRAY_ELEMENTS; index++) { if ((PVALUES[index] < 0) || (PVALUES[index] > 1)) { printf("array[%u] = %lf, which is outside the interval [0,1]\n", index, PVALUES[index]); printf("died at %s line %u\n", __FILE__, __LINE__); exit(EXIT_FAILURE); } }
unsigned int *restrict ro = order(o_double, NO_OF_ARRAY_ELEMENTS, false); if (ro == NULL) { printf("failed to malloc at %s line %u.\n", __FILE__, __LINE__); perror(""); exit(EXIT_FAILURE); } free(o_double); o_double = NULL; double *restrict cummin_input = malloc(sizeof(double) * NO_OF_ARRAY_ELEMENTS); if (TYPE == 0) {//BH method for (size_t index = 0; index < NO_OF_ARRAY_ELEMENTS; index++) { const double NI = (double)NO_OF_ARRAY_ELEMENTS / (NO_OF_ARRAY_ELEMENTS - index);// n/i simplified cummin_input[index] = NI * PVALUES[o[index]];//PVALUES[o[index]] is p[o] } } else if (TYPE == 1) {//BY method double q = 1.0; for (size_t index = 2; index < (1+NO_OF_ARRAY_ELEMENTS); index++) { q += (double) 1.0/index; } for (size_t index = 0; index < NO_OF_ARRAY_ELEMENTS; index++) { const double NI = (double)NO_OF_ARRAY_ELEMENTS / (NO_OF_ARRAY_ELEMENTS - index);// n/i simplified cummin_input[index] = q * NI * PVALUES[o[index]];//PVALUES[o[index]] is p[o] } } else if (TYPE == 3) {//Hochberg method for (size_t index = 0; index < NO_OF_ARRAY_ELEMENTS; index++) { // pmin(1, cummin((n - i + 1L) * p[o]))[ro] cummin_input[index] = (index + 1) * PVALUES[o[index]]; } } free(o); o = NULL; double *restrict cummin_array = NULL; cummin_array = cummin(cummin_input, NO_OF_ARRAY_ELEMENTS); free(cummin_input); cummin_input = NULL;//I don't need this anymore double *restrict pmin = pminx(cummin_array, NO_OF_ARRAY_ELEMENTS, 1); free(cummin_array); cummin_array = NULL; double *restrict q_array = malloc(NO_OF_ARRAY_ELEMENTS*sizeof(double)); for (unsigned int index = 0; index < NO_OF_ARRAY_ELEMENTS; index++) { q_array[index] = pmin[ro[index]]; }
free(ro); ro = NULL; free(pmin); pmin = NULL; return q_array; }
int main(void) { const double PVALUES[] = {4.533744e-01, 7.296024e-01, 9.936026e-02, 9.079658e-02, 1.801962e-01, 8.752257e-01, 2.922222e-01, 9.115421e-01, 4.355806e-01, 5.324867e-01, 4.926798e-01, 5.802978e-01, 3.485442e-01, 7.883130e-01, 2.729308e-01, 8.502518e-01, 4.268138e-01, 6.442008e-01, 3.030266e-01, 5.001555e-02, 3.194810e-01, 7.892933e-01, 9.991834e-01, 1.745691e-01, 9.037516e-01, 1.198578e-01, 3.966083e-01, 1.403837e-02, 7.328671e-01, 6.793476e-02, 4.040730e-03, 3.033349e-04, 1.125147e-02, 2.375072e-02, 5.818542e-04, 3.075482e-04, 8.251272e-03, 1.356534e-03, 1.360696e-02, 3.764588e-04, 1.801145e-05, 2.504456e-07, 3.310253e-02, 9.427839e-03, 8.791153e-04, 2.177831e-04, 9.693054e-04, 6.610250e-05, 2.900813e-02, 5.735490e-03};//just the pvalues const double CORRECT_ANSWERS[6][50] = {//each first index is type {6.126681e-01, 8.521710e-01, 1.987205e-01, 1.891595e-01, 3.217789e-01, 9.301450e-01, 4.870370e-01, 9.301450e-01, 6.049731e-01, 6.826753e-01, 6.482629e-01, 7.253722e-01, 5.280973e-01, 8.769926e-01, 4.705703e-01, 9.241867e-01, 6.049731e-01, 7.856107e-01, 4.887526e-01, 1.136717e-01, 4.991891e-01, 8.769926e-01, 9.991834e-01, 3.217789e-01, 9.301450e-01, 2.304958e-01, 5.832475e-01, 3.899547e-02, 8.521710e-01, 1.476843e-01, 1.683638e-02, 2.562902e-03, 3.516084e-02, 6.250189e-02, 3.636589e-03, 2.562902e-03, 2.946883e-02, 6.166064e-03, 3.899547e-02, 2.688991e-03, 4.502862e-04, 1.252228e-05, 7.881555e-02, 3.142613e-02, 4.846527e-03, 2.562902e-03, 4.846527e-03, 1.101708e-03, 7.252032e-02, 2.205958e-02},//Benjamini-Hochberg {1.000000e+00, 1.000000e+00, 8.940844e-01, 8.510676e-01, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 5.114323e-01, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.754486e-01, 1.000000e+00, 6.644618e-01, 7.575031e-02, 1.153102e-02, 1.581959e-01, 2.812089e-01, 1.636176e-02, 1.153102e-02, 1.325863e-01, 2.774239e-02, 1.754486e-01, 1.209832e-02, 2.025930e-03, 5.634031e-05, 3.546073e-01, 1.413926e-01, 2.180552e-02, 1.153102e-02, 2.180552e-02, 4.956812e-03, 3.262838e-01, 9.925057e-02},//Benjamini & Yekutieli {1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 7.019185e-01, 1.000000e+00, 1.000000e+00, 2.020365e-01, 1.516674e-02, 5.625735e-01, 1.000000e+00, 2.909271e-02, 1.537741e-02, 4.125636e-01, 6.782670e-02, 6.803480e-01, 1.882294e-02, 9.005725e-04, 1.252228e-05, 1.000000e+00, 4.713920e-01, 4.395577e-02, 1.088915e-02, 4.846527e-02, 3.305125e-03, 1.000000e+00, 2.867745e-01},//Bonferroni {9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 4.632662e-01, 9.991834e-01, 9.991834e-01, 1.575885e-01, 1.383967e-02, 3.938014e-01, 7.600230e-01, 2.501973e-02, 1.383967e-02, 3.052971e-01, 5.426136e-02, 4.626366e-01, 1.656419e-02, 8.825610e-04, 1.252228e-05, 9.930759e-01, 3.394022e-01, 3.692284e-02, 1.023581e-02, 3.974152e-02, 3.172920e-03, 8.992520e-01, 2.179486e-01},//Hochberg {1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 4.632662e-01, 1.000000e+00, 1.000000e+00, 1.575885e-01, 1.395341e-02, 3.938014e-01, 7.600230e-01, 2.501973e-02, 1.395341e-02, 3.052971e-01, 5.426136e-02, 4.626366e-01, 1.656419e-02, 8.825610e-04, 1.252228e-05, 9.930759e-01, 3.394022e-01, 3.692284e-02, 1.023581e-02, 3.974152e-02, 3.172920e-03, 8.992520e-01, 2.179486e-01},//Holm { 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.987624e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.595180e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 4.351895e-01, 9.991834e-01, 9.766522e-01, 1.414256e-01, 1.304340e-02, 3.530937e-01, 6.887709e-01, 2.385602e-02, 1.322457e-02, 2.722920e-01, 5.426136e-02, 4.218158e-01, 1.581127e-02, 8.825610e-04, 1.252228e-05, 8.743649e-01, 3.016908e-01, 3.516461e-02, 9.582456e-03, 3.877222e-02, 3.172920e-03, 8.122276e-01, 1.950067e-01}//Hommel }; //the following loop checks each type with R's answers const char *restrict TYPES[] = {"bh", "by", "bonferroni", "hochberg", "holm", "hommel"}; for (unsigned short int type = 0; type <= 5; type++) { double *restrict q = p_adjust(PVALUES, sizeof(PVALUES) / sizeof(*PVALUES), TYPES[type]); double error = fabs(q[0] - CORRECT_ANSWERS[type][0]); // printf("%e - %e = %g\n", q[0], CORRECT_ANSWERS[type][0], error); // puts("p q"); // printf("%g\t%g\n", pvalues[0], q[0]); for (unsigned int i = 1; i < sizeof(PVALUES) / sizeof(*PVALUES); i++) { const double this_error = fabs(q[i] - CORRECT_ANSWERS[type][i]); // printf("%e - %e = %g\n", q[i], CORRECT_ANSWERS[type][i], error); error += this_error; } double_say(q, sizeof(PVALUES) / sizeof(*PVALUES)); free(q); q = NULL; printf("\ntype %u = '%s' has cumulative error of %g\n", type, TYPES[type], error); }
return 0; } </lang>
- Output:
[1] 6.126681e-01 0.8521710465 0.1987205200 0.1891595417 0.3217789286 [5] 0.9301450000 0.4870370000 0.9301450000 0.6049730556 0.6826752564 [10] 0.6482628947 0.7253722500 0.5280972727 0.8769925556 0.4705703448 [15] 0.9241867391 0.6049730556 0.7856107317 0.4887525806 0.1136717045 [20] 0.4991890625 0.8769925556 0.9991834000 0.3217789286 0.9301450000 [25] 0.2304957692 0.5832475000 0.0389954722 0.8521710465 0.1476842609 [30] 0.0168363750 0.0025629017 0.0351608437 0.0625018947 0.0036365888 [35] 0.0025629017 0.0294688286 0.0061660636 0.0389954722 0.0026889914 [40] 0.0004502862 0.0000125223 0.0788155476 0.0314261300 0.0048465270 [45] 0.0025629017 0.0048465270 0.0011017083 0.0725203250 0.0220595769 [50] type 0 = 'bh' has cumulative error of 8.03053e-07 [1] 1.000000e+00 1.0000000000 0.8940844244 0.8510676197 1.0000000000 [5] 1.0000000000 1.0000000000 1.0000000000 1.0000000000 1.0000000000 [10] 1.0000000000 1.0000000000 1.0000000000 1.0000000000 1.0000000000 [15] 1.0000000000 1.0000000000 1.0000000000 1.0000000000 0.5114323399 [20] 1.0000000000 1.0000000000 1.0000000000 1.0000000000 1.0000000000 [25] 1.0000000000 1.0000000000 0.1754486368 1.0000000000 0.6644618149 [30] 0.0757503083 0.0115310209 0.1581958559 0.2812088585 0.0163617595 [35] 0.0115310209 0.1325863108 0.0277423864 0.1754486368 0.0120983246 [40] 0.0020259303 0.0000563403 0.3546073326 0.1413926119 0.0218055202 [45] 0.0115310209 0.0218055202 0.0049568120 0.3262838334 0.0992505663 [50] type 1 = 'by' has cumulative error of 3.64072e-07 [1] 1.000000e+00 1.0000000000 1.0000000000 1.0000000000 1.0000000000 [5] 1.0000000000 1.0000000000 1.0000000000 1.0000000000 1.0000000000 [10] 1.0000000000 1.0000000000 1.0000000000 1.0000000000 1.0000000000 [15] 1.0000000000 1.0000000000 1.0000000000 1.0000000000 1.0000000000 [20] 1.0000000000 1.0000000000 1.0000000000 1.0000000000 1.0000000000 [25] 1.0000000000 1.0000000000 0.7019185000 1.0000000000 1.0000000000 [30] 0.2020365000 0.0151667450 0.5625735000 1.0000000000 0.0290927100 [35] 0.0153774100 0.4125636000 0.0678267000 0.6803480000 0.0188229400 [40] 0.0009005725 0.0000125223 1.0000000000 0.4713919500 0.0439557650 [45] 0.0108891550 0.0484652700 0.0033051250 1.0000000000 0.2867745000 [50] type 2 = 'bonferroni' has cumulative error of 6.5e-08 [1] 9.991834e-01 0.9991834000 0.9991834000 0.9991834000 0.9991834000 [5] 0.9991834000 0.9991834000 0.9991834000 0.9991834000 0.9991834000 [10] 0.9991834000 0.9991834000 0.9991834000 0.9991834000 0.9991834000 [15] 0.9991834000 0.9991834000 0.9991834000 0.9991834000 0.9991834000 [20] 0.9991834000 0.9991834000 0.9991834000 0.9991834000 0.9991834000 [25] 0.9991834000 0.9991834000 0.4632662100 0.9991834000 0.9991834000 [30] 0.1575884700 0.0138396690 0.3938014500 0.7600230400 0.0250197306 [35] 0.0138396690 0.3052970640 0.0542613600 0.4626366400 0.0165641872 [40] 0.0008825610 0.0000125223 0.9930759000 0.3394022040 0.0369228426 [45] 0.0102358057 0.0397415214 0.0031729200 0.8992520300 0.2179486200 [50] type 3 = 'hochberg' has cumulative error of 2.7375e-07 [1] 1.000000e+00 1.0000000000 1.0000000000 1.0000000000 1.0000000000 [5] 1.0000000000 1.0000000000 1.0000000000 1.0000000000 1.0000000000 [10] 1.0000000000 1.0000000000 1.0000000000 1.0000000000 1.0000000000 [15] 1.0000000000 1.0000000000 1.0000000000 1.0000000000 1.0000000000 [20] 1.0000000000 1.0000000000 1.0000000000 1.0000000000 1.0000000000 [25] 1.0000000000 1.0000000000 0.4632662100 1.0000000000 1.0000000000 [30] 0.1575884700 0.0139534054 0.3938014500 0.7600230400 0.0250197306 [35] 0.0139534054 0.3052970640 0.0542613600 0.4626366400 0.0165641872 [40] 0.0008825610 0.0000125223 0.9930759000 0.3394022040 0.0369228426 [45] 0.0102358057 0.0397415214 0.0031729200 0.8992520300 0.2179486200 [50] type 4 = 'holm' has cumulative error of 2.8095e-07 [1] 9.991834e-01 0.9991834000 0.9991834000 0.9987623800 0.9991834000 [5] 0.9991834000 0.9991834000 0.9991834000 0.9991834000 0.9991834000 [10] 0.9991834000 0.9991834000 0.9991834000 0.9991834000 0.9991834000 [15] 0.9991834000 0.9991834000 0.9991834000 0.9991834000 0.9595180000 [20] 0.9991834000 0.9991834000 0.9991834000 0.9991834000 0.9991834000 [25] 0.9991834000 0.9991834000 0.4351894700 0.9991834000 0.9766522500 [30] 0.1414255500 0.0130434007 0.3530936533 0.6887708800 0.0238560222 [35] 0.0132245726 0.2722919760 0.0542613600 0.4218157600 0.0158112696 [40] 0.0008825610 0.0000125223 0.8743649143 0.3016908480 0.0351646120 [45] 0.0095824564 0.0387722160 0.0031729200 0.8122276400 0.1950066600 [50] type 5 = 'hommel' has cumulative error of 4.35302e-07
Kotlin
<lang scala>// version 1.1.51
import java.util.Arrays
typealias IAE = IllegalArgumentException
fun seqLen(start: Int, end: Int) =
when { start == end -> IntArray(end + 1) { it + 1 } start < end -> IntArray(end - start + 1) { start + it } else -> IntArray(start - end + 1) { start - it } }
var baseArr: DoubleArray? = null
fun compareIncrease(a: Int, b: Int): Int = baseArr!![b].compareTo(baseArr!![a])
fun compareDecrease(a: Int, b: Int): Int = baseArr!![a].compareTo(baseArr!![b])
fun order(array: DoubleArray, decreasing: Boolean): IntArray {
val size = array.size var idx = IntArray(size) { it } baseArr = array.copyOf() if (!decreasing) { idx = Arrays.stream(idx) .boxed() .sorted { a, b -> compareDecrease(a, b) } .mapToInt { it } .toArray() } else { idx = Arrays.stream(idx) .boxed() .sorted { a, b -> compareIncrease(a, b) } .mapToInt { it } .toArray() } baseArr = null return idx
}
fun cummin(array: DoubleArray): DoubleArray {
val size = array.size if (size < 1) throw IAE("cummin requires at least one element") val output = DoubleArray(size) var cumulativeMin = array[0] for (i in 0 until size) { if (array[i] < cumulativeMin) cumulativeMin = array[i] output[i] = cumulativeMin } return output
}
fun cummax(array: DoubleArray): DoubleArray {
val size = array.size if (size < 1) throw IAE("cummax requires at least one element") val output = DoubleArray(size) var cumulativeMax = array[0] for (i in 0 until size) { if (array[i] > cumulativeMax) cumulativeMax = array[i] output[i] = cumulativeMax } return output
}
fun pminx(array: DoubleArray, x: Double): DoubleArray {
val size = array.size if (size < 1) throw IAE("pmin requires at least one element") return DoubleArray(size) { if (array[it] < x) array[it] else x }
}
fun doubleSay(array: DoubleArray) {
print("[ 1] %e".format(array[0])) for (i in 1 until array.size) { print(" %.10f".format(array[i])) if ((i + 1) % 5 == 0) print("\n[%2d]".format(i + 1)) } println()
}
fun intToDouble(array: IntArray) = DoubleArray(array.size) { array[it].toDouble() }
fun doubleArrayMin(array: DoubleArray) =
if (array.size < 1) throw IAE("pAdjust requires at least one element") else array.min()!!
fun pAdjust(pvalues: DoubleArray, str: String): DoubleArray {
val size = pvalues.size if (size < 1) throw IAE("pAdjust requires at least one element") val type = when(str.toLowerCase()) { "bh", "fdr" -> 0 "by" -> 1 "bonferroni" -> 2 "hochberg" -> 3 "holm" -> 4 "hommel" -> 5 else -> throw IAE("'$str' doesn't match any accepted FDR types") } if (type == 2) { // Bonferroni method return DoubleArray(size) { val b = pvalues[it] * size when { b >= 1 -> 1.0 0 <= b && b < 1 -> b else -> throw RuntimeException("$b is outside [0, 1)") } } } else if (type == 4) { // Holm method val o = order(pvalues, false) val o2Double = intToDouble(o) val cummaxInput = DoubleArray(size) { (size - it) * pvalues[o[it]] } val ro = order(o2Double, false) val cummaxOutput = cummax(cummaxInput) val pmin = pminx(cummaxOutput, 1.0) return DoubleArray(size) { pmin[ro[it]] } } else if (type == 5) { // Hommel method val indices = seqLen(size, size) val o = order(pvalues, false) val p = DoubleArray(size) { pvalues[o[it]] } val o2Double = intToDouble(o) val ro = order(o2Double, false) val q = DoubleArray(size) val pa = DoubleArray(size) val npi = DoubleArray(size) { p[it] * size / indices[it] } val min = doubleArrayMin(npi) q.fill(min) pa.fill(min) for (j in size - 1 downTo 2) { val ij = seqLen(1, size - j + 1) for (i in 0 until size - j + 1) ij[i]-- val i2Length = j - 1 val i2 = IntArray(i2Length) { size - j + 2 + it - 1 } val pi2Length = i2Length var q1 = j * p[i2[0]] / 2.0 for (i in 1 until pi2Length) { val temp_q1 = p[i2[i]] * j / (2.0 + i) if(temp_q1 < q1) q1 = temp_q1 } for (i in 0 until size - j + 1) { q[ij[i]] = minOf(p[ij[i]] * j, q1) } for (i in 0 until i2Length) q[i2[i]] = q[size - j] for (i in 0 until size) if (pa[i] < q[i]) pa[i] = q[i] } for (index in 0 until size) q[index] = pa[ro[index]] return q } val ni = DoubleArray(size) val o = order(pvalues, true) val oDouble = intToDouble(o) for (index in 0 until size) { if (pvalues[index] !in 0.0 .. 1.0) { throw RuntimeException("array[$index] = ${pvalues[index]} is outside [0, 1]") } ni[index] = size.toDouble() / (size - index) } val ro = order(oDouble, false) val cumminInput = DoubleArray(size) if (type == 0) { // BH method for (index in 0 until size) { cumminInput[index] = ni[index] * pvalues[o[index]] } } else if (type == 1) { // BY method var q = 0.0 for (index in 1 until size + 1) q += 1.0 / index for (index in 0 until size) { cumminInput[index] = q * ni[index] * pvalues[o[index]] } } else if (type == 3) { // Hochberg method for (index in 0 until size) { cumminInput[index] = (index + 1) * pvalues[o[index]] } } val cumminArray = cummin(cumminInput) val pmin = pminx(cumminArray, 1.0) return DoubleArray(size) { pmin[ro[it]] }
}
fun main(args: Array<String>) {
val pvalues = doubleArrayOf( 4.533744e-01, 7.296024e-01, 9.936026e-02, 9.079658e-02, 1.801962e-01, 8.752257e-01, 2.922222e-01, 9.115421e-01, 4.355806e-01, 5.324867e-01, 4.926798e-01, 5.802978e-01, 3.485442e-01, 7.883130e-01, 2.729308e-01, 8.502518e-01, 4.268138e-01, 6.442008e-01, 3.030266e-01, 5.001555e-02, 3.194810e-01, 7.892933e-01, 9.991834e-01, 1.745691e-01, 9.037516e-01, 1.198578e-01, 3.966083e-01, 1.403837e-02, 7.328671e-01, 6.793476e-02, 4.040730e-03, 3.033349e-04, 1.125147e-02, 2.375072e-02, 5.818542e-04, 3.075482e-04, 8.251272e-03, 1.356534e-03, 1.360696e-02, 3.764588e-04, 1.801145e-05, 2.504456e-07, 3.310253e-02, 9.427839e-03, 8.791153e-04, 2.177831e-04, 9.693054e-04, 6.610250e-05, 2.900813e-02, 5.735490e-03 )
val correctAnswers = listOf( doubleArrayOf( // Benjamini-Hochberg 6.126681e-01, 8.521710e-01, 1.987205e-01, 1.891595e-01, 3.217789e-01, 9.301450e-01, 4.870370e-01, 9.301450e-01, 6.049731e-01, 6.826753e-01, 6.482629e-01, 7.253722e-01, 5.280973e-01, 8.769926e-01, 4.705703e-01, 9.241867e-01, 6.049731e-01, 7.856107e-01, 4.887526e-01, 1.136717e-01, 4.991891e-01, 8.769926e-01, 9.991834e-01, 3.217789e-01, 9.301450e-01, 2.304958e-01, 5.832475e-01, 3.899547e-02, 8.521710e-01, 1.476843e-01, 1.683638e-02, 2.562902e-03, 3.516084e-02, 6.250189e-02, 3.636589e-03, 2.562902e-03, 2.946883e-02, 6.166064e-03, 3.899547e-02, 2.688991e-03, 4.502862e-04, 1.252228e-05, 7.881555e-02, 3.142613e-02, 4.846527e-03, 2.562902e-03, 4.846527e-03, 1.101708e-03, 7.252032e-02, 2.205958e-02 ), doubleArrayOf( // Benjamini & Yekutieli 1.000000e+00, 1.000000e+00, 8.940844e-01, 8.510676e-01, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 5.114323e-01, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.754486e-01, 1.000000e+00, 6.644618e-01, 7.575031e-02, 1.153102e-02, 1.581959e-01, 2.812089e-01, 1.636176e-02, 1.153102e-02, 1.325863e-01, 2.774239e-02, 1.754486e-01, 1.209832e-02, 2.025930e-03, 5.634031e-05, 3.546073e-01, 1.413926e-01, 2.180552e-02, 1.153102e-02, 2.180552e-02, 4.956812e-03, 3.262838e-01, 9.925057e-02 ), doubleArrayOf( // Bonferroni 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 7.019185e-01, 1.000000e+00, 1.000000e+00, 2.020365e-01, 1.516674e-02, 5.625735e-01, 1.000000e+00, 2.909271e-02, 1.537741e-02, 4.125636e-01, 6.782670e-02, 6.803480e-01, 1.882294e-02, 9.005725e-04, 1.252228e-05, 1.000000e+00, 4.713920e-01, 4.395577e-02, 1.088915e-02, 4.846527e-02, 3.305125e-03, 1.000000e+00, 2.867745e-01 ), doubleArrayOf( // Hochberg 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 4.632662e-01, 9.991834e-01, 9.991834e-01, 1.575885e-01, 1.383967e-02, 3.938014e-01, 7.600230e-01, 2.501973e-02, 1.383967e-02, 3.052971e-01, 5.426136e-02, 4.626366e-01, 1.656419e-02, 8.825610e-04, 1.252228e-05, 9.930759e-01, 3.394022e-01, 3.692284e-02, 1.023581e-02, 3.974152e-02, 3.172920e-03, 8.992520e-01, 2.179486e-01 ), doubleArrayOf( // Holm 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 4.632662e-01, 1.000000e+00, 1.000000e+00, 1.575885e-01, 1.395341e-02, 3.938014e-01, 7.600230e-01, 2.501973e-02, 1.395341e-02, 3.052971e-01, 5.426136e-02, 4.626366e-01, 1.656419e-02, 8.825610e-04, 1.252228e-05, 9.930759e-01, 3.394022e-01, 3.692284e-02, 1.023581e-02, 3.974152e-02, 3.172920e-03, 8.992520e-01, 2.179486e-01 ), doubleArrayOf( // Hommel 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.987624e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.595180e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 4.351895e-01, 9.991834e-01, 9.766522e-01, 1.414256e-01, 1.304340e-02, 3.530937e-01, 6.887709e-01, 2.385602e-02, 1.322457e-02, 2.722920e-01, 5.426136e-02, 4.218158e-01, 1.581127e-02, 8.825610e-04, 1.252228e-05, 8.743649e-01, 3.016908e-01, 3.516461e-02, 9.582456e-03, 3.877222e-02, 3.172920e-03, 8.122276e-01, 1.950067e-01 ) ) val types = listOf("bh", "by", "bonferroni", "hochberg", "holm", "hommel") val f = "\ntype %d = '%s' has cumulative error of %g" for (type in 0 until types.size) { val q = pAdjust(pvalues, types[type]) var error = 0.0 for (i in 0 until pvalues.size) { error += Math.abs(q[i] - correctAnswers[type][i]) } doubleSay(q) println(f.format(type, types[type], error)) }
}</lang>
- Output:
[ 1] 6.126681e-01 0.8521710465 0.1987205200 0.1891595417 0.3217789286 [ 5] 0.9301450000 0.4870370000 0.9301450000 0.6049730556 0.6826752564 [10] 0.6482628947 0.7253722500 0.5280972727 0.8769925556 0.4705703448 [15] 0.9241867391 0.6049730556 0.7856107317 0.4887525806 0.1136717045 [20] 0.4991890625 0.8769925556 0.9991834000 0.3217789286 0.9301450000 [25] 0.2304957692 0.5832475000 0.0389954722 0.8521710465 0.1476842609 [30] 0.0168363750 0.0025629017 0.0351608438 0.0625018947 0.0036365888 [35] 0.0025629017 0.0294688286 0.0061660636 0.0389954722 0.0026889914 [40] 0.0004502862 0.0000125223 0.0788155476 0.0314261300 0.0048465270 [45] 0.0025629017 0.0048465270 0.0011017083 0.0725203250 0.0220595769 [50] type 0 = 'bh' has cumulative error of 8.03053e-07 [ 1] 1.000000e+00 1.0000000000 0.8940844244 0.8510676197 1.0000000000 [ 5] 1.0000000000 1.0000000000 1.0000000000 1.0000000000 1.0000000000 [10] 1.0000000000 1.0000000000 1.0000000000 1.0000000000 1.0000000000 [15] 1.0000000000 1.0000000000 1.0000000000 1.0000000000 0.5114323399 [20] 1.0000000000 1.0000000000 1.0000000000 1.0000000000 1.0000000000 [25] 1.0000000000 1.0000000000 0.1754486368 1.0000000000 0.6644618149 [30] 0.0757503083 0.0115310209 0.1581958559 0.2812088585 0.0163617595 [35] 0.0115310209 0.1325863108 0.0277423864 0.1754486368 0.0120983246 [40] 0.0020259303 0.0000563403 0.3546073326 0.1413926119 0.0218055202 [45] 0.0115310209 0.0218055202 0.0049568120 0.3262838334 0.0992505663 [50] type 1 = 'by' has cumulative error of 3.64072e-07 [ 1] 1.000000e+00 1.0000000000 1.0000000000 1.0000000000 1.0000000000 [ 5] 1.0000000000 1.0000000000 1.0000000000 1.0000000000 1.0000000000 [10] 1.0000000000 1.0000000000 1.0000000000 1.0000000000 1.0000000000 [15] 1.0000000000 1.0000000000 1.0000000000 1.0000000000 1.0000000000 [20] 1.0000000000 1.0000000000 1.0000000000 1.0000000000 1.0000000000 [25] 1.0000000000 1.0000000000 0.7019185000 1.0000000000 1.0000000000 [30] 0.2020365000 0.0151667450 0.5625735000 1.0000000000 0.0290927100 [35] 0.0153774100 0.4125636000 0.0678267000 0.6803480000 0.0188229400 [40] 0.0009005725 0.0000125223 1.0000000000 0.4713919500 0.0439557650 [45] 0.0108891550 0.0484652700 0.0033051250 1.0000000000 0.2867745000 [50] type 2 = 'bonferroni' has cumulative error of 6.50000e-08 [ 1] 9.991834e-01 0.9991834000 0.9991834000 0.9991834000 0.9991834000 [ 5] 0.9991834000 0.9991834000 0.9991834000 0.9991834000 0.9991834000 [10] 0.9991834000 0.9991834000 0.9991834000 0.9991834000 0.9991834000 [15] 0.9991834000 0.9991834000 0.9991834000 0.9991834000 0.9991834000 [20] 0.9991834000 0.9991834000 0.9991834000 0.9991834000 0.9991834000 [25] 0.9991834000 0.9991834000 0.4632662100 0.9991834000 0.9991834000 [30] 0.1575884700 0.0138396690 0.3938014500 0.7600230400 0.0250197306 [35] 0.0138396690 0.3052970640 0.0542613600 0.4626366400 0.0165641872 [40] 0.0008825610 0.0000125223 0.9930759000 0.3394022040 0.0369228426 [45] 0.0102358057 0.0397415214 0.0031729200 0.8992520300 0.2179486200 [50] type 3 = 'hochberg' has cumulative error of 2.73750e-07 [ 1] 1.000000e+00 1.0000000000 1.0000000000 1.0000000000 1.0000000000 [ 5] 1.0000000000 1.0000000000 1.0000000000 1.0000000000 1.0000000000 [10] 1.0000000000 1.0000000000 1.0000000000 1.0000000000 1.0000000000 [15] 1.0000000000 1.0000000000 1.0000000000 1.0000000000 1.0000000000 [20] 1.0000000000 1.0000000000 1.0000000000 1.0000000000 1.0000000000 [25] 1.0000000000 1.0000000000 0.4632662100 1.0000000000 1.0000000000 [30] 0.1575884700 0.0139534054 0.3938014500 0.7600230400 0.0250197306 [35] 0.0139534054 0.3052970640 0.0542613600 0.4626366400 0.0165641872 [40] 0.0008825610 0.0000125223 0.9930759000 0.3394022040 0.0369228426 [45] 0.0102358057 0.0397415214 0.0031729200 0.8992520300 0.2179486200 [50] type 4 = 'holm' has cumulative error of 2.80950e-07 [ 1] 9.991834e-01 0.9991834000 0.9991834000 0.9987623800 0.9991834000 [ 5] 0.9991834000 0.9991834000 0.9991834000 0.9991834000 0.9991834000 [10] 0.9991834000 0.9991834000 0.9991834000 0.9991834000 0.9991834000 [15] 0.9991834000 0.9991834000 0.9991834000 0.9991834000 0.9595180000 [20] 0.9991834000 0.9991834000 0.9991834000 0.9991834000 0.9991834000 [25] 0.9991834000 0.9991834000 0.4351894700 0.9991834000 0.9766522500 [30] 0.1414255500 0.0130434007 0.3530936533 0.6887708800 0.0238560222 [35] 0.0132245726 0.2722919760 0.0542613600 0.4218157600 0.0158112696 [40] 0.0008825610 0.0000125223 0.8743649143 0.3016908480 0.0351646120 [45] 0.0095824564 0.0387722160 0.0031729200 0.8122276400 0.1950066600 [50] type 5 = 'hommel' has cumulative error of 4.35302e-07
Perl
<lang perl>#!/usr/bin/env perl
use strict; use warnings; use Cwd; my $TOP_DIRECTORY = getcwd(); local $SIG{__WARN__} = sub {#kill the program if there are any warnings my $message = shift; my $fail_filename = "$TOP_DIRECTORY/$0.FAIL"; open my $fh, '>', $fail_filename or die "Can't write $fail_filename: $!"; printf $fh ("$message @ %s\n", getcwd()); close $fh; die "$message\n"; };#http://perlmaven.com/how-to-capture-and-save-warnings-in-perl
sub pmin { my $array_ref = shift; my $x = 1; unless ((ref $array_ref) =~ m/ARRAY/) { print "cummin requires an array.\n"; die; } my @pmin_array; my $n = scalar @$array_ref; for (my $index = 0; $index < $n; $index++) { if (@$array_ref[$index] < $x) { $pmin_array[$index] = @$array_ref[$index]; } else { $pmin_array[$index] = $x; } } return @pmin_array; }
sub cummin { my $array_ref = shift; unless ((ref $array_ref) =~ m/ARRAY/) { print "cummin requires an array.\n"; die; } my @cummin; my $cumulative_min = @$array_ref[0]; foreach my $p (@$array_ref) { if ($p < $cumulative_min) { $cumulative_min = $p; } push @cummin, $cumulative_min; } return @cummin; }
sub cummax { my $array_ref = shift; unless ((ref $array_ref) =~ m/ARRAY/) { print "cummin requires an array.\n"; die; } my @cummax; my $cumulative_max = @$array_ref[0]; foreach my $p (@$array_ref) { if ($p > $cumulative_max) { $cumulative_max = $p; } push @cummax, $cumulative_max; } return @cummax; }
sub order {#made to match R's "order" my $array_ref = shift; my $decreasing = 'false'; if (defined $_[0]) { my $option = shift; if ($option =~ m/true/i) { $decreasing = 'true'; } elsif ($option =~ m/false/i) { #do nothing, it's already set to false } else { print "2nd option should only be case-insensitive 'true' or 'false'"; die; } } unless ((ref $array_ref) =~ m/ARRAY/) { print "You should have entered an array.\n"; die; } my @array; my $max_index = scalar @$array_ref-1; if ($decreasing eq 'false') { @array = sort { @$array_ref[$a] <=> @$array_ref[$b] } 0..$max_index; } elsif ($decreasing eq 'true') { @array = sort { @$array_ref[$b] <=> @$array_ref[$a] } 0..$max_index; } }
sub min { my $array_ref = shift; unless ((ref $array_ref) =~ m/ARRAY/) { print "min requires an array.\n"; die; } my $min = @$array_ref[0]; foreach my $e (@$array_ref) { if ($e < $min) { $min = $e; } } return $min; }
sub min2 { my $n1 = shift; my $n2 = shift; if ($n1 < $n2) { return $n1; } else { return $n2; } }
sub p_adjust { my $pvalues_ref = shift; unless ((ref $pvalues_ref) =~ m/ARRAY/) { print "p_adjust requires an array.\n"; die; } my $method; if (defined $_[0]) { $method = shift; } else { $method = 'Holm'; } my %methods = ( 'bh' => 1, 'fdr' => 1, 'by' => 1, 'holm' => 1, 'hommel' => 1, 'bonferroni' => 1, 'hochberg' => 1 ); my $method_found = 'no'; foreach my $key (keys %methods) { if ((uc $method) eq (uc $key)) { $method = $key; $method_found = 'yes'; last; } } if ($method_found eq 'no') { if ($method =~ m/benjamini-?\s*hochberg/i) { $method = 'bh'; $method_found = 'yes'; } elsif ($method =~ m/benjamini-?\s*yekutieli/i) { $method = 'by'; $method_found = 'yes'; } } if ($method_found eq 'no') { print "No method could be determined from $method.\n"; die; } my $lp = scalar @$pvalues_ref; my $n = $lp; my @qvalues; if ($method eq 'hochberg') { my @o = order($pvalues_ref, 'TRUE'); my @cummin_input; for (my $index = 0; $index < $n; $index++) { $cummin_input[$index] = ($index+1)* @$pvalues_ref[$o[$index]];#PVALUES[$o[$index]] is p[o] } my @cummin = cummin(\@cummin_input); undef @cummin_input; my @pmin = pmin(\@cummin); undef @cummin; my @ro = order(\@o); undef @o; @qvalues = @pmin[@ro]; } elsif ($method eq 'bh') { my @o = order($pvalues_ref, 'TRUE'); my @cummin_input; for (my $index = 0; $index < $n; $index++) { $cummin_input[$index] = ($n/($n-$index))* @$pvalues_ref[$o[$index]];#PVALUES[$o[$index]] is p[o] } my @ro = order(\@o); undef @o; my @cummin = cummin(\@cummin_input); undef @cummin_input; my @pmin = pmin(\@cummin); undef @cummin; @qvalues = @pmin[@ro]; } elsif ($method eq 'by') { my $q = 0.0; my @o = order($pvalues_ref, 'TRUE'); my @ro = order(\@o); for (my $index = 1; $index < ($n+1); $index++) { $q += 1.0 / $index; } my @cummin_input; for (my $index = 0; $index < $n; $index++) { $cummin_input[$index] = $q * ($n/($n-$index)) * @$pvalues_ref[$o[$index]];#PVALUES[$o[$index]] is p[o] } my @cummin = cummin(\@cummin_input); undef @cummin_input; my @pmin = pmin(\@cummin); undef @cummin; @qvalues = @pmin[@ro]; } elsif ($method eq 'bonferroni') { for (my $index = 0; $index < $n; $index++) { my $q = @$pvalues_ref[$index]*$n; if ((0 <= $q) && ($q < 1)) { $qvalues[$index] = $q; } elsif ($q >= 1) { $qvalues[$index] = 1.0; } else { print "Failed to get Bonferroni adjusted p."; die; } } } elsif ($method eq 'holm') { my @o = order($pvalues_ref); my @cummax_input; for (my $index = 0; $index < $n; $index++) { $cummax_input[$index] = ($n - $index) * @$pvalues_ref[$o[$index]]; } my @ro = order(\@o); undef @o; my @cummax = cummax(\@cummax_input); undef @cummax_input; my @pmin = pmin(\@cummax); undef @cummax; @qvalues = @pmin[@ro]; } elsif ($method eq 'hommel') { my @i = 1..$n; my @o = order($pvalues_ref); my @p = @$pvalues_ref[@o]; my @ro = order(\@o); undef @o; my @pa; my @q; my $min = $n*$p[0]; for (my $index = 0; $index < $n; $index++) { my $temp = $n*$p[$index] / ($index + 1); if ($temp < $min) { $min = $temp; } } for (my $index = 0; $index < $n; $index++) { $pa[$index] = $min;#q <- pa <- rep.int(min(n * p/i), n) $q[$index] = $min;#q <- pa <- rep.int(min(n * p/i), n) } for (my $j = ($n-1); $j >= 2; $j--) {
- printf("j = %zu\n", j);
my @ij = 1..($n - $j + 1);#ij <- seq_len(n - j + 1) for (my $i = 0; $i < $n - $j + 1; $i++) { $ij[$i]--;#R's indices are 1-based, C's are 0-based } my $I2_LENGTH = $j - 1; my @i2; for (my $i = 0; $i < $I2_LENGTH; $i++) { $i2[$i] = $n-$j+2+$i-1;
- R's indices are 1-based, C's are 0-based, I added the -1
}
my $q1 = $j * $p[$i2[0]] / 2.0; for (my $i = 1; $i < $I2_LENGTH; $i++) {#loop through 2:j my $TEMP_Q1 = $j * $p[$i2[$i]] / (2 + $i); if ($TEMP_Q1 < $q1) { $q1 = $TEMP_Q1; } }
for (my $i = 0; $i < ($n - $j + 1); $i++) {#q[ij] <- pmin(j * p[ij], q1) $q[$ij[$i]] = min2( $j*$p[$ij[$i]], $q1); }
for (my $i = 0; $i < $I2_LENGTH; $i++) {#q[i2] <- q[n - j + 1] $q[$i2[$i]] = $q[$n - $j];#subtract 1 because of starting index difference }
for (my $i = 0; $i < $n; $i++) {#pa <- pmax(pa, q) if ($pa[$i] < $q[$i]) { $pa[$i] = $q[$i]; } }
- printf("j = %zu, pa = \n", j);
- double_say(pa, N);
}#end j loop @qvalues = @pa[@ro]; } else { print "$method doesn't fit my types.\n"; die; } return @qvalues; } my @pvalues = (4.533744e-01, 7.296024e-01, 9.936026e-02, 9.079658e-02, 1.801962e-01, 8.752257e-01, 2.922222e-01, 9.115421e-01, 4.355806e-01, 5.324867e-01, 4.926798e-01, 5.802978e-01, 3.485442e-01, 7.883130e-01, 2.729308e-01, 8.502518e-01, 4.268138e-01, 6.442008e-01, 3.030266e-01, 5.001555e-02, 3.194810e-01, 7.892933e-01, 9.991834e-01, 1.745691e-01, 9.037516e-01, 1.198578e-01, 3.966083e-01, 1.403837e-02, 7.328671e-01, 6.793476e-02, 4.040730e-03, 3.033349e-04, 1.125147e-02, 2.375072e-02, 5.818542e-04, 3.075482e-04, 8.251272e-03, 1.356534e-03, 1.360696e-02, 3.764588e-04, 1.801145e-05, 2.504456e-07, 3.310253e-02, 9.427839e-03, 8.791153e-04, 2.177831e-04, 9.693054e-04, 6.610250e-05, 2.900813e-02, 5.735490e-03);
my %correct_answers = ( 'Benjamini-Hochberg' => [6.126681e-01, 8.521710e-01, 1.987205e-01, 1.891595e-01, 3.217789e-01, 9.301450e-01, 4.870370e-01, 9.301450e-01, 6.049731e-01, 6.826753e-01, 6.482629e-01, 7.253722e-01, 5.280973e-01, 8.769926e-01, 4.705703e-01, 9.241867e-01, 6.049731e-01, 7.856107e-01, 4.887526e-01, 1.136717e-01, 4.991891e-01, 8.769926e-01, 9.991834e-01, 3.217789e-01, 9.301450e-01, 2.304958e-01, 5.832475e-01, 3.899547e-02, 8.521710e-01, 1.476843e-01, 1.683638e-02, 2.562902e-03, 3.516084e-02, 6.250189e-02, 3.636589e-03, 2.562902e-03, 2.946883e-02, 6.166064e-03, 3.899547e-02, 2.688991e-03, 4.502862e-04, 1.252228e-05, 7.881555e-02, 3.142613e-02, 4.846527e-03, 2.562902e-03, 4.846527e-03, 1.101708e-03, 7.252032e-02, 2.205958e-02], 'Benjamini-Yekutieli' => [1.000000e+00, 1.000000e+00, 8.940844e-01, 8.510676e-01, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 5.114323e-01, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.754486e-01, 1.000000e+00, 6.644618e-01, 7.575031e-02, 1.153102e-02, 1.581959e-01, 2.812089e-01, 1.636176e-02, 1.153102e-02, 1.325863e-01, 2.774239e-02, 1.754486e-01, 1.209832e-02, 2.025930e-03, 5.634031e-05, 3.546073e-01, 1.413926e-01, 2.180552e-02, 1.153102e-02, 2.180552e-02, 4.956812e-03, 3.262838e-01, 9.925057e-02], 'Bonferroni' => [1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 7.019185e-01, 1.000000e+00, 1.000000e+00, 2.020365e-01, 1.516674e-02, 5.625735e-01, 1.000000e+00, 2.909271e-02, 1.537741e-02, 4.125636e-01, 6.782670e-02, 6.803480e-01, 1.882294e-02, 9.005725e-04, 1.252228e-05, 1.000000e+00, 4.713920e-01, 4.395577e-02, 1.088915e-02, 4.846527e-02, 3.305125e-03, 1.000000e+00, 2.867745e-01],
'Hochberg' => [9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 4.632662e-01, 9.991834e-01, 9.991834e-01, 1.575885e-01, 1.383967e-02, 3.938014e-01, 7.600230e-01, 2.501973e-02, 1.383967e-02, 3.052971e-01, 5.426136e-02, 4.626366e-01, 1.656419e-02, 8.825610e-04, 1.252228e-05, 9.930759e-01, 3.394022e-01, 3.692284e-02, 1.023581e-02, 3.974152e-02, 3.172920e-03, 8.992520e-01, 2.179486e-01], 'Holm' => [1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 1.000000e+00, 4.632662e-01, 1.000000e+00, 1.000000e+00, 1.575885e-01, 1.395341e-02, 3.938014e-01, 7.600230e-01, 2.501973e-02, 1.395341e-02, 3.052971e-01, 5.426136e-02, 4.626366e-01, 1.656419e-02, 8.825610e-04, 1.252228e-05, 9.930759e-01, 3.394022e-01, 3.692284e-02, 1.023581e-02, 3.974152e-02, 3.172920e-03, 8.992520e-01, 2.179486e-01],
'Hommel' => [9.991834e-01, 9.991834e-01, 9.991834e-01, 9.987624e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.595180e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 9.991834e-01, 4.351895e-01, 9.991834e-01, 9.766522e-01, 1.414256e-01, 1.304340e-02, 3.530937e-01, 6.887709e-01, 2.385602e-02, 1.322457e-02, 2.722920e-01, 5.426136e-02, 4.218158e-01, 1.581127e-02, 8.825610e-04, 1.252228e-05, 8.743649e-01, 3.016908e-01, 3.516461e-02, 9.582456e-03, 3.877222e-02, 3.172920e-03, 8.122276e-01, 1.950067e-01]);
foreach my $method ('Hochberg','Benjamini-Hochberg','Benjamini-Yekutieli', 'Bonferroni', 'Holm', 'Hommel') {
print "$method\n";
my @qvalues = p_adjust(\@pvalues, $method);
my $error = 0.0;
foreach my $q (0..$#qvalues) {
$error += abs($qvalues[$q] - $correct_answers{$method}[$q]);
}
printf("type $method has cumulative error of %g.\n", $error);
}
</lang>
- Output:
Hochberg type Hochberg has cumulative error of 2.7375e-07. Benjamini-Hochberg type Benjamini-Hochberg has cumulative error of 8.03053e-07. Benjamini-Yekutieli type Benjamini-Yekutieli has cumulative error of 3.64072e-07. Bonferroni type Bonferroni has cumulative error of 6.5e-08. Holm type Holm has cumulative error of 2.8095e-07. Hommel type Hommel has cumulative error of 4.35302e-07.
Perl 6
<lang perl6>########################### Helper subs ###########################
sub adjusted (@p, $type) { "\n$type\n" ~ format adjust( check(@p), $type ) }
sub format ( @p, $cols = 5 ) {
my $i = -$cols; my $fmt = "%1.10f"; join "\n", @p.rotor($cols, :partial).map: { sprintf "[%2d] { join ' ', $fmt xx $_ }", $i+=$cols, $_ };
}
sub check ( @p ) { die 'p-values must be in range 0.0 to 1.0' if @p.min < 0 or 1 < @p.max; @p }
multi ratchet ( 'up', @p ) { my $m; @p[$_] min= $m, $m = @p[$_] for ^@p; @p }
multi ratchet ( 'dn', @p ) { my $m; @p[$_] max= $m, $m = @p[$_] for ^@p .reverse; @p }
sub schwartzian ( @p, &transform, :$ratchet ) {
my @pa = @p.map( {[$_, $++]} ).sort( -*.[0] ).map( {[transform(.[0]), .[1]]} ); @pa[*;0] = ratchet($ratchet, @pa»[0]); @pa.sort( *.[1] )»[0]
}
- The various p-value correction routines #############
multi adjust( @p, 'Benjamini-Hochberg' ) {
@p.&schwartzian: * * @p / (@p - $++) min 1, :ratchet('up')
}
multi adjust( @p, 'Benjamini-Yekutieli' ) {
my \r = ^@p .map( { 1 / ++$ } ).sum; @p.&schwartzian: * * r * @p / (@p - $++) min 1, :ratchet('up')
}
multi adjust( @p, 'Hochberg' ) {
my \m = @p.max; @p.&schwartzian: * * ++$ min m, :ratchet('up')
}
multi adjust( @p, 'Holm' ) {
@p.&schwartzian: * * ++$ min 1, :ratchet('dn')
}
multi adjust( @p, 'Šidák' ) {
@p.&schwartzian: 1 - (1 - *) ** ++$, :ratchet('dn')
}
multi adjust( @p, 'Bonferroni' ) {
@p.map: * * @p min 1
}
- Hommel correction can't be easily reduced to a one pass transform
multi adjust( @p, 'Hommel' ) {
my @s = @p.map( {[$_, $++]} ).sort( *.[0] ).values; # sorted my \z = +@p; # array si(z)e my @pa = @s»[0].map( * * z / ++$ ).min xx z; # p adjusted my @q; # scratch array for (1 ..^ z).reverse -> $i { my @L = 0 .. z - $i; # lower indices my @U = z - $i ^..^ z; # upper indices my $q = @s[@U]»[0].map( { $_ * $i / (2 + $++) } ).min; @q[@L] = @s[@L]»[0].map: { [min] $_ * $i, $q, @s[*-1][0] }; @pa = ^z .map: { [max] @pa[$_], @q[$_] } } @pa[@s[*;1].map( {[$_, $++]} ).sort( *.[0] ).values»[1]]
}
- The task ###########################
my @p-values =
4.533744e-01, 7.296024e-01, 9.936026e-02, 9.079658e-02, 1.801962e-01, 8.752257e-01, 2.922222e-01, 9.115421e-01, 4.355806e-01, 5.324867e-01, 4.926798e-01, 5.802978e-01, 3.485442e-01, 7.883130e-01, 2.729308e-01, 8.502518e-01, 4.268138e-01, 6.442008e-01, 3.030266e-01, 5.001555e-02, 3.194810e-01, 7.892933e-01, 9.991834e-01, 1.745691e-01, 9.037516e-01, 1.198578e-01, 3.966083e-01, 1.403837e-02, 7.328671e-01, 6.793476e-02, 4.040730e-03, 3.033349e-04, 1.125147e-02, 2.375072e-02, 5.818542e-04, 3.075482e-04, 8.251272e-03, 1.356534e-03, 1.360696e-02, 3.764588e-04, 1.801145e-05, 2.504456e-07, 3.310253e-02, 9.427839e-03, 8.791153e-04, 2.177831e-04, 9.693054e-04, 6.610250e-05, 2.900813e-02, 5.735490e-03
for < Benjamini-Hochberg Benjamini-Yekutieli Bonferroni Hochberg Holm Hommel Šidák > {
say adjusted @p-values, $_
}</lang>
- Output:
Benjamini-Hochberg [ 0] 0.6126681081 0.8521710465 0.1987205200 0.1891595417 0.3217789286 [ 5] 0.9301450000 0.4870370000 0.9301450000 0.6049730556 0.6826752564 [10] 0.6482628947 0.7253722500 0.5280972727 0.8769925556 0.4705703448 [15] 0.9241867391 0.6049730556 0.7856107317 0.4887525806 0.1136717045 [20] 0.4991890625 0.8769925556 0.9991834000 0.3217789286 0.9301450000 [25] 0.2304957692 0.5832475000 0.0389954722 0.8521710465 0.1476842609 [30] 0.0168363750 0.0025629017 0.0351608438 0.0625018947 0.0036365888 [35] 0.0025629017 0.0294688286 0.0061660636 0.0389954722 0.0026889914 [40] 0.0004502863 0.0000125223 0.0788155476 0.0314261300 0.0048465270 [45] 0.0025629017 0.0048465270 0.0011017083 0.0725203250 0.0220595769 Benjamini-Yekutieli [ 0] 1.0000000000 1.0000000000 0.8940844244 0.8510676197 1.0000000000 [ 5] 1.0000000000 1.0000000000 1.0000000000 1.0000000000 1.0000000000 [10] 1.0000000000 1.0000000000 1.0000000000 1.0000000000 1.0000000000 [15] 1.0000000000 1.0000000000 1.0000000000 1.0000000000 0.5114323399 [20] 1.0000000000 1.0000000000 1.0000000000 1.0000000000 1.0000000000 [25] 1.0000000000 1.0000000000 0.1754486368 1.0000000000 0.6644618149 [30] 0.0757503083 0.0115310209 0.1581958559 0.2812088585 0.0163617595 [35] 0.0115310209 0.1325863108 0.0277423864 0.1754486368 0.0120983246 [40] 0.0020259303 0.0000563403 0.3546073326 0.1413926119 0.0218055202 [45] 0.0115310209 0.0218055202 0.0049568120 0.3262838334 0.0992505663 Bonferroni [ 0] 1.0000000000 1.0000000000 1.0000000000 1.0000000000 1.0000000000 [ 5] 1.0000000000 1.0000000000 1.0000000000 1.0000000000 1.0000000000 [10] 1.0000000000 1.0000000000 1.0000000000 1.0000000000 1.0000000000 [15] 1.0000000000 1.0000000000 1.0000000000 1.0000000000 1.0000000000 [20] 1.0000000000 1.0000000000 1.0000000000 1.0000000000 1.0000000000 [25] 1.0000000000 1.0000000000 0.7019185000 1.0000000000 1.0000000000 [30] 0.2020365000 0.0151667450 0.5625735000 1.0000000000 0.0290927100 [35] 0.0153774100 0.4125636000 0.0678267000 0.6803480000 0.0188229400 [40] 0.0009005725 0.0000125223 1.0000000000 0.4713919500 0.0439557650 [45] 0.0108891550 0.0484652700 0.0033051250 1.0000000000 0.2867745000 Hochberg [ 0] 0.9991834000 0.9991834000 0.9991834000 0.9991834000 0.9991834000 [ 5] 0.9991834000 0.9991834000 0.9991834000 0.9991834000 0.9991834000 [10] 0.9991834000 0.9991834000 0.9991834000 0.9991834000 0.9991834000 [15] 0.9991834000 0.9991834000 0.9991834000 0.9991834000 0.9991834000 [20] 0.9991834000 0.9991834000 0.9991834000 0.9991834000 0.9991834000 [25] 0.9991834000 0.9991834000 0.4632662100 0.9991834000 0.9991834000 [30] 0.1575884700 0.0138396690 0.3938014500 0.7600230400 0.0250197306 [35] 0.0138396690 0.3052970640 0.0542613600 0.4626366400 0.0165641872 [40] 0.0008825611 0.0000125223 0.9930759000 0.3394022040 0.0369228426 [45] 0.0102358057 0.0397415214 0.0031729200 0.8992520300 0.2179486200 Holm [ 0] 1.0000000000 1.0000000000 1.0000000000 1.0000000000 1.0000000000 [ 5] 1.0000000000 1.0000000000 1.0000000000 1.0000000000 1.0000000000 [10] 1.0000000000 1.0000000000 1.0000000000 1.0000000000 1.0000000000 [15] 1.0000000000 1.0000000000 1.0000000000 1.0000000000 1.0000000000 [20] 1.0000000000 1.0000000000 1.0000000000 1.0000000000 1.0000000000 [25] 1.0000000000 1.0000000000 0.4632662100 1.0000000000 1.0000000000 [30] 0.1575884700 0.0139534054 0.3938014500 0.7600230400 0.0250197306 [35] 0.0139534054 0.3052970640 0.0542613600 0.4626366400 0.0165641872 [40] 0.0008825611 0.0000125223 0.9930759000 0.3394022040 0.0369228426 [45] 0.0102358057 0.0397415214 0.0031729200 0.8992520300 0.2179486200 Hommel [ 0] 0.9991834000 0.9991834000 0.9991834000 0.9987623800 0.9991834000 [ 5] 0.9991834000 0.9991834000 0.9991834000 0.9991834000 0.9991834000 [10] 0.9991834000 0.9991834000 0.9991834000 0.9991834000 0.9991834000 [15] 0.9991834000 0.9991834000 0.9991834000 0.9991834000 0.9595180000 [20] 0.9991834000 0.9991834000 0.9991834000 0.9991834000 0.9991834000 [25] 0.9991834000 0.9991834000 0.4351894700 0.9991834000 0.9766522500 [30] 0.1414255500 0.0130434007 0.3530936533 0.6887708800 0.0238560222 [35] 0.0132245726 0.2722919760 0.0542613600 0.4218157600 0.0158112696 [40] 0.0008825611 0.0000125223 0.8743649143 0.3016908480 0.0351646120 [45] 0.0095824564 0.0387722160 0.0031729200 0.8122276400 0.1950066600 Šidák [ 0] 0.9998642526 0.9999922727 0.9341844137 0.9234670175 0.9899922294 [ 5] 0.9999922727 0.9992955735 0.9999922727 0.9998642526 0.9998909746 [10] 0.9998642526 0.9999288207 0.9995533892 0.9999922727 0.9990991210 [15] 0.9999922727 0.9998642526 0.9999674876 0.9992955735 0.7741716825 [20] 0.9993332472 0.9999922727 0.9999922727 0.9899922294 0.9999922727 [25] 0.9589019598 0.9998137104 0.3728369461 0.9999922727 0.8605248833 [30] 0.1460714182 0.0138585952 0.3270159382 0.5366136349 0.0247164330 [35] 0.0138585952 0.2640282766 0.0528503728 0.3723753774 0.0164308228 [40] 0.0008821796 0.0000125222 0.6357389664 0.2889497995 0.0362651575 [45] 0.0101847015 0.0389807074 0.0031679962 0.5985019850 0.1963376344
R
<lang R>
- p.adjust's source code is below
- > p.adjust
- function (p, method = p.adjust.methods, n = length(p))
- {
- method <- match.arg(method)
- if (method == "fdr")
- method <- "BH"
- nm <- names(p)
- p <- as.numeric(p)
- p0 <- setNames(p, nm)
- if (all(nna <- !is.na(p)))
- nna <- TRUE
- p <- p[nna]
- lp <- length(p)
- stopifnot(n >= lp)
# if (n <= 1) # return(p0)
- if (n == 2 && method == "hommel")
- method <- "hochberg"
- p0[nna] <- switch(method, bonferroni = pmin(1, n * p), holm = {
- i <- seq_len(lp)
# o <- order(p)
- ro <- order(o)
- pmin(1, cummax((n - i + 1L) * p[o]))[ro]
- }, hommel = {
# if (n > lp) p <- c(p, rep.int(1, n - lp))
- i <- seq_len(n)
# o <- order(p)
- p <- p[o]
- ro <- order(o)
- q <- pa <- rep.int(min(n * p/i), n)
- for (j in (n - 1):2) {
# ij <- seq_len(n - j + 1)
- i2 <- (n - j + 2):n
- q1 <- min(j * p[i2]/(2:j))
- q[ij] <- pmin(j * p[ij], q1)
- q[i2] <- q[n - j + 1]
- pa <- pmax(pa, q)
- }
- pmax(pa, p)[if (lp < n) ro[1:lp] else ro]
- }, hochberg = {
# i <- lp:1L
- o <- order(p, decreasing = TRUE)
- ro <- order(o)
- pmin(1, cummin((n - i + 1L) * p[o]))[ro]
- }, BH = {
- i <- lp:1L
- o <- order(p, decreasing = TRUE)
- ro <- order(o)
- pmin(1, cummin(n/i * p[o]))[ro]
- }, BY = {
- i <- lp:1L
- o <- order(p, decreasing = TRUE)
- ro <- order(o)
- q <- sum(1L/(1L:n))
# pmin(1, cummin(q * n/i * p[o]))[ro]
- }, none = p)
- p0
- }
- <bytecode: 0x3a61b88>
- <environment: namespace:stats>
p <- c(4.533744e-01, 7.296024e-01, 9.936026e-02, 9.079658e-02, 1.801962e-01, 8.752257e-01, 2.922222e-01, 9.115421e-01, 4.355806e-01, 5.324867e-01, 4.926798e-01, 5.802978e-01, 3.485442e-01, 7.883130e-01, 2.729308e-01, 8.502518e-01, 4.268138e-01, 6.442008e-01, 3.030266e-01, 5.001555e-02, 3.194810e-01, 7.892933e-01, 9.991834e-01, 1.745691e-01, 9.037516e-01, 1.198578e-01, 3.966083e-01, 1.403837e-02, 7.328671e-01, 6.793476e-02, 4.040730e-03, 3.033349e-04, 1.125147e-02, 2.375072e-02, 5.818542e-04, 3.075482e-04, 8.251272e-03, 1.356534e-03, 1.360696e-02, 3.764588e-04, 1.801145e-05, 2.504456e-07, 3.310253e-02, 9.427839e-03, 8.791153e-04, 2.177831e-04, 9.693054e-04, 6.610250e-05, 2.900813e-02, 5.735490e-03)
p.adjust(p, method = 'BH') print("Benjamini-Hochberg") writeLines("\n")
p.adjust(p, method = 'BY') print("Benjamini & Yekutieli") writeLines("\n")
p.adjust(p, method = 'bonferroni') print("Bonferroni") writeLines("\n")
p.adjust(p, method = 'hochberg') print("Hochberg") writeLines("\n");
p.adjust(p, method = 'hommel') writeLines("Hommel\n") </lang>
- Output:
[1] 6.126681e-01 8.521710e-01 1.987205e-01 1.891595e-01 3.217789e-01 [6] 9.301450e-01 4.870370e-01 9.301450e-01 6.049731e-01 6.826753e-01 [11] 6.482629e-01 7.253722e-01 5.280973e-01 8.769926e-01 4.705703e-01 [16] 9.241867e-01 6.049731e-01 7.856107e-01 4.887526e-01 1.136717e-01 [21] 4.991891e-01 8.769926e-01 9.991834e-01 3.217789e-01 9.301450e-01 [26] 2.304958e-01 5.832475e-01 3.899547e-02 8.521710e-01 1.476843e-01 [31] 1.683638e-02 2.562902e-03 3.516084e-02 6.250189e-02 3.636589e-03 [36] 2.562902e-03 2.946883e-02 6.166064e-03 3.899547e-02 2.688991e-03 [41] 4.502862e-04 1.252228e-05 7.881555e-02 3.142613e-02 4.846527e-03 [46] 2.562902e-03 4.846527e-03 1.101708e-03 7.252032e-02 2.205958e-02 [1] "Benjamini-Hochberg" [1] 1.000000e+00 1.000000e+00 8.940844e-01 8.510676e-01 1.000000e+00 [6] 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 [11] 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 [16] 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 5.114323e-01 [21] 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 [26] 1.000000e+00 1.000000e+00 1.754486e-01 1.000000e+00 6.644618e-01 [31] 7.575031e-02 1.153102e-02 1.581959e-01 2.812089e-01 1.636176e-02 [36] 1.153102e-02 1.325863e-01 2.774239e-02 1.754486e-01 1.209832e-02 [41] 2.025930e-03 5.634031e-05 3.546073e-01 1.413926e-01 2.180552e-02 [46] 1.153102e-02 2.180552e-02 4.956812e-03 3.262838e-01 9.925057e-02 [1] "Benjamini & Yekutieli" [1] 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 [6] 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 [11] 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 [16] 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 [21] 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 [26] 1.000000e+00 1.000000e+00 7.019185e-01 1.000000e+00 1.000000e+00 [31] 2.020365e-01 1.516674e-02 5.625735e-01 1.000000e+00 2.909271e-02 [36] 1.537741e-02 4.125636e-01 6.782670e-02 6.803480e-01 1.882294e-02 [41] 9.005725e-04 1.252228e-05 1.000000e+00 4.713920e-01 4.395577e-02 [46] 1.088915e-02 4.846527e-02 3.305125e-03 1.000000e+00 2.867745e-01 [1] "Bonferroni" [1] 9.991834e-01 9.991834e-01 9.991834e-01 9.991834e-01 9.991834e-01 [6] 9.991834e-01 9.991834e-01 9.991834e-01 9.991834e-01 9.991834e-01 [11] 9.991834e-01 9.991834e-01 9.991834e-01 9.991834e-01 9.991834e-01 [16] 9.991834e-01 9.991834e-01 9.991834e-01 9.991834e-01 9.991834e-01 [21] 9.991834e-01 9.991834e-01 9.991834e-01 9.991834e-01 9.991834e-01 [26] 9.991834e-01 9.991834e-01 4.632662e-01 9.991834e-01 9.991834e-01 [31] 1.575885e-01 1.383967e-02 3.938014e-01 7.600230e-01 2.501973e-02 [36] 1.383967e-02 3.052971e-01 5.426136e-02 4.626366e-01 1.656419e-02 [41] 8.825610e-04 1.252228e-05 9.930759e-01 3.394022e-01 3.692284e-02 [46] 1.023581e-02 3.974152e-02 3.172920e-03 8.992520e-01 2.179486e-01 [1] "Hochberg" [1] 9.991834e-01 9.991834e-01 9.991834e-01 9.987624e-01 9.991834e-01 [6] 9.991834e-01 9.991834e-01 9.991834e-01 9.991834e-01 9.991834e-01 [11] 9.991834e-01 9.991834e-01 9.991834e-01 9.991834e-01 9.991834e-01 [16] 9.991834e-01 9.991834e-01 9.991834e-01 9.991834e-01 9.595180e-01 [21] 9.991834e-01 9.991834e-01 9.991834e-01 9.991834e-01 9.991834e-01 [26] 9.991834e-01 9.991834e-01 4.351895e-01 9.991834e-01 9.766522e-01 [31] 1.414256e-01 1.304340e-02 3.530937e-01 6.887709e-01 2.385602e-02 [36] 1.322457e-02 2.722920e-01 5.426136e-02 4.218158e-01 1.581127e-02 [41] 8.825610e-04 1.252228e-05 8.743649e-01 3.016908e-01 3.516461e-02 [46] 9.582456e-03 3.877222e-02 3.172920e-03 8.122276e-01 1.950067e-01 Hommel
zkl
<lang zkl>fcn bh(pvalues){ // Benjamini-Hochberg
psz,pszf := pvalues.len(), psz.toFloat(); n_i := psz.pump(List,'wrap(n){ pszf/(psz - n) }); # N/(N-0),N/(N-1),.. o,ro := order(pvalues,True),order(o,False); # sort pvalues, sort indices in := psz.pump(List,'wrap(n){ n_i[n]*pvalues[o[n]] }); pmin := cummin(in).apply((1.0).min); # (min(1,c[0]),min(1,c[1]),...) ro.apply(pmin.get); # (pmin[ro[0]],pmin[ro[1]],...)
}
fcn by(pvalues){ // Benjamini & Yekutieli
psz,pszf := pvalues.len(), psz.toFloat(); o,ro := order(pvalues,True),order(o,False); # sort pvalues, sort indices n_i := psz.pump(List,'wrap(n){ pszf/(psz - n) }); # N/(N-0),N/(N-1),.. q := [1..psz].reduce(fcn(q,n){ q+=1.0/n },0.0); in := psz.pump(List,'wrap(n){ q * n_i[n] * pvalues[o[n]] }); cummin(in).apply((1.0).min) : ro.apply(_.get);
}
fcn hochberg(pvalues){
psz,pszf := pvalues.len(), psz.toFloat(); o,ro := order(pvalues,True),order(o,False); # sort pvalues, sort indices n_i := psz.pump(List,'wrap(n){ pszf/(psz - n) }); # N/(N-0),N/(N-1),.. in := psz.pump(List,'wrap(n){ pvalues[o[n]]*(n + 1) }); cummin(in).apply((1.0).min) : ro.apply(_.get);
}
fcn cummin(pvalues){ // R's cumulative minima --> list of mins
out,m := List.createLong(pvalues.len()), pvalues[0]; foreach pv in (pvalues){ out.append(m=m.min(pv)) } out
} fcn order(list,downUp){ // True==increasing, --> List(int) sorted indices
f:=(downUp) and fcn(a,b){ a[1]>b[1] } or fcn(a,b){ a[1]<b[1] }; [0..].zip(list).pump(List()).sort(f).pump(List,T("get",0))
}
fcn bonferroni(pvalues){ // -->List
sz,r := pvalues.len(),List(); foreach pv in (pvalues){ b:=pv*sz; if(b>=1.0) r.append(1.0); else if(0.0<=b<1.0) r.append(b); else throw(Exception.ValueError(
"%g is outside of the interval I planned.".fmt(b)));
} r
}
fcn hommel(pvalues){
psz,indices := pvalues.len(), [1..psz].walk(); // 1,2,3,4... o,ro := order(pvalues,False),order(o,False); # sort pvalues, sort indices p := o.apply('wrap(n){ pvalues[n] }).copy(); // pvalues[*o] npi := [1..].zip(p).apply('wrap([(n,p)]){ p*psz/n }); min := (0.0).min(npi); // min value in npi pa,q := List.createLong(psz,min), pa.copy(); #(min,min,,,) foreach j in ([psz - 1..2,-1]){ ij:=[0..psz - j].walk(); i2:=(j - 1).pump(List,'+(psz - j + 1)); q1:=(0.0).min((j-1).pump(List,'wrap(n){ p[i2[n]]*j/(2 + n) })); foreach i in (psz - j + 1){ q[ij[i]] = q1.min(p[ij[i]]*j) } foreach i in (j - 1){ q[i2[i]] = q[psz - j] } foreach i in (psz){ pa[i] = pa[i].max(q[i]) } } psz.pump(List,'wrap(n){ pa[ro[n]] }); // Hommel q-values
}</lang> <lang zkl>pvalues:=T(
4.533744e-01, 7.296024e-01, 9.936026e-02, 9.079658e-02, 1.801962e-01, 8.752257e-01, 2.922222e-01, 9.115421e-01, 4.355806e-01, 5.324867e-01, 4.926798e-01, 5.802978e-01, 3.485442e-01, 7.883130e-01, 2.729308e-01, 8.502518e-01, 4.268138e-01, 6.442008e-01, 3.030266e-01, 5.001555e-02, 3.194810e-01, 7.892933e-01, 9.991834e-01, 1.745691e-01, 9.037516e-01, 1.198578e-01, 3.966083e-01, 1.403837e-02, 7.328671e-01, 6.793476e-02, 4.040730e-03, 3.033349e-04, 1.125147e-02, 2.375072e-02, 5.818542e-04, 3.075482e-04, 8.251272e-03, 1.356534e-03, 1.360696e-02, 3.764588e-04, 1.801145e-05, 2.504456e-07, 3.310253e-02, 9.427839e-03, 8.791153e-04, 2.177831e-04, 9.693054e-04, 6.610250e-05, 2.900813e-02, 5.735490e-03);
bh(pvalues) : format(_,"\nBenjamini-Hochberg"); by(pvalues) : format(_,"\nBenjamini & Yekutieli"); bonferroni(pvalues) : format(_,"\nBonferroni"); hochberg(pvalues) : format(_,"\nHochberg"); hommel(pvalues) : format(_,"\nHommel");
fcn format(list,title){
print(title,":"); foreach n in ([1..list.len(),5]){ print("\n[%2d]:".fmt(n)); foreach x in (list[n-1,5]){ print(" %.6e".fmt(x)) } } println();
}</lang>
- Output:
Benjamini-Hochberg: [ 1]: 6.126681e-01 8.521710e-01 1.987205e-01 1.891595e-01 3.217789e-01 [ 6]: 9.301450e-01 4.870370e-01 9.301450e-01 6.049731e-01 6.826753e-01 [11]: 6.482629e-01 7.253722e-01 5.280973e-01 8.769926e-01 4.705703e-01 [16]: 9.241867e-01 6.049731e-01 7.856107e-01 4.887526e-01 1.136717e-01 [21]: 4.991891e-01 8.769926e-01 9.991834e-01 3.217789e-01 9.301450e-01 [26]: 2.304958e-01 5.832475e-01 3.899547e-02 8.521710e-01 1.476843e-01 [31]: 1.683638e-02 2.562902e-03 3.516084e-02 6.250189e-02 3.636589e-03 [36]: 2.562902e-03 2.946883e-02 6.166064e-03 3.899547e-02 2.688991e-03 [41]: 4.502862e-04 1.252228e-05 7.881555e-02 3.142613e-02 4.846527e-03 [46]: 2.562902e-03 4.846527e-03 1.101708e-03 7.252032e-02 2.205958e-02 Benjamini & Yekutieli: [ 1]: 1.000000e+00 1.000000e+00 8.940844e-01 8.510676e-01 1.000000e+00 [ 6]: 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 [11]: 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 [16]: 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 5.114323e-01 [21]: 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 [26]: 1.000000e+00 1.000000e+00 1.754486e-01 1.000000e+00 6.644618e-01 [31]: 7.575031e-02 1.153102e-02 1.581959e-01 2.812089e-01 1.636176e-02 [36]: 1.153102e-02 1.325863e-01 2.774239e-02 1.754486e-01 1.209832e-02 [41]: 2.025930e-03 5.634031e-05 3.546073e-01 1.413926e-01 2.180552e-02 [46]: 1.153102e-02 2.180552e-02 4.956812e-03 3.262838e-01 9.925057e-02 Bonferroni: [ 1]: 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 [ 6]: 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 [11]: 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 [16]: 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 [21]: 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 [26]: 1.000000e+00 1.000000e+00 7.019185e-01 1.000000e+00 1.000000e+00 [31]: 2.020365e-01 1.516674e-02 5.625735e-01 1.000000e+00 2.909271e-02 [36]: 1.537741e-02 4.125636e-01 6.782670e-02 6.803480e-01 1.882294e-02 [41]: 9.005725e-04 1.252228e-05 1.000000e+00 4.713920e-01 4.395577e-02 [46]: 1.088915e-02 4.846527e-02 3.305125e-03 1.000000e+00 2.867745e-01 Hochberg: [ 1]: 9.991834e-01 9.991834e-01 9.991834e-01 9.991834e-01 9.991834e-01 [ 6]: 9.991834e-01 9.991834e-01 9.991834e-01 9.991834e-01 9.991834e-01 [11]: 9.991834e-01 9.991834e-01 9.991834e-01 9.991834e-01 9.991834e-01 [16]: 9.991834e-01 9.991834e-01 9.991834e-01 9.991834e-01 9.991834e-01 [21]: 9.991834e-01 9.991834e-01 9.991834e-01 9.991834e-01 9.991834e-01 [26]: 9.991834e-01 9.991834e-01 4.632662e-01 9.991834e-01 9.991834e-01 [31]: 1.575885e-01 1.383967e-02 3.938014e-01 7.600230e-01 2.501973e-02 [36]: 1.383967e-02 3.052971e-01 5.426136e-02 4.626366e-01 1.656419e-02 [41]: 8.825610e-04 1.252228e-05 9.930759e-01 3.394022e-01 3.692284e-02 [46]: 1.023581e-02 3.974152e-02 3.172920e-03 8.992520e-01 2.179486e-01 Hommel: [ 1]: 9.991834e-01 9.991834e-01 9.991834e-01 9.987624e-01 9.991834e-01 [ 6]: 9.991834e-01 9.991834e-01 9.991834e-01 9.991834e-01 9.991834e-01 [11]: 9.991834e-01 9.991834e-01 9.991834e-01 9.991834e-01 9.991834e-01 [16]: 9.991834e-01 9.991834e-01 9.991834e-01 9.991834e-01 9.595180e-01 [21]: 9.991834e-01 9.991834e-01 9.991834e-01 9.991834e-01 9.991834e-01 [26]: 9.991834e-01 9.991834e-01 4.351895e-01 9.991834e-01 9.766522e-01 [31]: 1.414256e-01 1.304340e-02 3.530937e-01 6.887709e-01 2.385602e-02 [36]: 1.322457e-02 2.722920e-01 5.426136e-02 4.218158e-01 1.581127e-02 [41]: 8.825610e-04 1.252228e-05 8.743649e-01 3.016908e-01 3.516461e-02 [46]: 9.582456e-03 3.877222e-02 3.172920e-03 8.122276e-01 1.950067e-01