Exploring Regression Models for Regression Analysis (1): Regression Analysis, Linear Regression, HC Standard Errors

들어가며

차라투 블로그 “Exploring Regression Models for Regression Analysis”에서는 세(네) 장에 걸쳐서 통계에서의 Regression Analysis를 위한 여러 Regression Model들을 수학적으로 깊게 탐구합니다. 1장에서는 Regression Analysis의 개념, (simple, multiple or general) linear regression에 대한 개념 및 Analysis 과정에 대한 소개로 시작해서 data의 특성 중 하나인 Homo, Heteroskedasticity의 개념, Heteroskedasticity로 추정되는 경우 linear regression 모델의 robust한 (co)variance를 구하는 방법 중 하나이자 필수로 고려해야 하는 Heteroskedasticity-consistent standard errors (heteroskedasticity-robust standard errors)와 Wild Bootstrap를 살펴볼 것입니다. 2장에서는 1장에서 다룬 기본 linear regression에서 link function을 도입하여 regression의 개념을 outcome of single yes/no, outcome of single K-way, count 등 non-normal한 종속변수로 확장한 Generalized linear model의 개념을 Exponential Family, Link Function와 같은 핵심 개념과 함께 깊게 살펴보며, model의 parameter를 estimate하는 알고리즘인 IRLS(Fisher scoring), HC standard errors의 clustered data 버전인 Cluster-robust standard error를 다룹니다. 3장에서는 GLM에서 여전히 남아있는 observations의 independent 조건(오차항의 독립성)을 극복하여 clustered(panel, longitudinal..) data를 고려할 수 있는 모델인 GEE, GLMM의 근본적인 원리를 M-estimator, robust (sandwich) estimator와 같은 수학적 개념들과 함께 탐구할 것입니다. {4장은 계획 중에 있습니다.}

결국 “Exploring Regression Models for Regression Analysis”는 Regression Analysis을 논리적으로 수행하기 위해 어떤 Regression Model 들이 있으며, 이들의 수학적인 원리가 무엇인지를 deep dive 하는 글입니다. 영문 글들 중에서도 이 내용들을 이어지게 다뤄주는 글은 거의 없고, 아주 친절하게 수식을 전개하였기 때문에 근본적인 Regression Models의 수학적 원리를 공부하기 좋은 글들이라고 생각하며, R 코드는 간단하게만 주석 형태로 제공하기 때문에 로컬에서 주석 해제 후 돌려보시고, 실제 의학 연구에 관련 분석이 필요하시다면 해당 연구를 차라투로 문의해주시길 바랍니다.

1. Regression Analysis

Regression Analysis는 통계학에서 가장 널리 사용되는 방법론 중 하나로, 어떤 종속변수(dependent variable)와 하나 이상의 독립변수(independent variables) 간의 관계를 추정하고 해석하는 분석 기법입니다. (종속변수 또한 벡터, 즉 여러 개일 수 있습니다.)

의학 분야에서는 예후를 예측하는 모델을 만들거나(예: 특정 약물 투여 후 혈압, outcome 등의 변화를 예측), 특정 위험인자(risk factor)가 결과에 유의미한 영향을 미치는지(예: 어떤 치료법이 환자의 생존율에 유의미한 차이를 주는지) 등을 살펴보기 위해 Regression Analysis가 필수적으로 활용됩니다. 이러한 Regression Analysis을 완벽하게 수행하기 위해서는 본인의 data가 어떤 특성을 갖고 있는지 파악하고, 이에 대해 어떠한 Regression Model을 고려해야 하는지 알고, 이 모델이 어떻게 data에 적합(fit)하는지에 대한 대략적인 수학 또는 알고리즘의 원리와, 적합된 모델이 통계적으로 유의미한지 검정(test)하고 예측 성능이나 해석력을 평가하는 방법까지 전 과정을 이해해야 합니다. 본 1장에서는 가장 간단한, 대신 많은 가정이 필요한 Regression Model인 (General) Linear Regression에 대해서 위 과정과 함께 살펴볼 것이며, 이후의 장들에서는 data의 가정이 조금씩 깨질 때 (data의 특성이 변할 때) 어떤 Regression Model을 고려해야 하는지와 이 Model들의 parameters는 어떻게 수학적으로 추정하는지 살펴볼 것입니다.

2. Linear Regression (LM)

2.1. (Simple, Multiple) Linear Regression 정의

---

(General) Linear Regression은 의학뿐만 아니라 다양한 분야에서 가장 basic한 Regression Model입니다. Linear Regression은 독립변수 $X$ 와 종속변수 $Y$ 간의 선형(linear) 관계를 가정하고, 이를 통해 종속변수를 설명하거나 예측합니다. 독립변수가 1개일 때는 Simple Linear Regression, 2개 이상일 때는 Multiple Linear Regression라고 부릅니다. 의학 연구에서 예를 들면, 혈압(종속변수)을 나이, 체중, 성별(독립변수) 등으로 설명하거나 예측하는 과정을 생각할 수 있습니다. 이러한 Linear regression은 아래 네 가지 가정을 기반으로 추정하며, 이 가정들에 대해서 잘 이해하는 것은 매우 중요하고, 이들 중 특정 가정들이 위배될 경우 이후에 다룰 Regression Model들을 고려해야 함을 명심해주시면 좋을 것 같습니다.

• 선형성(Linearity) 가정

  독립변수와 종속변수가 선형 관계에 있다고 가정합니다. 산점도(scatter plot)를 통해 대략적인 선형 관계 여부를 확인할 수 있습니다.

• 오차항의 정규성(Normality) 가정

  주어진 독립변수의 값에서 종속변수의 확률분포는 정규분포를 따른다고 가정합니다. 이는 오차항 ε이 정규분포를 따른다고도 표현할 수 있습니다.(종속변수가 정규분포를 따른다는 뜻은 Linear Model이 종속변수의 mean을 예측하므로, 이 둘의 차인 오차항은 평균이 0이고 분산은 종속변수와 같은 정규분포를 따르게 되기 때문입니다.) 이는 정규 P-P plot 혹은 Q-Q plot 등을 통해 가정 위배 여부를 대략적으로 확인할 수 있습니다.

• 오차항의 독립성(Independence) 가정

  각 관측치(Observations, Data(set)) 또는 잔차(residual)가 서로 독립이라고 가정합니다. 즉, data간의 상관관계가 없다고 가정하는 것이고, 잔차산점도(residual plot), Durbin-Watson 통계량 등을 통해 자기상관(autocorrelation)이 있는지 살펴볼 수 있습니다.

• 오차항의 등분산성(Homoscedasticity) 가정

  모든 독립변수의 값에서 종속변수의 분산이 동일하다고 가정합니다. 잔차산점도 등을 통해 잔차가 일정한 분산을 가지는지 대략적으로 확인할 수 있습니다.

이후 Linear Regression Model은 다음과 같은 과정으로 분석을 수행하게 됩니다; (1) 최소제곱법(Least Squares)를 통해 model parameters를 추정(estimate), (2) 결정계수($R^2$)를 통해 모델이 종속변수를 얼마나 잘 설명하는지 확인, (3) F-검정(F-test)으로 전체 회귀식의 유의성을 검정, (4) 검정(t-test)으로 각 회귀계수(regression coefficient)가 유의미한지 확인. 특히, 독립변수 각각의 Model coefficient(parameter)를 검정 함으로써 종속변수와의 상관성을 분석하는 (4) 과정은 유의성을 판단하는 가장 중요한 검정 과정으로, 고전적으로 Wald Test, Likelihood Ratio Test, Score Test가 자주 사용됩니다.

2.2. Linear Regression 수학적 표현 및 추정

---

Linear Regression은 다음과 같은 형태를 가집니다:

$$Y_i={\beta }_0+{\beta }_1{X}_{i1}+{\beta }_2{X}_{i2}+\cdots +{\beta }_k{X}_{ik}+{\epsilon }_i$$

혹은 Matrix 형태로 간단히 쓰면,

$$\mathbf{y}=\mathbf{X}\mathit{\beta }+\mathit{\epsilon }$$

$$where,\mathbf{y}\in {\mathbb{R}}^n,\mathbf{X}\in {\mathbb{R}}^{n\times p},\mathit{\beta }\in {\mathbb{R}}^p,\mathit{\epsilon }\in {\mathbb{R}}^n$$

• $\mathbf{y}\in {\mathbb{R}}^n$: 관측된 종속변수들의 벡터

• $\mathbf{X}\in {\mathbb{R}}^{n\times p}$: 설계 행렬, 각 행은 하나의 관측치(Observation)이며 ($n$), 각 열은 하나의 독립변수에 ($p$) 해당

• $\mathit{\beta }\in {\mathbb{R}}^p$: 추정하고자 하는 회귀 계수(or 모수 or parameters) 벡터 (여기선 intercept ${\beta }_0$ 제외)

• $\mathit{\epsilon }\in {\mathbb{R}}^n$: 랜덤 오차 항 벡터, 위 네 가정에 따라 랜덤 오차항은 평균이 $E(\mathit{\epsilon })=0$이고, 분산이 모든 관측치에서 같으며 독립인 정규분포

입니다. Multiple Linear Regression의 모델 $\mathit{\beta }$의 분산은 scalar variance가 아니라 (co)variance matrix 형태가 될 것입니다. (simple도 intercept를 고려하면 matrix) 앞으로 세 장에 걸쳐 이를 분산 행렬이라고 부르겠지만, 정확히는 분산-공분산 행렬이라고도 부릅니다. 또한, 위 matrix form 수식처럼 intercept를 제외하는 경우도 많은데, 항상 있다고 가정하며 가독성을 위해 제외한 것이니 혼동할 필요 없이 자연스럽게 읽으시면 됩니다. Covariance matrix에서는 회귀계수 서로 간의 공분산이 들어있으므로, 각 회귀계수(독립변수의 유의성)에 대한 검정은 Covariance Matrix에 있는 각 diagonal 원소들을 사용하여 수행하게 됩니다. (회귀계수 각각에서 스스로의 대한 공분산 = 회귀계수의 분산) 예를 들어 ${\beta }_2$의 분산은 $p$가 5 (독립변수가 4개, intercept 포함)일 때 Covariance Matrix의 3행 3열 값이 될 것입니다. (intercept가 없다면 2행 2열) 여기서 확인할 수 있듯 Regression Model의 parameter 추정 역시 중요하지만, parameter의 Covariance Matrix를 추정하는 것 또한 유의성 판단에서 매우 중요하며, 이후의 장들 또한 자연스럽게 Model의 소개, paramater 추정법, parameter의 (Co)variance Matrix를 추정법으로 이루어지게 될 것입니다.

또한, Linear Model 식은 선형대수학에서 Hyper Plane으로 정의하는 형태가 되는데, 쉽게 설명하자면 독립변수가 하나일 경우 종속변수와의 2차원 평면에서 직선(2차원의 Hyper Plane)을 띄고, 독립변수가 두 개일 경우 종속변수와의 3차원 공간에서 평면을 띄는(3차원의 Hyper Plane), 선형적으로 공간을 두 부분으로 가르는 공간이라고 생각하시면 될 것 같습니다. (마찬가지로 쉽게 생각하면, 비선형적일 경우 직선이 아니라 곡선, 평면이 아니라 곡면일 것입니다.)

이제 어떻게 Linear Model을 regression, 즉 fit(parameter를 추정)할 것인지 설명하겠습니다. 가장 기본적으로 알려진 Model의 통계학적 parameter 추정법으로는 Ordinary Least Squares(OLS), Maximum Likelihood Estimation(MLE), Method of Moment(MOM)가 있습니다. 이후의 모델들은 대부분 MLE로 모델을 추정하지만, Linear Model은 특별하게도 MLE와 OLS의 추정 결과가 같고, OLS의 추정 결과는 BLUE(Best Linear Unbiased Estimator) 입니다. 최소제곱법(OLS, Ordinary Least Squares) 추정량(estimator)은 오차 제곱합(SSR, Sum of Squared Residuals), 즉 모든 데이터에서 실제 종속변수 값과 모델의 예측 값의 차이를 제곱한 수들의 합을 최소화하는 $\mathit{\beta }$를 찾는 것이며, 다음의 closed-form solution을 갖습니다: $$\hat{\mathit{\beta }}={\left({\mathbf{X}}^{\mathrm{\top }}\mathbf{X}\right)}^{-1}{\mathbf{X}}^{\mathrm{\top }}\mathbf{y}$$

그리고, 고전적 가정(특히 등분산성, 독립성, 정상성)이 모두 충족된다고 할 때, 이 $\hat{\mathit{\beta }}$의 추정오차의 공분산행렬(covariance matrix)은 $$\mathrm{Var}(\hat{\mathit{\beta }})={\sigma }^2{\left({\mathbf{X}}^{\mathrm{\top }}\mathbf{X}\right)}^{-1}$$ 이며, 실제 계산할 때는 데이터에서 추정한 ${\hat{\sigma }}^2$를 통해 $$\widehat{\mathrm{Var}}(\hat{\mathit{\beta }})={\hat{\sigma }}^2{\left({\mathbf{X}}^{\mathrm{\top }}\mathbf{X}\right)}^{-1},where{\hat{\sigma }}^2=\frac{\sum e_i^2}{N-K},e_i=y_i-x_i\hat{\mathit{\beta }}$$

로 추정합니다. (분모에서 K를 뺀 이유는 degree of freedom에 의한 것이며, 에러 term이 정확히는 차원에 맞춰 $x_i^{\mathrm{\top }}\hat{\mathit{\beta }}$ 지만, 앞으로 이정도 표기는 가독성을 위해 넘어가겠습니다.) 이렇게 구한 모델 추정량과 추정오차의 공분산행렬을 통해 각 회귀계수에 대한 검정을 수행하거나 신뢰구간(CI)을 계산할 수 있게 됩니다.

다음으로 넘어가기 전에, Linear Regression에서 (1) $\hat{\beta }$이 closed-form solution $({\mathbf{X}}^{\mathrm{\top }}\mathbf{X}{)}^{-1}{\mathbf{X}}^{\mathrm{\top }}\mathbf{y}$을 갖게 되는 과정과, (2) 위에서 언급하였듯 model의 parameter ($\hat{\beta }$) OLS 추정량이 MLE(Maximum likelihood estimation) 추정량과 같음을 증명하겠습니다.

(1) Prove $\hat{\beta }=({\mathbf{X}}^{\mathrm{\top }}\mathbf{X}{)}^{-1}{\mathbf{X}}^{\mathrm{\top }}\mathbf{y}$.

우선, OLS 추정량은 오차 제곱합(SSR)을 최소화하는 식이므로 오차 SSR을 표현하는 식 $J(\beta )$를 적으면 다음과 같습니다.

$$J(\beta )=\frac{1}{2}(X\beta -y{)}^{\mathrm{\top }}(X\beta -y)=\frac{1}{2}\sum _{i=1}^n(x_i\beta -y^{(i)}{)}^2$$

이제 위 $J(\beta )$를 $\beta$에 대해 미분하였을 때 0이 나오는 $\hat{\beta }$이 SSR이 최소가 되는 OLS 추정량입니다.(이에 대한 증명은 안하겠지만, 간단하게 볼록한 2차 함수에서는 미분값이 0인 점이 최소점인 것과 같은 원리라고 생각하시면 되겠습니다.) 식을 풀어보면,

$$\begin{gathered} {\mathrm{\nabla }}_{\beta }J(\beta )={\mathrm{\nabla }}_{\beta }\frac{1}{2}(X\beta -y{)}^{\mathrm{\top }}(X\beta -y) \\ =\frac{1}{2}{\mathrm{\nabla }}_{\beta }((X\beta {)}^{\mathrm{\top }}X\beta -(X\beta {)}^{\mathrm{\top }}y-y^{\mathrm{\top }}(X\beta )+y^{\mathrm{\top }}y) \end{gathered}$$ 입니다. $\mathrm{\nabla }$는 미분하는 변수가 scalar가 아니라 vector이기 때문에 사용되는 기호이며(gradient), 그냥 $\beta$로 미분한다고만 생각하시면 됩니다. 이후 미분하는 과정 또한 크게 다르지 않습니다. 식을 계속 진행해보면, $y^{\mathrm{\top }}y$는 식에 $\beta$가 없고 $(X\beta {)}^{\mathrm{\top }}y$와 $y^{\mathrm{\top }}(X\beta )$ 둘다 단순히 두 벡터의 내적인 똑같은 식이므로 $$=\frac{1}{2}{\mathrm{\nabla }}_{\beta }({\beta }^{\mathrm{\top }}(X^{\mathrm{\top }}X)\beta -2(X^{\mathrm{\top }}y{)}^{\mathrm{\top }}\beta )$$이고, 이제 $\beta$로 미분하면, 우항은 그냥 미분, 좌항은 두 번 나오므로 곱미분을 수행하여
$$=\frac{1}{2}(2X^{\mathrm{\top }}X\beta -2X^{\mathrm{\top }}y)$$
$$=X^{\mathrm{\top }}X\beta -X^{\mathrm{\top }}y=0$$

입니다. 결국 OLS 추정량 $\hat{\beta }$는

$$X^{\mathrm{\top }}X\hat{\beta }=X^{\mathrm{\top }}y$$ 이고, 양변 앞에 역행렬을 곱해주면
$$\hat{\beta }=(X^{\mathrm{\top }}X{)}^{-1}{X}^{\mathrm{\top }}y$$가 됩니다.

(2) Prove OLS estimator is same as MLE estimator in Linear Model (Regression).

위 Linear Regression Model 식을 다음과 같이 작성할 수 있습니다. $$\mathbf{y}=\mathbf{X}\mathit{\beta }+\mathit{\epsilon },\mathit{\epsilon }\sim \mathcal{N}(\mathbf{0},{\sigma }^2\mathbf{I})$$

식이 어색할 수 있지만 정확히 이전의 수학적 표현과 동일하며, $\mathit{\epsilon }$ term이 복잡해 보이는 이유는, 이전의 네 가지 조건들을 수식으로 반영하여 모델을 확률 기반으로 해석하였기 때문입니다. 식을 설명드리자면, 선형성(Linearity) 가정에 의해 error의 mean이 (기댓값이) 0이 되고, 오차항의 정규성 가정에 의해 정규 분포를 따르며, 오차항의 독립성 및 등분산성 가정에 의해 분산(Covariance Matrix)이 값이 모두 ${\sigma }^2$인 대각행렬, 즉 다른 observations간의 correlation은 0이고(없고), 각 observations의 분산은 ${\sigma }^2$로 일정함을 표현한 것입니다.

최대가능도방법 또는 최대우도법(Maximum Likelihood Method)은 어떤 모수와 표집한 값들이 주어졌을 때, 표집한 값들이 나올 가능도(확률)을 최대로 만드는 모수를 선택하는 아주 general한 모델의 점 추정 방법이며, 처음 들어보셨다면 대표님께서 더욱 자세하게 설명하신 블로그 글을 읽으시면 좋을 것 같습니다. 철학적으로는 종속변수의 확률 분포를 데이터에 맞게 가정한 후, 모든 observations(data) 각각이 나올 확률(분포)을 모두 곱한 식이 최대가 되도록 하는 모수값이 바로 MLE를 통해 추정한 parameter라고 생각하시면 됩니다. 간단한 예시를 들어보자면, 동전을 던져 앞면이 7번, 뒷면이 3번 나온 데이터에서 동전이 앞면이 나올 확률을 모수로 두고 종속변수의 분포를 베르누이 분포로 가정하면, 각각의 동전을 던져 얻은 관측치들 10개의 data의 Likelihood(확률)를 모두 곱하면 모수가 하나였으니 변수가 하나인 하나의 식(함수)가 나오고, 이 식(함수)을 최대화 하는 모수를 찾는 method가 MLE라고 간단하게 정리할 수 있을 것 같습니다. (이때 추정된 결과는 0.7일 것이고, 모수로 식을 미분한 함수(score function)이 0인 값을 찾음으로써 모수를 추정하게 되며, MLE 이외에도 MAP, 베이지안 등 다른 추정 기법들도 있지만 여기에서는 MLE 정도를 이해하면 충분할 것 같습니다.)

이제 돌아와서 LM에서는 OLS estimator와 MLE estimator가 상응함을 보이겠습니다. 이후의 식 전개는 dim이 1일 때로 증명하겠습니다. Multi-dimension에서는 분포식과 term이 더 복잡하지만 meaning은 1차원과 정확히 동일하기 때문입니다. 오차 항 $\mathit{\epsilon }$는 정규분포를 따르므로, n개의 관측치에 대한 Likelihood는 $$\mathcal{L}(\mathit{\beta })=\prod _{i=1}^n\frac{1}{\sqrt{2\pi {\sigma }^2}}\exp (-\frac{(y_i-{\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta }{)}^2}{2{\sigma }^2})$$

이며, 로그를 씌운 log Likelihood는 $$\ell (\mathit{\beta })=\log \mathcal{L}(\mathit{\beta })$$

입니다. 식이 복잡해 보이지만 데이터 n개에 대하여 정규(가우시안)분포를 따르는 종속변수가 모델의 예측값을 가질 확률을 단순히 모두 곱한 식입니다. 결국 얻어진 표본 data 자체가 나올 확률에 대한 식을 세운 후, 이를 최대화 하는 parameter를 찾는 과정이 MLE estimation이라고 직관적으로 해석할 수 있습니다.

돌아와서, log를 씌운 이유를 생각해보겠습니다. log의 그래프를 생각해보시면 log함수는 직관적으로도 monotonically increase하기 때문에 likelihood를 최대화 하는 parameter와 log likelihood를 최대화 하는 parameter는 같으며, log는 곱셈을 덧셈으로, 지수 term을 아래로 만들어주어 likelihood 식이 매우 간단해지기 때문에 (아래에서 실제로 그런지 확인해보세요.) MLE에서는 항상 log likelihood를 최대화하는 parameter를 찾는 방향으로 parameter를 추정합니다. 이어서 전개하면,

$$=\log \prod _{i=1}^n\frac{1}{\sqrt{2\pi {\sigma }^2}}\exp (-\frac{(y_i-{\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta }{)}^2}{2{\sigma }^2})$$
$$=\sum _{i=1}^n\log \frac{1}{\sqrt{2\pi {\sigma }^2}}-\frac{(y_i-{\mathbf{x}}_{\mathbf{i}}^{\mathrm{\top }}\mathit{\beta }{)}^2}{2{\sigma }^2}$$
$$=n\log \frac{1}{\sqrt{2\pi {\sigma }^2}}-\frac{1}{2{\sigma }^2}\sum _{i=1}^n(y_i-{\mathbf{x}}_{\mathbf{i}}^{\mathrm{\top }}\mathit{\beta }{)}^2$$

입니다. Log likelihood를 최대화 하는 것은 negative Log likelihood를 최소화 하는 것과 같으므로 양변에 -를 곱하면
$$-\ell (\mathit{\beta })=-n\log \frac{1}{\sqrt{2\pi {\sigma }^2}}+\frac{1}{2{\sigma }^2}\sum _{i=1}^n(y^{(i)}-{\mathbf{x}}^{(i{)}^{\mathrm{\top }}}\mathit{\beta }{)}^2$$

이고, 상수 $-n\log \frac{1}{\sqrt{2\pi {\sigma }^2}}$는 $\mathit{\beta }$에 영향을 주지 않으므로 제거하면, $$-\log \mathcal{L}(\mathit{\beta })\approx \frac{1}{2{\sigma }^2}\sum _{i=1}^n(y^{(i)}-{\mathbf{x}}^{(i{)}^{\mathrm{\top }}}\mathit{\beta }{)}^2$$

가 됩니다. 아까 보았듯 OLS estimator는 SSR을 최소화하는 아래 문제로 정의되므로, $${\hat{\mathit{\beta }}}_{\text{OLS}}=\arg \underset{\mathit{\beta }}{min}\sum _{i=1}^n(y^{(i)}-{\mathbf{x}}^{(i{)}^{\mathrm{\top }}}\mathit{\beta }{)}^2$$

결국 위 두 식, OLS estimator와 MLE estimator의 목적식이 같음을 알 수 있고(종속변수와 모델 예측의 차이의 제곱을 모든 data point에서 더한 값을 최소화 하는 beta), 따라서 추정량 또한 같습니다. $${\hat{\mathit{\beta }}}_{\text{MLE}}={\hat{\mathit{\beta }}}_{\text{OLS}}$$

위 증명은 MLE로 LM의 parameter를 추정한 것이 아니라 식 전개를 통해 OLS와 목적식이 같음을 보인 것이고, 실제 MLE로 parameter를 추정하는 과정은 2장부터 자세히 보게 될 것입니다. 한 가지 첨언하자면, Linear model에서는 이 두 추정량이 같지만, 정규성 가정이 깨진 경우 두 추정량은 같지 않을 수 있습니다. 이유를 생각해보면, 정규분포는 확률이 가장 큰, 즉 mode가 mean이며 mean을 기준으로 양 옆이 symmetric한 분포이기 때문에 두 추정량이 같으며, 정규분포가 아닐 경우 mean을 추정하려고 하는 MLE는 절대적인 거리를 좁히려는 OLS와 다른 값을 추정하게 된다고 직관적으로 생각해볼 수 있습니다.

3. Linear Regression in Homo vs. Heteroskedasticity

앞서 언급한 Linear Regression 모형은 “오차항의 등분산성(homoscedasticity)”을 가정합니다. 그러나 실제 연구 상황에서는 이 가정이 위배되는 경우가 흔합니다. Heteroskedasticity는 이러한 오차항은 독립 변수에 따라 같지 않을 수 있다는 뜻으로, 등분산성 가정이 깨진 경우를 의미합니다. 즉,

• Homoscedasticity: 모든 관측값에 대해 오차항(또는 종속변수)의 분산이 동일
• Heteroskedasticity: 오차항(또는 종속변수)의 분산이 관측값에 따라 다름

입니다. 즉 이전 수식 중 오차 항을 $\mathit{\epsilon }$로 두었다면, observations의 분산은 더이상 ${\sigma }^2$로 일정하지 않고 observations마다 다른 값을 갖을 수 있다는 뜻입니다. 예를 들어, 임상 연구에서 나이(age)에 따라 혈압(blood pressure)의 분산이 달라질 수 있습니다. 이런 상황에선 등분산성을 가정하는 것이 부적절해집니다. 등분산성 가정이 깨진 상황에서 Linear regression을 수행하면, $\hat{\mathit{\beta }}$ 자체가 편향되진 않더라도(즉, 일관성(consistency)은 여전히 유지), $\hat{\mathit{\beta }}$ 의 분산 추정치 ${\hat{\sigma }}^2({\mathbf{X}}^{\mathrm{\top }}\mathbf{X}{)}^{-1}$ 가 bias를 갖게 되어 잘못된 표준오차, 잘못된 유의성 검정 결과로 이어질 수 있습니다. 따라서, Heteroskedasticity가 의심되는 상황에서는 모델의 분산을 안정적으로 추정해야 하며, 이때 1번으로 고려되는 방법 중 하나가 Heteroskedasticity-consistent(heteroskedasticity-robust) standard errors, 줄여서 HC standard errors를 통한 (Co)variance Matrix estimation입니다.

4. Heteroskedasticity-Consistent Standard Errors (HC Standard Errors)

4.1. HC Standard Errors 정의

---

HC(Heteroskedasticity-Consistent) standard errors는, 고전적 선형회귀 가정 중 등분산성만 깨졌을 때(나머지 선형성, 정규성, 독립성 가정 in LM 은 그대로 유지) 회귀계수 추정치의 분산을 “강건(robust)”하게 추정하기 위한 method입니다. 이는 heteroskedasticity-robust standard errors 또는 HCCME(Heteroskedasticity-Consistent Covariance Matrix Estimation) 라고도 부르며, 이전에로 구한 $\hat{\beta }$를 그대로 사용하되, 그 공분산행렬 추정만 새롭게(robust하게) 구하는 방식입니다. 다만, 독립성(오차항들이 서로 독립), 선형성, 정규성 등의 다른 가정이 또 깨져 있다면 HC SE만으로는 대응할 수 없음을 명심하시길 바랍니다. 예를 들어, 오차항이 서로 상관되어 있는 clustered data에서는 더 이상 추정량이 consist하지 않기 때문에, 이보다 한 단계 더 나아간 cluster-robust standard errors, 혹은 GEE, GLMM 등의 모델을 고려해야 합니다. 즉, “등분산성 가정만 깨졌을 때” 쓰는 표준오차(or covariance matrix) 추정량이라고 생각하면 됩니다. 또한, 모델이 실제로 크게 잘못 설정되어 있다면(선형성조차 안 맞는 경우, 모델의 estimator가 크게 편향, HC se가 모델에서 구한 se와 크게 차이가 나는 경우 등), 그때는 variance(standard error)만 “robust standard errors”로 바꾼다고 해서 문제가 해결되지 않고, 모델을 수정하거나 data를 다시 확인하고 분석을 다른 방식으로 수행해야 합니다.

즉 Greene의 말처럼,

“Simply computing a robust covariance matrix for an otherwise inconsistent estimator does not give it redemption.”

이라는 점을 늘 유의해야 합니다. 정리하자면, 모델을 잘 구성하고, 모델의 분산과 큰 차이가 없을 경우에는 모델의 parameter의 분산 효율성이 최대가 아니거나 MLE 추정을 통해 얻은 모델의 분산이 더 이상 불편향이 아닐 때도 heteroskedasticity-robust standard errors를 통해서 모델의 분산을 robust하게 estimate할 수 있으며, 분석 결과 편향이 너무 크지 않다면 이후의 검정에서 이 HC se를 사용함으로써 이후의 통계 분석에서 설득력을 높일 수 있습니다. 추가로 스포하자면, 아래에서 다룰 Heteroskedasticity-robust Standard Errors 식은 Linear Model에서 robust하게 분산을 추정하는 Robust(or Sandwich) Estimation이며, 이후의 모델에서도 해당 모델에 대한 Sandwich Estimator를 고려함으로써 안정적인 분산 추정 방법에 대해 다룰 것입니다.

4.2. HC(CME)0 수학적 표현

---

이제 위에서 설명한 Heteroskedasticity한 data에서 robust하게 standard errors, 혹은 Covariance Matrix를 estimate 하는 heteroskedasticity-consistent standard errors의 수식을 유도하겠습니다.

이전에 선형 회귀 모델을 다음과 같이 정의한 적이 있습니다.
$$\mathbf{y}=\mathbf{X}\mathit{\beta }+\mathit{\epsilon }$$ 이때, 위에서는 등분산성을 가정하였기 때문에 $\mathit{\epsilon }$를 원소값이 모두 ${\sigma }^2$ 인 대각행렬로 가정하여 각 observations의 분산은 ${\sigma }^2$로 일정함을 가정했다면, 이번에 정의한 $\mathit{\epsilon }$은 $E(\mathit{\epsilon })=0$, $E(\mathit{\epsilon }{\mathit{\epsilon }}^{\mathrm{\top }})=\mathrm{\Phi }$ 입니다. 이때 $\mathrm{\Phi }$는 $\mathrm{diag}({\epsilon }_i^2)$로 정의되는 대각행렬이며(대각 성분을 제외하면 모두 0인 행렬), 식을 통해서 여전히 error term의 mean이 0이고 observations들은 independent하지만(대각행렬이므로), 더 이상 observations의 분산이 서로 일치하지 않음을 표현한 것입니다. 즉, 이제 Homoscedasticity에서 Heteroskedasticity을 고려하기 시작했다는 것을 수식을 해석함으로써 알 수 있습니다.

Variance estimatation

위 식의 OLS 추정량은 여전히 $$\hat{\mathit{\beta }}=({\mathbf{X}}^{\mathrm{\top }}\mathbf{X}{)}^{-1}{\mathbf{X}}^{\mathrm{\top }}\mathbf{y}$$

입니다. (still consistent) 이제 모델의 variance 식을 전개해보면, $$\mathrm{Var}(\hat{\mathit{\beta }})=\mathrm{Var}(({\mathbf{X}}^{\mathrm{\top }}\mathbf{X}{)}^{-1}{\mathbf{X}}^{\mathrm{\top }}\mathbf{y})$$

이며, $({\mathbf{X}}^{\mathrm{\top }}\mathbf{X}{)}^{-1}{\mathbf{X}}^{\mathrm{\top }}$는 determinant한 값으로, variance는 식 내의 determinant한 salar는 Var 밖으로 제곱과 함께 나오는 것처럼 ($\mathrm{Var}(aX)=a^2\mathrm{Var}(X)$), 행렬 또는 벡터에 대해도 비슷하게 $$\begin{gathered} =({\mathbf{X}}^{\mathrm{\top }}\mathbf{X}{)}^{-1}{\mathbf{X}}^{\mathrm{\top }}\mathrm{Var}(\mathbf{y}){(({\mathbf{X}}^{\mathrm{\top }}\mathbf{X}{)}^{-1}{\mathbf{X}}^{\mathrm{\top }})}^{\mathrm{\top }} \\ =({\mathbf{X}}^{\mathrm{\top }}\mathbf{X}{)}^{-1}{\mathbf{X}}^{\mathrm{\top }}\mathrm{Var}(\mathbf{y})\mathbf{X}({\mathbf{X}}^{\mathrm{\top }}\mathbf{X}{)}^{-1} \end{gathered}$$ 가 됩니다. (벡터와 행렬을 다룰 때에는 위처럼 제곱한다고 생각하시면 될 것 같습니다.) 참고로, $${(({\mathbf{X}}^{\mathrm{\top }}\mathbf{X}{)}^{-1})}^{\mathrm{\top }}={(({\mathbf{X}}^{\mathrm{\top }}\mathbf{X}{)}^{\mathrm{\top }})}^{-1}=({\mathbf{X}}^{\mathrm{\top }}\mathbf{X}{)}^{-1}$$ 입니다.

따라서, 처음 설정한대로 $\mathrm{Var}(\mathbf{y})=\mathrm{\Phi }$를 넣으면, 최종적으로 추정된 consistent model $\hat{\mathit{\beta }}$의 (co)variance matrix $\mathrm{Var}(\hat{\mathit{\beta }})$는 다음과 같습니다. $$({\mathbf{X}}^{\mathrm{\top }}\mathbf{X}{)}^{-1}{\mathbf{X}}^{\mathrm{\top }}\mathrm{\Phi }\mathbf{X}({\mathbf{X}}^{\mathrm{\top }}\mathbf{X}{)}^{-1}$$

하나 확인할 수 있는 사실은, 등분산성 가정($\mathrm{\Phi }={\sigma }^2\mathbf{I}$) 하에선 위 식이 아래처럼 단순화되어 이전 (co)variance matrix가 나왔던, 가정하의 특수한 결과라는 것입니다. $$\begin{gathered} \mathrm{Var}(\hat{\mathit{\beta }})=({\mathbf{X}}^{\mathrm{\top }}\mathbf{X}{)}^{-1}{\mathbf{X}}^{\mathrm{\top }}\mathrm{\Phi }\mathbf{X}({\mathbf{X}}^{\mathrm{\top }}\mathbf{X}{)}^{-1} \\ =({\mathbf{X}}^{\mathrm{\top }}\mathbf{X}{)}^{-1}{\mathbf{X}}^{\mathrm{\top }}{\sigma }^2\mathbf{I}\mathbf{X}({\mathbf{X}}^{\mathrm{\top }}\mathbf{X}{)}^{-1} \\ ={\sigma }^2({\mathbf{X}}^{\mathrm{\top }}\mathbf{X}{)}^{-1}{\mathbf{X}}^{\mathrm{\top }}\mathbf{X}({\mathbf{X}}^{\mathrm{\top }}\mathbf{X}{)}^{-1} \\ ={\sigma }^2({\mathbf{X}}^{\mathrm{\top }}\mathbf{X}{)}^{-1} \end{gathered}$$ 또한, 이 때의 (OLS 추정량)의 분산 ${\sigma }^2$의 식도 그냥 보여드리고 넘어갔었는데, 이는 대각 성분이 모두 같기 때문에 단순히 잔차 제곱합 $\sum e_i^2$을 자유도($N-K$)로 나눈 값($s^2$)으로 추정하여 ${\hat{\sigma }}^2=\frac{\sum e_i^2}{N-K}$ 로 한번에 표현한 식이었다는 것도 확인할 수 있습니다.

HC0’s Robustness

이분산성이 존재할 때 $\mathrm{\Phi }$는 대각행렬이지만 대각원소가 서로 다릅니다. White는 1980년에 이 heteroskedasticity에 robust하게 분산을 추정하기 위해 위 식에서 $\mathrm{\Phi }={\sigma }^2\mathbf{I}$로 term이 소거되게 하는 대신, ${\mathrm{\Phi }}_{ii}$를 잔차 제곱 $e_i^2$로 추정하는 HC0을 제안했습니다: $${\hat{\mathrm{\Phi }}}_{ii}=e_i^2\Rightarrow \hat{\mathrm{\Phi }}=\mathrm{diag}(e_i^2)$$

이를 위 (co)variance 식에 대입하면 HC0 분산 추정량이 얻어집니다: $$\mathrm{HC}0=({\mathbf{X}}^{\mathrm{\top }}\mathbf{X}{)}^{-1}{\mathbf{X}}^{\mathrm{\top }}\mathrm{diag}(e_i^2)\mathbf{X}({\mathbf{X}}^{\mathrm{\top }}\mathbf{X}{)}^{-1}$$

이는 variance 식으로부터 시작해서 얻은 식이고, 데이터로부터 얻은 matrix $\hat{\mathrm{\Phi }}=\mathrm{diag}(e_i^2)$를 사용하여 이분산성의 구체적 형태를 가정하지 않았기 때문에 안정적으로 heteroskedasticity를 고려합니다. 이처럼 앞 뒤가 같은 모양 사이에 데이터로부터 얻은 분산 추정 식이 들어간 형태가 샌드위치 같다고 해서 Sandwich Estimator라고도 부르며, 여전히 가운데 matrix 항이 diagonal 함을 통해 오차(또는 관측치)간의 상관관계는 없다는 것도 기억해야 합니다.

4.3. HC(CME) 1~3 수학적 표현

---

이 HC Standard Error가 나온 이후, 위 식으로부터 구한 분산은 소표본에서는 모델의 overfitting(과적합, 딥러닝과 머신러닝에서도 자주 대두되는 문제로 데이터가 적어 학습한 데이터에서는 오차가 지만 보지 못한 데이터에 대해서는 오차가 크게 나오는 상황, 통계 분석은 test data를 따로 두지 않기 때문에 overfitting 대처가 더욱 중요.)에 의해 residual이 작게 나오고, 이에 따라 모델의 covariance가 과소평가될 가능성이 크다는 문제가 제시되었습니다. 이에, 소표본에서도 안정적으로 추정하기 위해 White(1980)의 HC Standard Errors를 HC0로 두고, 여러 버전의 HC Standard Errors가 최근까지 다양한 학자들에 의해 연구되고 있습니다. 보통 HC0~HC5까지를 HC Standard Errors로 부르며, HC 오른쪽 숫자가 버전을 뜻하여 이 수가 커질수록 최근에 나온 것이고, 각각은 simulation 실험을 통해 소표본에서 더욱 강건함을 보이는 식으로 논문이 나옵니다. 중요한 점은, 첫 번째로 이 HC Standard Errors들은 소표본에서는 값이 다르지만 표본의 크기가 무수히 커질수록 asymptotically하게 같으며 소표본에서의 고려를 위해 값을 더 크게 만드는 추가 term을 넣어준다는 점(물론 수식적으로 필요성을 증명하지만), 둘 째로 R의 sandwich 패키지에서는 4까지 지원하고, 3을 넘는 HC 시리즈가 쓰이는 경우는 거의 볼 수 없기 때문에 3까지의 식 정도만 다루어도 충분하다는 점, 마지막으로 위 직관적인 이유로나 수식적으로나 식을 해석해 보면 거의 모든 표본 cases에 대해 버전이 클수록 분산을 더욱 크게 추정한다는 점입니다. 따라서 유의성을 보이기 위해서는 버전이 작은게 유리할 것입니다.^^ 아래에서는 HC1부터 3까지의 수식을 살펴보겠습니다. 각각의 철학에 대한 자세한 증명 과정 또한 다룰까 했지만 중요하지 않기 때문에 간단히 소개한 후 로컬에서 돌려 보실 R 예시 코드를 보여드릴 것입니다.

HC1: 소표본 편향 보정

HC1은 소표본에서 자유도 조정 인자 $\frac{n}{n-k}$를 도입하여 편향을 줄인다는 철학에서 출발하여, 잔차 $e_i^2$의 기대값이 ${\sigma }_i^2(1-h_{ii})$이므로, $E(e_i^2)\approx {\sigma }_i^2$이 되도록 HC0에 $\frac{N}{N-K}$을 나누어서(1보다 크며, n이 커짐에 따라 1로 수렴하는 term 추가) 스케일링합니다. 즉, 식은 다음과 같습니다: $$\begin{gathered} \mathrm{HC}1=\frac{N}{N-K}\cdot \mathrm{HC}0 \\ =\frac{N}{N-K}\cdot ({\mathbf{X}}^{\mathrm{\top }}\mathbf{X}{)}^{-1}{\mathbf{X}}^{\mathrm{\top }}\mathrm{diag}(e_i^2)\mathbf{X}({\mathbf{X}}^{\mathrm{\top }}\mathbf{X}{)}^{-1} \end{gathered}$$

HC2: Leverage 보정

위 HC1에서 $E(e_i^2)={\sigma }_i^2(1-h_{ii})$ 라고 이야기 하였습니다. 일단 그렇다고 하면, 기존에는 ${\sigma }_i^2$를 $e_i^2$으로 추정하였다면, 사실 $\frac{e_i^2}{1-h_{ii}}$의 기대값이 ${\sigma }_i^2$이기 때문에 HC1처럼 앞에 scalar term을 추가하는 대신, HC2는 이를 직접적으로 반영하여 $\hat{\mathrm{\Phi }}$를 $\mathrm{diag}(e_i^2)$ 대신 $\mathrm{diag}(\frac{e_i^2}{1-h_{ii}})$로 추정합니다. 이 때 $h_{ii}$은 Leverage(래버리지), $H(orh)$는 Hat Matrix라고 불리고, 식은 $\mathbf{H}=\mathbf{X}({\mathbf{X}}^{\mathrm{\top }}\mathbf{X}{)}^{-1}{\mathbf{X}}^{\mathrm{\top }}$이며, 위 $h_{ii}$은 이 matrix의 대각원소를 뜻합니다. Leverage에 대한 이야기 또한 길지만 간단하게 설명하면, $\mathbf{H}$가 Hat Matrix라고 부르는 이유는 이 $\mathbf{H}$ term에 $\mathbf{y}$를 곱하면, $$\begin{gathered} \mathbf{Hy}=\mathbf{X}({\mathbf{X}}^{\mathrm{\top }}\mathbf{X}{)}^{-1}{\mathbf{X}}^{\mathrm{\top }}\mathbf{y} \\ =\mathbf{X}\hat{\mathit{\beta }}=\hat{\mathbf{y}} \end{gathered}$$ 이 되어 $\mathbf{y}$에 hat을(추정치) 씌운다는 의미에서 비롯되었으며, 이는 기하학적으로 더욱 복잡하게 설명될 수 있지만 간단한 의미는 각 관측치가 얼마나 모델의 estimation에 영향을 주었는지를 보여주는 matrix입니다. 이 matrix의 diagonal값인 레버리지를 보며 관측치 하나하나의 영향력을 볼 수 있고, 이 값이 큰 관측치에 대해서는 아래 식에서처럼 분모를 작게 함으로써 분산을 크게 추정한다고 직관적으로 생각할 수 있습니다. $$\mathrm{HC}2=({\mathbf{X}}^{\mathrm{\top }}\mathbf{X}{)}^{-1}{\mathbf{X}}^{\mathrm{\top }}\mathrm{diag}\left(\frac{e_i^2}{1-h_{ii}}\right)\mathbf{X}({\mathbf{X}}^{\mathrm{\top }}\mathbf{X}{)}^{-1}$$

넘어가기 전에, 증명 없이 사용해왔던 수식 $E(e_i^2)={\sigma }_i^2(1-h_{ii})$ 를 간단하게만 증명해보겠습니다. 예측 오차는 실제값 $y_i$와 예측값 ${\hat{y}}_i$의 차이입니다: $$e_i=y_i-{\hat{y}}_i.$$ 잔차 벡터 $\mathbf{e}$는 다음과 같이 표현됩니다: $$\mathbf{e}=\mathbf{y}-\hat{\mathbf{y}}=(\mathbf{I}-\mathbf{H})\mathbf{y}.$$

에러 제곱의 기댓값은 잔차의 분산을 의미합니다. 잔차 벡터 $\mathbf{e}$의 공분산 행렬은 다음과 같이 계산됩니다: $$\text{Cov}(\mathbf{e})=\text{Cov}((\mathbf{I}-\mathbf{H})\mathbf{y})=(\mathbf{I}-\mathbf{H})\text{Cov}(\mathbf{y})(\mathbf{I}-\mathbf{H}{)}^T.$$ $\text{Cov}(\mathbf{y})={\sigma }_i^2\mathbf{I}$이므로 (당연히 $\sigma$는 관측치 $i$마다 다를 수 있습니다.): $$\text{Cov}(\mathbf{e})=(\mathbf{I}-\mathbf{H}){\sigma }_i^2\mathbf{I}(\mathbf{I}-\mathbf{H}{)}^T={\sigma }_i^2(\mathbf{I}-\mathbf{H}).$$ 여기서 $\mathbf{H}$는 대칭 행렬($\mathbf{H}={\mathbf{H}}^T$)이고 멱등 행렬(${\mathbf{H}}^2=\mathbf{H}$)이므로 (이 부분은 기하학적인 설명이 많이 필요해서 증명하지 않겠다만, 이들은 단지 특정 행렬들의 성질이며 이를 만족한다고 하고 넘기겠습니다.): $$\text{Cov}(\mathbf{e})={\sigma }_i^2(\mathbf{I}-\mathbf{H}).$$

$i$번째 잔차 $e_i$의 분산은 공분산 행렬의 $i$번째 대각 성분이고, $$\text{Var}(e_i)=[\text{Cov}(\mathbf{e}){]}_{ii}={\sigma }_i^2(1-h_{ii}).$$ 에러 제곱의 기댓값은 잔차의 분산과 동일하므로 (에러의 mean이 0이므로), $$E[e_i^2]=\text{Var}(e_i)={\sigma }_i^2(1-h_{ii}).$$ 가 증명됩니다.

HC3: 강화된 Leverage 보정

HC3은 Jackknife 접근법에서 영감을 받아 HC2에서 $\mathrm{\Phi }$의 분모 term을 더 강하게 $(1-h_{ii}{)}^2$로 둡니다. 직관적으로, 같은 $h_{ii}$에 대해 더 크게 분산을 추정 ($(1-h_{ii}{)}^2<1-h_{ii}$)하여, 고레버리지, 즉 outlier에 대해 더 안정적으로 소표본의 분산을 추정할 수 있습니다. $$\mathrm{HC}3=({\mathbf{X}}^{\mathrm{\top }}\mathbf{X}{)}^{-1}{\mathbf{X}}^{\mathrm{\top }}\mathrm{diag}\left(\frac{e_i^2}{(1-h_{ii}{)}^2}\right)\mathbf{X}({\mathbf{X}}^{\mathrm{\top }}\mathbf{X}{)}^{-1}$$

HC(CME)s 비교

위 내용을 간단하게 표로 정리하면 다음과 같습니다.

유형	수식	보정 요소	사용 시나리오
HC0	$\mathrm{diag}(e_i^2)$	없음	대표본
HC1	$\frac{n}{n-k}\mathrm{HC}0$	자유도	소표본 기본
HC2	$\frac{e_i^2}{1-h_{ii}}$	레버리지 1차	고레버리지 존재 시
HC3	$\frac{e_i^2}{(1-h_{ii}{)}^2}$	레버리지 2차	소표본 + 고레버리지

정리하자면, 동분산성 검정을 검정하여 이분산성이 확인되면 HC se 시리즈의 적용을 고려해야 하고, 표본 크기가 $n<50$인 경우는 HC3 또는 아래에서 다룰 bootstrap기반의 (Clustered) Bootstrap Covariance Matrix Estimation를, $n>500$인 경우는 HC0/HC1로 충분하다고 간단하게 생각해볼 수 있습니다. (물론 이는 정해진 규칙이 아니고, 시작은 항상 HC0부터 검정해보는 것이 맞습니다.)

5. (Wild) Bootstrap

마지막으로, sandwich estimator는 아니지만, 모든 모델에서 data로부터 안정적으로 분산을 추정하여 의학 연구에서 인정하는 기법 중 하나인 Bootstraping을 통한 분산 term 추정 기법을 간단하게 소개드리겠습니다. 이 method는 특정 가정으로부터 수학적인 추정식을 세워 분산을 추정하는 대신, Bootstraping 통계 기법을 사용하여 data 만으로 경험적인 분산을 추정합니다. 때문에 이번에 소개하는 (Wild) Bootstrap 은 어떤 가정이나 모델이든 사용할 수 있는, 심지어 2장에서 다룰 clustered(grouped) data 에서도 사용할 수 있는 방법입니다.(사용 대상인 모델이 정해지지 않은, 데이터 만으로 구하는 분산이기 때문에 가정이 전혀 들어가지 않아서입니다. 예를 들어 LM에서 HC se를 사용하는 경우 나머지 세 가정이 만족해야 했지만, Bootstrap은 그렇지 않습니다.) 이후에 다룰 샌드위치 추정량(Sandwich Estimator)과, 그의 OLS 버전인 HC(Heteroskedasticity-Consistent) 시리즈가 개발되었지만, 이들은 결국 대표본에서의 점근적 성질로부터 얻은 수식이고, 극단적인 쇼포본 같은 작은 클러스터 수에서는 이러한 데이터 만을 고려하여 추정하는 부트스트랩이 더 강력한 대안이 될 수도 있습니다. 이를 소개하기 전에 간단히 중요한 통계학적 resampling 기법인 잭나이프(Jackknife)와 부트스트랩(Bootstrap)을 얘기해보겠습니다. 이 둘은 모두 표본을 계속해서 얻기 힘든 상황 (현실에서는 대부분 그럴 것입니다.)에서 모수를 추정하기 위해, 이미 갖고 있는 표본의 data로부터 resampling을 하여 추정하는 통계학적 방법론입니다.

잭나이프(Jackknife)

Jackknife method는 한 번에 하나의 클러스터를 제외하고 모델을 추정 하는 것을 모든 클러스터에 대해 반복하여 분산을 계산합니다:

$$\mathrm{Var}({\hat{\beta }}_c)=\frac{n-1}{n}\sum _{c=1}^C{({\hat{\beta }}_c-\overline{\hat{\beta }})}^2,$$

여기서 ${\hat{\beta }}_c$는 c번째 클러스터를 제외한 추정치, $\overline{\hat{\beta }}$는 추정치의 평균이고, 이는 resampling 기법이지만 위에서 설명드렸던 HC3 추정량이 이 식으로부터 도출된 결과입니다. 물론, 이는 이론적인 결과이므로 각각이 서로 다른 결과를 도출할 수도 있을 것입니다. 단, resampling 기법들은 그 특성상 클러스터나 그 data point의 수가 많을 경우 계산 부하가 큽니다.

부트스트랩(Bootstrap)

Bootstrap method는 데이터를 무작위로 resampling하여 여러 버전의 데이터셋을 생성하고, 각각에서 모델을 추정한 후 추정값의 변동성을 분산으로 사용합니다. 대표적인 예시로 와일드 부트스트랩 (Wild Bootstrap)은 잔차(residual)에 랜덤 가중치를 곱해 인위적 으로 데이터의 분산을 크게 보정하고, 케이스 기반 (XY/Pairs 부트스트랩)은 기본적인 방법으로 클러스터 자체를 복원 추출하여 새 data를 생성한 뒤 구하는 과정을 반복하며, 잔차 기반 (Residual 부트스트랩)은 잔차 만을 리샘플링하고 재추정하는 방법입니다다.

Wild Bootstrap

위에서 설명드렸듯, 잔차 $e_i$에 “랜덤 가중치” $w_i$를 곱해 새로운 반응 변수를 생성하는 방법으로 다음과 같을 것입니다:

$$y^{\ast }=\hat{y}+w_i\hat{e}.$$

이때 가중치 분포 로 사용되는 대표적인 예시는 Rademacher: $w_i\in \{-1,1\}$, Mammen: $w_i=\frac{\sqrt{5}\pm 1}{2}$, Webb: 6점 분포($w_i\in \{-\sqrt{\frac{3}{2}},-\sqrt{\frac{2}{2}},-\sqrt{\frac{1}{2}},\sqrt{\frac{1}{2}},\sqrt{\frac{2}{2}},\sqrt{\frac{3}{2}}\}$)로 클러스터별 정규분포 등 다른 적용도 가능합니다. 위에서 언급하였듯, 이외에도 케이스 기반 (XY) 부트스트랩, 잔차 부트스트랩, 프랙셔널 부트스트랩(Fractional Bootstrap)등의 방법론이 있지만, data를 크게 왜곡하지 않는 Wild가 이러한 분산 추정에서 디폴트하게 고려됩니다.

6. R 예시: HC0~3 및 부트스트랩

아래 R 코드를 복사하여 로컬 환경에서 돌려보세요. 위 내용 중 “시리즈가 클 수록 강건하게 추정한다”는 점을 기억하시고 해석하면 됩니다.

## mtcars 데이터에서 mpg(연비)를 wt(차량 무게), hp(마력)으로 설명하는 회귀모형
#data(mtcars)
#model <- lm(mpg ~ wt + hp, data = mtcars)
#
## 기본 OLS 표준오차 (등분산 가정)
#summary(model)  
#
## HC 표준오차 계산
#library(sandwich)
#library(lmtest)
#
## HC 유형별 분산-공분산 행렬
#cov_hc0 <- vcovHC(model, type = "HC0")  # 기본 White 추정량
#cov_hc1 <- vcovHC(model, type = "HC1")  # 자유도 보정
#cov_hc2 <- vcovHC(model, type = "HC2")  # 잔차 스케일링
#cov_hc3 <- vcovHC(model, type = "HC3")  # 작은 표본 보정
#
## 계수 테스트 결과 비교
#coeftest(model, vcov = cov_hc0)  
#coeftest(model, vcov = cov_hc1)  
#
## 부트스트랩 분산-공분산 행렬 (기본 100회 복제)
#cov_wild <- vcovBS(fit, cluster = ~cluster_id, type = "wild", R = 1000)
#cov_xy <- vcovBS(model, R = 200, type = "xy")  # xy-쌍 부트스트랩
#
## 결과 비교
#coeftest(model, vcov = cov_wild)  
#coeftest(model, vcov = cov_xy)  
#
## 주의: 부트스트랩 결과는 실행 횟수 R, 실행 시마다 다를 수 있음

마무리하며

이번 1장에서는 Regression Analysis의 기본 개념과 (simple, multiple or General) Linear Regression의 이론적 배경, 그리고 등분산성 가정이 깨진(Heteroskedasticity) 상황에서 유용하게 쓰이는 Heteroskedasticity-Consistent Standard Errors(HC standard errors)에 대해 살펴보았습니다.

HC standard errors를 사용하면, Linear Regression 모델에서 등분산성 가정이 위배되더라도 Standard Errors(or Covariance Matrix)를 좀 더 타당하게 추정할 수 있습니다. 하지만, “오차항의 독립성, 선형성, 정규성” 등 나머지 가정이 크게 어긋난다면, 단순히 HC SE로 분산을 보정하는 것만으로는 해결되지 않습니다. HC SE가 기본 OLS 분산추정치와 크게 다르다면, “정말 모델이 맞는지”를 다시 고민하고, 필요하다면 모델을 재설정하거나 다른 방법을 모색해야 합니다.

어쨌든, Heteroskedasticity가 의심되는 상황에서 가장 먼저 고려할 만한 접근인 HC(Heteroskedasticity-Consistent) standard errors는 모델의 유의성 검정을 위해 안정적으로 분산을 추정할 때 널리 쓰이며, 의학 연구에서도 일상적으로 활용되고 있습니다. 다음 2장에서는 이런 Linear Regression을 더 확장한 Generalized Linear Model(GLM)의 개념을 본격적으로 다루고, HC standard errors의 clustered data 버전인 Cluster-robust standard error에 대해서도 다룰 예정입니다.