그런데, 여기에는 각각의 분포가 정규성을 따른다는 중요한 가정이 필요하다. 하지만 정규성을 따르지 않는다면 어떻게 해야 할까? 그럴 때 사용하는 것이 크루스칼-왈리스 검정(Kruskal-Wallis Test)이다. 이번 포스팅에서는 크루스칼 왈리스 검정과 그 사후 분석 (post hoc analysis)까지 알아보고자 한다.
PROC UNIVARIATE DATA=hong.df NORMAL PLOT;
CLASS SMOK;
VAR NICOT_ADDICT;
HISTOGRAM NICOT_ADDICT/ NORMAL (MU=EST SIGMA=EST);
RUN;
정규성 - 결과
결과 - 1) 비흡연자 (SMOK=0)
표본의 수가 2,000이 되지 않으므로 Shapiro-Wilk 통계량의 p-value를 보아야 하고, 이는 0.05보다 작다. Q-Q plot에서는 직선에서 벗어나는 표본들이 많고, 히스토그램에서도 right-skewed 분포를 보이므로 정규성을 따르지 않다고 결론을 짓는다.
결과 - 2) 과거 흡연자 (SMOK=1)
표본의 수가 2,000이 되지 않으므로 Shapiro-Wilk 통계량의 p-value를 보아야 하고, 이는 0.05보다 작다. Q-Q plot에서는 직선에서 벗어나는 표본들이 많고, 히스토그램에서도 right-skewed 분포를 보이므로 정규성을 따르지 않다고 결론을 짓는다.
결과 - 3) 현재 흡연자 (SMOK=2)
표본의 수가 2,000이 되지 않으므로 Shapiro-Wilk 통계량의 p-value를 보아야 하고, 이는 0.05보다 작다. Q-Q plot에서는 직선에서 벗어나는 표본들이 많고, 히스토그램에서도 uniform distribution 분포를 보이므로 정규성을 따르지 않다고 결론을 짓는다.
따라서 일원 배치 분산 분석 (ANOVA)을 시행할 수 없고, 크루스칼 왈리스 분석 (Kruskal-Wallis Test)을 시행해야 한다.
지금까지 많은 분석을 할 때 분산이 같은지 여부 (등분산성)을 중요하게 여겼었다. 하지만 크루스칼 왈리스 분석은 분산이 달라도 어느 정도 괜찮다. 물론 분산에 아주 큰 차이가 나면 위험해지지만, 어느 정도는 robust 하므로 웬만하면 그냥 써도 괜찮다.
크루스칼 왈리스 (Kruskal Wallis) 검정 코드
PROC NPAR1WAY DATA=hong.df WILCOXON CORRECT=NO;
CLASS SMOK;
VAR NICOT_ADDICT;
RUN;
P-value가 <0.0001이므로 귀무가설을 기각하고 대립 가설을 채택한다. 그렇다면 여기에서 귀무가설 및 대립 가설은 무엇이었는가?
귀무가설: $H_0=$ 세 집단의 모평균은 "모두" 동일하다.
대립 가설: $H_1=$ 세 집단의 모평균이 모두 동일한 것은 아니다.
우리는 대립 가설을 채택해야 하므로 "세 집단의 모평균이 모두 동일한 것은 아니다."라고 결론 내릴 것이다.
이 말은, 세 집단 중 두 개씩 골라 비교했을 때, 적어도 한 쌍에서는 차이가 난다는 것이다. 따라서 세 집단 중 두 개씩 골라 비교를 해보아야 하며, 이를 사후 분석 (post hoc analysis)이라고 한다. 세 집단에서 두 개씩 고르므로 가능한 경우의 수는 $_3 C_2 =3$이다.
(1) 비흡연자 vs 과거 흡연자
(2) 비흡연자 vs 현재 흡연자
(3) 과거 흡연자 vs 현재 흡연자
비모수 검정인 Kruskal Wallis 검정의 사후 분석으로 쓰이는 방법은 다음 네 가지가 있다.
1) Pairwise Wilcoxon Test
- 가장 적절하지 않다. 왜냐하면 Kruskal wallis 검정은 기본적으로 순위를 이용한 검정 방법인데, 전체 집단에서 특정 피험자의 순위는 부분 집단에서의 순위와 다르기 때문이다. 예를 들어 전체 집단에서 흡연자인 어떤 피험자의 니코틴 중독 점수의 순위는 (비흡연자 + 현재 흡연자) 집단에서의 순위와 크게 달라진다. 동시에 분산도 크게 달라지기 때문에 적절하지 않다.
2) Conover Test
- 상기 문제점을 해결하였고, t 분포에 근거해 판단한다.
3) Dunn Test
-상기 문제점을 해결하였고, z 분포에 근거해 판단한다.
4) Nemenyi Test
하지만 SAS는 위 네 개 중 하나의 방법이 아닌 다른 방법을 사용한다.
Dwass, Steel, Critchlow-Fligner (DSCF) multiple comparisonprocedure를 사용하며 이는 pairwise Wilcoxon test를 기반으로 하나, 특수한 방법으로 그 familywise error rate를 보정하고 있다. 이 방법은 Conover Test보다 더 우수한 것으로 받아들여지고 있다.
DSCF 사후 검정 - 코드
PROC NPAR1WAY DATA=hong.df WILCOXON DSCF CORRECT=NO;
CLASS SMOK;
VAR NICOT_ADDICT;
RUN;
일반적인 Kruskal Wallis 검정을 시행하되, 첫 번째 줄에 "DSCF"를 추가해 준 것이 끝이다.
결과
여러 값이 제시되지만, 마지막 열인 "Pr>DSCF"인 p-value만 보면 된다.
(1) 과거 흡연자(SMOK=1) vs 현재 흡연자 (SMOK=2) : p-value < 0.0001
(2) 과거 흡연자(SMOK=1) vs 비흡연자(SMOK=0):p-value < 0.0001
(3) 현재 흡연자(SMOK=2) vs 비흡연자(SMOK=0):p-value < 0.0001
모든 p-value가 0.05보다 작으므로 다음과 같이 결론 내릴 수 있다.
비흡연자 (0)와 과거 흡연자 (1) 사이에는 NICOT_ADDICT 중앙값에 차이가 있다.
비흡연자 (0)와 현재 흡연자 (2) 사이에는 NICOT_ADDICT 중앙값에 차이가 있다.
과거 흡연자 (1)와 현재 흡연자 (2)) 사이에는 NICOT_ADDICT 중앙값에 차이가 있다.
어떤 사후 분석을 쓸 것인가
이 논의에 대해 정답이 따로 있는 것은 아니다. 적절한 방법을 사용하여 논문에 제시하면 되고, 어떤 것이 정답이라고 콕 집어 이야기할 수는 없다. 다만, 사후 분석 방법이 여러 가지가 있다는 것은 '사후 분석 방법에 따라 산출되는 결과가 달라질 수 있다.'는 것을 의미하고, 심지어는 '어떤 사후 분석 방법을 채택하냐에 따라 유의성 여부가 달라질 수도 있다.'는 것을 의미한다. 심지어, Kruskal-Wallis test에서는 유의한 결과가 나왔는데, 사후 분석을 해보니 유의한 차이를 보이는 경우가 없을 수도 있다. 따라서 어떤 사후 분석 방법 결과에 따른 결과인지 유의하여 해석할 필요가 있다.
왜 굳이 Kruskal Wallis test를 쓰는 것인가?
이 시점에서 이런 의문이 들 수 있다.
"각각의 그룹별로 평균을 비교하면 되지, 굳이 왜 Kruskal Wallis test라는 방법까지 사용하는 것인가?"
아주 논리적인 의문점이다. 하지만, 반드시 Kruskal-Wallis test를 사용해야 한다. 그 이유는 다음과 같다. 본 사례는 흡연 상태에 따른 조합 가능한 경우의 수가 3인데, 각각 유의성의 기준을 $\alpha=0.05$로 잡아보자. 이때 세 번의 비교에서 모두 귀무가설이 학문적인 진실인데(평균에 차이가 없다.), 세 번 모두 귀무가설을 선택할 확률은 $\left( 1-0.05 \right)^3 \approx 0.86$이다. (이해가 어려우면 p-value에 대한 설명 글을 읽고 오길 바란다. 2022.09.05 - [통계 이론] - [이론] p-value에 관한 고찰)
그런데, Kruskal-Wallis test의 귀무가설은 "모든 집단의 중앙값이 같다."이다. 따라서 모든 집단의 중앙값이 같은 것이 학문적 진실일 때, 적어도 한 번이라도 대립 가설을 선택하게 될 확률은 $1-0.86=0.14$가 된다. 즉, 유의성의 기준이 올라가게 되어, 덜 보수적인 결정을 내리게 되고, 다시 말하면 대립 가설을 잘 선택하는 쪽으로 편향되게 된다. 학문적으로는 '다중 검정 (multiple testing)을 시행하면 1종 오류가 발생할 확률이 증가하게 된다.'라고 표현한다.
따라서, 각각을 비교해 보는 것이 아니라 한꺼번에 비교하는 Kruskal-Wallis test를 시행해야 함이 마땅하다.
여기에서 한 번 더 의문이 들 수 있다.
"사후 분석을 할 때에는 1종 오류가 발생하지 않는가?"
그렇다. 1종 오류가 발생할 확률이 있으므로, p-value의 기준을 더 엄격하게 (0.05보다 더 작게) 잡아야 한다. P-value를 보정하는 방법은 일반적으로는 Bonferroni, holm, hochberg, hommel, BH, BY, FDR 등이 있고, Kruskal-Wallis test의 사후 분석 방법으로는 특별히 Pairwise Wilcoxon Test,Conover Test,Dunn Test,Nemenyi Test, DSCF을 사용하고 있다.
코드 정리
*라이브러리 지정하기;
LIBNAME hong "C:/Users/User/Documents/Tistory_blog";
*파일 불러오기;
PROC IMPORT
DATAFILE="C:\Users\user\Documents\Tistory_blog\Data.xlsx"
DBMS=EXCEL
OUT=hong.df
REPLACE;
RUN;
*정규성 검정하기;
PROC UNIVARIATE DATA=hong.df NORMAL PLOT;
CLASS SMOK;
VAR NICOT_ADDICT;
HISTOGRAM NICOT_ADDICT/ NORMAL (MU=EST SIGMA=EST);
RUN;
*크루스칼 왈리스 검정;
PROC NPAR1WAY DATA=hong.df WILCOXON CORRECT=NO;
CLASS SMOK;
VAR NICOT_ADDICT;
RUN;
*크루스칼 왈리스 검정 & 사후검정;
PROC NPAR1WAY DATA=hong.df WILCOXON DSCF CORRECT=NO;
CLASS SMOK;
VAR NICOT_ADDICT;
RUN;
[SAS] 크루스칼 왈리스 검정 (비모수 일원 배치 분산 분석: Kruskal-Wallis Test) 정복 완료!
[SAS] 윌콕슨 순위 합 검정, 맨 휘트니 U 검정 (비모수 독립 표본 중앙값 검정: Wilcoxon rank sum test, Mann-Whitney U test) - NPAR1WAY()
두 분포의 평균이 다른지 확인하는 방법을 이전에는 독립 표본 T검정 (Two-Sample T test)로 시행했었다. (2022.10.04 - [모평균 검정/SAS] - [SAS] 독립 표본 T검정 (Independent samples T-test) - PROC TTEST) 하지만 여기에는 중요한 가정이 필요한데, 각각의 분포가 정규성을 따르는 것이다. 하지만 분포가 정규성을 따르지 않는다면 어떻게 해야 할까? 그럴 때 사용하는 것이 Wilcoxon rank sum test (윌콕슨 순위 합 검정)이다. 윌콕슨 순위 합 검정 (Wilcoxon rank sum test)는 이후에 Mann과 Whitney가 개정을 하였고 그때 U 통계량을 사용하므로 맨 휘트니 U 검정(Mann-Whitney U test)이라고도 한다. 또한 세 사람의 이름을 붙여 Wilcoxon-Mann-Whitney Test (WMW test)라고도 한다.
이번 포스팅에서는 윌콕슨 순위 합 검정 (Wilcoxon rank sum test)에 대해 알아볼 것이다.
목표: 음주 여부에 따라 Gamma-glutamyl transferase(GGT)의 중앙값이 모집단 수준에서 서로 다르다고 말할 수 있는가?
정규성 - 코드
만약 GGT의 분포가 정규성을 따른다면, 독립 표본 T 검정으로 이를 확인할 수 있을 것이다. 독립 표본 (Indepent samples, Two tamples) T 검정의 전제조건은 각각의 독립 표본의 분포가 정규성을 따른다는 것이다. 즉, 여기에서는 음주 여부에 따라 GGT의 정규성을 검정하도록 하겠다.
PROC UNIVARIATE DATA=hong.df NORMAL PLOT;
CLASS ALCOHOL;
VAR GGT;
HISTOGRAM GGT/ NORMAL (MU=EST SIGMA=EST);
RUN;
정규성 검정 - 결과
1) ALCOHOL=0 (비음주자 군)
2) ALCOHOL=1 (음주자 군)
N수가 2,000개 미만이므로 Shapiro-Wilk 통계량의 p-value를 보면 0.05 이하이며, Q-Q Plot은 대부분의 데이터가 선상에 있지 않고, 히스토그램에서도 정규성을 따르지 않는 것처럼 보인다. 따라서 독립 표본 T검정 (Two-sample T-test)를 시행할 수 없고, 맨 휘트니 U 검정 (Mann-Whitney U test)을 시행해야 한다.
맨 휘트니 U 검정 (Mann-Whitney U test) 코드
PROC NPAR1WAY DATA=hong.df WILCOXON CORRECT=NO;
CLASS ALCOHOL;
VAR GGT;
RUN;
빨간 박스 안 양측 검정 p-value가 매우 작으므로 유의한 결과임을 알 수 있다. 따라서 귀무가설을 기각하고 대립 가설을 택해야 한다. 여기에서 귀무가설과 대립 가설은 무엇일까?
귀무가설 $H_0=$ 두 분포는 동일하다.
대립 가설$H_1=$ 두 분포는 평균만 다를 뿐, 분포는 똑같이 생겼다.
즉, 대립 가설이 의미하는 것은 두 분포는 평행이동 관계에 있다는 것이다. 그래야 중앙값이 달라진 것인지 확인할 수 있다.
평행이동 관계
이 말은, 맨 휘트니 U 검정에도 필요한 가정이 있다는 말이다. 두 분포는 동일하게 생겼어야 한다.
조금 어렵게 이야기하면, 두 분포의 first moment는 다르지만, 그보다 고차원의 central moment는 같아야 한다.
만약 생김새가 다르면 맨 휘트니 U 검정을 사용하면 안 되고, Robuts rank order test 등을 시행해야 한다. (물론 현실에서는 이런 내용을 무시한 채 '정규성을 따르지 않으니 맨 휘트니 U 검정을 사용한다'는 연구자가 대다수다.) 이에 관한 내용은 다음 링크를 확인하길 바란다. 링크 추가 예정
그런데 위 히스토그램을 보면 음주자와 비음주자의 GGT분포는 right skewed 된 분포로 어느 정도 비슷하게 생겼다. 모든 central moment를 조사하는 것은 현실적으로 힘들지만, 2nd cental moment인 분산의 차이를 검정해보면 다음과 같음을 알 수 있다. 등분산성 검정 방법은 다음 링크를 확인하길 바란다.
R에서 NP3_raw는 right skewed, NP4_raw는 left skewed로 서루 매우 다르게 생겼음을 알 수 있다.
분산, 중앙값이 완벽히 일치하며, 등분산성 검정에서도 분산에 차이가 없다고 나온다. 하지만 맨 휘트니 U 검정 결과 p-value는 매우 작게 나온다. 만약 맨 휘트니 U 검정이 "중앙값이 같은지 확인하는 검정"이라고만 알고 있다면 중앙값이 완전히 같은 이런 경우에도 "중앙값에 차이가 있다."라고 결론 내릴 것이다. 하지만 이 경우 귀무가설을 기각하고 "두 분포는 다르게 생겼다."라고 결론 내려야 옳은 결론에 다다른다.
따라서 만약 생긴 모양이 같고 중앙값만 다르다면, 유의미한 맨 휘트니 U 검정의 결과는 중앙값의 차이로 귀결될 수 있다. 하지만 모양이 다르다면 중앙값의 차이로 귀결될 수 없음을 유의하길 바란다.
코드 정리
*라이브러리 지정하기;
LIBNAME hong "C:/Users/User/Documents/Tistory_blog";
*파일 불러오기;
PROC IMPORT
DATAFILE="C:\Users\user\Documents\Tistory_blog\Data.xlsx"
DBMS=EXCEL
OUT=hong.df
REPLACE;
RUN;
*정규성 검정;
PROC UNIVARIATE DATA=hong.df NORMAL PLOT;
CLASS ALCOHOL;
VAR GGT;
HISTOGRAM GGT/ NORMAL (MU=EST SIGMA=EST);
RUN;
*맨 휘트니 U 검정;
PROC NPAR1WAY DATA=hong.df WILCOXON CORRECT=NO;
CLASS ALCOHOL;
VAR GGT;
RUN;
[SAS] 윌콕슨 순위 합 검정, 맨 휘트니 U 검정 (비모수 독립 표본 중앙값 검정: Wilcoxon rank sum test, Mann-Whitney U test) 정복 완료!
PROC UNIVARIATE DATA=hong.df NORMAL PLOT;
VAR NP1_raw;
HISTOGRAM NP1_raw/ NORMAL (MU=EST SIGMA=EST);
RUN;
정규성 검정 - 결과
N수가 2,000개 미만이므로 Shapiro-Wilk 통계량의 p-value를 보면 0.05 이하이며, Q-Q Plot은 대부분의 데이터가 선상에 있지 않고, 히스토그램에서도 정규성을 따르지 않는 것처럼 보인다. 따라서 일표본 T검정 (One-sample T-test)를 시행할 수 없고, 일표본 윌콕슨 부호 순위 검정 (One-Sample Wilcoxon Signed Rank Test)을 시행해야 한다.
일표본 윌콕슨 부호 순위 검정 (One-Sample Wilcoxon Signed Rank Test)의 귀무가설과 대립 가설
이번 일표본 윌콕슨 부호 순위 검정 (One-Sample Wilcoxon Signed Rank Test)의 귀무가설, 대립 가설은 다음과 같다.
귀무가설: $H_0=$ 모집단의 NP1_raw 중앙값은 55이다.
대립 가설: $H_1=$ 모집단의 중앙값은 55보다 크거나 작다. (양측 검정)
일표본 윌콕슨 부호 순위 검정 (One-Sample Wilcoxon Signed Rank Test) 코드
PROC UNIVARIATE DATA=hong.df MU0=55;
VAR NP_1_raw;
RUN;
PROC UNIVARIATE DATA=hong.df MU0=55; : hong 라이브러리의 df데이터를 사용하고, 검정할 중앙값은 55다. VAR NP1_raw; : NP1_raw에 대해 검정을 시행하라.
결과
우리가 보아야 할 내용은 부호 순위의 p값이다. 0.0001보다 작으므로 귀무가설을 기각하고 대립 가설을 선택하며 중앙값이 55는 아니라고 결론 내릴 수 있다.
알려진 중앙값에 대해 일표본 윌콕슨 부호 순위 검정 코드
위 정규성 검정에서 중앙값은 49.96328로 거의 50에 근사한다. 그렇다면 50에 대해 검정하면 어떤 결과를 내보일까?
PROC UNIVARIATE DATA=hong.df MU0=50;
VAR NP1_raw;
RUN;
결과
p-value가 0.6560>0.05로 귀무가설을 채택하는 것이 합리적이며, 중앙값이 50이 아니라고는 할 수 없는 상황이다.
중요한 가정: 대칭성
위 논의에서 빠진 정말 중요한 가정이 하나 있다. 신경심리검사 1 원점수 (NP1_raw)의 분포가 어떤 점수를 기점으로 좌우대칭이라는 것이다. 윌콕슨 부호 순위 검정은 좌우대칭일 때 그 검정력이 극대화되며, 좌우대칭이 아니면 정확도가 떨어진다. 신경심리검사 1 원점수 (NP1_raw)의 분포를 히스토그램을 통해 알아보자. 이는 위의 정규성 검정에서 나와있으며 결과는 다음과 같다.
봉우리가 세 개인 Trimodal shape을 하고 있고, 중앙값인 50을 기준으로 대칭인 것처럼 보인다. 따라서 위 분석에는 문제가 없는 듯하다.
가정이 성립하지 않는다면?
우리가 궁금하고, 문제가 되는 상황은 가정이 성립하지 않는 상황이다. 실제로 이럴 때가 많기 때문이다. 신경심리검사 2 원점수인 NP2_raw데이터의 분포를 확인해보자.
PROC UNIVARIATE DATA=hong.df NORMAL PLOT;
VAR NP2_raw;
HISTOGRAM NP2_raw/ NORMAL (MU=EST SIGMA=EST);
RUN;
이는 심하게 오른쪽으로 치우쳐진 right skewed data다. 이때 중위수(1.004532)에 대해 Wilcoxon signed rank test를 시행하면 어떻게 될까?
PROC UNIVARIATE DATA=hong.df MU0=1.004532;
VAR NP2_raw;
RUN;
참 중앙값으로 "NP2_raw의 중앙값이 참 중앙값과 같니?"라고 물어보았는데 "아니요"라고 대답하는 참사가 벌어진다. 이는 분포에 대칭성이 없기 때문이다. 이럴 때 선택할 수 있는 옵션이 많지 않긴 한데, 그중 하나는 변환 (transformation)이다. 이런 right skewed data는 로그 변환 (log-transformation)을 하면 대칭성이 갖추어지는 경우가 많다. 로그 변환을 하고 히스토그램을 그려 분포를 확인해보자. (로그 변환에 관한 글은 다음 링크에서 확인할 수 있다.2022.10.06 - [통계 프로그램 사용 방법/SAS] - [SAS] 변수 계산 (산술 연산) - DATA, SET)
DATA hong.df;
SET hong.df;
log_NP2_raw=log(NP2_raw);
RUN;
PROC UNIVARIATE DATA=hong.df NORMAL PLOT;
VAR log_NP2_raw;
HISTOGRAM log_NP2_raw/ NORMAL (MU=EST SIGMA=EST);
RUN;
충분히 대칭성을 띄고 있다고 할 수 있다. "혹시 정규성을 띄지 않을까?"라는 생각으로 정규성을 확인해보면 다음과 같음을 알 수 있다.
히스토그램에서는 정규분포를 따르는 것 같아 보이지만, QQ plot에서 직선 위에 있지 않은 점들이 많고, Shapiro-wilk test에서 p-value가 매우 작게 나오므로 정규성을 띠지 않는다고 판단하는 것이 맞다. 따라서 일표본 t검정(One-sample t-test)이 아닌 일표본 윌콕슨 부호 순위 검정 (One sample Wilcoxon signed rank test)을 시행해야 한다. 위와 같이 log_NP2_raw의 중앙값에 대해 검정해보자.
코드
PROC UNIVARIATE DATA=hong.df MU0=0.00452;
VAR log_NP2_raw;
RUN;
p-value가 0.7366으로 귀무가설을 기각할 수 없게 되어 중앙값은 0.00452가 아니라고 할 수 없다고 결론을 내리게 된다.
그 외 다룰 내용 1: 정규성을 안 따르면 꼭 Wilcoxon test를 시행해야 하나?
T-test의 전제조건은 표본 평균이 정규분포를 따르는 것이다. 그런데 중심 극한 정리에 따라 표본의 수 (우리 데이터에서는 n=1,000)가 커질수록 표본 평균의 분포는 정규성을 띠게 된다. 따라서 n수가 적당히 크기만 하면 표본이 정규분포를 따르는지와 관계없이 표본평균은 정규성을 따른다고 할 수 있다. 이런 이유로 정규성과 관계없이 n수가 크기만 하면 t-test의 이용이 정당화되기도 한다. 합리적인 말이다.
그런데 왜 정규분포를 따르는지 왜 다들 확인할까? 만약 모분포가 정규분포였다면 표본도 정규분포를 따를 확률이 높아진다. 물론 꼭 그런 건 아니지만 말이다. 따라서 표본이 정규분포를 따른다면, 모분포를 정규분포를 따른다고 할 수 있고, 그렇다면 n수가 작아도 표본 평균의 분포도 정규성을 따르니 t-test를 사용해도 되는 안전한 환경이 구축된다. 이런 이유로 정규성을 확인하게 된다.
그 외 다룰 내용 2: Wilcoxon signed rank test는 짝지어진 자료의 검정방법 아닌가?
1) 통계량: SPSS와 R은 계산된 통계량을 바로 제시한다. 하지만 SAS는 부호 순위 통계량의 기댓값인 $\frac {n(n+1)}{4}$를 뺀 값을 제시한다. 따라서 다를 수밖에 없다.
2) R에서는 중앙값 자체를 코드에 넣었지만, 본 포스팅에서는 중앙값이라고 제시된 반올림된 값을 코드에 넣었다.
3) (가장 중요한 내용) R과 SPSS는 윌콕슨 부호 순위 통계량을 정규분포에 근사시키지만, SAS는 t 분포에 근사 시키며 이용하는 분산도 아주 조금은 다르다.
4) 하지만, 표본의 수가 20이하일 때에는 근사 시키지 않은 정확한 p-value (exact p-value)을 사용한다. R에서도 50 미만인 표본에 대해서는 exact method를 사용하므로 20 미만의 표본의 수를 가진 데이터로 검정할 때에는 같은 결과를 내주게 된다.
코드 정리
*라이브러리 지정하기;
LIBNAME hong "C:/Users/User/Documents/Tistory_blog";
*파일 불러오기;
PROC IMPORT
DATAFILE="C:\Users\user\Documents\Tistory_blog\Data.xlsx"
DBMS=EXCEL
OUT=hong.df
REPLACE;
RUN;
**NP1_raw;
*정규성 검정;
PROC UNIVARIATE DATA=hong.df NORMAL PLOT;
VAR NP1_raw;
HISTOGRAM NP1_raw/ NORMAL (MU=EST SIGMA=EST);
RUN;
*Wilcoxon signed rank test;
PROC UNIVARIATE DATA=hong.df MU0=55;
VAR NP_1_raw;
RUN;
*중앙값 근사치(50)으로 검정하기;
PROC UNIVARIATE DATA=hong.df MU0=50;
VAR NP1_raw;
RUN;
**NP2_raw;
*NP2_raw 대칭성 확인;
PROC UNIVARIATE DATA=hong.df NORMAL PLOT;
VAR NP2_raw;
HISTOGRAM NP2_raw/ NORMAL (MU=EST SIGMA=EST);
RUN;
*NP2_raw 비대칭이지만 wilcoxon test 시행하기;
PROC UNIVARIATE DATA=hong.df MU0=1.004532;
VAR NP2_raw;
RUN;
*로그변환 후 히스토그램으로 대칭성 확인하기, 정규성 검정하기;
DATA hong.df;
SET hong.df;
log_NP2_raw=log(NP2_raw);
RUN;
PROC UNIVARIATE DATA=hong.df NORMAL PLOT;
VAR log_NP2_raw;
HISTOGRAM log_NP2_raw/ NORMAL (MU=EST SIGMA=EST);
RUN;
*로그변환된 변수로 Wilcoxon test 시행하기;
PROC UNIVARIATE DATA=hong.df MU0=0.00452;
VAR log_NP2_raw;
RUN;
[SAS] 일표본 윌콕슨 부호 순위 검정 (비모수 일표본 중앙값 검정: One-Sample Wilcoxon Signed Rank Test) 정복 완료!
그런데, 여기에는 각각의 분포가 정규성을 따른다는 중요한 가정이 필요하다. 하지만 정규성을 따르지 않는다면 어떻게 해야 할까? 그럴 때 사용하는 것이 크루스칼-왈리스 검정(Kruskal-Wallis Test)이다. 이번 포스팅에서는 크루스칼 왈리스 검정과 그 사후 분석 (post hoc analysis)까지 알아보고자 한다.
2) 분석하고자 하는 변수인 NICOT_ADDICT을 "종속변수"에 넣고, 흡연 상태에 따라 분석할 것이므로 SMOK을 "요인(F)"에 넣는다. 그 뒤 "도표(T)..."를 선택한다.
3) "히스토그램(H)", "검정과 함께 정규성 도표(O)" 체크박스를 클릭하고 "계속(C)"를 누르고, 돌아가 "확인"을 누른다.
결과
1) Q-Q Plot
2) 히스토그램
3) Shapiro-Wilk 검정
N수가 2,000개 미만이므로 Shapiro-Wilk 통계량의 p-value를 보면 0.05 이하이며, Q-Q Plot은 대부분의 데이터가 선상에 있지 않고, 히스토그램에서도 정규성을 따르지 않는 것처럼 보인다. 따라서 일원 배치 분산 분석 (ANOVA)을 시행할 수 없고, 크루스칼 왈리스 분석 (Kruskal-Wallis Test)을 시행해야 한다.
지금까지 많은 분석을 할 때 분산이 같은지 여부 (등분산성)을 중요하게 여겼었다. 하지만 크루스칼 왈리스 분석은 분산이 달라도 어느 정도 괜찮다. 물론 분산에 아주 큰 차이가 나면 위험해지지만, 어느 정도는 robust 하므로 웬만하면 그냥 써도 괜찮다.
크루스칼 왈리스 (Kruskal Wallis) 검정
1) 분석(A) > 비모수검정(N) > 독립표본(I)
2) 이때 나오는 창의 첫 페이지인 "목적"은 건들지 않는다.
3) "필드"를 누르고 분석하고자 하는 NICOT_ADDICT를 "검정 필드(T)"로 넘긴다. 흡연 상태에 따라 검정할 것이므로 SMOK을 "집단(G)"로 옮긴다.
4) "사용자 정의에 의한 검정(C)"를 누르고 "Kruskal-Wallis 일원배치 분산분석(k표본)(W)"을 체크한다.
결과
p-value가 0.000 (<0.001)으로 0.05보다 작으므로 귀무가설을 기각하고 대립 가설을 채택한다. 그렇다면 여기에서 귀무가설 및 대립 가설은 무엇이었는가?
귀무가설: $H_0=$ 세 집단의 모집단 수준에서 중앙값은 "모두" 동일하다.
대립가설:$H_1=$ 세 집단의 모집단 수준에서 중앙값이 모두 동일한 것은 아니다.
우리는 대립 가설을 채택해야 하므로 "세 집단의 모집단 수준에서 중앙값이 모두 동일한 것은 아니다."라고 결론 내릴 것이다.
사후 분석
그런데, 세 집단의 중앙값이 모두 동일하지 않다는 말은 세 집단 중 두 개씩 골라 비교했을 때, 적어도 한 쌍에서는 차이가 난다는 것이다. 따라서 세 집단 중 두 개씩 골라 비교를 해보아야 하며, 이를 사후 분석 (post hoc analysis)이라고 한다. 세 집단에서 두 개씩 고르므로 가능한 경우의 수 $_3C_2=3$이다.
(1) 비흡연자 vs 과거 흡연자
(2) 비흡연자 vs 현재 흡연자
(3) 과거 흡연자 vs 현재 흡연자
비모수 검정인 Kruskal Wallis 검정의 사후 분석으로 쓰이는 방법은 다음 네 가지가 있다.
1) Pairwise Wilcoxon Test
2) Conover Test
3) Dunn Test
4) Nemenyi Test
이중 SPSS에서 기본으로 제공하는 사후검정 방법은 Dunn Test이다. 이때 p value 보정 방법은 Bonferroni를 사용하고 있다.
그리고, Kruskal-Wallis에서 유의한 차이를 보이는 경우 SPSS는 사후분석을 자동으로 실행한다.
위 표에서 봐야 하는 값은 마지막 열인 "조정된 유의확률"이다.
0-1
비흡연자 (0)와 과거 흡연자 (1) 비교
0-2
비흡연자 (0)와 현재 흡연자 (2) 비교
1-2
과거 흡연자 (1)와 현재 흡연자 (2) 비교
모든 p-value가 0.05보다 작으므로 다음과 같이 결론 내릴 수 있다.
비흡연자 (0)와 과거 흡연자 (1) 사이에는 NICOT_ADDICT 중앙값에 차이가 있다.
비흡연자 (0)와 현재 흡연자 (2) 사이에는 NICOT_ADDICT 중앙값에 차이가 있다.
과거 흡연자 (1)와 현재 흡연자 (2)) 사이에는 NICOT_ADDICT 중앙값에 차이가 있다.
하지만 위 결과에서 볼 수 있듯이, 영가설 (귀무가설)은 "표본1 및 표본2에서 분포가 같다"이다. 만약 분포가 똑같이 생겼는데, 중앙값만 다른 경우라면, 즉 평행이동을 한 경우라면 귀무가설의 기각은 '중앙값에 차이가 난다'를 의미한다.
이 논의에 대해 정답이 따로 있는 것은 아니다. 적절한 방법을 사용하여 논문에 제시하면 되고, 어떤 것이 정답이라고 콕 집어 이야기할 수는 없다. 다만, 사후 분석 방법이 여러 가지가 있다는 것은 '사후 분석 방법에 따라 산출되는 결과가 달라질 수 있다.'는 것을 의미하고, 심지어는 '어떤 사후 분석 방법을 채택하냐에 따라 유의성 여부가 달라질 수도 있다.'는 것을 의미한다. 심지어, Kruskal-Wallis test에서는 유의한 결과가 나왔는데, 사후 분석을 해보니 유의한 차이를 보이는 경우가 없을 수도 있다. 따라서 어떤 사후 분석 방법 결과에 따른 결과인지 유의하여 해석할 필요가 있다.
왜 굳이 Kruskal Wallis test를 쓰는 것인가?
이 시점에서 이런 의문이 들 수 있다.
"각각의 그룹별로 평균을 비교하면 되지, 굳이 왜 Kruskal Wallis test라는 방법까지 사용하는 것인가?"
아주 논리적인 의문점이다. 하지만, 반드시 Kruskal-Wallis test를 사용해야 한다. 그 이유는 다음과 같다. 본 사례는 흡연 상태에 따른 조합 가능한 경우의 수가 3인데, 각각 유의성의 기준을 α=0.05α=0.05로 잡아보자. 이때 세 번의 비교에서 모두 귀무가설이 학문적인 진실인데(평균에 차이가 없다.), 세 번 모두 귀무가설을 선택할 확률은 $\left(1-0.05 \right)^3 \approx 0.86 $이다. (이해가 어려우면 p-value에 대한 설명 글을 읽고 오길 바란다. 2022.09.05 - [통계 이론] - [이론] p-value에 관한 고찰)
그런데, Kruskal-Wallis test의 귀무가설은 "모든 집단의 중앙값이 같다."이다. 따라서 모든 집단의 중앙값이 같은 것이 학문적 진실일 때, 적어도 한 번이라도 대립 가설을 선택하게 될 확률은 $1-0.86=0.14$ 가 된다. 즉, 유의성의 기준이 올라가게 되어, 덜 보수적인 결정을 내리게 되고, 다시 말하면 대립 가설을 잘 선택하는 쪽으로 편향되게 된다. 학문적으로는 '다중 검정 (multiple testing)을 시행하면 1종 오류가 발생할 확률이 증가하게 된다.'라고 표현한다.
따라서, 각각을 비교해보는 것이 아니라 한꺼번에 비교하는 Kruskal-Wallis test를 시행해야 함이 마땅하다.
[SPSS] 윌콕슨 순위 합 검정, 맨 휘트니 U 검정 (비모수 독립 표본 중앙값 검정: Wilcoxon rank sum test, Mann-Whitney U test)
두 분포의 평균이 다른지 확인하는 방법을 이전에는 독립 표본 T검정 (Two-Sample T test)로 시행했었다. (2022.11.12 - [모평균 검정/R] - [R] 독립 표본 T검정 (Independent samples T-test) - t.test(), var.test(), levene.test()) 하지만 여기에는 중요한 가정이 필요한데, 각각의 분포가 정규성을 따르는 것이다. 하지만 분포가 정규성을 따르지 않는다면 어떻게 해야 할까? 그럴 때 사용하는 것이 Wilcoxon rank sum test (윌콕슨 순위 합 검정)이다. 윌콕슨 순위 합 검정 (Wilcoxon rank sum test)는 이후에 Mann과 Whitney가 개정을 하였고 그때 U 통계량을 사용하므로 맨 휘트니 U 검정(Mann-Whitney U test)이라고도 한다. 또한 세 사람의 이름을 붙여 Wilcoxon-Mann-Whitney Test (WMW test)라고도 한다.
이번 포스팅에서는 윌콕슨 순위 합 검정 (Wilcoxon rank sum test)에 대해 알아볼 것이다.
2) 분석하고자 하는 변수인 GGT을 "종속변수"에 넣고, 음주 여부에 따라 분석할 것이므로 ALCOHOL을 "요인(F)"에 넣는다. 그 뒤 "도표(T)..."를 선택한다.
3) "히스토그램(H)", "검정과 함께 정규성 도표(O)" 체크박스를 클릭하고 "계속(C)"를 누르고, 돌아가 "확인"을 누른다.
결과
1) Q-Q Plot
2) 히스토그램
3) Shapiro-Wilk 검정
N수가 2,000개 미만이므로 Shapiro-Wilk 통계량의 p-value를 보면 0.05 이하이며, Q-Q Plot은 대부분의 데이터가 선상에 있지 않고, 히스토그램에서도 정규성을 따르지 않는 것처럼 보인다. 따라서 독립 표본 T검정 (Two-sample T-test)를 시행할 수 없고, 맨 휘트니 U 검정 (Mann-Whitney U test)을 시행해야 한다.
맨 휘트니 U 검정 (Mann-Whitney U test)
1) 분석(A) > 비모수검정(N) > 독립표본(I)
2) 이때 나오는 창의 첫 페이지인 "목적"은 건들지 않는다.
3) "필드"를 누르고 분석하고자 하는 GGT를 "검정 필드(T)"로 넘긴다. 음주 여부에 따라 검정할 것이므로 ALCOHOL을 "집단(G)"로 옮긴다.
4) "사용자 정의에 의한 검정(C)"를 누르고 "Mann-Whitney의 U(2표본)(H)"을 체크한다.
결과
유의확률이 0.000으로 0.001보다 작다. 따라서 유의한 결과임을 알 수 있다. 따라서 귀무가설을 기각하고 대립 가설을 택해야 한다. 여기에서 귀무가설과 대립 가설은 무엇일까?
귀무가설$H_0=$ 분포는 동일하다.
대립가설$H_1=$ 두 분포는 평균만 다를 뿐, 분포는 똑같이 생겼다.
즉, 대립 가설이 의미하는 것은 두 분포는 평행이동 관계에 있다는 것이다. 그래야 대립가설을 선택하더라도 중앙값이 달라진 것인지 확인할 수 있다.
평행이동 관계
이 말은, 맨 휘트니 U 검정에도 필요한 가정이 있다는 말이다. 두 분포는 동일하게 생겼어야 한다.
조금 어렵게 이야기하면, 두 분포의 first moment는 다르지만, 그보다 고차원의 central moment는 같아야 한다.
만약 맨 휘트니 U 검정이 "중앙값이 같은지 확인하는 검정"이라고만 알고 있다면 중앙값이 완전히 같은 이런 경우에도 "중앙값에 차이가 있다."라고 결론 내릴 것이다. 하지만 이 경우 귀무가설을 기각하고 "두 분포는 다르게 생겼다."라고 결론 내려야 옳은 결론에 다다른다. 그러므로 맨 휘트니 U 검정의 가정인 "모양이 똑같이 생겼다."를 반드시 준수하길 바란다.
[SPSS] 윌콕슨 순위 합 검정, 맨 휘트니 U 검정 (비모수 독립 표본 중앙값 검정: Wilcoxon rank sum test, Mann-Whitney U test) 정복 완료!
어떤 분포의 평균이 특정 값인지 확인하는 방법을 이전에는 일표본 T검정 (One-Sample T test)로 시행했었다. (2022.11.29 - [모평균 검정/SPSS] - [SPSS] 일표본 T검정 (One-sample T-test)) 하지만 여기에는 중요한 가정이 필요한데, 분포가 정규성을 따르는 것이다. 하지만 분포가 정규성을 따르지 않는다면 어떻게 해야 할까? 그럴 때 사용하는 것이 One-Sample Wilcoxon Signed Rank Test (일표본 윌콕슨 부호 순위 검정)이다.
이번 포스팅에서는 One-Sample Wilcoxon Signed Rank Test (일표본 윌콕슨 부호 순위 검정)에 대해서 알아볼 것이다.
2) 분석하고자 하는 변수인 NP1_raw을 "종속변수"에 넣어준다. 그 뒤 "도표(T)..."를 선택한다.
3) "히스토그램(H)", "검정과 함께 정규성 도표(O)" 체크박스를 클릭하고 "계속(C)"를 누르고, 돌아가 "확인"을 누른다.
결과
N수가 2,000개 미만이므로 Shapiro-Wilk 통계량의 p-value를 보면 0.05 이하이며, Q-Q Plot은 대부분의 데이터가 선상에 있지 않고, 히스토그램에서도 정규성을 따르지 않는 것처럼 보인다. 따라서 일표본 T검정 (One-sample T-test)를 시행할 수 없고, 일표본 윌콕슨 부호 순위 검정 (One-Sample Wilcoxon Signed Rank Test)을 시행해야 한다.
일표본 윌콕슨 부호 순위 검정 (One-Sample Wilcoxon Signed Rank Test)
이번 일표본 윌콕슨 부호 순위 검정 (One-Sample Wilcoxon Signed Rank Test)의 귀무가설과 대립 가설은 다음과 같다.
귀무 가설:$H_0=$ 모집단의 NP1_raw 중앙값은 55이다.
대립 가설:$H_1=$ 모집단의 중앙값은 55보다 크거나 작다. (양측 검정)
1) 분석(A) > 비모수검정(N) > 일표본(O)
2) 이때 나오는 창의 첫 페이지인 "목적"은 건들지 않는다.
3) "필드"를 누르고 분석하고자 하는 NP1_raw를 오른쪽으로 넘긴다.
4) "사용자 정의에 의한 검정(T)"를 누르고 "평균과 가설값 비교(Wilcoxon 부호 순위 검정)"을 체크하고, 가설 중위수에는 검정하고자 하는 값인 55를 적는다. 그리고 "실행"을 누른다.
결과
p-value=0.000으로 이는 0.001보다 작다는 뜻이며, 0.05보다 작으므로 귀무가설을 기각하고 대립 가설을 받아들인다. 대립 가설은 중위수가 55가 아니라고 설정했었으므로 적어도 중위수가 55는 아니라고 이야기할 수 있다. 그렇다면 중위수가 얼마였길래 그럴까?
중위수 확인
1) 분석(A) > 기술통계량(E) > 데이터 탐색 (E)
2) 분석하고자 하는 변수인 NP1_raw을 "종속변수"에 넣어준다. 그 뒤 "확인"을 누른다.
결과
중위수는 49.96328279로 거의 50에 가깝다. 그렇다면 Wilcoxon 부호 순위 검정을 50에 대해 시행하면 어떻게 될까? 가설 중위수에 50을 넣고 실행을 눌러보자
결과
p-value>0.05이므로 귀무가설을 채택하여 50과 다르다고 할 수 없다는 결론을 얻는다.
중요한 가정: 대칭성
위 논의에서 빠진 정말 중요한 가정이 하나 있다. 신경심리검사 1 원점수 (NP1_raw)의 분포가 어떤 점수를 기점으로 좌우대칭이라는 것이다. 윌콕슨 부호 순위 검정은 좌우대칭일 때 그 검정력이 극대화되며, 좌우대칭이 아니면 정확도가 떨어진다. 신경심리검사 1 원점수 (NP1_raw)의 분포를 히스토그램을 통해 알아보고자 하며 이는 위에 구해놓았다.
봉우리가 세 개인 Trimodal shape을 하고 있고, 중앙값인 50을 기준으로 대칭인 것처럼 보인다. 따라서 위 분석에는 문제가 없는 듯하다.
가정이 성립하지 않는다면?
우리가 궁금하고, 문제가 되는 상황은 가정이 성립하지 않는 상황이다. 실제로 이럴 때가 많기 때문이다. 신경심리검사 2 원점수인 NP2_raw데이터의 분포를 확인해보자.
1) 분석(A) > 기술통계량(E) > 데이터 탐색 (E)
2) 분석하고자 하는 변수인 NP2_raw을 "종속변수"에 넣어준다. 그 뒤 "도표(T)..."를 선택한다.
3) "히스토그램(H)", "검정과 함께 정규성 도표(O)" 체크박스를 클릭하고 "계속(C)"를 누르고, 돌아가 "확인"을 누른다.
결과
분포는 오른쪽으로 꼬리가 긴 right skewed 모양을 하고 있으며, 중위수는 1.004532153이다. 이 값에 대해 일표본 윌콕슨 부호 순위 검정을 해보자. 보통 우리는 "정규성을 만족하지 않으므로 윌콕슨 부호 순위 검정으로 모집단의 중위수를 검정하자."라고 생각한다. 만약 이런 논리에 따른다면 1.004532135에 대해 일표본 윌콕슨 부호 순위 검정을 시행한다면 당연히 유의하지 않은 결과가 나와야 한다.
1) 분석(A) > 비모수검정(N) > 일표본(O)
2) 이때 나오는 창의 첫 페이지인 "목적"은 건들지 않는다.
3) "필드"를 누르고 분석하고자 하는 NP2_raw를 오른쪽으로 넘긴다.
4) "사용자 정의에 의한 검정(T)"를 누르고 "평균과 가설값 비교(Wilcoxon 부호 순위 검정)"을 체크하고, 가설 중위수에는 검정하고자 하는 값인 1.004532135를 적는다. 그리고 "실행"을 누른다.
결과
p-value가 0.001보다도 작다는 결론에 도달한다. 참 중앙값으로 "NP2_raw의 중앙값이 참 중앙값과 같니?"라고 물어보았는데 "아니요"라고 대답하는 참사가 벌어진다. 이렇게 '잘못된' 결론을 내는 이유는 사실 우리가 윌콕슨 부호 순위 검정에 대해 '잘못된' 지식이자 선입견을 갖고 있기 때문이다. 윌콕슨 부호 순위 검정을 좌우 대칭인 데이터에만 사용할 수 있고, 그럴 때에만 중위수에 대한 검정을 의미한다. 이럴 때 선택할 수 있는 옵션이 많지 않긴 한데, 그중 하나는 변환 (transformation)이다. 이런 right skewed data는 로그 변환 (log-transformation)을 하면 대칭성이 갖추어지는 경우가 많다. 로그 변환을 하고 히스토그램을 그려 분포를 확인해보자. (로그 변환에 관한 글은 다음 링크에서 확인할 수 있다. 2022.11.30 - [통계 프로그램 사용 방법/SPSS] - [SPSS] 변수 계산 (산술 연산))
2) 분석하고자 하는 변수인 NP2_raw을 "종속변수"에 넣어준다. 그 뒤 "도표(T)..."를 선택한다.
3) "히스토그램(H)", "검정과 함께 정규성 도표(O)" 체크박스를 클릭하고 "계속(C)"를 누르고, 돌아가 "확인"을 누른다.
결과
대칭성을 띠는 듯 하지만 정규성을 따르지는 않는다. 따라서 일표본 윌콕슨 부호 순위 검정을 시행하도록 한다.
1) 분석(A) > 비모수검정(N) > 일표본(O)
2) 이때 나오는 창의 첫 페이지인 "목적"은 건들지 않는다.
3) "필드"를 누르고 분석하고자 하는 NP2_raw를 오른쪽으로 넘긴다.
4) "사용자 정의에 의한 검정(T)"를 누르고 "평균과 가설값 비교(Wilcoxon 부호 순위 검정)"을 체크하고, 가설 중위수에는 검정하고자 하는 값을 넣는데, 원래 중위수가 1.004532135이었으므로 로그 변환을 하면 0과 가까우니 0을 적는다. 그리고 "실행"을 누른다.
결과
이제야 유의하지 않다는 결론을 맞게 잘 내준다.
그 외 다룰 내용 1: 정규성을 안 따르면 꼭 Wilcoxon test를 시행해야 하나?
T-test의 전제조건은 표본 평균이 정규분포를 따르는 것이다. 그런데 중심 극한 정리에 따라 표본의 수 (우리 데이터에서는 n=1,000)가 커질수록 표본 평균의 분포는 정규성을 띠게 된다. 따라서 n수가 적당히 크기만 하면 표본이 정규분포를 따르는지와 관계없이 표본 평균은 정규성을 따른다고 할 수 있다. 이런 이유로 정규성과 관계없이 n수가 크기만 하면 t-test의 이용이 정당화되기도 한다. 합리적인 말이다.
그런데 왜 정규분포를 따르는지 왜 다들 확인할까? 만약 모분포가 정규분포였다면 표본도 정규분포를 따를 확률이 높아진다. 물론 꼭 그런 건 아니지만 말이다. 따라서 표본이 정규분포를 따른다면, 모분포를 정규분포를 따른다고 할 수 있고, 그렇다면 n수가 작아도 표본 평균의 분포도 정규성을 따르니 t-test를 사용해도 되는 안전한 환경이 구축된다. 이런 이유로 정규성을 확인하게 된다.
그 외 다룰 내용 2: Wilcoxon signed rank test는 짝지어진 자료의 검정방법 아닌가?
그런데, 여기에는 각각의 분포가 정규성을 따른다는 중요한 가정이 필요하다. 하지만 정규성을 따르지 않는다면 어떻게 해야 할까? 그럴 때 사용하는 것이 크루스칼-왈리스 검정(Kruskal-Wallis Test)이다. 이번 포스팅에서는 크루스칼 왈리스 검정과 그 사후 분석 (post hoc analysis)까지 알아보고자 한다.
Shapiro-Wilk normality test
data: dfne$NICOT_ADDICT
W = 0.88076, p-value < 2.2e-16
2) 과거 흡연자
Shapiro-Wilk normality test
data: dffo$NICOT_ADDICT
W = 0.87644, p-value = 1.122e-14
3) 현재 흡연자
Shapiro-Wilk normality test
data: dfcu$NICOT_ADDICT
W = 0.93791, p-value = 1.364e-08
N수가 2,000개 미만이므로 Shapiro-Wilk 통계량의 p-value를 보면 0.05 이하이며, Q-Q Plot은 대부분의 데이터가 선상에 있지 않고, 히스토그램에서도 정규성을 따르지 않는 것처럼 보인다. 따라서 일원 배치 분산 분석 (ANOVA)을 시행할 수 없고, 크루스칼 왈리스 분석 (Kruskal-Wallis Test)을 시행해야 한다.
지금까지 많은 분석을 할 때 분산이 같은지 여부 (등분산성)을 중요하게 여겼었다. 하지만 크루스칼 왈리스 분석은 분산이 달라도 어느 정도 괜찮다. 물론 분산에 아주 큰 차이가 나면 위험해지지만, 어느 정도는 robust 하므로 웬만하면 그냥 써도 괜찮다.
크루스칼 왈리스 (Kruskal Wallis) 검정 코드
kruskal.test(NICOT_ADDICT~SMOK, data=df)
kruskal.test(NICOT_ADDICT~SMOK, data=df) : 데이터 df의 SMOK에 따라 NICOT_ADDICT의 중앙값에 차이가 있는지 Kruskal-Wallis 검정을 시행하라.
결과
Kruskal-Wallis rank sum test
data: NICOT_ADDICT by SMOK
Kruskal-Wallis chi-squared = 560.98, df = 2, p-value < 2.2e-16
p-value가 0.05보다 작으므로 귀무가설을 기각하고 대립 가설을 채택한다. 그렇다면 여기에서 귀무가설 및 대립 가설은 무엇이었는가?
귀무가설: $H_0=$세 집단의 모집단 수준에서 중앙값은 "모두" 동일하다.
대립가설: $H_1=$ 세 집단의 모집단 수준에서 중앙값이 모두 동일한 것은 아니다.
우리는 대립 가설을 채택해야 하므로 "세 집단의 모집단 수준에서 중앙값이 모두 동일한 것은 아니다."라고 결론 내릴 것이다.
사후 분석
그런데, 세 집단의 중앙값이 모두 동일하지 않다는 말은 세 집단 중 두 개씩 골라 비교했을 때, 적어도 한 쌍에서는 차이가 난다는 것이다. 따라서 세 집단 중 두 개씩 골라 비교를 해보아야 하며, 이를 사후 분석 (post hoc analysis)이라고 한다. 세 집단에서 두 개씩 고르므로 가능한 경우의 수 $_3C_2=3$이다.
(1) 비흡연자 vs 과거 흡연자
(2) 비흡연자 vs 현재 흡연자
(3) 과거 흡연자 vs 현재 흡연자
비모수 검정인 Kruskal Wallis 검정의 사후 분석으로 쓰이는 방법은 다음 네 가지가 있다.
1) Pairwise Wilcoxon Test
2) Conover Test
3) Dunn Test
4) Nemenyi Test
코드
#1) Pairwise Wilcoxon Test
pairwise.wilcox.test(df$NICOT_ADDICT, df$SMOK, p.adjust.method="bonferroni", correct=FALSE)
#2) Conover Test
install.packages("DescTools")
library("DescTools")
ConoverTest(NICOT_ADDICT~SMOK, data=df)
#3) Dunn Test
DunnTest(NICOT_ADDICT~SMOK, data=df)
#4) Nemenyi Test
NemenyiTest(NICOT_ADDICT~SMOK, data=df)
ConoverTest(NICOT_ADDICT~SMOK, data=df) : df데이터의 SMOK에 따라 NICOT_ADDICT의 중앙값에 차이가 있는지 Conover test를 시행하라. Conover test, Dunn test, Nemenyi test모두 구조는 같다. 이 세 개의 함수는 "DescTools"라는 패키지에 포함되어 있으므로 설치를 먼저 해야 한다. 패키지 설치에 관한 내용은 다음 링크를 확인하길 바란다. 2022.08.05 - [통계 프로그램 사용 방법/R] - [R] 패키지 설치하기 - install.packages(), library()
각각의 pairwise로 robust rank order test를 시행하고, 각각의 p-value를 모아 보정을 하는 것이다.
p<-c(p1, p2, p3) : p1, p2, p3을 묶어 하나의 p에 저장하라 p.adjust(p, "bonferroni") : p에 있는 값들을 보정하되, bonferroni adjustment를 사용하라.
Bonferroni 외에도, holm, hochberg, hommel, BH, BY, fdr을 넣어 각각의 방법으로 보정을 시도할 수 있다.
어떤 사후 분석을 쓸 것인가
이 논의에 대해 정답이 따로 있는 것은 아니다. 적절한 방법을 사용하여 논문에 제시하면 되고, 어떤 것이 정답이라고 콕 집어 이야기할 수는 없다. 다만, 사후 분석 방법이 여러 가지가 있다는 것은 '사후 분석 방법에 따라 산출되는 결과가 달라질 수 있다.'는 것을 의미하고, 심지어는 '어떤 사후 분석 방법을 채택하냐에 따라 유의성 여부가 달라질 수도 있다.'는 것을 의미한다. 심지어, Kruskal-Wallis test에서는 유의한 결과가 나왔는데, 사후 분석을 해보니 유의한 차이를 보이는 경우가 없을 수도 있다. 따라서 어떤 사후 분석 방법 결과에 따른 결과인지 유의하여 해석할 필요가 있다.
왜 굳이 Kruskal Wallis test를 쓰는 것인가?
이 시점에서 이런 의문이 들 수 있다.
"각각의 그룹별로 평균을 비교하면 되지, 굳이 왜 Kruskal Wallis test라는 방법까지 사용하는 것인가?"
아주 논리적인 의문점이다. 하지만, 반드시 Kruskal-Wallis test를 사용해야 한다. 그 이유는 다음과 같다. 본 사례는 흡연 상태에 따른 조합 가능한 경우의 수가 3인데, 각각 유의성의 기준을 $\alpha=0.05$로 잡아보자. 이때 세 번의 비교에서 모두 귀무가설이 학문적인 진실인데(평균에 차이가 없다.), 세 번 모두 귀무가설을 선택할 확률은 $\left( 1-0.05 \right)^3 \approx 0.86$이다. (이해가 어려우면 p-value에 대한 설명 글을 읽고 오길 바란다. 2022.09.05 - [통계 이론] - [이론] p-value에 관한 고찰)
그런데, Kruskal-Wallis test의 귀무가설은 "모든 집단의 중앙값이 같다."이다. 따라서 모든 집단의 중앙값이 같은 것이 학문적 진실일 때, 적어도 한 번이라도 대립 가설을 선택하게 될 확률은 $1-0.86=0.14$가 된다. 즉, 유의성의 기준이 올라가게 되어, 덜 보수적인 결정을 내리게 되고, 다시 말하면 대립 가설을 잘 선택하는 쪽으로 편향되게 된다. 학문적으로는 '다중 검정 (multiple testing)을 시행하면 1종 오류가 발생할 확률이 증가하게 된다.'라고 표현한다.
따라서, 각각을 비교해보는 것이 아니라 한꺼번에 비교하는 Kruskal-Wallis test를 시행해야 함이 마땅하다.
여기에서 한 번 더 의문이 들 수 있다.
"사후 분석을 할 때에는 1종 오류가 발생하지 않는가?"
그렇다. 1종 오류가 발생할 확률이 있으므로, p-value의 기준을 더 엄격하게 (0.05보다 더 작게) 잡아야 한다. P-value를 보정하는 방법은 일반적으로는 Bonferroni, holm, hochberg, hommel, BH, BY, FDR 등이 있고, Kruskal-Wallis test의 사후 분석 방법으로는 특별히 Pairwise Wilcoxon Test, Conover Test, Dunn Test, Nemenyi Test을 사용하고 있다.
Shapiro-Wilk normality test
data: df_whtn$NP5_raw
W = 0.95661, p-value = 6.406e-11
2) 정상인
Shapiro-Wilk normality test
data: df_wohtn$NP5_raw
W = 0.9426, p-value = 4.956e-13
N수가 2,000개 미만이므로 Shapiro-Wilk 통계량의 p-value를 보면 0.05 이하이며, Q-Q Plot은 대부분의 데이터가 선상에 있지 않고, 히스토그램에서도 정규성을 따르지 않는 것처럼 보인다. 따라서 독립 표본 T검정 (Two-sample T-test)를 시행할 수 없고, 두 분포는 비슷하게 생기지 않았으므로 맨 휘트니 U 검정 (Mann-Whitney U test)도 시행할 수 없다. 그런데 다행히도 각각의 분포는 중앙을 기점으로 대칭성을 이루고 있는 듯하다. 따라서 Robust rank order test를 시행할 수 있겠다.
Robust Rank-Order Distributional Test
data: NP5_raw by HTN
z = -0.090036, p-value = 0.9283
alternative hypothesis: true location shift is not equal to 0
Robust Rank-Order Distributional Test
data: NP5_raw by HTN
z = -0.090036, p-value = 0.9283
alternative hypothesis: true location shift is not equal to 0
p-value를 보면 0.9283으로 0.05보다 크므로 귀무가설을 택한다. 즉, 중앙값에 차이가 없다고 결론 내린다.
로버스트 순위 순서 검정은 대칭성에 대한 전제만 요구할 뿐, 분포의 모양에 대해서는 어떠한 가정도 하지 않는다. 따라서 이런 경우에 아주 유용하다. 그러면 대칭성을 만족하지 못해 로버스트 순위 순서 검정도 하지 못할 때에는 어떻게 해야 할까? 적절한 변환을 통해 모양을 비슷하게 만들어주면 된다. 흔히 알려져 있는 대수적 변환 외에도 변환 함수를 커스터마이징하면 웬만한 분포는 대칭으로 만들 수 있다.
코드 정리
##작업 디렉토리 지정
setwd("C:/Users/user/Documents/Tistory_blog")
##데이터 불러오기
install.packages("readr")
library("readr")
df<-read_csv("Data.csv")
##데이터 나누기
df_whtn<-df[df$HTN==1,]
df_wohtn<-df[df$HTN==0,]
##정규성 검정
#고혈압 환자
# 1) Q-Q plot 그리기
qqnorm(df_whtn$NP5_raw)
qqline(df_whtn$NP5_raw)
# 2) 히스토그램 그리기
hist(df_whtn$NP5_raw, prob=TRUE, breaks=20) #자세하게 보기 위해 breaks=20을 추가함
NP5_rawrange<-seq(min(df_whtn$NP5_raw),max(df_whtn$NP5_raw),length=max(max(df_whtn$NP5_raw)-min(df_whtn$NP5_raw),100))
ND<-dnorm(NP5_rawrange,mean=mean(df_whtn$NP5_raw),sd=sd(df_whtn$NP5_raw))
lines(NP5_rawrange, ND, lwd=2)
# 3) Shapiro-Wilk test 시행하기
shapiro.test(df_whtn$NP5_raw)
# 1) Q-Q plot 그리기
qqnorm(df_wohtn$NP5_raw)
qqline(df_wohtn$NP5_raw)
# 2) 히스토그램 그리기
hist(df_wohtn$NP5_raw, prob=TRUE, breaks=20) #자세하게 보기 위해 breaks=20을 추가함
NP5_rawrange<-seq(min(df_wohtn$NP5_raw),max(df_wohtn$NP5_raw),length=max(max(df_wohtn$NP5_raw)-min(df_wohtn$NP5_raw),100))
ND<-dnorm(NP5_rawrange,mean=mean(df_wohtn$NP5_raw),sd=sd(df_wohtn$NP5_raw))
lines(NP5_rawrange, ND, lwd=2)
# 3) Shapiro-Wilk test 시행하기
shapiro.test(df_wohtn$NP5_raw)
#필요 패키지 설치
install.packages("trend")
library("trend")
#로버스트 순위 순서 검정
rrod.test(df_whtn$NP5_raw, df_wohtn$NP5_raw)
rrod.test(NP5_raw~HTN, data=df)
[R] 로버스트 순위 순서 검정 (비모수 독립 표본 중앙값 검정: Robust rank order test, Flinger-Pollicello test) 정복 완료!
