Rosetta Code

P-value correction

Revision as of 18:34, 1 November 2017 by rosettacode>Hailholyghost (Created page)
(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)

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.

P-value correction is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.

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, Holm, Hochberg, Hommel, Bonferroni. Each of which has its own advantages and disadvantages. This work is a translation of the R source code. In order to confirm that the new functions are working correctly, each value is compared to R's output and a cumulative absolute error is returned.

C

Works with: C99
Translation of: R

<lang C>

  1. include <stdio.h>//printf
  2. include <stdlib.h>//qsort
  3. include <math.h>//fabs
  4. include <stdbool.h>//bool data type
  5. 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; }

unsigned short int string2FDRtype (const char *restrict STRING) { //this functionshould be separate from p_adjust because this function will be //called in multiple places if (strcasecmp(STRING, "BH") == 0) { return 0; } else if (strcasecmp(STRING, "fdr") == 0) { return 0; } else if (strcasecmp(STRING, "by") == 0) { return 1; } else if (strcasecmp(STRING, "Bonferroni") == 0) { return 2; } else if (strcasecmp(STRING, "hochberg") == 0) { return 3; } else if (strcasecmp(STRING, "holm") == 0) { return 4; } else if (strcasecmp(STRING, "hommel") == 0) { return 5; } else { printf("%s doesn't match any accepted cases.\n", STRING); printf("Failed at %s line %u\n", __FILE__, __LINE__); exit(EXIT_FAILURE); } }

double double_array_min(const double *restrict ARRAY, const size_t NO_OF_ARRAY_ELEMENTS) { 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); } double min = ARRAY[0]; for (size_t index = 1; index < NO_OF_ARRAY_ELEMENTS; index++) { if (ARRAY[index] < min) { min = ARRAY[index]; } } return min; }

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); } const unsigned short int TYPE = string2FDRtype(STRING); //--------------------------------------------------------------------------- //--------------------------------------------------------------------------- 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) unsigned int *restrict indices = seq_len(NO_OF_ARRAY_ELEMENTS, NO_OF_ARRAY_ELEMENTS); //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 *restrict npi = malloc(sizeof(double) * NO_OF_ARRAY_ELEMENTS); if (npi == 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++) { npi[index] = (double)NO_OF_ARRAY_ELEMENTS * p[index] / indices[index]; } free(indices); indices = NULL; const double min = double_array_min(npi, NO_OF_ARRAY_ELEMENTS); free(npi); npi = NULL; 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 }

const size_t PI2_LENGTH = I2_LENGTH;

double q1 = j * p[i2[0]] / 2.0; for (size_t i = 1; i < PI2_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 double *restrict n_i = malloc(sizeof(double) * NO_OF_ARRAY_ELEMENTS);//R's n/i if (n_i == NULL) { printf("failed to malloc at %s line %u.\n", __FILE__, __LINE__); perror(""); exit(EXIT_FAILURE); } 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); } n_i[index] = (double)NO_OF_ARRAY_ELEMENTS / (NO_OF_ARRAY_ELEMENTS - index);// n/i simplified }

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++) { cummin_input[index] = n_i[index] * PVALUES[o[index]];//PVALUES[o[index]] is p[o] } } else if (TYPE == 1) {//BY method double q = 0.0; for (size_t index = 1; index < (1+NO_OF_ARRAY_ELEMENTS); index++) { q += (double) 1.0/index; } for (size_t index = 0; index < NO_OF_ARRAY_ELEMENTS; index++) { cummin_input[index] = q * n_i[index] * 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; free(n_i); n_i = 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

R

<lang R>

  1. p.adjust's source code is below
  2. > p.adjust
  3. function (p, method = p.adjust.methods, n = length(p))
  4. {
  5. method <- match.arg(method)
  6. if (method == "fdr")
  7. method <- "BH"
  8. nm <- names(p)
  9. p <- as.numeric(p)
  10. p0 <- setNames(p, nm)
  11. if (all(nna <- !is.na(p)))
  12. nna <- TRUE
  13. p <- p[nna]
  14. lp <- length(p)
  15. stopifnot(n >= lp)
#   if (n <= 1) 
#       return(p0)
  1. if (n == 2 && method == "hommel")
  2. method <- "hochberg"
  3. p0[nna] <- switch(method, bonferroni = pmin(1, n * p), holm = {
  4. i <- seq_len(lp)
#       o <- order(p)
  1. ro <- order(o)
  2. pmin(1, cummax((n - i + 1L) * p[o]))[ro]
  3. }, hommel = {
#       if (n > lp) p <- c(p, rep.int(1, n - lp))
  1. i <- seq_len(n)
#       o <- order(p)
  1. p <- p[o]
  2. ro <- order(o)
  3. q <- pa <- rep.int(min(n * p/i), n)
    1. for (j in (n - 1):2) {
#           ij <- seq_len(n - j + 1)
  1. i2 <- (n - j + 2):n
  2. q1 <- min(j * p[i2]/(2:j))
  3. q[ij] <- pmin(j * p[ij], q1)
  4. q[i2] <- q[n - j + 1]
  5. pa <- pmax(pa, q)
  6. }
  7. pmax(pa, p)[if (lp < n) ro[1:lp] else ro]
  8. }, hochberg = {
#       i <- lp:1L
  1. o <- order(p, decreasing = TRUE)
  2. ro <- order(o)
  3. pmin(1, cummin((n - i + 1L) * p[o]))[ro]
  4. }, BH = {
  5. i <- lp:1L
  6. o <- order(p, decreasing = TRUE)
  7. ro <- order(o)
  8. pmin(1, cummin(n/i * p[o]))[ro]
  9. }, BY = {
  10. i <- lp:1L
  11. o <- order(p, decreasing = TRUE)
  12. ro <- order(o)
  13. q <- sum(1L/(1L:n))
#       pmin(1, cummin(q * n/i * p[o]))[ro]
  1. }, none = p)
  2. p0
  3. }
  4. <bytecode: 0x3a61b88>
  5. <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
Retrieved from "https://rosettacode.org/wiki/P-value_correction?oldid=278799"