P-value correction: Difference between revisions

Line 3,340:

[40] 0.0008821796 0.0000125222 0.6357389664 0.2889497995 0.0362651575

[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, 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, 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, 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, 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.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>

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}}==