Add FDR and Bonferroni correction functions (#87)

* Add FDR and Bonferroni functions.

* Add references.

* Fix directive.

Author: tsalo

docs/api.rst

@@ -113,3 +113,5 @@ API
    stats.ensure_2d
    stats.q_profile
    stats.q_gen
+   stats.bonferroni
+   stats.fdr

docs/conf.py

@@ -168,6 +168,7 @@
     "sklearn": ("https://scikit-learn.org/stable", (None, "./_intersphinx/sklearn-objects.inv")),
     "matplotlib": ("https://matplotlib.org/", (None, "https://matplotlib.org/objects.inv")),
     "pandas": ("https://pandas.pydata.org/pandas-docs/stable/", None),
+    "statsmodels": ("http://www.statsmodels.org/stable/", None),
 }
 
 sphinx_gallery_conf = {

docs/references.bib

@@ -1,3 +1,38 @@
+@article{benjamini1995controlling,
+  title={Controlling the false discovery rate: a practical and powerful approach to multiple testing},
+  author={Benjamini, Yoav and Hochberg, Yosef},
+  journal={Journal of the Royal statistical society: series B (Methodological)},
+  volume={57},
+  number={1},
+  pages={289--300},
+  year={1995},
+  publisher={Wiley Online Library},
+  url={https://doi.org/10.1111/j.2517-6161.1995.tb02031.x},
+  doi={10.1111/j.2517-6161.1995.tb02031.x}
+}
+
+@article{benjamini2001control,
+  author={Yoav Benjamini and Daniel Yekutieli},
+  title={The control of the false discovery rate in multiple testing under dependency},
+  volume={29},
+  journal={The Annals of Statistics},
+  number={4},
+  publisher={Institute of Mathematical Statistics},
+  pages={1165 -- 1188},
+  year={2001},
+  doi={10.1214/aos/1013699998},
+  url={https://doi.org/10.1214/aos/1013699998}
+}
+
+@article{bonferroni1936teoria,
+  title={Teoria statistica delle classi e calcolo delle probabilita},
+  author={Bonferroni, Carlo},
+  journal={Pubblicazioni del R Istituto Superiore di Scienze Economiche e Commericiali di Firenze},
+  volume={8},
+  pages={3--62},
+  year={1936}
+}
+
 @article{brockwell2001comparison,
   title={A comparison of statistical methods for meta-analysis},
   author={Brockwell, Sarah E and Gordon, Ian R},
@@ -88,6 +123,17 @@ @article{sangnawakij2019meta
   doi={10.1177/0962280217718867}
 }
 
+@article{shaffer1995multiple,
+  title={Multiple hypothesis testing},
+  author={Shaffer, Juliet Popper},
+  journal={Annual review of psychology},
+  volume={46},
+  number={1},
+  pages={561--584},
+  year={1995},
+  publisher={Annual Reviews 4139 El Camino Way, PO Box 10139, Palo Alto, CA 94303-0139, USA}
+}
+
 @article{stouffer1949american,
   title={The american soldier: Adjustment during army life.(studies in social psychology in world war ii), vol. 1},
   author={Stouffer, Samuel A and Suchman, Edward A and DeVinney, Leland C and Star, Shirley A and Williams Jr, Robin M},

pymare/stats.py

@@ -146,3 +146,90 @@ def q_gen(y, v, X, tau2):
     beta = weighted_least_squares(y, v, X, tau2)
     w = 1.0 / (v + tau2)
     return (w * (y - X.dot(beta)) ** 2).sum(0)
+
+
+def bonferroni(p_values):
+    """Perform Bonferroni correction on p values.
+
+    This correction is based on the one described in :footcite:t:`bonferroni1936teoria` and
+    :footcite:t:`shaffer1995multiple`.
+
+    .. versionadded:: 0.0.4
+
+    Parameters
+    ----------
+    p_values : :obj:`numpy.ndarray`
+        Uncorrected p values.
+
+    Returns
+    -------
+    p_corr : :obj:`numpy.ndarray`
+        Corrected p values.
+
+    References
+    ----------
+    .. footbibliography::
+    """
+    p_corr = p_values * p_values.size
+    p_corr[p_corr > 1] = 1
+    return p_corr
+
+
+def fdr(p_values, q=0.05, method="bh"):
+    """Perform FDR correction on p values.
+
+    .. versionadded:: 0.0.4
+
+    Parameters
+    ----------
+    p_values : :obj:`numpy.ndarray`
+        Array of p values.
+    q : :obj:`float`, optional
+        Alpha value. Default is 0.05.
+    method : {"bh", "by"}, optional
+        Method to use for correction.
+        Either "bh" (Benjamini-Hochberg :footcite:p:`benjamini1995controlling`) or
+        "by" (Benjamini-Yekutieli :footcite:p:`benjamini2001control`).
+        Default is "bh".
+
+    Returns
+    -------
+    p_adjusted : :obj:`numpy.ndarray`
+        Array of adjusted p values.
+
+    Notes
+    -----
+    This function is adapted from ``statsmodels``, which is licensed under a BSD-3 license.
+
+    References
+    ----------
+    .. footbibliography::
+
+    See Also
+    --------
+    statsmodels.stats.multitest.fdrcorrection
+    """
+    sort_idx = np.argsort(p_values)
+    revert_idx = np.argsort(sort_idx)
+    p_sorted = p_values[sort_idx]
+
+    n_tests = p_values.size
+
+    # empirical cumulative density function
+    ecdf = np.linspace(0, 1, n_tests + 1)[1:]
+    if method == "by":
+        # NOTE: I don't know what cm stands for
+        cm = np.sum(1 / np.arange(1, n_tests + 1))
+        ecdffactor = ecdf / cm
+    else:
+        ecdffactor = ecdf
+
+    p_adjusted = p_sorted / ecdffactor
+    p_adjusted = np.minimum.accumulate(p_adjusted[::-1])[::-1]
+    # NOTE: Why not this?
+    # p_adjusted = np.maximum.accumulate(p_adjusted)
+
+    p_adjusted[p_adjusted > 1] = 1
+    p_adjusted = p_adjusted[revert_idx]
+
+    return p_adjusted