회귀에서의 분산 검정¶

(1) 참값은 설계행렬만으로 정해진다. \(\hat{\boldsymbol\beta} = (\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'\mathbf{Y}\)에 \(\mathbf{Y} = \mathbf{X}\boldsymbol\beta + \boldsymbol\varepsilon\)을 넣으면

이다. \(E[\boldsymbol\varepsilon \mid \mathbf{X}] = \mathbf{0}\)이므로 이 줄에서 이미 불편성이 끝난다. 오차의 분산 구조가 어떻든 상관없다. 분산을 재면

이고, 이 자료에서는 \(\varepsilon_i = x_i z_i\)이므로

으로 완전히 알려져 있다. 그러므로 참 표준오차를 추정 없이 수로 적을 수 있다. 단순회귀에서는 손으로도 쓸 수 있다.

OLS 공식이 겨누는 값은 다르다. OLS는 \(\boldsymbol\Omega = \sigma^2\mathbf{I}\)로 믿고 \(\sigma^2\)을 \(s^2 = \mathbf{e}'\mathbf{e}/(n-2)\)로 추정한다. \(\mathbf{e} = (\mathbf{I}-\mathbf{H})\boldsymbol\varepsilon\)이므로

이다. 곧 \(s^2\)은 \(\sigma_i^2\)들을 지렛대로 가중한 평균 하나로 수렴한다. 등분산이면 \(\operatorname{tr}(\mathbf{I}-\mathbf{H}) = n-2\)라 정확히 \(\sigma^2\)이 되지만, 이분산이면 스칼라 하나로 대각이 다른 \(\boldsymbol\Omega\)를 흉내 내야 한다.

여기에 문제의 핵심이 있다. OLS 표준오차는 참 표준오차에 공통 배수를 곱한 꼴

인데 참값은 \(\sqrt{\left[(\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'\boldsymbol\Omega\mathbf{X}(\mathbf{X}'\mathbf{X})^{-1}\right]_{jj}}\)이다. 고를 수 있는 수가 하나뿐이므로 여러 계수를 동시에 맞힐 수 없고, 어느 쪽으로 틀릴지도 정해져 있지 않다. 연습문제 2 가 지렛대와 잔차의 상관으로 설명한 것이 이 대수의 내용이다.

(2) 수치로 확인한다. 설계행렬 \(\mathbf{X}\)는 보기의 것을 그대로 고정하고 오차만 4000번 다시 뽑는다. \(\mathbf{X}\)를 고정하는 것은 위의 모든 식이 \(\mathbf{X}\)에 조건을 건 식이기 때문이다.

import numpy as np
from scipy import stats
import statsmodels.api as sm
from statsmodels.stats.diagnostic import het_breuschpagan, het_white

rng = np.random.default_rng(42)
n = 100
X = rng.uniform(1, 10, size=n)

# 잡음에 X 를 곱해 분산이 X 에 비례해 커지도록 만든다. 전형적인 이분산이다.
epsilon = rng.normal(0, 1, size=n) * X
Y = 3 + 2 * X + epsilon

X_with_const = sm.add_constant(X)
model = sm.OLS(Y, X_with_const).fit()
residuals = model.resid

# Breusch-Pagan 은 분산이 설명변수의 선형함수로 커지는 경우를 잘 잡는다.
bp_stat, bp_p, bp_f, bp_fp = het_breuschpagan(residuals, X_with_const)
print(f"Breusch-Pagan: LM = {bp_stat:.4f}, p = {bp_p:.6f}")

# White 는 제곱항과 교차항까지 넣어 비선형 형태의 이분산도 잡는다.
w_stat, w_p, w_f, w_fp = het_white(residuals, X_with_const)
print(f"White:         LM = {w_stat:.4f}, p = {w_p:.6f}")

# 이분산이 있어도 계수 추정값 자체는 여전히 불편이다. 망가지는 것은
# 표준오차다. HC3 로버스트 표준오차는 이분산을 셈에 넣어 계산하므로,
# 모형을 바꾸지 않고도 추론을 바로잡을 수 있다.
# 아래 출력에서 계수는 그대로이고 표준오차와 t 값만 달라지는 것을 본다.
robust_model = model.get_robustcov_results(cov_type='HC3')
print(f"\ncoefficients:      {np.round(model.params, 4)}")
print(f"OLS std errors:    {np.round(model.bse, 4)}")
print(f"Robust std errors: {np.round(robust_model.bse, 4)}")
print(f"OLS t-values:      {np.round(model.tvalues, 3)}")
print(f"Robust t-values:   {np.round(robust_model.tvalues, 3)}")

출력:

Breusch-Pagan: LM = 14.2174, p = 0.000163
White:         LM = 14.7165, p = 0.000637

coefficients:      [3.1523 1.9672]
OLS std errors:    [1.5014 0.254 ]
Robust std errors: [1.0664 0.2719]
OLS t-values:      [2.1   7.746]
Robust t-values:   [2.956 7.235]

두 검정 모두 이분산을 강하게 탐지한다(\(p = 0.00016\)과 \(p = 0.00064\)). 자료를 \(\operatorname{Var}(\varepsilon_i) \propto X_i^2\)이 되도록 생성했으므로 당연한 결과이다. 그런데 이 한 벌만 보면 어느 표준오차가 옳은지 알 수 없다. \(1.5014\)와 \(1.0664\) 가운데 무엇이 참값에 가까운가. (1)의 식과 반복추출이 그 답을 준다.

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

# 위 보기의 설계행렬을 그대로 쓴다. X 는 고정하고 오차만 다시 뽑는다.
rng = np.random.default_rng(42)
n = 100
X = rng.uniform(1, 10, size=n)
Xc = sm.add_constant(X)

# (1) 참 공분산행렬을 설계행렬만으로 정확히 계산한다.
#     Var(eps_i) = X_i^2 이므로 Omega = diag(X^2) 다. 모수 추정이 전혀 없다.
XtX_inv = np.linalg.inv(Xc.T @ Xc)
Omega = np.diag(X ** 2)
V_true = XtX_inv @ Xc.T @ Omega @ Xc @ XtX_inv
se_true = np.sqrt(np.diag(V_true))

# OLS 공식이 겨누는 값:  sigma^2 (X'X)^-1 에서 sigma^2 을 s^2 으로 바꾼 것.
# E[s^2] = tr((I-H) Omega)/(n-2) 이므로 그것을 넣는다.
H = Xc @ XtX_inv @ Xc.T
Es2 = np.trace((np.eye(n) - H) @ Omega) / (n - 2)
se_ols_target = np.sqrt(Es2 * np.diag(XtX_inv))

print(f"참 표준오차 (샌드위치)      : {np.round(se_true, 4)}")
print(f"OLS 공식이 겨누는 값        : {np.round(se_ols_target, 4)}   (E[s^2] = {Es2:.3f})")
print(f"비 (OLS 겨냥 / 참)          : {np.round(se_ols_target / se_true, 4)}")

# (2) 4000 번 반복해 실제 표준편차·표준오차·포함률을 잰다.
R = 4000
beta = np.array([3.0, 2.0])
rng2 = np.random.default_rng(7)
b = np.empty((R, 2))
se_o = np.empty((R, 2))
se_r = np.empty((R, 2))
for r in range(R):
    Y = beta[0] + beta[1] * X + rng2.normal(0, 1, n) * X
    m = sm.OLS(Y, Xc).fit()
    b[r] = m.params
    se_o[r] = m.bse
    se_r[r] = m.get_robustcov_results(cov_type="HC3").bse

print(f"\n반복 {R} 회 (X 고정, 오차만 다시 뽑음)")
print(f"계수의 평균                 : {np.round(b.mean(axis=0), 4)}   참값 [3. 2.]")
print(f"계수의 평균 - 참값          : {np.round(b.mean(axis=0) - beta, 4)}"
      f"   (몬테카를로 표준오차 {np.round(b.std(axis=0, ddof=1) / np.sqrt(R), 4)})")
print(f"계수의 실제 표준편차        : {np.round(b.std(axis=0, ddof=1), 4)}")
print(f"  참 표준오차 (1)           : {np.round(se_true, 4)}")
print(f"OLS 표준오차의 평균         : {np.round(se_o.mean(axis=0), 4)}")
print(f"HC3 표준오차의 평균         : {np.round(se_r.mean(axis=0), 4)}")

t = stats.t(n - 2).ppf(0.975)
for name, se in [("OLS", se_o), ("HC3", se_r)]:
    cov = np.mean(np.abs(b - beta) <= t * se, axis=0)
    print(f"명목 95% 구간의 포함률 ({name})  : {np.round(cov, 4)}")

출력:

참 표준오차 (샌드위치)      : [1.0367 0.2543]
OLS 공식이 겨누는 값        : [1.4257 0.2412]   (E[s^2] = 34.913)
비 (OLS 겨냥 / 참)          : [1.3752 0.9481]

반복 4000 회 (X 고정, 오차만 다시 뽑음)
계수의 평균                 : [2.9915 2.0019]   참값 [3. 2.]
계수의 평균 - 참값          : [-0.0085  0.0019]   (몬테카를로 표준오차 [0.0165 0.0041])
계수의 실제 표준편차        : [1.0404 0.2565]
  참 표준오차 (1)           : [1.0367 0.2543]
OLS 표준오차의 평균         : [1.4176 0.2398]
HC3 표준오차의 평균         : [1.0434 0.2547]
명목 95% 구간의 포함률 (OLS)  : [0.9922 0.937 ]
명목 95% 구간의 포함률 (HC3)  : [0.952  0.9452]

① 계수는 불편이다. 평균이 \((2.9915,\ 2.0019)\)로 참값 \((3, 2)\)에서 \((-0.0085,\ 0.0019)\)만큼 떨어져 있고, 몬테카를로 표준오차가 \((0.0165,\ 0.0041)\)이다. 두 편차 모두 한 표준오차 안이라 \(0\)과 구별되지 않는다. 이분산이 뚜렷한 자료인데도 그렇다. (1)의 한 줄 유도가 말한 대로다.

② 참값을 맞히는 것은 HC3 뿐이다. 반복추출로 잰 \(\hat\beta\)의 실제 표준편차가 \((1.0404,\ 0.2565)\)이고 (1)이 설계행렬만으로 계산한 참 표준오차가 \((1.0367,\ 0.2543)\)이다. 모수 추정 없이 적은 식이 모의실험과 소수 셋째 자리까지 맞는다. 그리고

HC3 는 두 계수 모두 맞히고 OLS 는 두 계수 모두 틀린다. (1)에서 예측한 OLS 의 겨냥값 \((1.4257,\ 0.2412)\)도 관측된 평균 \((1.4176,\ 0.2398)\)과 맞는다(제곱근을 먼저 취해 평균하므로 Jensen 부등식으로 조금 작게 나온다).

③ 어긋나는 방향이 서로 반대다. 비 \((1.3752,\ 0.9481)\)이 그 말이다. OLS 는 절편의 표준오차를 \(38\%\) 과대, 기울기의 표준오차를 \(5\%\) 과소평가한다. 그 결과 포함률이

• OLS: 절편 \(0.9922\)(과대포함, 구간이 너무 넓다), 기울기 \(0.9370\)(과소포함, 구간이 너무 좁다)
• HC3: 절편 \(0.9520\), 기울기 \(0.9452\) — 둘 다 명목 \(0.95\)

가 된다. 같은 자료, 같은 적합인데 한 계수는 지나치게 보수적이고 다른 계수는 지나치게 자신만만하다. (1)에서 본 대로 OLS 에는 고를 수 있는 수가 \(E[s^2] = 34.913\) 하나뿐이므로 두 계수를 동시에 맞힐 길이 없다. 아래 상자와 연습문제 2 가 그 방향이 왜 계수마다 다른지를 지렛대로 설명한다.

그러므로 "로버스트 표준오차는 보수적"이 아니라 "로버스트 표준오차는 옳다"가 맞는 말이다. 기울기에서는 커지고 절편에서는 작아졌는데, 두 변화 모두 참값 쪽으로의 이동이었다.

일반적 결과. 보조회귀 \(Z_i = \gamma_0 + \boldsymbol{\gamma}'\mathbf{X}_i + u_i\)에서 \(H_0: \boldsymbol{\gamma} = \mathbf{0}\)을 검정할 때, \(n\)개 관측값에 대한 라그랑주 승수(스코어) 검정통계량은

이다. 여기서 \(R^2\)은 보조회귀의 결정계수이고 \(p\)는 검정하는 계수의 개수이다.

직관. \(R^2\)은 \(Z_i\)의 변동 중 \(\mathbf{X}_i\)가 설명하는 비율이다. \(H_0\)이 참이면 \(\mathbf{X}_i\)가 아무것도 설명하지 못하므로 \(R^2 \approx 0\)이고, \(nR^2\)이 작다. 관계가 있으면 \(R^2\)이 커지고 \(nR^2\)도 커진다.

\(n\)을 곱하는 이유는 \(R^2\) 자체가 표본크기와 무관한 비율이기 때문이다. 같은 \(R^2 = 0.1\)이라도 \(n = 20\)에서는 우연일 수 있지만 \(n = 1000\)에서는 확실한 신호이다. \(nR^2\)이 그 차이를 반영한다.

다른 LM 검정과의 연결. 같은 \(nR^2\) 형태가 여러 진단검정에 등장한다.

모두 "잔차(또는 그 제곱)에 아직 설명되지 않은 구조가 남아 있는가"를 묻는 같은 질문의 변형이다. \(\square\)

관찰. 절편은 \(1.501 \to 1.066\)(29% 감소), 기울기는 \(0.254 \to 0.272\)(7% 증가)이다.

원인. 샌드위치 추정량

은 각 관측값의 실제 잔차 크기 \(e_i^2\)로 가중한다. OLS는 이를 공통 \(s^2\)으로 대체한다.

두 값의 차이는 \(e_i^2\)과 관측값의 지렛대(leverage) 사이의 상관으로 결정된다.

• \(e_i^2\)이 큰 관측값이 그 계수에 대해 지렛대도 크면 → 로버스트 표준오차가 커진다.
• \(e_i^2\)이 큰 관측값이 그 계수에 대해 지렛대가 작으면 → 로버스트 표준오차가 작아진다.

이 예에서. \(\operatorname{Var}(\varepsilon_i) \propto X_i^2\)이므로 \(X\)가 큰 관측값의 잔차가 크다.

• 기울기에 대해서는 \(X\)가 극단적인(작거나 큰) 관측값의 지렛대가 크다. \(X\)가 큰 쪽에서 잔차도 크므로 양의 상관이 생겨 로버스트 표준오차가 커진다.
• 절편에 대해서는 \(X\)가 작은 관측값의 영향력이 크다(외삽 거리가 짧으므로). 그런데 그 관측값들의 잔차는 작다. 음의 상관이 생겨 로버스트 표준오차가 작아진다.

일반화. 로버스트 표준오차는 "OLS보다 보수적인 값"이 아니라 올바른 값이다. OLS 표준오차가 참값보다 클 수도 작을 수도 있으며, 로버스트 추정량은 어느 방향이든 바로잡는다.

실무적 함의: 이분산이 의심되면 로버스트 표준오차를 쓰되, 그것이 자동으로 "더 안전한" 결과를 준다고 기대하지 말라. 유의성이 오히려 강해질 수 있다. \(\square\)

자유도 계산. White 검정의 보조회귀에 들어가는 항(상수 제외)은

• 원래 설명변수: \(p\)개
• 제곱항 \(X_j^2\): \(p\)개
• 교차항 \(X_j X_l\) (\(j < l\)): \(\binom{p}{2} = \frac{p(p-1)}{2}\)개

합계는

(이진 설명변수가 있으면 \(X_j^2 = X_j\)이므로 그만큼 항이 줄어들고, statsmodels는 이런 완전공선 항을 자동으로 제거한다. 본문 보기에서 \(p=1\)인데 White 자유도가 2인 것이 표와 일치한다.)

검정력 문제. 카이제곱 검정에서 자유도가 커지면 임계값이 커진다.

이분산의 신호가 소수의 항에만 집중되어 있다면, 나머지 무의미한 항들이 \(R^2\)을 거의 올리지 못하면서 임계값만 높인다. 곧 신호를 잡음으로 희석한다.

구체적으로 \(p = 10\)이면 White 검정은 65개 항을 쓴다. \(n = 200\)인 자료에서 65개 설명변수의 보조회귀는 그 자체로 과적합 위험이 있고, 검정력이 Breusch-Pagan(자유도 10)보다 크게 떨어진다.

실무 지침.

• \(p\)가 작고(3 이하) 이분산의 형태를 모르면 White 검정.
• \(p\)가 크거나 이분산이 특정 변수와 관련될 것으로 예상되면 Breusch-Pagan, 또는 의심되는 변수만 넣은 보조회귀.
• 어느 쪽이든 탐지 여부와 무관하게 로버스트 표준오차를 기본으로 쓰는 것이 가장 간단한 해법이다. 이분산이 없으면 로버스트 표준오차가 OLS와 거의 같으므로 잃는 것이 거의 없다. \(\square\)

왜 제거하는가. 검정의 목적은 \(X\)가 작은 구간과 큰 구간에서 오차분산이 다른지 보는 것이다. 자료를 정렬하여 절반씩 나누면 두 집단의 경계 근처 관측값들은 \(X\) 값이 비슷하므로 분산도 비슷하다. 이 관측값들이 두 집단의 대비를 흐린다.

가운데를 잘라내면 두 집단의 \(X\) 범위가 확실히 분리되어 분산 차이가 뚜렷해진다. 검정력이 올라간다.

\(c = n/5\)의 근거. Goldfeld와 Quandt의 원논문은 모의실험으로 \(c \approx n/3\)까지 검토했고, 이후 문헌에서 \(n/5\) 정도가 대비 강화와 표본 손실 사이의 절충으로 자리 잡았다. 너무 많이 자르면 각 집단의 자유도가 줄어 오히려 검정력이 떨어진다.

한계.

1. 정렬 변수를 미리 알아야 한다. 이분산이 어느 변수와 관련되는지 모르면 쓸 수 없다. 여러 변수를 시도하면 다중검정 문제가 생긴다.
2. 단조 이분산만 탐지한다. 나비넥타이 모양(가운데가 좁고 양끝이 넓은)처럼 비단조 패턴은 놓친다. 오히려 가운데를 잘라내는 것이 역효과를 낸다.
3. F 검정의 정규성 민감도를 물려받는다. 4단계에서 두 잔차제곱합의 비를 F 분포와 비교하는데, 15.3절에서 보았듯 F 검정은 비정규성에 극도로 민감하다. 오차가 두꺼운 꼬리를 가지면 크기가 심하게 부풀려진다.
4. 관측값을 버린다. \(n/5\)를 버리는 것은 정보의 낭비이다. Breusch-Pagan은 모든 관측값을 쓴다.
5. 다변량으로 확장되지 않는다. 설명변수가 여러 개일 때 하나만 골라 정렬해야 한다.

결론. Goldfeld-Quandt 검정은 역사적 의의와 교육적 직관은 있지만, 오늘날 실무에서는 Breusch-Pagan(Koenker 판)이나 White 검정이 거의 모든 면에서 낫다. 이분산의 원인 변수가 명확하고 관계가 단조이며 오차가 정규에 가깝다고 확신하는 특수한 경우에만 고려할 만하다. \(\square\)