P-value correction: Difference between revisions
Content deleted Content added
Added Go |
|||
Line 3,340: | Line 3,340: | ||
[40] 0.0008821796 0.0000125222 0.6357389664 0.2889497995 0.0362651575 |
[40] 0.0008821796 0.0000125222 0.6357389664 0.2889497995 0.0362651575 |
||
[45] 0.0101847015 0.0389807074 0.0031679962 0.5985019850 0.1963376344</pre> |
[45] 0.0101847015 0.0389807074 0.0031679962 0.5985019850 0.1963376344</pre> |
||
=={{header|Phix}}== |
|||
Translation of Kotlin (version 2), except for the Hommel part, which is translated from Go.<br> |
|||
Note that sq_min(), extract(), and custom_sort() as used below require 0.8.0+ |
|||
<lang Phix>enum UP, DOWN |
|||
function ratchet(sequence p, integer direction) |
|||
atom m = p[1] |
|||
for i=1 to length(p) do |
|||
if iff(direction=UP?p[i]>m:p[i]<m) then p[i] = m end if |
|||
m = p[i] |
|||
end for |
|||
return sq_min(p,1) |
|||
end function |
|||
function schwartzian(sequence p, mult, integer direction) |
|||
sequence order = custom_sort(p,tagset(length(p))) |
|||
if direction=UP then order = reverse(order) end if |
|||
sequence pa = ratchet(sq_mul(mult,extract(p,order)), direction) |
|||
return extract(pa,order,invert:=true) |
|||
end function |
|||
function adjust(sequence p, string method) |
|||
integer size = length(p) |
|||
sequence mult = tagset(size) |
|||
switch method |
|||
case "Benjamini-Hochberg": |
|||
mult = sq_div(size,sq_sub(size+1,mult)) |
|||
return schwartzian(p, mult, UP) |
|||
case "Benjamini-Yekutieli": |
|||
atom q = sum(sq_div(1,mult)) |
|||
mult = sq_div(q*size,sq_sub(size+1,mult)) |
|||
return schwartzian(p, mult, UP) |
|||
case "Bonferroni": |
|||
return sq_min(sq_mul(p,size),1) |
|||
case "Hochberg": |
|||
return schwartzian(p, mult, UP) |
|||
case "Holm": |
|||
mult = sq_sub(size+1,mult) |
|||
return schwartzian(p, mult, DOWN) |
|||
case "Hommel": |
|||
sequence ivdx = repeat(0,size) |
|||
for i=1 to size do ivdx[i] = {p[i],i} end for |
|||
ivdx = sort(ivdx) |
|||
sequence s = vslice(ivdx,1), |
|||
m = sq_div(sq_mul(s,size),mult), |
|||
{q,pa} @= repeat(min(m),size), |
|||
order = vslice(ivdx,2) |
|||
for j=size-1 to 2 by -1 do |
|||
sequence lwr = tagset(size-j+1), |
|||
upr = sq_add(size-j+1,tagset(j-1)) |
|||
atom qmin = j*s[upr[1]]/2 |
|||
for i=2 to length(upr) do |
|||
qmin = min(s[upr[i]]*j/(i+1),qmin) |
|||
end for |
|||
for i=1 to length(lwr) do |
|||
q[lwr[i]] = min(s[lwr[i]]*j, qmin) |
|||
end for |
|||
for i=1 to length(upr) do |
|||
q[upr[i]] = q[size-j+1] |
|||
end for |
|||
pa = sq_max(pa,q) |
|||
end for |
|||
return extract(pa,order,invert:=true) |
|||
case "Sidak": |
|||
for i=1 to length(p) do |
|||
p[i] = 1 - power(1-p[i],size) |
|||
end for |
|||
return p |
|||
else |
|||
return {} -- (unknown method) |
|||
end switch |
|||
return p |
|||
end function |
|||
constant {types,correct_answers} = columnize({ |
|||
{"Benjamini-Hochberg", |
|||
{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-Yekutieli", |
|||
{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}}, |
|||
{"Bonferroni", |
|||
{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}}, |
|||
{"Hochberg", |
|||
{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}}, |
|||
{"Holm", |
|||
{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}}, |
|||
{"Hommel", |
|||
{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}}, |
|||
{"Sidak", |
|||
{1.0000000000, 1.0000000000, 0.9946598274, 0.9914285749, 0.9999515274, |
|||
1.0000000000, 0.9999999688, 1.0000000000, 1.0000000000, 1.0000000000, |
|||
1.0000000000, 1.0000000000, 0.9999999995, 1.0000000000, 0.9999998801, |
|||
1.0000000000, 1.0000000000, 1.0000000000, 0.9999999855, 0.9231179729, |
|||
0.9999999956, 1.0000000000, 1.0000000000, 0.9999317605, 1.0000000000, |
|||
0.9983109511, 1.0000000000, 0.5068253940, 1.0000000000, 0.9703301333, |
|||
0.1832692440, 0.0150545753, 0.4320729669, 0.6993672225, 0.0286818157, |
|||
0.0152621104, 0.3391808707, 0.0656206307, 0.4959194266, 0.0186503726, |
|||
0.0009001752, 0.0000125222, 0.8142104886, 0.3772612062, 0.0430222116, |
|||
0.0108312558, 0.0473319661, 0.0032997780, 0.7705015898, 0.2499384839}}}) |
|||
-- {"Unknown",{1}}}) |
|||
constant 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} |
|||
if length(pValues)=0 or min(pValues)<0 or max(pValues)>1 then |
|||
crash("p-values must be in range 0.0 to 1.0") |
|||
end if |
|||
for i=1 to length(types) do |
|||
string ti = types[i] |
|||
sequence res = adjust(pValues,ti) |
|||
if res={} then |
|||
printf(1,"\nSorry, do not know how to do %s correction.\n"& |
|||
"Perhaps you want one of these?:\n %s\n", |
|||
{ti,join(types[1..$-1],"\n ")}) |
|||
exit |
|||
end if |
|||
-- printf(1,"%s\n",{ti}) |
|||
-- res = correct_answers[i] -- (for easier comparison only) |
|||
-- pp(res,{pp_FltFmt,"%13.10f",pp_IntFmt,"%13.10f",pp_Maxlen,75,pp_Pause,0}) |
|||
atom error = sum(sq_abs(sq_sub(res,correct_answers[i]))) |
|||
printf(1,"%s has cumulative error of %g\n", {ti,error}) |
|||
end for</lang> |
|||
{{out}} |
|||
Matches Kotlin (etc) when some of those lines just above are uncommented. |
|||
<pre> |
|||
Benjamini-Hochberg has cumulative error of 8.03052e-7 |
|||
Benjamini-Yekutieli has cumulative error of 3.64071e-7 |
|||
Bonferroni has cumulative error of 6.5e-8 |
|||
Hochberg has cumulative error of 2.7375e-7 |
|||
Holm has cumulative error of 2.8095e-7 |
|||
Hommel has cumulative error of 4.35302e-7 |
|||
Sidak has cumulative error of 7.26897e-10 |
|||
</pre> |
|||
=={{header|Python}}== |
=={{header|Python}}== |