오차(Errors)

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

회귀분석의 가정: \(\Pr(Y|x)\)의 조건 분포는 동일한 분산을 가진 정규분포를 따르며, \(Y\)의 조건 평균(여기서는 \(\mu_1, \mu_2, \dots, \mu_5\))은 모두 일직선 상에 놓이게 됩니다 (Fox 2016, 17).

• 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\)가 됩니다.

계수(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)에 대해 짚고 넘어갈 필요가 있습니다.

Best Linear Unbiased Estimator

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

모델의 특정(specification)

다양한 변수들 (Various variables)

## 필요한 패키지를 불러오겠습니다.
library(ezpickr)
library(ggplot2)
library(futurevisions)
library(patchwork)
library(tidyverse)

## 먼저 가상의 데이터를 만들어 보겠습니다.

fake.data1 <- tibble(Consumption = seq(45, 94.5, 0.5) + 
                       rnorm(100, 5, 4),
                    Gender = rep("male", 100),
                    region = gl(n = 2, k = 50,
                               length = 2*50,
                               labels = c("PK", "TK"),
                               ordered = F),
                    Income = seq(50, 149, 1) + 
                      rnorm(100, 2, 10))

fake.data2 <- tibble(Consumption = seq(10, 59.6, 0.5) +
                       rnorm(100, 2, 8),
                    Gender = rep("female", 100),
                    region = gl(n = 2, k = 50,
                               length = 2*50,
                               labels = c("Seoul", "Honam"),
                               ordered = F),
                    Income = seq(30, 129, 1) + 
                      rnorm(100, 2, 10))

fake.data <- bind_rows(fake.data1, fake.data2) %>% 
  mutate(
  Gender = Gender %>% parse_factor(., levels = c("female", "male"),
                                   include_na = F, ordered = T),
  region = region %>% as.character() %>%
    parse_factor(., levels = c("PK", "TK", "Honam", "Seoul"),
                 include_na = F, ordered = T)
)
fake.data

## 각 성별에 따라서 기울기는 같지만 절편은 다른 회귀모델을 만들겠습니다.

model1 <- lm(Consumption ~ Gender + Income, data = fake.data)
summary(model1)

## 
## Call:
## lm(formula = Consumption ~ Gender + Income, data = fake.data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -20.0571  -5.8740   0.4099   4.4216  22.5627 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 16.75799    1.84282   9.094   <2e-16 ***
## Gender.L    21.83555    0.85002  25.688   <2e-16 ***
## Income       0.42513    0.01916  22.186   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.123 on 197 degrees of freedom
## Multiple R-squared:  0.8921, Adjusted R-squared:  0.891 
## F-statistic: 814.7 on 2 and 197 DF,  p-value: < 2.2e-16

## 회귀분석의 예측값을 별도의 열로 추가한 데이터를 만듭니다.
predict.fake <- bind_cols(fake.data, pred.gender = predict(model1))

predict.fake %>% 
  ggplot(aes(x = Income, y = Consumption, color=Gender)) + 
  geom_point(alpha = 1/3) + 
   geom_smooth(aes(y = pred.gender, color = Gender),
              method="lm", se=FALSE, fullrange=TRUE, 
              alpha = 1/3) + 
  scale_color_manual(values=futurevisions("mars")) + 
  theme_bw() + theme(legend.position = "bottom")

model2 <- lm(Consumption ~ region + Income, data = fake.data)
summary(model2)

## 
## Call:
## lm(formula = Consumption ~ region + Income, data = fake.data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -21.6209  -5.2047  -0.2184   4.2042  21.8285 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  27.59874    3.22962   8.545 3.67e-15 ***
## region.L    -29.24316    1.25653 -23.273  < 2e-16 ***
## region.Q     -8.40970    2.08702  -4.030 8.01e-05 ***
## region.C     15.46734    1.13612  13.614  < 2e-16 ***
## Income        0.30649    0.03482   8.802 7.20e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.832 on 195 degrees of freedom
## Multiple R-squared:  0.9007, Adjusted R-squared:  0.8987 
## F-statistic: 442.3 on 4 and 195 DF,  p-value: < 2.2e-16

predict.fake <- bind_cols(fake.data, pred.region = predict(model2))

predict.fake %>% 
  ggplot(aes(x = Income, y = Consumption, color=region)) + 
  geom_point(alpha = 1/3) + 
  geom_smooth(aes(y = pred.region, color = region),
              method="lm", se=FALSE, fullrange=TRUE,
              alpha = 1/3) + 
  scale_color_manual(values=futurevisions("titan")) + 
  theme_bw() + theme(legend.position = "bottom")

간단한 정리

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

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

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()

영가설 유의성 검정 (Null Hypothesis Significance Testing; NHST)

우리가 생각하는 모집단에 대한 \(\beta_k\)가 사실(true)라고 가정해보겠습니다. \(\beta_k\)가 보여주고자 하는 모수를 \(w\)라고 하겠습니다. \(\beta_k = w\)라는 가정 하에서, 우리는 관측가능한 표본들로부터 추출한 \(\hat{\beta}_k\)들로부터 나타나는 분포가 얼마나 그 \(w\)를 보여주는 \(\beta_k\)를 포함할 가능성을 가지는지를 생각해볼 수 있습니다. 나아가 \(\beta_k = w\)라고 가정한다면 우리는 사실상 알 수 없는 모수를 아는 것과 다름없게 됩니다. 왜냐하면 \(\beta_k\)의 표집분포의 기대값이 \(w\)를 제대로 보여준다면, 이후에는 충분한 수의 표집과정을 통해 \(\beta_k\)를 충분히 확보만 하면 되는 것이니까요.

기술적으로 우리는 이것을 \((n-k-1)\)의 자유도를 가진 Student t 분포라고 합니다. 그리고 상당한 수의 관측치를 확보하게 된다면, 이 분포는 가우시안(Gaussian) 분포라고 할 수 있습니다. 아래와 같은 공식으로 나타낼 수 있겠네요.

\[ \frac{\hat{\beta_k} - \beta_k}{se(\hat{\beta})} \sim t_{n-k-1} \]

그렇다면 이제 영가설 유의성검정의 절차를 한 번 살펴보도록 하겠습니다. 영가설을 \(H_0: \beta_k = w\)이라고 설정해보도록 하겠습니다. 그리고 이때 t 통계치는 \(\frac{\hat{\beta_k} - w}{se(\hat{\beta_k}-w)} \sim t_{n-k-1}\)로 나타낼 수 있겠네요. 그렇다면 여기서 문제입니다. \(\beta_k\)에 대한 영가설 혹은 연구가설을 기각하기 위해서 우리에게 필요한 t 통계치는 얼마일까요? 즉 \(\hat{\beta_k}\)은 얼마나 \(\beta_k\)로부터 떨어져 있어야할까요?

먼저 영가설 \(H_0: \beta_k = w\)에 대해 그것과 비교하기 위한 대안가설(연구가설)을 \(H_1: \beta_k > w\)이라고 해보겠습니다. 즉, 우리가 관심있는 것은 \(\beta_k > w\)이지만 이것을 직접적으로 검증하거나 혹은 확증할 수는 없기 때문에 \(\beta_k \neq w\)라는 것을 통해 영가설을 기각함으로써 영가설이 기각될 확률, 연구가설이 유의미할 ’확률’을 구하게 되는 것입니다.

이것이 우리가 흔히 말하는 유의수준(significance level)입니다. 기술적으로 서술하자면 “모집단에서 추출한 100개의 표본 중에서 영가설이 사실일 경우를 기각하는 것을 몇 번 관측할 수 있느냐”라고 하는 것입니다. 만약 100개의 표본 중에서 12개의 표본이 영가설의 기대대로 \(\beta_k = w\)라는 것을 보여줬다면, 영가설이 사실일 확률은 0.12 이 될 것입니다. 이때 우리는 이 영가설이 사실일 확률 0.12를 \(\alpha\)라고 표현합니다. \(\alpha\)의 값이 더 작을 수록 영가설을 기각하는 \(\hat{\beta_k}\)가 더 많다는 것을 의미합니다.

위의 플롯은 영가설(\(H_0\))이 \(\beta_k = 0\)이라고 했을 때, 연구가설(\(H_1\))이 \(\beta_k > 0\)인 경우, 관측치가 10,000개인 표본에서의 기각역(rejection region)을 보여주고 있습니다. 일단 전체 곡선 면적의 합은 1입니다. 당연하겠죠? 밀도함수로 나타낸 분포이니까 전체의 총합은 1입니다. 그리고 플롯에서 기각역이라고 나타난 선 우측의 면적의 총합은 0.05가 됩니다. 통상적으로 그것은 우리가 상정한 영가설대로 \(\beta_k = 0\)일 확률을 의미합니다.

이번에는 영가설(\(H_0\))이 \(\beta_k = 0\)이라고 했을 때, 연구가설(\(H_1\))이 \(\beta_k < 0\)인 경우, 관측치가 10,000개인 표본에서의 기각역(rejection region)을 살펴보겠습니다. 마찬가지로 전체 곡선 면적의 합은 1입니다. 앞서의 플롯과 다른 점은 기각역의 위치입니다. 연구가설이 상정한 계수의 부호가 달라졌기 때문에 기각역의 위치도 달라진 것이죠. 따라서 이때는 선 좌측의 면적의 총합이 0.05가 되고, 영가설대로 \(\beta_k = 0\)일 확률을 의미합니다.

우리는 t 통계치가 우리가 상정한 100개의 표본 중에서 n개의 표본만이 영가설에 부합할 것이다라고 할 수 있는 일종의 결정적 기준값(critical value)보다 클 경우에 영가설을 기각할 수 있습니다. 이 결정적 기준값은 유의수준(\(\alpha\)) 혹은 전체 확률에서 유의수준을 제한 값(\(1-\alpha\))을 Student t 혹은 정규누적밀도함수에 대입하면 구할 수 있습니다.

영가설의 기각 이후

영가설을 기각했다고 해보겠습니다. 그렇다면 우리는 그 기각 결과를 가지고 어떻게 연구가설에 대해 설명할 수 있을까요?

사실 \(\alpha\)라고 하는 특정한 값을 사전에 미리 설정하고 그것을 선긋듯이 어떠한 결과를 결정하는 도구로 사용한다는 것은 문제의 소지가 있어 보입니다. 예를 들어 \(p\) 값이 0.051인 경우와 0,049일 때, 우리는 이들의 통계적 의미를 어떻게 이해해야할까요?

\(p\) 값은 주어진 t 통계치 하에서 영가설을 기각할 수 있는 가장 작은 \(\alpha\)를 의미합니다. 공식을 통해서 살펴보자면, \(T\)가 우리가 얻을 수 있는 모든 가능한 검정 통계량(test statistics)라고 해보겠습니다. 예를 들어 Student t나 Gaussian이 있겠네요. 이때 만약 가설이 특정한 관계의 방향을 설정하지 않고 있다면, 단지 영가설 \(\beta_k = 0\)만을 기각하면 되기 때문에 이른바 양측검정을 상정한 연구가설을 설정하게 될 것입니다. 따라서 이때는 관계의 방향을 상정하지 않는 \(Pr(|T|>|t|) = 2Pr(T>|t|)\)가 성립하게 됩니다.

이제까지는 \(\beta_k = 0\)인 경우를 중심으로 살펴보았습니다. 왜냐하면 사실 회귀분석 같은 경우 영가설은 우리가 기대한 변수가 종속변수에 대해 ’효과가 없을 것’을 주로 기대하기 때문입니다. 따라서 \(\beta_k = 0\)이 기각되었다는 것은 우리가 관심을 가지고 있는 변수가 종속변수에 대해 유의미한 영향이 있을 수 있다는 것을 시사합니다.

하지만 사실 영가설은 반드시 0, 효과가 없다는 것만을 기대할 필요는 없습니다. 예를 들어서 \(H_0: \beta_k = g\)라는 특정한 값을 가지는 영가설에 대해서 \(H_1: \beta_k > g\)와 같이 계수값이 그 특정한 값보다 클 경우를 상정할 수 있습니다. 이때 유의수준 \(\alpha\)를 설정한다면, \(t = \frac{\hat{\beta_k} - g}{se(\hat{\beta_k})}\)가 되므로, 우리는 \(t > c\)이기만 하면 영가설을 기각할 수 있습니다.⁵ 조금 더 풀어서 말하자면 일반적인 t 통계량은 다음과 같이 쓸 수 있습니다.

\[ t = \frac{\text{추정치 - 가설로 기대하는 값}}{\text{추정치의 표준오차}} \]

이제 정치학 분야에서 실제로 사용할법한 데이터를 이용해서 위의 내용을 한 번 더 살펴보도록 하겠습니다. Quality of Government 데이터셋을 이용해서 \(y = a_0 + \sum_k \beta_k x_k + e\)의 형태를 갖는 선형회귀모델을 한 번 만들어 보겠습니다.

# 단순선형회귀모델과 다중선형회귀모델
model.simple <- 
  lm(wdi_gdpcapcon2010 ~ wdi_trade, data=QOG.sub)
result1 <- broom::tidy(model.simple) %>%
  mutate(model = "SLR")
model.multiple <- 
  lm(wdi_gdpcapcon2010 ~ wdi_trade + p_polity2 + 
       wdi_pop1564, data=QOG.sub)
result2 <- broom::tidy(model.multiple) %>% 
  mutate(model = "MLR")
results <- bind_rows(result1, result2) %>% 
   mutate_if(is.numeric, round, 3)

이렇게 만들어진 두 선형회귀모델을 가지고 영가설이 \(\beta_k = 0\)인 경우와 \(\beta_k = c\)인 경우를 한 번 살펴보도록 하겠습니다. 물론 이 경우에는 \(c\)가 합당한 비교 기준이라고 정당화할 수 있어야할 것입니다. 다중선형회귀모델을 기준으로 모델을 수식으로 표현한다면 아래와 같습니다.

\[ \begin{aligned} \text{Economy} = &\beta_1\text{Trade} + \beta_2\text{Level of Democracy}\\ &+ \beta_3\text{Working Population} + e. \end{aligned} \]

results %>% knitr::kable()

여기서 \(x_1\)에 대한 계수, \(\beta_1\)을 중심으로 일단 영가설을 특정해보도록 하겠습니다.

• 먼저 \(H_0: \beta_1 = 0\)이라고 할때, 민주주의 변수, 노동가능 인구비율이 통제되었을 때, GDP 대비 재화와 서비스의 수출 및 수입의 총합으로 측정된 무역개방성의 한 단위 증가는 2010년 고정 달러로 측정된 1인당 GDP, 즉 그 국가의 경제규모에 미치는 효과가 없다라는 것이 영가설임을 이해할 수 있습니다.

• 그 다음으로는 결정적 기준값으로 영가설을 수립해보도록 하겠습니다. 이번의 영가설은 \(\beta_1 = 164.355\)라고 하겠습니다. 즉, 무역 개방성의 한 단위 변화는 1인당 GDP로 측정된 경제규모가 163.48 증가하는 만큼의 효과를 가지고 있다는 것을 의미합니다.

• 여기서 결정적 기준값을 \(\beta_1 = 164.355\)라고 설정한 것은 단순선형회귀 모델에서의 무역 개방성의 계수값입니다. 따라서 이 영가설을 통해 우리는 단순선형회귀모델에서 무역 개방성과 경제 규모 간의 관계가 다중선형회귀 모델에서 다른 통제변수들이 통제된 가운데 나타나는 무역 개방성과 경제 규모 간의 관계와 다른지, 그리고 그러한 차이가 통계적으로 유의미한지를 살펴볼 수 있습니다.

그렇다면 이번에는 우리가 관심을 가지고 있는, 연구가설에 대해 살펴보겠습니다. 효과의 방향(부호)를 생각할 필요가 없는 경우(양측, two-sided)와 방향도 고려해야 하는 경우(단측, one-sided) 모두를 특정해보도록 하겠습니다.

• 효과가 없음이 영가설일 때(\(H_0: \beta_1 = 0\))

  1. 단측 연구가설(\(H_A: \beta_1 > 0\), 또는 \(H_A: \beta_1 < 0\)): 무역 개방성은 경제규모에 대해 긍정적 또는 부정적 효과를 가지고 있다. 이때의 연구가설은 긍정적 효과 또는 부정적 효과가 별개의 것이다.

  2. 양측 연구가설(\(H_A: \beta_1 \neq 0\)): 무역 개방성은 경제규모에 대해 ‘효과가 있다.’

• 연구가설이 \(H_A: \beta_1 = 164.355\)일 경우

  1. 단측 연구가설(\(H_A: \beta_1 > 164.355\), 또는 \(H_A: \beta_1 < 164.355\)): 경제규모에 영향을 미칠 수 있는 다른 요인들을 고려하지 않을 때보다 고려했을 경우(통제했을 경우), 무역 개방성이 경제규모에 미치는 효과가 더 크다/작다.

  2. 양측 연구가설(\(H_A: \beta_1 \neq 164.355\)): 경제규모에 영향을 미칠 수 있는 다른 요인들을 고려하지 않을 때와 고려했을 경우(통제했을 경우), 무역 개방성이 경제규모에 미치는 효과는 다르다(같지 않다).

이번에는 각각의 영가설과 연구가설에 대한 검정을 시행하고 그 결과를 제시 및 해석해보도록 하겠습니다.

# 바로 이전의 분석을 이용해보도록 하겠습니다.
# 단순선형회귀모델의 무역 개방성에 대한 계수와 표준오차, 자유도
b1.simple <- results[2,2]
se1.simple <- results[2,3]
simple.df <- model.simple$df

# 다중선형회귀모델의 무역 개방성에 대한 계수와 표준오차, 자유도
b1.multiple <- results[4,2]
se1.multiple <- results[4,3]
multiple.df <- model.multiple$df

# 첫 번째 연구가설 중 단측 가설입니다.
# HA: beta_1 > 0
pt(as.numeric((b1.simple-0)/se1.simple), 
   simple.df, lower = FALSE)

## [1] 4.131871e-08

# 0.05보다 훨씬 작은 값을 확인할 수 있습니다.
# 즉, 영가설의 기대대로 표본에서 결과를 얻을 확률이 매우 낮다는 것이고
# 영가설을 기각하기에 충분한 경험적 근거를 확보했다고 할 수 있습니다.

# 양측 가설을 한 번 살펴보겠습니다. 방향을 고려할 필요가 없죠.
# HA: beta_1 != 0
2*pt(-abs(as.numeric((b1.simple-0)/se1.simple)), simple.df)

## [1] 8.263742e-08

# 방향을 고려할 필요가 없기 때문에 방향을 고려했을 때의 
# 확률에 2를 곱해줍니다 (좌측 + 우측)

# 두 번째 연구가설(단순 vs. 다중) 중 단측 가설입니다.
# HA: beta_1_Multi > beta_1_Simple
pt(as.numeric((b1.multiple-b1.simple)/(se1.multiple-b1.simple)), 
   multiple.df, lower = FALSE)

## [1] 0.3761995

# 0.05보다 훨씬 큰 값, 영가설은 기각되지 않습니다.
# 당연합니다. 다중회귀분석에서의 무역 개방성의 계수값과
# 단순회귀분석의 무역 개방성의 계수값은 일단 단순회귀분석의
# 계수값이 더 컸습니다. 따라서 위의 연구가설, 다중회귀분석의
# 계수값이 단순회귀분석의 그것보다 클 것이라는 기대는 기각되는 것이
# 우리의 상식에 부합합니다.

# HA: beta_1_Multi != beta_1_Simple
2*pt(-abs(as.numeric((b1.multiple-0)/(se1.multiple-b1.simple))), multiple.df)

## [1] 0.3731572

# 방향을 고려할 필요가 없기 때문에 방향을 고려했을 때의 
# 확률에 2를 곱해줍니다 (좌측 + 우측). 마찬가지로 영가설 기각.
# 둘이 서로 같지 않다는 기대죠? 하지만 두 계수 값의 차이는 각 계수값이
# 서로 다른 표본에서 계산되어 가지게 될 표집분포에 따른 편차,
# 즉, 표준오차에 따라 고려했을 때, 계수의 차이가 유의미하기 보다는
# 표본의 차이에 따라 나타날 수 있는 차이라고 볼 수도 있습니다.
# 그 결과 영가설은 기각되지 않습니다. 
# 두 효과는 통계적으로 차이가 있다고 말하기에는 영가설을 
# 기각할 수 있는 충분한 경험적 근거를 확보하지 못한 것입니다.

\(\beta_k\)와 \(\hat{\beta_k}\)에 대한 또 다른 접근

서로 다른 데이터를 보여주는 간단한 분포 세 개를 보도록 하겠습니다.

세 분포의 공통점과 차이점은 무엇일까요? 우선 세 분포는 평균(mean)과 표준편차(standard deviation)가 매우 비슷합니다. 하지만 실제 관측치들의 분포 양상은 매우 다르죠. 첫 번째 분포는 정규분포를 따른다면, 두 번째 분포는 지수분포, 마지막 분포는 크게 양극화된 양봉분포(bimodal distribution)의 형태를 띄고 있습니다.

단순히 데이터를 평균과 표준편차만으로 섣불리 이해할 경우 그 데이터를 이용한 통계적 추론에서 오류를 범할 수 있습니다. 따라서 우리는 ‘보다 구조화된’ (more structured) 접근법을 취할 필요가 있습니다.

우선 가장 클래식한 선형회귀모델의 가정들을 다시 한 번 생각해보겠습니다. 총 다섯 개의 세부 가정으로 살펴볼 수 있는데, 순서대로 다중선형회귀분석(Multiple linear regression)에 대한 가정이라는 의미로 MLR라고 라벨링을 해보도록 하겠습니다. 어디까지나 여기서 클래식한 다중선형회귀분석의 가정들은 교차사례(cross-sectional) 자료에 적용가능합니다. 시계열(time-series)이 포함되면 또 고려해야할 문제들이 있기 때문인데, 이는 아마 고급 통계파트의 자료에서 다루게 될 것 같습니다.

• MLR1. 모집단에서의 모수들의 관계가 선형이다.

• MLR2. 무작위로 추출된 표본이다.

• MLR3. 설명변수들 간의 완벽한 다중공선성은 존재하지 않는다.

• MLR4. 설명변수와 오차는 독립적이다(Zero-conditional mean)

  • 주어진 설명변수라는 조건 하에서 오차항(\(u\))의 기대값은 0이다.

  • 만약 이 가정이 성립되지 않는다면 우리는 모델 특정에 문제가 있었다고 생각할 수 있다. 예를 들면, 종속변수에 영향을 미치는 적절한 변수를 모델에 포함시키지 않아 오차항이 설명변수와 상관성을 가진다고 볼 수 있다.

  • 이 MLR4는 선형회귀분석의 추정치, 계수의 편향성 문제와 매우 밀접한 관계를 가진다.

• MLR5. (오차항의) 동분산성

  • 주어진 설명변수들에 대해 오차항의 조건 분포가 일정해야 한다.

  • 만약 오차항이 이분산성(heteroskedasticity)을 띈다고 하더라도 그것이 OLS 추정치가 편향되었다는 것을 의미하지는 않는다.

  • 다만 OLS 추정치의 ’효율성’이 낮다는 것을 의미한다. 즉, OLS 중 least 의 조건 최소 분산(least variance)이 아니라는 것을 의미할 뿐이다.

• MLR1부터 MLR5까지가 충족되었을 때, 우리는 OLS의 결과를 BLUE (Best Linear Unbiased Estimates)라고 한다.

자, 이 다섯 가정을 우리는 한 데 묶어서 Gauss-Markov 가정이라고 합니다. 이 클래식한 선형회귀모델의 가정을 조금 단순하게 묶어보자면 다음과 같이 표현할 수 있을 것 같습니다.

• \(E(\hat{\beta}_k) = \beta_k\). 당연하겠죠? 만약 BLUE라면 우리가 표본들로부터 얻는 평균들의 평균이 모집단의 기대값과 일치할 거라는 가정 또한 충족될 것입니다.

• \(Var(\hat{\beta}_k)\)는 표본들로부터 얻은 평균들의 분산으로 어떠한 값을 가지게 될 것입니다. 뭐, 놀랄 건 아니죠. 모집단에서 표본을 추출하는 것은 본연적으로 일정한 오류를 내재하게 되어있기 때문에 평균들도 항상 모집단의 기대값과 같지는 않을 것이고, 그런 표본들의 평균들 간의 차이는 일종의 분산으로 나타날 것이니까요.

그런데 여기서 우리는 위의 두 발견으로부터 하나의 식을 도출해낼 수 있습니다. 만약 우리가 구한 OLS 추정치가 BLUE라면? 그래서 위의 두 조건이 충족된다면? 우리는 표본들로부터 구한 평균들이 어떠한 분포를 보일 것이라고 생각해볼 수 있습니다. 바로 모집단의 기대값을 보여줄 \(E(\hat{\beta}_k)\)를 중심으로 \(Var(\hat{\beta}_k)\)라는 분산을 가진 정규분포입니다.

\[ \hat{\beta}_k \sim N(\hat{\beta}_k, Var(\hat{\beta}_k)) \]

앞서 살펴보았던 \(\hat{\beta}_k \sim N(\hat{\beta}_k, Var(\hat{\beta}_k))\)가 성립한다고 할 때, 우리가 “알 수 있는 것”과 “알 수 없는 것”은 무엇일까요?

• 먼저 우리가 알 수 있는 것은 크게 세 가지라고 할 수 있습니다. 우리가 가지고 있는 예측변수이자 확률변수인 \(x\)가 어떠한 평균과 분산을 가진 정규분포를 따를 것이라는 것을 알 수 있습니다: \(x\text{에 대한 분포} \sim N(\mu, \tau^2): \frac{1}{2\sqrt{\tau^2\pi}}\text{exp}[-\frac{(x-\mu)^2}{2\tau^2}]\). 그리고 실제로 표본을 통해서 \(\hat{\beta_k}\)와 \(Var{\hat{\beta_k}}\)도 구할 수 있습니다.

• 반면에 모집단의 모수, \(\beta_k\)는 알 수 없습니다.

만약 우리가 \(\beta_k\)에 대한 가정을 세울 수 있고, 이론적 분석이 아니라 실제 경험적 데이터를 통해 그러한 가정이 충족되는지 아닌지를 확인할 수 있다면 어떻게 될까요?

신뢰구간(confidence interval)

영가설 유의성 검정, 즉 NHST에서 우리는 이미 \(\alpha\)에 대해서 살펴보았습니다. 알 수 없는 모수 \(\beta_k\)에 대한 알고 있는 통계치 \(\hat{\beta_k}\)를 바탕으로 한 분포의 가정에서 우리는 다시 이 \(\alpha\)로 돌아옵니다. \(\alpha\)라는 개념의 바탕에는 “이론적으로 무수히 많은 표본”이 전제되어 있습니다. 적어도 이론적 수준에서 우리는 하나의 모집단으로부터 수많은 표본을 반복하여 추출할 수 있습니다. 그렇게 무수히 많은 표본들로부터 얻은 통계치의 분포를 표집분포(sampling distribution)이라고 할 수 있고, 그 표집분포를 바탕으로 우리는 얼마나 많은 표본이 영가설의 기대에 부합하는 통계치를 가지는지, 전체 표본에서 그러한 표본의 비율을 \(\alpha\)라고 나타냅니다. 따라서 흔히 사용하는 유의수준 0.05라는 기준은 이런 맥락에서 보자면 “100번의 반복 추출된 표본 중에서 5개의 표본의 표본평균이 영가설의 기대를 포함하고 있을 경우”를 말한다고 할 수 있습니다. 그러면 \(\alpha\)는 조금 더 직관적으로 이해할 수 있습니다.

• \(\alpha\)는 표집분포에 있어서 일종의 꼬리확률(tail probability)이라고 생각할 수 있습니다.

• 신뢰구간은 우리가 일반적으로 모수가 그 사이에 존재할 것이라고 기대하는 구간입니다.

• 반복해서 추출한 표본들에 대해 95%의 신뢰구간은 100번 추출한 표본의 \(\hat{\beta_k}\) 중에서 95개의 표본이 모수 \(\beta_k\)를 포함한다고 기대할 수 있습니다.

행렬? 벡터? 전치행렬? 역행렬?

저도 OLS 기본을 마무리하고 좀 더 깊게 들어갔을 때, 갑자기 [오랜만에] 등장한 이 행렬 때문에 당황했습니다. 하지만 R로 행렬들이 계산되면서 변화하는 것을 추적하며 공부하면 예전에 수학의 정석 붙잡고 노트에 필기하며 낑낑대었던 것보다는 훨씬 편하게 이해하실 수 있을 겁니다.

대체 행렬(matrix), 벡터(vector), 전치행렬(transposed matrix), 역행렬(inversed matrix)라는 이놈들이 왜 OLS를 공부하는 데 등장하는 것일까요? 그리고 OLS에서 \(X\beta\)가 제대로 돌아가기 위해서 가정되어야 하는 것들은 무엇이 있을까요? 마지막으로 이것들을 R로 어떻게 보여줄 수 있을까요? 참고로 \(X\beta\)는 OLS를 통해 얻은 일련의 계수값을 \(\beta\)라는 하나의 집단으로 나타내고 \(X\)도 마찬가지로 모델에 포함된 설명변수 일체를 표시한 것입니다.

행렬과 벡터를 정의하기에 앞서, 우리는 스칼라(scalar)라는 것에 대해서 알아야만 합니다. 저도 정치학을 전공한 사람이고 방법론은 연구문제를 풀기 위한 도구로 공부하는 사람이기 때문에 엄청 깊은 증명이나 디테일에 집착하지는 않겠습니다. 저도 지금 큰 문제없이 연구방법을 통해 연구문제들을 풀어가고 있으니, 제가 이해한 정도와 비슷한 정도로만 이해하셔도 연구에는 큰 지장은 없으실겁니다. 물론 주정공이나 부전공이 방법론이라던가 하는 경우에는 조금 얘기가 달라지겠죠? 어디까지나 부외자에 한정된 이야기입니다.

스칼라는 하나의 차원에서 측정된 수치화된 결과를 말합니다. 따라서 스칼라는 방향에 대한 정보는 가지지 않고 크기(magnitude)에 대한 정보만 가지고 있다고 할 수 있습니다. 벡터는 이런 스칼라들의 집합입니다. 스칼라와는 다르게 벡터는 크기뿐아니라 방향에 대한 정보도 가지고 있습니다. 예를 들면, \(n\)개의 스칼라를 가진 \(a\)라는 벡터가 있다고 해보겠습니다. 우리는 이 벡터를 \(n\)차원의 벡터라고 부를 수 있습니다. 이런 벡터는 아래의 수식 과 같이 열 또는 행으로 나타낼 수 있습니다.

\[\begin{equation} \label{eq:1} a = \begin{bmatrix}a_1 \\ a_2 \\ a_3 \\ \vdots \\ a_n\end{bmatrix} \text{ or } a = \begin{bmatrix}a_1, a_2, a_3, \cdots, a_n\end{bmatrix} \end{equation}\]

그렇다면 행렬은 무엇일까요? 행렬은 행과 열로 요소들(숫자일 수도 있고 일종의 함수들일 수도 있습니다)이 사각형의 형태를 띄게끔 배열된 결과물입니다. 행렬이라는 이름이 행(raw)과 열(column)이 합쳐진 것이니까 이 정도는 쉽게 유추할 수 있습니다. 행렬을 구성하는 행과 열의 수를 가지고 우리는 행렬의 차원(dimension)이 어떻게 되는지를 말할 수 있습니다. 예를 들면, 아래의 수식 는 요소들의 배열, 즉 행렬의 모습을 보여주고 있습니다. \(j\)는 행을, \(k\)는 열을 의미하는데요, \(a_{jk}\)는 \(j\)-행과 \(k\)-열에 위치한 요소를 나타냅니다.

\[\begin{equation} \label{eq:2} A = \begin{bmatrix}a_{jk}\end{bmatrix} = \begin{bmatrix} a_{11} & a_{12} & \cdots & a_{1n} \\ a_{21} & a_{21} & \cdots & a_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{m1} & a_{m2} & \cdots & a_{mn} \\ \end{bmatrix} \end{equation}\]

행렬의 정의에 따르면, 우리는 벡터 역시도 일종의 행렬이라고 이해할 수 있습니다. 단 하나의 행(1\(\times n\)) 또는 열(\(m \times1\))로 구성된 행렬인 거죠. 그리고 행렬을 전치(transposition)한다는 것은 말 그대로 위치(position)를 뒤바꾸는 것(trans)이기 때문에 행렬의 행과 열을 서로 뒤바꾼다는 것을 의미합니다(수식 ). 예를 들면 2개의 행과 3개의 열을 가지고 있던 \(2\times3\) 행렬이 3개의 행과 2개의 열을 가진 \(3\times 2\)의 행렬로 전치되는 것입니다.

\[\begin{equation} \label{eq:3} A = \begin{bmatrix}a_{jk}\end{bmatrix} \Rightarrow A^T = \begin{bmatrix}a_{kj}\end{bmatrix} = \begin{bmatrix} a_{11} & a_{21} & \cdots & a_{m1} \\ a_{12} & a_{22} & \cdots & a_{m2} \\ \vdots & \vdots & \ddots & \vdots \\ a_{1n} & a_{2n} & \cdots & a_{mn} \\ \end{bmatrix} \end{equation}\]

주어진 행렬의 곱셈에 대한 역원을 역행렬이라고 합니다.모든 행렬이 반드시 그 역행렬을 가지는 것은 아닙니다. A라는 행렬이 역행렬을 가진다고 한다면, 우리는 A라는 행렬이 “역행렬을 가질 수 있는 행렬”이라는 의미에서 가역행렬(invertible matrix)라고 부릅니다. 반면에 역행렬을 가지지 않는 행렬을 특이행렬(singular matrix)라고 합니다. 만약에 어떤 행렬이 역행렬을 가진다면, 이때 그 역행렬은 가역행렬에 대해 유일한 역행렬입니다. 1:1 관계가 성립한다는 거죠. A라는 행렬에 대한 역행렬을 우리는 A\(^{-1}\)라고 표기합니다.

우리는 행렬을 이용해서 단선순형회귀모델을 표현할 수 있습니다. 이때 필요한 것은 \(x\)와 \(y\)로, 실제 관측치 \(x_i\)를 요소로 하는 예측변수의 행렬이 모델의 계수값이라고 할 수 있는 행렬을 거치게 되면 \(\hat{y}_i\)의 행렬을 산출하는 논리입니다. 마찬가지로 다중선형회귀모델도 행렬로 표현할 수 있습니다. 그러나 이때에는 \(y_i\)의 변화를 설명하기 위해 두 개 이상의 예측변수들을 나타내는 행렬이 필요합니다. 즉, 최소 세 차원 이상의 행렬이 필요하다는 것입니다.

앞서 정리한 것처럼 벡터는 스칼라들의 집합입니다. 따라서 우리는 \(k\)개의 예측변수들과 \(n\)개의 관측치를 이용한 행렬들로 수식 와 같이 다중선형회귀모델을 나타낼 수 있습니다.

\[\begin{equation} \label{eq:4} y = \begin{bmatrix}y_1 \\ y_2 \\ \vdots \\ y_n \end{bmatrix} X = \begin{bmatrix}1 & x_{21} & x_{31} & \cdots & x_{k1}\\ 1 & x_{22} & x_{32} & \cdots & x_{k2}\\ \vdots & \vdots & \vdots & \ddots & \vdots\\ 1 & x_{2n} & x_{3n} & \cdots & x_{kn}\\\end{bmatrix} \beta = \begin{bmatrix} \beta_1 \\ \beta_2 \\ \vdots \\ \beta_n \\\end{bmatrix} u = \begin{bmatrix} u_1 \\ u_2 \\ \vdots \\ u_n \\\end{bmatrix} \end{equation}\]

수식 에 입각해 우리는 선형모델을 각각의 관측된 \(y\)값과 \(X\), 그리고 추정된 \(\beta\)와 \(u\)를 가지로 수식 로 보여줄 수 있습니다.

\[\begin{equation} \label{eq:5} \begin{split} &y_1 = \beta_1 + \beta_2x_{21} +\cdots+\beta_kx_{k1} + u_1\\ &y_2 = \beta_1 + \beta_2x_{22} +\cdots+\beta_kx_{k2} + u_2\\ &\vdots\\ &y_n = \beta_1 + \beta_2x_{2n} +\cdots+\beta_kx_{kn} + u_n\\ \end{split} \end{equation}\]

여기서 \(X\)와 \(y\)는 실제 관측된 데이터로부터 나온 예측변수와 종속변수의 개별 관측값을 의미합니다. 따라서 현실의 데이터는 항상 우리가 알지 못하는 요인들에 의해 영향을 받을 수 있기 때문에 우리는 이 ‘관측하지 못한’ 그리고 ‘비체계적인’ 요인들을 나타내는 오차항을 모델에 포함시킵니다. 그러므로 다중선형회귀모델을 보여주는 행렬은 아래의 수식 와 같이 정리할 수 있습니다.

\[\begin{equation} \label{eq:6} y = \begin{bmatrix}y_1 \\ y_2 \\ \vdots \\ y_n \end{bmatrix} = \begin{bmatrix}1 & x_{21} & x_{31} & \cdots & x_{k1}\\ 1 & x_{22} & x_{32} & \cdots & x_{k2}\\ \vdots & \vdots & \vdots & \ddots & \vdots\\ 1 & x_{2n} & x_{3n} & \cdots & x_{kn}\\\end{bmatrix} \begin{bmatrix} \beta_1 \\ \beta_2 \\ \vdots \\ \beta_n \\\end{bmatrix} + \begin{bmatrix} u_1 \\ u_2 \\ \vdots \\ u_n \\\end{bmatrix} \end{equation}\]

그리고 예측변수 \(X\)의 1은 절편값을 나타냅니다. 행렬 계산 공식을 떠올리지 못하더라도 이 수식 을 직관적으로 이해해보자면, 관측된 종속변수 \(y_i\)는 관측된 \(k\)개의 예측변수의 관측치들인 \(x_{ki}\)가 모델에 의해 추정된 \(\beta\)들과 결합하여 계산되는 예측값 \(\hat{y}_i\)에 현실에서 관측하지 못한 오차들이 더해진 결과라는 것입니다. 우리는 위의 행렬을 계산해서 OLS, 즉 오차의 최소제곱합(least sum of squared errors)을 구할 수 있습니다. 이건 그냥 참고로 알아두세요.

\[ \begin{split} S&=\Sigma^n_{i=1}\epsilon^{2}_{i} = \Sigma^n_{i=1}(Y_i - \beta_1 - \beta_2X_i1 - \cdots - \beta_kX_{ik})^2\\ &=\epsilon'\epsilon = (Y - X\beta)'(Y-X\beta) = Y'Y - Y'X\beta - \beta'X'Y + \beta'X'X\beta \end{split} \]

이렇게 계산된 \(S\)에 편미분을 취하면?

\[ \begin{split} \frac{\partial S}{\partial \beta} =&-2X'Y + 2X'Y\beta = 0\\ \rightarrow&X'X\hat{\beta}= X'Y\\ \rightarrow&\hat{\beta} = (X'X)^{-1}X'Y \end{split} \]

이렇게 됩니다. 이 행렬계산을 통해서 \(\hat{\beta}_k\)를 얻을 수 있고(이게 계수값이죠!), R은 이 행렬계산을 구현할 수 있습니다.

우리가 수많은 변수들과 관측치들을 가지고 있다고 해봅시다. \(i \leq n\), 즉 \(n\)개의 관측치를 가지고 있다고 할 때, 우리는 \[ y_i = \alpha + \beta_1x_{i1} + \cdots + \beta_{K}x_{iK} + \epsilon_i \]

라는 다중선형회귀모델을 구축할 수 있게 됩니다. 이는 예측변수들의 수에 절편까지를 고려한 \(K+1\)개의 제한된 변수를 가진 함수를 가지 모델이라는 것을 의미하며, 우리는 이 모델에서 오차의 제곱합이 최소가 되는 계수값들을 구할 수가 있게 되는 것입니다. 조금 더 직관적이게 행렬을 보여드리자면,

\[ y = \begin{bmatrix}y_1 \\ y_2 \\ \vdots \\ y_n \end{bmatrix}, x_k = \begin{bmatrix}x_{1k} \\ x_{2k} \\ \vdots & x_{Nk}\end{bmatrix}, \text{ 그리고 }, \epsilon = \begin{bmatrix} \epsilon_1 \\ \epsilon_2 \\ \vdots \\ \epsilon_n \\\end{bmatrix} \]

이라고 할 때, 각각의 \(y\), \(x_k\), \(\epsilon\)은 \(N\times 1\)의 벡터입니다. 관측된 종속변수, 예측변수, 그리고 오차항이죠. 이를 선형회귀모델에 집어넣으면 아래와 같이 다시 써볼 수 있습니다.

\[ y = \begin{bmatrix}y_1 \\ y_2 \\ \vdots \\ y_n \end{bmatrix} = \alpha + \beta_1 \begin{bmatrix}x_{11} \\ x_{21} \\ \vdots & x_{N1}\end{bmatrix} + \cdots + \beta_k\begin{bmatrix}x_{1K} \\ x_{2K} \\ \vdots & x_{NK}\end{bmatrix} + v \]

이렇게 보면 조금 더 명확하죠? 우리가 지겹도록 보았던 \(y = \alpha + \beta_1x_1 + \cdots + \beta_kx_k + \epsilon\)과 형태가 같습니다. 그런데 이렇게 길게 쓰는 건 너무 비효율적이니 계수값들의 벡터를 \(\beta = (\alpha\;\;\beta_k\;\;\cdots\;\; \beta_k)'\)라고 나타내고 예측변수들을 \(\textbf{X} = (1\;\;x_1\;\;\cdots\;\; x_k)\)라고 표현합니다. 이때 \(\textbf{X}\)는 \(N\times(K+1)\)의 행렬이며, \(\beta\)는 \(1\times(K+1)\)의 벡터입니다. 이해하기는 어렵지 않겠죠? 왜냐면 예측변수들이야 변수들의 개수(\(K\))에 절편을 표현할 항이 하나 더 추가되서 (\(K+1\))에 실제 관측된 관측치의 수(\(N\))을 고려해준 것이고, 계수값은 변수들에 절편을 포함한 개수만큼 계산되는데 위에서 살펴본 것처럼 열로 표현되니까 벡터라고 할 수 있는 것입니다. 그렇다면 우리는 행렬로 나타내는 OLS를 다음과 같이 심플하게 보여줄 수 있습니다.

\[ \textbf{y} = \textbf{X}\beta+\epsilon \]

이 하나의 식이 모든 \(N\)개의 관측치와 \(K\)개의 예측변수들에 대한 정보를 함축하고 있는 것입니다.

행렬 OLS에서의 가우스-마르코프

당연히 행렬 OLS도 우리가 이전 챕터들에서 보았던 OLS와 동일하니만큼 OLS의 기본가정, 가우스-마르코프 가정을 \(\textbf{y} = \textbf{X}\beta+\epsilon\)을 이용해 보여줄 수 있습니다. 먼저 우리의 추정치로부터 잔차의 벡터를 얻어냅니다. 그냥 식의 좌변과 우변을 요리조리 옮겨주면 됩니다.

\[ \hat{\epsilon} = \textbf{y}-\textbf{X}\hat{\beta} \]

개별 종속변수의 관측치에서 우리가 추정한 계수—모델에 일련의 예측변수들의 관측치들을 하나하나 대입하여 계산된 예측값(\(\hat{\textbf{y}} = \textbf{X}\hat{\beta}\))을 제해주면 나오는 것이 우리가 모델로 설명하지 못한 종속변수의 변동(variations), 오차(errors)가 될 것입니다. 물론 표본을 가지고 하는 것이니 잔차(residuals)라고 하는 게 정확한 표현입니다. 그리고 우리는 이렇게 얻은 잔차의 제곱합을 최소한으로 하고 싶은거죠. 이게 OLS의 핵심이니까요. 최소로 만드는 값을 찾는 것은 여러 가지 방법이 있겠습니다만, 여기서는 편미분을 통해 제곱합이 0이 되는 해를 구하는 방식을 취하겠습니다.

\[ \frac{\partial \hat{\epsilon}'\hat{\epsilon}}{\partial\hat{\beta}} = \frac{\partial(\textbf{y}-\textbf{X}\hat{\beta})'(\textbf{y}-\textbf{X}\hat{\beta})}{\partial \hat{\beta}} = 2\textbf{X}'\textbf{X}\hat{\beta}-2\textbf{X}'\textbf{y} = 0 \]

여기서 \(\hat{\epsilon}'\hat{\epsilon}\)은 제곱합, \(\sum^{N}_{i=1}\epsilon^{2}_{i}\)를 벡터로 보여준 것입니다. 아무튼 \(2\textbf{X}'\textbf{X}\hat{\beta}-2\textbf{X}'\textbf{y} = 0\)라고 할 때, 우리는 이 식을 다시 \(\hat{\beta} = (\textbf{X}'\textbf{X})^{-1}\textbf{X}'\textbf{y}\)로 나타낼 수 있습니다. 이것이 회귀계수(물론 절편까지!)를 구하는 보다 간명한 해(solution)가 되겠습니다. 그리고 이때 회귀계수의 분산은 \(Var(\hat{\beta}) = \sigma^2(X'X)^{-1}\)가 됩니다. \(Var(\hat{\beta})\)는 \((K+1)\times(K+1)\)의 차원을 갖는 행렬이 되고, 이를 우리는 분산-공분산행렬(variance-covariance matrix)이라고 합니다. 이 분산-공분산행렬은 무엇이고, 왜 중요한지를 이제부터 살펴보겠습니다.

표집분포 (Sampling distribution)

지겹게 들어왔던 표집분포로 다시 돌아왔습니다. 그런데 이번에는 조금 다른 방식으로 표집분포를 나타내보겠습니다. 행렬에 대해서 살펴보았으니 표집분포를 행렬의 형태로 나타내 보겠습니다.

\[ \hat{\beta}\sim N(\beta, \hat{\sigma}^2(\textbf{X}'\textbf{X})^{-1}) \]

수없이 많은 표본을 뽑아서 그 표본들에 동일한 선형회귀모델을 적용해 계산한 계수값 \(\hat{\beta_k}\)들은 표집에 문제가 없다면 모집단의 모수, \(\beta\)를 평균으로 하고 설명변수들 간에 나타날 수 있는 분산-공분산, 일종의 표집과정의 본연적 한계로 나타나는 오차에 따라 분산(\(\hat{\sigma}^2(\textbf{X}'\textbf{X})^{-1}\))을 가지게 될 것이라는 이야기입니다. 이때, 우리는 모집단의 분산을 알 수 없으므로 \(\sigma^2\) 대신에 표본의 분산, \(\hat{\sigma}^2\)를 사용합니다.

행렬 OLS의 유용성

행렬 OLS는 어떠한 장점을 가지고 있길래 사용하는 걸까요? 무엇보다도 예측값에 있어서 분석적 접근(analytical approach)보다 유용하다고 할 수 있습니다. 분석적 접근이란 쉽게 말해서 통계 이론적으로 OLS에 접근하는 방식을 의미하는데요, 모집단과 표본의 관계를 바탕으로 오차와 잔차의 관계, 계수값과 영가설을 바탕으로 t와 p의 값을 계산하고, 표준편차와 관측치의 수를 바탕으로 표준오차를 계산하는 등의 접근법을 말합니다. 이제까지 해오던 방식입니다.

반면에 행렬 OLS는 우리가 직접적으로 관측한 데이터를 바탕으로 모델을 통해 일종의 예측값들을 얻을 수 있도록 해줍니다. 정확하게는 우리가 관심을 가지고 있는 행렬 \(\bar{\textbf{X}}\)를 가지고 모델을 거쳐 나올 예측값을 계산할 수 있게 됩니다. 다시 말하면, \(\bar{\textbf{X}}\)는 우리가 임의로 만들 수 있습니다. 이 \(\bar{\textbf{X}}\)는 관측된 실제 변수 \(\text{X}\)와 마찬가지로 모든 \(K\)개의 변수들에 대한 가정을 포함하고 있습니다. 즉, 원래 변수는 1, 4, 11, 5, 6 이런 관측치의 배열을 가지고 있다고 할 때, 우리는 이 변수가 1부터 11까지 1씩 증가할 때 나타날 \(\textbf{y}\), 예측값 \(\hat{\textbf{y}}\)를 얻을 수 있습니다. \(\hat{\textbf{y}} = \bar{\textbf{X}}\hat{\beta}\)라는 과정을 통해서 우리는 변수들의 관계를 쉽게 이해하고 실질적인 의미를 파악할 수 있습니다. 그리고 무엇보다 이 \(\hat{\textbf{y}} = \bar{\textbf{X}}\hat{\beta}\)는 R을 통해 구현될 수 있습니다.

행렬 OLS, R로 보기

행렬 OLS가 과연 분석적 접근법으로 계산해낸 OLS 추정 결과와 같은 결과를 가질지를 알아보기 위해서 실제의 데이터를 이용한 R 분석을 수행해보겠습니다. 언제나와 같이 QoG 데이터셋에서 2015년 기준의 POLITY2, 무역개방성, 2010년 기준 1인당 GDP를 따로 서브셋을 만들어 사용할 것입니다. 먼저 모델이 어떻게 생겼는지를 보겠습니다.

# 먼저 model1이 우리가 추정할 모델입니다.
## 종속변수는 로그를 취한 1인당 GDP
## 1은 절편을 의미
model1 <- lm(log2(wdi_gdpcapcon2010) ~ 
               1 + p_polity2 + wdi_trade, data=QOG2) 

위의 모델은 로그값을 취한 1인당 GDP를 종속변수로 하고, POLITY2 즉, 민주주의와 무역개방성을 예측변수로 합니다. 민주주의 수준과 무역개방성을 가지고 1인당 GDP를 설명하고자 하는 모델인 것이죠. 그럼 행렬 OLS로 이 모델을 재구성해보겠습니다.

X <- model.matrix(object = log2(wdi_gdpcapcon2010) ~ 1 +  
                    p_polity2 + wdi_trade, data=QOG2)
head(X)

##   (Intercept) p_polity2 wdi_trade
## 1           1        -1  55.80373
## 2           1         9  71.80102
## 3           1         2  59.73290
## 4           1        -2  70.32683
## 5           1        -7  72.60151
## 6           1         9  22.48623

먼저 우리가 필요로 하는 예측변수들의 행렬, \(\textbf{X}\)를 만들어줍니다. 로그화된 1인당 GDP의 개별 관측치에 대응하는 절편과 실제 관측된 민주주의, 무역개방성의 관측치가 행렬의 형태로 저장됩니다.

Y <- log2(QOG2$wdi_gdpcapcon2010)
head(Y)

## [1]  9.276256 12.143589 12.216623 11.864574 12.578616 13.357908

마찬가지로 로그값을 취한 1인당 GDP도 \(\textbf{Y}\)라는 행렬에 저장해줍니다. 벡터도 행렬이니까요. 이제 저 둘을 계산해보겠습니다.

beta.Hat <- solve(t(X) %*% X) %*% t(X) %*% Y

이게 우리가 수식으로 썼던 \(\hat{\beta} = (\textbf{X}'\textbf{X})^{-1}\textbf{X}'\textbf{y}\)에 대응하는 R 코드입니다. 순서대로 한 번 보겠습니다. 결과가 어떻게 변화하는지 보세요.

step1 <- solve(t(X) %*% X) 
step1

##               (Intercept)     p_polity2     wdi_trade
## (Intercept)  0.0286305657 -7.398553e-04 -2.205416e-04
## p_polity2   -0.0007398553  1.843584e-04 -9.382902e-07
## wdi_trade   -0.0002205416 -9.382902e-07  2.629664e-06

step2 <- solve(t(X) %*% X) %*% t(X)
step2 %>% as.data.frame() %>% dplyr::select(1:4) %>% head() 

# 이런게 153개 열이 있습니다
step3 <- solve(t(X) %*% X) %*% t(X) %*% Y
step3

##                    [,1]
## (Intercept) 10.67984565
## p_polity2    0.08517992
## wdi_trade    0.01423212

자, 위의 저 결과가 바로 행렬 OLS로 구한 \(\hat{\beta}\)가 됩니다. 그럼 이번에는 표준오차를 구해볼까요? 이미 위에서 수식으로 다 보여드린 내용을 R로 다시 풀어놓은 것과 다름이 없습니다. 여기서 처음보는 건 아마 chol 함수일텐데요, 이는 촐레스키 분해(Cholesky Factorization) 함수입니다.⁶ 저는 여기서 자세한 설명은 하지 않을 것이기 때문에 참고할만한 블로그 링크를 공유하겠습니다.

sig.sq <- sum((Y - X%*%beta.Hat)^2)/(nrow(X)-ncol(X)) ## 분산이죠?
VCV <- sig.sq*chol2inv(chol(t(X)%*%X))
SE <- sqrt(diag(VCV))              
# 행렬 OLS 결과
cbind(beta.Hat, SE) %>% knitr::kable()  

# R lm()을 이용한 OLS 결과
coef(summary(model1)) %>% knitr::kable()

두 결과의 계수값과 표준오차 추정치가 일치하는 것을 확인하실 수 있습니다.