scipy.stats: 点双列相関係数 pointbiserialr

特別な名前がついているが，ピアソンの積率相関係数において，片方の変数が 2 値データの場合ということである。

pointbiserialr は point biserial correlation coefficient r で，訳すと，点双列相関係数ということである。

2 値変数は連続変数なので（知らない人も多いかもしれないが），当たり前なのだが，その昔，計算環境が劣悪だった頃は，特別な場合に簡単な計算式で計算できるということで誕生したものかな。

普通に（数値に変換しなければならないかもしれないが） scipy.stats.pearsonr を使えばよい。

import numpy as np

a = np.array([0, 0, 0, 1, 1, 1, 1])
b = np.arange(7)

1． 点双列相関係数

pointbiserialr(x, y)

戻り値はピアソンの積率相関係数と無相関検定（$H_0:\rho=0$）の両側検定の $p$ 値である。

from scipy.stats import pointbiserialr

pointbiserialr(a, b)

PointbiserialrResult(correlation=0.8660254037844387, pvalue=0.011724811003954626)

ちなみに，相関係数が $r$，サンプルサイズが $n$ のとき無相関検定の両側検定の $p$ 値は以下のように求められる。

from scipy.stats import t

def pvalue(r, n):
    return t.sf(abs(r) * np.sqrt(n - 2) / np.sqrt(1 - r**2), n-2) * 2 if r**2 != 0 else 1

pvalue(0.8660254037844387, len(a))

0.011724811003954614

2.　ピアソンの積率相関係数

pearsonr(x, y)

戻り値はピアソンの積率相関係数と無相関検定（$H_0: \rho = 0$）の両側検定の $p$ 値である。

from scipy.stats import pearsonr

pearsonr(a, b)

(0.8660254037844387, 0.011724811003954626)

3.　numpy の corrcoef

numpy には corrcoef がある。

corrcoef(x, y=None, rowvar=True, bias=<no value>, ddof=<no value>, *, dtype=None)

np.corrcoef(a, b)[0,1]

0.8660254037844385

ちなみに，np.corrcoef() は相関係数行列を求めることができる。無相関検定の結果は得られない。

x = [[2,1,2,3,4], [3,2,1,2,3], [5,4,3,4,5], [2,1,3,4,5]]
r = np.corrcoef(x)
print(r)

[[1.         0.41931393 0.41931393 0.97072534]
 [0.41931393 1.         1.         0.18898224]
 [0.41931393 1.         1.         0.18898224]
 [0.97072534 0.18898224 0.18898224 1.        ]]

以下のような関数で $p$ 値の行列を求めることができる（前出の pvalue() を使用する）。

def r(x):
    k, n = np.array(x).shape
    p = np.ones((k, k))
    r = np.corrcoef(x)
    for i in range(k):
        for j in range(i+1, k):
            p[i, j] = p[j, i] = pvalue(r[i, j], n)
    return p

r(x)

array([[1.00000000e+00, 4.82198952e-01, 4.82198952e-01, 5.98625570e-03],
       [4.82198952e-01, 1.00000000e+00, 1.40426542e-24, 7.60820376e-01],
       [4.82198952e-01, 1.40426542e-24, 1.00000000e+00, 7.60820376e-01],
       [5.98625570e-03, 7.60820376e-01, 7.60820376e-01, 1.00000000e+00]])

4.　φ 係数

ファイ係数は，2変数がともに 2 値データのときのピアソンの積率相関係数の絶対値に該当する。

さらに，2×2 分割表の場合には $\phi$ 係数は scipy.stats.contingency のデフォルトで得られる クラメールの $V$ に一致する。

import numpy as np
np.random.seed(1234)
from scipy.stats import uniform
x = uniform.rvs(size=100) > 0.5
y = uniform.rvs(size=100) > 0.5
import pandas as pd
df = pd.DataFrame({'x': x, 'y': y})
tbl = pd.crosstab(df['x'], df['y'])
print(tbl)

y      False  True
x                 
False     17    26
True      27    30

from scipy.stats.contingency import association

association(tbl)

0.07812845233116676

from scipy.stats import pearsonr

pearsonr(x, y)

(-0.07812845233116661, 0.43973315252059714)

なお，$\phi$ 係数に適切な符号を付けてピアソンの積率相関係数と同じにするには，分割表の 4 つのマス目の数値を a, b, c, d としたとき，$a\cdot d - b \cdot c$ の符号をつければよい。

例の場合だと 17*30 - 26*27 = -192 なので，0.07812845233116676 の符号は - である。

同じ分割表 tbl に基づく検定は，相関係数の検定の他に独立性の検定（いわゆる $\chi^2$ 検定）もある。同じ結果になるだろうか？

独立性の検定は scipy.stats.chi2_contingency で行える。
連続性の補正は correction で指定できる（デフォルトで　　True）。
戻り値は，$\chi^2$ 値，$p$ 値，自由度，期待値行列である。

from scipy.stats import chi2_contingency
chi2_contingency(tbl)

(0.3338817445304601,
 0.5633822668015476,
 1,
 array([[18.92, 24.08],
        [25.08, 31.92]]))

chi2_contingency(tbl, correction=False)

(0.6104055063663398,
 0.43463517068036106,
 1,
 array([[18.92, 24.08],
        [25.08, 31.92]]))

同じ $p$ 値になるわけではないということか。