RMANOVA (Repeated Measured ANOVA, 반복측정 ANOVA)의 우아함

이전에 특정 시점을 기준으로 그 시점의 전후로 2개 집단의 차이를 검정했던 것이 있었는데, 그것이 무엇이냐면 Paired t Test입니다. 변화가 있었는지 보는 검정 방법이고, 왠지 있어 보이는 분석 방법이라서 애정이 가는 분석/검정입니다. 이 분석의 Application으로는 보통 주어진 데이터에서 패턴을 찾아내서 분석하여 의미를 찾기보다는 의학 쪽에서 제한된 환자들을 대상으로 시간별 변화 추이들을 볼 때 많이 사용합니다.

Paired t Test는 변화 한번만 보니까 집단이 2개인 경우인데요, Repeated Measured ANOVA는 변화를 여러 번 보는 것입니다. 그러니까 변화에 대해 비교하는 시점이 여러 개가 되고, 그러니까 비교하는 집단이 여러 개가 되는 거죠. 그러니까 ANOVA입니다. 반복 측정이라는 말이 붙었을 뿐이고요.

확인 시점을 2회로 전후 비교는 paired t Test
확인 시점을 3회 이상으로 하여 매 순간의 비교는 RMANOVA
라고 정리하면 좋겠습니다.

paired t Test에서 가장 중요하게 생각했던 것을 상기해 보면 Paired니까 한 개의 개체에 대해서 전/후의 데이터가 쌍이 꼭 맞아야 했었듯이, RMANOVA도 당연히 모든 시점에서 각 개체에 대한 데이터가 있어야 하겠죠. 모든 시점에서 데이터가 마련되어 있지 않은 개체는 아까비~라도 버려야 합니다.

어쨌든 paired t Test에서 들었던 예를 확장해서 접근해 볼까 합니다. 어떤 정신과 의사가 환자들을 상대로 스트레스 호르몬인 코르티솔을 줄이기 위해 코끼리를 타고 바흐의 음악을 계속 듣게 했다고 합니다. 이런 치료 기법을 아라불라샤라고 하면 좋겠군요. - 이런 치료 기법은 실제로는 없습니다. 주의 -
현재 코르티솔을 측정하고, 1주일 후에 측정하고, 2주일 후에 측정했다고 한다면 말이죠. 즉 세 시점이 되겠군요.

ANOVA에서 F검정에 대한 산수적인 접근은 충분히 했다고 생각하고, 이번에는 직접 파이썬의 도움을 얻어서 곧바로 분석하는 법을 따라가 볼까 합니다. 이제까지 본 통계 분석 방법을 한번 해보는 거죠.

일단 세 시점을 세 개의 집단으로 보았다면, 이 세 집단이 과연 다를까 미리보면 좋으니까 평균을 시각화해서 확인해 본다면,

시점을 구분하여 현재 - 1주뒤 - 2주 뒤이고, 데이터 분포를 보면 차이가 많이 나서 딱! 봐도 평균들이 차이가 나 보이는군요? 최대 최솟값도 보면 표본 분산이 많이 차이가 나는 것 같네요. 그러면 일단 차이가 나서 다르겠구먼 음.. 정도의 결론을 마음속에 품고 있으면 결과가 나왔을 때 이게 뭐지? 라며 놀라지 않겠습니다.

이제와서는 눈을 감고 잠시 생각해 보면, 당연한 이야기이겠지만, - 그래야 차이를 0으로 둔 그렇다~치고 분포를 가정할 수 있겠죠! -
Null Hypothesis는 세 시점에서 변화가 없다. 이고,
Alternative Hypothesis는 적어도 한 시점에서는 변화가 있다고 할 수 있겠습니다.

그러면, 준비도 되었겠다 시작해 보시죠. 파이썬으로 RMANOVA를 하는 방법은! statsmodel의 AnovaRM 모듈을 이용하면 아주 간단하게 할 수 있습니다.

사용법은 AnovaRM(data=판다스 데이터프레임, depvar=종속변수, subject=개체구분자, within=[시점]).fit()) 입니다.

간단하쥬? 그렇긴 한데, 데이터를 우리가 보기 편한 형식으로는 AnovaRM을 사용할 수 없는 것이 조금 불편하긴 한데요, 어떤 식이냐면, 각 시점에 대한 데이터를 한 개의 컬럼에 모두 넣은 후 Subject와 시점을 붙여 넣어서 데이터를 만들어야 합니다.

이렇게 하고 나면, df_data가 데이터를 담고 있는 dataframe가 된다고 하면,

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)는 엄청 작은 값이 나오겠죠. 그러니까 0이나 다름없으니까 0으로 표시됩니다 - 얼마나~ 작으면 0이겠어요. - . 결론적으로 F 검정의 Null Hypothesis를 기각하게 되는데 그것이 3개의 시점 중에 차이가 나는 다른 것이 적어도 하나 있다고 결론 맺을 수 있겠습니다. 어허.

간단하게 Null Hypothesis가 기각되지 않았다면, 쓱싹 어떻게든 넘어갈 요량이었는데, 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으로 엄청나게 작은 결과를 내뱉는군요. 그러니까, 어쨌든 Null Hypothesis가 기각될 수 있겠습니다.

혹시나 모듈이 잘 동작하는지 확인차 correction=False로 주면?

import pingouin as pg
pg.rm_anova(dv='value', within='viewPoint', subject='id', data=df_data, correction=False)

>
Source    ddof1 ddof2   F           p-unc ...
viewPoint 2     38      58.888199   2.282217e-12 ...

F통계량이 58.8882, p value는 거의 0. 아까의 결과와 같군요. 휴. 일단 Null Hypothesis를 기각하고 나면, 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

아이고 다 다르네요. 결국 코르티솔은 계속 줄어서 코끼리를 타고 바흐를 듣는 것은 개선의 효과가 있었다고 볼 수 있겠습니다. 후.