P-value correction

########################### 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] ; # 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] )»[1]]
}

multi adjust ( @p, $unknown ) {
    note "\nSorry, do not know how to do $unknown correction.\n" ~
    "Perhaps you want one of these?:\n" ~
    <Benjamini-Hochberg Benjamini-Yekutieli Bonferroni Hochberg
    Holm Hommel Šidák>.join("\n");
    exit
}

########################### 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, $_
}

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

Last updated