Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

point-biserial / biserial

点双列相関係数（point-biserial correlation）

点双列相関係数 （point-biserial correlation） は 連続変数と自然な二値変数（連続変数を離散化したわけではない二値変数）の間の相関を測るための相関係数。

参考文献
import numpy as np

def point_biserial_correlation(x: np.ndarray, y: np.ndarray) -> float:
    """
    Compute the point-biserial correlation between a continuous variable x
    and a binary variable y (0 or 1), assuming y is a true categorical variable.

    Parameters
    ----------
    x : np.ndarray
        Continuous variable.
    y : np.ndarray
        Binary variable (0 and 1), true categories.

    Returns
    -------
    float
        Point-biserial correlation coefficient.
    """
    x = np.asarray(x)
    y = np.asarray(y)
    
    assert set(np.unique(y)).issubset({0, 1}), "y must be binary (0/1)"
    
    x1 = x[y == 1]
    x0 = x[y == 0]
    
    M1 = np.mean(x1)
    M0 = np.mean(x0)
    s = np.std(x, ddof=0)
    
    p = np.mean(y)
    q = 1 - p

    return (M1 - M0) / s * np.sqrt(p * q)

np.random.seed(0)
x = np.random.normal(size=100)
x_ = x + abs(x.min())
p = (x_ - x_.min()) / x_.max()
y = np.random.binomial(n=1, p=p, size=100)

point_biserial_correlation(x, y)
np.float64(0.42540375845000344)

scipy.stats にも実装がある

point-biserial と Pearsonの積率相関係数は等しい

Y∈{0,1}Y\in\{0,1\}のときのピアソンの積率相関係数はpoint-biserialと等しい

Proof of Point-Biserial Correlation being a special case of Pearson Correlation - Cross Validated

YYが二値変数のため、回帰直線を描くとY=0Y=0の点のXXの平均M0M_0とY=1Y=1の点のXXの平均M1M_1の2点の直線になる。

この回帰直線の傾きはβ=M1−M0/(1−0)=M1−M0\beta = M_1 - M_0 / (1 - 0) = M_1 - M_0

ピアソンの相関係数の定義は

r=Cov⁡(X,Y)sXsYr=\frac{\operatorname{Cov}(X, Y)}{s_X s_Y}
(2)

であり、回帰係数の定義から

β=Cov⁡(X,Y)Var⁡(Y)=sXsYr  ⟹  r=β⋅sYsX\beta=\frac{\operatorname{Cov}(X, Y)}{\operatorname{Var}(Y)}=\frac{s_X}{s_Y} r \implies r = \beta \cdot \frac{s_Y}{s_X}
(3)

であるため

r=M1−M0⋅sYsX=M1−M0sX⋅p(1−p)r = M_1 - M_0 \cdot \frac{s_Y}{s_X} =\frac{M_1 - M_0}{s_X} \cdot \sqrt{p(1-p)}
(4)

これはpoint-biserialに等しい

from scipy.stats import pointbiserialr, pearsonr
pointbiserialr(x, y) == pearsonr(x, y)
True

双列相関係数（biserial correlation）

バイシリアル相関係数 （biserial correlation, 双列相関係数とも） は、連続変数と人工的に二値化した変数（連続変数を閾値で分けたもの）の間の相関係数。

仮定：

  • YY が自然なカテゴリ（二値）ではなく、連続変数を人工的にしきい値で切ったものという仮定が必要

  • 連続変数 XX のほうは正規分布に近いことが望ましい

import numpy as np
from scipy.stats import norm

def biserial_correlation(x: np.ndarray, y: np.ndarray) -> float:
    """
    Compute the biserial correlation coefficient between a continuous variable x
    and a dichotomized variable y (0 or 1), assuming y was split from a latent normal variable.

    Parameters
    ----------
    x : np.ndarray
        Continuous variable.
    y : np.ndarray
        Dichotomous variable (0 and 1), assumed to be derived from a latent normal variable.

    Returns
    -------
    float
        Biserial correlation coefficient.
    """
    x = np.asarray(x)
    y = np.asarray(y)

    assert set(np.unique(y)).issubset({0, 1}), "y must be binary (0/1)"
    
    x1 = x[y == 1]
    x0 = x[y == 0]
    M1 = np.mean(x1)
    M0 = np.mean(x0)
    s = np.std(x, ddof=1)

    p = np.mean(y)
    q = 1 - p
    z = norm.ppf(p)
    phi = norm.pdf(z)

    return (M1 - M0) / s * (p * q) / phi
from scipy.stats import multivariate_normal

rho = 0.05
cov = np.array([[1, rho], [rho, 1]])
X = multivariate_normal.rvs(cov=cov, size=100, random_state=0)
x = X[:, 0]
y = 1 * (X[:, 1] > 0.5)
biserial_correlation(x,  y)
np.float64(0.06162599756903869)

背景

Biserial correlationは Pearson (1909) が「変数 AA の階級ごとに、潜在変数 BB が閾値を超える割合だけがある」という状況での相関の議論を開始したあたりから研究が始まっている

biserialとpoint-biserialのイメージの違い

  • biserialは人工的な二値変数が対象なので1つの連続値ylatenty_{\text{latent}}の分布をある閾値で切断したものを扱っている

  • point-biserialは自然な二値変数なので2つのクラス{0,1}\{0, 1\}の分布はそれぞれ分かれており、重なることもありうるイメージ。

    • 例えば潜在的な能力ylatenty_{\text{latent}}が高い人が正答する（y=1y=1になる）確率は高いが、100％ではなく偶然誤答することもありうる

Source
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
np.random.seed(42)
n = 500
x = np.random.normal(loc=0, scale=1, size=n)
threshold = 0.3
y_biserial = (x > threshold).astype(int)
y_point = (x + np.random.normal(0, 0.5, n) > 0).astype(int)

fig, axs = plt.subplots(1, 2, figsize=(8, 3), sharey=True)

axs[0].hist(x[y_biserial == 0], bins=20, alpha=0.6, label="y=0")
axs[0].hist(x[y_biserial == 1], bins=20, alpha=0.6, label="y=1")
axs[0].axvline(threshold, color='k', linestyle='--', label='Threshold')
axs[0].set_title("Biserial: Artificial Dichotomization")
axs[0].set_xlabel(r"$y_{\text{latent}}$")
axs[0].legend()

axs[1].hist(x[y_point == 0], bins=20, alpha=0.6, label="y=0")
axs[1].hist(x[y_point == 1], bins=20, alpha=0.6, label="y=1")
axs[1].set_title("Point-Biserial: Natural Categories")
axs[1].set_xlabel(r"$y_{\text{latent}}$")
axs[1].legend()

plt.suptitle("Comparison of Biserial vs Point-Biserial Correlation")
plt.tight_layout()
plt.show()
<Figure size 800x300 with 2 Axes>
Source
# Create scatter plots with color-coded categories for better intuition
fig, axs = plt.subplots(1, 2, figsize=(8, 3), sharey=True)

# Biserial: artificial thresholding
axs[0].scatter(x, y_biserial + np.random.normal(0, 0.02, size=n), alpha=0.4, c=y_biserial, cmap='coolwarm')
axs[0].axvline(threshold, color='k', linestyle='--', label='Threshold')
axs[0].set_title("Biserial: Artificial Dichotomy")
axs[0].set_xlabel("X (continuous)")
axs[0].set_ylabel("Y (binary)")
axs[0].legend()

# Point-biserial: natural binary
axs[1].scatter(x, y_point + np.random.normal(0, 0.02, size=n), alpha=0.4, c=y_point, cmap='coolwarm')
axs[1].set_title("Point-Biserial: Natural Categories")
axs[1].set_xlabel("X (continuous)")

plt.suptitle("Scatter Plot: Biserial vs Point-Biserial Correlation")
plt.tight_layout()
plt.show()
<Figure size 800x300 with 2 Axes>

Point-Biserialに補正係数をかけたのがBiserial

point-biserial （=Pearson’s r）

rpbi=Xˉ1−Xˉ0sX⋅pqr_{\text{pbi}}=\frac{\bar{X}_1-\bar{X}_0}{s_X} \cdot \sqrt{p q}
(7)

に対して補正係数pqϕ(z)\frac{\sqrt{pq}}{\phi(z)}をかけたのがBiserial

rbi=rpbi⋅pqϕ(z)=Xˉ1−Xˉ0sX⋅pqϕ(z)r_{\text{bi}} = r_{\text{pbi}} \cdot \frac{\sqrt{pq}}{\phi(z)} = \frac{\bar{X}_1-\bar{X}_0}{s_X} \cdot \frac{p q}{\phi(z)}
(8)

Peters & Van Voorhis (1940) は相関係数がρ\rhoの 2 変量正規分布に従う確率変数 X,YX, Yがあるとき、XXを平均値や中央値で二値化した確率変数をXdX_dとすると、XdX_dとYYの間の相関係数は0.798ρ0.798 \rhoになる、つまり真の相関係数の約 0.8 倍へと過小評価する問題があることを報告している。

Cohen (1983)によれば、Peters and Van Voorhis (1940)が報告しているような二値化の希薄化の係数は次のように一般化できる

e=ϕ(z)p(1−p)e = \frac{\phi(z)}{\sqrt{p(1 - p)}}
(9)

  • ϕ(z)\phi(z)：標準正規分布の密度関数

  • pp：二値化した変数の比率

この希薄化誤差eeの逆数をpoint-biserialに掛けているのがbiserial

import numpy as np
from scipy.stats import pearsonr

# 2変量標準正規分布に従うデータ
rho = 0.75
size = 10000
np.random.seed(0)
data = np.random.multivariate_normal(mean=[0, 0], cov=[[1, rho], [rho, 1]], size=size)
x, y = data[:, 0], data[:, 1]

# 平均値で二値化する
xd = 1 * (x >= np.mean(x))  # 平均で2値化
rho_d = pearsonr(xd, y)[0]
print(f"真のρ={rho:.2f}, 離散化後のpearson={rho_d:.3f}, 比率={rho_d / rho:.3f}")


from scipy.stats import norm
p = xd.mean()
mean_z = 0  # 「標準正規分布」 & 「平均値で二分」の仮定より、閾値=平均値は0
e = norm.pdf(mean_z) / np.sqrt(p * (1 - p))
print(f"{e=:.3f}, 補正後（biserial）={rho_d / e:.3f}")
真のρ=0.75, 離散化後のpearson=0.601, 比率=0.801
e=0.798, 補正後（biserial）=0.753

rho = 0.5
cov = np.array([[1, rho], [rho, 1]])
X = multivariate_normal.rvs(cov=cov, size=100, random_state=0)
x = X[:, 0]
y = 1 * (X[:, 1] > 0.5)
print(f"biserial: {biserial_correlation(x, y):.5f}")
print(f"point_biserial: {point_biserial_correlation(x, y):.5f}")
biserial: 0.42521
point_biserial: 0.32762

from ordinalcorr import polyserial
polyserial(x, y)
np.float64(0.4539866448381744)

補正係数はどういう関数になっているのか

相関の希薄化（attenuation）の効果は、二値化の分割点以上の値の比率p=E⁡[1(Y≥τ)]p = \operatorname{E}[\mathbb{1}(Y \geq \tau)]の関数である

（閾値の標準正規分布上の位置zzが閾値以上の値の比率ppの関数であるため）

e(p)=ϕ(z)p(1−p)e(p) = \frac{\phi(z)}{\sqrt{p(1 - p)}}
(10)

e(p)e(p)はppが0あるいは1に近いときに極端に小さくなる。

biserial相関係数で使用する補正項は1/e(p)1/e(p)なので、補正係数が極端に大きくなることに相当する。

Source
import numpy as np
from scipy.stats import norm
c = 0.000001
p = np.linspace(0 + c, 1 - c, 500)

def attenuation_effect(p: np.ndarray) -> np.ndarray:
    z = norm.ppf(p)
    e = norm.pdf(z) / np.sqrt(p * (1 - p))
    return e

e = attenuation_effect(p)

import matplotlib.pyplot as plt
# fig, ax = plt.subplots(figsize=[4, 2])
fig, axes = plt.subplots(figsize=[9, 2], ncols=2)
ax = axes[0]
ax.plot(p, e)
ax.set(
    title=r"Attenuation effect $e(p) = \frac{\phi(z)}{\sqrt{p(1-p)}}$",
    ylabel=r"$e(p) = \frac{\phi(z)}{\sqrt{p(1-p)}}$",
    xlabel=r"proportion of population above split point $p$",
    ylim=[0,1]
)

ax = axes[1]
ax.plot(p, 1/e)
ax.set(
    title=r"Inverse of Attenuation effect $e(p)$",
    ylabel=r"$1 / e(p)$",
    xlabel=r"proportion of population above split point $p$",
)

fig.show()
<Figure size 900x200 with 2 Axes>

そのため、例えばppが極めて0や1に近い値のときは補正係数が大きくなりすぎてbiserial相関係数の絶対値が1を超えることもありうる。

import numpy as np
from scipy.stats import pearsonr

rho = -0.999
n = 100
np.random.seed(0)
data = np.random.multivariate_normal(mean=[0, 0], cov=[[1, rho], [rho, 1]], size=n)
x, y = data[:, 0], data[:, 1]

threshold = 0.999
yd = 1 * (y >= threshold)
p = yd.mean()
e = attenuation_effect(p)

r_pbi = pearsonr(x, yd).statistic
r_bi = r_pbi / e
print(f"{rho=:.3f}, {p=:.3f}, {r_pbi=:.3f}, {r_bi=:.3f}")
rho=-0.999, p=0.170, r_pbi=-0.714, r_bi=-1.060

なお、離散化前の潜在変数に正規分布を仮定しない場合は異なる補正係数となる（Bedrick, 1995）。

参考文献

統計学
ポリシリアル相関係数
統計学
スピアマンの順位相関係数（Spearman’s rho）