[Statistics] Analysis

데이터 사이언티스트 되기 책: https://recipesds.tistory.com/

Analysis

이제부턴 분석에 대해 다룰 것이다. 검정에서 봐왔던 $t$-test도 사실 분석의 일종이라 할 수 있다. 하지만 2개 집단 차이를 볼 때 차이 분석이라 부르지 않고 보통 $t$ 검정이라 한다. 이제 다룰 ANOVA는 분석이라 하는데 $t$ 검정이라 불리는 이유는 검정을 하기 전에 표본 분석을 통한 통계량을 구하기 위해 적당한 Grouping이나 데이터 변환등을 동원해서 데이터 처리를 하는 단계가 있는데, 차이 분석의 경우는 표본을 구하기만 하면 곧바로 검정이 가능하므로 특별하게 분석이라 부르지 않고 그냥 차이 검정 또는 $t$ 검정이라 부른다. 분석의 예를 들면, ANOVA는 분산의 아이디어로부터 표본(원시 (Raw) 데이터)을 분산의 형태로 변환해서, F검정을 하고, 교차분석은 교차표라는 형태로 데이터를 Agrregation 한 후에 카이제곱 검정을 한다.

ANOVA

ANOVA(Analysis of Variance)는 분산분석으로 3개 이상의 집단 간 평균의 차이를 검정할 때 사용한다. 평균의 차이를 검정하는데 이름이 분산분석인 이유는 분산의 성질과 원리를 이용한 평균차이 분석이기 때문이다. ANOVA는 평균을 직접 비교하지 않고, 집단내의 분산과 각 집단 간의 평균의 분산을 이용해 평균이 다른지 확인하는 방법이다.

3개의 집단이 있다고 할 때, 위 그림을 통해 알 수 있듯이 집단 평균 간의 분산이 크고, 집단 내의 분산이 작으면 평균이 서로 확실히 다르다는 것을 알 수 있다. 위 그림에서는 3번 케이스를 제외하고는 평균이 다르다고 단정 지을 수 없다. ANOVA의 기준은 집단간 평균의 분산이 커지고, 집단 내 분산이 작아지는 3번의 경우를 기준으로 찾아낼 수 있고, 이러한 경향이 커질수록 통계량 값은 커진다. 따라서 평균 제곱 간의 비를 (집단 간 평균의 분산 / 집단 내 분산) 으로 정의할 수 있다. 이때의 검정 통계량 $F$는 $F$ 분포를 따르고, 이 차이가 통계적으로 유의한 지를 분석해서 평균이 모두 같다는 귀무가설을 검증하면 된다.

분산 자체는 $\chi^2$ 카이스퀘어 분포를 따르는데, 분자 분모가 모두 카이스퀘어 분포를 따르는 경우에 이런 비율이 따르는 분포는 $F$ 분포를 따른다. 따라서 ANOVA에서는 $F$ 분포를 사용한다. 우리가 이용하는 F분포는 Null Hypothesis가 집단간 평균이 같을 경우의 F값(통계량)의 분포이다. 만약 F분포에서 관측된 F값(통계량)으로 p값을 계산해 보았는데 p값이 너무 작은 경우에는 평균이 서로 같은 환경에서 관측하기 어려운 것이므로 집단 중에 최소한 1개는 평균이 다르다라고 결론을 낼 수 있다.

통계량은 앞서 (집단 간 평균의 분산 / 집단 내 분산)를 (Between Mean Variance/ Within Variance)로 표현한다. 이 의미는 Between Mean Variance가 Within Variance에 비해 크면 최소한 1개 그룹은 차이가 난다고 할 수 있다. 그러므로 $F = \frac{Mean_{Between}}{Mean_{Within}}$ 을 (설명된 분산 / 설명되지 않은 분산) 으로 표기하는데, 설명된 분산이란 이미 계산할 값들이 정해지는 경우이고, 설명되지 않은 분산이란 계산 시 경우에 따라 Random 하게 다르게 정의하고, 자유도는 Between의 경우 k group-1, Within의 경우 n-k (sample-group)이 된다. 두 자유도의 합은 n-1이 되며, 이유는 k 그룹으로 mean을 구했으니 k-1 자유도가 되고, Within의 경우 전체 n 개에서 k개 그룹의 mean을 이미 구했으므로 n-k가 된다.

예를 들어서 남학생, 여학생, 외계인 학생의 발 길이에 대한 데이터가 있다고 하자.

boy = [27.3, 28.5, 29.7]
girl = [23.5, 24.2, 22.2, 25.7]
alien = [19.2, 18.6, 17.1]

이 데이터를 기반으로 각 집단의 평균에 대한 분산(Between), 각 집단 내 평균에 대한 분산(Within)을 구하면 다음과 같다.

\[F = \frac{⓵/dof_{between}}{⓶/dof_{within}} = \frac{156.66/2}{11.6/7} = 47.268\]

결국 집단별 표본평균의 분산은 중심극한정리에 의해 $\sigma_x = \frac{\sigma}{\sqrt{n}}$ 이므로, $\sigma^2 = n\sigma_x^2$이 되고, 표본분산을 이용해 추정한 모분산은 $\sigma^2 = ns_x^2$ 가 되므로, ⓵은 각각의 집단의 표본평균의 분산을 각각 구해서 더한 것 즉, $3 \times (28.5-23.6)^2 + 4 \times (23.9 - 23.6)^2 + 3 \times (18.3 - 23.6)^2$이고, 결국 집단 수에 대한 자유도로 평균을 낸 값이다. ⓶는 집단내 편차들의 전체 제곱합으로, 단순하게 각 집단 내의 모든 편차 제곱을 다 더해서 자유도로 평균낸 값이다.

5% 유의수준으로 검정해보면, 유의 수준 5% 일 때의 F값을 F테이블을 이용해서 구해보면 4.74이다.

결과적으로 p value가 매우 작으므로 귀무가설은 기각되어 3개의 그룹은 차이가 있다고 결론지을 수 있다.

파이썬을 이용하면 다음과 같다.

from scipy import stats
 
boy = [27.3, 28.5, 29.7]
girl = [23.5, 24.2, 22.2, 25.7]
alien = [19.2, 18.6, 17.1] 
 
stats.f_oneway(boy, girl, alien)

> F_onewayResult(statistic=47.26810344827576, pvalue=8.603395069970194e-05)

여기서 one way는 일원분산분석이라고 해서 one way AVONVA라 한다. 이 의미는 우리가 비교하려는 요인의 개수가 1개라는 의미이다. 그룹(집단)을 군 또는 수준 이라고도 하고, 요인은 각 그룹(집단)을 구분짓게 하는 것인데, 실험요인/인자/독립변수(Factor) 라고도 한다. 발길이에 대한 예를 다시 한번 살펴보면 group(boy, girl, alien)이 독립변수, 발길이가 종속변수. 즉, group에 따라 발길이가가 달라지는지 1개의 독립변수와 1개의 종속변수를 검정한 것이라고 할 수 있다.

다른 예를 들어 한 마케팅연구에서 30명의 표본을 선정한 후 이들을 임의로 세 가지 형태의 콜라 광고 중 하나를 시청하게 하였다. 1시간 동안 시청하게 한 후, 광고에 반응하는 것(구매욕구)을 측정하여 세 가지 광고의 차이가 있는지 알고자 할때,

측정값 : 구매욕구 (종속변수)
요인 : 광고 (독립변수)
수준 : 3가지 광고

가 된다. 광고에 따라 구매욕구가 달라지는가를 연구할 떄, one way ANOVA를 한다고 할 수 있다. 여기에 광고에 따라 그리고, 남/여에 따라 어떻게 되는지 확인하게 되면 two way ANOVA가 된다.

3개 이상의 집단끼리 차이를 분석할 때, 2개씩 짝이어서 t 검정을 하지 않고 ANOVA를 쓰는 이유는, 예를 들어 3개 그룹을 비교할 때는 3가지 짝지움을 할 수 있다. 이때, 1개의 짝지움에 대해 t Test를 하는 경우에는 α=5%로 보았을 때에는 잘못 판단하지 않을 확률이 95%이다. 이제 3번을 연속해서 검정하게 되는데, 유의하지 않을 확률 95%, 유의할 확률 5%의 동전을 던진다고 할 때, 3번 연속해서 동전을 던진 후에 한 번도 유의하지 않을 확률은 95% x 95% x 95%이다. 따라서 최소한 한 번이라도 유의할 확률은 1 - 95% x 95% x 95%가 된다. 이렇게 되면 원래 α는 5%로 시작했지만 3번 연속 검정을 하게 되면 1 - 95% x 95% x 95% = 0.142625가 되어 약 14.2%정도로 α가 커진다. 이렇게 되면 원래 5%로 판단하려고 했던 처음 의도와 다르게 14.2%를 기준으로 유의성을 판단하는 것과 같은 효과가 된다. 이를 일반화 하면, 전체 그룹에서 2개를 뽑는 경우만큼 t test를 하고, 신뢰도는 1-α이므로 $1 - (1 - \alpha)^{_{n}C _{2}}$ 가 된다. 이렇게 되었을 때 유의확률α가 커지니까 False Positive 오류 (Type I Error)를 범할 확률이 커지게 된다.

추가적으로 2개 그룹의 평균을 비교할 때 t test 대신 ANOVA를 써도 된다. 등분산의 2개 집단에 대해 Independent two Samples t-test결과와 ANOVA결과는 p value가 똑같이 나온다. 게다가 두 분석의 통계량은 제곱 관계로서, t값을 제곱하면 F값이 나온다.

ANOVA는 분산의 동질성(모분포에 대해)이 매우 중요하다. 이 분산의 동질성 가정을 만족하지 못하는 경우에는, Welch 검정을 하게된다.

Post Hoc

Post Hoc(사후분석)은 ANOVA에서 유의(Significant)하다는 검정 결과가 나왔을 때, 집단 중 어느 집단이 다른 것인가를 찾아내는 것이다. ANOVA의 Alternative Hypothesis는 최소한 1개의 집단은 평균이 다르다이므로 어느 것이 다른지를 찾아야 한다. 앞서 말했듯이 t test를 $_{n} C _{2}$번 해서 찾아내는 것은 $1 - (1 - \alpha)^{ _{n}C _{2}}$ 로 유의수준이 커지므로 Type I Error가 늘어난다.

사후분석을 위한 방법에는 Fisher’s LSD(피셔의 LSD), Tukey’ HSD 투키의 HSD), Bonferroni correction (봉페로니교정), Duncan (던칸의 방법), Scheffe (셰페의 방법), Games Howell (게임즈 하웰) 등이 있고, 이런 분석들을 통틀어 Post Hoc이라 한다. 이러한 사후분석들은 모두 등분산을 가정하고 차이는 다음과 같다.

1. Fisher’s LSD: 집단을 짝지어 t test를 아무 보정 없이 하는 방법으로 좋은 방법이 아니다.
2. Tuckey’ HSD: 비교 집단간 표본크기가 동일한 경우 사용하고, t분포를 사용한다.
3. Bonferroni: t검정을 짝지어서 검정한 후 유의수준 $\alpha$를 5%로 보정한다. 보정된 유의수준은 FWE라 부른다.
4. Duncan: Duncan은 작은 차이에도 차이가 난다고 결과를 낸다. 집단 간 차이가 꼭 드러나야 하는 경우 사용하는 것이 좋고, 사회과학 등에서 설문을 할 때 많이 사용한다.
5. Sheffe: 큰 차이가 나야 차이가 난다고 결과를 낸다. F분포를 활용해 검정한다.

위는 일반적인 ANOVA의 Post Hoc이고, 등분산이 아닌 경우 Welch 검정을 하는데 이 경우에는 Games Howell를 사용한다. Games Howell는 Welch의 방법으로 t분포의 자유도를 바꿔서 검정한다.

만약 ANOVA에서 유의하다고 결과가 나왔는데, Post Hoc에서 유의하지 않다고 나올 수도 있다. 즉, ANOVA에선 차이가 있다고 해놓고 막상 어떤 집단이 다른지 모르는 경우엔 집단 간 차이가 유의한 것으로만 단순 해석하는 것이 합리적이다.

앞서 발길이 예시로 실제 검정을 해보면 다음과 같다.

boy = [27.3, 28.5, 29.7]
girl = [23.5, 24.2, 22.2, 25.7]
alien = [19.2, 18.6, 17.1]

1.정규성 검정(Shapiro-Wilk)

from scipy.stats import shapiro
 
boy = [27.3, 28.5, 29.7]
girl = [23.5, 24.2, 22.2, 25.7]
alien = [19.2, 18.6, 17.1] 
 
print(shapiro(boy))
print(shapiro(girl))
print(shapiro(alien))
 
>
ShapiroResult(statistic=1.0, pvalue=0.9999986886978149)
ShapiroResult(statistic=0.9968306422233582, pvalue=0.9891670346260071)
ShapiroResult(statistic=0.9423076510429382, pvalue=0.5367357134819031)

세 집단 모두 정규성 검정을 통과한다. 사실 표본이 작기 때문에 큰 의미는 없고, 모집단이 정규성을 띄므로 표본의 평균도 정규적이라고 가정할 수 있다.

2.등분산성 검정(Levene, Bartlett)

from scipy.stats import levene, bartlett
 
print(levene(boy, girl, alien))
print(bartlett(boy, girl, alien))
 
>
LeveneResult(statistic=0.19816176470588293, pvalue=0.8246838520550716)
BartlettResult(statistic=0.19083867589274103, pvalue=0.9089916798322615)

Levene과 Bartlett에서 모두 등분산 검정을 통과하므로 3개의 집단이 모두 모분산이 동질하다고 할 수 있다.

3.One way ANOVA

from scipy.stats import f_oneway
 
f_oneway(boy, girl, alien)
 
>
F_onewayResult(statistic=47.26810344827576, pvalue=8.603395069970194e-05)

유의하다는 결과가 나오므로 세 집단에 차이가 있다.

이제 어느 집단이 다른지 확인할 수 있도록 Post Hoc을 해보면 표본의 크기가 서로 다르므로 Bonferroni를 사용한다. Bonferroni를 하기 전에 데이터의 모양새가 값과 그룹의 종류에 대한 데이터를 순서를 맞추어 넣으면 pandas dataframe로 다음과 같다.

import pandas as pd
 
lst_value = boy + girl + alien
lst_group = ["boy" for i in range(len(boy))] + ["girl" for i in range(len(girl))] + ["alien" for i in range(len(alien))]
df_total = pd.DataFrame({'Group':lst_group, 'Value':lst_value}).reset_index(drop=True)
print(df_total)
 
>
    Group  Value
0    boy   27.3
1    boy   28.5
2    boy   29.7
3   girl   23.5
4   girl   24.2
5   girl   22.2
6   girl   25.7
7  alien   19.2
8  alien   18.6
9  alien   17.1

이렇게 데이터를 만들어야 Bonfrroni를 할 수 있다.

from statsmodels.sandbox.stats.multicomp import MultiComparison
from scipy.stats import ttest_ind
 
comp = MultiComparison(df_total['Value'], df_total['Group'])
ret = comp.allpairtest(scipy.stats.ttest_ind, method='bonf')
print(ret[0])
 
>
Test Multiple Comparison ttest_ind 
FWER=0.05 method=bonf
==============================================
group1 group2   stat    pval  pval_corr reject
----------------------------------------------
 alien    boy -10.9355 0.0004    0.0012   True
 alien   girl  -5.5521 0.0026    0.0078   True
   boy   girl   4.4257 0.0069    0.0206   True
----------------------------------------------

pval_corr이 집단의 개수로 곱해진 보정된 p value이고, 모두 유의하다는 것을 알 수 있다. 따라서 서로 모두 차이가 있다.

RMANOVA

특정 시점을 기준으로 그 시점의 전후 2개 집단의 차이를 검정할 때 Paired t Test를 했었다. Paired t Test는 변화 한 번만 봐서 집단이 2개인 경우인데, RMANOVA(Repeated Measured ANOVA)는 변화를 여러번 보는 것이다. 따라서 변화에 대해 비교 시점이 여러 개가 되고, 비교하는 집단이 여러 개가 된다.

Paired t Test때와 비슷한 예시로 어떤 정신과 의사가 환자들을 상대로 스트레스 호르몬인 코르티솔을 줄이기 위해 코끼리를 타고 바흐의 음악을 계속 듣게 했다고 하자. 현재 코르티솔을 측정하고, 1주일 후에 측정하고, 2주일 후에 측정했다고 하면 세 시점이 된다.

산수적인 계산은 ANOVA와 동일하므로 파이썬으로 분석을 해보자. Null Hypothesis는 “세 시점에 변화가 없다.”이고, Alternative Hypothesis는 “적어도 한 시점에서는 변화가 있다.”이다. 파이썬으로 RMANOVA를 하기 위해서 statsmodel의 AnovaRM 모듈을 이용한다.

AnovaRM를 사용하려면 각 시점에 대한 데이터를 한 개의 컬럼에 모두 넣은 후 Subject와 시점을 붙여 넣어야 한다.

from statsmodels.stats.anova import AnovaRM
 
print(AnovaRM(data=df_data, depvar='value', subject='id', within=['viewPoint']).fit())
 
>
 Anova
                Anova
=======================================
          F Value Num DF  Den DF Pr > F
---------------------------------------
viewPoint 58.8882 2.0000 38.0000 0.0000

이런 식으로 AnovaRM 분석을 할 수 있고, F검정 결과까지 한 번에 볼 수 있다. F Value가 58.8882로 매우 큰 통계량이 나왔고 p value(Pr>F)는 매우 작다. 따라서 F 검정의 Null Hypothesis를 기각하여 3개의 시점 중에 차이가 나는 다른 것이 적어도 하나 있다고 할 수 있다.

이제 Null Hypothesis가 기각되었으니 어떤 시점이 다른 시점과 다른지 알아야 한다. RMANOVA를 할 때는 전제조건이 있는데, 정규성 전제와 구형성전제가 있다. 정규성 전제는 기존에 했던 것과 동일하고 구형성(Sphericity)은 다음과 같다.

ANOVA의 등분산성에 대응되어서 RMANOVA는 차이의 분산들이 등분산성을 가져야 한다는 의미이다. 이를 pingouin 패키지를 통해 테스트할 수 있고, Mauchy’s Test of Sphericity라 한다.

import pingouin as pg
pg.sphericity(data=df_data, dv='value', within='viewPoint', subject='id')
 
>
SpherResults(spher=False, W=0.6111256376806817, chi2=8.864148863341935, dof=2, pval=0.011889799495965416)

결과를 확인해보니 구형성을 만족하지 못한다. 이럴 때 t test에서 등분산이 아닐 때 Welch’s Testing을 위해서 correction을 줄 수 있던 것과 같이 할 수 있다.

import pingouin as pg
pg.rm_anova(dv='value', within='viewPoint', subject='id', data=df_data, correction=True)
 
> 
Source    ddof1 ddof2  F            p-unc	...
viewPoint  2     38    58.888199    2.282217e-12 ...

이런 경우에 correction을 True로 주면 구형성을 만족하지 못하는 경우에도, Greenhouse-Geisser correction과 Huynh-Feldt correction을 적용하여 결과를 낸다. correction을 한 결과의 F값은 58.888199이고, p value는 2.282217e-12로 귀무가설이 기각될 수 있다. correction=False로 주면 F통계량이 58.8882, p value는 거의 0. 으로 아까의 결과와 동일한 결과를 낸다.

이제 Post Hoc을 할 차례다. RMANOVA도 Post Hoc이 있는데, Benjamini/Hochberg FDR correction이라는 방법이다. 각 시점별로 pair를 만들어서 t Testing을 할 때 FWER를 보정해서 결과를 알려준다. 아이디어는 Bonferroni와 같고 이 방법은 pingouin에서 제공한다.

import pingouin as pg
print(pingouin.__version__)
posthoc = pg.pairwise_ttests(dv='value', within='viewPoint', subject='id', data=df_data)
print(posthoc)
 
>
A             B            Paired T         dof alternative p-unc
1 week later  2weeks later  TRUE  4.152317  19  two-sided   5.41E-04
1 week later  present TRUE  TRUE  -6.997791 19  two-sided   1.15E-06
2 weeks later present TRUE  TRUE  -9.38871  19  two-sided   1.44E-08

결과를 확인하면 다 다른다는 것을 알 수 있다. 결국 코르티솔은 계속 줄어서 코끼리를 타고 바흐를 듣는 것은 개선의 효과가 있었다고 할 수 있다.

Cross tabulation

이번엔 교차분석에 대해 알아볼 것이다. cross table(contingency table)은 빈도에 관련한 교차표이다. 범주형 자료에 대해서 빈도를 분석한다고 생각할 수 있다. 다음과 같이 교차집계표를 만들 수 있다면, 이미 검정에서 살펴본 독립성(연관성), 동질성 검정이 가능하다.

이런 식의 데이터가 있다고 할 때, cross table은 두 개의 컬럼을 교차하여 만든다.

cross table은 Grouping and Agggregation이다. row와 column을 Grouping을 한 후에 그에 대한 집계를 하는 것이다. 따라서 사과, 수박, 오렌지가 아무리 많이 나와도 이 세 가지로 Grouping이 가능하고, 출하일은 9월 1일, 9월 2일, 9월 3일로 모두 Grouping 해서 몰아넣을 수 있다.

cross tabulation은 pandas의 crosstab을 이용해 만들 수 있다.

pd.crosstab(df_data['과일'], df_data['출하일'])

이 결과는 과일의 unique value들 (Grouping)과 출하일의 unique value만(Grouping)을 index와 column으로 만들어서 그 수를 count(Aggregation)하는 결과이다.

raw data가 pandas의 Dataframe일 때, pandas의 pivot으로도 만들 수 있다.

pd.pivot_table(df_data, index='과일', columns='출하일', aggfunc='size')

이렇게 하면 똑같은 결과가 나온다. pivot은 pandas dataframe을 통째로 데이터로 넣은 다음에 그중에 index와 column이 무엇인지 정해주면 된다. pivot을 이용해서 집계를 할 수 있는데, 과일별 - 출하일에 대해서 당도의 (각각의 Grouping) 평균(Aggregation)을 보고 싶다고 한다면 데이터 개수를 count를 하는 대신에, count자리에 다른 column값을 넣을 수 있다.

pd.pivot_table(df_data, index='과일', columns='출하일', values='당도')

$\chi^2$ Testing, Cramer V

교차분석은 관심 있는 표본으로부터 Categorical(명목형) 데이터에 대하여 pivot table이나 cross table (contingency table)을 만들어 행과 열을 교차하여 교차빈도집계표를 구한 후에 카이스퀘어 검정을 이용하여 동질성, 독립성을 검정하는 것이다.

예시를 통해 분석, 검정, 연관도 확인까지 해보자. 아래와 같인 데이터가 있다고 하자.

최종적으로 카이스퀘어 검정을 하니까 가설을 설정하면 다음과 같다.

$H_0$: 두 변수는 독립이다.
$H_1$: 두 변수는 독립이 아니다. (연관이 있다.)

일단 첫 번째로 거주지에 따른 대학 진학에 차이가 있는지 확인해보자.

df_cross_table = pd.crosstab(df_data['거주지'], df_data['대학진학'])
 
df_cross_table = pd.pivot_table(df_data, index='거주지', columns='대학진학', aggfunc="size")
 
>
대학진학  미진학   진학
거주지           
광역시    90  122
시구군    42   50
시군구    16   30
특별시   110  150

이렇게 만들어진 교차집계표로 바로 카이스퀘어 검정을 할 수 있다. 교차표의 각 row가 각각 list로 들어가야 하니까 각 row에 대한 list를 만들어 [[90, 122],[42, 50], [16, 30], [110, 150]]의 형태가 되도록 넣어주면 된다.

from scipy.stats import chi2_contingency
rows = [row.to_list() for i, row in df_cross_table.iterrows()]
chi2_contingency(rows, correction=False)
 
>
(1.490709512981629,
 0.6844162240071692,
 3,
 array([[ 89.66557377, 122.33442623],
        [ 38.91147541,  53.08852459],
        [ 19.4557377 ,  26.5442623 ],
        [109.96721311, 150.03278689]]))

p value가 0.68이므로 거주지에 따른 대학진학결과가 독립이라 할 수 있다. 결과로 나오는 첫 번째는 카이제곱 통계량, 두 번째는 p value, 세 번째는 자유도, 마지막은 서로 독립일 때의 expected 빈도의 row list이다.

이번에는 부모의 학력이 대학 진학에 영향을 끼치는지 확인해보자.

df_cross_table = pd.crosstab(df_data['부모학력'], df_data['대학진학'])
>
대학진학  미진학   진학
부모학력          
고졸    112  126
대졸     70  140
대학원졸   76   86
 
 
rows = [row.to_list() for i, row in df_cross_table.iterrows()]
chi2_contingency(rows, correction=False)
>
(10.539167067134844,
 0.005145753160303171,
 2,
 array([[100.66229508, 137.33770492],
        [ 88.81967213, 121.18032787],
        [ 68.51803279,  93.48196721]]))

부모의 학력과 진학은 독립이 아님을 확인할 수 있다. 서로 연관이 있으므로 얼마나 연관이 있는지 알아봐야 하는데, 명목형 변수-명목형 변수(Categorical - Categorical)를 분석하는 경우에는 크래머V (Cramer V)라는 것을 이용해서 얼마나 서로 연관되어 있는지 알 수 있다. 명목형 변수는 대상에 대해 측정하면 대상을 일정한 범주에 속하게 하며 대상에 이름이 붙여지지만 각 범주 간에 순위는 없는 항목을 말한다. 예를 들어 성별(남, 여), 인종(황인종, 흑인종, 백인종), 혈액형(AB, A, B, O), 검사결과여부(양성, 음성) 등이 있다.

\[CramerV = \sqrt{\frac{\chi_{stat}^2 \div n}{min(r-1, c-1)}}\]

위 식을 Cramer V - 연관 계수라 한다. 카이스퀘어의 결과값은 고정된 같은 확률에 대해 전체 표본의 개수가 크면 큰 값이 나오는 경향이 있기 때문에 이를 없애기 위해 $n$으로 나눈다. 그리고 이 나눈 값의 최댓값이 $min(r-1, c-1)$이 되는데 이 값으로 다시 나누어 0~1사이의 값으로 정규화 한다. 이때 이 최댓값을 Cramer의 phi $\phi$라 한다.

쉽게 말해 카이스퀘어값과 관련되어 있으니까 서로 독립일 때 기대했던 값과 차이가 많이 날 수록 더 연관되어 있다고 볼 수 있다. 결과를 계산해보면,

import numpy as np
 
x2 = chi2_contingency(rows, correction=False)[0]
n = np.sum(rows)
minDimension = min(np.array(rows).shape)-1
 
V = np.sqrt((x2/n) / minDimension)
 
print(V)
 
>
0.13144323132393237

0.131 정도로 서로 영향을 끼치긴 하는데 그렇게 큰 영향은 아니다. 이 정도는 통계적으로 유의하고, 약한 연관관계를 갖는다고 할 수 있다. 보통 0.6이상이 되면 강한 연관관계가 있다고 한다.

Correlation Analysis

Covariance, Correlation Coefficient, Spearman

교차분석에서 Categorical - Categorical 데이터에 대한 연관도를 계산해 보았다면 이번엔 연속형-연속형 변수들의 관계를 살펴볼 것이다. 먼저 공분산에 대해 알아야 하는데, 1개 변수의 흩어짐에 대한 정도를 나타내는 것이 분산이라면, 공분산은 2개의 확률변수의 차이에 대한 분산이 아닌 같이 변화하는 양을 나타내는 상관 정도를 나타내는 값이다. 따라서 Scattering 관점에서 분산은 값이 커질수록 더 흩어지는데, 공분산은 값이 커질수록 더 흩어지는 지표가 아니다.

표본 데이터에 대해 공분산은 다음과 같이 구해진다.

\[S_{XY} = \frac{\sum_{i=1}^n(x_i-\bar{x})(y_i-\bar{y})}{(n-1)}\]

비교를 위해 분산의 공식을 보면,

\[s^2 = \frac{\sum_{i=1}^n(x_i - \bar{x})^2}{n-1} = \frac{\sum_{i=1}^n(x_i - \bar{x})(x_i - \bar{x})}{n-1}\]

분산은 자기 자신과의 관계라 볼 수 있는데, 공분산 Covariance은 각각의 평균으로부터 얼마나 같이 흩어져 있는지를 본다.

상관관계는 (x, y) 쌍의 데이터의 좌표를 나타내는 점들이 오른쪽 위 또는 오른쪽 아래를 향하는 경향성을 보인다면 양의 또는 음의 상관관계가 있다고 표현하는데, 오른쪽 위를 향한 경향성을 띌때 공분산을 계산해보면 Positive가 나와서 양의 상관, 오른쪽 아래를 향하는 경향성이 보이면 공분산 계산 결과가 Negative가 나와서 음의 상관으로 볼 수 있다. 양의 상관이 나온다면 같이 증가하는 형태, 음의 상관이 나온다면 하나가 증가할 때 나머지 하나가 감소하는 형태이다.

x, y 각각의 평균으로부터 차이의 부호가 $(x-\bar{x})(y-\bar{y})$의 부호이고, 4분면으로 표시해서 따지면, 공분산이 Positive인 경우에는 +를 갖는 값이 -를 갖는 값보다 많으므로 합이 양수이고, 음수일 때는 반대이다. 만약 +와 -값이 비등하게 있다면 다 더했을 때 0에 가까운 값을 가질 테고 그럴 때는 서로 상관도가 작다고 볼 수 있다. 따라서 Positive이든, Negative이든 공분산이 크다면 선형적 연관성이 높다고 할 수 있으며 공분산의 범위는 $-\infty \sim \infty$ 이다.

공분산은 값의 크기 보다는 Sign(부호)에 의해서 선형 관계의 방향을 나타내기 때문에 문제가 있는데, x, y의 단위에 영향을 많이 받는다는 것이다. 상관관계 강도와 상관없이 단위가 크면 큰 값이 나오게 된다. 예를 들면 1000점 만점의 영어성적과 수학성적의 공분산과 50점 만점의 영어성적과 수학 성적을 비교하게 되면 1000점 만점의 영어성적과 수학 성적의 상관관계가 더 적더라도 1000점 만점의 공분산이 더 큰 값을 가질 수 있다.

이를 해결하기 위해서 정규화를 한다. Correlation Coefficient통해 비교가 가능하도록 하는데, 데이터 정규화처럼 각각의 표준편차로 정규화하고, 이렇게 하면 $-1 \sim +1$ 사이의 값을 갖게 된다.

\[r = \frac{S_{XY}}{S_X \cdot S_Y} = \frac{\frac{\sum_{i=1}^n(x_i-\bar{x})(y_i-\bar{y})}{(n-1)}}{\sqrt{\frac{\sum_{i=1}^n(x_i - \bar{x})^2}{n-1}} \cdot \sqrt{\frac{\sum_{i=1}^n(y_i - \bar{y})^2}{n-1}}}\]

상관계수의 수식을 보면 공분산을 각각의 분산으로 나눈 값이 된다. 이 값이 자기 자신과의 상관이면 최대값이 되는데 그 값이 1이다.

상관계수도 검정이 가능하다. 그렇다면 상관계수가 어떤 분포를 따를지를 알아야 Null Hypothesis를 세울 수 있는데, 여기서 모집단의 상관계수를 $\rho$로 표본의 상관계수를 $r$로 표현한다.

자세한 유도는 너무 복잡하고 결론적으로 $r$분포는 $t$분포를 따르고, 통계량은 다음과 같다.

\[t_{stat} = \frac{r -\rho_0}{\frac{\sqrt{1-r^2}}{\sqrt{n-2}}}\]

주어진 데이터(표본)에서 r을 계산하고 나면 모상관계수 $\rho_0$를 안다는 가정하에 위와 같이 t분포의 검정 통계량을 구할 수 있다.

추가적으로 표본평균분포의 자유도 n-1과 상관계수분포의 자유도 n-2에 차이가 있는데, 표본평균을 따질 때는 평균을 구할 때 데이터를 쓰면서 $s$를 구한다. 따라서 $t(n-1)$을 따르고, 상관계수는 $r$을 구할 때 분모의 분산을 구하는 과정에서 이미 $X$의 평균과 $Y$의 평균을 구하기 위해 데이터를 하나씩 쓰므로 최종 자유도가 n-2가 되어 $t(n-2)$를 따른다.

1 Sample t Test를 했을 때의 분포 그림을 비교해보면 다음과 같다.

이제 Null Hypothesis와 Alternative Hypothesis를 설정하기 위해서 $\rho_0$를 특정해야 한다. $\rho_0=0$로 특정하면 가설이 다음과 같다.

$H_0$: 모상관계수는 0이다. $\rho_0 = 0$
$H_1$: 모상관계수는 0이 아니다. $\rho_0 \neq 0$

$\rho_0$는 0이 아닌 어떤 특정값으로도 t testing이 가능하다. 예를 들면 $\rho_0 = 0.5$도 가능하다.

이제 검정을 위해 실제 예를 들어 데이터가 다음과 같이 있다고 하자.

import numpy as np
import pandas as pd
 
x_data = [1, 2, 4, 8, 10, 12, 13, 8, 17, 20]
y_data =  [2, 3, 7, 7, 9, 11, 12, 11, 15, 17]

Covariance공분산은 numpy에서 직접 계산할 수 있고, Pearson Correlation Coefficient상관계수는 pandas의 Dataframe에서 계산할 수 있다.

import math
 
### 공분산 ; numpy로 곧바로 계산 가능하다.
dof = len(x_data)-1  # 자유도
np.cov(x_data, y_data, ddof=dof) # covariace
>
array([[348.5, 260. ],
       [260. , 208.4]])
 
### 상관계수 ; pandas로 곧바로 계산 가능하다. 
df_corr = pd.DataFrame({'x_data':x_data, 'y_data':y_data})
>
x_data  y_data
0       1       2
1       2       3
2       4       7
3       8       7
4      10       9
....
 
df_corr.corr() # 상관계수를 계산!
>
          x_data    y_data
x_data  1.000000  0.964768
y_data  0.964768  1.000000

두 결과 모두 x_data와 y_data를 행과 열로 순서대로 둔 행렬 형태로 나온다. x와 y사이의 공분산은 260이다. x_data의 분산은 348.5, y의 분산은 208.4이다. Correlation Coefficient는 0.964768이다. 당연히 자기 자신과의 Correlation은 1.000000이 나온다.

Correlation Coefficient의 t 검정은 scipy의 pearsonr()를 통해 할 수 있다.

from scipy.stats import pearsonr
pearsonr(x_data, y_data)
>
(0.9647684650479604, 6.459870025582558e-06)

Correlation Coefficient는 0.9647684650479604, p value는 6.459870025582558e-06로 상관도가 매우크다. t 검정 결과 역시 유의하다. 따라서 $\rho_0 = 0$라는 귀무가설을 기각할 수 있으며, 모상관계수는 0이 아니라 주장할 수 있다. 즉, x와 y는 상관이 있다고 할 수 있다.

상관관계는 두 개의 변수끼리 선형관계가 있는지에 대한 조사 방법이다. 서로 독립이라면 선형관계가 없다고 할 수 있다. 즉, 독립 $\to$ 공분산 = 0 이 된다. 하지만 반대로 공분산이 0이라 해서 항상 독립인건 아니다. 이유는 공분산은 선형 관계에 대한 의존성을 측정하고, 공분산이 0이 되면 두 변수의 선형 관계가 없다는 것을 의미하지만 독립이라고 볼 수는 없다.

앞서 공분산 식을 볼 때 공분산의 자유도가 n-2가 아닌 이유는 x와 y의 평균을 구해서 n-2가 되는 것이 아니라 x, y로 이루어진 2차원의 $(x_i, y_i)$ 데이터 중에서 쌍으로 된 $(\bar{x}, \bar{y})$로 이루어진 표본평균을 구했기 때문이다. 따라서 $(x, y)$ 한 쌍을 사용했으므로 n-1이 된다. 상관계수에 대한 통계량의 자유도가 n-2가 나오는 이유는 Normalization을 위해서 분모에 나눌 때 x와 y의 분산을 따로 사용했기 때문에 각각 데이터가 하나씩 나가서 n-2가 된다.

상관계수를 해석핼 때 주의할 점이 있는데, 상관계수가 크다고 해서 데이터가 더 조밀하게 모여 있다고 할 수 없다. 그저 4사분면 중 어디에 더 (x, y) 쌍이 많은지만 나타낼 뿐이다. 즉, 더 밀집도보다는 선형관계인가 정도만 가늠할 수 있다. 또한 상관관계를 볼 때 개별 데이터가 아닌 그룹별 평균 데이터를 사용하여 상관관계에 대한 결론을 내면 각각의 데이터에 대한 상관관계가 과장될 수 있다. 아래 그림과 같이 개별 데이터에 대한 상관관계는 실제로 훨씬 작을 수 있다.

마지막으로 상관관계를 볼 때 Outlier가 있거나 비선형관계일 때는 상관계수가 의미를 퇴색한다. 비선형관계일 때는 의미를 찾기 어렵고 Outlier는 큰 cov값을 갖는 점이 있을 테니 이것이 상관관계의 왜곡을 만들어 낸다. 만약 Outlier가 있는 경우에 Robust 하게 상관관계를 측정할 수 있는 방법이 있는데 그것이 Spearman 상관계수이다. Spearman 상관계수를 계산할 때는 숫자 데이터를 그냥 다 쓰는 것이 아니라 숫자 데이터에 순위를 다시 매긴 후에 사용하기 때문에 갑자기 큰 숫자가 나오더라도 문제가 되지 않는다.

Spearman 상관계수는 Pearson 상관계수로부터 시작하여 $d_i = x_i - y_i$로 정의하여 $d$에 대해 정리하면 $r = 1 - \frac{6\sum d_i^2}{n(n^2-1)}$이 된다. d는 x의 순위와 y 순위의 차를 의미한다.

Correlation Analysis

지금까지 상관계수에 대해 다뤘으니 이제 실제 상관분석에 대해 알아보자. 우선 상관계수의 값에 대해서 해당 값이 어느 정도의 상관이 있는지는 다음과 같다.

상관계수 범위	상관의 정도
$\pm 0.2$ 미만	거의 상관 관계 없음
$\pm 0.2 \sim \pm 0.4$	조금 상관 있음
$\pm 0.4 \sim \pm 0.7$	상관 있음
$\pm 0.7 \sim \pm 0.9$	분명히 상관 관계가 있음
$\pm 0.9$ 이상	매우 높은 상관 관계가 있음
$\pm 1.0$	자기 자신

실무적으로는 0.4 이상만 되도 어느 정도 관계가 있다 할 수 있고, 0.9이상이면 사실상 같은 변수라 판단해도 무리가 없다.

이제 예를 들어 다음과 같은 데이터가 있다고 하자. 카페의 매출에 상관관계가 있는 것이 무엇인지 알아내기 위해 고객의 평가 데이터를 모았다고 한다.

이 데이터들에 대해 Correlation Analysis를 하면 다음과 같다.

import padas as pd
df_data.corr() # df_data는 pandas dataframe임

대각 성분은 자기 자신과의 성분이기 때문에 1이 나온다. 나머지 값들은 다른 요소들 간의 상관관계고 당연히 대각 성부을 기준으로 위와 아래는 대칭이다. 값을 보면 1인당 매출에 커피맛과 알바외모가 큰 양의 상관관계가 있음을 알 수 있다.

이러한 상관관계 값에 대한 검정결과 p value도 확인할 수 있는데 다음과 같다.

df_data.corr(method=lambda x, y : pearsonr(x,y)[1])

p value table을 보면 커피맛-1인당매출, 알바외모-1인당매출의 값이 Significant함을 알 수 있다.

이렇게 확인한 상관분석은 보통 회귀분석을 할 때, 회귀를 하기 전 상관분석을 한다. 상관분석 후 얻을 수 있는 시사점은 다음 2가지가 있다.

1. 독립변수들끼리 너무 서로 상관관계가 큰 경우에 어떤 현상이 일어나는데, 그 어떤 현상이라는 것이 서로 상관관계가 큰 변수끼리 다중공선성을 일으킬 가능성이 매우 높다. 다중공선성이란 독립변수끼리 같이 움직이긴 하지만 당연히 비슷하게 같이 움직이는 독립변수 중에 대표적인 1개 독립변수 하나만 있으면 대충 어떤 식으로 회귀 모델에 그 변수가 기여하는지 알 수 있는데, 이것들이 모두 모델에 포함되면 독립변수 하나만 있을 때 보다 그 두 개 변수의 합 때문에 그 변수의 합의 분산이 1개 있을 때 보다 커져버려서 회귀 모델의 계수의 분산을 크게 할 수 있다는 뜻이다.
2. 독립변수들이 종속변수와 너무 관계가 없어도 어려워서 어느정도는 관계가 있어줘야 회귀결과에도 반영될 수 있다.

여기에 추가적으로 상관계수는 큰데 검정을 했을 때 유의하지 않다면 보통 표본의 크기가 작아서 그런 것일 가능성이 있다.

정리해서 상관분석을 통해 할 수 있는 것은 Feature Extraction, Feature Selection이다. Feature Selection은 적당한 Feature를 고르는 것이고, Feature Extraction은 다르게 말하면 차원 축소이다. Feature Selection은 데이터중에서 우리가 관심 있는 종속변수 y가 있는 경우 y와 상관이 적당히 있는 것들만 골라서 관련이 적은 컬럼을 삭제하여 독립변수의 수를 줄여 차원을 축소한다는 것이고, Feature Extraction의 경우 종속변수 y가 없는 데이터라면 PCA 주성분 분석을 이용해 어떤 것이 이 데이터중에서 주 성분인지를 알아내어 독립변수를 합치는 방법으로 차원을 축소하기도 한다. 여기에서 Feature는 데이터의 컬럼을 의미한다.

다중공선성: https://blog.naver.com/vnf3751/220833952857