점이연 상관과 파이 계수¶

import numpy as np
from scipy import stats

# 점이연상관: 한쪽이 이분변수(0/1), 다른 쪽이 연속변수인 경우다.
group = np.array([0, 0, 0, 0, 0, 1, 1, 1, 1, 1])
scores = np.array([65, 70, 68, 72, 66, 78, 82, 85, 80, 76])

r_pb, p_val = stats.pointbiserialr(group, scores)
print(f"Point-biserial r = {r_pb:.4f}, p-value = {p_val:.4f}")

# 사실 점이연상관은 0/1 에 그대로 피어슨을 쓴 것과 같다. 이름만 다를 뿐이다.
r_pearson, _ = stats.pearsonr(group, scores)
print(f"Pearson r        = {r_pearson:.4f}")

# 파이계수: 양쪽이 모두 이분변수인 경우다. 2x2 표에서 카이제곱을 전체
# 도수로 나눈 뒤 제곱근을 취하면 나온다. 이 또한 0/1 에 대한 피어슨과 같다.
table = np.array([[60, 30], [40, 70]])
chi2, p, dof, expected = stats.chi2_contingency(table, correction=False)
phi = np.sqrt(chi2 / table.sum())
print(f"Phi coefficient  = {phi:.4f}")

# (1) 점이연 공식과 직접 비교
n0, n1, N = 5, 5, 10
m0, m1 = scores[:5].mean(), scores[5:].mean()
S_Y = scores.std(ddof=0)
print(f"\n(1) (m1-m0)/S_Y * sqrt(n0*n1/n^2)"
      f" = ({m1:.1f}-{m0:.1f})/{S_Y:.6f} * {np.sqrt(n0 * n1 / N**2):.1f}"
      f" = {(m1 - m0) / S_Y * np.sqrt(n0 * n1 / N**2):.6f}")

# (2) 파이계수는 0/1 자료에 그대로 피어슨을 쓴 것이다
(a, b), (c, d) = table
xs = np.r_[np.ones(a + b), np.zeros(c + d)]
ys = np.r_[np.ones(a), np.zeros(b), np.ones(c), np.zeros(d)]
print(f"\n(2) 0/1 자료의 피어슨 {stats.pearsonr(xs, ys)[0]:.6f}")
print(f"    (ad-bc)/sqrt(...)  "
      f"{(a * d - b * c) / np.sqrt((a + b) * (c + d) * (a + c) * (b + d)):.6f}")
print(f"    chi2 = {chi2:.6f}   n*phi^2 = {table.sum() * phi**2:.6f}"
      f"   200/11 = {200 / 11:.6f}")

# (3) 모집단 꼴과 천장
sigma = np.sqrt((((scores[:5] - m0)**2).sum() + ((scores[5:] - m1)**2).sum()) / N)
d_eff = (m1 - m0) / sigma
print(f"\n(3) sigma = {sigma:.6f}  d = {d_eff:.6f}"
      f"  d/sqrt(d^2+4) = {d_eff / np.sqrt(d_eff**2 + 4):.6f}")
pp, qq = (a + b) / table.sum(), (a + c) / table.sum()
ceiling = np.sqrt(pp * (1 - qq) / (qq * (1 - pp)))
print(f"    phi = 1/sqrt(11) = {1 / np.sqrt(11):.6f}"
      f"   천장 3/sqrt(11) = {ceiling:.6f}   비 {phi / ceiling:.6f}")

Point-biserial r = 0.9029, p-value = 0.0003
Pearson r        = 0.9029
Phi coefficient  = 0.3015

(1) (m1-m0)/S_Y * sqrt(n0*n1/n^2) = (80.2-68.2)/6.645299 * 0.5 = 0.902894

(2) 0/1 자료의 피어슨 0.301511
    (ad-bc)/sqrt(...)  0.301511
    chi2 = 18.181818   n*phi^2 = 18.181818   200/11 = 18.181818

(3) sigma = 2.856571  d = 4.200840  d/sqrt(d^2+4) = 0.902894
    phi = 1/sqrt(11) = 0.301511   천장 3/sqrt(11) = 0.904534   비 0.333333

import numpy as np

print("두 집단이 아무리 떨어져도 (d 를 키워도) |r_pb| 의 한계")
print(f"{'p (한쪽 비율)':>13s} {'d=2':>8s} {'d=5':>8s} {'d=20':>8s} {'d→∞':>8s}")
for p in [0.5, 0.4, 0.3, 0.2, 0.1, 0.05]:
    q = 1 - p
    vals = [d * np.sqrt(p * q) / np.sqrt(d**2 * p * q + 1) for d in [2, 5, 20]]
    print(f"{p:13.2f} {vals[0]:8.4f} {vals[1]:8.4f} {vals[2]:8.4f} {1.0:8.4f}")

두 집단이 아무리 떨어져도 (d 를 키워도) |r_pb| 의 한계
    p (한쪽 비율)      d=2      d=5     d=20      d→∞
         0.50   0.7071   0.9285   0.9950   1.0000
         0.40   0.6999   0.9258   0.9948   1.0000
         0.30   0.6757   0.9165   0.9941   1.0000
         0.20   0.6247   0.8944   0.9923   1.0000
         0.10   0.5145   0.8321   0.9864   1.0000
         0.05   0.3996   0.7368   0.9747   1.0000

rng = np.random.default_rng(23001)
print("\n두 집단이 50 표준편차 떨어진 (사실상 완전 분리) 자료")
print(f"{'n1':>4s} {'n2':>4s} {'p':>6s} {'r_pb':>8s}")
for n1, n2 in [(50, 50), (40, 60), (20, 80), (10, 90), (5, 95)]:
    x = np.r_[np.zeros(n1), np.ones(n2)]
    y = np.r_[rng.normal(0, 1, n1), rng.normal(50, 1, n2)]
    print(f"{n1:4d} {n2:4d} {n2 / (n1 + n2):6.2f} {np.corrcoef(x, y)[0, 1]:8.4f}")

두 집단이 50 표준편차 떨어진 (사실상 완전 분리) 자료
  n1   n2      p     r_pb
  50   50   0.50   0.9992
  40   60   0.60   0.9991
  20   80   0.80   0.9985
  10   90   0.90   0.9980
   5   95   0.95   0.9961

import numpy as np
from scipy import stats

rng = np.random.default_rng(23001)
print(f"{'n1':>4s} {'n2':>4s} {'r_pb':>11s} {'t/√(t²+df)':>13s} {'차이':>10s}")
for _ in range(5):
    n1, n2 = rng.integers(10, 40, 2)
    a = rng.normal(0, 1, n1)
    b = rng.normal(0.8, 1, n2)
    x = np.r_[np.zeros(n1), np.ones(n2)]
    y = np.r_[a, b]
    r = np.corrcoef(x, y)[0, 1]
    t = stats.ttest_ind(a, b, equal_var=True).statistic
    df = n1 + n2 - 2
    ident = -t / np.sqrt(t**2 + df)          # b 가 x=1 이므로 부호를 맞춘다
    print(f"{n1:4d} {n2:4d} {r:+11.6f} {ident:+13.6f} {abs(r - ident):10.2e}")

  n1   n2        r_pb    t/√(t²+df)         차이
  21   19   +0.448066     +0.448066   1.11e-16
  22   17   +0.523646     +0.523646   2.22e-16
  21   37   +0.273406     +0.273406   1.11e-16
  39   26   +0.280764     +0.280764   1.11e-16
  37   15   +0.356779     +0.356779   0.00e+00

import numpy as np
from scipy import stats

print("행 비율 p, 열 비율 q 일 때 φ 의 최댓값")
print(f"{'p':>6s} {'q':>6s} {'φ_max':>8s}")
for p, q in [(0.5, 0.5), (0.5, 0.3), (0.5, 0.1),
             (0.3, 0.3), (0.2, 0.8), (0.1, 0.9)]:
    phi = (min(p, q) - p * q) / np.sqrt(p * (1 - p) * q * (1 - q))
    print(f"{p:6.1f} {q:6.1f} {phi:8.4f}")

print("\n실제 표에서 확인 (n=1000, 최대한 겹치게 배치)")
for p, q in [(0.5, 0.5), (0.5, 0.2), (0.2, 0.8)]:
    n = 1000
    a = int(min(p, q) * n)
    T = np.array([[a, int(p * n) - a],
                  [int(q * n) - a, n - a - (int(p * n) - a) - (int(q * n) - a)]])
    chi = stats.chi2_contingency(T, correction=False).statistic
    print(f"  p={p}, q={q}: 표={T.tolist()}, φ=√(χ²/n)={np.sqrt(chi / n):.4f}")

행 비율 p, 열 비율 q 일 때 φ 의 최댓값
     p      q    φ_max
   0.5    0.5   1.0000
   0.5    0.3   0.6547
   0.5    0.1   0.3333
   0.3    0.3   1.0000
   0.2    0.8   0.2500
   0.1    0.9   0.1111

실제 표에서 확인 (n=1000, 최대한 겹치게 배치)
  p=0.5, q=0.5: 표=[[500, 0], [0, 500]], φ=√(χ²/n)=1.0000
  p=0.5, q=0.2: 표=[[200, 300], [0, 500]], φ=√(χ²/n)=0.5000
  p=0.2, q=0.8: 표=[[200, 0], [600, 200]], φ=√(χ²/n)=0.2500

p=0.1, q=0.9 이면 (n=1000)
  행 합:  100 / 900
  열 합:  900 / 100

  "완벽한 연관"을 만들려 해도 100 과 900 을 맞출 수 없다
  → 최선이 φ = 0.111

import warnings
warnings.filterwarnings("ignore")

import numpy as np

rng = np.random.default_rng(23002)
print(f"{'ρ':>6s} {'원래 r':>8s} {'x 만 분할':>10s} {'이론 ρ√(2/π)':>13s} "
      f"{'둘 다 분할':>11s} {'이론 (2/π)arcsin ρ':>19s}")
for rho in [0.2, 0.4, 0.6, 0.8]:
    z = rng.standard_normal((300_000, 2))
    x = z[:, 0]
    y = rho * z[:, 0] + np.sqrt(1 - rho**2) * z[:, 1]
    xd = (x > np.median(x)).astype(float)
    yd = (y > np.median(y)).astype(float)
    print(f"{rho:6.2f} {np.corrcoef(x, y)[0, 1]:8.4f} "
          f"{np.corrcoef(xd, y)[0, 1]:10.4f} {rho * np.sqrt(2 / np.pi):13.4f} "
          f"{np.corrcoef(xd, yd)[0, 1]:11.4f} "
          f"{2 / np.pi * np.arcsin(rho):19.4f}")

     ρ     원래 r     x 만 분할    이론 ρ√(2/π)      둘 다 분할    이론 (2/π)arcsin ρ
  0.20   0.2003     0.1584        0.1596      0.1268              0.1282
  0.40   0.4005     0.3219        0.3192      0.2644              0.2620
  0.60   0.5970     0.4754        0.4787      0.4064              0.4097
  0.80   0.8004     0.6383        0.6383      0.5911              0.5903

from scipy import stats

rng = np.random.default_rng(23002)
B = 6_000
a = b = c = 0
for _ in range(B):
    z = rng.standard_normal((100, 2))
    x = z[:, 0]
    y = 0.3 * z[:, 0] + np.sqrt(0.91) * z[:, 1]
    a += stats.pearsonr(x, y).pvalue < 0.05
    xd = (x > np.median(x)).astype(float)
    b += stats.pearsonr(xd, y).pvalue < 0.05
    yd = (y > np.median(y)).astype(float)
    c += stats.pearsonr(xd, yd).pvalue < 0.05
print("검정력 (ρ=0.3, n=100, 명목 0.05)")
print(f"  연속 그대로:   {a / B:.4f}")
print(f"  x 만 이분화:   {b / B:.4f}")
print(f"  둘 다 이분화:  {c / B:.4f}")

검정력 (ρ=0.3, n=100, 명목 0.05)
  연속 그대로:   0.8647
  x 만 이분화:   0.6760
  둘 다 이분화:  0.5563

n = 100 (연속) 의 정보를
이분화하면 n ≈ 250 이 있어야 회복된다
→ 자료의 60% 를 버린 셈

import warnings
warnings.filterwarnings("ignore")

import numpy as np
from scipy import stats
from scipy.optimize import brentq

def tetrachoric(T):
    """2x2 표에서 사분상관을 수치적으로 구한다."""
    n = T.sum()
    h = stats.norm.ppf((T[0, 0] + T[0, 1]) / n)
    k = stats.norm.ppf((T[0, 0] + T[1, 0]) / n)
    target = T[0, 0] / n

    def f(r):
        return stats.multivariate_normal.cdf(
            [h, k], mean=[0, 0], cov=[[1, r], [r, 1]]) - target

    return brentq(f, -0.999, 0.999)

rng = np.random.default_rng(23002)
print(f"{'참 ρ':>7s} {'φ (0 에서 분할)':>15s} {'사분상관':>10s}")
for rho in [0.3, 0.5, 0.7]:
    z = rng.standard_normal((200_000, 2))
    x = z[:, 0]
    y = rho * z[:, 0] + np.sqrt(1 - rho**2) * z[:, 1]
    xd, yd = (x > 0).astype(int), (y > 0).astype(int)
    T = np.array([[((xd == 0) & (yd == 0)).sum(), ((xd == 0) & (yd == 1)).sum()],
                  [((xd == 1) & (yd == 0)).sum(), ((xd == 1) & (yd == 1)).sum()]])
    print(f"{rho:7.2f} {np.corrcoef(xd, yd)[0, 1]:15.4f} {tetrachoric(T):10.4f}")

    참 ρ     φ (0 에서 분할)       사분상관
   0.30          0.1916     0.2964
   0.50          0.3380     0.5063
   0.70          0.4962     0.7029

이분화로 잃은 것
  ├─ 점추정의 축소  → 사분상관이 되돌린다  ✓
  └─ 정보(정밀도)   → 되돌릴 수 없다       ✗

파이썬:  from semopy import ...  또는
         factor_analyzer.factor_analyzer.corr  등
R:       psych::tetrachoric(), psych::polychoric()

자료가 이진인가?
    │
    ├─ 원래부터 이진 (합격/불합격, 생존/사망)
    │     ├─ 상대 변수도 이진  → φ 또는 오즈비
    │     │     └─ 잠재 연속이 자연스러운가?  → 사분상관
    │     └─ 상대 변수가 연속  → 점이연 r_pb (+ 코언 d 병기)
    │
    └─ 연속인데 이분화하려 한다
          ├─ 임상적 절단점이 미리 정해져 있다  → 허용
          └─ 그렇지 않다                        → 하지 말 것

처치군(n=45)과 대조군(n=55)의 회복 점수

  점이연 r = 0.42,  95% CI [0.24, 0.57]
  코언 d   = 0.93   (집단 비율 p = 0.45)
  t(98)    = 4.60,  p < 0.001

집단 비율이 균형에 가까우므로 r 과 d 의 변환이 안정적이다.
다른 비율의 연구와 비교할 때는 d 를 사용했다.