P-value correction: Difference between revisions

m (→‎{{header|Phix}}: syntax coloured, made p2js compatible)
 
(3 intermediate revisions by 3 users not shown)
Line 53:
 
Link with <code>-lm</code>
<langsyntaxhighlight Clang="c">#include <stdio.h>//printf
#include <stdlib.h>//qsort
#include <math.h>//fabs
Line 601:
return 0;
}
</syntaxhighlight>
</lang>
 
{{out}}
Line 693:
{{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 Raku entry. If using gcc, you need to link to the math library (-lm).
<langsyntaxhighlight Clang="c">#include <stdio.h>
#include <stdlib.h>
#include <math.h>
Line 887:
each_i(0, 8) adjusted(p_values, types[i]);
return 0;
}</langsyntaxhighlight>
 
{{output}}
Line 896:
=={{header|C sharp|C#}}==
{{trans|Java}}
<langsyntaxhighlight lang="csharp">using System;
using System.Collections.Generic;
using System.Linq;
Line 1,227:
}
}
}</langsyntaxhighlight>
{{out}}
<pre>[ 1] 6.126681E-001 8.521710E-001 1.987205E-001 1.891595E-001 3.217789E-001
Line 1,309:
=={{header|C++}}==
{{trans|Java}}
<langsyntaxhighlight lang="cpp">#include <algorithm>
#include <functional>
#include <iostream>
Line 1,630:
 
return 0;
}</langsyntaxhighlight>
{{out}}
<pre>[ 1] 0.6126681081 0.8521710465 0.1987205200 0.1891595417 0.3217789286
Line 1,718:
{{trans|Kotlin}}
''This work is based on R source code covered by the '''GPL''' license. It is thus a modified version, also covered by the GPL. See the [https://www.gnu.org/licenses/gpl-faq.html#GPLRequireSourcePostedPublic FAQ about GNU licenses]''.
<langsyntaxhighlight Dlang="d">import std.algorithm;
import std.conv;
import std.math;
Line 2,057:
writefln("\ntype %d = '%s' has a cumulative error of %g", type, types[type], error);
}
}</langsyntaxhighlight>
{{out}}
<pre>[ 1] 6.126681e-01 0.8521710465 0.1987205200 0.1891595417 0.3217789286
Line 2,140:
=={{header|Go}}==
{{trans|Kotlin (Version 2)}}
<langsyntaxhighlight lang="go">package main
 
import (
Line 2,443:
fmt.Println(s)
}
}</langsyntaxhighlight>
 
{{out}}
Line 2,546:
{{works with|Java|8}}
''This work is based on R source code covered by the '''GPL''' license. It is thus a modified version, also covered by the GPL. See the [https://www.gnu.org/licenses/gpl-faq.html#GPLRequireSourcePostedPublic FAQ about GNU licenses]''.
<langsyntaxhighlight Javalang="java">import java.util.Arrays;
import java.util.Comparator;
 
Line 2,895:
}
}
}</langsyntaxhighlight>
{{out}}
<pre>[ 1] 6.126681e-01 0.8521710465 0.1987205200 0.1891595417 0.3217789286
Line 2,975:
 
type 5 = 'hommel' has a cumulative error of 4.35302e-07</pre>
 
=={{header|jq}}==
'''Adapted from [[#Wren|Wren]]'''
 
'''Works with jq, the C implementation of jq'''
 
'''Works with gojq, the Go implementation of jq'''
 
'''Works with jaq, the Rust implementation of jq'''
 
The def of `_nwise` is included for the sake of gojq; it may be omitted if using jq or jaq.
<syntaxhighlight lang="jq">
### For gojq
# Require $n > 0
def nwise($n):
def _n: if length <= $n then . else .[:$n] , (.[$n:] | _n) end;
if $n <= 0 then "nwise: argument should be non-negative" else _n end;
 
### Generic functions
 
def array($n): . as $in | [range(0;$n)|$in];
 
def lpad($len): tostring | ($len - length) as $l | (" " * $l) + .;
 
def rpad($len): tostring | ($len - length) as $l | . + (" " * $l);
 
def round($ndec): pow(10;$ndec) as $p | . * $p | round / $p;
 
# tabular print
def tprint($columns; $width):
reduce _nwise($columns) as $row ("";
. + ($row|map(lpad($width)) | join(" ")) + "\n" );
 
# Emit the permutation p such that [range(0;length) as $i | .[$p[$i]]] is sorted
def sort_index:
[range(0;length) as $i | [$i, .[$i]]]
| sort_by(.[1])
| map(.[0]);
 
 
### p-value Corrections
 
def types: [
"Benjamini-Hochberg", "Benjamini-Yekutieli", "Bonferroni", "Hochberg",
"Holm", "Hommel", "Šidák"
];
 
######################################
# The functions in this section expect
# an array of p-values as input.
######################################
 
def pFormat($cols):
map(round(10) | rpad(12)) | tprint($cols; 12);
 
def check:
if (length == 0 or min < 0 or max > 1)
then "p-values must be in the range 0 to 1 inclusive" | error
else .
end;
 
# $dir should be "UP" or "DOWN"
def ratchet($dir):
{ m: .[0], p: .}
| if $dir == "UP"
then reduce range(1; .p|length) as $i (.;
if (.p[$i] > .m) then .p[$i] = .m end
| .m = .p[$i])
else reduce range(1; .p|length) as $i (.;
if (.p[$i] < .m) then .p[$i] = .m end
| .m = .p[$i] )
end
| .p
| map( if . < 1 then . else 1 end);
 
# If $dir is "UP" then reverse is called
def schwartzian($mult; $dir):
length as $size
| (sort_index | if $dir == "UP" then reverse else . end) as $order
| ([range(0;$size) as $i | $mult[$i] * .[$order[$i]] ]
| ratchet($dir)) as $pa
| ($order | sort_index) as $order2
| [ range(0; $size) as $i | $pa[$order2[$i]]] ;
 
# $type should be one of `types`
def adjust($type):
length as $size
| if $size == 0 then "The array of p-values cannot be empty." | error end
| if $type == "Benjamini-Hochberg"
then
[range(0;$size) as $i | $size / ($size - $i)] as $mult
| schwartzian($mult; "UP")
elif $type == "Benjamini-Yekutieli"
then (reduce range(1; 1+$size) as $i (0; . + (1/$i))) as $q
| [range(0; $size) as $i | $q * $size / ($size - $i)] as $mult
| schwartzian($mult; "UP")
 
elif $type == "Bonferroni"
then map( [(. * $size), 1] | min)
 
elif $type == "Hochberg"
then
[range(0;$size) as $i | $i + 1] as $mult
| schwartzian($mult; "UP")
 
elif $type == "Holm"
then
[range(0; $size) as $i | $size - $i] as $mult
| schwartzian($mult; "DOWN")
 
elif $type == "Hommel"
then
sort_index as $order
| [range(0; $size) as $i | .[$order[$i]]] as $s
| [range(0; $size) as $i | $s[$i] * $size / ($i + 1)] as $m
| ($m | min) as $min
| { q: ($min | array($size)),
pa: ($min | array($size)) }
| reduce range($size-1; 1; -1) as $j (.;
.lower = (0 | array($size - $j + 1)) # lower indices
| reduce range(0; .lower|length) as $i (.; .lower[$i] = $i)
| .upper = (0|array($j - 1))
| reduce range(0; .upper|length) as $i (.; .upper[$i] = $size - $j + 1 + $i)
| .qmin = ($j * $s[.upper[0]] / 2)
| reduce range(1; .upper|length) as $i (.;
($s[.upper[$i]] * $j / (2 + $i)) as $temp
| if $temp < .qmin then .qmin = $temp end )
| reduce range(0; .lower|length) as $i (.;
.q[.lower[$i]] = ([.qmin, ($s[.lower[$i]] * $j)] | min) )
| reduce range(0; .upper|length) as $i (.; .q[.upper[$i]] = .q[$size - $j])
| reduce range(0; $size) as $i (.; if (.pa[$i] < .q[$i]) then .pa[$i] = .q[$i] end)
)
| ($order | sort_index) as $order2
| [range(0; $size) as $i | .pa[$order2[$i] ]]
 
elif $type == "Šidák"
then map(1 - pow(1 - .; $size) )
 
else
"\nSorry, do not know how to do '\($type)' correction.\n" +
"Perhaps you want one of the following?\n" +
(types | map( " \(.)" ) | join("\n") )
end;
 
def adjusted($type):
"\n\($type)",
(check | adjust($type) | pFormat(5));
 
### Example
 
def 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
];
 
pValues | adjusted( types[] )
</syntaxhighlight>
{{output}}
The output shown here is from a run using jq. The output using gojq
is the same except that numbers are presented without using scientific notation.
<pre>
Benjamini-Hochberg
0.6126681081 0.8521710465 0.19872052 0.1891595417 0.3217789286
0.930145 0.487037 0.930145 0.6049730556 0.6826752564
0.6482628947 0.72537225 0.5280972727 0.8769925556 0.4705703448
0.9241867391 0.6049730556 0.7856107317 0.4887525806 0.1136717045
0.4991890625 0.8769925556 0.9991834 0.3217789286 0.930145
0.2304957692 0.5832475 0.0389954722 0.8521710465 0.1476842609
0.016836375 0.0025629017 0.0351608437 0.0625018947 0.0036365888
0.0025629017 0.0294688286 0.0061660636 0.0389954722 0.0026889914
0.0004502863 1.25223e-05 0.0788155476 0.03142613 0.004846527
0.0025629017 0.004846527 0.0011017083 0.072520325 0.0220595769
 
 
Benjamini-Yekutieli
1 1 0.8940844244 0.8510676197 1
1 1 1 1 1
1 1 1 1 1
1 1 1 1 0.5114323399
1 1 1 1 1
1 1 0.1754486368 1 0.6644618149
0.0757503083 0.0115310209 0.1581958559 0.2812088585 0.0163617595
0.0115310209 0.1325863108 0.0277423864 0.1754486368 0.0120983246
0.0020259303 5.63403e-05 0.3546073326 0.1413926119 0.0218055202
0.0115310209 0.0218055202 0.004956812 0.3262838334 0.0992505663
 
 
Bonferroni
1 1 1 1 1
1 1 1 1 1
1 1 1 1 1
1 1 1 1 1
1 1 1 1 1
1 1 0.7019185 1 1
0.2020365 0.015166745 0.5625735 1 0.02909271
0.01537741 0.4125636 0.0678267 0.680348 0.01882294
0.0009005725 1.25223e-05 1 0.47139195 0.043955765
0.010889155 0.04846527 0.003305125 1 0.2867745
 
 
Hochberg
0.9991834 0.9991834 0.9991834 0.9991834 0.9991834
0.9991834 0.9991834 0.9991834 0.9991834 0.9991834
0.9991834 0.9991834 0.9991834 0.9991834 0.9991834
0.9991834 0.9991834 0.9991834 0.9991834 0.9991834
0.9991834 0.9991834 0.9991834 0.9991834 0.9991834
0.9991834 0.9991834 0.46326621 0.9991834 0.9991834
0.15758847 0.013839669 0.39380145 0.76002304 0.0250197306
0.013839669 0.305297064 0.05426136 0.46263664 0.0165641872
0.0008825611 1.25223e-05 0.9930759 0.339402204 0.0369228426
0.0102358057 0.0397415214 0.00317292 0.89925203 0.21794862
 
 
Holm
1 1 1 1 1
1 1 1 1 1
1 1 1 1 1
1 1 1 1 1
1 1 1 1 1
1 1 0.46326621 1 1
0.15758847 0.0139534054 0.39380145 0.76002304 0.0250197306
0.0139534054 0.305297064 0.05426136 0.46263664 0.0165641872
0.0008825611 1.25223e-05 0.9930759 0.339402204 0.0369228426
0.0102358057 0.0397415214 0.00317292 0.89925203 0.21794862
 
 
Hommel
0.9991834 0.9991834 0.9991834 0.99876238 0.9991834
0.9991834 0.9991834 0.9991834 0.9991834 0.9991834
0.9991834 0.9991834 0.9991834 0.9991834 0.9991834
0.9991834 0.9991834 0.9991834 0.9991834 0.959518
0.9991834 0.9991834 0.9991834 0.9991834 0.9991834
0.9991834 0.9991834 0.43518947 0.9991834 0.97665225
0.14142555 0.0130434007 0.3530936533 0.68877088 0.0238560222
0.0132245726 0.272291976 0.05426136 0.42181576 0.0158112696
0.0008825611 1.25223e-05 0.8743649143 0.301690848 0.035164612
0.0095824564 0.038772216 0.00317292 0.81222764 0.19500666
 
 
Šidák
1 1 0.9946598274 0.9914285749 0.9999515274
1 0.9999999688 1 1 1
1 1 0.9999999995 1 0.9999998801
1 1 1 0.9999999855 0.9231179729
0.9999999956 1 1 0.9999317605 1
0.9983109511 1 0.506825394 1 0.9703301333
0.183269244 0.0150545753 0.4320729669 0.6993672225 0.0286818157
0.0152621104 0.3391808707 0.0656206307 0.4959194266 0.0186503726
0.0009001752 1.25222e-05 0.8142104886 0.3772612062 0.0430222116
0.0108312558 0.0473319661 0.003299778 0.7705015898 0.2499384839
</pre>
 
=={{header|Julia}}==
<langsyntaxhighlight lang="julia">using MultipleTesting, IterTools, Printf
 
p = [4.533744e-01, 7.296024e-01, 9.936026e-02, 9.079658e-02, 1.801962e-01,
Line 3,001 ⟶ 3,261:
println("\n", corr)
printpvalues(adjust(p, corr))
end</langsyntaxhighlight>
 
{{out}}
Line 3,058 ⟶ 3,318:
''This work is based on R source code covered by the '''GPL''' license. It is thus a modified version, also covered by the GPL. See the [https://www.gnu.org/licenses/gpl-faq.html#GPLRequireSourcePostedPublic FAQ about GNU licenses]''.
 
<langsyntaxhighlight lang="scala">// version 1.1.51
 
import java.util.Arrays
Line 3,339 ⟶ 3,599:
println(f.format(type, types[type], error))
}
}</langsyntaxhighlight>
 
{{out}}
Line 3,427 ⟶ 3,687:
 
To avoid licensing issues, this version follows the approach of the Raku entry of which it is a partial translation. However, the correction routines themselves have been coded independently, common code factored out into separate functions (analogous to Raku) and (apart from the Šidák method) agree with the Raku results.
<langsyntaxhighlight lang="scala">// version 1.2.21
 
typealias DList = List<Double>
Line 3,572 ⟶ 3,832:
 
types.forEach { println(adjusted(pValues, it)) }
}</langsyntaxhighlight>
 
{{out}}
Line 3,604 ⟶ 3,864:
=={{header|Nim}}==
{{trans|Kotlin (Version 2)}}
<langsyntaxhighlight Nimlang="nim">import algorithm, math, sequtils, strformat, strutils, sugar
 
type
Line 3,758 ⟶ 4,018:
 
for ctype in CorrectionType:
echo adjusted(PVals, ctype)</langsyntaxhighlight>
 
{{out}}
Line 3,849 ⟶ 4,109:
{{trans|C}}
''This work is based on R source code covered by the '''GPL''' license. It is thus a modified version, also covered by the GPL. See the [https://www.gnu.org/licenses/gpl-faq.html#GPLRequireSourcePostedPublic FAQ about GNU licenses]''.
<langsyntaxhighlight lang="perl">#!/usr/bin/env perl
 
use strict;
Line 4,159 ⟶ 4,419:
printf("type $method has cumulative error of %g.\n", $error);
}
</syntaxhighlight>
</lang>
 
{{out}}
Line 4,179 ⟶ 4,439:
=={{header|Phix}}==
Translation of Kotlin (version 2), except for the Hommel part, which is translated from Go.
<!--<langsyntaxhighlight Phixlang="phix">(phixonline)-->
<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span>
<span style="color: #008080;">enum</span> <span style="color: #000000;">UP</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">DOWN</span>
Line 4,373 ⟶ 4,633:
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"%s has cumulative error of %g\n"</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">ti</span><span style="color: #0000FF;">,</span><span style="color: #000000;">error</span><span style="color: #0000FF;">})</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<!--</langsyntaxhighlight>-->
{{out}}
Matches Kotlin (etc) when some of those lines just above are uncommented.
Line 4,389 ⟶ 4,649:
{{trans|Perl}}
''This work is based on R source code covered by the '''GPL''' license. It is thus a modified version, also covered by the GPL. See the [https://www.gnu.org/licenses/gpl-faq.html#GPLRequireSourcePostedPublic FAQ about GNU licenses]''.
<langsyntaxhighlight lang="python">from __future__ import division
import sys
 
Line 4,624 ⟶ 4,884:
error += abs(q[i] - correct_answers[key][i])
print '%s error = %g' % (key.upper(), error)
</syntaxhighlight>
</lang>
 
{{out}}
Line 4,639 ⟶ 4,899:
The '''p.adjust''' function is built-in, see [https://stat.ethz.ch/R-manual/R-devel/library/stats/html/p.adjust.html R manual].
 
<langsyntaxhighlight Rlang="rsplus">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,
Line 4,667 ⟶ 4,927:
 
p.adjust(p, method = 'hommel')
writeLines("Hommel\n")</langsyntaxhighlight>
 
{{out}}
Line 4,739 ⟶ 4,999:
{{works with|Rakudo|2019.03.1}}
 
<syntaxhighlight lang="raku" perl6line>########################### Helper subs ###########################
 
sub adjusted (@p, $type) { "\n$type\n" ~ format adjust( check(@p), $type ) }
Line 4,832 ⟶ 5,092:
{
say adjusted @p-values, $_
}</langsyntaxhighlight>
 
{{out}}
Line 4,921 ⟶ 5,181:
=={{header|Ruby}}==
{{trans|Perl}}
<langsyntaxhighlight lang="ruby">def pmin(array)
x = 1
pmin_array = []
Line 5,185 ⟶ 5,445:
puts "total error for #{method} = #{error}"
end
</syntaxhighlight>
</lang>
{{out}}
<pre>Benjamini-Yekutieli
Line 5,202 ⟶ 5,462:
 
=={{header|Rust}}==
<langsyntaxhighlight lang="rust">
use std::iter;
 
Line 5,435 ⟶ 5,695:
}
}
</syntaxhighlight>
</lang>
{{out}}
<pre style="height:60ex;overflow:scroll;">
Line 5,525 ⟶ 5,785:
=={{header|SAS}}==
 
<langsyntaxhighlight lang="sas">data pvalues;
input raw_p @@;
cards;
Line 5,542 ⟶ 5,802:
 
proc multtest pdata=pvalues bon sid hom hoc holm;
run;</langsyntaxhighlight>
 
'''output'''
Line 5,619 ⟶ 5,879:
First, install the package with:
 
<syntaxhighlight lang ="stata">ssc install qqvalue</langsyntaxhighlight>
 
Given a dataset containing the p-values in a variable, the qqvalue command generates another variable with the adjusted p-values. Here is an example showing the result with all implemented methods:
 
<langsyntaxhighlight lang="stata">clear
 
#delimit ;
Line 5,645 ⟶ 5,905:
}
 
list</langsyntaxhighlight>
 
'''output'''
Line 5,720 ⟶ 5,980:
{{libheader|Wren-math}}
{{libheader|Wren-sort}}
<langsyntaxhighlight ecmascriptlang="wren">import "./dynamic" for Enum
import "./fmt" for Fmt
import "./seq" for Lst
import "./math" for Nums
import "./sort" for Sort
 
var Direction = Enum.create("Direction", ["UP", "DOWN"])
Line 5,874 ⟶ 6,134:
2.177831e-04, 9.693054e-04, 6.610250e-05, 2.900813e-02, 5.735490e-03
]
types.each { |type| System.print(adjusted.call(pValues, type)) }</langsyntaxhighlight>
 
{{out}}
Line 5,885 ⟶ 6,145:
''This work is based on R source code covered by the '''GPL''' license. It is thus a modified version, also covered by the GPL. See the [https://www.gnu.org/licenses/gpl-faq.html#GPLRequireSourcePostedPublic FAQ about GNU licenses]''.
 
<langsyntaxhighlight 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),..
Line 5,949 ⟶ 6,209:
}
psz.pump(List,'wrap(n){ pa[ro[n]] }); // Hommel q-values
}</langsyntaxhighlight>
<langsyntaxhighlight 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,
Line 5,975 ⟶ 6,235:
}
println();
}</langsyntaxhighlight>
{{out}}
<pre style="height:45ex">
2,515

edits