상관 시각화¶

import numpy as np

np.random.seed(7)
n = 200

# 공통 요인 z1 을 섞는 정도를 달리해 상관 구조를 만든다.
# X2 는 X1 과 강한 양, X3 은 X1 과 중간 음, X4 는 어느 것과도 무관하다.
z1 = np.random.randn(n)
z2 = np.random.randn(n)

x1 = z1
x2 = 0.7 * z1 + 0.3 * z2
x3 = -0.5 * z1 + np.random.randn(n) * 0.8
x4 = np.random.randn(n)

data = np.column_stack([x1, x2, x3, x4])
labels = ['X1', 'X2', 'X3', 'X4']

# (1) 에서 손으로 구한 모집단 상관행렬
s2, s3 = np.sqrt(0.58), np.sqrt(0.89)
theory = np.array([
    [1.0,        0.7 / s2,          -0.5 / s3,          0.0],
    [0.7 / s2,   1.0,               -0.35 / (s2 * s3),  0.0],
    [-0.5 / s3, -0.35 / (s2 * s3),   1.0,               0.0],
    [0.0,        0.0,                0.0,               1.0],
])

sample = np.corrcoef(data, rowvar=False)

print("쌍      이론      표본      SE      z")
for i in range(4):
    for j in range(i + 1, 4):
        rho, r = theory[i, j], sample[i, j]
        se = (1 - rho ** 2) / np.sqrt(n)   # 표본상관의 표준오차
        print(f"{labels[i]}-{labels[j]}  {rho:+.4f}  {r:+.4f}  {se:.4f}  {(r - rho) / se:+.3f}")

쌍      이론      표본      SE      z
X1-X2  +0.9191  +0.9192  0.0110  +0.001
X1-X3  -0.5300  -0.5115  0.0508  +0.365
X1-X4  +0.0000  -0.0156  0.0707  -0.220
X2-X3  -0.4871  -0.4688  0.0539  +0.339
X2-X4  +0.0000  +0.0272  0.0707  +0.384
X3-X4  +0.0000  +0.0295  0.0707  +0.417

# rowvar=False 는 "행이 관측, 열이 변수"라는 뜻이다. 기본값은 그 반대이므로
# 자료행렬을 그대로 넣으면 엉뚱한 행렬이 나온다.
corr_matrix = np.corrcoef(data, rowvar=False)

print(f"rowvar=False : {corr_matrix.shape}")
print(f"대각합       : {np.trace(corr_matrix):.1f}")
print(f"고윳값       : {np.round(np.linalg.eigvalsh(corr_matrix)[::-1], 4)}")

wrong = np.corrcoef(data)          # 기본값 rowvar=True
print(f"\nrowvar 생략  : {wrong.shape}")
print(f"대각 성분    : {wrong[0, 0]:.1f}, 성분 범위 [{wrong.min():.2f}, 1.00], 대칭 {np.allclose(wrong, wrong.T)}")
print(f"계수(rank)   : {np.linalg.matrix_rank(wrong)}")

rowvar=False : (4, 4)
대각합       : 4.0
고윳값       : [2.2914 1.0035 0.6266 0.0786]

rowvar 생략  : (200, 200)
대각 성분    : 1.0, 성분 범위 [-1.00, 1.00], 대칭 True
계수(rank)   : 3

import matplotlib.pyplot as plt

# 발산형 색지도를 쓰고 vmin/vmax 를 -1 과 1 로 못박는다. 이래야 흰색이
# 정확히 0 에 놓여, 색만 보고도 부호와 세기를 읽을 수 있다.
k = data.shape[1]
fig, ax = plt.subplots(figsize=(6, 5))
im = ax.imshow(corr_matrix, cmap='RdBu_r', vmin=-1, vmax=1)
ax.set_xticks(range(k))
ax.set_yticks(range(k))
ax.set_xticklabels(labels)
ax.set_yticklabels(labels)

for i in range(k):
    for j in range(k):
        ax.text(j, i, f'{corr_matrix[i, j]:.2f}',
                ha='center', va='center', fontsize=11,
                color='white' if abs(corr_matrix[i, j]) > 0.5 else 'black')

fig.colorbar(im, ax=ax, fraction=0.046, pad=0.04)
ax.set_title('Correlation Heatmap')
plt.tight_layout()
plt.show()

# 그림에서 읽을 수치를 찍어 둔다.
off = corr_matrix[np.triu_indices(k, 1)]
print(f"비대각 6칸: {np.round(off, 4)}")
print(f"|r| < 0.05 인 칸 {np.sum(np.abs(off) < 0.05)}개, 그 범위 "
      f"[{off[np.abs(off) < 0.05].min():+.4f}, {off[np.abs(off) < 0.05].max():+.4f}]")
print(f"가장 진한 비대각 칸 |r| = {np.abs(off).max():.4f}, 대각 칸 = 1.0000")

비대각 6칸: [ 0.9192 -0.5115 -0.0156 -0.4688  0.0272  0.0295]
|r| < 0.05 인 칸 3개, 그 범위 [-0.0156, +0.0295]
가장 진한 비대각 칸 |r| = 0.9192, 대각 칸 = 1.0000

# 열지도는 숫자 하나로 요약하지만 산점도 행렬은 관계의 모양을 보여 준다.
# 상관계수가 같아도 모양이 다를 수 있으므로 둘을 함께 본다.
# 대각선에는 그 변수 자신의 분포를 그린다.
fig, axes = plt.subplots(k, k, figsize=(10, 10))
for i in range(k):
    for j in range(k):
        ax = axes[i, j]
        if i == j:
            ax.hist(data[:, i], bins=20, edgecolor='k', alpha=0.7)
        else:
            ax.scatter(data[:, j], data[:, i], s=8, alpha=0.5)
        if j == 0:
            ax.set_ylabel(labels[i])
        if i == k - 1:
            ax.set_xlabel(labels[j])
        if j != 0:
            ax.set_yticklabels([])
        if i != k - 1:
            ax.set_xticklabels([])

fig.suptitle('Scatter Matrix (Pair Plot)', fontsize=14, y=1.01)
plt.tight_layout()
plt.show()

# 표준화한 점구름의 축 비를 직접 재어 공식과 맞춘다.
print("쌍      |r|      공식      주성분으로 잰 값")
for i, j in [(0, 1), (0, 2), (1, 2), (2, 3)]:
    Z = np.column_stack([data[:, i], data[:, j]])
    Z = (Z - Z.mean(axis=0)) / Z.std(axis=0, ddof=1)
    lam = np.linalg.eigvalsh(np.cov(Z, rowvar=False))     # 오름차순
    r = abs(corr_matrix[i, j])
    print(f"{labels[i]}-{labels[j]}  {r:.4f}  {np.sqrt((1 + r) / (1 - r)):.4f}    "
          f"{np.sqrt(lam[1] / lam[0]):.4f}")

쌍      |r|      공식      주성분으로 잰 값
X1-X2  0.9192  4.8724    4.8724
X1-X3  0.5115  1.7589    1.7589
X2-X3  0.4688  1.6629    1.6629
X3-X4  0.0295  1.0300    1.0300

import numpy as np

np.random.seed(7)
n = 200
z1 = np.random.randn(n)
x1 = z1
x5 = x1 ** 2 + np.random.normal(0, 0.1, n)

r = np.corrcoef(x1, x5)[0, 1]
print(f"Pearson r(x1, x5) = {r:.4f}")

Pearson r(x1, x5) = -0.0493

import numpy as np
import matplotlib.pyplot as plt

def lower_triangle_heatmap(corr, labels):
    k = corr.shape[0]
    mask = np.triu(np.ones_like(corr, dtype=bool))
    masked = np.ma.array(corr, mask=mask)

    fig, ax = plt.subplots(figsize=(6, 5))
    im = ax.imshow(masked, cmap='RdBu_r', vmin=-1, vmax=1)
    ax.set_xticks(range(k))
    ax.set_yticks(range(k))
    ax.set_xticklabels(labels, rotation=45)
    ax.set_yticklabels(labels)

    for i in range(k):
        for j in range(i):
            ax.text(j, i, f'{corr[i, j]:.2f}',
                    ha='center', va='center', fontsize=11)

    fig.colorbar(im, ax=ax)
    ax.set_title('Lower Triangle Correlation Heatmap')
    plt.tight_layout()
    plt.show()

# 사용 예
np.random.seed(7)
data = np.random.randn(200, 4)
corr = np.corrcoef(data, rowvar=False)
lower_triangle_heatmap(corr, ['X1', 'X2', 'X3', 'X4'])

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(42)
n = 100
# 1번 무리: (2, 2) 둘레
c1 = np.random.normal(loc=[2, 2], scale=0.5, size=(n, 2))
# 2번 무리: (-2, 2) 둘레
c2 = np.random.normal(loc=[-2, 2], scale=0.5, size=(n, 2))
data = np.vstack([c1, c2])

r = np.corrcoef(data[:, 0], data[:, 1])[0, 1]
print(f"Pearson r = {r:.4f}")   # about -0.01

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 4))
ax1.scatter(data[:, 0], data[:, 1], s=10, alpha=0.6)
ax1.set_title(f'Scatter Plot (r = {r:.3f})')
ax1.set_xlabel('X1')
ax1.set_ylabel('X2')

corr = np.corrcoef(data, rowvar=False)
ax2.imshow(corr, cmap='RdBu_r', vmin=-1, vmax=1)
ax2.set_title('Heatmap')
plt.tight_layout()
plt.show()

Pearson r = -0.0117

import warnings
warnings.filterwarnings("ignore")

import numpy as np
from scipy.cluster.hierarchy import linkage, leaves_list
from scipy.spatial.distance import squareform

rng = np.random.default_rng(28001)
p = 12
groups = np.repeat([0, 1, 2], 4)              # 3개 블록
R = np.eye(p)
for i in range(p):
    for j in range(p):
        if i != j:
            R[i, j] = 0.75 if groups[i] == groups[j] else 0.05
X = rng.standard_normal((3_000, p)) @ np.linalg.cholesky(R).T
Rh = np.corrcoef(X, rowvar=False)

print(f"  참 구조: 변수 {p}개, 3개 블록 (블록 내 0.75, 블록 간 0.05)")
print(f"  무작위 순서: {rng.permutation(p).tolist()}")

Z = linkage(squareform(1 - np.abs(Rh), checks=False), method="average")
order = leaves_list(Z)
print(f"  계층 군집 순서: {order.tolist()}")
print(f"  각 위치의 참 블록: {groups[order].tolist()}")
runs = 1 + sum(groups[order][i] != groups[order][i - 1] for i in range(1, p))
print(f"  블록이 바뀌는 횟수 = {runs}  (완벽히 정렬되면 3)")

iu = np.triu_indices(p, 1)
same = groups[iu[0]] == groups[iu[1]]
print(f"\n  블록 내 평균 상관: {Rh[iu][same].mean():.4f} (참 0.75)")
print(f"  블록 간 평균 상관: {Rh[iu][~same].mean():.4f} (참 0.05)")

  참 구조: 변수 12개, 3개 블록 (블록 내 0.75, 블록 간 0.05)
  무작위 순서: [5, 4, 1, 0, 2, 8, 10, 11, 9, 7, 6, 3]
  계층 군집 순서: [9, 8, 10, 11, 3, 0, 1, 2, 4, 7, 5, 6]
  각 위치의 참 블록: [2, 2, 2, 2, 0, 0, 0, 0, 1, 1, 1, 1]
  블록이 바뀌는 횟수 = 3  (완벽히 정렬되면 3)

  블록 내 평균 상관: 0.7542 (참 0.75)
  블록 간 평균 상관: 0.0486 (참 0.05)

from scipy.cluster.hierarchy import linkage, leaves_list
from scipy.spatial.distance import squareform

d = 1 - np.abs(R)                       # 거리로 변환
Z = linkage(squareform(d, checks=False), method="average")
order = leaves_list(Z)
R_sorted = R[np.ix_(order, order)]      # 행과 열을 함께 재배열

import numpy as np

rng = np.random.default_rng(28001)
R2 = np.eye(6)
R2[np.triu_indices(6, 1)] = rng.uniform(0.02, 0.08, 15)
R2 = R2 + R2.T - np.eye(6)
off = R2[np.triu_indices(6, 1)]
print("모든 상관이 0.02~0.08 인 자료 (사실상 무관)")
print(f"  실제 범위: [{off.min():.4f}, {off.max():.4f}]")
print(f"  vmin/vmax 자동  → 가장 진한 칸이 {off.max():.2f} 인데 최대 색으로 표시")
print(f"  vmin=-1, vmax=1 → 모든 칸이 거의 흰색 (올바름)")

모든 상관이 0.02~0.08 인 자료 (사실상 무관)
  실제 범위: [0.0221, 0.0744]
  vmin/vmax 자동  → 가장 진한 칸이 0.07 인데 최대 색으로 표시
  vmin=-1, vmax=1 → 모든 칸이 거의 흰색 (올바름)

plt.imshow(R, cmap="RdBu_r", vmin=-1, vmax=1)
                 └─ 발산형   └─ 고정 범위

1. cmap 은 발산형 (RdBu_r, coolwarm, bwr)
2. vmin=-1, vmax=1 고정
3. 0 이 중앙(흰색)에 오도록
4. 색막대(colorbar)를 반드시 붙인다

전체 남성의 약 8% 가 적록색각 이상이다
  → RdBu (빨강-파랑) 는 비교적 안전
  → RdYlGn (빨강-노랑-초록) 은 위험

대안: 타원 그림 (모양으로 부호화, 다음 절)
      숫자 주석 병기

□ vmin=-1, vmax=1 로 고정했는가 (또는 범위를 명시했는가)
□ 발산형 색상표를 썼는가
□ 색막대가 있는가
□ 대각선(모두 1)을 가렸는가
□ 위쪽 삼각형을 가려 중복을 없앴는가

import warnings
warnings.filterwarnings("ignore")

import numpy as np

rng = np.random.default_rng(28002)
p, n = 5, 300_000
X = np.empty((n, p))
X[:, 0] = rng.standard_normal(n)
for j in range(1, p):
    X[:, j] = 0.8 * X[:, j - 1] + np.sqrt(1 - 0.64) * rng.standard_normal(n)

R = np.corrcoef(X, rowvar=False)
P = np.linalg.inv(R)                     # 정밀도 행렬
D = np.sqrt(np.diag(P))
PC = -P / np.outer(D, D)                 # 부분상관 행렬
np.fill_diagonal(PC, 1.0)

print("주변 상관행렬:")
print(np.round(R, 3))
print("\n부분 상관행렬 (나머지 전부 통제):")
print(np.round(PC, 3))
print(f"\n  주변: 1-5 도 {R[0, 4]:.3f} 로 강해 보인다")
print(f"  부분: 이웃만 남는다 (1-2={PC[0, 1]:.3f}), 1-5={PC[0, 4]:.3f}")

주변 상관행렬:
[[1.    0.799 0.64  0.512 0.41 ]
 [0.799 1.    0.801 0.641 0.513]
 [0.64  0.801 1.    0.8   0.64 ]
 [0.512 0.641 0.8   1.    0.8  ]
 [0.41  0.513 0.64  0.8   1.   ]]

부분 상관행렬 (나머지 전부 통제):
[[ 1.     0.623 -0.002  0.     0.   ]
 [ 0.623  1.     0.491 -0.001 -0.   ]
 [-0.002  0.491  1.     0.488  0.   ]
 [ 0.    -0.001  0.488  1.     0.624]
 [ 0.    -0.     0.     0.624  1.   ]]

  주변: 1-5 도 0.410 로 강해 보인다
  부분: 이웃만 남는다 (1-2=0.623), 1-5=0.000

import warnings
warnings.filterwarnings("ignore")

import numpy as np
from scipy import stats
from statsmodels.stats.multitest import multipletests

rng = np.random.default_rng(28002)
B, p, n = 3_000, 15, 40
m = p * (p - 1) // 2
raw = adj = 0
for _ in range(B):
    X = rng.standard_normal((n, p))              # 모두 독립
    r = np.corrcoef(X, rowvar=False)[np.triu_indices(p, 1)]
    t = r * np.sqrt((n - 2) / (1 - r**2))
    pv = 2 * stats.t.sf(np.abs(t), n - 2)
    raw += (pv < 0.05).sum()
    adj += multipletests(pv, alpha=0.05, method="fdr_bh")[0].sum()
print("독립 자료 p=15, n=40 에서 * 가 몇 개 찍히는가")
print(f"  쌍의 수 = {m}")
print(f"  보정 없이:        {raw / B:.2f} 개")
print(f"  FDR(BH) 보정 후:  {adj / B:.2f} 개")

독립 자료 p=15, n=40 에서 * 가 몇 개 찍히는가
  쌍의 수 = 105
  보정 없이:        5.21 개
  FDR(BH) 보정 후:  0.05 개

n = 1000 인 연구:  r = 0.07 도 p < 0.05  → **  찍힌다
n = 20 인 연구:    r = 0.40 도 p = 0.08  → 안 찍힌다

별표가 많은 그림 = 표본이 큰 연구
                 ≠ 관계가 강한 자료

mask = np.abs(R) < 0.2
R_show = np.where(mask, np.nan, R)
plt.imshow(R_show, cmap="RdBu_r", vmin=-1, vmax=1)
  → NaN 칸은 배경색으로 비어 보인다

□ vmin=-1, vmax=1 로 고정
□ 발산형 색상표 (RdBu_r, coolwarm)
□ 색막대 표시
□ 계층 군집으로 변수 재배열 (+ 순서를 그림 설명에 명시)
□ 위쪽 삼각형과 대각선 가리기
□ 변수가 20개 이하면 각 칸에 r 값 표기
□ 별표는 FDR 보정 후에만
□ 표본 크기 n 을 그림 설명에

열지도는 "어디를 볼지" 알려 주는 도구이지
"무엇이 있는지" 알려 주는 도구가 아니다
  → 흥미로운 칸을 찾으면 반드시 그 쌍의 산점도를 그린다

그림 3. 12개 생리 지표의 상관 열지도 (n = 248)

피어슨 상관. 변수는 평균 연결 계층 군집(거리 = 1-|r|)으로
재배열했으며 왼쪽 덴드로그램에 군집 구조를 표시했다.
색은 -1(파랑)에서 +1(빨강)으로 고정했다. 아래쪽 삼각형만
표시했고, |r| < 0.2 인 칸은 비워 두었다.
별표는 BH 절차로 FDR 5% 를 통제한 결과다.