pingouin
pingouin copied to clipboard
raphaelvallat
divide by zero encountered to calculate p-value when r==1
In correlation.py the function _correl_pvalue(r, nx, k=0) does not handle the case when r=1.0.
This results in:
lib/python3.12/site-packages/pingouin/effsize.py:152): RuntimeWarning: divide by zero encountered in arctanh
z = np.arctanh(stat) # R-to-z transform
lib/python3.12/site-packages/pingouin/correlation.py:68): RuntimeWarning: divide by zero encountered in scalar divide
tval = r * np.sqrt(dof [/](https://file+.vscode-resource.vscode-cdn.net/) (1 - r**2))
Try it with a dataset like:
0 X Y 3
0 0.281 -0.281 0.281 0.281
1 -0.171 0.171 -0.171 -0.171
2 -0.517 0.517 -0.517 -0.517
3 0.667 -0.667 0.667 0.667
4 0.835 -0.835 0.835 0.835
5 0.887 -0.887 0.887 0.887
6 -0.851 0.851 -0.851 -0.851
7 -0.315 0.315 -0.315 -0.315
8 -0.343 0.343 -0.343 -0.343
9 -0.017 0.017 -0.017 -0.017
10 -0.543 0.543 -0.543 -0.543
11 -0.749 0.749 -0.749 -0.749
12 0.643 -0.643 0.643 0.643
13 -0.371 0.371 -0.371 -0.371
14 -0.443 0.443 -0.443 -0.443
15 0.949 -0.949 0.949 0.949
16 -0.027 0.027 -0.027 -0.027
17 -0.489 0.489 -0.489 -0.489
18 0.845 -0.845 0.845 0.845
19 0.999 -0.999 0.999 0.999
This doesn't seem critical, since the user does get the warning. But addressing it would be straightforward, and should probably be handled at the same time as some other correlation edge-cases (e.g., when r ~0 in #435).
Notably, the zero division only occurs with the percentage bend correlation method (method = 'percbend').
import pandas as pd
import pingouin as pg
df = pd.DataFrame(
columns=["0", "X", "Y", "3"],
data=[
[ 0.281, -0.281, 0.281, 0.281],
[-0.171, 0.171, -0.171, -0.171],
[-0.517, 0.517, -0.517, -0.517],
[ 0.667, -0.667, 0.667, 0.667],
[ 0.835, -0.835, 0.835, 0.835],
[ 0.887, -0.887, 0.887, 0.887],
[-0.851, 0.851, -0.851, -0.851],
[-0.315, 0.315, -0.315, -0.315],
[-0.343, 0.343, -0.343, -0.343],
[-0.017, 0.017, -0.017, -0.017],
[-0.543, 0.543, -0.543, -0.543],
[-0.749, 0.749, -0.749, -0.749],
[ 0.643, -0.643, 0.643, 0.643],
[-0.371, 0.371, -0.371, -0.371],
[-0.443, 0.443, -0.443, -0.443],
[ 0.949, -0.949, 0.949, 0.949],
[-0.027, 0.027, -0.027, -0.027],
[-0.489, 0.489, -0.489, -0.489],
[ 0.845, -0.845, 0.845, 0.845],
[ 0.999, -0.999, 0.999, 0.999]
],
)
x = df["X"].to_numpy()
y = df["Y"].to_numpy()
correlation_methods = [
"pearson",
"spearman",
"kendall",
"bicor",
"percbend",
"shepherd",
"skipped",
]
print("Pingouin version:", pg.__version__)
print("=" * 40)
for method in correlation_methods:
print("Correlation method:", method)
correlation = pg.corr(x, y, method=method)
rval = correlation.at[method, "r"]
print("r value:", rval)
print("-" * 40)
Running the above code gives me:
Pingouin version: 0.5.4
========================================
Correlation method: pearson
r value: -1.0
----------------------------------------
Correlation method: spearman
r value: -1.0
----------------------------------------
Correlation method: kendall
r value: -1.0
----------------------------------------
Correlation method: bicor
r value: -0.9999999999999999
----------------------------------------
Correlation method: percbend
pingouin\correlation.py:59: RuntimeWarning: divide by zero encountered in scalar divide
tval = r * np.sqrt(dof / (1 - r**2))
r value: -1.0
----------------------------------------
Correlation method: shepherd
r value: -1.0
----------------------------------------
Correlation method: skipped
sklearn\covariance\_robust_covariance.py:748: UserWarning: The covariance matrix associated to your dataset is not full rank
warnings.warn(
pingouin\correlation.py:133: UserWarning: The skipped correlation relies on the Minimum Covariance Determinant algorithm, which gives slightly different results in Python (scikit-learn) than in the original Matlab library (LIBRA). As such, the skipped correlation may be different from the Matlab robust correlation toolbox (see issue 164 on Pingouin's GitHub). Make sure to double check your results or use another robust correlation method.
warnings.warn(
r value: -1.0
----------------------------------------
Same thing with pcorr():
[lib/python3.12/site-packages/pingouin/effsize.py:152](lib/python3.12/site-packages/pingouin/effsize.py:152): RuntimeWarning:
divide by zero encountered in arctanh
z = np.arctanh(stat) # R-to-z transform
[lib/python3.12/site-packages/pingouin/correlation.py:59](lib/python3.12/site-packages/pingouin/correlation.py:59): RuntimeWarning:
divide by zero encountered in scalar divide
And there is no percbend-ing: pg.pcorr(dataset[[variable, 'X', 'Y']])
The problem with the warning is that the output is undefined. Should the value be 0, 1 or -1?