11. 기초통계이론 2단계

Chapter11. 기초통계이론 2단계

1. 통계 모형 Preview

본격적으로 모델링을 시작하기에 앞서 간단한 주의사항 및 분석모형에 대한 소개를 하고 넘어가도록 하겠습니다. 흔히 사람들이 분석과정에서 실수하는 경우는 다음과 같은 상황입니다. 힘들게 어려워 보였던 예측 모형 및 알고리즘을 공부하고 이를 사용하기 위하여 바로 모델에 적합시키려고 하는 경우입니다. 물론 복잡한 모형을 적용시키면 결과가 잘 나올거 같고 멋도 있어 보일 수 있지만, 전혀 그렇지 않습니다. 데이터 분석은 요리와 비슷합니다. 어떤 요리를 할지 결정하고, 가져온 재료를 보지는 않죠. 데이터도 마찬가지입니다. 결국 분석 모형은 데이터에 맞는 분석모형을 적용시켜야지, 무작정 어려운 분석모형 적용시킨다고해서 결과가 다 나오는 것이 아닙니다.

그러므로 항상 데이터에 대한 특성을 충분히 이해하고 있어야 하며, 이런 상황에는 어떤 분석 모형을 적용시켜야 하는지 바로 생각이 떠올라야 할 것입니다. 해당 부분의 이해를 위하여 간단하게 통계 모형들에 대하여 다루고 넘어가도록 하겠습니다.

2. 분석모형을 선택하는 기준

1. 분석하고자 하는 바를 명확히 규명해야 합니다.

이는 가설을 세우는 과정과 같습니다.

2. 가설에 맞는 데이터들에 대한 변수 척도 구분해야 합니다.

변수의 척도에 따라 적용해야 되는 모형이 정해집니다. 모형은 우리가 정하는 것이 아닌, 데이터가 정해줍니다.

3. 분석하고자 하는 주제가 ’차이’를 검정하고 싶은 것인지, ’관계’를 검정하고 싶은 것인지에 따라 갈립니다.

차이를 보는 검정은 흔히 집단 간에 평균 차이를 검정하고 싶은 경우 입니다. 예를 들어 인사관리 데이터에서 이직 여부(0: 안함, 1: 이직)에 따라 직무 만족도가 다른지 검정하는 경우가 있을 수 있습니다. 반대로 관계를 보고자 하는 검정은변수 간의 함수적 관계를 보고자하는 것과 같습니다. 예를 들어 마케팅 투자비용이 마케팅 효과에 미치는 영향을 검정하고 싶은 경우가 있을 수 있습니다.

4. 다음으로는 Response Variable(종속 변수)가 연속형이냐, 이산형이냐에 따라 결정됩니다.

• 종속 변수가 연속형일 때
  • 차이를 보고자 할 때 : T 검정(T-test), 분산분석(Anova)
  • 관계를 보고자 할 때 : 회귀분석(Regression)
• 종속 변수가 이산형일 때
  • 연관을 보고자 할 때 : 카이제곱 독립성 검정(Chi square Independent Test)
  • 관계를 보고자 할 때 : 로지스틱 회귀분석(Logisitc Regression)

3. 분석모형별 가설 검정

통계 분석에서 차이가 없는(즉 의미가 없는) 사실은 알아내고 싶은 사실이 아닙니다. 차이가 있는(즉 의미가 있는, 유의하다) 사실에 관심이 있으며, 의미가 있는 사실을 대립가설에 배치합니다.

4. \(t\)검정

일반적으로 \(t\)검정이라하면 독립표본 \(t\)검정을 의미합니다. 독립적인 두 집단에서 추출된 표본들의 평균이 차이가 있는지 확인하기 위해 실시되며 혹시 두 집단이 상황적으로 독립이 아닌 집단이라면 대응표본 \(t\)검정(paired t-test)등을 포함한 다른 분석을 진행하게 됩니다.

원리는 위에서 가설검정 파트에서 했던 일표본 \(t\)검정과 동일합니다. 두 개의 독립적인 정규모집단에서 추출된 표본을 통해 계산된 두 집단의 표본평균 차를 이용합니다. 두 표본평균 차의 분포를 확인 후 우리가 가지고 있는 표본에서 계산된 차이가 두 평균이 같다고 가정했을 때(귀무가설이 사실일 때) 충분히 나올 법한 차이인가를 검정하는 과정입니다.

독립표본 t검정의 특성

• 분석의 목적은 두 집단의 평균이 차이가 있는지 확인하는 것이다.
• 분석의 대상은 독립적인 두 정규모집단에서 추출된 표본들이다.
• 분석의 수단은 \(t\)분포다.
• 분석의 원리는 추출된 표본으로 계산된 차이가 우연히 나올만한 차이인지 아닌지를 \(t\)분포(두 표본평균 차의 분포)를 이용하여 확인하는 것이다.

하지만 여기서 하나의 문제가 생깁니다. 바로 분산에 대한 처리입니다. 가설검정 파트에서 했던 일표본 \(t\)검정의 경우 집단이 하나밖에 없으니 상관 없지만 집단이 두 개인 이 경우에서는 분산을 어떻게 처리해야 할까요.

독립표본 \(t\)검정에서는 두 집단의 분산이 같은 경우와 다른 경우를 나누어 생각합니다. 이를 흔히 등분산 가정이라 합니다. 일반적으로 \(t\)검정에서는 두 집단의 등분산이 같은 경우, 즉, 등분산 가성이 성립되었을 때는 정확하게(exact) \(t\)분포를 이용하고 그렇지 않을 경우에는 근사적인(approximate) 검정을 실시합니다. 그런 관계로, \(t\)검정 전에는 반드시 두 집단의 분산이 같은지에 대한 검정을 진행해야 합니다. 참고로 등분산을 확인하는 방법은 F분포를 이용한 Levene 등분산 검정(\(H_0: 두 \; 집단의\;분산은\; 동일하다\))을 비롯한 몇 가지의 검정이 있으나 여기서 자세히 설명하진 않겠습니다.

만약 두 집단의 등분산 가정이 성립한다고 판단이 되면 집단에 상관없이 전체의 분산을 계산해서 사용합니다. 이를 합동 표본분산(pooled sample variance)이라고 합니다. 그 형태는 다음과 같습니다. 여기서 \(n_1\)은 X표본(집단 1)에서의 표본 수고 \(n_2\) 는 Y표본(집단 2)에서의 표본수입니다.
\[ s_p ^2 = \frac 1 {n_1 + n_2 -2 } \left\{ \sum_{i=1}^{n_1 }(X_i -\overline X )^2 + \sum_{j=1}^{n_2 }(Y_j -\overline Y )^2 \right\} \]

합동 표본분산의 형태를 살펴보면 두 집단의에서의 변동을 각각 계산하고 더한 후 전체 자유도로 나누어준 것을 확인할 수 있습니다. 자유도는 역시 기존 다른 추정량들과 마찬가지로 합동 표본분산을 추정하는데 온전히 사용된 데이터 수로 결정됩니다. 이와 같은 것을 pooling 한다고 표현하며 단일 집단에서의 표본평균과 비교해보면 pooling이 어떤 구조로 이루어지는 것인지 어렵지 않게 이해하실 수 있을 것입니다.
\[ (단일\;집단에서의 \; 표본평균) \qquad s^2 = \frac 1 {n-1} \sum_{i=1}^{n} (X_i -\overline X)^2 \]

이제 이를 이용해 두 집단의 표본평균 차에 대한 분포를 도출해봅시다.

두 독립적인 정규모집단에서 추출된 표본으로 계산된 표본평균들은 정규분포를 따르기 때문에 정규분포 파트에서 학습한대로 그들의 선형결합 \((\overline X - \overline Y)\) 역시 정규분포를 따를 것입니다. 이것에 착안하여 일표본 \(t\)검정에서와 같이 스튜던트 정리를 이용합니다. \[ \frac {(\overline X-\overline Y) - (\mu_1 - \mu_2) } { \sqrt{\frac {s_p^2} n_1 + \frac {s_p^2} n_2}} \; \sim \; t(n_1+n_2 -2) \]

여기서 \(\mu_1\) 과 \(\mu_2\)는 각 두 집단으 모평균이고 그 차이인 \((\mu_1 -\mu_2)\)이 우리의 최종 목적입니다. 만약 두 집단의 사이의 유의한 차이가 없다면 0이 될 것이고 그렇지 않다면 충분히 유의한 차이를 보일 것입니다.

이 역시 일표본 \(t\)검정에서의 검정통계량과 비교해보시면 같은 형태를 띠는 것을 알 수 있습니다. \[ (단일 \; 집단에서의 \; 검정통계량) \qquad \frac { \overline X -\mu} {s/\sqrt n } \;=\; \frac {\overline X -\mu}{\sqrt {\frac {s^2} n}} \; \sim \;t(n-1) \]

이제 이 검정통계량을 이용하여 우리가 원하는 \(t\)검정을 실시할 수 있습니다.

위에서 말씀드린대로 우리의 최종 목적은 두 집단의 평균이 같은지 다른지를 확인하는 것입니다. 이를 귀무가설과 대립가설로 표현하면 다음과 같습니다. \[ H_0 \; : \; \mu_1 = \mu_2 \\ H_1 \; : \; \mu_1 \neq \mu_2 \]

만약 귀무가설이 사실이라고 가정하면 위 검정통계량에서 \((\mu_1 -\mu_2)\)은 0이 되므로 검정통계량은 좀 더 간소화가 가능합니다. 이는 귀무가설 하에서의 검정통계량이니 \(T_0\) 라고 하겠습니다. \[ T_0\;=\; \frac {(\overline X-\overline Y) } {\sqrt{\frac {s_p^2} n_1 + \frac {s_p^2} n_2}} \; \sim \; t(n_1+n_2 -2) \]

판단은 원리는이 \(T_0\) 를 이용하여, 우리가 가지고 있는 표본들로 계산된 검정통계량이 과연 (주어진 \(t\)분포하에서) 자연스럽게 나올만한 값인지, 아니면 우연이라고 보기에는 너무 극단에 가까운 값인지 확률로써 확인하는 것입니다. 그 방식은 일표본 \(t\)검정과 정확히 일치합니다.

만약 계산된 검정통계량보다 극단적인 값이 나올 확률(유의확률)이 사전에 우리가 정한 유의수준 \(\alpha\) 보다 작다면 귀무가설을 기각하고 그렇지 않다면 기각하지 못하게 될 것입니다.

다음은 두 집단의 분산이 다를 경우입니다. 이 때는 일반적인 등분산 \(t\)검정에서의 검정통계량과 형태는 같게 하되 두 집단의 분산을 따로 추정합니다. 여기서 \(s_1^{2}\)과 \(s_2^{2}\) 는 각각 집단1과 집단2의 표본분산입니다. \[ \frac {(\overline X-\overline Y) - (\mu_1 - \mu_2) } { \sqrt{\frac {s_1^2} n_1 + \frac {s_2^2} n_2}} \; \sim \; t(degree\; of \; freedom) \]

그런데 이 경우 \(t\)분포의 자유도를 계산하기가 애매합니다. 그래서 대부분의 경우는 세터스웨이트 공식과 같은 자유도의 근사치를 얻는 방법을 이용하게 됩니다. R을 비롯한 통계 프로그램에서는 자동으로 자유도 근사치를 얻어 검정을 진행합니다.

5. \(t\)검정(R Code)

데이터는 앞단에서 다루었던 HR데이터를 이용하겠습니다.

이직 여부에 따라 직원들의 직무만족도에 차이가 있는지 검정을 하고자 합니다. 이 경우, 이직 여부(left)는 0 : 이직 안함, 1 : 이직이기 때문에 이직 여부는 2개의 수준을 가지고 있는 명목형 변수이고직무만족도(satisfaction_level)는 0 ~ 1 사이에 있는 연속형 변수입니다.

즉, 연속형 변수를 두 수준을 지니고 있는 명목형 변수에 따라 차이가 있는지 검정하고 싶기에 T 검정을 진행하는 것이 적합한 상황입니다.

T 검정을 R에서 진행하는 방법은 다음과 같습니다.

• 등분산 검정

비교하고자 하는 두 잡단의 분산이 같은지 검정하기 위함입니다.

\[ H_0 : 두\;집단의\;분산이 \;동일하다.\\ H_1 :두\;집단의\;분산이\;다르다.\\ \]

여기서 두 집단은 left변수의 수준인 0(이직 안함)과 1(이직)을 의미합니다.

# 라이브러리 불러오기
library(car)
# 등분산검정 실행
HR = read.csv('F:/Dropbox/DATA SET/HR_comma_sep.csv')
HR$left = as.factor(HR$left)
leveneTest(satisfaction_level ~ left , data = HR)

Levene's Test for Homogeneity of Variance (center = median)
         Df F value    Pr(>F)    
group     1   122.4 < 2.2e-16 ***
      14997                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Levene’s 등분산 검정은 두 집단의 분산비를 비교하기 때문에 F분포를 따르게 됩니다.

R 실행 결과값에 나오는 \(F\ value = 122.4\)는 위 분포에서의 검정통계량입니다. 앞에 가설검정단계에서 검정통계량을 기준으로 계산 된 유의확률을 통해 귀무가설 기각 여부를 판단할 수 있다고 했습니다. 유의확률을 나타내는 \(Pr(>F)\)은 \(< 2.2e-16\)으로 표시가 되어 있습니다. R에서 \(e-16\)은 \(\frac{1}{10^{16}}\)을 의미합니다. 그러니 유의확률은 0에 매우 근접한 값을 가지는 것을 뜻하며 유의수준 \(\alpha=0.05\)보다 훨씬 작기에 두 집단의 분산이 동일하다.라는 귀무가설을 기각할 수 있습니다. 즉, 이직을 한 직원들과 이직을 하지 않은 직원들 간의 직무만족도의 분산은 동일하지 않다라는 것을 알 수 있습니다.

• \(t\)검정

앞서 등분산 검정의 결과가 어떻게 나왔는가에 따라서 t 검정의 옵션이 변합니다. 검정의 귀무가설과 대립가설은 다음과 같습니다. \[ H_0:두\;집단의\;평균이\;같다.\\ H_1:두\;집단의\;평균이\;다르다. \]

# 등분산이 동일할 경우
t.test(satisfaction_level ~ left , data = HR,var.equal = TRUE)

    Two Sample t-test

data:  satisfaction_level by left
t = 51.613, df = 14997, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 0.2181017 0.2353215
sample estimates:
mean in group 0 mean in group 1 
      0.6668096       0.4400980 

# 등분산이 동일하지 않을 경우
t.test(satisfaction_level ~ left , data = HR,var.equal = FALSE)

    Welch Two Sample t-test

data:  satisfaction_level by left
t = 46.636, df = 5167, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 0.2171815 0.2362417
sample estimates:
mean in group 0 mean in group 1 
      0.6668096       0.4400980 

현재 데이터는 두 집단의 분산이 동일하지 않으므로 t.test에서 옵션을 var.equal = FALSE로 설정해줍니다.

분석 결과 검정통계량 \(t\ value = 46.636\)이며 유의확률(\(p-value\))은 \(<2.2e-16\)으로 0과 매우 가깝습니다. 이로써 두 집단(이직 여부)의 평균이 같다라는 귀무가설을 기각할 수 있으며, 이직 여부에 따라 직무만족도의 평균차이가 존재한다는 것을 알 수 있습니다. R 결과표에 따르면 이직을 하지 않은 집단의 평균(mean in group 0)은 0.66이고 이직을 한 집단의 평균(mean in group 1)은 0.44로 이직을 하지 않은 집단의 직무만족도가 더 높다고 할 수 있습니다.

6. 분산분석

분산분석은 \(t\)검정과 마찬가지로 차이를 보는 분석입니다. \(t\)검정과 다른 점은 분산분석은 두 집단은 물론이고 세 집단 이상에서도 집단 간 평균 차이를 볼 수 있는 점입니다. 두 집단에서 진행하면 \(t\)검정과 같은 결과를 줍니다.

ANOVA(Analysis of Variance)의 의미를 조금 풀어서 해석해보자면 다음과 같습니다.

분산(변동)을 분석하여 평균을 비교한다.정도로 해석할 수 있습니다. 보통 차이를 비교할 때, 3개 이상의 개체에 대한 동시비교는 불가능합니다. 예를 들어, 사고 싶은 제품이 3개가 있는데, 그 중 1개만 살 수 있는 상황이라면 본능적으로 3개 중에 2개를 먼저 비교하고, 그 중 선택된 1개가 나머지 1개를 비교하여 구매할 제품을 선택하는 과정과 같습니다. 분산분석은 이러한 한계점을 극복하고자 직접적으로 평균을 비교하는 것이 아닌, 분산을 이용하여 평균을 비교하는 방법입니다. 그래서 명칭도 ’분산을 분석한다’의 의미로 Analysis of Variance라고 부릅니다.

예를 들어 HR 데이터 셋에서 satisfaction_level(직무 만족도) 변수를 salary(연봉 수준, Low, Mid, High 3개의 수준을 가지고 있는 순서형 변수) 변수에 따라 평균차이가 있는지 없는지 통계적으로 검정하고자 합니다. 그렇다면 연속형 변수를 3 집단 간의 평균차이를 하고자 하는 것이니, \(t\)검정이 아닌 분산분석을 적용해야합니다.

분산분석의 아이디어는 다음과 같습니다.

1. 먼저 전체집단의 평균(\(\overline{X}{\cdot\cdot}\)) , 그리고 각 집단의 평균(\(\overline{X}{i\cdot}\))을 구합니다. 여기서 \(i\)는 각 집단을 의미합니다.(salary변수의 low, mid, high)

2. 임의의 데이터(\(X{ij}\))를 기준으로 총 변동을 계산합니다.

총 변동 : \(\overline{X}{\cdot\cdot}-X{ij}\)

3. 계산 된 총 변동을 집단 간(between) 변동과 집단 내(within) 변동으로 분리를 합니다.

집단 간(between) 변동 : \(\overline{X}{\cdot\cdot} - \overline{X}{i\cdot}\)
집단 내(within) 변동 : \(\overline{X}{i\cdot}-X{ij}\)

총 변동을 집단 간 및 집단 내로 분리한다는 것은 총 변동 중에 집단 간(between) 변동이 집단 내(within) 변동에 비해 얼마나 많은 부분을 차지하고 있는가를 보고자 하는 것입니다.

4. 모든 임의의 데이터(\(X_{ij}\))에 대해서 변동의 제곱합을 구해 줍니다.

제곱합을 구해주는 이유는 임의의 데이터(\(X_{ij}\)) 위치에 따라서 총 변동의 값이 양수(+)가 계산 될 수가 있고 혹은 음수(-)가 계산 될 수가 있기 때문입니다. 부호간 덧셈을 하면 변동의 합이 상쇄가 되기 때문에 제곱합을 함으로써 문제를 해결해주고자 하는겁니다.

5. 계산된 제곱합을 평균제곱합으로 보정을 해준 다음, F 통계량을 구해 유의확률(p-value)을 계산합니다.

집단 간 변동이 집단 내 변동에 비해 유의한 크기를 가지고 있다면, 귀무가설을 기각할 수 있습니다.

이제부터, 분산분석의 원리와 관련된 몇 가지 키 아이디어를 말씀드리겠습니다.

• 전체 평균을 구하는 이유는 3집단 이상의 변동을 비교하기 위한 기준선을 구하는 것입니다.

• 분산분석의 핵심은 총 변동 중에서 집단 간(between) 변동이 집단 내(within) 변동에 비해 얼마나 큰지 비교하고자 하는 것입니다.

• 집단 간 변동에 집중하는 이유는 분산분석에서 분석하고자 하는 것이 전체 집단의 평균과 특정 집단의 평균의 차이 즉, 다른 집단 간의 평균 차이이기 때문입니다.

• 반대로 집단 내(within) 변동은 동일 집단 내 차이를 나타내기 때문에 집단 간 평균차이를 비교하고자 하는 분산분석에서는 Error취급을 합니다. 만약 집단 내 차이가 크면 총 변동 내에서 집단 간 차이가 차지하는 부분이 줄어들기 때문에, 총 변동이 고정된 상태에서 Error가 작다는 것은 곧 집단 간의 차이가 크다는 것으로 의미가 연결됩니다.

그럼 이제 위에서 언급한 세 가지 변동을 관찰된 데이터를 통해 추정해보도록 하겠습니다. 여기서 변동을 표현할 때 사용되는 SS는 제곱합이란 의미로 Sum of Squares를 의미합니다.

\[ 총\; 변동\;: \;SST = \sum_{i=1}^k \sum_{j=1}^{n_i} (X_{ij}-\overline X_{\cdot \cdot})^2 \qquad\\ i=집단의\; 첨자\;j= 표본의\; 첨자 \] \(SST\)는 보시는바와 같이 어떤 그룹이든 신경쓰지 않고 전부 편차의 제곱합을 시켜 더합니다. \(\overline X_{\cdot \cdot}\) 역시 그룹과 상관없이 모든 표본들의 전체 평균입니다. 만약 위의 예처럼 집단 변수(salary 변수)가 low, mid, high 세 그룹이라면 \(k\) 는 3이 되겠고 \(n_i\) 들은 low, mid, high각 집단에서의 표본 수입니다. 각 집단에서의 표본 수는 달라도 상관은 없으나 되도록 비슷한 것이 분석의 신뢰를 높여줍니다.

\[ 집단 \;간(between) 변동\;:SSG = \sum_{i=1}^k \sum_{j=1}^{n_i} (\overline X_{i\cdot} - \overline X_{\cdot \cdot })^2\; =\;\sum_{i=1}^k n_i(\overline X_{i\cdot} - \overline X_{\cdot \cdot })^2 \qquad \]

집단 간 제곱합은 처리제곱합이라고도 많이 불립니다. 여기서 \(X_{i\cdot}\) 은 \(i\) 그룹의 평균입니다. 예를 들어 \(\overline X _{2\cdot}\)은 salary변수에서 mid에 속하는 표본만의 평균이라고 보시면 됩니다. 즉, 전체 평균과 그룹평균의 차이에 대한 제곱합이죠. 이것이 그룹의 효과입니다. \(SST\)가 그룹과 상관 없는 전체 변동이라면 이는 그룹들이 가지고 있는 변동인 것이죠.

\[ 집단\;내(within) 변동\;:SSE = \sum_{i=1}^k \sum_{j=1}^{n_i}(X_{ij}-\overline X_{i\cdot})^2 \qquad \] 마지막으로 \(SSE\)는 그룹 내 변동을 표현하는 제곱합입니다. 형태를 보시면 개별 관찰값과 그 관측값이 속한 그룹의 그룹 평균의 차이를 이용해 제곱합을 계산하는 구조입니다. 앞에서 집단 내 변동은 보고자 하는 집단 간의 차이와는 반대되는 성향을 가지고 있는 Error라고 하였습니다. 만약 이 \(SSE\)가 \(SSG\)에 비해 상대적으로 작다면 우리는 그룹이 가지고 있는 효과가 유의하다고 판단하 것입니다.

\[ \sum_{i=1}^k \sum_{j=1}^{n_i} (X_{ij}-\overline X_{\cdot \cdot})^2\;=\;\sum_{i=1}^k n_i(\overline X_{i\cdot} - \overline X_{\cdot \cdot })^2\;+\;\sum_{i=1}^k \sum_{j=1}^{n_i}(X_{ij}-\overline X_{i\cdot})^2\\ SST\qquad \quad \quad =\quad \quad \qquad SSG \quad \quad \qquad +\quad \quad \qquad SSE \]

만약 기본적인 가정들이 성립한다면 수리적으로 이 두 제곱합(\(SSG\), \(SSE\))은 직교분해(orthogonal decomposition)되었다고 할 수 있습니다. 직교분해는 교차항이 없이 각 변동을 표현하는 제곱합항만으로 분해되는 것이며 그 의미는 두 제곱합이 통계적으로 독립이라는 것입니다. 그리고 이 변동의 추정치인 제곱합 역시 표본으로부터 만들어진 통계량이므로 확률분포를 갖습니다.

그럼 우리는 서로 독립적인 집단 간 제곱합과 집단 내 제곱합을 갖게 됩니다. 이 둘을 조합하여 상대적 비교가 가능한 F분포를 따르는 새로운 통계량을 얻어낼 수 있습니다. (주어진 표본들을 이용해 이 통계량을 계산한 값을 흔히 F값이라고 부릅니다)

분산분석에서의 가설검정 절차

\[ \frac {SSG/(SSG의 \; 자유도)} {SSE/(SSE의 \; 자유도) } \; \sim \; F(SSG의 \; 자유도, SSE의 \;자유도) \] \(SSG\)와 \(SSE\)의 자유도 역시 예전에 했던 자유도의 논리와 크게 다르지 않습니다.

\(SSG\)는 K개의 그룹의 평균(\(\overline X_{i\cdot}\))이 제곱합을 만드는데 선택될 수 있고, 만드는 과정에서 전체 평균( \(\overline X _{\cdot \cdot}\)) 하나를 추정해서 사용하였으니 1개의 자유도를 잃습니다. 그래서 \(SSG\)의 자유도는 \(K-1\)이 됩니다.

\(SSE\)의 경우, 전체 표본 수가 (여기서는 N이라고 하겠습니다) 제곱합을 만드는데 선택될 수 있고 만드는 과정에서 \(k\)개의 그룹평균(\(\overline X_{i\cdot})\)을 추정해서 사용하였으니 k개의 자유도를 잃습니다. 즉, \(SSE\)의 자유도는 \(N-k\)입니다.

분산분석의 검정통계량은 다음과 같습니다. 이 검정통계량을 가지고 유의확률(p-value)을 계산합니다. \[ F = \frac {SSG / df_{G} } { SSE /df_{E}} \; \sim \; F(df_{G}, df_{E})\\ k= number\; of\; group, \qquad n_i= number\; of\; samples\; for\; each \; group,\qquad N = \sum_{i=1}^k n_i\\ df_{G}=k-1\qquad df_{E}=N-k \]

이렇게 진행된 분산분석의 결과를 정리한 것을 분산분석표(Anova Table)라고 합니다. 분산분석의 결과해석은 항상 분산분석표를 이용하여 해석하니 꼭 기억해두시길 바랍니다

**
분산분석의 원리를 이용하여 분석 모형을 설정해보면 다음과 같습니다.**

\[ Y{ij} =\mu{i} +\epsilon{ij} \; , \; \epsilon \; \sim \; iid \; N (0,\sigma^2)\\ \Leftrightarrow \quad Y{ij} = \mu + \alpha i +\epsilon{ij} \;, \; \epsilon{ij} \sim \; iid\; N(0,\sigma^2)\\ \\\;\\ i = 1,2, \cdots ,k \qquad (각\;집단 )\\ j=1,2,\cdots ,ni \quad (각\;집단에서의\;개별\;관찰치)\\ ai =\mu_{i} -\mu \quad (각\;집단의\;효과) \]

여기서 각 집단의 효과 \(\alpha_{i}\) 는 각 집단의 평균과 전체평균의 차이입니다. 우리가 위에서 변동을 이야기하면서 사용했던 \(\overline X{i\cdot} - \overline X {\cdot \cdot}\) 이 추정하려는 것이 바로 이 각 집단의 효과 \(\alpha_{i}\) 라고 할 수 있습니다. 또한 \(iid\)란 independent and identically distributed의 줄임말로, 어떤 확률변수가 독립적이며 동일한 분포를 보인다는 의미입니다. 위 모형에서는 오차항 \(\epsilon\)이 독립적으로 평균이 0이고 분산이 \(\sigma^2\) 인 정규분포를 따른다는 것으로 해석할 수 있습니다.

귀무가설이 기각이 되었을 때, 즉 평균이 같지 않을 경우

분산분석에서 대립가설의 뜻은 ‘모든 집단별 평균은 같지 않다’ 가 아닌 ‘적어도 하나의 집단 평균은 다르다’ 이라는 점입니다. \(K-1\)개의 집단별 평균이 같다고해도 하나만 유의하게 다르다면 귀무가설은 기각될 것입니다. 이는 분산분석이 가지고있는 구조적 한계입니다. 지금까지 언급했던 검정 원리를 보아도 아노바는 집단 간 변동의 총합을 이용한 검정을 하고, 그 말은 특정 집단의 개별 효과를 전혀 확인할 수 없다는 뜻입니다. 다시 말해서 귀무가설이 기각되었을 때, 검정의 결과를 보면 각 집단별 평균이 다르다는 것은 알 수 있지만 어떤 집단이 얼만큼 큰지는 알 수 없다는 것이죠. 그래서 이 경우 대게 사후검정(post hoc test)라 부르는 집단별 비교 검정을 진행하게 됩니다. 원리는 아주 간단합니다.

위의 예에서 low, mid, high 세 집단의 효과가 같은지에 대한 분산분석 결과 귀무가설이 기각되어 세 집단의 효과가 모두 같지는 않다는 결과를 얻었다면 사후검정으로 (low-medium), (low-high), (medium-high) 이렇게 두개씩 비교하면서 개별 비교를 하는 것입니다. 그 분석 방법으로는 위에서 배운 독립 \(t\)검정을 이용할 수도 있겠고 같은 원리를 공유하는 Duncan이나 Scheffe라고 불리는 개별 비교법을 이용할 수도 있습니다.

7. 분산분석(R Code)

분산분석을 R에서 실행하는 방법은 다음과 같습니다. 분석하고자 하는 것은 satisfaction_level(직무 만족도)의 평균이 salary(연봉 수준, low, medium, high) 집단에 따라 차이가 있는지 통계적으로 검정하고자 합니다.

귀무가설은 다음과 같습니다. \[ H_0:salary(연봉) \; 수준별로\;satisfaction\; level(직무만족도)의 \;평균이\;같을\;것이다.\\ H_1:not\;H_0 \]

ANOVA = aov(satisfaction_level ~ salary, data = HR)
summary(ANOVA)

               Df Sum Sq Mean Sq F value   Pr(>F)    
salary          2    2.3  1.1693   18.96 5.97e-09 ***
Residuals   14996  924.8  0.0617                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

분석결과를 보시면 분산분석표가 나옵니다. Df는 자유도, Sum sq는 제곱합을 Mean Sq는 평균제곱합을 의미합니다. 위에서 언급했던 분산분석표를 그대로 가져와서 값을 채워넣으면 다음과 같습니다.

정해진 순서대로 제곱합과 자유도를 구한 후 평균제곱을 계산한뒤 검정통계량인 F value를 구했습니다. 마지막으로 유의확률 p-value를 확인할 수 있습니다. p-value는 0과 매우 가까운 값이기 때문에 당연히 유의 수준 0.05보다 작습니다. 그러므로 귀무가설 \(H_0\)를 기각할 수 있습니다. 즉 연봉 수준별로 직무만족도의 평균은 동일하지 않습니다.

하지만 이 상태로는 위에서 말했듯이 ’평균이 같지 않다.’만을 알 수 있지, 정확히 어떻게 다른지 알 수가 없기 때문에 사후검정을 진행해야됩니다.

TUKEY = TukeyHSD(ANOVA)
TUKEY

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = satisfaction_level ~ salary, data = HR)

$salary
                   diff         lwr          upr     p adj
low-high    -0.03671654 -0.05461098 -0.018822102 0.0000046
medium-high -0.01565305 -0.03372130  0.002415191 0.1049520
medium-low   0.02106349  0.01111999  0.031006988 0.0000021

plot(TUKEY)

사후검정은 여러 방법이 있습니다. 지금 사용한 방법은 Tukey 사후검정입니다. 해석을 하는 방법은 다음과 같습니다. 3개의 집단을 2개씩 따로 분석을 한 것이며, 비교하는 두 집단 간의 차이가 같은지 다른지 다시 한번 검정하는 것입니다. plot은 사후검정 결과를 나타낸 것이며 보는 법은 간단합니다. 각 신뢰구간 안에 0이 포함되는지 아닌지만 보시면 됩니다.

• low - high : low 집단의 satisfaction_level의 평균과 high 집단의 satisfaction_level의 평균차이를 비교하니 신뢰구간이 \([-0.054(lwr),-0.018(upr)]\)사이에 있는 것을 알 수 있습니다. 해당 구간 안에 0이 포함되어있지 않으며, 0보다 작기 때문에 low집단의 satisfaction_level 평균은 high 집단의 satisfaction_level에 비해 작다는 것을 의미합니다.

• medium - high : medium 집단과 high 집단의 satisfaction_level의 비교를 진행한 결과 신뢰구간이 \([-0.033,0.0024]\)입니다. 0을 포함하기 때문에 medium집단과 high집단의 satisfaction_level의 평균은 같다고 볼 수 있습니다.

• medium - low : 신뢰구간이 0을 포함하지 않고 0보다 크기 때문에 medium이 low집단보다 satisfaction_level의 평균이 높다고 볼 수 있습니다.

정리를 하면 satisfaction_level의 평균은 \(low < medium = high\) 순서인 것을 알 수 있습니다. 즉 low집단의 평균이 낮기 때문에 귀무가설을 기각하게 된 것입니다.

지금까지 진행한 이론과 분석을 정확히는 일원배치 분산분석(One way Anova)이라고 합니다. 만약 평균을 비교할 때 하나의 요인만 사용하지 않고 여러개의 요인을 사용하여 분석을 하고 싶을 때가 있습니다. 만약 2개의 요인을 통해 평균을 비교하고자 하면 이원배치 분산분석(Two way Anova)이라고 합니다.

T_ANOVA  = aov(satisfaction_level ~ salary + left +
               salary:left, data = HR)
summary(T_ANOVA)

               Df Sum Sq Mean Sq  F value   Pr(>F)    
salary          2    2.3    1.17   22.276 2.19e-10 ***
left            1  137.8  137.79 2624.960  < 2e-16 ***
salary:left     2    0.0    0.01    0.155    0.856    
Residuals   14993  787.0    0.05                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

요인을 하나 더 넣어주면 바로 이원배치 분산분석을 진행할 수 있습니다. 다만 신경써야할 부분은 요인을 두개를 사용할 경우, 각 요인 간의 교호작용이 있는지 확인을 해야됩니다. 위 분석 결과에서는 단일 효과(salary, left)는 유의한 것을 알 수 있지만 교호작용효과(salary:left)는 유의하지 않다는 것을 알 수 있습니다. 이런 경우 교호작용항을 제거하고 모형을 다시 분석하면 됩니다.

8. 상관분석

만약 두개의 변수 간 관계를 분석하고 싶은 경우, 우리는 종종 상관계수(correlation)이란 것을 구하고는 합니다. 너무 유명한 용어라서, 상관계수의 정확한 의미는 알지 못하더라도, 상관계수가 대충 어떤 것인지는 많은 사람들이 알고 있습니다. 한번 개념을 정립하고 넘어가도록 하겠습니다. 먼저, 상관분석이란 두 변수의 관계에서 하나의 변수가 증가하면, 다른 하나의 변수도 증가하는지 혹은 감소하는 경향이 있는지 확인을 하기 위해 분석을 진행합니다. 우리는 그러한 경향을 확인하기 위해 공분산(Covariance)이라는 값을 계산합니다.

공분산과 상관계수

\[ COV[X,Y] = E[(X- \overline X)(Y- \overline Y)] \] 이렇게 계산을 하면, X와 Y의 상관관계를 계산할 수가 있습니다. 하지만, 한가지 문제가 존재합니다. 공분산은 변수의 단위에 따라 범위가 무한대까지 확장됩니다. 즉, 변수가 바뀌면 공분산의 단위도 바뀌기 때문에 비교하는 값으로 확인하기에는 문제가 존재합니다. 이러한 문제점을 해결하기 위하여 공분산에 두 변수의 분산을 나누어줍니다. 그러면 공분산은 -1 ~ 1의 변수의 단위에 상관 없이 일정한 범위를 가지는 상관계수(Correlation)로 변환이 됩니다.

\[ Corr[X,Y] = \frac{COV[X,Y]}{VAR[X]\ VAR[Y]} \\ -1 \leq Corr[X,Y] \leq 1 \]

해석방법은 매우 간단합니다. 상관계수가 1에 가까울 수록 강한 긍정관계를 가지고 있는 것이고 -1에 가까울 수록 강한 부정관계를 가지고 있습니다. 만약 0에 가까울 경우, 두 변수는 관계가 없다고 해석할 수가 있습니다.

상관분석

상관분석이란, 두 변수 간의 상관관계가 0인지 아닌지 확인을 하는 통계적 검정방법입니다. 따라서 귀무가설과 대립가설은 다음처럼 세울 수가 있습니다.

\[ H_0: \rho=0\\ H_1: \rho \neq 0 \]

9. 상관분석 (R code)

상관분석에서 주의할 점은 상관분석은 단순히 두 변수 간의 선형 관계를 파악하는 것뿐입니다. 즉 이 말은 비선형관계는 상관계수로 잡아내기 힘들 수가 있습니다. 다음의 예시를 통해 알아보도록 하겠습니다.

library(ggplot2)
x1 = runif(n = 100,min = -10, max = 10)
y = x1 * 10 + rnorm(n = 100, mean = 3, sd = 5)

ggplot() +
  geom_point(aes(x = x1, y= y),size = 3) +
  geom_text(aes(x = 5, y = -30),label = round(cor(x1,y),4)) +
  theme_bw() 

위 두 변수는 산점도로 보나, 상관계수로 보나 거의 1에 가까운 상관계수를 가지고 있습니다. 이는 두 변수의 관계를 하나의 선으로 표현할 수 있다는 의미와 같습니다. 그렇다면 상관분석 검정을 해보도록 하겠습니다.

cor.test(x1,y)

    Pearson's product-moment correlation

data:  x1 and y
t = 113.5, df = 98, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.9943742 0.9974582
sample estimates:
      cor 
0.9962179 

\(p-value\)가 매우 낮은 것을 확인할 수 있습니다. 귀무가설을 기각할 수 있으며, 두 변수의 상관계수는 0이 아니라는 결론을 내릴 수가 있습니다.

y = x1^2 + x1 * 10 + rnorm(n = 100, mean = 3, sd = 5)

ggplot() +
  geom_point(aes(x = x1, y= y),size = 3) +
  geom_text(aes(x = 0, y = 100),label = round(cor(x1,y),4)) +
  theme_bw() 

두번째 산점도에서는 선형보다는 2차 함수꼴의 관계를 가지고 있습니다. 또한 상관계수가 0.8 ~ 0.9로 낮아진 것을 확인할 수 있습니다. 여전히 높은 값인 것은 맞지만 선형관계가 아니기에 상관계수가 감소를 하였습니다. 하지만, 그렇다고 저 두 변수의 상관도가 낮아졌다고 판단하기는 힘듭니다. 선형 관계가 아닐 뿐이지, 비선형 관계는 그대로 유지하고 있습니다.

y = sin(x1) + rnorm(n = 100, mean = 3, sd = 0.3)

ggplot() +
  geom_point(aes(x = x1, y= y),size = 3) +
  geom_text(aes(x = 0, y = 5),label = round(cor(x1,y),4)) +
  theme_bw() 

이번 산점도는 상관계수가 매우 낮게 나오는 것을 확인할 수 있습니다. 근데, 그렇다고 상관이 없다고 할 수 없습니다. 산점도에는 뚜렷히 삼각함수 \(sin(x)\)꼴의 관계를 가지고 있는 것을 확인할 수 있기 때문입니다. 이렇게 상관분석은 정말 편한 분석이지만, 그렇다고 만능인 분석이 아닙니다. 저는 개인적으로 상관계수를 구하기 보다는 산점도를 직접 봄으로써 변수의 상관정도를 파악하고, 이러한 방법이 더 정확합니다.

10. 단순 선형 회귀분석의 추정

• 회귀분석 : 인과관계를 가지고 있는 두 변수간의 함수관계를 통계적으로 규명하고자 하는 분석

  분산분석과 회귀분석은 선형모형이라는 큰 줄기에서 같은 방법론이라는 말씀을 드린 바 있습니다. 두 모형 모두 예측자(predictor)에 따른 평균 반응값을 추정 혹은 예측하는 모형으로, 주어진 데이터를 통해서 선형 모형을 설정하고 새로운 값에 대한 반응값을 예측하는 것에 그 목적이 있습니다. 예측자란 반응 값을 예측하기 위해 사용되는 것으로 설명 변수(explanatory variable)와 혼용되는 개념으로 이해하시면 될 것 같습니다. 또한 아노바와 회귀모형의 차이점이 있다면 아노바와 달리 회귀분석은 일반적으로 연속형 예측자를 가지고 있는 경우에 사용된다는 점입니다. 회귀분석은 설명 예측자와 반응 값의 수에 따라 단순 회귀분석, 다중 회귀분석 등으로 분류되며 설명 예측자와 반응 값의 관계에 따라 선형 회귀분석, 다항 회귀분석(비선형 회귀분석) 등으로 나누어 집니다. 먼저 선형 회귀분석에 대해 알아보도록 하겠습니다.

  일반적인 회귀분석의 가장 핵심 관점은 두 변수간의 관계를 선형적으로 바라보는 것입니다. 그렇기에 관계를 확인하려는 두 변수가 선형적이지 않다면 애초에 성립될 수 없는 것이 회귀분석입니다. 이러한 선형성은 분석 전에 산점도와 같은 그래프를 통해 확인하는 것이 일반적이며 선형 회귀분석 전체를 아우르는 중요한 가정이라고 할 수 있습니다.

회귀분석의 최종 목표는 두 변수 사이의 선형성이 존재한다는 가정하에 그 선형관계를 대표할 수 있는 하나의 직선을 구하는 것입니다. 예를 들어 Y와 X라는 두 변수 사이의 관계가 \(Y= a+bX\) 라는 직선으로 대표될 수 있다면, 두 변수 사이의 관계를 한 눈에 알 수 있는 것은 물론이고 새로운 \(X\)값이 주어졌을 때 \(Y\)값을 손쉽게 예측할 수 있을 것입니다. 여기서 고민해 보아야 하는 것은 \(X ,Y\) 관찰치가 주어졌을 때, 그 관계를 표현하는 \(a\) 와 \(b\) 를 어떤 방법으로 추정해야 가장 효율적이고 두 변수간의 관계를 정확히 나타낼 수 있을까에 대한 방법론적 문제입니다.

회귀분석에서는 잔차의 제곱합을 최소로 만들어주는 \(a\)와 \(b\)를 추정하는 방법을 택하였고 이를 최소제곱 추정법(method of least squares estimation)이라 부릅니다. 여기서 잔차(residual)는 실제 관측값과 예측된 값의 차이를 말하며 그 제곱합을 최소로 한다는 것은 해당 선과 실제 관측값의 차이의 총량을 최소로 하겠다는 것을 뜻합니다. 다음 그림을 참고하면 어렵지 않게 이해하실 수 있을 것입니다.

파란색 직선이 적합된 회귀선이고 검은 점들이 각 관찰치입니다. 이 경우 실제 관찰치와 적합된 회귀선의 차이가 바로 잔차이고 그 제곱은 위와 같이 사격형의 면적으로 표현될 것입니다.

잔차를 그냥 합하면 + / - 가 상쇄되어 의미가 없어지지만, 이런식으로 제곱 후에 총합 이용하면 실제 관찰치와 적합된 회귀선 사이의 차이를 표현해줄 수 있고 이것은 회귀선이 얼마나 관찰치들을 잘 대표하고 있는지에 대한 척도로 생각할 수 있습니다. 이런 관점에서 봤을 때 최소제곱법, 즉, 이 잔차제곱합을 최소화하는 회귀선은 충분히 합리적인 회귀선이라고 볼 수 있을 것입니다.

회귀선의 평가

최소제곱법을 이용하면 어떤 자료에서든 회귀선을 추정할 수 있고, 최소제곱법이 아무리 합리적인 방법이라고 하여도 관찰치들을 절대적으로 잘 대표하는 직선을 항상 제시하지는 못합니다. 그래서 회귀선을 적합시킨 후에는 꼭 그 회귀선이 얼마나 정확한지에 대한 검정을 해주어야 합니다. 이를 회귀선의 적합도 혹은 정도를 평가한다고 하며 분산분석의 아이디어를 그대로 가져와서 진행하게 됩니다.

• \(Y_i (관찰치)\), \(\overline{Y}\)(관찰치의 평균), \(\hat{Y_i}\)(회귀선에 의한 추정값)을 계산한다.

• 총 편차(\(Y_i-\overline{Y}\))를 회귀선을 기준으로 두 영역으로 나누어준다.

• 회귀선에 의해 설명이 가능한 영역 : \(\hat{Y_i}-\overline{Y}\)

• 회귀선에 의해 설명이 되지 않은 영역 : \(Y_i-\hat{Y_i}\)

• 모든 관찰치에 대해 제곱합을 계산한 후, 직교분해를 한다.

\[ \sum_{i=1 }^n (Y_i - \overline Y)^2 \; = \; \sum_{i=1 }^n (\hat Y_i - \overline Y)^2 \; + \; \sum_{i=1 }^n (Y_i - \hat Y_i)^2\\ SST \qquad = \qquad SSR \qquad + \qquad SSE \]

• 분산분석과 마찬가지로 분산분석표를 작성하여 회귀선의 유의성 검정을 한다.

분산분석과 마찬가지로 회귀선으로 설명하는 상대적인 비중을 나타내는 \(MSR\)이 모형의 잔차(Error)를 설명하는 \(MSE\)보다 상대적으로 커야지 회귀선이 유의할 수 있습니다. 이 제곱합을 비교해주기 위해 마찬가지로 \(F\)분포를 이용합니다.

여기서 많은 분들이 SSE의 자유도를 이해하기 어려워 하십니다. 이는 결정계수(\(R^2\))와도 관련이 있는 부분인데, 이 부분에 대하여 간단하게 설명하고 넘어가도록 하겠습니다.

자유도는 자유롭게 선택될 수 있는 자료 수로, 온전히 해당 통계량에 대한 정보를 생성하는 데 사용되는 자료 수라고 했습니다. 만약 관찰치가 2개 밖에 없다면 회귀선은 무조건 그 두 점을 지나야 합니다. 그래야 직선이 만들어지기 때문입니다. 이 경우 두 선을 지나는 회귀선을 추정하기만 할 뿐, 그 회귀선을 평가할 자유도는 없습니다. 같은 관점으로 관찰치가 3개라면 2개는 위와 같은 논리로 사용되고 회귀선의 정도를 평가하는 데에는 하나의 관찰치만을 사용하게 될 겁니다. 즉, 기본적으로 회귀선은 2개의 자유도를 차지합니다. 그렇기 때문에 회귀선 정도 평가를 위해 사용되는 \(SSE\)는 \(n-2\)의 자유도를 갖게 됩니다.

설명 예측자가 2개 이상인 다중회귀분석에도 그대로 적용됩니다. 예를 들어 2개인 경우 회귀식은 3차원의 좌표공간에서 어떤 평평한 면으로 표현할 수 있습니다. 이 면은 최소한 3개의 점이 있어야 형성될 수 있습니다. 그렇기 때문에 관찰치가 3개인 경우,적합된 회귀식을 무조건 그 세 점을 지나갈 것이며, 회귀식의 정도를 파악하고 싶다면 최소한 4개 이상의 관찰치가 필요합니다.(3개의 점만 있다면 회귀면은 모든 점을 지나갈테고 이는 회귀선 평가에 관해서 아무 의미가 없습니다.) 즉, 이 경우 회귀식은 3의 자유도를 차지합니다. 회귀식 평가를 위한 \(SSE\)는 n-3의 자유도를 갖게 됩니다.

이를 일반화 하면 SSE는 n-예측자 수-1 가 됩니다. 일반적으로 예측자 수를 k로 표현하니, n-k-1이라고 할 수 있습니다. 이와 같은 이유로 회귀선을 평가하는 데 사용되는 가설과 검정통계량은 다음과 같습니다.

\[ H_0 : 추정된 \; 회귀식은\; 유의하지\;않다 \\ \\ H_1 : 추정된\; 회귀식은 \;유의하다 \\ \; \\ \]

\[ F_0=\frac {SSR/df_R } { SSE/df_E} \sim F (df_R,df_E) \\ \; \\ \]

\[ df_R=k \qquad df_E=n-k-1\\ \; \\ n:number\;of\; observations \qquad k:number\;of\;predictors \]

만약 MSE에 비해 MSR이 상대적으로 큰 값을 가지고 있다면 \(F_0\)는 높은 값을 갖겠고, 이는 회귀식으로 설명할 수 있는 변동이 그렇지 않은 변동에 비해 크다는 것을 의미합니다. 이 경우 귀무가설을 기각합니다.

\(R^2\) 계산

회귀선이 전체 변동 중 설명하는 비중을 회귀선의 설명력, \(R^2\)라고 합니다. \[ R^2 = \frac {SSR} {SST} = 1- \frac {SSE} { SST} \]

이는 결정계수(coefficient of determination, \(R^2\))라는 척도로 전체 변동중 회귀선에 의해 설명되는 변동이 차지하는 비중을 표현하는 값입니다. 0과 1 사이의 값을 가지며 ‘회귀식의 기여율’ 이라는 관점으로 해석될 수 있습니다. 이 역시 위의 F검정과 같은 원리이며, 같은 방법으로 회귀선을 평가하는 척도입니다.

회귀모형 설정

회귀모형을 설정해보도록 하겠습니다. 여기서 \(X\) 값은 주어진 값으로 취급하므로 오차항 \(\epsilon\) 의 분포는 곧 반응값 \(Y\)의 분포입니다. \[ Yi = \beta 0 + \beta1 Xi + \epsilon_{i} \; \\ \; \\\quad \epsilon_{i} \sim iid \; N (0,\sigma^2) \]

분산분석에서와 같이 오차항 \(\epsilon\)에 대한 가정은 곧 회귀분석을 하기 위한 가정이고 \(iid \; N\) 에 따라 독립성, 정규성, 등분산성 세 가지 성질을 요합니다.

이 세 가정은 오차항의 추정량인 잔차를 통해 진단되며 잔차가 관찰치에 따라 특정한 패턴을 보이거나 커졌다 작아졌다 하는 양상을 보인다면 독립성, 등분산성을 의심해보아야 하며 잔차의 히스토그램을 그려보거나 관련 검정을 통해 정규성을 만족하는지에 대해 확인해 보아야 합니다.

또한 실제 데이터를 기지고 추정한 회귀식은 다음과 같이 표현합니다.

\[ \hat Y_i = b_0 + b_1 X_i \]

절편 \(b_0\) 는 예측자가 0의 값을 가질 때의 반응 값을 의미하며 초기값을 나타냅니다. 기울기 \(b_1\)은 예측자가 가지고 있는 반응값에 대한 단위 당 영향도를 의미합니다. 기울기가 양의 값을 갖고 있다면 예측자가 증가할수록 반응 값 역시 증가할 것이며 음의 값을 갖는다면 예측자의 증가에 따라 반응 값은 감소할 것입니다. 또한 절댓값이 클수록 그 영향이 크다는 것을 의미합니다.

이 영향도 역시 유의한지 검정을 진행할 수 있습니다. 이를 회귀계수의 검정 혹은 회귀식의 기울기 검정이라고 합니다. 기울기 검정은 \(t\)분포를 이용하게 됩니다. 만약 정규성의 가정이 만족하면, 회귀계수는 \(t\)분포를 따르기 때문이죠. 회귀계수가가 0인가를 확인하는 검정은 해당 설명변수가 반응 값에 영향을 미치는지 아닌지를 확인 하는 것과 같습니다. 여러 설명변수가 포함된 모형에서도 진행할 수 있으며, 마찬가지로 각 변수의 독립적인 영향력을 평가하게 됩니다.

\[ H_0:\beta_1=0 \qquad 회귀계수는\; 0이다.\\ H_1:\beta_1\neq0\qquad 회귀계수는\;0이\; 아니다. \]

11. 회귀분석(R Code)

회귀분석은 제가 만들어 둔 데이터로 진행을 하도록 하겠습니다.

# 데이터 불러오기
Regression = read.csv("F:\\Dropbox\\DATA SET(Dropbox)/Regression.csv")

**
산점도**

회귀분석은 우선적으로 산점도를 그려보고 선형성을 판단해야됩니다.

library(ggplot2)

ggplot(Regression,aes(x = X , y = y)) +
 geom_point() +
 geom_smooth(method = 'lm') +
 theme_classic()

산점도를 그려본 결과 X와 y는 선형 관계를 가지고 있는 것을 알 수 있습니다. 이런 경우 단순 선형 회귀분석을 진행하면 됩니다.

회귀분석

Reg = lm(y ~ X  , data = Regression)
anova(Reg)

Analysis of Variance Table

Response: y
           Df Sum Sq Mean Sq F value    Pr(>F)    
X           1 143757  143757  3882.2 < 2.2e-16 ***
Residuals 148   5480      37                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(Reg)

Call:
lm(formula = y ~ X, data = Regression)

Residuals:
    Min      1Q  Median      3Q     Max 
-9.5537 -5.7116  0.2738  5.0961 10.2835 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   1.9839     0.9975   1.989   0.0486 *  
X            10.0924     0.1620  62.308   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 6.085 on 148 degrees of freedom
Multiple R-squared:  0.9633,    Adjusted R-squared:  0.963 
F-statistic:  3882 on 1 and 148 DF,  p-value: < 2.2e-16

\(p-value\)는 0과 매우 가까운 값입니다. 그러므로 회귀분석의 귀무가설(\(H_0:회귀식은\; 유의하지\; 않다\))을 기각하게 됩니다. 즉, 회귀선은 유의합니다.

R분석 결과를 자세히 보도록 하겠습니다.

• Residuals는 \(y_i-\hat{y_i}\)의 요약 값을 나타냅니다.

• Coefficients는 추정된 회귀식을 나타냅니다.

• Estimate는 회귀식의 절편과 기울기를 나타냅니다.

• (Intercept) : \(b_0\)의 값은 1.9839입니다.

• \(X\) : \(b_1\)의 값은 10.0924입니다.

• Pr(>|t|)는 각각 \(b_0\), \(b_1\)에 대한 가설검정 결과를 나타냅니다.
  만약 설명 변수를 여러 개 사용하는 다중 회귀분석을 진행할 경우, 각 설명 변수에 대한 가설 검정을 진행하게 될 것입니다.

• Multiple R-squared는 회귀모형의 \(R^2\)를 나타냅니다.

• Adjusted R-squared는 조금은 다른 \(R^2\)를 나타냅니다.

• 일반적으로 회귀 모형은 설명 변수가 많을 수록 설명력은 높아지는 구조를 가지고 있습니다.

• 만약 다른 변수로 구성 된 두 회귀 모형의 설명력을 비교하고자 할 때 Adjusted R-squared를 사용합니다.

\[ Model_1:\hat{y_i}=b_0+b_1x_{1i}+b_2x_{2i}\\ Model_2:\hat{y_i}=b_{0^*}+b_1x_{1i}+b_2x_{2i}+b_3x_{3i} \]

예를 들어 회귀식이 2개가 있고 두 모형 중 \(R^2\)가 더 높은 모형을 택하고자 합니다.
\(Model2\)는 3개의 설명 변수가 투입되었기에 설명 변수를 2개를 투입한 \(Model1\)에 비해 우위에 있습니다. \(R^2\) 의 구조상 설명변수가 많아지면 항상 오를 수 밖에 없는 구조이기 때문입니다. 그렇지만 설명 변수가 많은 것은 항상 바람직한 것이 아닙니다. 그만큼 많은 계수를 추정해야하기 때문이죠. 그렇기에 설명 변수가 더 많이 투입된 부분을 보정해주어 더 바람직한 모형의 기준을 제시해 주는 \(R^2\)값이 Adjusted R-square입니다.

• p-value : 회귀모형에 대한 유의성을 나타냅니다.

회귀분석의 결과를 정리하면 다음과 같습니다.

• 적합된 회귀선은 \(\hat{y_i}=1.9839 +10.0924x_i\)입니다. \(x_i\)가 한 단위 증가하면 \(\hat{y_i}\)는 10.0924만큼 증가합니다.
• \(R^2\)는 0.963으로 96%의 설명력을 가진다고 볼 수 있습니다.

잔차 진단

회귀분석을 진행한 후에는 항상 가정 검토를 위해 잔차진단을 해야합니다. 회귀분석에서의 가정은 정규성, 등분산성, 독립성입니다. 그런데 독립성은 데이터가 수집되는 과정에서 확인하여야 할 가정이므로 수집된 데이터로 분석이 진행된 상황에서는 확인할 수 없습니다. 그렇기에 여기서는 정규성, 등분산성을 다루고, 추가적으로 영향점이라는 것을 다루도록 하겠습니다.

plot(Reg)

)))

• 정규성

qqnorm(Reg$residuals)

회귀분석에서 잔차의 정규성을 진단하는 이유는 신뢰구간 추정과 가설검을 정확하게 하기 위해서 진행합니다. 우측 상단에 배치되어 있는 Normal Q-Q를 보면 됩니다. x축은 Theoretical Quantiles이며, y축은 Standardized Residuals입니다. QQ plot은 역확률을 이용하여 두 분포를 비교하는 plot입니다. 여기서는 정규분포와 잔차의 분포를 비교하고 있습니다. 만약 선이 대각선을 따라가는 일직선이라면 그만큼 잔차의 분포는 정규분포와 비슷하다고 볼 수 있습니다. 반대로 선이 대각선과는 멀어지는 곡선 형태라면 그만큼 잔차의 분포는 정규분포와 다르다는 것을 의미합니다.

• 등분산성

  등분산성은 회귀분석에서 매우 중요한 가정 중 하나 입니다. 여기서 등분산의 주체는 오차입니다. 그렇지만 실제 오차를 정량화 할 수 없으니 오차의 추정치로써 잔차를 사용하게 됩니다. 잔차는 추정된 회귀선과 실제 값의 차이입니다. 즉, 등분성을 보는 것은 선과 점 사이의 거리가 패턴이 없이 일정한가를 보는 것과 같습니다.

이를 확인하기 위해 조금 극단적인 예를 보겠습니다.

좌측 산점도와 우측 산점도가 있을 때, 두 산점도에 회귀선을 적합시켜보도록 하겠습니다.

우측의 회귀선은 직관적으로 판단해도 회귀선에 문제가 없습니다. 하지만 좌측 회귀선은 그러지 못합니다. 그 이유는 회귀선과 데이터의 차이(잔차)가 x가 커지면서 같이 늘어나고 있기 때문입니다. 즉, 등분산성을 만족한다고 할 수 없습니다. 이 때는, 회귀선이 데이터를 잘 설명한다고 보기 어렵습니다. 이처럼 회귀분석에서 등분산성이 위배되면 회귀분석은 데이터를 설명하지 못한다고 판단하기 때문에 좋은 회귀선이라고 할 수 없습니다. 좌측의 회귀선에 대한 잔차의 등분산 진단 그래프로 보면 다음과 같은 플롯을 확인할 수 있습니다.

x축은 Fitted value입니다. 즉, \(\hat{y_i}\)를 의미합니다. y축은 Resiuals(\(y_i-\hat{y_i}\))을 의미합니다. 잔차가 추정값이 증가함에 따라 증가하는 패턴을 보이고 있습니다. 위에서 말씀드린 것처럼 등분산성을 만족하지 못한다는 것을 알 수 있습니다.
다시 데이터로 돌아와서 위에서 그린 잔차진단 그래프를 확인해보겠습니다. 등분산성은 좌측 상,하 그래프를 보면 됩니다.

좌측은 y축이 Residuals, 우측은 y축이 표준화 잔차값입니다. 이는 또한 설명변수와 반응 값과의 관계를 확인 할 수도 있습니다. 만약 패턴이 존재하지 않는다면 상관 없지만 Fitted value 에 따라 잔차가 특정한 패턴을 보인다면 다른 분석을 고려하거나 설명변수의 변환을 생각해 보아야합니다.

• 독립성

ggplot(NULL) +
  geom_point(aes(x = 1:nrow(Regression), y = Reg$residuals)) +
  geom_hline(yintercept = 0, linetype = "dashed",col = 'red') + 
  xlab("Index") + ylab("Residuals") + 
  theme_bw()

잔차의 독립성이란, 잔차가 ’자기상관(Auto correlation)’이 있는지 없는지 판단하는 것입니다. 우리는 확률변수 및 확률분포를 다룰 때, 모든 확률실험은 독립시행을 했다고 가정을 하였습니다. 이 말은 전 시점에서의 확률실험이 현재의 확률실험에 영향을 주지 못하며, 지금의 확률실험은 미래의 결과에 영향을 주지 못한다는 것입니다. 이를 ’자기상관’이라고 합니다. 대표적으로 주식, 날씨 등이 해당이 됩니다. 만약, 독립성이 위배가 된다면 이 떄는 시계열 분석(Time Series)방법론을 이용하여 회귀분석을 해야 합니다. 검정방법은 매우 간단합니다. Durbin-Watson검정을 하는 방법도 있지만, 잔차들을 시점 순서대로 그래프를 그린다 다음, 패턴이 없다면 독립성을 충족한다고 판단할 수가 있습니다.

• 영향점

회귀분석과 같은 통계모형에는 영향점이라는 것이 존재합니다. 영향점은 쉽게 설명하면 회귀선을 자기 자신한테 당겨오는 데이터입니다. 즉 회귀선의 무게중심과 멀리 있어, 이에 영향을 준다고 볼 수 있습니다. Cook’s distance라는 값을 계산하여 영향점을 판단합니다. plot에 찍힌 92, 134, 144번째 데이터가 회귀선에 영향을 주고 있음을 나타내고 있습니다. 이 점들은 회귀선을 기울기에 크게 영향을 주어 예측에 왜곡을 줄 수 있습니다. 영향점을 발견하였을 경우, 먼저 영향점에 대한 확인이 필요합니다. 영향점은 일반적인 x와 먼 값을 갖고 있으면서 높은 잔차를 보이므로 우선적으로 x값이 정당한 값인가를 확인해야 합니다. 측정 자체가 잘못되었을 가능성이 있기 때문입니다. 또한 예측에 목적이 있다면 해당 데이터를 제거하고 회귀분석을 진행하는 방법도 고려해야 합니다.

12. 다중 회귀분석(Multiple Regression)

• 다중 회귀분석 : 단순 선형회귀분석의 확장판으로 예측자가 2개 이상 쓰이는 경우

  다중 회귀분석은 예측자를 2개 이상 쓰는 경우로, 회귀분석과 거의 동일하다고 볼 수 있습니다. 식 표현은 행렬식을 이용해 표현을 하는데, 이 책의 취지와는 맞지 않으므로 간단하게 다중 회귀분석을 진행할 때 주의해야할 점들에 대해 다루면서 진행하겠습니다.

\[ \hat{y_i}=b_0+b_1x_{1i}+b_2x_{2i}\\ \]

회귀식이 위 식처럼 구해져 있을 때, 회귀식의 해석은 다음과 같이 진행합니다.

\(x_{1i}\)가 1 단위 증가하면 \(\hat{y_i}\)는 \(b_1\)만큼 변한다.(단, \(x2_i\)는 고정)
\(x_{2i}\)가 1 단위 증가하면 \(\hat{y_i}\)는 \(b_2\)만큼 변한다.(단, \(x1_i\)는 고정)

여기서 중요한 점은 서로 다른 예측자는 상관성이 적어야 한다는 것입니다. 만약 \(x_1\)과 \(x_2\)의 상관성이 높다면, \(x_1\)이 증가하였을 때, \(x_2\)는 고정되지 못하고 함께 증가해버리는 문제가 발생해버립니다. 이런 상황을 다중공선성(Multicolinearity)라고 합니다. 회귀식에 다중공선성이 존재할 경우, 회귀식의 분산이 팽창을 하게 되며, \(b_1\), \(b_2\)에 대한 추정이 불확실하게 됩니다. 그러므로 다중공선성이 존재할 경우, 회귀식에 투입되는 예측자를 다시 조정하여 회귀식을 구해야합니다. 이 부분은 후에 주성분 분석에서 다루도록 하겠습니다.

x1 = runif(n = 100, min = -10, max = 10)
x2 = 0.7*x1 + rnorm(n = 100,mean = 0,sd = 1)

y = 1.3*x1 + rnorm(n = 100, mean = 6, sd = 3)

DF = data.frame(
  x1 = x1,
  x2 = x2,
  y = y
)

library(GGally)

ggpairs(DF)

위 그림을 보시면 \(x_1\), \(x_2\)와의 상관관계가 매우 높은 것을 확인할 수가 있습니다. 이러한 경우, 회귀분석을 구했을 때, 다중공선성이 매우 의심되는 상황이 발생합니다.

다중 회귀분석

M_Reg = lm(y ~ x1 + x2)
summary(M_Reg)

Call:
lm(formula = y ~ x1 + x2)

Residuals:
    Min      1Q  Median      3Q     Max 
-7.6871 -1.8648  0.2094  1.6262  6.5584 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   6.2825     0.2965  21.192  < 2e-16 ***
x1            1.1548     0.2202   5.244 9.18e-07 ***
x2            0.2368     0.3083   0.768    0.444    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.937 on 97 degrees of freedom
Multiple R-squared:  0.8793,    Adjusted R-squared:  0.8768 
F-statistic: 353.4 on 2 and 97 DF,  p-value: < 2.2e-16

언뜻 보기에는 회귀분석에 큰 문제는 없어보입니다. 다중회귀분석에서는 \(x_1\), \(x_2\)의 기울기 \(b_1\), \(b_2\)에 대한 가설 검정이 진행이 됩니다. \(x_1\)의 회귀계수 \(b_1\)은 유의한 것으로 나왔지만, \(x_2\)의 회귀계수 \(b_2\)는 유의하지 않은 것으로 나왔습니다.

다중공선성 진단

library(car)

vif(M_Reg)

      x1       x2 
19.66709 19.66709 

다중공선성은 분산팽창지수(VIF)를 측정함으로써 판단이 가능합니다. 기본적으로 5보다 작으면 다중공선성이 없다고 판단할 수 있습니다. 지금은 다중공선성이 매우 높은 상황이므로 다중공선성이 존재하고 있다고 확인할 수 있습니다.

13. 다항 회귀분석(Polynomial Regression)

• 다항 회귀분석 : 예측자들이 1차항으로 구성된 것이 아닌, 2차항, 3차항 등으로 구성되어 있는 회귀식
  \[ \hat y = b_0+b_1x_i+b_2x_{i}^2+\cdots+b_px_p^{p} \]

  다항 회귀분석은 위 식처럼 구성이 될 수 있습니다. 다항회귀분석에서는 매우 중요한 개념이 하나 따라오는데, 이를 확인하고 다항 회귀분석을 진행하도록 하겠습니다.

분산-편차의 Trade off 관계

• Trade off : 두 개의 목표 중에서 하나를 달성하려고 하면 다른 목표가 희생되어야 하는 관계를 의미합니다.

기계학습에서 예측 모형을 만드는 것은 항상 Trade off 관계를 생각해야 됩니다. 기본적으로 통계학에서는 모형의 Target Variable(종속 변수)이 연속형(Continuous)일 때는 MSE 와 Bias에 주목합니다. 만약 Target Variable이 범주형(Categorical)일 경우에는 모형의 Error Rate에 주목합니다. 그 이유는 모형의 정확성은 MSE 혹은 Bias가 얼마나 작은지에 따라 결정되기 때문입니다.

MSE(Mean Squared Error)

앞서 단순선형회귀에서 MSE를 다루었지만, 다시 한번 복습하면서 재차 다루어보도록 하겠습니다.

\(MSE\)를 이해하기 위해서는 위 그림의 의미를 제대로 숙지하고 있어야 합니다.

\[ Y_i = 실제 관측값 \, \]

\[ \hat{Y_i} = 예측값 \]

\[ \overline{Y} = 평균 \]

여기서 예측값 \(\hat{Y_i}\)은 추정된 회귀식 \(\hat{Y_i} = \beta_0 + \beta_1X_i\)으로부터 추정된 예측 값입니다. 평균값 \(\overline{Y}\)은 \(\overline{Y} = \frac{1}{n}\Sigma(Y_i)\) 평균 산술식으로부터 계산된 표본평균입니다.

보라색 간격에 해당되는 \(\hat{Y_i}-\overline{Y}\)는 간격의 차이를 설명할 수 있기 때문에(\(\beta_0 + \beta_1X_i - \frac{1}{n}\Sigma(Y_i)\)) 추정된 회귀식이설명이 가능한 영역이 입니다. 하지만 초록색 간격에 해당되는 \(Y_i - \hat{Y_i}\)는 실제로 관측된 \(Y_i\)값이 왜 저기에 찍혔는지에 대해서는 설명을 할 수 없습니다. 그런관계로 해당 영역을 설명이 불가능한 영역입니다.

모든 모형은 설명력이 높으며 (혹은 설명 못하는 영역이 적은) 예측이 잘 되는 모형이 좋은 모형입니다. 결국 모형의 결함은 설명을 하지 못하는 \(Y_i - \hat{Y_i}\)은 오차로 계산하게 됩니다. 이를 잔차(Residuals, 혹은 Error)라고 합니다.

추정된 회귀식 \(\hat{Y_i} = \beta_0 + \beta_1X_i\)이 변수들 간의 인과관계를 제대로 설명하는지 측정하기 위하여 잔차의 합을 계산하게 됩니다. 여기서 잔차는 상황에 따라 양수가 될 수도, 음수가 될 수도 있습니다. 그렇기에 부호 간 계산으로 잔차의 합이 상쇄되는 것을 방지하기 위하여 잔차를 제곱하여 합을 구하게 됩니다. 이를 오차의 제곱합(Sum Squred Error, \(SSE\))라고 부릅니다. 그런 다음 계산된 \(SSE\)를 보정하기 위하여 \(SSE\)를 오차의 자유도(\(df_e\))로 나눕니다. 그렇게 계산 된 값을 오차 평균 제곱합(Mean Squared Error, \(MSE\))라고 부릅니다.

\[ SSE = \Sigma(Y_i - \hat{Y_i})^2 \]

\[ MSE = \frac{1}{df_E}\Sigma{(Y_i-\hat{Y_i})^2} \] SSE를 오차의 자유도로 나누어주는 이유

• 제곱 된 값은 항상 양수입니다
• 양수를 모두 더하게 되면, 데이터가 많을 수록 값은 커지게 됩니다.
• 그 의미는 SSE자체가 데이터가 많을 수록 단순히 커지는 의미이기 때문에, 정말 오차가 높은가? 에 대한 평가기준이 잘못 해석될 수가 있습니다.
• 이를 자유도로 나눔으로써 평균이 계산되고, 보정된 평균오차를 모형의 Error 수준으로 판단하게 됩니다.

위와 같은 계산으로 회귀식이 설명 가능한 영역인, \(\hat{Y_i} - \overline{Y}\)은 각각 다음과 같이 계산이 됩니다. \[ SSR = \Sigma(\hat{Y_i}-\overline{Y})^2 \]

\[ MSR = \frac{1}{df_R}\Sigma(\hat{Y_i}-\overline{Y})^2 \]

그렇다면 추정되는 회귀식이 데이터의 인과관계를 얼마나 잘 설명하는지 계산하기 위해 설명을 하지 못하는 영역 대비, 설명을 할 수 있는 영역을 비교하게 됩니다. \[ F \, value = \frac{MSR}{MSE} = \frac{\frac{1}{df_R}\Sigma(\hat{Y_i}-\overline{Y})^2}{\frac{1}{df_E}\Sigma{(Y_i-\hat{Y_i})^2}} \] \(\frac{MSR}{MSE}\)는 두 집단의 분산을 비교하는 F 분포를 따르게 됩니다.

추정된 회귀식의 설명하는 영역이 설명하지 못하는 영역에 비해 얼마나 큰지 나타내는 검정통계량(Test Statistics)\(F \, value\)는 값이 클수록 회귀식의 귀무가설(\(H_0\,: 회귀식의 \, 기울기가 \,0이다\))를 기각할 가능성이 커지게 됩니다. 그렇다면 \(F\ value\)값이 커질려면 다음과 같습니다.

• \(MSR\)이 증가
• \(MSE\)가 감소

분석모형의 성능 평가를 \(MSE\)로 하는 이유입니다. \(MSE\)가 작은 모형일수록 회귀식의 오차가 줄기때문에 그만큼 현상을 잘 설명한다고 할 수 있습니다.

Variation & Bias

회귀식으로 추정된 \(\hat{Y_i}\)는 얼핏보면 단일 값인 점 추정(Point Estimation)으로 생각할 수 있지만, 모든 통계분석 모형은 구간 추정(Interval Estimation)입니다. 구간 추정이란 소리는 추정값에 대한 신뢰구간을 계산한다는 의미입니다.

신뢰구간의 의미는 똑같은 \(X\)값이 주어졌을 때, 추정값 \(\hat{y} = \beta_0 + \beta_1 X\) 의 값이 \([\hat{y} - \alpha , \hat{y} +\alpha]\)의 범위에 속한다는 의미입니다. 만약 분산이 크다면, 이 신뢰구간의 길이는 길어지게 되고, 추정의 신뢰성이 떨어지는 문제가 발생합니다. 반대로 편의(Bias)는 추정된 값이 모집단의 특성, 즉 모수를 반영하지 못한다는 의미입니다.

위 그림을 보시면 분산과 편향이 크고 작을 때에 따라 모형의 정확성이 어떻게 변하는지 알 수 있습니다. 모수의 True Value가 원 정중앙에 있다고 하였을 때, Variance 가 크다는 것은 추정값의 범위가 넓은 것을 의미하고, Bias가 크다는 것은 영점조준 사격 훈련 때 탄집군은 생겼지만 영점이 잘못잡혔다와 비슷하다고 생각하시면 됩니다.

선형 & 비선형 Modeling

Linear Regression(선형 회귀분석)과 Non - Linear Regression(비선형 회귀분석)을 잠깐 다루고 가겠습니다.

library(ggplot2)

ggplot(Regression) +
  geom_point(aes(x = X, y = y),col = 'royalblue',alpha = 0.4) +
  geom_smooth(aes(x = X, y = y),col = 'red') +
  theme_bw() +
  xlab("") + ylab("")


ggplot(Regression) +
  geom_point(aes(x = X, y = y2),col = 'royalblue',alpha = 0.4) +
  geom_smooth(aes(x = X, y = y2),col = 'red') +
  theme_bw() +
  xlab("") + ylab("")

)

사람들이 회귀분석을 돌릴 때, 가장 실수하는 부분은 단순하게 변수들 간의 상관관계만을 파악해서 분석하는 경우입니다. 상관관계는 두 변수의 관계가 선형성을 띄는지를 판단하는 것일 뿐입니다. 만약 두 변수가 비선형 관계에 있을 경우, 상관관계는 낮게 잡힐 수도 있습니다. 하지만 상관계수가 낮게 잡힌다고 해서 이 두 변수 간에 관계가 존재하지 않는 것은 아닙니다. 비선형으로 회귀식을 잡으면 충분히 관계를 설명할 수가 있게 됩니다.

ggplot(Regression) +
  geom_point(aes(x = X, y = y3),col = 'royalblue',alpha = 0.4) +
  geom_smooth(aes(x = X, y = y3),col = 'red') +
  theme_bw() +
  xlab("") + ylab("")

특히 이런 경우, 데이터의 상관관계수가 매우 낮게 계산이 됩니다. 상관계수만 보면 매우 낮기때문에 일반적으로 모델링 할 생각부터 안하게 됩니다. 하지만 분석모형에서는 이러한 비선형 관계들로 관계식을 추정할 수 있습니다.

• 선형 회귀분석

두 변수의 관계가 선형인 경우에 대해서 회귀분석을 추정해보겠습니다.

LINEAR = lm(y ~ X,data = Regression)
summary(LINEAR)

Call:
lm(formula = y ~ X, data = Regression)

Residuals:
    Min      1Q  Median      3Q     Max 
-9.5537 -5.7116  0.2738  5.0961 10.2835 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   1.9839     0.9975   1.989   0.0486 *  
X            10.0924     0.1620  62.308   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 6.085 on 148 degrees of freedom
Multiple R-squared:  0.9633,    Adjusted R-squared:  0.963 
F-statistic:  3882 on 1 and 148 DF,  p-value: < 2.2e-16

ggplot(Regression) +
  geom_smooth(aes(x = X, y = predict(LINEAR, newdata = Regression)),col = "red",  method = 'lm') +
  geom_point(aes(x = X , y = y),col = 'royalblue') + ylab("") + xlab("") + ggtitle("Linear Regression") +
  theme_bw()

일반적인 선형회귀분석, 즉 선형을 완벽하게 띄고 있는 변수 간의 관계는 간단하게 선형으로 적합시키면 문제가 없습니다.

• Polynomial Regression

변수 간의 관계가 \(y=x^2\)형태를 가지는 데이터에 대해 2차항 회귀분석(다항 회귀분석)을 적용시켜 보겠습니다.

# 제곱꼴 관계를 선형으로 적합
NonLinear = lm(y2 ~ X, data = Regression)
summary(NonLinear)

Call:
lm(formula = y2 ~ X, data = Regression)

Residuals:
    Min      1Q  Median      3Q     Max 
-154.19  -84.79  -46.28   46.35  446.95 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -208.872     21.225  -9.841   <2e-16 ***
X            103.788      3.447  30.113   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 129.5 on 148 degrees of freedom
Multiple R-squared:  0.8597,    Adjusted R-squared:  0.8587 
F-statistic: 906.8 on 1 and 148 DF,  p-value: < 2.2e-16

2차항의 관계를 선형으로 적합하였을 때의 설명력은 85 ~ 86%가 나온 것을 알 수가 있습니다.

다음으로 다항 회귀분석(2차항)을 적용시켜보도록 하겠습니다.

NonLinear2 = lm(y2 ~ poly(X,2), data = Regression)
summary(NonLinear2)

Call:
lm(formula = y2 ~ poly(X, 2), data = Regression)

Residuals:
    Min      1Q  Median      3Q     Max 
-65.180 -25.879   5.213  29.693  39.436 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  345.341      2.582  133.74   <2e-16 ***
poly(X, 2)1 3899.149     31.625  123.29   <2e-16 ***
poly(X, 2)2 1527.866     31.625   48.31   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 31.63 on 147 degrees of freedom
Multiple R-squared:  0.9917,    Adjusted R-squared:  0.9916 
F-statistic:  8767 on 2 and 147 DF,  p-value: < 2.2e-16

회귀식을 \(\hat{y_i} = \beta_0 + \beta_1x_i +\beta_2x_i^2\)형태로 적합한 결과 설명력은 99%로 상승한 것을 볼 수 있습니다.

ggplot(Regression) +
  geom_smooth(aes(x = X, y = predict(NonLinear2,newdata = Regression)),col = "red") +
  geom_point(aes(x = X, y = y2),col = 'royalblue') + 
  ylab("") + xlab("") + 
  ggtitle("Polynomial Regression") +
  theme_bw()

• 유연성이 있는 회귀분석

다음 회귀분석은 \(y=sin(x)\)꼴을 가지는 두 변수 간의 관계를 회귀식으로 추정해보도록 하겠습니다. 워낙 형태가 괴이하기 때문에 몇차 항을 적합시켜야할지 모르겠습니다. 그러하니 변수항의 차수(Degree of Polynomial)을 2 ~ 10까지 주고 Testing을 해보도록 하겠습니다.

# Train Set & Test Set 형성

TRAIN = Regression[1:100,]
TEST = Regression[101:150,]

# 저장공간 생성

DEGREE = c()
TEST_MSE = c()
TRAIN_MSE = c()
Adj_R = c()
TEST_VAR = c()
TRAIN_VAR = c()

# 적합 모델 찾기 

for( degree in 2:10){

  FLEXIBLE_MODEL = lm(y3 ~ poly(X,degree),data = TRAIN)

  ## Summary Save

  SUMMARY = summary(FLEXIBLE_MODEL)
  ANOVA = anova(FLEXIBLE_MODEL)
  DEGREE = c(DEGREE,degree)

  ## R_SQUARE

  Adj_R = c(Adj_R,SUMMARY$adj.r.squared)

  ## Train Set

  TRAIN_MSE = c(TRAIN_MSE, ANOVA$`Mean Sq`[2])
  TRAIN_VAR = c(TRAIN_VAR, var(FLEXIBLE_MODEL$fitted.values))

  ## Test Set

  Pred = predict(FLEXIBLE_MODEL, newdata = TEST)
  TEST_RESIDUALS = (Pred - TEST$y3)
  TEST_MSE_VALUE = sum(TEST_RESIDUALS^2)/(nrow(TEST))
  TEST_MSE = c(TEST_MSE,TEST_MSE_VALUE)
  TEST_VAR = c(TEST_VAR, var(Pred))
}

## Test 결과 데이터프레임 생성

F_DATA = data.frame(
  DEGREE = DEGREE,
  Adj_R = Adj_R,
  TRAIN_MSE = TRAIN_MSE,
  TRAIN_VAR = TRAIN_VAR,
  TEST_MSE = TEST_MSE,
  TEST_VAR = TEST_VAR
)

항차를 2차항부터 10차항까지 차례대로 추정해본 결과는 다음과 같습니다.

library(dplyr)
library(reshape)

ggplot(TRAIN) +
  geom_point(aes(x= X, y = y3) , col = 'royalblue', alpha = 0.8) +
  geom_smooth(aes(x = X, y = predict(FLEXIBLE_MODEL,newdata = TRAIN)),col = 'red') +
  xlab("") + ylab("") + ggtitle("Flexible Regression") +
  theme_bw()


ggplot(F_DATA) +
  geom_point(aes(x = DEGREE, y = Adj_R * 100)) +
  geom_line(aes(x = DEGREE, y = Adj_R * 100)) +
  geom_text(aes(x = DEGREE, y = Adj_R * 100 + 5,
                label = paste(round(Adj_R*100,2),"%",sep="")),size = 3)+ 
  scale_x_continuous(breaks = seq(2,10,by = 1)) + 
  ylab("Adj_R2") +
  theme_bw()


F_DATA %>%
  select(DEGREE,TRAIN_MSE,TEST_MSE) %>%
  melt(id.vars = c("DEGREE")) %>%
  ggplot() +
  geom_point(aes(x = DEGREE, y = value, col = variable)) +
  geom_line(aes(x = DEGREE, y = value, col = variable)) +
  labs(col = "") +
  theme_bw() +
  theme(legend.position = "bottom") +
  xlab("Degree") + ylab("MSE")

F_DATA %>%
  select(DEGREE,TRAIN_VAR,TEST_VAR) %>%
  melt(id.vars = c("DEGREE")) %>%
  ggplot() +
  geom_point(aes(x = DEGREE, y = value, col = variable)) +
  geom_line(aes(x = DEGREE, y = value, col = variable)) +
  labs(col = "") +
  theme_bw() +
  theme(legend.position = "bottom") +
  xlab("Degree") + ylab("Variation")

)))

• \(R^2\)는 6차항부터 급격하게 올라가는 것을 볼 수가 있습니다. 따라서 6차항은 되어야 \(y=sin(x)\)형태의 관계를 잘 설명할 수 있는 편이라고 생각할 수 있습니다.
• \(MSE\)는 Train Set과 Test Set에 따라 추세가 다릅니다. 차수가 높아질수록 Train Set의 \(MSE\)는 감소하는 것을 알 수있습니다. 하지만 Test Set의 MSE는 감소하다가 증가하는 것을 확인할 수 있습니다. 이는 Train Set은 기가막히게 잘 맞추지만 새로운 데이터인 Test Set은 맞추지 못하는 OverFitting이 발생하였다고 볼 수 있습니다.
• \(Variance\)는 항차가 올라갈수록 대체로 증가하는 추세에 있는 것을 볼 수 있습니다. 여기서 고차항의 회귀모형의 단점이 제대로 드러납니다. 분석 모형이 유연할수록(항차가 높을수록) 회귀추정값의 분산은 높게 뛰기 마련입니다. 이는 회귀식에 의한 추정값에 대한 신뢰구간이 길어진다는 의미이며, 결과에 대한 신뢰도가 떨어진다는 것을 의미합니다.

여기까지 비선형 회귀분석에 대해 알아보았습니다.