P-value correction: Difference between revisions
→{{header|jq}}
m (→{{header|R}}) |
|||
(2 intermediate revisions by 2 users not shown) | |||
Line 53:
Link with <code>-lm</code>
<
#include <stdlib.h>//qsort
#include <math.h>//fabs
Line 601:
return 0;
}
</syntaxhighlight>
{{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).
<
#include <stdlib.h>
#include <math.h>
Line 887:
each_i(0, 8) adjusted(p_values, types[i]);
return 0;
}</
{{output}}
Line 896:
=={{header|C sharp|C#}}==
{{trans|Java}}
<
using System.Collections.Generic;
using System.Linq;
Line 1,227:
}
}
}</
{{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}}
<
#include <functional>
#include <iostream>
Line 1,630:
return 0;
}</
{{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]''.
<
import std.conv;
import std.math;
Line 2,057:
writefln("\ntype %d = '%s' has a cumulative error of %g", type, types[type], error);
}
}</
{{out}}
<pre>[ 1] 6.126681e-01 0.8521710465 0.1987205200 0.1891595417 0.3217789286
Line 2,140:
=={{header|Go}}==
{{trans|Kotlin (Version 2)}}
<
import (
Line 2,443:
fmt.Println(s)
}
}</
{{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]''.
<
import java.util.Comparator;
Line 2,895:
}
}
}</
{{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}}==
<
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</
{{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]''.
<
import java.util.Arrays
Line 3,339 ⟶ 3,599:
println(f.format(type, types[type], error))
}
}</
{{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.
<
typealias DList = List<Double>
Line 3,572 ⟶ 3,832:
types.forEach { println(adjusted(pValues, it)) }
}</
{{out}}
Line 3,604 ⟶ 3,864:
=={{header|Nim}}==
{{trans|Kotlin (Version 2)}}
<
type
Line 3,758 ⟶ 4,018:
for ctype in CorrectionType:
echo adjusted(PVals, ctype)</
{{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]''.
<
use strict;
Line 4,159 ⟶ 4,419:
printf("type $method has cumulative error of %g.\n", $error);
}
</syntaxhighlight>
{{out}}
Line 4,179 ⟶ 4,439:
=={{header|Phix}}==
Translation of Kotlin (version 2), except for the Hommel part, which is translated from Go.
<!--<
<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>
<!--</
{{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]''.
<
import sys
Line 4,624 ⟶ 4,884:
error += abs(q[i] - correct_answers[key][i])
print '%s error = %g' % (key.upper(), error)
</syntaxhighlight>
{{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].
<
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")</
{{out}}
Line 4,739 ⟶ 4,999:
{{works with|Rakudo|2019.03.1}}
<syntaxhighlight lang="raku"
sub adjusted (@p, $type) { "\n$type\n" ~ format adjust( check(@p), $type ) }
Line 4,832 ⟶ 5,092:
{
say adjusted @p-values, $_
}</
{{out}}
Line 4,921 ⟶ 5,181:
=={{header|Ruby}}==
{{trans|Perl}}
<
x = 1
pmin_array = []
Line 5,185 ⟶ 5,445:
puts "total error for #{method} = #{error}"
end
</syntaxhighlight>
{{out}}
<pre>Benjamini-Yekutieli
Line 5,202 ⟶ 5,462:
=={{header|Rust}}==
<
use std::iter;
Line 5,435 ⟶ 5,695:
}
}
</syntaxhighlight>
{{out}}
<pre style="height:60ex;overflow:scroll;">
Line 5,525 ⟶ 5,785:
=={{header|SAS}}==
<
input raw_p @@;
cards;
Line 5,542 ⟶ 5,802:
proc multtest pdata=pvalues bon sid hom hoc holm;
run;</
'''output'''
Line 5,619 ⟶ 5,879:
First, install the package with:
<syntaxhighlight lang
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:
<
#delimit ;
Line 5,645 ⟶ 5,905:
}
list</
'''output'''
Line 5,720 ⟶ 5,980:
{{libheader|Wren-math}}
{{libheader|Wren-sort}}
<
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)) }</
{{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]''.
<
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
}</
<
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();
}</
{{out}}
<pre style="height:45ex">
|