Renamed 'r' with 'pointbiserialr' in convert_effsize (#325)

* Renamed 'r' with 'pointbiserialr' in convert_effsize

* Updated changelog + unit tests

* Temp fix for bug in plot_paired

Author: raphaelvallat

docs/changelog.rst

@@ -14,6 +14,7 @@ v0.6.0.dev
 **Bugfixes**
 
 - Fixed a bug where the boolean value returned by :py:func:`pingouin.anderson` was inverted. It returned True when the data was NOT coming from the tested distribution, and vice versa. `PR 308 <https://github.com/raphaelvallat/pingouin/pull/308>`_.
+- Fixed misleading documentation and ``input_type`` in the :py:func:`convert_effsize` function. When converting from a Cohen's d effect size to a correlation coefficient, the resulting correlation is **not** a Pearson correlation but instead a `point-biserial correlation <https://en.wikipedia.org/wiki/Point-biserial_correlation_coefficient>`_. To avoid any confusion, ``input_type='r'`` has been deprecated and replaced with ``input_type='pointbiserialr'``. For more details, see `issue 302 <https://github.com/raphaelvallat/pingouin/issues/302>`_.
 
 **New functions**
 

notebooks/03_EffectSizes.ipynb

@@ -32,13 +32,13 @@ def compute_esci(
     Parameters
     ----------
     stat : float
-        Original effect size. Must be either a correlation coefficient or a
-        Cohen-type effect size (Cohen d or Hedges g).
+        Original effect size. Must be either a correlation coefficient or a Cohen-type effect size
+        (Cohen d or Hedges g).
     nx, ny : int
         Length of vector x and y.
     paired : bool
-        Indicates if the effect size was estimated from a paired sample.
-        This is only relevant for cohen or hedges effect size.
+        Indicates if the effect size was estimated from a paired sample. This is only relevant for
+        cohen or hedges effect size.
     eftype : string
         Effect size type. Must be "r" (correlation) or "cohen" (Cohen d or Hedges g).
     confidence : float
@@ -57,9 +57,8 @@ def compute_esci(
 
     Notes
     -----
-    To compute the parametric confidence interval around a
-    **Pearson r correlation** coefficient, one must first apply a
-    Fisher's r-to-z transformation:
+    To compute the parametric confidence interval around a **Pearson r correlation** coefficient,
+    one must first apply a Fisher's r-to-z transformation:
 
     .. math:: z = 0.5 \\cdot \\ln \\frac{1 + r}{1 - r} = \\text{arctanh}(r)
 
@@ -69,14 +68,12 @@ def compute_esci(
 
     where :math:`n` is the sample size.
 
-    The lower and upper confidence intervals - *in z-space* - are then
-    given by:
+    The lower and upper confidence intervals - *in z-space* - are then given by:
 
     .. math:: \\text{ci}_z = z \\pm \\text{crit} \\cdot \\text{SE}
 
-    where :math:`\\text{crit}` is the critical value of the normal distribution
-    corresponding to the desired confidence level (e.g. 1.96 in case of a 95%
-    confidence interval).
+    where :math:`\\text{crit}` is the critical value of the normal distribution corresponding to
+    the desired confidence level (e.g. 1.96 in case of a 95% confidence interval).
 
     These confidence intervals can then be easily converted back to *r-space*:
 
@@ -85,10 +82,9 @@ def compute_esci(
         \\text{ci}_r = \\frac{\\exp(2 \\cdot \\text{ci}_z) - 1}
         {\\exp(2 \\cdot \\text{ci}_z) + 1} = \\text{tanh}(\\text{ci}_z)
 
-    A formula for calculating the confidence interval for a
-    **Cohen d effect size** is given by Hedges and Olkin (1985, p86).
-    If the effect size estimate from the sample is :math:`d`, then it follows a
-    T distribution with standard error:
+    A formula for calculating the confidence interval for a **Cohen d effect size** is given by
+    Hedges and Olkin (1985, p86). If the effect size estimate from the sample is :math:`d`, then
+    it follows a T distribution with standard error:
 
     .. math::
 
@@ -107,15 +103,14 @@ def compute_esci(
 
     .. math:: \\text{ci}_d = d \\pm \\text{crit} \\cdot \\text{SE}
 
-    where :math:`\\text{crit}` is the critical value of the T distribution
-    corresponding to the desired confidence level.
+    where :math:`\\text{crit}` is the critical value of the T distribution corresponding to the
+    desired confidence level.
 
     References
     ----------
     * https://en.wikipedia.org/wiki/Fisher_transformation
 
-    * Hedges, L., and Ingram Olkin. "Statistical models for meta-analysis."
-      (1985).
+    * Hedges, L., and Ingram Olkin. "Statistical models for meta-analysis." (1985).
 
     * http://www.leeds.ac.uk/educol/documents/00002182.htm
 
@@ -211,8 +206,7 @@ def compute_bootci(
     y : 1D-array, list, or None
         Second sample. Required only for bivariate functions.
     func : str or custom function
-        Function to compute the bootstrapped statistic.
-        Accepted string values are:
+        Function to compute the bootstrapped statistic. Accepted string values are:
 
         * ``'pearson'``: Pearson correlation (bivariate, paired x and y)
         * ``'spearman'``: Spearman correlation (bivariate, paired x and y)
@@ -501,12 +495,13 @@ def convert_effsize(ef, input_type, output_type, nx=None, ny=None):
     ef : float
         Original effect size.
     input_type : string
-        Effect size type of ef. Must be ``'r'`` or ``'cohen'``.
+        Effect size type of ef. Must be ``'cohen'`` or ``'pointbiserialr'``.
     output_type : string
         Desired effect size type. Available methods are:
 
         * ``'cohen'``: Unbiased Cohen d
         * ``'hedges'``: Hedges g
+        * ``'pointbiserialr'``: Point-biserial correlation
         * ``'eta-square'``: Eta-square
         * ``'odds-ratio'``: Odds ratio
         * ``'AUC'``: Area Under the Curve
@@ -527,15 +522,17 @@ def convert_effsize(ef, input_type, output_type, nx=None, ny=None):
 
     Notes
     -----
-    The formula to convert **r** to **d** is given in [1]_:
+    The formula to convert from a`point-biserial correlation
+    <https://en.wikipedia.org/wiki/Point-biserial_correlation_coefficient>`_ **r** to **d** is
+    given in [1]_:
 
-    .. math:: d = \\frac{2r}{\\sqrt{1 - r^2}}
+    .. math:: d = \\frac{2r_{pb}}{\\sqrt{1 - r_{pb}^2}}
 
-    The formula to convert **d** to **r** is given in [2]_:
+    The formula to convert **d** to a point-biserial correlation **r** is given in [2]_:
 
     .. math::
 
-        r = \\frac{d}{\\sqrt{d^2 + \\frac{(n_x + n_y)^2 - 2(n_x + n_y)}
+        r_{pb} = \\frac{d}{\\sqrt{d^2 + \\frac{(n_x + n_y)^2 - 2(n_x + n_y)}
         {n_xn_y}}}
 
     The formula to convert **d** to :math:`\\eta^2` is given in [3]_:
@@ -584,35 +581,35 @@ def convert_effsize(ef, input_type, output_type, nx=None, ny=None):
     >>> pg.convert_effsize(.45, 'cohen', 'hedges', nx=10, ny=10)
     0.4309859154929578
 
-    3. Convert Pearson r to Cohen d
+    3. Convert a point-biserial correlation to Cohen d
 
-    >>> r = 0.40
-    >>> d = pg.convert_effsize(r, 'r', 'cohen')
+    >>> rpb = 0.40
+    >>> d = pg.convert_effsize(rpb, 'pointbiserialr', 'cohen')
     >>> print(d)
     0.8728715609439696
 
-    4. Reverse operation: convert Cohen d to Pearson r
+    4. Reverse operation: convert Cohen d to a point-biserial correlation
 
-    >>> pg.convert_effsize(d, 'cohen', 'r')
+    >>> pg.convert_effsize(d, 'cohen', 'pointbiserialr')
     0.4000000000000001
     """
     it = input_type.lower()
     ot = output_type.lower()
 
     # Check input and output type
-    for input in [it, ot]:
-        if not _check_eftype(input):
-            err = "Could not interpret input '{}'".format(input)
+    for inp in [it, ot]:
+        if not _check_eftype(inp):
+            err = "Could not interpret input '{}'".format(inp)
             raise ValueError(err)
-    if it not in ["r", "cohen"]:
-        raise ValueError("Input type must be 'r' or 'cohen'")
+    if it not in ["pointbiserialr", "cohen"]:
+        raise ValueError("Input type must be 'cohen' or 'pointbiserialr'")
 
     # Pass-through option
     if it == ot or ot == "none":
         return ef
 
-    # Convert r to Cohen d (Rosenthal 1994)
-    d = (2 * ef) / np.sqrt(1 - ef**2) if it == "r" else ef
+    # Convert point-biserial r to Cohen d (Rosenthal 1994)
+    d = (2 * ef) / np.sqrt(1 - ef**2) if it == "pointbiserialr" else ef
 
     # Then convert to the desired output type
     if ot == "cohen":
@@ -627,7 +624,7 @@ def convert_effsize(ef, input_type, output_type, nx=None, ny=None):
                 "Hedges g. Returning Cohen's d instead"
             )
             return d
-    elif ot == "r":
+    elif ot == "pointbiserialr":
         # McGrath and Meyer 2006
         if all(v is not None for v in [nx, ny]):
             a = ((nx + ny) ** 2 - 2 * (nx + ny)) / (nx * ny)
@@ -640,6 +637,12 @@ def convert_effsize(ef, input_type, output_type, nx=None, ny=None):
     elif ot == "odds-ratio":
         # Borenstein et al. 2009
         return np.exp(d * np.pi / np.sqrt(3))
+    elif ot == "r":
+        # https://github.com/raphaelvallat/pingouin/issues/302
+        raise ValueError(
+            "Using effect size 'r' in `pingouin.convert_effsize` has been deprecated. "
+            "Please use 'pointbiserialr' instead."
+        )
     else:  # ['auc']
         # Ruscio 2008
         from scipy.stats import norm
@@ -666,7 +669,8 @@ def compute_effsize(x, y, paired=False, eftype="cohen"):
         * ``'none'``: no effect size
         * ``'cohen'``: Unbiased Cohen d
         * ``'hedges'``: Hedges g
-        * ``'r'``: correlation coefficient
+        * ``'r'``: Pearson correlation coefficient
+        * ``'pointbiserialr'``: Point-biserial correlation
         * ``'eta-square'``: Eta-square
         * ``'odds-ratio'``: Odds ratio
         * ``'AUC'``: Area Under the Curve
@@ -684,8 +688,8 @@ def compute_effsize(x, y, paired=False, eftype="cohen"):
 
     Notes
     -----
-    Missing values are automatically removed from the data. If ``x`` and ``y``
-    are paired, the entire row is removed.
+    Missing values are automatically removed from the data. If ``x`` and ``y`` are paired, the
+    entire row is removed.
 
     If ``x`` and ``y`` are independent, the Cohen :math:`d` is:
 
@@ -702,22 +706,19 @@ def compute_effsize(x, y, paired=False, eftype="cohen"):
         d_{avg} = \\frac{\\overline{X} - \\overline{Y}}
         {\\sqrt{\\frac{(\\sigma_1^2 + \\sigma_2^2)}{2}}}
 
-    The Cohen’s d is a biased estimate of the population effect size,
-    especially for small samples (n < 20). It is often preferable
-    to use the corrected Hedges :math:`g` instead:
+    The Cohen's d is a biased estimate of the population effect size, especially for small samples
+    (n < 20). It is often preferable to use the corrected Hedges :math:`g` instead:
 
     .. math:: g = d \\times (1 - \\frac{3}{4(n_1 + n_2) - 9})
 
-    The common language effect size is the proportion of pairs where ``x`` is
-    higher than ``y`` (calculated with a brute-force approach where
-    each observation of ``x`` is paired to each observation of ``y``,
-    see :py:func:`pingouin.wilcoxon` for more details):
+    The common language effect size is the proportion of pairs where ``x`` is higher than ``y``
+    (calculated with a brute-force approach where each observation of ``x`` is paired to each
+    observation of ``y``, see :py:func:`pingouin.wilcoxon` for more details):
 
     .. math:: \\text{CL} = P(X > Y) + .5 \\times P(X = Y)
 
-    For other effect sizes, Pingouin will first calculate a Cohen :math:`d` and
-    then use the :py:func:`pingouin.convert_effsize` to convert to the desired
-    effect size.
+    For other effect sizes, Pingouin will first calculate a Cohen :math:`d` and then use the
+    :py:func:`pingouin.convert_effsize` to convert to the desired effect size.
 
     References
     ----------
@@ -822,7 +823,7 @@ def compute_effsize_from_t(tval, nx=None, ny=None, N=None, eftype="cohen"):
     N : int, optional
         Total sample size (will not be used if nx and ny are specified)
     eftype : string, optional
-        desired output effect size
+        Desired output effect size.
 
     Returns
     -------

pingouin/effsize.py

@@ -1,8 +1,9 @@
-import pandas as pd
-import numpy as np
 import pytest
-
+import numpy as np
+import pandas as pd
 from unittest import TestCase
+from scipy.stats import pearsonr, pointbiserialr
+
 from pingouin.effsize import compute_esci, compute_effsize, compute_effsize_from_t, compute_bootci
 from pingouin.effsize import convert_effsize as cef
 
@@ -216,29 +217,48 @@ def test_compute_boot_esci(self):
 
     def test_convert_effsize(self):
         """Test function convert_effsize.
+
         Compare to https://www.psychometrica.de/effect_size.html
         """
         # Cohen d
         d = 0.40
         assert cef(d, "cohen", "none") == d
-        assert round(cef(d, "cohen", "r"), 4) == 0.1961
-        cef(d, "cohen", "r", nx=10, ny=12)  # When nx and ny are specified
-        assert np.allclose(cef(1.002549, "cohen", "r"), 0.4481248)  # R
+        assert round(cef(d, "cohen", "pointbiserialr"), 4) == 0.1961
+        cef(d, "cohen", "pointbiserialr", nx=10, ny=12)  # When nx and ny are specified
+        assert np.allclose(cef(1.002549, "cohen", "pointbiserialr"), 0.4481248)  # R
         assert round(cef(d, "cohen", "eta-square"), 4) == 0.0385
         assert round(cef(d, "cohen", "odds-ratio"), 4) == 2.0658
         cef(d, "cohen", "hedges", nx=10, ny=10)
-        cef(d, "cohen", "r")
+        cef(d, "cohen", "pointbiserialr")
         cef(d, "cohen", "hedges")
 
-        # Correlation coefficient
-        r = 0.65
-        assert cef(r, "r", "none") == r
-        assert round(cef(r, "r", "cohen"), 4) == 1.7107
-        assert np.allclose(cef(0.4481248, "r", "cohen"), 1.002549)
-        assert round(cef(r, "r", "eta-square"), 4) == 0.4225
-        assert round(cef(r, "r", "odds-ratio"), 4) == 22.2606
+        # Point-biserial correlation
+        rpb = 0.65
+        assert cef(rpb, "pointbiserialr", "none") == rpb
+        assert round(cef(rpb, "pointbiserialr", "cohen"), 4) == 1.7107
+        assert np.allclose(cef(0.4481248, "pointbiserialr", "cohen"), 1.002549)
+        assert round(cef(rpb, "pointbiserialr", "eta-square"), 4) == 0.4225
+        assert round(cef(rpb, "pointbiserialr", "odds-ratio"), 4) == 22.2606
+        # Using actual values
+        np.random.seed(42)
+        x1, y1 = np.random.multivariate_normal(mean=[1, 2], cov=[[1, 0.5], [0.5, 1]], size=100).T
+        xy1 = np.hstack((x1, y1))
+        xy1_bool = np.repeat([0, 1], 100)
+        # Let's calculate the ground-truth point-biserial correlation
+        r_biserial = pearsonr(xy1_bool, xy1)[0]  # 0.50247
+        assert np.isclose(r_biserial, pointbiserialr(xy1_bool, xy1)[0])
+        # Now the Cohen's d
+        d = abs(compute_effsize(x1, y1, paired=True, eftype="cohen"))  # 1.15651
+        # And now we can convert point-biserial r <--> d
+        r_convert = cef(abs(d), "cohen", "pointbiserialr", nx=100, ny=100)  # 0.50247
+        assert np.isclose(r_convert, r_biserial)
+        d_convert = cef(r_biserial, "pointbiserialr", "cohen", nx=100, ny=100)  # 1.162
+        assert abs(d - d_convert) < 0.1
 
         # Error
+        with pytest.raises(ValueError):
+            # DEPRECATED - https://github.com/raphaelvallat/pingouin/issues/302
+            cef(d, "cohen", "r")
         with pytest.raises(ValueError):
             cef(d, "coucou", "hibou")
         with pytest.raises(ValueError):
@@ -252,6 +272,7 @@ def test_compute_effsize(self):
         compute_effsize(x=x, y=y, eftype="odds-ratio", paired=False)
         compute_effsize(x=x, y=y, eftype="eta-square", paired=False)
         compute_effsize(x=x, y=y, eftype="cles", paired=False)
+        compute_effsize(x=x, y=y, eftype="pointbiserialr", paired=False)
         compute_effsize(x=x, y=y, eftype="none", paired=False)
         # Unequal variances
         z = np.random.normal(2.5, 3, 30)

pingouin/tests/test_effsize.py

@@ -325,6 +325,7 @@ def _check_eftype(eftype):
         "hedges",
         "cohen",
         "r",
+        "pointbiserialr",
         "eta-square",
         "odds-ratio",
         "auc",

pingouin/utils.py

@@ -3,3 +3,4 @@ codecov
 pytest-cov
 openpyxl
 mpmath
+numpy<=1.23