범주형 설명변수¶

import numpy as np
import pandas as pd
from sklearn.linear_model import LinearRegression

rng = np.random.default_rng(0)
n = 900

# 세 도시. 효과에 순서가 없다: 부산이 가장 높고 서울이 가장 낮다.
city = rng.integers(0, 3, n)
CITY = ["서울", "부산", "대구"]
true = np.array([10.0, 16.0, 12.0])
y = true[city] + rng.normal(0, 2, n)

print("집단별 표본평균")
for k, name in enumerate(CITY):
    print(f"  {name}  n={np.sum(city == k):3d}  평균 {y[city == k].mean():7.4f}"
          f"   (참값 {true[k]:.1f})")

# 정수 코딩: 도시 이름을 0, 1, 2 로 바꿔 넣는다.
X_int = city.reshape(-1, 1).astype(float)
m_int = LinearRegression().fit(X_int, y)
print(f"\n정수 코딩 (서울=0, 부산=1, 대구=2)")
print(f"  기울기 {m_int.coef_[0]:+.4f}, 절편 {m_int.intercept_:.4f}, "
      f"R^2 {m_int.score(X_int, y):.4f}")
print(f"  예측값: 서울 {m_int.predict([[0]])[0]:.4f}, "
      f"부산 {m_int.predict([[1]])[0]:.4f}, "
      f"대구 {m_int.predict([[2]])[0]:.4f}")

집단별 표본평균
  서울  n=287  평균  9.9962   (참값 10.0)
  부산  n=286  평균 15.8841   (참값 16.0)
  대구  n=327  평균 11.9519   (참값 12.0)

정수 코딩 (서울=0, 부산=1, 대구=2)
  기울기 +0.8759, 절편 11.6630, R^2 0.0547
  예측값: 서울 11.6630, 부산 12.5389, 대구 13.4148

# 집단평균과 집단크기만으로 정수 코딩의 기울기를 계산한다.
counts = np.array([np.sum(city == k) for k in range(3)])
means = np.array([y[city == k].mean() for k in range(3)])
codes = np.arange(3.0)
x_bar, y_bar = (codes * counts).sum() / n, y.mean()
S_xy = (counts * (codes - x_bar) * (means - y_bar)).sum()
S_xx = (counts * (codes - x_bar) ** 2).sum()
print(f"x-bar = {x_bar:.6f},  y-bar = {y_bar:.4f}")
print(f"S_xy = {S_xy:.4f},  S_xx = {S_xx:.4f}")
print(f"S_xy / S_xx     = {S_xy / S_xx:.6f}   (sklearn 기울기 {m_int.coef_[0]:.6f})")
print(f"y-bar - b*x-bar = {y_bar - S_xy / S_xx * x_bar:.6f}   (sklearn 절편 {m_int.intercept_:.6f})")
print()
print(f"R^2          = {m_int.score(X_int, y):.6f}")
print(f"corr(x, y)^2 = {np.corrcoef(city, y)[0, 1] ** 2:.6f}")
print()
total_ss = ((y - y_bar) ** 2).sum()
within_ss = sum(((y[city == k] - means[k]) ** 2).sum() for k in range(3))
between_ss = (counts * (means - y_bar) ** 2).sum()
rss_int = ((y - m_int.predict(X_int)) ** 2).sum()
print(f"총 제곱합 TSS        = {total_ss:.4f}")
print(f"집단간 BSS           = {between_ss:.4f}   BSS/TSS = {between_ss / total_ss:.6f}")
print(f"집단내 WSS           = {within_ss:.4f}   (집단평균 모형의 RSS)")
print(f"정수 코딩의 RSS      = {rss_int:.4f}")
print(f"등간격 제약이 낸 손실 = {rss_int - within_ss:.4f}  "
      f"(= TSS 의 {(rss_int - within_ss) / total_ss:.4f})")

x-bar = 1.044444,  y-bar = 12.5778
S_xy = 536.2542,  S_xx = 612.2222
S_xy / S_xx     = 0.875914   (sklearn 기울기 0.875914)
y-bar - b*x-bar = 11.662986   (sklearn 절편 11.662986)

R^2          = 0.054742
corr(x, y)^2 = 0.054742

총 제곱합 TSS        = 8580.4087
집단간 BSS           = 5167.1933   BSS/TSS = 0.602208
집단내 WSS           = 3413.2155   (집단평균 모형의 RSS)
정수 코딩의 RSS      = 8110.6959
등간격 제약이 낸 손실 = 4697.4805  (= TSS 의 0.5475)

# 첫 열(서울)을 빼서 서울을 기준범주로 삼는다.
D = np.eye(3)[city][:, 1:]
m_d = LinearRegression().fit(D, y)

print(f"더미 부호화 (기준 = 서울)")
print(f"  절편 {m_d.intercept_:.4f}  = 서울 평균")
for k, name in enumerate(CITY[1:], 1):
    print(f"  {name} 더미 계수 {m_d.coef_[k - 1]:+.4f}"
          f"  = {name} 평균 - 서울 평균 "
          f"= {y[city == k].mean() - y[city == 0].mean():+.4f}")
print(f"  R^2 {m_d.score(D, y):.4f}")

더미 부호화 (기준 = 서울)
  절편 9.9962  = 서울 평균
  부산 더미 계수 +5.8878  = 부산 평균 - 서울 평균 = +5.8878
  대구 더미 계수 +1.9557  = 대구 평균 - 서울 평균 = +1.9557
  R^2 0.6022

# 더미 모형의 예측값이 집단평균인지, R^2 가 BSS/TSS 인지 확인한다.
from scipy import stats

fitted = m_d.predict(D)
print("예측값이 집단평균과 같은가")
for k, name in enumerate(CITY):
    print(f"  {name}  예측 {fitted[city == k][0]:.6f}   평균 {y[city == k].mean():.6f}"
          f"   차이 {abs(fitted[city == k][0] - y[city == k].mean()):.2e}")
print(f"  예측값이 집단 안에서 모두 같은가? "
      f"{all(np.ptp(fitted[city == k]) < 1e-12 for k in range(3))}")
print()
rss_dummy = ((y - fitted) ** 2).sum()
print(f"더미 모형의 RSS = {rss_dummy:.4f}   WSS = {within_ss:.4f}"
      f"   차이 {abs(rss_dummy - within_ss):.2e}")
print(f"R^2 = 1 - RSS/TSS = {1 - rss_dummy / total_ss:.6f}")
print(f"      BSS/TSS     = {between_ss / total_ss:.6f}")
print()
F = (between_ss / 2) / (within_ss / (n - 3))
print(f"일원배치 F = (BSS/2) / (WSS/(n-3)) = {F:.4f}")
print(f"  p-값 = {stats.f.sf(F, 2, n - 3):.4g}   (자유도 2, {n - 3})")

예측값이 집단평균과 같은가
  서울  예측 9.996228   평균 9.996228   차이 1.60e-14
  부산  예측 15.884073   평균 15.884073   차이 4.62e-14
  대구  예측 11.951942   평균 11.951942   차이 2.66e-14
  예측값이 집단 안에서 모두 같은가? True

더미 모형의 RSS = 3413.2155   WSS = 3413.2155   차이 4.55e-13
R^2 = 1 - RSS/TSS = 0.602208
      BSS/TSS     = 0.602208

일원배치 F = (BSS/2) / (WSS/(n-3)) = 678.9745
  p-값 = 2.79e-180   (자유도 2, 897)

full = np.eye(3)[city]                        # 더미 3개 전부
X_bad = np.column_stack([np.ones(n), full])   # 절편 + 더미 3개
X_ok = np.column_stack([np.ones(n), full[:, 1:]])

print(f"절편 + 더미 3개 : 열 {X_bad.shape[1]}개, 계수(rank) "
      f"{np.linalg.matrix_rank(X_bad)}")
print(f"절편 + 더미 2개 : 열 {X_ok.shape[1]}개, 계수(rank) "
      f"{np.linalg.matrix_rank(X_ok)}")
print(f"\n더미 3개의 합이 절편 열과 같은가?  "
      f"{np.allclose(full.sum(axis=1), 1.0)}")

# 해가 유일하지 않다는 것을 직접 보인다.
b1, *_ = np.linalg.lstsq(X_bad, y, rcond=None)
shift = np.array([1.0, -1.0, -1.0, -1.0])     # 절편 +1, 더미 전부 -1
b2 = b1 + shift

print(f"\n해 1: {np.round(b1, 4)}")
print(f"해 2: {np.round(b2, 4)}")
print(f"두 해의 예측이 같은가? {np.allclose(X_bad @ b1, X_bad @ b2)}")
print(f"잔차제곱합  해1 {np.sum((y - X_bad @ b1) ** 2):.4f}   "
      f"해2 {np.sum((y - X_bad @ b2) ** 2):.4f}")

절편 + 더미 3개 : 열 4개, 계수(rank) 3
절편 + 더미 2개 : 열 3개, 계수(rank) 3

더미 3개의 합이 절편 열과 같은가?  True

해 1: [9.4581 0.5382 6.426  2.4939]
해 2: [10.4581 -0.4618  5.426   1.4939]
두 해의 예측이 같은가? True
잔차제곱합  해1 3413.2155   해2 3413.2155

# 영공간을 확인하고, lstsq 가 어느 해를 고르는지 본다.
print(f"X_bad @ shift 의 절대값 최댓값 = {np.abs(X_bad @ shift).max():.2e}"
      f"   (0 이면 shift 가 영공간에 있다)")
print(f"영공간의 차원 = 열 {X_bad.shape[1]}개 - rank {np.linalg.matrix_rank(X_bad)} = "
      f"{X_bad.shape[1] - np.linalg.matrix_rank(X_bad)}")
print()
print(f"해 1 의 노름 = {np.linalg.norm(b1):.6f}")
print(f"해 2 의 노름 = {np.linalg.norm(b2):.6f}")
print(f"해 1 과 shift 의 내적 = {b1 @ shift:.2e}   (0 이면 최소노름해다)")
norms = [(t, np.linalg.norm(b1 + t * shift)) for t in np.linspace(-2, 2, 401)]
t_best = min(norms, key=lambda p: p[1])[0]
print(f"노름을 최소로 하는 t = {t_best:.4f}")
print()
print("식별되지 않는 것: 개별 계수.  식별되는 것: 절편 + 더미 계수의 합")
for k, name in enumerate(CITY):
    print(f"  {name}  해1 {b1[0] + b1[k + 1]:.6f}   해2 {b2[0] + b2[k + 1]:.6f}"
          f"   집단평균 {y[city == k].mean():.6f}")
print()
print("차이도 식별된다 (기준범주 부호화가 쓰는 양)")
print(f"  부산 - 서울  해1 {b1[2] - b1[1]:+.6f}   해2 {b2[2] - b2[1]:+.6f}")
print(f"  대구 - 서울  해1 {b1[3] - b1[1]:+.6f}   해2 {b2[3] - b2[1]:+.6f}")

X_bad @ shift 의 절대값 최댓값 = 0.00e+00   (0 이면 shift 가 영공간에 있다)
영공간의 차원 = 열 4개 - rank 3 = 1

해 1 의 노름 = 11.715699
해 2 의 노름 = 11.885185
해 1 과 shift 의 내적 = 1.55e-14   (0 이면 최소노름해다)
노름을 최소로 하는 t = 0.0000

식별되지 않는 것: 개별 계수.  식별되는 것: 절편 + 더미 계수의 합
  서울  해1 9.996228   해2 9.996228   집단평균 9.996228
  부산  해1 15.884073   해2 15.884073   집단평균 15.884073
  대구  해1 11.951942   해2 11.951942   집단평균 11.951942

차이도 식별된다 (기준범주 부호화가 쓰는 양)
  부산 - 서울  해1 +5.887845   해2 +5.887845
  대구 - 서울  해1 +1.955714   해2 +1.955714

import numpy as np
from sklearn.linear_model import LinearRegression

rng = np.random.default_rng(1)
n_edu = 6000
level = rng.integers(0, 5, n_edu)                 # 교육수준 5단계
true_edu = np.array([0.0, 0.2, 0.5, 3.0, 3.4])    # 임금 효과 — 등간격이 아니다
y_edu = true_edu[level] + rng.normal(0, 1, n_edu)

X_int = level.reshape(-1, 1).astype(float)        # 정수 인코딩
X_hot = np.eye(5)[level][:, 1:]                   # 원-핫 (기준 범주 제외)

for name, X in [("정수 인코딩", X_int), ("원-핫 인코딩", X_hot)]:
    m = LinearRegression().fit(X, y_edu)
    print(f"{name:>14}: R^2 {m.score(X, y_edu):.4f}"
          f"   계수 {np.round(m.coef_, 3)}")
print(f"{'참 효과':>14}: (기준 대비) {np.round(true_edu - true_edu[0], 3)}")

        정수 인코딩: R^2 0.5975   계수 [0.977]
       원-핫 인코딩: R^2 0.6933   계수 [0.226 0.552 3.09  3.48 ]
          참 효과: (기준 대비) [0.  0.2 0.5 3.  3.4]

print("기준범주를 바꾸면")
for ref in range(3):
    keep = [k for k in range(3) if k != ref]
    D = np.eye(3)[city][:, keep]
    m = LinearRegression().fit(D, y)
    parts = ", ".join(f"{CITY[k]} {c:+.4f}" for k, c in zip(keep, m.coef_))
    print(f"  기준={CITY[ref]:3s} 절편 {m.intercept_:8.4f} | {parts} "
          f"| R^2 {m.score(D, y):.4f}")

print("\n어느 기준을 쓰든 예측 집단평균은 같다")
for ref in range(3):
    keep = [k for k in range(3) if k != ref]
    D = np.eye(3)[city][:, keep]
    m = LinearRegression().fit(D, y)
    mu = {ref: m.intercept_}
    for k, c in zip(keep, m.coef_):
        mu[k] = m.intercept_ + c
    print(f"  기준={CITY[ref]:3s} " +
          ", ".join(f"{CITY[k]} {mu[k]:.4f}" for k in range(3)))

기준범주를 바꾸면
  기준=서울  절편   9.9962 | 부산 +5.8878, 대구 +1.9557 | R^2 0.6022
  기준=부산  절편  15.8841 | 서울 -5.8878, 대구 -3.9321 | R^2 0.6022
  기준=대구  절편  11.9519 | 서울 -1.9557, 부산 +3.9321 | R^2 0.6022

어느 기준을 쓰든 예측 집단평균은 같다
  기준=서울  서울 9.9962, 부산 15.8841, 대구 11.9519
  기준=부산  서울 9.9962, 부산 15.8841, 대구 11.9519
  기준=대구  서울 9.9962, 부산 15.8841, 대구 11.9519

# 효과 부호화: 부산 -> (1,0), 대구 -> (0,1), 서울(기준) -> (-1,-1)
E = np.zeros((n, 2))
E[city == 1, 0] = 1
E[city == 2, 1] = 1
E[city == 0, :] = -1

m = LinearRegression().fit(E, y)
group_means = [y[city == k].mean() for k in range(3)]

print(f"절편 {m.intercept_:.4f}")
print(f"  세 집단 평균의 평균 {np.mean(group_means):.4f}")
print(f"계수 {np.round(m.coef_, 4)}   (부산, 대구)")
print(f"  서울 효과 = -(계수의 합) = {-m.coef_.sum():+.4f}")
print(f"\n각 집단 평균 - 전체(집단평균의 평균)")
for k, name in enumerate(CITY):
    print(f"  {name} {group_means[k] - np.mean(group_means):+.4f}")

절편 12.6107
  세 집단 평균의 평균 12.6107
계수 [ 3.2733 -0.6588]   (부산, 대구)
  서울 효과 = -(계수의 합) = -2.6145

각 집단 평균 - 전체(집단평균의 평균)
  서울 -2.6145
  부산 +3.2733
  대구 -0.6588

import numpy as np
from sklearn.linear_model import LinearRegression

rng = np.random.default_rng(3)
N = 2000
SIGMA = 3.0

# 참 효과가 전부 0 이므로, 계수의 제곱평균제곱근이 곧 '0 에서 벗어난 정도'다.
print(f"관측 {N}개를 K개 수준에 고르게 나눈다 (참 효과는 전부 0)")
print(f"{'K':>5s}{'수준당 n':>10s}{'모수':>6s}"
      f"{'훈련 R^2':>10s}{'계수 RMS':>10s}")
for K in [2, 5, 20, 100, 400]:
    g = rng.integers(0, K, N)
    y0 = rng.normal(0, SIGMA, N)          # 집단 효과가 전혀 없는 자료
    D = np.eye(K)[g][:, 1:]
    m = LinearRegression().fit(D, y0)
    rms = np.sqrt((m.coef_ ** 2).mean())
    print(f"{K:>5d}{N / K:>10.1f}{K:>6d}"
          f"{m.score(D, y0):>10.4f}{rms:>10.4f}")

관측 2000개를 K개 수준에 고르게 나눈다 (참 효과는 전부 0)
    K     수준당 n    모수    훈련 R^2    계수 RMS
    2    1000.0     2    0.0002    0.0913
    5     400.0     5    0.0007    0.1731
   20     100.0    20    0.0091    0.3379
  100      20.0   100    0.0518    0.7033
  400       5.0   400    0.2049    1.5411