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.




package main
 
import (
    "fmt"
    "log"
    "math"
    "os"
    "sort"
    "strconv"
    "strings"
)
 
type pvalues = []float64
 
type iv1 struct {
    index int
    value float64
}
type iv2 struct{ index, value int }
 
type direction int
 
const (
    up direction = iota
    down
)
 
// Test also for 'Unknown' correction type.
var ctypes = []string{
    "Benjamini-Hochberg", "Benjamini-Yekutieli", "Bonferroni", "Hochberg",
    "Holm", "Hommel", "Šidák", "Unknown",
}
 
func minimum(p pvalues) float64 {
    m := p[0]
    for i := 1; i < len(p); i++ {
        if p[i] < m {
            m = p[i]
        }
    }
    return m
}
 
func maximum(p pvalues) float64 {
    m := p[0]
    for i := 1; i < len(p); i++ {
        if p[i] > m {
            m = p[i]
        }
    }
    return m
}
 
func adjusted(p pvalues, ctype string) (string, error) {
    err := check(p)
    if err != nil {
        return "", err
    }
    temp := pformat(adjust(p, ctype), 5)
    return fmt.Sprintf("\n%s\n%s", ctype, temp), nil
}
 
func pformat(p pvalues, cols int) string {
    var lines []string
    for i := 0; i < len(p); i += cols {
        fchunk := p[i : i+cols]
        schunk := make([]string, cols)
        for j := 0; j < cols; j++ {
            schunk[j] = strconv.FormatFloat(fchunk[j], 'f', 10, 64)
        }
        lines = append(lines, fmt.Sprintf("[%2d]  %s", i, strings.Join(schunk, " ")))
    }
    return strings.Join(lines, "\n")
}
 
func check(p []float64) error {
    cond := len(p) > 0 && minimum(p) >= 0 && maximum(p) <= 1
    if !cond {
        return fmt.Errorf("p-values must be in range 0.0 to 1.0")
    }
    return nil
}
 
func ratchet(p pvalues, dir direction) {
    size := len(p)
    m := p[0]
    if dir == up {
        for i := 1; i < size; i++ {
            if p[i] > m {
                p[i] = m
            }
            m = p[i]
        }
    } else {
        for i := 1; i < size; i++ {
            if p[i] < m {
                p[i] = m
            }
            m = p[i]
        }
    }
    for i := 0; i < size; i++ {
        if p[i] > 1.0 {
            p[i] = 1.0
        }
    }
}
 
func schwartzian(p pvalues, mult pvalues, dir direction) pvalues {
    size := len(p)
    order := make([]int, size)
    iv1s := make([]iv1, size)
    for i := 0; i < size; i++ {
        iv1s[i] = iv1{i, p[i]}
    }
    if dir == up {
        sort.Slice(iv1s, func(i, j int) bool {
            return iv1s[i].value > iv1s[j].value
        })
    } else {
        sort.Slice(iv1s, func(i, j int) bool {
            return iv1s[i].value < iv1s[j].value
        })
    }
    for i := 0; i < size; i++ {
        order[i] = iv1s[i].index
    }
    pa := make(pvalues, size)
    for i := 0; i < size; i++ {
        pa[i] = mult[i] * p[order[i]]
    }
    ratchet(pa, dir)
    order2 := make([]int, size)
    iv2s := make([]iv2, size)
    for i := 0; i < size; i++ {
        iv2s[i] = iv2{i, order[i]}
    }
    sort.Slice(iv2s, func(i, j int) bool {
        return iv2s[i].value < iv2s[j].value
    })
    for i := 0; i < size; i++ {
        order2[i] = iv2s[i].index
    }
    pa2 := make(pvalues, size)
    for i := 0; i < size; i++ {
        pa2[i] = pa[order2[i]]
    }
    return pa2
}
 
func adjust(p pvalues, ctype string) pvalues {
    size := len(p)
    if size == 0 {
        return p
    }
    fsize := float64(size)
    switch ctype {
    case "Benjamini-Hochberg":
        mult := make(pvalues, size)
        for i := 0; i < size; i++ {
            mult[i] = fsize / float64(size-i)
        }
        return schwartzian(p, mult, up)
    case "Benjamini-Yekutieli":
        q := 0.0
        for i := 1; i <= size; i++ {
            q += 1.0 / float64(i)
        }
        mult := make(pvalues, size)
        for i := 0; i < size; i++ {
            mult[i] = q * fsize / (fsize - float64(i))
        }
        return schwartzian(p, mult, up)
    case "Bonferroni":
        p2 := make(pvalues, size)
        for i := 0; i < size; i++ {
            p2[i] = math.Min(p[i]*fsize, 1.0)
        }
        return p2
    case "Hochberg":
        mult := make(pvalues, size)
        for i := 0; i < size; i++ {
            mult[i] = float64(i) + 1
        }
        return schwartzian(p, mult, up)
    case "Holm":
        mult := make(pvalues, size)
        for i := 0; i < size; i++ {
            mult[i] = fsize - float64(i)
        }
        return schwartzian(p, mult, down)
    case "Hommel":
        order := make([]int, size)
        iv1s := make([]iv1, size)
        for i := 0; i < size; i++ {
            iv1s[i] = iv1{i, p[i]}
        }
        sort.Slice(iv1s, func(i, j int) bool {
            return iv1s[i].value < iv1s[j].value
        })
        for i := 0; i < size; i++ {
            order[i] = iv1s[i].index
        }
        s := make(pvalues, size)
        for i := 0; i < size; i++ {
            s[i] = p[order[i]]
        }
        m := make(pvalues, size)
        for i := 0; i < size; i++ {
            m[i] = s[i] * fsize / (float64(i) + 1)
        }
        min := minimum(m)
        q := make(pvalues, size)
        for i := 0; i < size; i++ {
            q[i] = min
        }
        pa := make(pvalues, size)
        for i := 0; i < size; i++ {
            pa[i] = min
        }
        for j := size - 1; j >= 2; j-- {
            lower := make([]int, size-j+1) // lower indices
            for i := 0; i < len(lower); i++ {
                lower[i] = i
            }
            upper := make([]int, j-1) // upper indices
            for i := 0; i < len(upper); i++ {
                upper[i] = size - j + 1 + i
            }
            qmin := float64(j) * s[upper[0]] / 2.0
            for i := 1; i < len(upper); i++ {
                temp := s[upper[i]] * float64(j) / (2.0 + float64(i))
                if temp < qmin {
                    qmin = temp
                }
            }
            for i := 0; i < len(lower); i++ {
                q[lower[i]] = math.Min(s[lower[i]]*float64(j), qmin)
            }
            for i := 0; i < len(upper); i++ {
                q[upper[i]] = q[size-j]
            }
            for i := 0; i < size; i++ {
                if pa[i] < q[i] {
                    pa[i] = q[i]
                }
            }
        }
        order2 := make([]int, size)
        iv2s := make([]iv2, size)
        for i := 0; i < size; i++ {
            iv2s[i] = iv2{i, order[i]}
        }
        sort.Slice(iv2s, func(i, j int) bool {
            return iv2s[i].value < iv2s[j].value
        })
        for i := 0; i < size; i++ {
            order2[i] = iv2s[i].index
        }
        pa2 := make(pvalues, size)
        for i := 0; i < size; i++ {
            pa2[i] = pa[order2[i]]
        }
        return pa2
    case "Šidák":
        p2 := make(pvalues, size)
        for i := 0; i < size; i++ {
            p2[i] = 1.0 - math.Pow(1.0-float64(p[i]), fsize)
        }
        return p2
    default:
        fmt.Printf("\nSorry, do not know how to do '%s' correction.\n", ctype)
        fmt.Println("Perhaps you want one of these?:")
        temp := make([]string, len(ctypes)-1)
        for i := 0; i < len(temp); i++ {
            temp[i] = fmt.Sprintf("  %s", ctypes[i])
        }
        fmt.Println(strings.Join(temp, "\n"))
        os.Exit(1)
    }
    return p
}
 
func main() {
    p := 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,
    }
    for _, ctype := range ctypes {
        s, err := adjusted(p, ctype)
        if err != nil {
            log.Fatal(err)
        }
        fmt.Println(s)
    }
}