Diagnostics

여기서는 먼저 선형회귀(linear regression) 모델에 대해 간단하게 리뷰하고 어떤 경우에 선형회귀분석을 사용하는 것이 바람직한지를 진단하는 기준들에 대해 살펴보도록 하겠습니다. 그리고 나아가 간단하게 종속변수가 연속형이 아닌, 0과 1일 경우에 선형회귀모델을 적용한 경우, 선형확률모형(linear probability models)에 대해서도 살펴보도록 하겠습니다.

  • 선형확률모형을 통해 우리는 종속변수가 연속형이 아닐 때, 즉 예측변수와 종속변수 간 선형 관계를 가정할 수 없을 때, OLS를 사용하는 것이 왜 ‘비효율적’(inefficient)이고 그 추정치가 BLUE1일 수 없는지를 알 수 있습니다.

필요한 패키지를 로드하기

library(ezpickr)   # 다른 확장자의 파일을 R로 불러오기 위한 패키지
library(here)      # 현재 작업디렉토리를 R-스크립트가 위치한 디렉토리로 자동설정하는 패키지
library(broom)     # 통계분석 함수의 결과를 정연하게 정리해주는 패키지
library(ggfortify) # ggplot2 시각화를 도와주는 패키지
library(estimatr)  # 로버스트 표준오차나 신뢰구간 등과 같은 추정치들을 계산해주는 패키지
library(patchwork)
library(futurevisions) # 시각화를 위한 컬러패키지
library(tidyverse) # 데이터 관리 및 전처리를 위한 주요 패키지
library(showtext)
library(extrafont)
loadfonts()
font_add("NanumGothic", "NanumGothic.ttf")
theme_set(
  theme_bw() +
    theme(text = element_text(family = "NanumGothic"))
)
Sys.setlocale("LC_CTYPE","korean")
## [1] "Korean_Korea.949"

선형회귀모델의 간단한 리뷰

선형회귀모델은 몇 가지 가정들에 바탕을 두고 있습니다. 이 가정들을 통틀어 가우스 마르코프 가정(Gaus-Markov assumptions)라고 합니다. 여기서는 편의상 다중선형회귀모델의 가정이라는 의미로 MLR(Multiple linear regression)의 가정, MLR1부터 MLR5로 표현하겠습니다.

  • MLR1. 선형회귀분석은 다차항(polynomial terms) 등을 포함하여 비선형성(non-linearity)도 모델에 반영할 수 있지만 기본적으로는 종속변수와 예측변수 간의 관계가 선형일 것 등을 가정합니다. 여기서 언급한 선형회귀모델이 표현할 수 있는 비선형성의 문제는 이 다음 주차의 상호작용항에 관한 내용에서 자세히 살펴보실 수 있습니다.

  • MLR2. 잔차(residuals)는 정규적으로 분포되어 있어야 합니다.—잔차의 정규성

  • MLR3-4. 데이터는 iid 가정(assumed indendently and identically distributed)을 충족시켜야 합니다.

    • MLR3. 데이터의 각 확률변수는 상호독립적이며(independent),

    • MLR4. 모두 동일한 확률분포를 가지고(등분산성, identical) 정규분포를 따른다(normally distributed)는 가정입니다. 즉, 데이터의 각 확률변수들과 잔차는 서로 상관하지 않아야 합니다.

    • MLR5. 등분산성(homoscedasticity)

      • 모든 잔차(오차)는 같은 분산을 가진다는 가정입니다.

이처럼 선형회귀모델은 강한 가정들을 위에 성립합니다. 또한 OLS, 제곱합 최소화 방법을 통한 추정치는 자료에 민감하다는 특징을 가집니다. 소수의 극단적인 이탈치들이 잔차의 제곱합을 크게 변화시킬 수 있기 때문입니다. 따라서 우리에게는 두 가지 선택지가 존재합니다.

  1. OLS를 통한 추정을 포기하는 것

  2. OLS로 자료를 분석하되, 데이터를 “진단”(diagnose)해보는 것

Fox (2016, 266)은 회귀모델을 수립하기에 앞서 데이터에 대한 주의 깊은 분석을 통해 여러 문제들을 예상하고 미리 대처할 수 있다고 밝히고 있습니다. 대표적인 문제들로는,

  • 이탈치 (outliers)

  • 정규적으로 분포되지 않은 잔차(non-normally distributed errors)

  • 일정하지 않은 분산(non-constant variance)

  • 비선형성(non-linearity)

  • 공선성(collinearity)

등이 있다고 할 수 있습니다.

이탈치

이탈치(outliers)는 설명변수의 값 중 다른 값들과는 달리 유독 튀는, 극단적인 값을 갖는 관측치를 의미합니다. 단, 단변량 분석에서의 이탈치와 양변량, 혹은 다변량에서의 이탈치는 차이가 있습니다.

그리고 이탈치는 회귀모델의 추정에 큰 영향을 미칠 수 있습니다.

회귀선의 적합도(fit)에 영향 (Influence)을 미치기 때문에, 이탈치와 같이 “이례적인” 데이터는 제곱합(least square)에 문제가 될 수 있습니다.

  • \(\text{Influence} = \text{Leverage}\times\text{Discrepancy}\)

    • \(\text{Leverage} = X\text{에서 관측치들이 떨어져 있는 정도}\)

      • \(h_i\)는 레버리지를 측정하는 대표적인 지표로, \(x_i\)와 \(\bar{x}\) 간의 거리를 보여줍니다.

      • \(h_i = \frac{1}{n}+\frac{(X_i-\bar{X})^2}{\sum^{n}_{j=1}(X_j-\bar{X})^2}\)

    • \(\text{Discrepancy} = Y\text{에서 관측치들이 떨어져 있는 정도}\)

      • 큰 레버리지를 가진 관측치는 작은 잔차를 가지는 경향이 있습니다.

      • \(V(E_i) = \sigma^{2}_{\epsilon} (1-h_i)\)

      • 표준화된 잔차: \(E_i = \frac{\sigma^{2}_{\epsilon} (1-h_i)}{S_E\sqrt{1-h}}\)

      • 스튜던트화 잔차: \(E_i^{*} = \frac{\sigma^{2}_{\epsilon} (1-h_i)}{S_E(-i)\sqrt{1-h}}\)

        • 스튜던트화 잔차 = 표준화 잔차에서 \(i\)번째 관측치를 제외한 결과

        • 특정한 관측치가 아니라 1부터 \(n\) 사이의 아무 관측치를 제외하고 \(E_i\)와 \(E_i^{*}\)의 비교

      • 가설검정의 문제처럼 이탈치를 식별(outlier identification)하여, 다른 관측치들과 유의미한 차이의 거리를 가진 값을 관측치로 “결정”하는 기법입니다.

  • \(\text{Influence} = \text{hat-value}\times\text{residuals}\)

이같은 이탈치의 영향력을 측정하는 대표적인 지표는 쿡스의 거리(Cook’s distance)가 있습니다.

\[ \frac{E^{\prime2}}{k+1} \times \frac{h_i}{1-h_i} = \text{Residual}\times\text{Leverage} \]

  1. 한번에 하나씩 관측치를 표본에서 제거합니다.

  2. \(n-1\)의 표본을 대상으로 회귀모델을 다시 분석합니다.

  3. \(i\)번째 표본이 제거되었을 때의 전체 예측값이 얼마나 변화하였는지를 분석합니다.

쿡스의 거리는 평균적으로 데이터셋에서 문제가 될 것이라고 생각되는 관측치가 제거되었을 때, 예측된 \(y\)의 값이 어떻게 변할지를 살펴보는 지표입니다.

비정규적으로 분포된 오차

Fox (2016, 297)에 따르면 “최소제곱합 추정법의 타당성(validity)이 강력하다 할지라도 그 효율성(efficiency)은 그렇지 않”습니다. 여러 개의 정점을 가진 오차의 분포(multimodal error distribution)는 데이터를 몇 가지 집단으로 나누는 설명변수가 누락되었을 가능성을 시사합니다. 즉, 데이터에 체계적인 영향을 미치는 변수를 빼먹었을 수 있다는 신호입니다.

오차의 비정규성을 검정하기 위해서 일단 여기서는 단변량 관계를 살펴보도록 하겠습니다. 분위비교플롯(quantile-cimparison plot)은 스튜던트화 잔차의 표본 분포를 단위정규분포의 분위와 비교 및 시각화하는 플롯입니다. 이른바 Q-Q 플롯을 통해 우리는 잔차 분포의 꼬리 부분이 어떻게 나타나는지를 효과적으로 보여줄 수 있습니다.

일정하지 않은 오차의 분산

일정한 오차의 분산을 동분산성 혹은 등분산성(homoscedasticity)라고 하며, 일정하지 않은 분산을 이분산성(heteroscadasticity)라고 합니다. OLS 추정치의 결과가 가지는 효율성과 타당성은 (1) 표본 크기, (2) 변동량(variations)의 정도, (3) \(X\) 값이 어떻게 이루어져 있는지, 그리고 (4) 오차 분산과 \(X\)의 관계 등에 의해 영향을 받습니다. 당연히 표본의 크기가 크면 클수록 트렌드가 존재한다면 그 패턴이 더 선명하게 나타날 것이니 보다 정확하게 추정할 수 있을 것이구요, 변수의 변화량이 많을수록 설명할 수 있는 것도 많을 것입니다. 이분산성은 오차의 분산이 일정하지 않아 \(X\)와 체계적인 관계를 맺을 수 있는 확률이 더 높다는 점에서 OLS 결과에 영향을 미칠 수 있습니다.

이분산성을 확인하기 위해서는 잔차 플롯(residual plot)을 살펴보아야 합니다. 대개 스튜던트화된 잔차와 예측값 간의 관계를 보거나, 혹은 절대값을 취한 스튜던트화 잔차와 예측값의 관계를 살펴봅니다. 논리는 간단합니다. 예측값은 우리가 이미 가지고 있는 \(X\)로 \(Y\)를 추정한 결과입니다. 잔차는 그 예측값과 실제 관측치 간의 차이를 보여주죠. 만약 예측값과 잔차 간에 체계적 관계가 존재한다면? 우리는 \(Y\)를 설명하기 위한 일련의 설명변수, \(X\)에 적절한 변수를 다 포함시키지 못한 것입니다. 통계적으로는 가장 큰 오차의 분산이 가장 작은 잔차의 분산보다 10배 이상 클 때 심각할 정도의 결과의 편향이 나타날 수 있다고 봅니다 (보수적으로는 4배를 기준으로 삼기도 함).

비선형성

\(Y\)와 \(X\) 간의 관계가 비선형일 때, 선형성을 가정한 선형회귀모델은 결과를 제대로 있습니다. 선형회귀모델의 회귀계수는 결과적으로 \(X\)의 한 단위 변화가 \(Y\)에 미치는 한계효과를 보여주는 것인데, 이 한계효과는 비선형적 관계를 단편적으로밖에 보여주지 못합니다. Fox (2016)은 성분잔차플롯(component-plus-residual plot)을 통해 변수들 간 관계의 비선형성 여부를 확인할 것을 주문합니다.

성분잔차플롯은 \(e + \hat{\beta}_{j} X_{j} \sim X_{j}\)의 관계를 보여주는 산포도입니다. \(\hat{\beta}_{j} X_{j}\)은 \(j\)번째 설명변수가 예측값에 기여하는, 계수값에 개별 설명변수의 관측치를 곱한 값입니다. 즉, 예측값에 오차를 더한 값과 개별 설명변수 간의 관계를 보여줌으로써, 만약 설명변수가 선형관계에 있다면 플롯은 \(\hat{beta}\)라는 기울기에 밀접하게 모인 관측치의 분포를 보여줄 것입니다. 비선형적이라면 선을 통과하기는 하지만 모여있지는 않은 휘어진 분포 등을 보여줄 것입니다.

공선성

공선성은 다른 문제들에 비해서는 상대적으로 문제시되지는 않습니다. 현실 세계에서 우리가 관측할 수 있는 통계치들, 표본의 값들이 공선성에서 자유로울 수는 없는데 만약 그렇다면 공선성이 어느 정도여야 선형회귀모델의 추정치가 편향되어 신뢰하기 어렵다고 판단할 수 있을까요?

먼저 설명변수들 간에 완벽한 공선성이 존재한다고 생각해보겠습니다. 이 경우, 최소제곱합으로 만들어진 모델은 유일한 해(solution), 계수를 가지지 못합니다. 왜냐하면 계수값에 대한 표집 분산이 무한대가 되기 때문입니다. 단순하게 설명하자면, \(Y\)를 설명하기 위한 각 \(X\)의 변량이 모두 동일하게 겹치므로 중복…중복… 그 결과 무한대나 다름 없는 해가 생기는 것입니다.

만약 공선성이 완전하지는 않지만 꽤 높은 수준으로 존재한다면 어떻게 될까요? 계수(_j$에 대한 표집 분산은 다음과 같이 추정되는데,

\[ Var(\beta_j) = \frac{1}{1-R^2_j}\times\frac{\sigma^2_{\epsilon}}{(n-1)S^2_j} \]

Fox (2016, 342)의 Figure 13.1에서 확인할 수 있듯이, “\(R_j\)가 0.9까지 도달하지 않는 한, 추정치의 정확성(precision)은 훼손되지 않”습니다.

다중선형회귀모델

선형회귀모델을 진단하기에 앞서 간단한 모델을 하나 만들어보도록 하겠습니다. ggplot2 패키지에서 제공하는 diamonds 데이터로 다중선형회귀모델을 만들도록 하겠습니다. 이 예제에서 사용할 종속변수는 바로 다이아몬드 데이터의 가격(price) 변수입니다. 따라서 가격 변수의 분포를 살펴보겠습니다.

df <- diamonds

df %>% 
  ggplot(aes(price)) + geom_histogram()

일단 전체 다이아몬드 데이터에서 가격 변수의 분포는 좌측으로 치우친(left-skewed) 것을 확인할 수 있습니다. 즉, OLS의 기본 가정에서 언급되었던 종속변수가 정규적으로 분포되어 있을 것이라는 가정이 위배된다 충분히 의심해볼 수 있습니다.

자, 이제는 이 df 데이터의 변수들을 필요에 따라 조작해주도록 하겠습니다.

  • sold의 경우는 예제이니까 30퍼센트의 확률로 팔린다 (1)가 나타나도록 이항분포에서 무작위 추출을 한 변수입니다.

  • price_bin은 원래 숫자형 변수인 price를 조건에 따라 low, medium, high의 레이블을 갖는 요인형 변수로 조작하였습니다.

df <- df %>%
  mutate(
    sold = rbinom(nrow(df), 1, 0.3),
    price_bins = case_when(    # ifelse와 비슷한 조건문
      price <= 2500L ~ 'low',  # price가 2500 이하이면 'low'라고 레이블
      price > 2500L & price <= 7500L ~ 'medium', # 2500 초과 7500 이하면 'medium'
      price > 7500L  ~ 'high', # 7500 초과면 'high'
      TRUE ~ NA_character_     # 위의 세 조건을 만족시키지 못하는 관측치는 NA 처리
    ) %>% parse_factor(., levels = c('low', 'medium', 'high'), include_na = F)
  ) # 그리고 이렇게 분류된 결과를 요인형 변수로 변경, NA는 제외

분석하게 용이하게 변수의 순서도 조금 손을 봐주도록 하겠습니다.

df$index <- 1:nrow(df) # df 데이터에 행번을 표시하는 변수 추가
df <- df %>% select(
  ## 행을 식별하는 row_id 추가
  row_id = index, sold, carat, cut, color, clarity, depth, 
  table, price, price_bins, x, y, z
  ## row_id, sold, carat, cut, color, clarity, depth, table, 
  ## price, price_bins, x, y, z 순으로 변수의 순서를 변경.
)
head(df)

이번에는 데이터 안에서 가격이 high인 다이아몬드가 몇 개나 있을지 세어보겠습니다. 이미 price_bin이라는 변수를 새로 만들어서 숫자형인 price가 얼마를 기준으로 초과할 때 high라고 할 수 있는지를 지정해 두었습니다.

df %>% count(price_bins)

그렇다면 각각의 가격 범주에서의 평균 가격은 어떻게 될까요?

df %>% select(price, price_bins) %>%
  group_by(price_bins) %>%
  summarise_all(mean, na.rm = T)

선형관계 살펴보기

여기까지가 대략적인 데이터의 기술분석이라고 한다면, 이제는 이번에는 다이아몬드의 무게(carot)와 가격 간의 관계를 살펴보도록 하겠습니다. 이번에는 통계치보다는 가시화된 플롯으로 그 결과를 구성해보도록 하겠습니다. x축에는 무게, y축에는 가격이 오도록 플롯을 그립니다.

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)) +
  theme_bw() +
  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 값이 커질수록 관측치가 옅어지는 것을 확인할 수 있습니다. 아마도 이차항(quadratic term)을 포함하는 것이 더 관계를 잘 보여줄 수 있을 것으로 보입니다.

  • \(y = \alpha + \beta_1x_1 + \epsilon\)과 같은 선형 모델에서

  • \(y = \alpha + \beta_1{x_1}^2 + \epsilon\)과 같은 관계일 것이라고 상정해보는 것이죠.

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]) + 
  # 2차항을 포함하였을 때 변화한 회귀선
  geom_smooth(method = 'lm', formula = y ~ x + I(x^2), 
              se = F, color = futurevisions("mars")[4]) + 
  coord_cartesian(ylim = c(-1000, 22000)) +
  theme_bw() +
  labs(
    x = 'Carat of Diamond', y = 'Price of Diamond',
    title = "Effect of Carat on Diamond Price"
  ) + 
  scale_y_continuous(labels = scales::dollar_format())

어떤 것이 데이터를 더 잘 보여주는 것 같나요? 보라색 회귀선, 이차항으로 그 관계를 표현해준 것이 좀 더 carat과 price의 관계를 잘 보여주는 것으로 보이네요. 그럼 이번에는 다이아몬드 데이터에 관측치가 너무 많기 때문에 딱 1,000의 배수가 되는 관측치들만 뽑아서 살펴보도록 하겠습니다.—천번째, 이천번째, 삼천번째…

temp_df <- df %>%
  dplyr::filter(row_id %% 1000 == 0) # %%는 나머지를 구하는 연산자이다.
glimpse(temp_df)
## Rows: 53
## Columns: 13
## $ row_id     <int> 1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000, 9000, 10...
## $ sold       <int> 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0,...
## $ carat      <dbl> 1.12, 0.72, 0.90, 0.74, 1.05, 1.01, 1.00, 1.01, 0.91, 1....
## $ cut        <ord> Premium, Ideal, Good, Ideal, Very Good, Premium, Premium...
## $ color      <ord> J, D, G, D, I, H, F, G, E, D, D, D, G, H, F, I, E, H, J,...
## $ clarity    <ord> SI2, VS2, SI2, SI1, SI2, SI2, SI2, SI2, SI1, SI1, VS2, S...
## $ depth      <dbl> 60.6, 61.8, 63.8, 61.8, 62.3, 62.8, 62.6, 61.9, 63.5, 64...
## $ table      <dbl> 59, 57, 56, 56, 59, 61, 58, 59, 57, 59, 59, 62, 60, 59, ...
## $ price      <int> 2898, 3099, 3303, 3517, 3742, 3959, 4155, 4327, 4512, 47...
## $ price_bins <fct> medium, medium, medium, medium, medium, medium, medium, ...
## $ x          <dbl> 6.68, 5.76, 6.13, 5.84, 6.42, 6.33, 6.37, 6.46, 6.11, 6....
## $ y          <dbl> 6.61, 5.73, 6.16, 5.88, 6.46, 6.25, 6.31, 6.39, 6.14, 6....
## $ z          <dbl> 4.03, 3.55, 3.92, 3.62, 4.01, 3.95, 3.97, 3.98, 3.89, 4....

새로운 데이터셋, temp_df의 row_id를 보면 천의 자리의 배수로 관측치들이 추출되어 있다는 것을 확인할 수 있습니다. 이제, 이 데이터를 가지고 선형회귀모델을 다시 만들어보겠습니다.

# 종속변수는 price, 예측변수는 carat
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 패키지를 이용하면 baseR의 summary 패키지보다 조금 더 세련된 방식으로 원하는 결과를 얻을 수 있습니다.

## 선형모델로부터 예측값과 그 외 기타 통계치들을 fitted_values 라는 객체에 저장
fitted_values <-  broom::augment(model)  
head(fitted_values)

.fitted는 각 관측치에 따른 예측값으로 주어진 선형회귀모델에 개별 price와 carat, 즉 \(\text{price}_i\)와 \(\text{carat}_i\)를 대입하여 얻어낸 \(\hat{Y_i}\)입니다. 그리고 .resid는 실제 종속변수의 관측치, \(Y_i\)와 예측값의 차이인 잔차를 보여줍니다.

fitted_values %>%  
  # 각 1,000번째 관측치로 필터링하여 총 53개의 값을 가지는 temp_df를
  ggplot(aes(x = carat, y = price)) + 
  # 기준으로 종속: price/예측: carot 
  geom_point() +                     
  geom_smooth(method = 'lm', se = F, color = futurevisions("mars")[3]) +
  geom_segment( # 수직선을 그리는 옵션
    aes(x = carat, y = price, xend = carat, yend = .fitted), 
    # xend는 carat, yend는 예측값 
    color = futurevisions("mars")[1], size = 0.3 
    # 즉, 관측치로부터 예측값에 수직인 선을 긋는 코드
  ) +
  geom_curve(  # 커브 화살표를 그리는 옵션
    aes(x = 2.75, y = 10000, xend = 2.55, yend = 15000), 
    color = futurevisions("mars")[1], # 커브 명령어
    ## 시작인 x가 2.75, y가 10000인데 끝은 2.55, 15000
    ## 즉, x는 -0.2만큼, y는 +5000만큼 좌상 방향의 곡선을 그리라는 명령어
    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())

이 문서는 기본적으로 OLS 회귀분석—최소자승법을 이용한 선형회귀분석에 대한 개념을 이해하고 있다는 전제를 바탕으로 작성되었습니다. 이 그림은 carat과 price 간의 관계를 보여줄 때, 선형회귀선을 긋는다는 것이 어떠한 의미인지를 보여줍니다.

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

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

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

    • OLS는 각 점들로부터 수직거리(즉, 실제 관측치와 예측값 간의 차이), 잔차(residuals)의 제곱 합이 가장 적은 선을 긋는 기법입니다.

      • 모집단에 대해서는 오차(errors)지만, 표본에 적용할 때에는 잔차(residuals)라고 표현합니다.

    • 따라서 Ordinary least square란 결국 잔차의 제곱합이 최소가 되는 선을 긋는 방법을 말합니다.

    • 제곱을 하는 이유는 가정 상 잔차의 총합은 0(\(\sum_{i=1}^n\epsilon_i = 0\))이기 때문에 잔차의 총합을 그냥 구하게 되면 서로 상쇄되어 0이 되므로 잔차의 총합이 가장 작아지는 선을 긋기 어렵습니다.

그렇다면 이번에는 다중선형회귀모델로 기존의 단순선형회귀모델을 확장해보도록 하겠습니다.

full_model <- lm(
  price ~ carat + cut + color + clarity + depth + x + y + z, 
  data = temp_df
)
summary(full_model) # BaseR의 summary() 함수
## 
## Call:
## lm(formula = price ~ carat + cut + color + clarity + depth + 
##     x + y + z, data = temp_df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1668.2  -540.8  -132.6   330.2  2610.6 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  39476.37   49735.10   0.794  0.43320    
## carat        12483.92    1742.57   7.164 3.92e-08 ***
## cut.L          964.37     605.16   1.594  0.12086    
## cut.Q         -922.27     527.28  -1.749  0.08986 .  
## cut.C          427.08     441.22   0.968  0.34033    
## cut^4          611.97     374.80   1.633  0.11232    
## color.L      -2199.29     493.21  -4.459 9.49e-05 ***
## color.Q       -548.75     421.33  -1.302  0.20207    
## color.C         82.06     456.69   0.180  0.85853    
## color^4        -13.17     451.26  -0.029  0.97690    
## color^5       -317.53     434.29  -0.731  0.47001    
## color^6       -519.72     399.58  -1.301  0.20266    
## clarity.L    14308.21    1974.52   7.246 3.12e-08 ***
## clarity.Q   -10561.87    1712.62  -6.167 6.70e-07 ***
## clarity.C     7119.07    1187.18   5.997 1.10e-06 ***
## clarity^4    -3150.77     607.94  -5.183 1.17e-05 ***
## clarity^5     1405.62     404.47   3.475  0.00149 ** 
## depth         -660.91     798.16  -0.828  0.41378    
## x              188.06    5272.89   0.036  0.97177    
## y            -7612.02    5874.17  -1.296  0.20430    
## z             9615.52   13534.97   0.710  0.48259    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 971.5 on 32 degrees of freedom
## Multiple R-squared:  0.965,  Adjusted R-squared:  0.9431 
## F-statistic: 44.13 on 20 and 32 DF,  p-value: < 2.2e-16
tidy_model <- broom::tidy(full_model) %>% 
  mutate_if(is.numeric, ~ round(., 2))  # 분석 결과를 티블로 나타냅니다.
tidy_conf <- broom::confint_tidy(full_model) # 신뢰구간을 추정
tidy_conf <- tidy_conf %>%
  ## 신뢰구간을 나타내기 위한 변수들의 이름을 변경
  rename(ll = conf.low, ul = conf.high) 

# 한 번에 티블로 보여줄 수 있도록 tidy_model과 tidy_conf의 열을 결합
tidy_model <- bind_cols(tidy_model, tidy_conf) %>%
  mutate_if(is.numeric, ~ round(., 2)) # 자, 결과는?

만들어진 티블은 선형회귀분석 결과를 잘 정리하여 보여줍니다. 먼저 상수항(절편)을 포함한 변수들의 계수값(estimate)이 나타나있고, 이어서 표준오차(std.error)를 보여주고 있습니다. t값(statistic)을 바탕으로 유의수준(p.value)도 보여주고 있고, 신뢰구간을 위해 최저값(ll)과 최고값(ul)도 나타납니다.

그렇다면 이제는 이 모델이 선형회귀모델의 기본적인 가정들, Gaus-Markov 가정을 위해하는지를 살펴보겠습니다. 먼저 선형회귀모델의 결과를 broom 패키지의 augment() 함수를 이용하여 별도의 객체인 fitted_model에 저장하였습니다.

fitted_model <- augment(full_model)
head(fitted_model)

fitted_model을 보면 변수 중에 표준화된 잔차(.std.resid)가 있습니다. 이 변수를 히스토그램으로 나타내보겠습니다.

fitted_model %>%
  ggplot(aes(x = .std.resid)) + geom_histogram()

표준화된 잔차의 값이 넓게 분포하여 있다는 것은 우리의 예측값과 개별 관측치 간의 차이가 큰 경우가 존재한다는 것입니다. 회귀분석의 가정을 보다 자세히 살펴보기 위해서 다음과 같은 플롯을 그려보았습니다.

autoplot(full_model)

autoplot()은 ggplot2 패키지에서 주로 사용되는 함수들을 데이터에 맞게 적용시켜주는 ggfortify 패키지에 속한 함수입니다. 순서대로 보자면,

  1. 잔차와 예측값 간의 관계 (Residuals vs Fitted)

  • 이 플롯은 잔차와 예측값의 관계가 선형인지 아닌지를 보여주는 플롯입니다.

  • 선형회귀분석의 가정이 만족되기 위해서는 이 플롯은 0의 값을 갖는 수평선을 보여주어야 합니다.

  • 여기서는 심하지는 않지만 약간의 선형성을 보여주고는 있습니다.

  • 따라서 이 선형 모형에 제곱항(squared terms)을 포함해주는 등의 조치를 취해줄 필요가 있다는 것을 보여줍니다.

  1. QQ 플롯

  • QQ 플롯은 잔차가 정규분포를 따르고 있는지를 보여줍니다.

  • 점선을 따라 분포되어 있어야 정규분포인데, 이 데이터의 결과는 역시 잔차의 정규성 또한 위반의 소지가 있음을 보여줍니다.

  1. 잔차의 동분산성 확인(Scale-Location or Spread-Location)

  • 잔차의 분산의 동질성을 확인하기 위해 그리는 플롯입니다(등분산성).

  • 선들이 평평하게 확산되는 수평선이 등분산성을 보여주는 좋은 지표입니다.

  • 이 예제 데이터는 일정한 경향성을 보여주기 때문에 이분산성의 문제(heteroscedasticity)가 존재할 수 있음을 보여줍니다.

    • 이분산성은 로버스트 표준오차(robust standard errors)를 통해서 쉽게 해결할 수 있습니다.

    • 로버스트 표준오차란 표준오차를 계산할 때, 단위(units) 간 차이를 고려해 보정해준 표준오차를 의미합니다.

  1. 잔차 대 레버리지(Residuals vs Leverage)

  • 극단적인 값들이 분석에 포함/배제될 경우에는 회귀분석 결과에 영향을 미칠 수 있습니다.

  • 따라서 이 플롯은 그러한 극단적인 이상치들이 있는지를 확인하는 데 사용됩니다.

  • 우측 극단에 값들이 나타나는 것을 확인할 수 있습니다.

  • 확실히 몇몇 극단적인 이탈치들이 제거될 필요가 있음을 보여줍니다.2

연속형이 아닐 경우에 선형회귀분석을 쓰면 어떻게 될까?

그렇다면 종속변수가 이항변수(binary)일때는 선형회귀모형(OLS)을 사용할 수 없을까요? 그건 아닙니다. 사용할 ‘수는’ 있습니다. 종속변수가 이항변수일 때, 선형회귀보델을 사용할 경우, 우리는 그것을 선형확률모형(linear probability model, LPM)이라고 부릅니다. 그러나 LPM을 사용하는 데 있어서는 여러 가지 제약들이 존재합니다.

  1. LPM은 이분산성이 거의 확정적입니다.

  • 왜냐하면 종속변수는 0에서 1 중 하나의 값만 가지는 데 비해 설명변수는 연속형일 경우 큰 변량(variations)을 가지므로 당연히 오차항의 이분산성이 존재할 가능성이 커집니다.

  • 이 문제는 표본의 크기가 크면 로버스트 표준오차를 취함으로써 쉽게 해결해줄 수 있습니다.

  1. LPM을 통해 우리는 0…1의 확률을 넘어서는 예측을 할 수 있습니다.

  • 만약 단위 구간(0…1)을 넘어서는 예측확률이 없거나, 거의 없다면 LPM은 대개 불편(unbiased)하고 일관될 것으로 기대됩니다.

  • 하지만 선형회귀모델은 기본적으로 종속변수가 0과 1 사이에 ’갇혀있다’고 가정하지는 않습니다. 즉, 경우에 따라 예측값이 0과 1을 넘는, 실존할 수 없는 값을 가지게 될 수 있습니다. 이를 외삽의 문제(extrapolation)라고 합니다.

그렇다면 LPM을 사용하는 이유는 무엇일까요?

  1. 빅데이터 분석을 할 때, 더 빠르게 계산이 가능합니다.

  • 물론 속도 문제는 이제 컴퓨터 기술이 발전해서 큰 문제는 아닐 수 있습니다.

  • 그러나 아무리 컴퓨터가 좋아도 확률론이 아니라 최대가능도를 반복계산하게 되면 시간이 더 걸리는 것은 사실이기 때문에, 실무에서는 LPM이 종종 사용됩니다.

  1. 해석하기가 정말 쉽습니다.

  • 계수 값은 결국 x의 한 단위 증가와 관련된 y의 증가분을 보여주는 “한계효과”입니다.

  • LPM 역시 y를 일종의 확률로 이해한다는 것에 차이가 있을 뿐, 기본적인 OLS 해석의 결을 같이하기 때문에 직관적입니다.

그럼 한 번 LPM 모형을 만들어보겠습니다. 종속변수를 다이아몬드의 판매 여부(sold)로 설정하겠습니다. 이분산성은 거의 확실하기 때문에 표준오차를 로버스트 표준오차로 설정해주도록 하겠습니다(se_type = "HC1).

lmp <- lm_robust(
  sold ~ price + carat + cut + color + clarity + depth + x + y + z, 
  data = df, se_type = "HC1"
)

# 만약에 위의 se_type을 지정해주지 않으면? 계수값은 변하지 않지만 표준오차가 달라지게 됩니다.

non_lmp <- lm(
  sold ~ price + carat + cut + color + clarity + depth + x + y + z, 
  data = df
) 

# 결과를 한 번 비교하여 살펴보겠습니다.
lmp_tidy <- tidy(lmp) %>% select(1:5)
nonlmp_tidy <- tidy(non_lmp) %>% select(1:5)
lmp.table <- lmp_tidy %>% 
  mutate_if(is.numeric, ~ round(., 3))
knitr::kable(lmp.table, digits = 3)
term estimate std.error statistic p.value
(Intercept) 0.234 0.113 2.062 0.039
price 0.000 0.000 1.480 0.139
carat 0.007 0.028 0.252 0.801
cut.L -0.019 0.009 -2.141 0.032
cut.Q 0.010 0.007 1.292 0.196
cut.C -0.002 0.006 -0.400 0.689
cut^4 0.004 0.005 0.818 0.414
color.L -0.013 0.008 -1.698 0.090
color.Q 0.000 0.006 -0.055 0.956
color.C -0.006 0.006 -1.020 0.308
color^4 0.000 0.005 0.064 0.949
color^5 0.006 0.005 1.109 0.267
color^6 0.007 0.005 1.489 0.136
clarity.L -0.020 0.014 -1.393 0.164
clarity.Q 0.005 0.012 0.431 0.667
clarity.C -0.011 0.010 -1.079 0.280
clarity^4 0.010 0.008 1.307 0.191
clarity^5 0.001 0.006 0.148 0.882
clarity^6 0.000 0.006 0.072 0.942
clarity^7 0.005 0.005 1.078 0.281
depth 0.002 0.002 1.428 0.153
x -0.008 0.010 -0.829 0.407
y -0.006 0.002 -2.738 0.006
z -0.002 0.009 -0.165 0.869 
non_lmp.table <- nonlmp_tidy %>% 
  mutate_if(is.numeric, ~ round(., 3))
knitr::kable(non_lmp.table, digits = 3)
term estimate std.error statistic p.value
(Intercept) 0.234 0.118 1.981 0.048
price 0.000 0.000 1.484 0.138
carat 0.007 0.028 0.254 0.799
cut.L -0.019 0.009 -2.166 0.030
cut.Q 0.010 0.007 1.310 0.190
cut.C -0.002 0.006 -0.403 0.687
cut^4 0.004 0.005 0.820 0.412
color.L -0.013 0.008 -1.682 0.093
color.Q 0.000 0.007 -0.055 0.957
color.C -0.006 0.006 -1.015 0.310
color^4 0.000 0.005 0.063 0.949
color^5 0.006 0.005 1.108 0.268
color^6 0.007 0.005 1.488 0.137
clarity.L -0.020 0.014 -1.398 0.162
clarity.Q 0.005 0.012 0.433 0.665
clarity.C -0.011 0.010 -1.086 0.277
clarity^4 0.010 0.008 1.311 0.190
clarity^5 0.001 0.006 0.148 0.882
clarity^6 0.000 0.006 0.073 0.942
clarity^7 0.005 0.005 1.080 0.280
depth 0.002 0.002 1.354 0.176
x -0.008 0.013 -0.627 0.531
y -0.006 0.008 -0.726 0.468
z -0.002 0.014 -0.114 0.909

표준오차가 미묘하게 다른 것을 확인할 수 있습니다. 좀 더 가시적으로 보기 위해서는 이 결과를 autoplot()으로 살펴보도록 하겠습니다.

## 먼저 로버스트 표준오차를 취하지 않은 경우.
lm(
  sold ~ price + carat + cut + color + clarity + depth + x + y + z, data = df
) %>% autoplot()

## 계수값의 비교
round(lmp_tidy$estimate, 4) %>% unname()
##  [1]  0.2340  0.0000  0.0071 -0.0186  0.0096 -0.0025  0.0041 -0.0132 -0.0004
## [10] -0.0061  0.0003  0.0058  0.0070 -0.0199  0.0052 -0.0108  0.0103  0.0009
## [19]  0.0004  0.0053  0.0023 -0.0084 -0.0057 -0.0016
round(nonlmp_tidy$estimate, 4)
##  [1]  0.2340  0.0000  0.0071 -0.0186  0.0096 -0.0025  0.0041 -0.0132 -0.0004
## [10] -0.0061  0.0003  0.0058  0.0070 -0.0199  0.0052 -0.0108  0.0103  0.0009
## [19]  0.0004  0.0053  0.0023 -0.0084 -0.0057 -0.0016
## 표준오차의 비교
round(lmp_tidy$std.error, 7) %>% unname()
##  [1] 0.1134970 0.0000018 0.0280851 0.0086673 0.0074233 0.0062419 0.0049800
##  [8] 0.0077550 0.0064699 0.0059627 0.0054670 0.0051890 0.0047140 0.0142877
## [15] 0.0120223 0.0100260 0.0078793 0.0064157 0.0055760 0.0049304 0.0016327
## [22] 0.0101919 0.0020832 0.0094339
round(nonlmp_tidy$std.error, 7)
##  [1] 0.1181076 0.0000017 0.0278586 0.0085666 0.0073208 0.0061931 0.0049657
##  [8] 0.0078289 0.0065164 0.0059892 0.0054953 0.0051931 0.0047191 0.0142304
## [15] 0.0119511 0.0099611 0.0078607 0.0064120 0.0055719 0.0049197 0.0017221
## [22] 0.0134810 0.0078529 0.0136041

계수는 한계효과를 보여주기 때문에 이 경우에는 해석이 꽤 쉬운 편입니다. 이번에는 가격에 대한 선형회귀모형을 분석하고 그 결과를 좀 더 보기 좋은 플롯으로 바꾸어 보겠습니다. 가격을 종속변수로 쓰는 이유는 다이아몬드 데이터는 실제 데이터기 때문에 무작위 값을 가지지 않기 때문입니다 (관계가 나타나면 잘 보일 것).

# 먼저 변수들을 코딩해주도록 하겠습니다. 기존의 다이아몬드 데이터에서 n = 1000인 표본을 추출합니다.
new_df <- df %>% sample_n(1000)
full_model <- lm(
  price ~ carat + cut + color + clarity + depth + x + y + z, data = new_df
)

summary(full_model)
## 
## Call:
## lm(formula = price ~ carat + cut + color + clarity + depth + 
##     x + y + z, data = new_df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7031.3  -499.9  -101.7   362.5  7491.1 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 10094.484   8214.436   1.229  0.21942    
## carat       13273.717    355.421  37.346  < 2e-16 ***
## cut.L         419.188    145.131   2.888  0.00396 ** 
## cut.Q        -122.962    125.758  -0.978  0.32843    
## cut.C         -28.513    102.127  -0.279  0.78015    
## cut^4         -66.788     79.632  -0.839  0.40184    
## color.L     -1904.263    115.041 -16.553  < 2e-16 ***
## color.Q      -765.285    106.421  -7.191 1.28e-12 ***
## color.C      -164.296     96.657  -1.700  0.08949 .  
## color^4      -123.002     86.268  -1.426  0.15424    
## color^5       -49.348     81.882  -0.603  0.54686    
## color^6       -32.431     75.668  -0.429  0.66831    
## clarity.L    3880.457    220.947  17.563  < 2e-16 ***
## clarity.Q   -1853.309    206.181  -8.989  < 2e-16 ***
## clarity.C    1005.449    175.756   5.721 1.41e-08 ***
## clarity^4    -284.212    139.859  -2.032  0.04241 *  
## clarity^5       4.263    111.340   0.038  0.96946    
## clarity^6    -165.885     93.486  -1.774  0.07630 .  
## clarity^7     -20.511     79.257  -0.259  0.79585    
## depth        -104.195    132.518  -0.786  0.43190    
## x           -3751.610    856.755  -4.379 1.32e-05 ***
## y             973.672    695.218   1.401  0.16167    
## z            1428.193   2203.711   0.648  0.51708    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 995.4 on 977 degrees of freedom
## Multiple R-squared:  0.939,  Adjusted R-squared:  0.9376 
## F-statistic: 683.8 on 22 and 977 DF,  p-value: < 2.2e-16

full_model은 가격을 종속변수로 하고, carat, cut, color, clarity, depth, x, y, z를 설명변수로 하는 다중선형회귀모델입니다. 이제 full_model의 결과를 broom 패키지의 함수들을 이용해 깔끔하게 정리해보겠습니다.

new_tidy_model <- broom::tidy(full_model) %>% print(n = Inf) 
## # A tibble: 23 x 5
##    term        estimate std.error statistic   p.value
##    <chr>          <dbl>     <dbl>     <dbl>     <dbl>
##  1 (Intercept) 10094.      8214.     1.23   2.19e-  1
##  2 carat       13274.       355.    37.3    2.31e-190
##  3 cut.L         419.       145.     2.89   3.96e-  3
##  4 cut.Q        -123.       126.    -0.978  3.28e-  1
##  5 cut.C         -28.5      102.    -0.279  7.80e-  1
##  6 cut^4         -66.8       79.6   -0.839  4.02e-  1
##  7 color.L     -1904.       115.   -16.6    1.94e- 54
##  8 color.Q      -765.       106.    -7.19   1.28e- 12
##  9 color.C      -164.        96.7   -1.70   8.95e-  2
## 10 color^4      -123.        86.3   -1.43   1.54e-  1
## 11 color^5       -49.3       81.9   -0.603  5.47e-  1
## 12 color^6       -32.4       75.7   -0.429  6.68e-  1
## 13 clarity.L    3880.       221.    17.6    3.20e- 60
## 14 clarity.Q   -1853.       206.    -8.99   1.27e- 18
## 15 clarity.C    1005.       176.     5.72   1.41e-  8
## 16 clarity^4    -284.       140.    -2.03   4.24e-  2
## 17 clarity^5       4.26     111.     0.0383 9.69e-  1
## 18 clarity^6    -166.        93.5   -1.77   7.63e-  2
## 19 clarity^7     -20.5       79.3   -0.259  7.96e-  1
## 20 depth        -104.       133.    -0.786  4.32e-  1
## 21 x           -3752.       857.    -4.38   1.32e-  5
## 22 y             974.       695.     1.40   1.62e-  1
## 23 z            1428.      2204.     0.648  5.17e-  1
new_tidy_conf <- broom::confint_tidy(full_model)
new_tidy_conf <- new_tidy_conf %>%
  rename(ll = conf.low, ul = conf.high)
new_tidy_model <- bind_cols(new_tidy_model, new_tidy_conf) 
print(new_tidy_model, n = Inf)
## # A tibble: 23 x 7
##    term        estimate std.error statistic   p.value     ll       ul
##    <chr>          <dbl>     <dbl>     <dbl>     <dbl>  <dbl>    <dbl>
##  1 (Intercept) 10094.      8214.     1.23   2.19e-  1 -6025. 26214.  
##  2 carat       13274.       355.    37.3    2.31e-190 12576. 13971.  
##  3 cut.L         419.       145.     2.89   3.96e-  3   134.   704.  
##  4 cut.Q        -123.       126.    -0.978  3.28e-  1  -370.   124.  
##  5 cut.C         -28.5      102.    -0.279  7.80e-  1  -229.   172.  
##  6 cut^4         -66.8       79.6   -0.839  4.02e-  1  -223.    89.5 
##  7 color.L     -1904.       115.   -16.6    1.94e- 54 -2130. -1679.  
##  8 color.Q      -765.       106.    -7.19   1.28e- 12  -974.  -556.  
##  9 color.C      -164.        96.7   -1.70   8.95e-  2  -354.    25.4 
## 10 color^4      -123.        86.3   -1.43   1.54e-  1  -292.    46.3 
## 11 color^5       -49.3       81.9   -0.603  5.47e-  1  -210.   111.  
## 12 color^6       -32.4       75.7   -0.429  6.68e-  1  -181.   116.  
## 13 clarity.L    3880.       221.    17.6    3.20e- 60  3447.  4314.  
## 14 clarity.Q   -1853.       206.    -8.99   1.27e- 18 -2258. -1449.  
## 15 clarity.C    1005.       176.     5.72   1.41e-  8   661.  1350.  
## 16 clarity^4    -284.       140.    -2.03   4.24e-  2  -559.    -9.75
## 17 clarity^5       4.26     111.     0.0383 9.69e-  1  -214.   223.  
## 18 clarity^6    -166.        93.5   -1.77   7.63e-  2  -349.    17.6 
## 19 clarity^7     -20.5       79.3   -0.259  7.96e-  1  -176.   135.  
## 20 depth        -104.       133.    -0.786  4.32e-  1  -364.   156.  
## 21 x           -3752.       857.    -4.38   1.32e-  5 -5433. -2070.  
## 22 y             974.       695.     1.40   1.62e-  1  -391.  2338.  
## 23 z            1428.      2204.     0.648  5.17e-  1 -2896.  5753.

마지막으로 이렇게 ’표’로 보여주는 것이 아니라 가독성이 좋게 계수플롯(coef.plot)의 형태로 나타내 보겠습니다.

library(wesanderson)
new_tidy_model %>%
  dplyr::filter(term != "(Intercept)") %>% # 절편, 상수항을 제외한 나머지만
  mutate(                                  
    term = str_replace(term, 'cut', 'cut -- '),  
    term = str_replace(term, 'color', 'color -- '),
    term = str_replace(term, 'clarity', 'clarity -- ')
  ) %>%
  rename(variable = term) %>% 
  ggplot(aes(x = reorder(variable, estimate), y = estimate, 
             # x축에 변수, y축에 추정치, 나타나는 순서는 변수명, 추청치 순
      color = p.value <= 0.05 # 유의수준이 0.05 이하인 경우와 아닌 경우의 색을 구분
    ) 
  ) +
  geom_point(show.legend = F) +
  geom_linerange(aes(ymin = ll, ymax = ul), show.legend = F) +
  geom_hline(yintercept = 0, color = 'purple', linetype = 'dashed') + 
  coord_flip() + # x축과 y축 그림에서 서로 바꾸기
  labs(x = 'Independent Variables', y = 'Estimate with Confidence Intervals') +
  scale_color_manual(values = wes_palette("Royal1"))

이 결과는 예측변수들이 가격과 어떠한 관계를 가지는지 다변량(multivariate) 관계를 바탕으로 보여주고 있습니다.

문제의 진단과 해결?

일반적으로 이탈치와 오차의 비정규적 분포의 문제는 변수를 제곱 또는 제곱근을 취하는 등의 변환(transform)을 통해 통계적으로 해결하는 방법이 있습니다. 이러한 방법들의 기본 논리는 비선형적 관계를 수리적 조치를 이용해 일견 선형성의 틀에 가두는 것과 같습니다. 예를 들어 \(\frac{1}{2}, 1, 2\)의 관계는 비선형적입니다. 하지만 로그화할 경우, \(\text{log}_2(\frac{1}{2}), \text{log}_2(1), \text{log}_2(2)\)는 \(-1, 0, 1\)로 선형적 틀에 가두어 보여줄 수 있습니다.

일정하지 않은 분산의 경우 일종의 가중치를 가중최소제곱 기법(weighted-least-squares)을 사용해볼 수 있습니다. 가중치를 부여한 제곱합을 최소화하는 방법으로 더 작은 분산을 가진 관측치에 가중치를 주어 일정하지 않은 분산에 대한 표준오차를 보정(correct)하는 것입니다. 또는, Huber-White 로버스트 표준오차라는 방법으로 공분산행렬의 대각행렬에 제곱근을 취하는 방밥이 있습니다. 자세한 내용은 King and Roberts (2015)를 참고하시기 바랍니다.

마지막으로 공선성의 문제를 해결하는 방법으로는 변수에 새로운 정보를 더하여 새로운 측정지표를 만들어내는 것이 있습니다. 그리고 새로운 정보를 포함한 지표와 기존 변수들 간에 어떤 것이 더 \(Y\)를 잘 설명하는지 순위를 매겨보는 것이죠. 공선성의 문제는

  • 변수의 수를 줄이고

  • \(X\)들 간의 독립성을 확보하고

  • 변수 간 상호이동성(interpretability, 서로 거의 비슷하면 사실 이동성이 매우 높은 것이라고 할 수 있습니다)을 줄이는

것을 통해 해결할 수 있습니다. 따라서 구성개념의 다른 측면을 보여줄 수 있는 경험지표를 설계하는 모델의 재구성(model respecification)과 변수 선택, 그리고 일종의 편향된 추정방법(e.g., ridge regression) 등을 사용하여 해결하고는 합니다.

  1. 불편추정량을 의미한다. Best Unbiased Linear Estimator.↩︎

  2. 다만 이것은 통계적 진단일 뿐입니다. 극단적인 값을 지니는 이른바 이탈치(outliers)에 대한 처치를 어떻게 할 것인가는 연구자의 관점에 좌우됩니다. 자세한 논의는 김웅진. 2016. “Deviant Case in the Social Scientific Generalization: A Methodological Discourse on the Analytic Conventions for Detecting Deviance.” 21세기정치학회보. 26권 4호, pp. 49-66. 또는 김웅진. 2016. “사회과학적 개념의 방법론적 경직성 : 국소성과 맥락성의 의도적 훼손.” 국제지역연구. 19권 4호, pp. 3-22.를 읽어보시기 바랍니다.↩︎