관측의 독립성 확인¶

import numpy as np
import pandas as pd
from statsmodels.formula.api import ols

# 이 페이지의 모든 진단은 아래 모형 하나를 놓고 수행한다.
# 집단마다 표준편차를 1.0, 1.3, 1.6으로 다르게 주어 진단이 무엇을 잡아내는지
# (그리고 무엇을 못 잡아내는지) 볼 수 있게 했다.
rng = np.random.default_rng(42)
n = 20
data = pd.DataFrame({
    "group": np.repeat(["A", "B", "C"], n),
    "response": np.concatenate([
        rng.normal(10.0, 1.0, n),
        rng.normal(10.8, 1.3, n),
        rng.normal(12.0, 1.6, n),
    ]),
})
model = ols("response ~ C(group)", data=data).fit()

print(data.groupby("group").response.agg(["count", "mean", "std"]).round(3))
print(f"\nF = {model.fvalue:.4f}, p = {model.f_pvalue:.4f}")

e = model.resid
print(f"\n잔차 전체의 합 = {e.sum():+.3e}")
print("집단별 잔차 합:", "  ".join(f"{g}: {v.sum():+.1e}"
                                for g, v in e.groupby(data["group"])))

# (2) 에서 유도한 잔차 사이의 상관을 모의로 확인한다.
print(f"\n이론  같은 집단 corr(e_i, e_j) = -1/(n-1) = {-1 / (n - 1):+.6f}")
print("      다른 집단 corr                        = +0.000000")

rng_sim = np.random.default_rng(123)
B = 2000
w1, w2, b1, b2 = [], [], [], []
for _ in range(B):
    eps = rng_sim.normal(0, 1, (3, n))          # 서로 독립인 오차
    r = eps - eps.mean(axis=1, keepdims=True)   # 집단평균만 뺀 잔차
    for g in range(3):
        w1.append(r[g, 0]); w2.append(r[g, 1])
    b1.append(r[0, 0]); b2.append(r[1, 0])
print(f"모의  같은 집단 corr = {np.corrcoef(w1, w2)[0, 1]:+.6f}"
      f"   ({3 * B}쌍, 오차 {1 / np.sqrt(3 * B):.4f})")
print(f"      다른 집단 corr = {np.corrcoef(b1, b2)[0, 1]:+.6f}"
      f"   ({B}쌍, 오차 {1 / np.sqrt(B):.4f})")

# 이 쪽에서 쓸 양: 잔차가 놓인 순서. 자료프레임은 집단별로 묶여 있다.
print("\n행 순서대로 본 집단 (처음 25행)")
print("  " + "".join(data["group"].values[:25]))

       count    mean    std
group                      
A         20   9.967  0.870
B         20  10.942  1.034
C         20  12.191  1.145

F = 23.7708, p = 0.0000

잔차 전체의 합 = -2.469e-13
집단별 잔차 합: A: -7.1e-14  B: -5.3e-15  C: -1.7e-13

이론  같은 집단 corr(e_i, e_j) = -1/(n-1) = -0.052632
      다른 집단 corr                        = +0.000000
모의  같은 집단 corr = -0.058448   (6000쌍, 오차 0.0129)
      다른 집단 corr = -0.004901   (2000쌍, 오차 0.0224)

행 순서대로 본 집단 (처음 25행)
  AAAAAAAAAAAAAAAAAAAABBBBB

import numpy as np
from statsmodels.stats.stattools import durbin_watson

# durbin_watson은 잔차를 **주어진 순서 그대로** 본다.
# 그래서 자료가 수집 순서대로 정렬되어 있어야 의미가 있다.
# 집단별로 정렬된 자료에 그냥 적용하면 집단 효과를 자기상관으로 오인할 수 있다.
dw_stat = durbin_watson(model.resid)
print(f"Durbin-Watson Statistic: {dw_stat:.4f}")

# (1) 의 정확한 분해 d = 2 - 2*rho - (e_1^2 + e_n^2)/S 을 확인한다.
e = model.resid.values
S = (e ** 2).sum()
rho = (e[1:] * e[:-1]).sum() / S
edge = (e[0] ** 2 + e[-1] ** 2) / S
print(f"rho-hat                        = {rho:+.8f}")
print(f"2(1 - rho-hat)                 = {2 * (1 - rho):.8f}   <- 흔히 쓰는 근사")
print(f"끝점 보정 (e_1^2 + e_n^2)/S    = {edge:.8f}")
print(f"2 - 2*rho - 끝점 보정          = {2 - 2 * rho - edge:.8f}")
print(f"durbin_watson 과의 차이        = {abs(2 - 2 * rho - edge - dw_stat):.3e}")

# (2) 완전한 등차수열 잔차에서는 d = 12 / (n(n+1)) 인가
print(f"\n{'n':>5s} {'직접 계산':>13s} {'12/(n(n+1))':>13s} {'차이':>9s}")
for m in [5, 10, 20, 50, 101]:
    a = np.arange(1, m + 1) - (m + 1) / 2     # 합이 0 인 등차수열
    d_direct, d_formula = durbin_watson(a), 12 / (m * (m + 1))
    print(f"{m:5d} {d_direct:13.10f} {d_formula:13.10f} {d_direct - d_formula:9.1e}")

Durbin-Watson Statistic: 2.1101
rho-hat                        = -0.07149872
2(1 - rho-hat)                 = 2.14299745   <- 흔히 쓰는 근사
끝점 보정 (e_1^2 + e_n^2)/S    = 0.03287376
2 - 2*rho - 끝점 보정          = 2.11012369
durbin_watson 과의 차이        = 0.000e+00

    n         직접 계산   12/(n(n+1))        차이
    5  0.4000000000  0.4000000000   0.0e+00
   10  0.1090909091  0.1090909091   0.0e+00
   20  0.0285714286  0.0285714286   0.0e+00
   50  0.0047058824  0.0047058824   0.0e+00
  101  0.0011648224  0.0011648224   0.0e+00

import numpy as np

# 집단별로 묶인 순서에서 d 의 중심이 정말 2(1 + (k-1)/N) 인가.
# 오차를 서로 독립으로 만들므로 귀무가설은 참이다.
rng_sim = np.random.default_rng(2026)
B = 20_000
print(f"{'k':>3s} {'n':>4s} {'모의 E[분자]/E[분모]':>20s} {'2(1+(k-1)/N)':>14s} {'차이':>9s}")
for k, m in [(2, 10), (3, 20), (5, 8), (4, 25)]:
    N = k * m
    num = np.empty(B)
    den = np.empty(B)
    for b in range(B):
        eps = rng_sim.normal(0, 1, (k, m))
        r = (eps - eps.mean(axis=1, keepdims=True)).ravel()   # 집단 순서로 늘어놓는다
        num[b] = (np.diff(r) ** 2).sum()
        den[b] = (r ** 2).sum()
    sim = num.mean() / den.mean()
    th = 2 * (1 + (k - 1) / N)
    print(f"{k:3d} {m:4d} {sim:20.5f} {th:14.5f} {sim - th:+9.5f}")

  k    n       모의 E[분자]/E[분모]   2(1+(k-1)/N)        차이
  2   10              2.10028        2.10000  +0.00028
  3   20              2.06733        2.06667  +0.00066
  5    8              2.20555        2.20000  +0.00555
  4   25              2.05902        2.06000  -0.00098

import numpy as np
import pandas as pd
from statsmodels.formula.api import ols
from statsmodels.stats.stattools import durbin_watson

# 자기상관이 **정말로 있는** 자료를 만든다. 측정은 수집 순서대로 AR(1) 오차를
# 받지만, 자료프레임에는 집단별로 묶여 저장된다. 실제 자료가 흔히 이 꼴이다.
rng_ar = np.random.default_rng(7)
N, rho_true = 60, 0.8
grp = np.repeat(["A", "B", "C"], 20)
mu = {"A": 10.0, "B": 10.8, "C": 12.0}

order = rng_ar.permutation(N)      # t 번째로 수집된 관측이 들어갈 행 번호
eps, x = np.empty(N), 0.0
for t in range(N):
    x = rho_true * x + rng_ar.normal(0, 1)
    eps[t] = x

y = np.empty(N)
for t in range(N):
    y[order[t]] = mu[grp[order[t]]] + eps[t]

df = pd.DataFrame({"group": grp, "response": y})
df["t"] = np.argsort(order)         # 각 행이 몇 번째로 수집되었는가
m_ar = ols("response ~ C(group)", data=df).fit()
r = m_ar.resid.values

print(f"참 자기상관 rho = {rho_true}")
print(f"집단순 잔차의 d = {durbin_watson(r):.4f}   <- 자료프레임 행 순서")
print(f"수집순 잔차의 d = {durbin_watson(r[np.argsort(df['t'].values)]):.4f}"
      f"   <- 측정 순서")

참 자기상관 rho = 0.8
집단순 잔차의 d = 1.8443   <- 자료프레임 행 순서
수집순 잔차의 d = 0.7927   <- 측정 순서

import matplotlib.pyplot as plt
import numpy as np
from scipy import stats
from statsmodels.stats.stattools import durbin_watson

# 그림에서 읽으려는 세 가지를 먼저 수로 적어 둔다.
e = model.resid.values
idx = np.arange(len(e))

# 추세: 잔차를 행 번호에 회귀한 기울기
sl = stats.linregress(idx, e)
print(f"추세   기울기 {sl.slope:+.6f} (표준오차 {sl.stderr:.6f}), p = {sl.pvalue:.4f}")

# 자기상관: 시차 1 상관과 d
rho1 = (e[1:] * e[:-1]).sum() / (e ** 2).sum()
print(f"자기상관  시차 1 rho-hat = {rho1:+.4f},  d = {durbin_watson(e):.4f}"
      f"  (집단순 중심 {2 * (1 + 2 / 60):.4f})")

# 군집: 부호의 런 수. 기대값 2*n1*n2/N + 1
sign = e > 0
runs = 1 + (sign[1:] != sign[:-1]).sum()
n1, n2 = sign.sum(), (~sign).sum()
N = n1 + n2
mu_r = 2 * n1 * n2 / N + 1
sd_r = np.sqrt(2 * n1 * n2 * (2 * n1 * n2 - N) / (N ** 2 * (N - 1)))
print(f"군집   양수 {n1}개, 음수 {n2}개,  런 {runs}개"
      f"  (기대 {mu_r:.2f}, 표준편차 {sd_r:.2f}, z = {(runs - mu_r) / sd_r:+.2f})")

# 잔차를 관측 순서대로 찍는다. 이웃한 점들이 같은 쪽으로 몰려 다니면
# 독립이 아니라는 신호다. 순서가 뜻을 갖는 자료(시간·공간)에서만 쓸 수 있다.
plt.scatter(range(len(model.resid)), model.resid, alpha=0.6)
plt.axhline(y=0, color='r', linestyle='--')
plt.xlabel("Observation Order")
plt.ylabel("Residuals")
plt.title("Residuals vs. Observation Order")
plt.show()

추세   기울기 -0.000381 (표준오차 0.007556), p = 0.9600
자기상관  시차 1 rho-hat = -0.0715,  d = 2.1101  (집단순 중심 2.0667)
군집   양수 32개, 음수 28개,  런 31개  (기대 30.87, 표준편차 3.82, z = +0.03)

import numpy as np
from scipy import stats
from statsmodels.stats.stattools import durbin_watson

# 보기 2 의 AR(1) 자료(참 rho = 0.8)를 두 순서로 같은 세 가지 수치로 재 본다.
def read_three(v, label):
    i = np.arange(len(v))
    sl = stats.linregress(i, v)
    sign = v > 0
    runs = 1 + (sign[1:] != sign[:-1]).sum()
    n1, n2 = sign.sum(), (~sign).sum()
    N = n1 + n2
    mu_r = 2 * n1 * n2 / N + 1
    sd_r = np.sqrt(2 * n1 * n2 * (2 * n1 * n2 - N) / (N ** 2 * (N - 1)))
    print(f"{label}  기울기 p = {sl.pvalue:.4f},  d = {durbin_watson(v):.4f},  "
          f"런 {runs}개 (기대 {mu_r:.1f}, z = {(runs - mu_r) / sd_r:+.2f})")

read_three(r, "집단순")
read_three(r[np.argsort(df["t"].values)], "수집순")

집단순  기울기 p = 0.4481,  d = 1.8443,  런 24개 (기대 30.9, z = -1.80)
수집순  기울기 p = 0.5023,  d = 0.7927,  런 14개 (기대 30.9, z = -4.41)

import numpy as np
from scipy import stats

rng = np.random.default_rng(1111)

def ar1(n, rho, rng):
    """주변분산이 1 인 AR(1) 계열."""
    e = rng.normal(0, np.sqrt(1 - rho**2), n)
    x = np.empty(n)
    x[0] = rng.normal()
    for t in range(1, n):
        x[t] = rho * x[t - 1] + e[t]
    return x

B = 8_000
print("k=3, 집단당 n=20, 모든 평균 0, 명목 0.05")
print(f"{'ρ':>6s} {'F 오류율':>9s} {'유효 n = 20(1-ρ)/(1+ρ)':>22s} "
      f"{'분산 팽창 (1+ρ)/(1-ρ)':>22s}")
for rho in [0.0, 0.2, 0.4, 0.6, 0.8]:
    a = sum(stats.f_oneway(*[ar1(20, rho, rng) for _ in range(3)]).pvalue < 0.05
            for _ in range(B))
    print(f"{rho:6.1f} {a / B:9.4f} {20 * (1 - rho) / (1 + rho):22.2f} "
          f"{(1 + rho) / (1 - rho):22.2f}")

k=3, 집단당 n=20, 모든 평균 0, 명목 0.05
     ρ     F 오류율   유효 n = 20(1-ρ)/(1+ρ)      분산 팽창 (1+ρ)/(1-ρ)
   0.0    0.0471                  20.00                   1.00
   0.2    0.1355                  13.33                   1.50
   0.4    0.2727                   8.57                   2.33
   0.6    0.4780                   5.00                   4.00
   0.8    0.7300                   2.22                   9.00

import numpy as np
from scipy import stats

rng = np.random.default_rng(1111)
B = 8_000
print("3집단 × 군집 5개 × 군집당 4명 (집단당 20명), 명목 0.05")
print(f"{'ICC':>6s} {'F 오류율':>9s} {'설계효과 1+(m-1)ICC':>18s}")
for icc in [0.0, 0.05, 0.1, 0.2, 0.4]:
    a = 0
    for _ in range(B):
        gs = []
        for g in range(3):
            vals = [rng.normal(0, np.sqrt(icc))
                    + rng.normal(0, np.sqrt(1 - icc), 4) for _ in range(5)]
            gs.append(np.concatenate(vals))
        a += stats.f_oneway(*gs).pvalue < 0.05
    print(f"{icc:6.2f} {a / B:9.4f} {1 + 3 * icc:18.2f}")

3집단 × 군집 5개 × 군집당 4명 (집단당 20명), 명목 0.05
   ICC     F 오류율    설계효과 1+(m-1)ICC
  0.00    0.0466               1.00
  0.05    0.0790               1.15
  0.10    0.1017               1.30
  0.20    0.1581               1.60
  0.40    0.2659               2.20

import numpy as np
from scipy import stats
from statsmodels.stats.stattools import durbin_watson

rng = np.random.default_rng(2222)
B = 4_000
N = 60          # 잔차 개수

def ar1(n, rho, rng):
    e = rng.normal(0, np.sqrt(1 - rho**2), n)
    x = np.empty(n)
    x[0] = rng.normal()
    for t in range(1, n):
        x[t] = rho * x[t - 1] + e[t]
    return x

print("(가) 시간적 자기상관을 잡는가  (기준: |d-2| > 2·1.96/√N)")
print(f"{'ρ':>6s} {'평균 d':>8s} {'탐지율':>8s}")
for rho in [0.0, 0.2, 0.4, 0.6, 0.8]:
    ds, hit = [], 0
    for _ in range(B):
        gs = [ar1(20, rho, rng) for _ in range(3)]
        r = np.concatenate([g - g.mean() for g in gs])
        d = durbin_watson(r)
        ds.append(d)
        hit += abs(d - 2) > 2 * 1.96 / np.sqrt(N)
    print(f"{rho:6.1f} {np.mean(ds):8.4f} {hit / B:8.4f}")

print("\n(나) 군집 자료(ICC=0.3)는 잡는가 — 수집 순서가 군집 순서와 무관할 때")
ds, hit = [], 0
for _ in range(B):
    r = []
    for g in range(3):
        vals = [rng.normal(0, np.sqrt(0.3))
                + rng.normal(0, np.sqrt(0.7), 4) for _ in range(5)]
        x = np.concatenate(vals)
        r.append(x - x.mean())
    res = np.concatenate(r)
    d = durbin_watson(res[rng.permutation(len(res))])
    ds.append(d)
    hit += abs(d - 2) > 2 * 1.96 / np.sqrt(N)
print(f"  평균 d = {np.mean(ds):.4f}, 탐지율 = {hit / B:.4f}   (ICC=0.3 인데도)")

(가) 시간적 자기상관을 잡는가  (기준: |d-2| > 2·1.96/√N)
     ρ     평균 d      탐지율
   0.0   2.0581   0.0512
   0.2   1.7149   0.1998
   0.4   1.3814   0.6855
   0.6   1.0426   0.9617
   0.8   0.7271   0.9990

(나) 군집 자료(ICC=0.3)는 잡는가 — 수집 순서가 군집 순서와 무관할 때
  평균 d = 2.0000, 탐지율 = 0.0473   (ICC=0.3 인데도)

import warnings
warnings.filterwarnings("ignore")

import numpy as np
import pandas as pd
from scipy import stats
import statsmodels.formula.api as smf

def run(ICC, C, M, B, seed):
    """C = 집단당 군집 수, M = 군집당 인원."""
    rng = np.random.default_rng(seed)
    a = b = c = 0
    for _ in range(B):
        rows, cid = [], 0
        for g in range(3):
            for _c in range(C):
                u = rng.normal(0, np.sqrt(ICC))
                rows.append(pd.DataFrame(
                    {"y": u + rng.normal(0, np.sqrt(1 - ICC), M),
                     "g": f"G{g}", "cl": cid}))
                cid += 1
        df = pd.concat(rows, ignore_index=True)

        a += stats.f_oneway(*[v.y.values
                              for _, v in df.groupby("g")]).pvalue < 0.05
        cm = df.groupby(["g", "cl"]).y.mean().reset_index()
        b += stats.f_oneway(*[v.y.values
                              for _, v in cm.groupby("g")]).pvalue < 0.05
        m = smf.mixedlm("y ~ C(g)", df, groups=df["cl"]).fit(reml=True)
        idx = list(m.params.index)
        names = [x for x in idx if x.startswith("C(g)")]
        R = np.zeros((len(names), len(idx)))
        for r, nm in enumerate(names):
            R[r, idx.index(nm)] = 1
        c += float(m.f_test(R).pvalue) < 0.05
    return a / B, b / B, c / B

print("ICC=0.3, 명목 0.05")
print(f"{'설계':>26s} {'개별 분산분석':>12s} {'군집 평균':>10s} {'혼합효과':>10s}")
for lab, C, M, B, seed in [("군집 5개/집단, 군집당 4", 5, 4, 800, 3333),
                           ("군집 15개/집단, 군집당 4", 15, 4, 600, 4444)]:
    a, b, c = run(0.3, C, M, B, seed)
    print(f"{lab:>26s} {a:12.4f} {b:10.4f} {c:10.4f}")

ICC=0.3, 명목 0.05
                        설계      개별 분산분석      군집 평균       혼합효과
           군집 5개/집단, 군집당 4       0.2188     0.0575     0.0862
          군집 15개/집단, 군집당 4       0.2183     0.0533     0.0583

import numpy as np
from scipy import stats

rng = np.random.default_rng(5555)
B = 20_000
print("제1종 오류율 (참 평균차 0, n=15 쌍)")
print(f"{'개체내 상관 φ':>12s} {'대응 t':>10s} {'독립 t':>10s}")
for phi in [0.0, 0.3, 0.6, 0.9]:
    a = b = 0
    for _ in range(B):
        u = rng.normal(0, np.sqrt(phi), 15)
        x = u + rng.normal(0, np.sqrt(1 - phi), 15)
        y = u + rng.normal(0, np.sqrt(1 - phi), 15)
        a += stats.ttest_rel(x, y).pvalue < 0.05
        b += stats.ttest_ind(x, y).pvalue < 0.05
    print(f"{phi:12.1f} {a / B:10.4f} {b / B:10.4f}")

rng = np.random.default_rng(5556)
print("\n검정력 (참 평균차 0.5, n=15 쌍)")
print(f"{'개체내 상관 φ':>12s} {'대응 t':>10s} {'독립 t':>10s}")
for phi in [0.0, 0.3, 0.6, 0.9]:
    a = b = 0
    for _ in range(B):
        u = rng.normal(0, np.sqrt(phi), 15)
        x = u + rng.normal(0, np.sqrt(1 - phi), 15)
        y = u + rng.normal(0.5, np.sqrt(1 - phi), 15)
        a += stats.ttest_rel(x, y).pvalue < 0.05
        b += stats.ttest_ind(x, y).pvalue < 0.05
    print(f"{phi:12.1f} {a / B:10.4f} {b / B:10.4f}")

제1종 오류율 (참 평균차 0, n=15 쌍)
    개체내 상관 φ       대응 t       독립 t
         0.0     0.0522     0.0508
         0.3     0.0507     0.0217
         0.6     0.0500     0.0034
         0.9     0.0534     0.0000

검정력 (참 평균차 0.5, n=15 쌍)
    개체내 상관 φ       대응 t       독립 t
         0.0     0.2455     0.2589
         0.3     0.3316     0.2275
         0.6     0.5255     0.1788
         0.9     0.9810     0.0888

import warnings
warnings.filterwarnings("ignore")

import numpy as np
from scipy import stats
import statsmodels.api as sm

def ar1(n, rho, rng):
    e = rng.normal(0, np.sqrt(1 - rho**2), n)
    x = np.empty(n)
    x[0] = rng.normal()
    for t in range(1, n):
        x[t] = rho * x[t - 1] + e[t]
    return x

rng = np.random.default_rng(6667)
B = 2_000
R = np.array([[1.0, -1, 0], [0, 1, -1]])
print("3집단, 시간순 관측, 명목 0.05")
print(f"{'n':>5s} {'ρ':>5s} {'표준 F':>9s} {'HAC 왈드 F':>11s} {'maxlags':>8s}")
for n in [20, 100, 400]:
    L = int(4 * (n * 3 / 100)**(2 / 9)) + 1      # 흔히 쓰는 경험 규칙
    for rho in [0.0, 0.6]:
        a = b = 0
        X = np.zeros((3 * n, 3))
        for i in range(3):
            X[i * n:(i + 1) * n, i] = 1
        for _ in range(B):
            gs = [ar1(n, rho, rng) for _ in range(3)]
            y = np.concatenate(gs)
            a += stats.f_oneway(*gs).pvalue < 0.05
            fit = sm.OLS(y, X).fit(cov_type="HAC",
                                   cov_kwds={"maxlags": L})
            b += float(fit.f_test(R).pvalue) < 0.05
        print(f"{n:5d} {rho:5.1f} {a / B:9.4f} {b / B:11.4f} {L:8d}")

3집단, 시간순 관측, 명목 0.05
    n     ρ      표준 F    HAC 왈드 F  maxlags
   20   0.0    0.0545      0.1535        4
   20   0.6    0.4925      0.3295        4
  100   0.0    0.0460      0.0720        6
  100   0.6    0.4475      0.1470        6
  400   0.0    0.0460      0.0535        7
  400   0.6    0.4860      0.1115        7

관측이 서로 독립인가?
    │
    ├─ 같은 개체를 여러 번 쟀는가 ──→ 반복측정/혼합효과
    ├─ 자연스러운 군집이 있는가  ──→ 혼합효과 또는 군집 평균
    ├─ 시간 순서가 있는가        ──→ 시계열 모형, HAC, GLS
    ├─ 공간적으로 가까운가       ──→ 공간 모형, 블록 설계
    └─ 모두 아니오               ──→ 표준 분산분석 가능