Given a list of p-values, adjust the p-values for multiple comparisons. This is done in order to control the false positive, or Type 1 error rate.

This is also known as the "false discovery rate" (FDR). After adjustment, the p-values will be higher but still inside [0,1].

The adjusted p-values are sometimes called "q-values".



Task

Given one list of p-values, return the p-values correcting for multiple comparisons

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



There are several methods to do this, see:

 Yoav Benjamini, Yosef Hochberg "Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing", Journal of the Royal Statistical Society. Series B, Vol. 57, No. 1 (1995), pp. 289-300, JSTOR:2346101
 Yoav Benjamini, Daniel Yekutieli, "The control of the false discovery rate in multiple testing under dependency", Ann. Statist., Vol. 29, No. 4 (2001), pp. 1165-1188, DOI:10.1214/aos/1013699998 JSTOR:2674075
 Sture Holm, "A Simple Sequentially Rejective Multiple Test Procedure", Scandinavian Journal of Statistics, Vol. 6, No. 2 (1979), pp. 65-70, JSTOR:4615733
 Yosef Hochberg, "A sharper Bonferroni procedure for multiple tests of significance", Biometrika, Vol. 75, No. 4 (1988), pp 800–802, DOI:10.1093/biomet/75.4.800 JSTOR:2336325
 Gerhard Hommel, "A stagewise rejective multiple test procedure based on a modified Bonferroni test", Biometrika, Vol. 75, No. 2 (1988), pp 383–386, DOI:10.1093/biomet/75.2.383 JSTOR:2336190


Each method has its own advantages and disadvantages.




with javascript_semantics
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),
                     qh = repeat(min(m),size),
                     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
                    qh[lwr[i]] = min(s[lwr[i]]*j, qmin)
                end for
                for i=1 to length(upr) do
                    qh[upr[i]] = qh[size-j+1]
                end for
                pa = sq_max(pa,qh)
            end for
            return extract(pa,order,invert:=true)
 
        case "Sidak":
            p = deep_copy(p)
            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
