P-value correction: Difference between revisions
Content deleted Content added
→{{header|Kotlin}}: Added second version to deal with licensing issues of first version. |
→{{header|C}}: Added second version to avoid licensing issues of first version. |
||
Line 26: | Line 26: | ||
=={{header|C}}== |
=={{header|C}}== |
||
===Version 1=== |
|||
{{works with|C99}} |
{{works with|C99}} |
||
{{trans|R}} |
{{trans|R}} |
||
Line 677: | Line 678: | ||
type 5 = 'hommel' has cumulative error of 4.35302e-07</pre> |
type 5 = 'hommel' has cumulative error of 4.35302e-07</pre> |
||
===Version 2=== |
|||
{{works with|C89}} |
|||
{{trans|Kotlin}} |
|||
To avoid licensing issues, this version is a translation of the Kotlin entry (Version 2) which is itself a partial translation of the Perl 6 entry. If using gcc, you need to link to the math library (-lm). |
|||
<lang C>#include <stdio.h> |
|||
#include <stdlib.h> |
|||
#include <math.h> |
|||
#define SIZE 50 |
|||
#define each_i(start, end) for (i = start; i < end; ++i) |
|||
typedef enum { UP, DOWN } direction; |
|||
typedef struct { int index; double value; } iv1; |
|||
typedef struct { int index; int value; } iv2; |
|||
/* test also for 'Unknown' correction type */ |
|||
const char *types[8] = { |
|||
"Benjamini-Hochberg", "Benjamini-Yekutieli", "Bonferroni", "Hochberg", |
|||
"Holm", "Hommel", "Šidák", "Unknown" |
|||
}; |
|||
int compare_iv1(const void *a, const void *b) { |
|||
double aa = ((iv1 *)a) -> value; |
|||
double bb = ((iv1 *)b) -> value; |
|||
if (aa > bb) return 1; |
|||
if (aa < bb) return -1; |
|||
return 0; |
|||
} |
|||
int compare_iv1_desc(const void *a, const void *b) { |
|||
return -compare_iv1(a, b); |
|||
} |
|||
int compare_iv2(const void *a, const void *b) { |
|||
return ((iv2 *)a) -> value - ((iv2 *)b) -> value; |
|||
} |
|||
void ratchet(double *pa, direction dir) { |
|||
int i; |
|||
double m = pa[0]; |
|||
if (dir == UP) { |
|||
each_i(1, SIZE) { |
|||
if (pa[i] > m) pa[i] = m; |
|||
m = pa[i]; |
|||
} |
|||
} |
|||
else { |
|||
each_i(1, SIZE) { |
|||
if (pa[i] < m) pa[i] = m; |
|||
m = pa[i]; |
|||
} |
|||
} |
|||
each_i(0, SIZE) if (pa[i] > 1.0) pa[i] = 1.0; |
|||
} |
|||
void schwartzian(const double *p, double *pa, direction dir) { |
|||
int i; |
|||
int order[SIZE]; |
|||
int order2[SIZE]; |
|||
iv1 iv1s[SIZE]; |
|||
iv2 iv2s[SIZE]; |
|||
double pa2[SIZE]; |
|||
each_i(0, SIZE) { iv1s[i].index = i; iv1s[i].value = p[i]; } |
|||
if (dir == UP) |
|||
qsort(iv1s, SIZE, sizeof(iv1s[0]), compare_iv1_desc); |
|||
else |
|||
qsort(iv1s, SIZE, sizeof(iv1s[0]), compare_iv1); |
|||
each_i(0, SIZE) order[i] = iv1s[i].index; |
|||
each_i(0, SIZE) pa[i] *= p[order[i]]; |
|||
ratchet(pa, dir); |
|||
each_i(0, SIZE) { iv2s[i].index = i; iv2s[i].value = order[i]; } |
|||
qsort(iv2s, SIZE, sizeof(iv2s[0]), compare_iv2); |
|||
each_i(0, SIZE) order2[i] = iv2s[i].index; |
|||
each_i(0, SIZE) pa2[i] = pa[order2[i]]; |
|||
each_i(0, SIZE) pa[i] = pa2[i]; |
|||
} |
|||
void adjust(const double *p, double *pa, const char *type) { |
|||
int i; |
|||
if (!strcmp(type, "Benjamini-Hochberg")) { |
|||
each_i(0, SIZE) pa[i] = (double)SIZE / (SIZE - i); |
|||
schwartzian(p, pa, UP); |
|||
} |
|||
else if (!strcmp(type, "Benjamini-Yekutieli")) { |
|||
double q = 0.0; |
|||
each_i(1, SIZE + 1) q += 1.0 / i; |
|||
each_i(0, SIZE) pa[i] = q * SIZE / (SIZE - i); |
|||
schwartzian(p, pa, UP); |
|||
} |
|||
else if (!strcmp(type, "Bonferroni")) { |
|||
each_i(0, SIZE) pa[i] = (p[i] * SIZE > 1.0) ? 1.0 : p[i] * SIZE; |
|||
} |
|||
else if (!strcmp(type, "Hochberg")) { |
|||
each_i(0, SIZE) pa[i] = i + 1.0; |
|||
schwartzian(p, pa, UP); |
|||
} |
|||
else if (!strcmp(type, "Holm")) { |
|||
each_i(0, SIZE) pa[i] = SIZE - i; |
|||
schwartzian(p, pa, DOWN); |
|||
} |
|||
else if (!strcmp(type, "Hommel")) { |
|||
int i, j; |
|||
int order[SIZE]; |
|||
int order2[SIZE]; |
|||
iv1 iv1s[SIZE]; |
|||
iv2 iv2s[SIZE]; |
|||
double s[SIZE]; |
|||
double q[SIZE]; |
|||
double pa2[SIZE]; |
|||
int indices[SIZE]; |
|||
each_i(0, SIZE) { iv1s[i].index = i; iv1s[i].value = p[i]; } |
|||
qsort(iv1s, SIZE, sizeof(iv1s[0]), compare_iv1); |
|||
each_i(0, SIZE) order[i] = iv1s[i].index; |
|||
each_i(0, SIZE) s[i] = p[order[i]]; |
|||
double min = s[0] * SIZE; |
|||
each_i(1, SIZE) { |
|||
double temp = s[i] / (i + 1.0); |
|||
if (temp < min) min = temp; |
|||
} |
|||
each_i(0, SIZE) q[i] = min; |
|||
each_i(0, SIZE) pa2[i] = min; |
|||
for (j = SIZE - 1; j >= 2; --j) { |
|||
each_i(0, SIZE) indices[i] = i; |
|||
int upper_start = SIZE - j + 1; /* upper indices start index */ |
|||
int upper_size = j - 1; /* size of upper indices */ |
|||
int lower_size = SIZE - upper_size; /* size of lower indices */ |
|||
double qmin = j * s[indices[upper_start]] / 2.0; |
|||
each_i(1, upper_size) { |
|||
double temp = s[indices[upper_start + i]] * j / (2.0 + i); |
|||
if (temp < qmin) qmin = temp; |
|||
} |
|||
each_i(0, lower_size) { |
|||
double temp = s[indices[i]] * j; |
|||
q[indices[i]] = (temp < qmin) ? temp : qmin; |
|||
} |
|||
each_i(0, upper_size) q[indices[upper_start + i]] = q[SIZE - j]; |
|||
each_i(0, SIZE) if (pa2[i] < q[i]) pa2[i] = q[i]; |
|||
} |
|||
each_i(0, SIZE) { iv2s[i].index = i; iv2s[i].value = order[i]; } |
|||
qsort(iv2s, SIZE, sizeof(iv2s[0]), compare_iv2); |
|||
each_i(0, SIZE) order2[i] = iv2s[i].index; |
|||
each_i(0, SIZE) pa[i] = pa2[order2[i]]; |
|||
} |
|||
else if (!strcmp(type, "Šidák")) { |
|||
each_i(0, SIZE) pa[i] = 1.0 - pow(1.0 - p[i], SIZE); |
|||
} |
|||
else { |
|||
printf("\nSorry, do not know how to do '%s' correction.\n", type); |
|||
printf("Perhaps you want one of these?:\n"); |
|||
each_i(0, 7) printf(" %s\n", types[i]); |
|||
exit(1); |
|||
} |
|||
} |
|||
void adjusted(const double *p, const char *type) { |
|||
int i; |
|||
double pa[SIZE] = { 0.0 }; |
|||
if (check(p)) { |
|||
adjust(p, pa, type); |
|||
printf("\n%s", type); |
|||
each_i(0, SIZE) { |
|||
if (!(i % 5)) printf("\n[%2d] ", i); |
|||
printf("%1.10f ", pa[i]); |
|||
} |
|||
printf("\n"); |
|||
} |
|||
else { |
|||
printf("p-values must be in range 0.0 to 1.0\n"); |
|||
exit(1); |
|||
} |
|||
} |
|||
int check(const double* p) { |
|||
int i; |
|||
each_i(0, SIZE) { |
|||
if (p[i] < 0.0 || p[i] > 1.0) return 0; |
|||
} |
|||
return 1; |
|||
} |
|||
int main() { |
|||
int i; |
|||
double p_values[SIZE] = { |
|||
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 |
|||
}; |
|||
each_i(0, 8) adjusted(p_values, types[i]); |
|||
return 0; |
|||
}</lang> |
|||
{{output}} |
|||
<pre> |
|||
Same as Kotlin (Version 2) output. |
|||
</pre> |
|||
=={{header|D}}== |
=={{header|D}}== |