Ordinary Least Squares I

The Fundamentals of Linear Regression and Dealing with Qualitative Explanatory Variables

Ordinary Least Squares I

The Fundamentals of Linear Regression

여러 가지 통계기법들과 선형회귀모델에 관해 이야기하기에 앞서, 정량적 사회과학연구를 이해하기 위한 몇 가지 핵심적 개념들을 다루고, 이해해보는 시간을 가져보고자 합니다. 먼저, 대체 “함수(function)”란 무엇일까요? 함수는 도대체 왜 유용하며, 왜 사용하는 걸까요? 그리고 이 자료는 선형회귀에 관한 것인데, 그에 앞서 함수에 대한 질문을 하는 것일까요?

함수는 일군의 특성이 다른 일군의 특성과 함께 변화하는(covary) 체계적 변화를 보여주는 방식입니다.

  • 다른 말로 하자면, 하나의 함수는 어떻게 우리가 관측한 일련의 특성들이 서로 연관되어 있는지를 보여준다고 할 수 있습니다.

  • 즉, 우리가 어떤 함수를 가지고 있고, 그와 관련된 일련의 특성들을 관측했다면, 우리는 관측된 특성들이 함수적 관계 속에서 어떠한 다른 특성들과 이어지게 될지를 기대할 수 있게 됩니다.

사회과학 연구에서는 이론적 검토라는 작업이 선행되는데, 선행연구들을 바탕으로 이론적인 관계를 기대하고 경험적 데이터들이 그 기대에 따라 배열되어 나타날 것이라고 주장(가설 수립)합니다. 함수는 그 관계를 간명하게 보여주는 방식입니다. 함수는 우리가 관측한 현상들의 특성들 간 체계적 관계를 논리적으로 보여주는 지도의 역할을 하기 때문에 유용하다고 할 수 있습니다.

선형회귀분석이라는 기법에 본격적으로 들어가기에 앞서 함수의 정의를 묻는 것이 중요한 이유는 선형회귀모델이 둘, 또는 그 이상의 변수들 간의 관계를 서술하고 추론하는 데 함수를 사용하기 때문입니다. 선형회귀모델은 변수들 간의 관계를 하나의 수리적 식(equation)으로 표현하고, 그것이 바로 함수입니다.

함수적 관계를 통해 우리가 파악할 수 있는 것은

  • 첫째 우리는 원인과 결과를 상정하고 그 관계를 살펴보고자 하며,

  • 둘째 이때 원인과 결과는 서로 변화하는 관계에 있다는 것입니다.

즉, 우리는 변수와 변수라는 공변하는(covarying) 인자들 간의 관계를 보고자 하는 것입니다. 그리고 선형회귀분석은 이렇게 변화하는 관계의 모습이 선형(linear)일 것이라는 기대를 가지고 있습니다. 우리는 기대하는 선형관계를 수식으로 나타내고 그 도함수(derivative, 혹은 편미분)을 취함으로써 변수들 간의 관계를 간략하게 보여주고자 합니다.

그렇다면 도함수 또는 미분이란 무엇일까요? 우리는 \(x\)에 대한 함수 \(f(x)\)를 도함수를 취해 \(f'(x)\)의 수리적 형태로 보여줄 수 있습니다.

\[f'(x) = \lim_{h\to\infty}\frac{f(x+h)-f(x)}{h}\]

조금 더 풀어서 설명하자면, 도함수는 원래의 함수로부터 도출된 또 다른 함수로 다른 모든 변수들의 값이 일정할 때, 한 변수의 부분적 변화가 종속변수의 얼마만큼의 변화와 관련이 있는지를 보여줍니다.

  • 따라서 도함수는 변화율을 보여주는 것으로, 직관적으로는 선형관계의 기울기를 보여준다고 할 수 있습니다.

  • 도함수를 가지고 우리는 원함수(original function)가 비선형 관계일 때에도 매 순간의 \(x\)의 변화에 따른 기울기의 변화, 또는 차이를 알 수 있습니다.

매우 간략한 회귀분석의 이론

먼저, 간단한 함수적 관계를 나타내는 식이 있습니다. 이를 범죄의 모델(Model of Crime)이라고 해보겠습니다.

\[\text{범죄(Crime)} = f(\text{임금 수준(LW)})\]

위의 함수적 관계에서 각각의 항이 나타내는 것은 무엇일까요? 그리고 위의 관계는 우리에게 어떠한 것을 말해주고 있을까요? 위의 함수적 관계는 범죄(Crime)의 변화량(variation)을 설명하는 데 있어서 임금 수준(LW)이 중요하다는 것 밖에는 말해주지 않습니다.

이 모델을 가지고 우리는 그저 “주로 임금 수준이 범죄를 설명하는 데 있어서 중요합니다’’라고밖에 말할 수 없습니다. 어떻게, 얼마나 중요한지 위의 함수적 관계는 알려주지 않습니다. 다만 그 둘이 관계가 있다는 것만 시사합니다. 그렇다면 위의 함수적 관계를 조금 더 자세하게 풀어서 표현해보겠습니다.

\[\text{범죄(Crime)} = \beta_0 + \beta_1\text{임금 수준(LW)} + u\]

이 모델은 어떤 정보를 가지고 있을까요? 하나씩 뜯어보면, \(\beta_0 + u\)는 임금 수준이 0일 경우에 나타날 범죄율을 의미합니다. 그리고 임금 수준이 한 단위 증가할 때마다 범죄는 \(\beta_1\)만큼 증가하겠죠. 즉, 임금이 아예 없었을 경우의 범죄가 \(\beta_0 + u\)였다면, 임금이 한 단위 증가했을 때에는 \(\beta_0 + \beta_1\times1 + u = \beta_0 + \beta_1 + u\)가 되는 것입니다. 조금 더 자세히 뜯어보면,

  • 범죄(Crime): 종속변수(dependent variable)

  • 임금(LW): 설명변수(explanatory variable)1

  • \(\beta_0\): 임금이 0일 때의 범죄 수준

  • \(\beta_1\): 임금이 한 단위 변화할 때 나타나는 범죄의 변화

  • \(u\): 오차항(error term)

그렇다면 이 선형모델을 가지고 우리는 무엇을 할 수 있을까요? 바로 모델을 통해 우리가 관심을 가지고 있는 변수의 한 단위 변화와 연관되는 종속변수의 변화를 포착하는 것입니다. 이를 한계효과(marginal effect)라고 하는데요, 단순선형회귀모델, 혹은 다중선형회귀모델에서는 각 변수들이 +로 결합되어 있습니다. 논리적으로 이 +는 서로가 독립적이라는 OR의 의미를 지니고 있습니다. 따라서 이와 같이 +로 변수들이 연결된 모델을 가산모형(additive models)이라고 하고, 이때의 한계효과는 각 변수들의 \(\beta_i\)라고 보시면 됩니다.

주어진 모델에서 한계효과를 수리적으로 나타내면 다음과 같습니다. \[\text{Marginal effects:} \frac{\partial\text{(범죄)}}{\partial\text{(임금 수준)}} = \beta_1\] 그리고 이때 우리는 한계효과를 통해 특정한 임금 수준에서의 범죄 수준을 예측할 수 있게 됩니다. \[\text{Predictions:} \widehat{\text{범죄}} = \beta_0 + \beta_1\widehat{\text{임금 수준}}\]

자, 그럼 차근차근 하나씩 구체적인 의미들을 살펴보겠습니다.

\(\beta_1\)에 대하여

다른 모든 조건들이 일정할 때(Ceteris paribus)라는 말이 있습니다. 통계방법에서는 다른 변수들의 효과를 실험실에서처럼 완벽하게 통제할 수 없기 때문에, 다른 변수들의 값을 동일하게 고정시켜둔 상태에서 우리가 관심을 가진 변수의 변화가 종속변수의 변화와 어떠한 관계에 놓여있는지를 살펴보게 됩니다. 따라서 \(\beta_1\)은 다른 모든 변수들이 동일한 값으로 고정되어 있을 때, 임금 수준의 값만 변화할 때 나타나는 범죄의 변화를 의미합니다.

그렇다면 \(u\)는 대체 무엇일까요? 이 모델에서 임금 수준과 \(u\)는 모두 확률변수입니다. 그리고 임금 수준이 변화하더라도 \(u\)는 변화하지 않습니다. \(u\)는 임금 수준으로 설명되지 않는 범죄의 나머지 부분으로, 우리가 관측하지 못하는(unobserved) 요소들입니다.

이론적으로는 KKV를 다시 한 번 불러와 그 개념대로 설명하면 이해가 쉽습니다. 우리는 종속변수, 범죄에 대해 체계적으로 연계되어 있을 것이라고 생각하고 체계적 요인(systematic factor)으로 임금 수준이라는 변수를 모델에 투입하였습니다. 그렇다면 우리의 이론적 기대가 다른 체계적 요인의 가능성을 제시하지 않는다면, 이론적으로는 범죄의 체계적 변화는 임금 수준으로 다 설명이 되고 체계적으로 설명되지 않는 비체계적 요인들만이 남게 될 것입니다. 이 비체계적 요인은 어떠한 경향성을 보이지 않는, 무작위(randamized)할 것이며 우리는 이를 확률변수(random variable)이라고 합니다. 마찬가지로 이 오차항—\(u\)는 설명변수와 +로 연결되어 있기 때문에 모델 안에서는 적어도 그들은 서로 관계가 없습니다(없어야만 합니다). 이를 수리적 공식으로 표현해보면 다음과도 같습니다. \[E(u|x) = E(u)\]

주어진 \(x\)라는 조건에서 \(u\)의 기대값은 그냥 \(u\)의 기대값과도 같다, 즉 \(x\)는 \(u\)와 어떤 관계도 맺고 있지 않다는 것입니다. Wooldridge (2016: 22)는 “\(\cdots\)우리는 관측할 수 없는 \(u\)가 설명변수 \(x\)와 관계되지 않도록 제한하는 가정을 만들었을 때만, 데이터의 무작위 표본(random sample)로부터 신뢰할 수 있는(reliable) 추정치, \(\beta_0, \beta_1\)을 확보할 수 있다.”고 언급하고 있습니다. 또한, “그러한 제한이 없이 우리는 다른 조건들이 일정할 때의 효과 \(\beta_1\)를 추정할 수 없다”고 서술하고 있습니다.

여기서 하나 더 중요한 것이 바로 “무작위 표본”이라는 것입니다. 무작위화(randomization)은 우리가 설명변수들과 오차항 간의 독립성 가정을 수립할 수 있도록 해주는 핵심적인 개념입니다. 만약 어떤 표본이 무작위로 뽑힌 것이 아니라 어떠한 기준에 의해서 뽑혔다면, 우리는 그 표본으로부터 얻은 결과가 실로 모집단에서의 관계를 보여주는지, 혹은 표본 선정에 관여한 기준으로 인해 나타난 경향성인지를 구분할 수 없을 것입니다.

\(x\)와 \(u\)에 대한 독립성 가정에 더해 Wooldridge (2016: 22)와 같이, 여기서도 \(u\)와 \(x\)가 어떻게 관계를 가지는지에 대한 가정을 별도로 수립하기 전까지는, 항상 \(u\)에 대하여 “모집단에서의 \(u\)의 평균값은 0이다”라는 가정을 수립하도록 합니다. 이 가정은 \(E(u) = 0\)이라고 표현할 수 있습니다.

\(E(u|x)\)의 의미

\(E(u|x)\)에 다시 한 번 살펴보겠습니다. \(E(u|x)\)는 \(x\)라는 설명변수가 주어진 조건 하에서의 \(u\), 오차항의 기대값을 보여주는 수리적 표현입니다. \(u, x\)가 확률변수이기 때문에, 우리는 \(x\)의 어떤 값이 주어졌다고 할 때의 \(u\)의 조건분포를 정의할 수 있습니다 (Wooldridge 2016: 22). 다른 표현으로는 \(x\)가 존재할 때의 \(u\)의 조건확률이라고 합니다. 단순선형회귀모델에서 우리는 이 조건확률을 0, \(E(u|x)=0\)이라고 가정합니다. 이론적으로 다시 한 번 살펴보았다면, 이번에는 좀 더 실질적인 맥락에서 이 가정을 살펴보겠습니다. 설명변수에 대한 오차항의 조건확률이 0이라는 것은 현실 속에서는 성립(혹은 충족)되기 어렵습니다. 왜냐하면 현실 속에서는 우리가 관측가능한 종속변수에 영향을 미칠 ‘수도’ 있는 여러 ‘관측 불가능한’ 혹은 ‘우리가 생각하지 못한’ 잠재적 요인들이 존재할 수 있기 때문입니다.

또한, 우리가 모집단에 대한 회귀 함수(PRF)를 수립한다고 하더라도 우리가 확보 가능한 자료는 어디까지나 관측가능한 표본에 지나지 않습니다. 표본은 우리가 모집단에서 추출할 때마다 달라질 수 있고, 그것은 표본이 가지는 본연의 한계입니다. 따라서 어디까지나 표본을 사용한다는 점에서 우리는 본연적 불확실성에서 자유로울 수 없고, 이 오차항에 대한 가정은 충족되기 힘듭니다. 만약 이 오차항에 대한 가정이 심각하게 위배된다면, 우리는 그러한 표본으로부터 얻은 추정량들에 대해 “편향되었다”(biased)라고 할 수 있을 것입니다.

논리적으로 \(\times\)는 AND로 두 가지 사건이 동시에 발생하였다는 것을 의미합니다. 그러나 이미 우리는 앞서 \(E(u)\)라고 가정한 바 있습니다. 따라서 어떤 \(x\)의 값이 주어지던 간에 상관없이, \(u\)의 기대값은 0이 됩니다. 이는 \(x\)와 \(u\)가 서로 독립적인 사건으로 \(u\)가 \(y\)와 \(x\)와 공변하지 않는다는 것을 의미합니다.

\[\begin{equation*} \begin{aligned} Cov(x, u)& = E((x-\mu)(u-\nu)) = E(xu - x\nu - u\mu + \mu\nu)\\ & = E(x\cdot u) - E(x\nu) - E(u\mu) + E(\mu\nu)\\ & = E(x\cdot u) - \nu\mu - \mu\nu + \mu\nu\\ & = E(x\cdot u) - \nu\mu\\ \end{aligned} \end{equation*}\]

위의 수리적 유도에 따라 논리적으로 \(E(u|x)=Cov(x,u)=E(x\cdot u)\neq 0\)가 된다면, 설명변수는 종속변수에 영향을 미치는 오차항과 연관성이 있게 됩니다. 이 경우, 우리는 순수한 \(x\)와 \(y\)의 관계와 \(y\)에 대한 \(u\)의 관계를 구별해내기 어렵습니다. 이것이 \(x\)와 \(u\)가 서로 상관성이 있다고 할 때, 그 추정 결과가 편향된다고 말하는 근거입니다. 이렇게 여러 번 반복하는 이유는 그만큼 중요해서겠죠?

모집단 회귀 함수(Population Regression Function; PRF)

이전 자료에서 추론(inference) 파트에서 살펴보았듯, 우리가 관심있는 것은 표본에서의 변수들의 관계가 아닙니다. 우리는 어디까지나 그 표본이 모집단에 대해 대표성 있는 관측가능한 사례라고 간주하고, 그 표본에서의 관계가 관측불가능한 모집단에서도 유의미하게 나타날 것이라고 기대하는 것입니다. 따라서 어떠한 모델을 구성할 때, 우리는 두 가지 수준에서 모델을 생각해볼 수 있습니다. 바로 모집단 수준에서의 모델과 표본 수준에서의 모델입니다. PRF는 모집단 수준에서의 회귀모델을 의미합니다. \[\begin{equation*} \begin{aligned} E(\text{범죄}|x, u)&= \beta_0 + \beta_1\text{임금 수준} + E(u|x)\\ &=\beta_0 + \beta_1\text{임금 수준}\\ &=E(\text{범죄}|x) \end{aligned} \end{equation*}\] 문제는 우리는 관측할 수 없지만 분명 존재할 \(\beta_0\)와 \(\beta_1\)에 대해 알고자 한다는 것입니다. 따라서 우리는 가지고 있는 데이터로부터 모집단의 \(\beta_0\)와 \(\beta_1\)을 추론할 필요가 있습니다. 그리고 표본에서 이 모집단의 \(\beta_0\)와 \(\beta_1\)을 보여주는 통계치(statistics)이자 예측치(estimates)는 \(\hat{\beta_0}\)와 \(\hat{\beta_1}\)라고 할 수 있습니다. 하지만 데이터만 가지고 잇을 때에는 \(\text{범죄}_i\)와 \(\text{임금 수준}_i\)라는 각 변수들의 개별 관측치들은 가지고 있는 것이지만 아직 그 두 변수의 관계를 보여주는 \(\beta_0\)와 \(\beta_1\)의 표본통계치, \(\hat{\beta_0}\)와 \(\hat{\beta_1}\)는 알지 못합니다.

추정-조건들들로부터 도출

여기에서는 두 가지 가정들을 살펴봅니다. 이 \(\beta_0\)과 \(\beta_1\)을 얻기 위해서는 이 두 가정들이 충족되어야 합니다. 하지만 관측된 자료들을 이용하는 연구에서는 사실 이 가정들이 모두 완벽하게 충족되기란 어렵습니다. 따라서 우리는 확보할 수 있는 데이터셋에 관해 매우 조심스러운 접근을 취할 필요가 있습니다. 과연 내가 가진 데이터는 이 두 가정을 충족시키는지, 그리고 충족시키지 않는다면 어떻게, 얼마만큼 위배하고 있는지를 면밀히 살펴보아야 합니다.

가정1.과 제약1.

제일 먼저 수립하는 가정과 그에 따른 제약은 바로 오차항의 기대값이 0으로 수렴한다는 것입니다. \(E(u) = 0\)으로 표현할 수 있겠죠? 그리고 오차항이 종속변수를 설명하기 위한 설명변수(\(x\))와 그 효과(\(\beta_1\)), 그리고 설명변수가 0일 경우의 값(\(\beta_0\))이라는 것을 감안할 때, 위의 가정은 다음과 같이 표현할 수 있습니다. \[\begin{equation*} \begin{aligned} E(u)&= 0\\ &= E(y-\beta_0-\beta_1x) \end{aligned} \end{equation*}\] 모집단 수준에서의 오차항에 대한 가정이 다음과 같고, 이것이 참이라고 한다면, 개별 관측치로 나타나는 표본에서의 오차항에 대한 대응—잔차(residuals)에도 그와 같은 가정을 적용할 수 있을 것입니다. \[\begin{equation*} \begin{aligned} E(u_i)&= 0\\ &= E(y_i-\hat{\beta_0}-\hat{\beta_1}x_i) &= \frac{1}{N}\sum^{N}_{i=1}(y_i-\hat{\beta_0}-\hat{\beta_1}x_i) = 0 \end{aligned} \end{equation*}\] 모집단에서 우리는 기대값(expected values)이라는 표현을 사용할 수 있지만, 표본 수준에서는 관측가능한 값, 평균(mean)으로 이를 보여줄 수 있습니다.

가정2.와 제약2.

두 번째 가정이란 설명변수와 오차항 간의 공변(covariation)이 존재하지 않는다는 것입니다. 즉, \(Cov(x, u)=0\)이라는 것이죠. 이것이 의미하는 게 무엇일까요?

\[\begin{equation*} \begin{aligned} E(u|x)&= Cov(x, u) = E(x\cdot u)=0\\ &=E(x(y=\beta_0-\beta_1x)) \end{aligned} \end{equation*}\]

모집단 수준에서 주어진 설명변수가 있을 때의 오차항의 기대값은 곧 둘 간의 상관성을 묻는 것과 같습니다. 그리고 논리적 기호로써 이는 \(\times\)이며, AND고 둘이 서로 동시에 작동한다는 것을 의미합니다. 그럼 표본 수준에서 살펴볼까요?

\[\begin{equation*} \begin{aligned} E(u|x)&= 0\\ &= E(x_i(y_i=\hat{\beta_0}-\hat{\beta_1}x_i))\\ &= \frac{1}{N}\sum^{N}_{i=1}(x_i(y_i=\hat{\beta_0}-\hat{\beta_1}x_i)) = 0 \end{aligned} \end{equation*}\]

이제 위의 두 가정들로부터 선형회귀모델의 \(\beta_0\)과 \(\beta_1\)을 추정하기 위한 \(\hat{\beta_0}\)과 \(\hat{\beta_1}\)을 얻는 과정을 유도해보도록 하겠습니다.

제1단계

가정 1(\(\frac{1}{N}\sum^{N}_{i=1}(y_i-\hat{\beta_0}-\hat{\beta_1}x_i)=0\))은 \(\bar{y}=\hat{\beta_0}+\hat{\beta_1}\bar{x}\)를 유도한다. 이때, \(\bar{x}\)와 \(\bar{y}\)는 각각 \(x\)와 \(y\)의 평균을 의미합니다.2 즉, 선형회귀모델은 “종속변수와 설명변수의 평균을 지나는” 선을 그리게 됩니다. 위의 식에 따라 우리는 \(\hat{\beta_0} = \bar{y}-\hat{\beta_1}\bar{x}\)라는 결과를 얻을 수 있습니다.

제2단계

가정 2는 표본에서의 설명변수와 오차항 간 공변이 없을, \(\sum^{N}_{i=1}(x_i(y_i=\hat{\beta_0}-\hat{\beta_1}x_i)) = 0\)를 의미합니다. 그리고 이 경우, 식을 풀어서 \(\sum^{N}_{i=1}x_i(y_i-\bar{y})=\hat{\beta_1}\sum^{N}_{i=1}x_i(x_i-\bar{x})\)로 표현할 수 있습니다. 여기서부터는 수리적으로 계속 유도해내는 거라서 굳이 하나하나 다 이해하실 필요가 있다고는 하지 않겠습니다만, 다시 쓰면\(\cdots\) \(\sum^{N}_{i=1}x_i(y_i-\bar{y})=\sum^{N}_{i=1}(x_i-\bar{x})(y_i-\bar{y})\)로 보여줄 수 있고, \(\sum^{N}_{i=1}x_i(_i-\bar{x})=\sum^{N}_{i=1}(x_i-\bar{x})^2\)라고 할 수 있습니다. 결론적으로 이 일련의 복잡한 식들을 통해서 우리는 앞의 \(\hat{\beta_0}, \hat{\beta_1}\)에 대입, 이 두 통계치를 \(x\)와 \(y\)로 보여줄 수 있습니다. \[\begin{equation*} \begin{aligned} &\hat{\beta_1}=\frac{\sum^{N}_{i=1}(x_i-\bar{x})(y_i-\bar{y})}{\sum^{N}_{i=1}(x_i-\bar{x})^2}\\ &\hat{\beta_0} = \bar{y}-\hat{\beta_1}\bar{x} \end{aligned} \end{equation*}\]

언제 \(\beta_1\)의 값이 커질까?

\(\beta_1\)은 간단히 얘기하면 기울기, \(x\)와 \(y\) 간의 관계 양상을 보여주는 추정치입니다. 그런데 이 \(\beta_1\)을 추론하기 위하여 우리가 구하는 표본의 통계치, \(\hat{\beta_1}\)은 자료의 척도(scale)와 밀접하게 관련되어 있습니다. Wooldridge (2016: 26)에서 \(\hat{\beta_1}\)를 구하기 위한 공식을 \(\hat{\beta_1} = \hat{\rho}_xy\cdot \frac{\hat{\sigma}_y}{\hat{\sigma}_x}\)이라고 보여주고 있습니다. 이것이 의미하는 것은 무엇일까요? 만약 \(y\)의 척도, 단위가 \(x\)의 척도, 단위에 비해 훨씬 크다면, \(\hat{\beta_1}\)는 상대적으로 매우 커진다는 것을 의미합니다. 분수니까요. \(x\)의 분산이 \(y\)의 분산보다 훨씬 작을 경우에도 마찬가지로 상대적으로 큰 \(\hat{\beta_1}\) 값을 가질 수 있다는 것을 보여줍니다.

추정-잔차의 최소화

위에는 \(\beta_0, \beta_1\)을 추정하기 위한 \(\hat{\beta_0}, \hat{\beta_1}\)을 구하기 위해 필요한 두 가지 가정을 살펴보았습니다. 이번에는 또 다른 관점에서 \(\hat{\beta_0}, \hat{\beta_1}\)을 살펴보겠습니다.

  • 예측값(fitted values): 주어진 설명변수 \(x\)를 우리가 구한 \(\hat{\beta_0}, \hat{\beta_1}\)에 대입하여 종속변수를 예측한 값입니다. \(\hat{y_i} = \hat{\beta_0} + \hat{\beta_1}x_i + \hat{u}_i\)라고 할 수 있습니다.

  • 잔차(Residuals, \(\hat{u}_i\))는 모집단의 오차항 \(u\)와 대응되는 표본통계치입니다.

  • 잔차/오차항이란 우리가 예측하지 못한, 혹은 관측하지 못한 것—불확실성입니다. 불확실성을 가능한 한 줄이는 것은 어찌보면 사회과학의 본령일지도 모릅니다.

    • 잔차의 평균이 0이라고 했기 때문에 단순하게 이 잔차들을 더한 값을 최소화하는 것은 의미가 없습니다. 비체계적 요인들이 서로 상쇄(offsets)된다고 할 때, 잔차를 최소화하는 방법은 크게 두 가지가 있을 것입니다. 절대값을 구해서 그 모든 절대값들의 합이 최소가 되게 하는 것, 그리고 제곱의 합이 최소가 되게 하는 것입니다. 이 둘의 차이는 개념적으로는 없습니다. 다만 후자—제곱합이 수리적으로 더 계산이 용이하기 때문에 이 방법을 택하는 것일 뿐입니다.

  • 잔차의 제곱합을 최소화하는 것은 앞서 가정들에서 \(\hat{\beta_0}, \hat{\beta_1}\)을 도출해낸 것과 같은 결과로 이어집니다.

표본 회귀 함수 (Sample regression function, SRF)

표본 회귀 함수는 모집단의 모집단 회귀 함수와 동일한 형태를 가집니다(\(\hat{y}=\hat{\beta_0}+\hat{\beta_1}x\)). 함수적 관계라는 측면에서는 \(\hat{y} = f(x, \hat{\beta_1}, \hat{\beta_0})\)라고도 표현할 수 있습니다. 우리는 이 표본 회귀 함수를 가지고 도함수(편미분) 등을 통해 우리가 원하는 결과들을 확인할 수 있습니다.

한 번 앞의 범죄-임금 수준 함수를 가지고 살펴보겠습니다. \(y\)가 범죄, \(x\)가 임금 수준이라고 하고 \(\hat{\beta_0}=11, \hat{\beta_1}=-0.5\)라고 해봅시다.

  • 우리는 \(x\)가 20일때와 2일 때의 예측값의 차이도 살펴볼 수 있고; \(f(x=20, \cdot) - f(x = 2, \cdot)\)

  • 도함수를 통해서 한계효과를 계산할 수도 있으며; \(\frac{\partial f(x, \cdot)}{\partial x}\)

  • \(x\)가 0일 때의 \(y\)의 값을 구할 수도 있고; \(f(x = 0, \cdot)\)

  • 변화율을 살펴볼 수도 있습니다; \(\frac{f(x=20, \cdot) - f(x = 2, \cdot)}{f(x = 2, \cdot)}\times100\)

모집단 회귀 함수와 표본 회귀 함수의 이해

이쯤에서 한 번 내용들을 정리하고 넘어가겠습니다. \(\beta_0, \beta_1, \hat{\beta_0}, \hat{\beta_1}, x, x_i, u, \hat{u_i}\), 그리고 \(\bar{x}\)와 같이 여러 가지 수리적 표현들이 등장했는데요, 이들 각각은 대체 무슨 의미를 지니고 있는 걸까요? 그리고 이 각각의 값들은 확정적(fixed)인걸까요 아니면 확률적(random)일까요? 이를 설명하기 위해서는 모집단과 표본의 관계에 대한 명확한 관점이 필요합니다.

우리는 모집단과 표본을 구별해야만 합니다. 그리고 모집단 회귀 모델과 표본 회귀 모델 역시 구분해야 합니다. 비록 우리가 모집단에 관심을 가지고 있다고 하더라도 현실 속에서 우리가 가질 수 있는 것은 오직 표본입니다. 따라서 우리는 관측불가능한 모집단의 특성들을 관측가능한 표본들을 통해서 추론하고자 합니다. 이는 이론적으로 모집단에 대해 우리가 관심을 가지는 값들—모수(parameters)는 이미 주어진 것들로 확정정인 값들이라는 것을 의미하며, 표본의 통계치들은 모집단에서 추출하는 과정 속에서 불확실성을 내재하게 되므로 표본에 따라 다르게 나타날 수 있습니다. 즉, 확률적입니다.

  • PRF는 모집단의 수준에서 어떠한 함수적 관계를 가지고 하나의 설명변수가 종속변수와 관계되어 있다고 가정할 때, 두 가지 모수를 구할 수 있다는 것을 보여줍니다. 바로 \(\beta_0\)와 \(\beta_1\)입니다.

    • \(\beta_0\)는 절편에 대한 모수값입니다. 우리가 상정한 설명변수의 값이 0일 때의 종속변수의 값을 보여주는 것이죠. 모집단의 속성을 보여주는 이 \(\beta_0\)은 모집단 수준에서는 확정적인 값입니다.

    • \(\beta_1\)은 기울기에 대한 모수값입니다. \(\beta_1\)은 다른 모든 조건들이 일정할 때, \(x\)와 \(y\) 간의 관계를 보여줍니다. \(\beta_0\)와 마찬가지로 확정적인 값입니다.

    • \(u\)은 관측되지 않았지만 종속변수에 영향을 미칠 것으로 기대되는 요인들입니다. 따라서 \(u\)는 \(x\)에 의해 설명되지 않는 \(y\)의 일부라고 표현할 수 있습니다. 하지만 이 \(u\)는 비체계적이고 확률적이기 때문에 서로 상쇄되어 0으로 수렴할 것이라고 가정할 수 있습니다. 모집단의 수준에서는 확률변수이지만 동시에 이 오차항도 고정된 값입니다. 0으로 수렴하는 고정된 값이지요.

    • \(x\)는 종속변수를 설명할 것으로 기대되는 변수로 마찬가지로 모집단 수준에서는 확정적인 값을 가집니다. 종속변수 \(y\)도 마찬가지겠죠?

  • SRF는 표본 회귀 함수로 PRF와 같은 형태를 취하지만 다만 표본 수준의 논의라는 점이 다를 뿐입니다.

    • \(\hat{\beta_0}, \hat{\beta_1}\)은 모수값, \(\beta_0, \beta_1\)을 추정하기 위해 구한 표본에서의 추정치입니다. \(\hat{\beta_0}\)은 절편에 대한 통계치이며 \(\hat{\beta_1}\)은 기울기에 대한 표본의 통계치입니다. 따라서 이 값들은 표본을 어떻게 뽑느냐에 따라서 매번 달라질 수 있으므로 “확률적”인 값들입니다.

    • \(x_i, \hat{u_i}\)은 각각 설명변수의 개별 관측값들과 주어진 설명변수 값에 대한 예측된 종속변수의 값과 실제 관측된 종속변수 값의 차이를 의미합니다. \(u_i\)는 \(y_i - \hat{y}\)라고 할 수 있습니다. 두 값 모두 표본 수준의 논의이므로 확률적인 값들입니다.

    • \(\bar{x}\)는 표본에서 관측된 설명변수의 평균으로 표본에 따라 달라집니다. \(\bar{y}\)도 마찬가지로 표본의 종속변수 평균이며 표본에 따라 달라지므로, 이 두 값 모두 확률적인 값들입니다.

일종의 “총합들” (sums)

우리는 다양한 종류의 제곱합들을 통해 회귀모델의 특성을 보여줄 수 있습니다.

  • 총 제곱합(Total sums of squares): SST = \(\sum^{N}_{i=1}(y_i-\bar{y})^2\)

  • 설명된 변동량의 제곱합(Explained sums of squares): SSE = \(\sum^{N}_{i=1}(\hat{y_i}-\bar{y})^2\)

  • 잔차의 제곱합(Residual sums of squares): SSR = \(\sum^{N}_{i=1}(\hat{y_i}-y_i)^2\)

  • 모델이 얼마나 잘 들어맞는지를 보여주는 측정지표인 결정계수, \(R^2 \equiv \frac{SSE}{SST}\)

이같은 변동량에 대한 논의는 \(y\), \(x\), 그리고 \(u\)에 대한 벤다이어그램을 통해 도식적으로도 보여줄 수 있습니다. 물론 변수가 많아질수록 이렇게 도식화하는 것은 점점 어려워지기는 합니다. 그러나 여기서는 단순선형회귀모델을 가정하고 간단한 그림을 그려보겠습니다.

R로 계산해보기

패키지 불러오기

library(ezpickr)
library(ggplot2)
library(futurevisions)
library(tidyverse)
theme_set(theme_bw())

예제 1. Quality of Governance Basic Dataset

앞에서의 이론적 논의들을 한 번 R을 통해 경험적 데이터를 가지고 계산해보도록 하겠습니다. 구체적으로는 주어진 데이터로 \(\hat{\beta_0}, \hat{\beta_1}\), 그리고 \(R^2\)를 구해볼 겁니다.

예제로 사용할 데이터는 2016년도 국가 수준의 집합자료들입니다. 사실 교차사례-시계열 자료로 구축할 수도 있지만 그렇게 되면 너무 복잡하니까 2016년도를 따로 떼어내서 교차사례 자료로 만들었습니다. 이 자료에서 국가의 노령화 정도와 경제 간의 관계에 관심을 가지고 있다고 하겠습니다.

분석단위는 따라서 국가가 되겠고, 노령화 수준이 높은 국가일수록 생산성이 낮아져서 경제 수준도 낮을 것이라는 이론적 기대를 가지고 모델을 만들어보겠습니다. 물론, 원래 모델을 구성하려면 더 세심한 고려가 필요하겠습니다만, 여기서는 R 예제를 보이는 것이 목표이니만큼 단순화된 모델 구축으로 넘어가겠습니다. 여하튼, 이 예제에서 상정된 이론적 모델은 다음과 같습니다: \(\text{경제수준} = \beta_0 + \beta_1\text{(노령화 수준)} + u\).

분석을 위해 사용한 자료는 The Quality of Governance Basic Dataset 입니다. 주요 변수로는 2016년도의 2010년 미국 달러 고정으로 계산된 1인당 GDP(wdi_gdpcapcon2010)와 노동가능인구 대비 노인 인구로 계산된 노령화 지수(wdi_agedr)입니다.

  • 먼저, 1인당 GDP는 는 국내총생산을 그 해 중간의 인구 수로 나눈 값이기 때문에 단순 GDP 보다는 국가의 경제 규모를 더 잘 보여주는 지표라고 판단하여 사용하였습니다. 인구 100만명인데 100만 달러 버는 나라와 1000만명인데 300만 달러를 버는 나라를 비교해보았을 때, 후자가 더 부유하다고 단순하게 말하기는 어려울 겁니다.

  • 그 다음으로 생산성에 있어서 노령 인구가 젊은 인구에 비해서 더 낮을 것이라고 볼 수 있기 때문에 노령인구 수준이 더 높을수록 경제 수준이 더 낮을 것이라고 기대합니다. 따라서 노령인구 비율은 노동 가능인구에 대해 측정됩니다. 노동가능인구는 약 15세에서 64세 사이의 인구이며, 노령인구는 64세 초과의 인구입니다.

QOG <- ezpickr::pick(file = "http://www.qogdata.pol.gu.se/data/qog_bas_ts_jan22.dta")
QOG_summary <- 
  QOG |> select(ccodecow, cname, year, 
                 wdi_agedr, wdi_gdpcapcon2010) |>
  dplyr::filter(year==2016) |> drop_na()

y <- QOG_summary$wdi_gdpcapcon2010
x <- QOG_summary$wdi_agedr
model <- lm(y ~ x, data=QOG_summary)
b1.hat <- (sd(y) / sd(x)) * cor(y, x)
b0.hat <- mean(y) - b1.hat * mean(x)
y.hat <- b0.hat + b1.hat * x # Predicted y given x
SSE <- sum((y.hat - mean(y))^2)
SSR <- sum((y - y.hat)^2)
R2 <- SSE / (SSE + SSR) 
round(b1.hat, 3)
[1] -429.507
round(b0.hat, 3)
[1] 38953.91
round(R2, 3)
[1] 0.155

각각의 결과들을 한 번 살펴보고 해석해볼까요? 편의상 소수점 없이 정수로만 나타내겠습니다.

  • 먼저 \(\hat{\beta_1}\)은 약 -418.258로 노령화 비율이 한 단위 증가할수록 \(\hat{\beta_1}\)만큼 1인당 GDP가 감소할 것이라는 것을 보여줍니다.

  • \(\hat{\beta_0}\)는 38207.65로 \(x\), 노령화 비율이 0일 때의 1인당 GDP의 규모를 보여줍니다.

  • 그렇다면 모델의 설명력은? \(R^2\)로 볼 수 있습니다. 약 0.187이네요. 비율이니까 다시 말하면, 노령화 비율로는 1인당 GDP의 변화를 약 14.8% 설명할 수 있다, 이렇게 이해하실 수 있을 것 같습니다. 한편, 이것은 단순하게 1인당 GDP의 평균으로 1인당 GDP의 변화를 설명하는 것보다 더 나은 설명력을 모델이 제공한다는 것을 의미합니다.

또, 척도에 따라 이 \(\beta_1\)의 값이 민감하게 변화한다는 것을 언급했었는데요. 이번에는 주어진 SRF, \(\hat{y} = \hat{\beta_0} + \hat{\beta_1}x + \epsilon\)에서 설명변수 \(x\)에 로그값을 취한 다른 측정척도의 \(z\)를 모델에 투입해보겠습니다.

  • 앞서 1인당 GDP와 노령화 비율로 구성한 SRF에 따르면, \(x\) 한 단위의 변화는 \(\hat{\beta_1}\)만큼의 종속변수 변화와 관계가 있습니다. 그리고 종속변수의 예측값은 \(x=1\)일 때, \(\hat{y} = \hat{\beta_0 + \hat{\beta_1}}\)로 나타낼 수 있습니다. 마찬가지로 \(x=2\)일 때는 \(\hat{y} = \hat{\beta_0} + 2\hat{\beta_1}\)로 보여줄 수 있고, \(x=1\)일 때와 \(x=2\)일 때의 차이가 \(\hat{\beta_1}\)입니다.

  • 그렇다면 만약 \(z\)를 투입하게 되면, \(z\) 한 단위의 변화가 종속변수의 변화와는 어떤 관계를 맺게 될까요? 도식적으로 보면 다음과 같을 겁니다. \(y = \hat{\beta_0} + \hat{\beta_1}\text{log}x\).

    • 이 경우에는, \(y + \partial y = \hat{\beta_0} + \hat{\beta_1}\text{log}(x + \partial x)\)라 할 수 있고,

    • 좀 더 자세하게는 아래와 같이 표현할 수 있습니다.

\[\begin{equation*} \begin{aligned} \partial y&= \hat{\beta_1}(\text{log}(x + \partial x) - \text{log}x)\\ &= \hat{\beta_1}\text{log}(1 + \frac{\partial x}{x}) \end{aligned} \end{equation*}\]

만약 \(\frac{\partial x}{x}\) is 작다면, t\(\text{log}(1 + \frac{\partial x}{x}) \approx \frac{\partial x}{x}\)일 거라고 생각해볼 수 있습니다.

여기서 우리는 \(x\)의 1퍼센트 변화가 \(y\)의 \(\frac{\hat{\beta_1}}{100}\)만큼의 변화와 관계가 있을 것이라고 이해할 수 있습니다. 즉, 변수의 측정척도는 연구에 있어서 중요합니다. 실질적으로 단순히 변화량을 보여주는 것만이 아니라 때로는 변화율을 보여주는 것이 실질적으로 더 많은 정보를 제공할 수 있습니다.

\(x\)에 로그값을 취하면 우리는 \(x\)의 한 단위 변화가 아니라 1% 변화가 \(y\)에 미치는 효과를 볼 수 있게 됩니다: \(\frac{\hat{\beta_1}}{100}\). 예를 들어, \(x\)가 연간 소득이라고 할 때, 사실 우리는 연간 소득의 1달러 변화가 궁금하다기 보다는 연간 소득의 비율 변화가 \(y\)와 어떠한 관계가 있는지 궁금할 것입니다. 따라서 연구목적에 따라 우리는 측정척도들을 변화시킬 수 있으며, 주의해야할 것은 그 경우 변화된 척도에 맞춰 해석해줘야 한다는 것입니다.

예제 2. ggplot2의 diamonds dataset

먼저 이 예제에서 사용할 종속변수는 바로 다이아몬드 데이터의 가격(price) 변수입니다. 그리고 이 예제에서는 다이아몬드의 무게(carot)와 가격 간의 관계를 살펴볼 것입니다. 이번에는 통계치보다는 시각화된 플롯으로 살펴보겠습니다. x축에는 무게, y축에는 가격을 설정하였습니다.

df <- diamonds 
df |>
  ggplot(aes(x = carat, y = price)) +
  geom_point(alpha = .05, color = futurevisions("mars")[3]) +
  geom_smooth(method = 'lm', se = F, color = futurevisions("mars")[1]) +
  coord_cartesian(ylim = c(-1000, 22000)) +
  labs(
    x = 'Carat of Diamond', y = 'Price of Diamond',
    title = "Effect of Carat on Diamond Price"
  ) + 
  scale_y_continuous(labels = scales::dollar_format())

결과를 보면 선형이라기 보다는 지수함수와도 같은 형태로 우측으로 가파르게 상승하는 관계가 존재하는 것을 확인할 수 있습니다. alpha = .05는 투명도를 의미하는 명령어로 관측치들의 집중도를 확인할 수 있습니다. 이 플롯을 통해 우리는 carat 값이 커질수록 관측치가 옅어지는 것을 확인할 수 있습니다.

그럼 이번에는 다이아몬드 데이터에 관측치가 너무 많기 때문에 딱 1,000의 배수가 되는 관측치들만 뽑아서 살펴보겠습니다—천번째, 이천번째, 삼천번째…

temp_df <- df |> rowid_to_column() |>
  dplyr::filter(rowid %% 1000 == 0)
glimpse(temp_df)
Rows: 53
Columns: 11
$ rowid   <int> 1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000, 9000, 10000, 1…
$ carat   <dbl> 1.12, 0.72, 0.90, 0.74, 1.05, 1.01, 1.00, 1.01, 0.91, 1.00, 0.…
$ cut     <ord> Premium, Ideal, Good, Ideal, Very Good, Premium, Premium, Prem…
$ color   <ord> J, D, G, D, I, H, F, G, E, D, D, D, G, H, F, I, E, H, J, E, E,…
$ clarity <ord> SI2, VS2, SI2, SI1, SI2, SI2, SI2, SI2, SI1, SI1, VS2, SI1, SI…
$ depth   <dbl> 60.6, 61.8, 63.8, 61.8, 62.3, 62.8, 62.6, 61.9, 63.5, 64.9, 62…
$ table   <dbl> 59, 57, 56, 56, 59, 61, 58, 59, 57, 59, 59, 62, 60, 59, 61, 56…
$ price   <int> 2898, 3099, 3303, 3517, 3742, 3959, 4155, 4327, 4512, 4704, 49…
$ x       <dbl> 6.68, 5.76, 6.13, 5.84, 6.42, 6.33, 6.37, 6.46, 6.11, 6.20, 6.…
$ y       <dbl> 6.61, 5.73, 6.16, 5.88, 6.46, 6.25, 6.31, 6.39, 6.14, 6.13, 6.…
$ z       <dbl> 4.03, 3.55, 3.92, 3.62, 4.01, 3.95, 3.97, 3.98, 3.89, 4.00, 3.…
model <- lm(price ~ carat, data = temp_df)

summary(model)

Call:
lm(formula = price ~ carat, data = temp_df)

Residuals:
    Min      1Q  Median      3Q     Max 
-6874.0  -727.7  -413.6   121.2  5785.7 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   -836.3      471.6  -1.774   0.0821 .  
carat         5719.0      440.7  12.977   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1983 on 51 degrees of freedom
Multiple R-squared:  0.7676,    Adjusted R-squared:  0.763 
F-statistic: 168.4 on 1 and 51 DF,  p-value: < 2.2e-16

조금 세련된 방법으로 {broom} 패키지의 함수를 이용해 선형모델로부터 예측값과 그 외 기타 통계치들을 추출해보도록 하겠습니다.

fitted_values <-  broom::augment(model) 
head(fitted_values)
# A tibble: 6 × 8
  price carat .fitted .resid   .hat .sigma   .cooksd .std.resid
  <int> <dbl>   <dbl>  <dbl>  <dbl>  <dbl>     <dbl>      <dbl>
1  2898  1.12   5569. -2671. 0.0219  1966. 0.0207       -1.36  
2  3099  0.72   3281.  -182. 0.0200  2003. 0.0000882    -0.0929
3  3303  0.9    4311. -1008. 0.0189  1998. 0.00254      -0.513 
4  3517  0.74   3396.   121. 0.0197  2003. 0.0000384     0.0617
5  3742  1.05   5169. -1427. 0.0204  1993. 0.00550      -0.727 
6  3959  1.01   4940.  -981. 0.0198  1998. 0.00252      -0.500 
fitted_values |>  
  ggplot(aes(x = carat, y = price)) + 
  geom_point(color = futurevisions("mars")[4]) +                      
  geom_smooth(method = 'lm', se = F, color = futurevisions("mars")[1]) + 
  geom_segment(
    aes(x = carat, y = price, xend = carat, yend = .fitted),
    color = futurevisions("mars")[1], size = 0.3) +
  geom_curve(
    aes(x = 2.75, y = 10000, xend = 2.55, yend = 15000), 
    color = futurevisions("mars")[1],
    size = 1, arrow = arrow(length = unit(.1, "inches"))) +
  annotate(
    "text", x = 2.75, y = 9500, label = "Residual / Error", 
    color = futurevisions("mars")[1]
  ) +
  geom_curve(
    aes(x = 1.5, y = 17500, xend = 2.0, yend = 11000), 
    color = futurevisions("mars")[3],
    size = 1, arrow = arrow(length = unit(.1, "inches")),
    curvature = -0.05
  ) +
  annotate(
    "text", x = 1.5, y = 18500, 
    label = "Regression line / Fitted Values", 
    color = futurevisions("mars")[3]
  ) +
  labs(x = 'Carat', y = "Price") +
  scale_y_continuous(labels = scales::dollar_format())

이 그림은 carat과 price 간의 관계를 보여줄 때, 선형회귀선을 긋는다는 것이 어떠한 의미인지를 보여줍니다.

  • 우리가 가진 관측치는 개별 점들입니다.

  • 우리는 각 점들로부터 얻은 정보를 바탕으로 그 추세를 파악하기 위한 예측을 선형관계로 보여주게 됩니다.

    • 이때 중요한 것이, 어떤 기준으로 선을 그어야 하냐는 것입니다.

    • 선형회귀, 잔차의 제곱합을 사용할 경우는 각 점들로부터 수직거리(즉, 실제 관측치와 예측값 간의 차이), 잔차(residuals)의 제곱 합이 가장 적은 선을 그리도록 합니다.

      • 모집단에 대해서는 실제 관측값과 예측값의 차이는 오차(errors)지만, 표본에 적용할 때에는 잔차(residuals)라고 표현합니다.

    • 이러한 방법을 일반화 최소자승법(Ordinary least squares), OLS라고도 하는데, 결국 잔차의 제곱합이 최소가 되는 선을 긋는 방법을 말합니다.

    • 제곱을 하는 이유는 가정 상 잔차의 총합은 0 \((\sum_{i=1}^n\epsilon_i = 0)\)이라 절대값 혹은 제곱값을 구해서 더해줘야 모델의 잔차를 추정할 수 있기 때문인데, 절대값은 수리적으로 더 계산이 까다롭기 때문입니다.

Dealing with Qualitative Explanatory Variables

다양한 변수들 (Various variables)

질적 변수(Qualitative variables)

제3주차에서 변수의 종류에 대해서 다루어본 적이 있습니다. 간략하게 말하면 일정한 간격을 가진 연속적인 수로 이루어진 연속형 변수(continuous variables)와 각 부류 간에 서로 배타적인 명목형 변수(혹은 분류형 변수, nominal or categorical variables), 그리고 그 사이에 서로 다른 부류임에도 불구하고 순위를 매길 수 있는 순위형 변수(ordinal variables)라고 구분할 수 있습니다.3

일단 변수가 연속형—양적 변수라고 할 때, 우리는 설명변수 \(x\)가 \(\mathbb{R}\), 실수에 속한다고 할 수 있습니다. 즉 \(x\in\mathbb{R}\)이라고 한다면, \(y = \alpha + \beta x + u\)에서 \(\beta\)가 의미하는 것은 무엇일까요?

한편, 변수가 질적 변수—예를 들어, 명목형이라고 하겠습니다. 따라서 이때 새로운 변수 \(w\)는 어떤 분류목에 속한다는 의미에서 \(w\in\{A, B, C\}\)라고 하겠습니다. 이때의 \(z = \gamma + \delta w + \nu\)에서 \(\delta\)가 의미하는 것은 무엇일까요?

위의 두 모델은 같은 형태를 가지지만 설명변수의 유형이 다릅니다. 양적 변수는 \(x \in \mathbb{R}\)이었고, 질적 변수는 \(w \in \{A,B,C\}\)로 나타낼 수 있었습니다. 바로 위에서는 모델의 특정(specification)에 대하여 살펴보았다면, 이번에는 설명변수의 유형에 따라 우리가 얻는 \(\beta_i\)의 의미가 어떻게 달라지는지를 살펴보고자 하는 것입니다.

질적 변수로 우리가 모델에 종종 포함하는 것은 이른바 더미변수(dummy variables, dummies)라고 불리는 것들입니다. 더미변수는 분류형 변수를 각 카테고리별로 쪼개어 각각 변수로 만든 결과로서, 그 카테고리에 속할 경우 1, 속하지 않을 경우 0으로 보여줍니다. 더미변수라고 부르는 이유는 이 변수가 해당 관측치가 특정한 카테고리에 속해있다 혹은 속해있지 않다는 것만을 말해주나는 점에서 그 정보량이 연속형이나 다른 유형의 변수들에 비해 상대적으로 적다는 것 때문입니다. \(w\) 변수를 예로 들자면, 우리는 \(w\)를 각각 \(\{A,B,C\}\)라는 카테고리별로 더미변수 3개로 쪼갤 수 있습니다. \[\begin{equation*} w = \begin{cases} A & \text{iff } wA = 1\\ B & \text{iff } wB = 1\\ C & \text{iff } wC = 1 \end{cases} \end{equation*}\]

그렇다면 이 더미변수들을 질적 변수로 구성했던 \(z = \gamma + \delta w + v\)에 투입하여 이 회귀모델을 실제 분석가능한 모델로 재구성해보도록 하겠습니다. \[\begin{equation*} \begin{aligned} &z = \gamma_1 I(w = A) + \gamma_2 I(w = B) + \gamma_3 I(w = C) + \nu\\ &z = \gamma_1 wA + \gamma_2 wB + \gamma_3 wC + \nu \end{aligned} \end{equation*}\] 자, 이렇게 분류형 변수를 더미변수들로 나누어 모델에 집어넣어 보면, 이론적으로는 문제가 없어보입니다. 왜냐하면 각각의 더미변수들을 모두 합쳐놓은 것이 원래의 분류형 변수일테니까요. 그럼 이대로 분석이 가능할까요? 아닙니다. 왜냐하면 위의 회귀분석모델은 완전다중공선성(Perfect Mulitcollinearity)로 인해 분석이 이루어지지 않게 됩니다. 세 더미변수의 총합, 즉 절편값이 1이 되기 때문입니다(\(wA+wB+wC = 1\)). 따라서 모든 조건이 일정하다고 할 때, \(z = \gamma + \delta_1 wA + \delta_2 wB + \delta_3 wC + \nu\)는 사실상 아무 것도 분석하지 못합니다.

더미변수들을 가지고 분석모델을 구성할 때는 전체 분류에 속하는 더미변수들 중 하나를 제외하고 모델에 투입해야 합니다. 그 제외된 분류가 바로 전체 더미들에 대한 기준변수가 됩니다. 예를 들어보면, 종속변수에 대한 계절의 효과를 보고 싶다고 합시다. 만약 봄, 여름, 가을, 겨울을 각각 더미로 측정할 경우 그 중 봄을 제외한 나머지 세 계절의 더미변수를 모델에 투입합니다. 이때, 여름/가을/겨울의 계수는 봄이라는 계절에 비교하여 해석할 수 있습니다. 일단 수리적 모형으로 살펴보겠습니다. 이번에는 사계절이 아니라 앞서의 \(w \in \{A,B,C\}\)의 분류형 변수를 기준으로 보겠습니다. \[\begin{equation*} \begin{aligned} &z = \gamma + \delta_2 wB + \delta_3 wC + v\\ &z = \gamma + \delta_1 wA + \delta_3 wC + v\\ &z = \gamma + \delta_1 wA + \delta_2 wB + v\\ \end{aligned} \end{equation*}\] 첫 번째 식의 경우 \(\delta_2, \delta_3\)는 \(wA\)에 비하여 \(wB\)와 \(wC\)가 종속변수에 미치는 효과를 보여줍니다. \(\delta_2 < 0\)이라면 우리는 \(wB\)가 \(wA\)에 비해 종속변수에 미치는 효과가 \(\delta_2\)만큼 ‘덜’ 하다는 의미입니다. 간략히 말하자면, 더미변수들의 효과는 항상 생략된 카테고리에 비교하여 해석되어야 합니다.

한번 R로 데이터를 통해 살펴보겠습니다. 사용할 데이터는 언제나와 같이 Quality of Government에서 가져오며 2016년 데이터입니다.

library(ezpickr)
library(tidyverse)
QOG <- ezpickr::pick(file = "http://www.qogdata.pol.gu.se/data/qog_bas_ts_jan19.dta")

QOG.s <-QOG |> 
  dplyr::select(ccode, cname, year, 
                wdi_agedr, wdi_trade, 
                wdi_gdpcapcon2010, ht_region) |>
  dplyr::filter(year==2016) |> drop_na()

rm(QOG)

QOG 데이터셋에서 국가코드, 국가명, 연도, 무역 개방성, 1인당 GDP, 지역에 해당하는 변수들을 따로 선별하여 서브셋을 만들고, 결측치를 제외하였습니다. 그리고 로그값을 취한 1인당 GDP가 무역 개방성 간의 관계를 보여주는 단순선형회귀모델 하나, 로그값을 취한 1인당 GDP와 지역 간 관계를 보여주는 단순선형회귀모델 하나를 각각 구축하였습니다.

QOG.s |> mutate(
  reg = case_when(
    ht_region == 1 ~ 
      "1. Eastern Europe and post Soviet Union",
    ht_region == 2 ~ 
      "2. Latin America",
    ht_region == 3 ~ 
      "3. North Africa & the Middle East",
    ht_region == 4 ~ 
      "4. Sub-Saharan Africa",
    ht_region == 5 ~ 
      "5. Western Europe and North America",
    ht_region == 6 ~ 
      "6. East Asia",
    ht_region == 7 ~ 
      "7. South-East Asia",
    ht_region == 8 ~ 
      "8. South Asia",
    ht_region == 9 ~ 
      "9. The Pacific",
    ht_region == 10 ~ 
      "10. The Caribbean",
    T ~ NA_character_),
  reg = factor(reg, levels = c(
       "1. Eastern Europe and post Soviet Union",
       "2. Latin America",
       "3. North Africa & the Middle East",
       "4. Sub-Saharan Africa",
       "5. Western Europe and North America",
       "6. East Asia",
       "7. South-East Asia",
       "8. South Asia",
       "9. The Pacific",
       "10. The Caribbean"))) -> QOG.s

table(QOG.s$reg) |> 
  knitr::kable(col.names = c("Region", "Counts"))
Region Counts
1. Eastern Europe and post Soviet Union 27
2. Latin America 18
3. North Africa & the Middle East 18
4. Sub-Saharan Africa 43
5. Western Europe and North America 23
6. East Asia 4
7. South-East Asia 11
8. South Asia 8
9. The Pacific 5
10. The Caribbean 9 
model_ch7_1 <- lm(
  log(wdi_gdpcapcon2010) ~
    wdi_trade, data = QOG.s)

model_ch7_2 <- lm(
  log(wdi_gdpcapcon2010) ~
    reg, data = QOG.s)

texreg::htmlreg(list(model_ch7_1, model_ch7_2),
                single.row = T, digits = 3)
Statistical models
  Model 1 Model 2
(Intercept) 7.700 (0.218)*** 8.848 (0.181)***
wdi_trade 0.010 (0.002)***  
reg2. Latin America   -0.265 (0.287)
reg3. North Africa & the Middle East   0.346 (0.287)
reg4. Sub-Saharan Africa   -1.743 (0.231)***
reg5. Western Europe and North America   1.911 (0.267)***
reg6. East Asia   0.658 (0.505)
reg7. South-East Asia   -0.436 (0.337)
reg8. South Asia   -1.385 (0.379)***
reg9. The Pacific   -0.809 (0.459)
reg10. The Caribbean   0.190 (0.363)
R2 0.116 0.623
Adj. R2 0.110 0.601
Num. obs. 166 166
***p < 0.001; **p < 0.01; *p < 0.05

연속형 변수와 질적 변수가 선형회귀모델에 포함되었을 때, 어떠한 차이가 있는지 살펴보겠습니다. 먼저, 첫 번째 모델은 다음과 같이 표현할 수 있습니다.

\[ \ln(\mathrm{GDPpc}) = \beta_0 + \beta_1\times \text{Trade openness} + \epsilon. \]

첫 번째 모델을 해석하자면,

  • 다른 조건들이 모두 일정할 때, 무역개방성이 한 단위 증가할 때, 1인당 GDP는 0.01%p 증가한다.

  • 그리고 이 모델은 로그값을 취한 1인당 GDP의 약 12%를 설명한다.

정도로 요약할 수 있습니다. 이 관계를 시각화하면 다음과 같습니다.

tidy_model1 <- model_ch7_1 |> broom::augment()
head(tidy_model1)
# A tibble: 6 × 8
  `log(wdi_gdpcapcon2010)` wdi_trade .fitted  .resid    .hat .sigma     .cooksd
                     <dbl>     <dbl>   <dbl>   <dbl>   <dbl>  <dbl>       <dbl>
1                     6.43      55.9    8.27 -1.85   0.00821   1.40 0.00718    
2                     8.45      74.7    8.46 -0.0111 0.00633   1.41 0.000000201
3                     8.48      56.3    8.28  0.207  0.00816   1.41 0.0000899  
4                     8.18      59.4    8.31 -0.125  0.00773   1.41 0.0000311  
5                     9.50      88.7    8.61  0.890  0.00605   1.41 0.00122    
6                     8.68      90.1    8.62  0.0557 0.00607   1.41 0.00000482 
# ℹ 1 more variable: .std.resid <dbl>
tidy_model1 |> 
  ggplot(aes(x = wdi_trade, y = `log(wdi_gdpcapcon2010)`)) + 
  geom_point(size = 2, alpha = 0.75) + 
  geom_smooth(method = "lm", se = F) + 
  labs(x = "Trade openness", y = "Logged GDPpc") + 
  theme_bw()

두 번째 모델은 다음과 같이 표현할 수 있습니다.

\[ \begin{aligned} \operatorname{log(wdi\_gdpcapcon2010)} &= \alpha + \beta_{1}(\operatorname{reg}_{\operatorname{2.\ Latin\ America}})\\ &\quad + \beta_{2}(\operatorname{reg}_{\operatorname{3.\ North\ Africa\ \&\ the\ Middle\ East}}) + \beta_{3}(\operatorname{reg}_{\operatorname{4.\ Sub-Saharan\ Africa}})\\ &\quad + \beta_{4}(\operatorname{reg}_{\operatorname{5.\ Western\ Europe\ and\ North\ America}}) + \beta_{5}(\operatorname{reg}_{\operatorname{6.\ East\ Asia}})\\ &\quad + \beta_{6}(\operatorname{reg}_{\operatorname{7.\ South-East\ Asia}}) + \beta_{7}(\operatorname{reg}_{\operatorname{8.\ South\ Asia}})\\ &\quad + \beta_{8}(\operatorname{reg}_{\operatorname{9.\ The\ Pacific}}) + \beta_{9}(\operatorname{reg}_{\operatorname{10.\ The\ Caribbean}})\\ &\quad + \epsilon \end{aligned} \]

두 번째 모델은 어떻게 해석할까요? 첫 번째와는 조금 다릅니다. 왜냐하면 일단 “동유럽과 구소비에트 연방(Eastern Europe and post Soviet Union)”에 해당하는 카테고리가 기준변수(reference variable)로 다른 변수들은 그 기준변수와 비교한 효과를 보여주기 때문입니다. 단적으로 편미분을 취한다고 했을 때, 라틴 아메리카 지역이 로그값을 취한 1인당 GDP에 미치는 효과를 보고 싶다면 다음과 같이 나타납니다.

\[ \frac{\partial \ln(GDPpc)}{\partial\text{Latin America}} = \beta_1 \] 사실상 \(\beta_1\)은 다음과 같은 식의 결과라고 할 수 있습니다.

\[ \begin{aligned} \ln(\mathrm{GDPpc})\: &=\beta_0 + \beta_1\times \text{Latin America = 1}\\ & + \beta_2\times 0 \\ & + \beta_3\times 0 \\ & + \beta_4\times 0 \\ & + \beta_5\times 0 + \beta_6\times 0\\ & + \beta_8\times 0 + \beta_9\times 0\\ & + \beta_{10}\times 0 + \epsilon. \end{aligned} \]

그렇다면 \(\beta_0\)은 무엇일까요? 네. 나머지 모든 지역이 0인 경우, 즉 기준변수인 동유럽과 구소비에트 연방 지역의 평균 1인당 GDP 값을 보여주는 것입니다. 즉, 질적 변수–더미변수를 모델에 집어넣으면 그 결과는 연속형 변수들과 같은 기울기의 차이를 보여주는 것이 아니라 기울기는 같지만 서로 다른 “절편”의 값을 보여준다고 할 수 있습니다. 라틴 아메리카 지역일 경우 평균 1인당 GDP 값을 모델에서 어떻게 확인할까요? 동유럽 지역일 경우 \(\beta_0\)이고 라틴 아메리카를 포함할 경우 \(\beta_0 + \beta_1\)이 되니 이 둘의 차이를 구하면 라틴 아메리카의 평균 1인당 GDP의 값을 구하게 되는 것입니다. 따라서 해석은 다음과 같습니다.

  • 다른 조건들이 모두 일정할 때, 라틴 아메리카 지역일 경우 1인당 GDP는 기준변수인 동유럽 및 구소비에트 연방지역에 비하여 평균적으로 -26%p 낮다.

  • 그리고 이 모델은 로그값을 취한 1인당 GDP의 약 62%를 설명한다.

정도로 요약할 수 있습니다. 이 관계를 시각화하면 다음과 같습니다.

tidy_model2 <- model_ch7_2 |> broom::augment()
head(tidy_model2)
# A tibble: 6 × 8
  `log(wdi_gdpcapcon2010)` reg   .fitted .resid   .hat .sigma .cooksd .std.resid
                     <dbl> <fct>   <dbl>  <dbl>  <dbl>  <dbl>   <dbl>      <dbl>
1                     6.43 8. S…    7.46 -1.04  0.125   0.941 1.98e-2     -1.18 
2                     8.45 1. E…    8.85 -0.397 0.0370  0.945 7.07e-4     -0.429
3                     8.48 3. N…    9.19 -0.712 0.0556  0.944 3.55e-3     -0.777
4                     8.18 4. S…    7.11  1.08  0.0233  0.941 3.18e-3      1.16 
5                     9.50 10. …    9.04  0.459 0.111   0.945 3.33e-3      0.516
6                     8.68 1. E…    8.85 -0.172 0.0370  0.945 1.33e-4     -0.186
tidy_model2 |> 
  ggplot(aes(x = reorder(reg, `log(wdi_gdpcapcon2010)`), 
             y = `log(wdi_gdpcapcon2010)`,
             color = reg,
             fill = reg)) + 
  geom_jitter(size = 2, alpha = 0.75, show.legend = F) + 
  geom_boxplot(alpha = 0.2, show.legend = F) + 
  labs(x = "Region", y = "Logged GDPpc") + 
  coord_flip() + 
  theme_bw() + theme(legend.title = element_blank(),
                     legend.position = "bottom")

아직 배우지는 않았지만 다중회귀분석을 통해서 보면 더 확연하게 살펴볼 수 있습니다.

model_ch7_3 <- lm(
  log(wdi_gdpcapcon2010) ~
    reg + wdi_trade, data = QOG.s)

tidy_model3 <- model_ch7_3 |> broom::augment()
head(tidy_model3)
# A tibble: 6 × 9
  `log(wdi_gdpcapcon2010)` reg   wdi_trade .fitted  .resid   .hat .sigma .cooksd
                     <dbl> <fct>     <dbl>   <dbl>   <dbl>  <dbl>  <dbl>   <dbl>
1                     6.43 8. S…      55.9    7.43 -1.01   0.125   0.905 1.83e-2
2                     8.45 1. E…      74.7    8.67 -0.213  0.0400  0.909 2.19e-4
3                     8.48 3. N…      56.3    9.06 -0.575  0.0572  0.908 2.36e-3
4                     8.18 4. S…      59.4    7.03  1.15   0.0237  0.904 3.64e-3
5                     9.50 10. …      88.7    9.01  0.482  0.111   0.908 3.62e-3
6                     8.68 1. E…      90.1    8.75 -0.0767 0.0378  0.909 2.66e-5
# ℹ 1 more variable: .std.resid <dbl>
library(RColorBrewer)
# Classic palette BuPu, with 4 colors
coul <- brewer.pal(4, "PuOr") 
# Add more colors to this palette :
coul <- colorRampPalette(coul)(10)

tidy_model3 |> 
  ggplot(aes(x = wdi_trade, y = `log(wdi_gdpcapcon2010)`,
             color = reg, fill = reg,
             group = reg, label = reg)) + 
  geom_point(size = 2, alpha = 0.75) + 
  geom_smooth(method = "lm", se = F, 
              aes(y = predict(model_ch7_3, QOG.s)),
              show.legend = F, fullrange = T) + 
  scale_color_manual(values =coul) + 
  labs(x = "Trade openness", y = "Logged GDPpc") + 
  theme_bw() + 
  theme(legend.title = element_blank(),
        legend.position = "bottom") + 
  guides(color=guide_legend(ncol=2),
         fill=guide_legend(ncol=2))

변수에 대한 심화

우리가 \(y_i = \delta_0 + \delta_1 x_i\)인 모델을 가지고, 이때 \(i \leq n\)이라고 합시다. 그리고 \(a = \frac{1}{n}\sum_i x_i\)이고 \(b=2sd(x)\)라고 하겠습니다. 이때, \(y_i = \lambda_0 + \lambda_1\frac{x_i -a}{b}\)를 어떻게 추정할 수 있을까요? Gelman (2008) 의 내용에서 표준오차에 대해 조금 더 자세히 이해하기 위해 가져와 본 내용입니다.

단순하게 보겠습니다. 만약 \(b = 2sd(x)\)가 아니라 \(b=sd(x)\)였다면, 두 번째 수식에서 \(\frac{x_i-a}{b}\)는 각 관측치들에서 평균을 제외하고 그것을 표준편차로 나누는, 일종의 표준화(standardization) 작업이라고 이해할 수 있을 것입니다. 표준화 작업은 우리로 하여금 변수의 측정척도의 영향에서 자유롭게 합니다. 그런데 만약 \(b=2sd(x)\)라고 한다면, 이 작업이 근본적으로 무언가 차이가 있을까요? 왜 우리는 꼭 \(sd(x)\)로 표준화를 해줘야만 하는 것일까요?

\(y = \delta_0 + \delta_1 x_i\)는 단순회귀모델입니다. 이때 \(\delta_0\)는 절편값이고 \(\delta_1\)은 모델의 회귀계수, 즉 \(x_i\)가 \(y\)와 맺는 관계 양상을 보여주는 것이죠. 그런데 우리가 두 번째로 수립한 수식에서 \(a\)는 \(x\)의 평균입니다. \(b\)는 \(x\)의 표준편차에 2를 곱한 것이죠. 따라서 \(a\), \(b\)로 되어 있는 이 식을 \(x\)와 관련하여 다시 구성하여주면 다음과 같습니다: \(y = \lambda_0 + \lambda_1(\frac{x_i - \bar{x}}{2\times\sigma_x})\). 앞서 말했다시피 일종의 표준화 작업입니다. 이론적으로 \(sd(x)\)와 \(2sd(x)\)의 차이를 논의하기 전에 실제 데이터를 가지고 어떠한 편화가 있는지를 살펴보겠습니다. 이미 만들어놓은 QOG.s 데이터셋의 1인당 GDP와 무역 개방성 변수를 이용해보겠습니다.4

QOG.s |> select(wdi_trade) |> mutate(
  wdi_trade = wdi_trade,
  wdi_trade1sd = (wdi_trade - mean(wdi_trade))/sd(wdi_trade),
  wdi_trade2sd = (wdi_trade - mean(wdi_trade))/2*sd(wdi_trade)
) |> summary() |> knitr::kable()
wdi_trade wdi_trade1sd wdi_trade2sd
Min. : 20.72 Min. :-1.3093 Min. :-1616.3
1st Qu.: 56.10 1st Qu.:-0.5973 1st Qu.: -737.3
Median : 77.26 Median :-0.1714 Median : -211.6
Mean : 85.78 Mean : 0.0000 Mean : 0.0
3rd Qu.:102.75 3rd Qu.: 0.3415 3rd Qu.: 421.5
Max. :407.43 Max. : 6.4735 Max. : 7991.1

그럼 이번에는 세 변수 각각으로 종속변수를 선형회귀분석으로 추정했을 때의 결과를 비교해보도록 하겠습니다. 일단 이 예제를 진행하는 데 있어서 기술적으로(technically) 조작한(manipulated) 변수들의 내용은 다음과 같습니다.

  • 첫째, 1인당 GDP는 달러 기준이므로 단위가 너무 커서 가시적으로 보여주기가 어렵기 때문에 로그를 취한 값을 사용했습니다. 무역 개방성은 비율이기 때문에 단위차이가 너무 크면 표준화한 결과와 원변수의 결과를 한 그래프 안에서 보여주기 어려울 테니까요.

  • 둘째, 무역 개방성은 원변수(Original), 1표준편차로 표준화한 값(One_std), 그리고 2표준편차로 표준화한 값(Two_std)로 조작하였습니다.

  • 절편값에는 크게 관심이 없으므로 제외하고 무역 개방성과 1인당 GDP 간의 선형회귀모델 결과만 보여드리도록 하겠습니다.

간단하게 정리하자면 예측변수의 통계적 유의성 여부, 방향성은 표준화 여부를 떠나서 같습니다. 그리고 여기서는 보여드리지 않았지만 \(R^2\)도 동일합니다. 그 얘기는 세 모델 모두 동일한 \(y\)의 변동량의 일부를 설명하고 있다는 것이겠죠? 그러나 \(\hat{\beta_x}, \hat{\beta_x^{\sigma}}, \hat{\beta_x^{2\sigma}}\)와 \(se_{\hat{\beta_x}}, se_{\hat{\beta_x^{\sigma}}}, se_{\hat{\beta_x^{2\sigma}}}\)는 서로 다릅니다.

만약 우리가 변수를 적절하게 표준화한다면, \(x\)와 \(y\)의 관계는 왜곡되지 않게 나타날 것입니다. 또한 변수들은 변수들의 평균에 영향을 받지만 동시에 변수들의 분산—표준편차에도 영향을 받습니다. 왜냐하면 회귀모델은 통계적 기법으로 얘기하자면 이른바 “평균 차이”를 검정하는 방법이기 때문입니다.5 그러나 변수를 표준화한 다음에는 \(x\)의 한 단위 변화라고 해석할 때, 그 한 단위가 무엇을 의미하는지 고민해볼 필요가 있습니다.

일반적인 표준화에 대한 함의가 위의 논의였다면, 다음은 Gelman (2008) 의 제안대로 \(2sd(x)\)를 이용해 표준화한 결과가 무슨 의미가 있는지 살펴보겠습니다. 일반 표준화와 비교해서 무슨 차이가 있나요? 보시면 그래프에서 원변수의 계수의 양상과 \(2sd(x)\) 표준화된 결과의 모습이 유사한 것을 확인할 수 있습니다. 즉, 측정척도의 영향에서 자유로우면서도 상대적으로 원변수의 경향성과 비슷한 모습을 보여줌으로써 우리는 단순 표준화 변수를 통한 결과보다 직관적인 이해가 가능하다고 할 수 있겠습니다.

다중선형회귀모델과 BLUE

단순선형회귀모델은 단순해보이지만 결코 단순하지 않습니다. 앞에서 우리는 +가 각 변수들 간의 독립적인 관계를 보여준다는 것을 확인하였습니다. 이번에는 범죄에 대한 회귀모델을 구축하되, 임금 수준 이외의 또 다른 변수—경찰의 노력(Police effort)를 포함해봅시다. 선형회귀모델이되 둘 이상의 설명변수를 포함하였으니 다중선형회귀모델(multiple linear regression model)이라고 할 수 있을 것입니다. 통계적 모델로 나타내자면 \(\text{범죄} = \beta_0 + \beta_1\text{임금 수준} + \beta_2\text{경찰의 노력} + u\)가 될 것입니다.

우리는 일단 다중선형회귀모델의 각 계수, \(\beta_i\)들을 어떻게 구할 수 있는지를 살펴볼 것입니다. 그리고 설명변수로서의 제곱항(square term)이 포함된 다중선형회귀모델에서 \(\beta_i\)가 가지는 의미를 분석하고자 합니다.다중선형회귀모델 통계적 특성은 다음과 같습니다.

  • MLR.1: 모수들은 선형관계를 이루고 있다.

  • MLR.2: 표본은 무작위로 추출된 것이다(random sampling)

  • MLR.3: 완전한 다중공선성은 존재하지 않아야 한다(No perfect collinearity)

  • MLR.4: 오차항과 설명변수들 간의 상관관계가 존재하지 않아야 한다(\(E(u|x_i,\cdots,x_k) = 0\)).

우리는 MLR.1부터 MLR.4까지를 가정하고, \(y_i = \hat{\beta_0} + \hat{\beta_1}x_{i1} + \hat{\beta_2}x_{i2}\)를 가정합니다. 따라서 다중선형회귀모델에서 \(\hat{\beta_1}\)을 구하기 위해서는 다음과 같은 공식이 성립합니다. \[\hat{\beta_1} = \frac{\sum^{n}_{i=1}\hat{r}_{i1}\cdot y_i}{\sum^{n}_{i=1}\hat{r}^2_{i1}}\] 이때, \(r_{i1}\)은 잔차로 \(x_{i1} = \hat{\alpha}_0 + \hat{\alpha}_1x_{i2} + r_{i1}\)에서 도출된 것입니다. 간단히 말하자면, 단순선형회귀모델에서와는 달리 다중선형회귀모델에서는 \(\hat{\beta_1}\)을 구할 때, \(x_1\)과 \(y\) 간의 관계뿐 아니라 \(x_1\)과 \(x_2\) 간의 관계, 그리고 \(x_2\)와 \(y\)의 관계도 고려해주어야 한다는 것입니다. 왜? \(x_1\)이 설명하고 있는 부분의 일부는 사실 \(x_1\)과 \(x_2\)가 같이 설명하는 “교집합”, “공분산”(covariance)일 수 있기 때문입니다.

공분산에 대한 내용은 뒤에서 더 자세하게 다루기로 하고, 우선 MLR.1부터 MLR.4까지의 가정들이 모두 충족되었을 때, 우리는 수없이 많은 표본들로부터 얻어내는 \(\hat{\beta}_k\)의 기대값이 모집단에서의 \(\beta_k\)와 같을 것이라고 생각할 수 있습니다: \(E(\hat{\beta}_k) = \beta_k\). 이를 최소자승법으로 구한 선형회귀모델의 계수값의 비편향성(unbiasedness)라고 합니다.

오차(Errors)

이번에는 오차에 대해서 배워보겠습니다. 앞서 우리는 모집단 수준에서의 오차의 기대값은 0으로 수렴한다(\(E(u) = 0\))는 것을 가정하였습니다. 그렇다면 오차의 분산(\(Var(u)\))은 어떨까요? 이 오차의 분산에 관한 내용이 바로 다중선형회귀모델의 다섯번째 통계적 특성이자 가정, MLR.5입니다.

  • MLR.5: 오차의 분산은 동질적이다(Homosckedasticity).

    • \(Var(u|x_i, \cdots, x_k)=\sigma^2\)

    • \(E(y|x_1, x_2) = \beta_0 + \beta_1x_1 + \beta_2x_2\)

    • \(Var(y|x_1, x_2) = Var(y) = Var(u) = \sigma^2 = E(u^2)\)

    • 오차의 표준편차는 \(\sqrt{\sigma^2} = \sigma\)가 됩니다.

  • 위에서부터 순서대로 저 가정의 내용을 풀어가보면, 먼저 우리는 주어진 설명변수들 하에서 오차항의 분산을 모집단 수준에서 \(\sigma^2\)라고 합니다. \(\sigma\)가 표준편차고 분산은 표준편차의 제곱이라는 수리적 정리를 여기서 딱히 증명할 필요는 없을 거 같으니 이렇게 넘어가겠습니다.

  • 그리고 \(x_1, x_2\)라는 다중선형회귀모델의 두 설명변수가 있다고 할 때, 그 설명변수들이 주어졌을 때의 종속변수를 예측할 수 있는 기대값은 오차를 제외한 두 설명변수의 다중선형회귀모델의 PRF로 결정됩니다.

  • 그리고 이러한 설명변수들이 주어졌을 때의 종속변수의 분산은 MLR.5에 따라 오차항의 분산과 같고, 오차항의 분산은 \(\sigma^2\)입니다. 이는 \(u^2\)의 기대값과도 같습니다.

  • 표준편차는 분산의 제곱근이므로 결과적으로 오차의 표준편차는 \(\sigma\)가 됩니다.

\(\hat{\sigma}^2\) 추정하기

표본에서의 오차—즉, 잔차의 분산(\(\hat{\sigma}^2\))은 표본의 크기의 영향을 받습니다. 정확히는 \[\hat{\sigma}^2=\frac{1}{n-K-1}\sum\hat{u}^2_i\]이며, 이는 곧 잔차의 제곱합을 n-k-1, 전체 관측치의 개수에서 변수의 개수-1을 하는 것과 같습니다(\(SSR/(n-K-1)\)).

MLR.1부터 MLR.5까지의 가정이 모두 충족된다면 모집단 오차의 분산, \(\sigma^2\)에 대한 불편추정량인 \(\hat{\sigma}^2\)을 얻을 수 있습니다.

계수(Coefficients)

최소자승법(ordinary least square)는 예측값과 실제 관측치 간의 차이인 잔차의 제곱합을 최소로 하는 \(\hat{\beta}_k\)를 구하는 방법을 말합니다. 따라서 최소자승법에 따라 구한 \(\beta_i\)가 편향되어 있지 않다는 것은 모집단에서 추출한 표본들로부터 구한 각 \(\hat{\beta}_k\)의 기대값이 모집단의 \(\beta_k\)와 같다는 것(\(E(\hat{\beta}_k)=\beta_k\))을 의미합니다.

그러나 \(\hat{\beta}_k\)는 하나의 추정치에 불과하므로 필연적으로 불확실성을 내포하고 있습니다. 따라서 표본의 통계치로 모집단의 모수를 추론할 때에는 단지 \(\hat{\beta}_k\)를 제시하는 것만으로는 부족합니. 예를 들어, 우리는 “모집단에서의 계수, \(\beta_k\)가 0보다 큰가?”라는 질문에 대답하기 위해서는 다음과 같은 것들을 필요로 합니다.

  • \(\beta_k\)에 대한 (최선의) 추정치인 \(\hat{\beta}_k\)

  • 우리가 그 최선의 추정치에 대해 얼마만큼 “신뢰할 수 있느냐”에 관한 \(\hat{\beta}_k\)의 분산(\(Var(\hat{\beta}_k)\)). 이 분산은 이론적으로 모집단에서 수많은 표본들을 뽑고 그 표본들에 대한 SRF에서 도출된 \(\hat{\beta}_k\)들이 얼마나 서로 다른지를 보여줍니다.

  • 또, \(\hat{\beta}_k\)에 대해 신뢰하기 위해서 t-statistic, p-값 등과 같은 통계치들을 확인하고는 합니다.

분산(Variance)

분산에 관한 내용들이 계속해서 나오는데, 과연 분산이 크다는 것과 작다는 것은 무엇을 의미하는 것일까요? 다중선형회귀모델에서 분산이 중요하게 논의되는 맥락을 이해하기 위해서는 먼저 표집분포(sampling distribution)에 대해 짚고 넘어갈 필요가 있습니다.

표집분포(Sampling distribution)

\(\hat{\beta}_k\)을 확률변수라고 생각해봅시다. 주어진 하나의 표본에 대해서 \(\hat{\beta}_k\)는 고정된 값입니다. 그 하나의 표본에서 \(\hat{\beta}_k\)는 정해져 있기 때문입니다. 그러나 사실 우리는 모집단으로부터 또 다른 표본을 확보할 수 있습니다. MLR.2의 무작위 표집(혹은 확률표집) 가정에 따라서 우리는 하나의 모집단에서 뽑아낸 수없이 많은 표본들로부터 \(\hat{\beta}_k\)들을 일종의 확률변수로서 분포로 보여줄 수 있게 됩니다. 우리는 이 표본들로부터의 얻어낸 \(\hat{\beta}_k\)들의 분포를 표집분포(sampling distribution)라고 합니다. 그리고 우리는 그 표집분포가 모집단에 대한 \(\beta_k\)의 기대값을 포함하고 있을 것이라고 기대하고, \(\hat{\beta}_k\)의 분산—\(Var(\hat{\beta}_k)\)이 그 분포가 모집단의 \(\beta_k\)을 포함할 “불확실성”을 보여주는 것입니다.

표집분포의 표준오차(Standard errors)

여기서 우리는 수리통계적으로 \(\hat{\beta}_k\)에 대한 표준편차의 추정치에 대한 정리를 생각해볼 수 있습니다. 즉, \(\hat{\beta}_k\)의 표준오차에 대한 추정치는 표본에 따라 조건적이며, 동시에 MLR.1과 MLR.5 가정에 기초하고 있습니다. \[se(\hat{\beta}_k) = \frac{\hat{\sigma}}{[SST_k\cdot (1-R^2_k)]^{1/2}}\] 이때, \(SST_k\)는 \(x_k\)의 변동량을 의미하며, \(R^2_k\)는 다른 모든 설명변수들로 \(x_k\)를 회귀분석한 \(R^2\) 결과라고 할 수 있습니다.

표집분포의 표준오차는 또 다른 방식으로도 보여줄 수 있는데, 수리적으로만 유도해보도록 하겠습니다. \[\begin{equation*} \begin{aligned} sd(\hat{\beta_k})&= \sqrt{Var(\hat{\beta_k})}\\ &= \frac{\sigma}{[SST_k\cdot(1-R^2_k)]^{1/2}}\\ se(\hat{\beta_k})&= \frac{\hat{\sigma}}{[SST_k\cdot(1-R^2_k)]^{1/2}}\\ &= \frac{\hat{\sigma}}{\sqrt{n}\cdot sd(x_k)\cdot \sqrt{1-R^2_k}} \end{aligned} \end{equation*}\]

즉, 표준오차가 작을수록 우리가 표본을 통해 얻어낸 표본의 \(\hat{\beta}_k\)들이 만들어내는 분포가 모집단의 \(\beta_k\)을 포함하고 있을 가능성이 높다고 할 수 있습니다.

Best Linear Unbiased Estimator

그렇다면 이제는 이른바 BLUE, 불편추정량에 대해 알아볼 차례입니다. 대체 이 불편추정량이라는 것이 무엇을 의미하는 것일까요? 바로 모든 표본들을 통틀어 ‘평균적으로’ 그 추정량이 최선의(효율적이고 편향되지 않은) 추정량이라는 것을 의미합니다. 그리고 표준오차가 작다는 것은 그만큼 우리의 추정량이 “효율적”(efficient)이라는 것을 말합니다. MLR.1부터 MLR.5까지의 가정이 충족되는 하에서 OLS 추정량은 불편추정량입니다.

모델의 특정(specification)

부적절한(irrelevant) 변수의 포함

만약 우리가 추정해야할 모델이 \(y = \beta_0 + \beta_1x + u\)라고 합시다. 그런데 모델을 \(y = \beta_0 + \beta_1x + \beta_2z + e\)로 수립해 추정하였다고 합시다. 이 경우에 부적절한 \(z\) 변수를 모델에 포함하여 모델을 잘못 특정하였다고 말합니다(misspecified).

잘못 특정된 모델: \(E(\hat{\beta_1})\)과 \(Var(\hat{\beta_1})\)에 미치는 영향

좀 더 자세히 이 잘못 특정된 모델이 \(E(\hat{\beta_1})\)에 미치는 효과를 살펴보도록 하겠습니다. 우리는 MLR.1로부터 MLR.4까지의 가정 하에서 \(E(\hat{\beta_k})=\beta_k\)가 된다는 사실을 알고 있습니다. 그렇다면 \(\beta_2\), 즉 잘못 집어넣은 이 변수의 계수값이 0이라면 그것이 MLR.1부터 MLR.4까지의 가정에 영향을 미칠까요?

  • 부적절한 변수를 모델에 포함하게 되면, 이 모델에서 우리는 \(\tilde{\beta_1}\)를 얻게 됩니다. \(\hat{\beta_1}\)이 제대로 된 모델에서 얻을 수 있는 OLS 추정치라고 합시다. 그러면 \(E(\tilde{\beta_1}|x_k)=\beta_1\)라고 표현할 수 있습니다.

  • 이 경우에 분산은 다음과 같은 관계를 가지게 됩니다.

\[\begin{equation*} \begin{aligned} Var(\hat{\beta_1}|x_k)&=\frac{\sigma^2}{\sum^N_{i=1}(x_{i1}-\bar{x_1})^2}\\ &\leq \frac{\sigma^2}{(1-R^2_1)\sum^N_{i=1}(X_{i1}-\bar{x_1})^2} = Var(\tilde{\beta_1}|x_k). \end{aligned} \end{equation*}\]

  • \(se(\hat{\beta_1}) = \frac{\sigma^2}{[SST_1\cdot(1-R^2_1)]^{1/2}}\)라는 것을 다시 한 번 떠올려봅시다.

    • 부적절한 변수를 모형에 포함시키더라도 \(SST_1\)는 변하지 않습니다.

    • 이때, \(R^2_1\)은 \(x_1\)을 종속변수로 하는 \(x_2\)의 회귀분석으로부터 얻어낸 \(R^2\)입니다.

      • 정확히는 \(x_1 = \alpha_0 + \alpha_1z + e\)

    • 만약 \(z\)가 \(x\)를 잘 설명한다면, \(R^2_1\)는 높아질 것이고, \((1-R^2_1)\)는 작아질 것입니다. 따라서 분모가 작아지므로 결과적으로 \(se(\hat{\beta_1})\)은 커지게 됩니다.

  • 따라서 표본에서 \(x_1\)과 \(x_2\)가 서로 상관되어 있다면, \(x_2\)를 포함하는 것은 \(\beta_1\)에 대한 추정량의 분산을 증가시킬 수 있습니다.

결론적으로 부적절한 변수가 모델에 포함될 경우, 모수 추정에 편향성(bias)는 나타나지 않지만 모수에 대한 추정치의 표준오차가 증가하는 결과가 나타납니다.

적절한(relevant) 변수의 제외

한편, 우리가 원래 추정해야할 모델이 \(y = \beta_0 + \beta_1x + \beta_2z + u\)라고 해보겠습니다. 그런데 모델을 잘못 특정해서 \(y = \beta_0 + \beta_1x + e\)로 추정했다고 할 때, \(e = \beta_2z + u\)라고 할 수 있습니다. 이 경우에는 무엇이 문제일까요? MLR.1부터 MLR.5까지의 가정들이 위배되었다고 할 수 있을까요?

일단 \(x\)와 \(z\)의 상관관계가 없다, 즉 \(Cov(x, z) = 0\)이라고 생각해보겠습니다. 이 경우에는 아무 문제가 없습니다. 왜냐하면 우리가 추정해야 하는 정확한 모델의 \(\beta_1\)와 잘못 특정한 모델의 \(\beta_1\)가 동일하기 때문입니다. 다만 잘못 특정한 모델의 오차항의 크기가 클 뿐입니다. 왜냐하면 오차항으로부터 \(y\)를 설명할 수 있는 \(z\)라는 요인을 추출해내지 못했기 때문입니다.

문제는 바로 \(Cor(x, z)\neq0\)일 경우입니다. Wooldridge(2016: 89-90)에서도 살펴볼 수 있듯이, 우리는 네 가지 누락변수로 인한 편의(Omitted variable bias, OVB)를 생각해볼 수 있습니다.

  • 만약 \(Cov(x, z) > 0, \beta_2 > 0\)이면, \(\hat{\beta_1} > \beta_1\)이므로 이 경우는 Positive OVB라고 합니다.

  • 만약 \(Cov(x, z) > 0, \beta_2 < 0\)이면, \(\hat{\beta_1} < \beta_1\)이므로 이 경우는 Negative OVB라고 합니다.

  • 만약 \(Cov(x, z) < 0, \beta_2 > 0\)이면, \(\hat{\beta_1} > \beta_1\)이므로 이 경우는 Negative OVB라고 합니다.

  • 만약 \(Cov(x, z) < 0, \beta_2 < 0\)이면, \(\hat{\beta_1} < \beta_1\)이므로 이 경우는 Positive OVB라고 합니다.

만약 적절한 변수를 제외하고 만든 모델로 OLS 추정을 하게될 경우에 우리는 \[\begin{equation*} \tilde{\beta_1}=\hat{\beta_1} + \frac{\sum^N_{i=1}(x_{i1}-\bar{x}_1)e}{\sum^N_{i=1}(x_{i1}-\bar{x}_1)^2} \end{equation*}\] 를 얻게 됩니다. 즉, 원래 추정하고자 했던 \(\hat{\beta_1}\)에 비해 뭔가가 더 붙는 거슬 확인할 수 있습니다. 그렇다면 이 우변의 기대값을 구해보도록 하겠습니다.

  • 일단, 위에서도 살펴보았던 것처럼 \(e = \beta_2z + u\)입니다.

  • 이걸 방금 전의 OLS 추정치, 우변에다가 대압해보겠습니다. 정말 끔찍한 수식이 나오지만 원래 추정하고자 했던 \(\hat{\beta_1}\)에 뭔가 점점 잡다한게 더해지고 있다는 것에서부터 이게 문제가 있다는 게 짐작이 가시겠죠?

\[\begin{equation*} \tilde{\beta_1}=\hat{\beta_1} + \frac{\hat{\beta_2}\sum^N_{i=1}(x_{i1}-\bar{x}_1)x_{i2} + \sum^N_{i=1}(x_{i1}-\bar{x}_1)e}{\sum^N_{i=1}(x_{i1}-\bar{x}_1)^2} \end{equation*}\]

이 수식으로 주어진 설명변수들에 대한 기대값을 구해보면,

\[\begin{equation*} \begin{aligned} E(\tilde{\beta_1}|x_k)&= \hat{\beta}_1\\ &\:\:\:+\frac{\hat{\beta_2}\sum^N_{i=1}(x_{i1}-\bar{x}_1)x_{i2} + \sum^N_{i=1}(x_{i1}-\bar{x}_1)E(e|x_k)}{\sum^N_{i=1}(x_{i1}-\bar{x}_1)^2}\\ &= \hat{\beta}_1 + \hat{\beta}_2\frac{\sum^N_{i=1}(x_{i1}-\bar{x}_1)x_{i2}}{\sum^N_{i=1}(x_{i1}-\bar{x}_1)^2}\\ &= \hat{\beta}_1 + \hat{\beta}_2\widehat{Cov(x_1, x_2)}/\widehat{Var(x_1)} \end{aligned} \end{equation*}\] 라는 결과를 도출하게 됩니다. 따라서, 적절한 변수를 모델에서 제외할 때 생기는 편향, 누락변수의 편향(OVB)의 크기는 \[\begin{equation*} \text{Bias}(\tilde{\beta_1}) = E(\tilde{\beta_1})|x_k - \beta_1 = \beta_2\frac{\widehat{Cov(x_1, x_2)}}{\widehat{Var(x_1)}} \end{equation*}\] 가 됩니다. 그렇다면, 위에서 살펴본 것 같이 두 경우를 생각해볼 수 있겠죠?

  • \(\hat{\beta_2}=0\)

  • \(\widehat{Cov(x_1, x_2)}=0\)

이 경우들에서는 편향성이 0이 됩니다. 따라서 일반적으로 종속변수에 영향을 미치는 변수를 누락시키는 문제는 누락변수가 포함된 변수들과 독립적이지 않은 한에야 포함된 변수들의 OLS 추정값들의 편향시키는 결과를 가져옵니다. 그리고 그 편향의 방향성과 크기는 위에서 살펴본 네 가지 경우의 수에 따라서 나타날 수 있으며, \(\hat{\beta_2}\)와 \(\widehat{Cov(x_1, x_2)}\)의 부호와 크기에 따라 결정됩니다.

Footnotes

  1. 저는 종속변수를 예측하기 위한 예측변수(predictor), 혹은 설명하기 위한 설명변수(explanatory variable)이라고 하고 독립변수(independent variable)라는 표현은 지양합니다. 왜냐하면 이 변수들이 외생적(exogeneous)이고 독립적이라는 근거는 어디에도 없기 때문입니다.↩︎

  2. \(\bar{x}\)와 \(\bar{y}\)는 주어진 표본에서의 특정한 “수”, 평균입니다.↩︎

  3. 구분하기 편하게 연속형 변수는 양적 변수(quatitative variables)로, 명목형과 순위형 변수는 질적 변수(qualitative variables)로 사용하겠습니다.↩︎

  4. 여기서 염두에 두어야할 것은 본래 QOG.s는 1인당 GDP, 무역 개방성, 그리고 노령화 지수를 대상으로 서브셋으로 만든 후에 결측치를 제거했기 때문에 노령화 지수가 결측치일 경우, 그에 해당하는 id의 1인당 GDP와 무역 개방성도 제외되었을 수 있습니다. 편의상 이전에 사용한 서브셋을 다시 사용하는 것이므로 이 점을 감안해야 합니다.↩︎

  5. 어떤 변수의 영향력 하에서 종속변수의 평균이 그 변수의 영향력 외에 있을 때의 평균과 차이가 있는지 여부, 이것은 영가설 기각 논리의 배경이기도 합니다.↩︎