반응형

[R] 크루스칼 왈리스 검정 (비모수 일원 배치 분산 분석: Kruskal-Wallis Test) - kruskal.test(), pairwise.wilcox.test(), ConoverTest(), DunnTest(), NemenyiTest()

 

셋 이상의 분포에서 평균이 다른지 확인하는 방법은 ANOVA (ANalysis Of VAriance)였다. 

2022.11.23 - [모평균 검정/R] - [R] 일원 배치 분산 분석 (One-way ANOVA, ANalysis Of VAriance) - aov(), pairwise.t.test(), ScheffeTest(), scheffe.test(), duncan.test(), TukeyHSD()

2022.11.30 - [모평균 검정/R] - [R] 등분산 가정이 성립하지 않는 일원 배치 분산 분석 (Welch's ANOVA, Brown-Forsythe ANOVA) - oneway.test(), bf.test(), tamhaneT2Test(), dunnettT3Test(), gamesHowellTest()

그런데, 여기에는 각각의 분포가 정규성을 따른다는 중요한 가정이 필요하다. 하지만 정규성을 따르지 않는다면 어떻게 해야 할까? 그럴 때 사용하는 것이 크루스칼-왈리스 검정(Kruskal-Wallis Test)이다.  이번 포스팅에서는 크루스칼 왈리스 검정과 그 사후 분석 (post hoc analysis)까지 알아보고자 한다.

 

*실습용 데이터는 아래 링크를 클릭하면 다운로드할 수 있습니다.

2022.08.04 - [공지사항 및 소개] - 분석용 데이터 (update 22.12.01)

 

분석용 데이터 (update 22.12.01)

2022년 12월 01일 버전입니다. 변수는 계속하여 추가될 예정입니다. 다음 카테고리에 있는 글에서 이용된 데이터입니다. - 기술 통계 - 범주형 자료 분석 - 모평균 검정 - 반복 측정 자료 분석 - 통계

medistat.tistory.com

 

코드를 보여드리기에 앞서 워킹 디렉토리부터 지정하겠다.

워킹 디렉토리에 관한 설명은 다음 링크된 포스트에서 볼 수 있다.

2022.08.05 - [통계 프로그램 사용 방법/R] - [R] 작업 디렉토리 (Working Directory) 지정 - getwd(), setwd()

setwd("C:/Users/user/Documents/Tistory_blog")

 

데이터를 불러와 df에 객체로 저장하겠다.

데이터 불러오는 방법은 다음 링크에서 볼 수 있다.

2022.08.05 - [통계 프로그램 사용 방법/R] - [R] 데이터 불러오기 : EXCEL - read_excel(), read.xlsx()

2022.08.10 - [통계 프로그램 사용 방법/R] - [R] 데이터 저장하기 : CSV 파일 - write.csv(), write_csv()

2022.08.10 - [통계 프로그램 사용 방법/R] - [R] 데이터 불러오기 : SAS file (.sas7bdat) - read.sas7bdat(), read_sas()

install.packages("readr")
library("readr")
df<-read_csv("Data.csv")

 

목표:  흡연 상태(SMOK)에 따라 니코틴 중독 점수 (NICOT_ADDICT) 중앙값이 모집단 수준에서 서로 다르다고 말할 수 있는가?

 

 만약 NICOT_ADDICT의 분포가 정규성을 따른다면, ANOVA로 이를 확인할 수 있을 것이다. 따라서 정규성 여부를 먼저 확인해보자.

1) 정규성 - 코드

흡연 상태 별로 NICOT_ADDICT의 정규성을 확인하기 위해서는 다음 두 가지의 일을 해야 한다.

1) 흡연 상태에 따라 데이터를 나누기

2) 각각 정규성을 확인하기.

 

1)데이터 나누기

조건에 따라 데이터를 나누는 방법에 대한 설명은 다음 링크에서 볼 수 있다. 2022.11.10 - [통계 프로그램 사용 방법/R] - [R] 조건에 맞는 자료 추출하기

위 링크에서 확인할 수 있듯이 여러 가지 방법으로 나눌 수 있지만 indexing을 이용하여 나누도록 하겠다.

dfne<-df[df$SMOK==0,]
dffo<-df[df$SMOK==1,]
dfcu<-df[df$SMOK==2,]

dfne: 비흡연자의 데이터

dffo: 과거 흡연자의 데이터

dfcu: 현재 흡연자의 데이터

 

2) 정규성 확인

정규성 검정에 관한 내용은 다음 링크에서 확인할 수 있다.

2022.08.11 - [기술 통계/R] - [R] 정규성 검정 (1) : Q-Q plot - qqnorm(), qqline()

2022.08.11 - [기술 통계/R] - [R] 정규성 검정 (2) : 히스토그램 - hist(), dnorm()

2022.08.11 - [기술 통계/R] - [R] 정규성 검정 (3) : 정량적 검정 (Shapiro-Wilk, Kolmogorov-Smirnov) - shapiro.test(), ks.test()

2022.08.12 - [기술 통계/R] - [R] 정규성 검정 (4) : 정량적 검정 (Lilliefors test) - lillie.test()

2022.08.16 - [기술 통계/R] - [R] 고급 Q-Q Plot - Van der Waerden, Rankit, Tukey, Blom

 

#비흡연자
# 1) Q-Q plot 그리기
qqnorm(dfne$NICOT_ADDICT)
qqline(dfne$NICOT_ADDICT)
# 2) 히스토그램 그리기
hist(dfne$NICOT_ADDICT, prob=TRUE)
NICOT_ADDICTrange<-seq(min(dfne$NICOT_ADDICT),max(dfne$NICOT_ADDICT),length=max(max(dfne$NICOT_ADDICT)-min(dfne$NICOT_ADDICT),100))
ND<-dnorm(NICOT_ADDICTrange,mean=mean(dfne$NICOT_ADDICT),sd=sd(dfne$NICOT_ADDICT))
lines(NICOT_ADDICTrange, ND, lwd=2)
# 3) Shapiro-Wilk test 시행하기
shapiro.test(dfne$NICOT_ADDICT)

#과거 흡연자
# 1) Q-Q plot 그리기
qqnorm(dffo$NICOT_ADDICT)
qqline(dffo$NICOT_ADDICT)
# 2) 히스토그램 그리기
hist(dffo$NICOT_ADDICT, prob=TRUE)
NICOT_ADDICTrange<-seq(min(dffo$NICOT_ADDICT),max(dffo$NICOT_ADDICT),length=max(max(dffo$NICOT_ADDICT)-min(dffo$NICOT_ADDICT),100))
ND<-dnorm(NICOT_ADDICTrange,mean=mean(dffo$NICOT_ADDICT),sd=sd(dffo$NICOT_ADDICT))
lines(NICOT_ADDICTrange, ND, lwd=2)
# 3) Shapiro-Wilk test 시행하기
shapiro.test(dffo$NICOT_ADDICT)

#현재 흡연자
# 1) Q-Q plot 그리기
qqnorm(dfcu$NICOT_ADDICT)
qqline(dfcu$NICOT_ADDICT)
# 2) 히스토그램 그리기
hist(dfcu$NICOT_ADDICT, prob=TRUE)
NICOT_ADDICTrange<-seq(min(dfcu$NICOT_ADDICT),max(dfcu$NICOT_ADDICT),length=max(max(dfcu$NICOT_ADDICT)-min(dfcu$NICOT_ADDICT),100))
ND<-dnorm(NICOT_ADDICTrange,mean=mean(dfcu$NICOT_ADDICT),sd=sd(dfcu$NICOT_ADDICT))
lines(NICOT_ADDICTrange, ND, lwd=2)
# 3) Shapiro-Wilk test 시행하기
shapiro.test(dfcu$NICOT_ADDICT)

 

결과

1) 비흡연자

	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

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)

pairwise.wilcox.test(df$NICOT_ADDICT, df$SMOK, p.adjust.method="bonferroni", correct=FALSE) : df 데이터의 SMOK에 따라 df 데이터의 NICOT_ADDICT의 중앙값에 차이가 있는지 Wilcoxon rank sum test를 시행하는데, p-value의 보정은 Bonferroni 방법으로 시행하고, 연속성 수정은 하지 말아라. (연속성 수정에 관한 내용은 다음 링크를 확인하길 바란다.2022.08.30 - [통계 이론] - [이론] 연속성을 수정한 카이 제곱 검정 (Chi-squared test with Yates's correction for continuity))

 

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 comparisons using Wilcoxon rank sum test 

data:  df$NICOT_ADDICT and df$SMOK 

  0       1      
1 8.3e-08 -      
2 < 2e-16 < 2e-16

P value adjustment method: bonferroni



 Conover's test of multiple comparisons : holm  

    mean.rank.diff    pval    
1-0       91.06188 2.6e-10 ***
2-0      532.38681 < 2e-16 ***
2-1      441.32493 < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1



 Dunn's test of multiple comparisons using rank sums : holm  

    mean.rank.diff    pval    
1-0       91.06188 2.3e-05 ***
2-0      532.38681 < 2e-16 ***
2-1      441.32493 < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1



 Nemenyi's test of multiple comparisons for independent samples (tukey)  

    mean.rank.diff    pval    
1-0       91.06188 6.9e-05 ***
2-0      532.38681 < 2e-16 ***
2-1      441.32493 < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons using Wilcoxon rank sum test 

data:  df$NICOT_ADDICT and df$SMOK 

  0       1      
8.3e-08 -      
< 2e-16 < 2e-16

P value adjustment method: bonferroni

 

빨간색 글씨가 p-value를 의미하는데, 이는 다음과 같이 보면 된다.

  0 1
1 비흡연자 (0)와 과거 흡연자 (1) 비교  
2 비흡연자 (0)와 현재 흡연자 (2) 비교 과거 흡연자 (1)와 현재 흡연자 (2) 비교

모든 p-value가 0.05보다 작으므로 다음과 같이 결론 내릴 수 있다.

비흡연자 (0)와 과거 흡연자 (1)  사이에는 NICOT_ADDICT 중앙값에 차이가 있다.

비흡연자 (0)와 현재 흡연자 (2) 사이에는 NICOT_ADDICT 중앙값에 차이가 있다.

과거 흡연자 (1)와 현재 흡연자 (2)) 사이에는 NICOT_ADDICT 중앙값에 차이가 있다.




 Conover's test of multiple comparisons : holm  

    mean.rank.diff    pval    
1-0       91.06188 2.6e-10 ***
2-0      532.38681 < 2e-16 ***
2-1      441.32493 < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1



 Dunn's test of multiple comparisons using rank sums : holm  

    mean.rank.diff    pval    
1-0       91.06188 2.3e-05 ***
2-0      532.38681 < 2e-16 ***
2-1      441.32493 < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1



 Nemenyi's test of multiple comparisons for independent samples (tukey)  

    mean.rank.diff    pval    
1-0       91.06188 6.9e-05 ***
2-0      532.38681 < 2e-16 ***
2-1      441.32493 < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

 

Conover Test, Dunn Test, Nemenyi Test 결과를 보는 방법은 같다. 

1-0 비흡연자 (0)와 과거 흡연자 (1) 비교
2-0 비흡연자 (0)와 현재 흡연자 (2) 비교
1-2 과거 흡연자 (1)와 현재 흡연자 (2) 비교

모든 p-value가 0.05보다 작으므로 다음과 같이 결론 내릴 수 있다.

비흡연자 (0)와 과거 흡연자 (1)  사이에는 NICOT_ADDICT 중앙값에 차이가 있다.

비흡연자 (0)와 현재 흡연자 (2) 사이에는 NICOT_ADDICT 중앙값에 차이가 있다.

과거 흡연자 (1)와 현재 흡연자 (2)) 사이에는 NICOT_ADDICT 중앙값에 차이가 있다.

 

즉, 어떤 방법을 사용하더라도 유의미하게 나온다.

 

하지만, 위 분석 방법들은 서로 다른 모양의 분포, 서로 다른 분산에 대한 고려가 부족하므로, 필자는 여기에 아직 쓰이고 있지는 않지만 한 가지 방법을 추가로 제안한다. 바로, Robust rank order test의 p-value를 보정하여 쓰는 것이다. (Robust rank order test에 대한 내용은 다음 링크를 확인하길 바란다.2022.12.01 - [모평균 검정/R] - [R] 로버스트 순위 순서 검정 (비모수 독립 표본 중앙값 검정: Robust rank order test, Flinger-Pollicello test) - rrod.test())

 

rrod1<-rrod.test(dfne$NICOT_ADDICT, dffo$NICOT_ADDICT)
rrod2<-rrod.test(dfne$NICOT_ADDICT, dfcu$NICOT_ADDICT)
rrod3<-rrod.test(dffo$NICOT_ADDICT, dfcu$NICOT_ADDICT)
p1<-rrod1$p.value
p2<-rrod2$p.value
p3<-rrod3$p.value
p<-c(p1, p2, p3)
p.adjust(p, "bonferroni")
p.adjust(p, "fdr")

각각의 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을 사용하고 있다. 

 

코드 정리

##워킹 디렉토리 지정
setwd("C:/Users/user/Documents/Tistory_blog")


##데이터 불러오기
install.packages("readr")
library("readr")
df<-read_csv("Data.csv")

##데이터 나누기
dfne<-df[df$SMOK==0,]
dffo<-df[df$SMOK==1,]
dfcu<-df[df$SMOK==2,]


##정규성 검정
#비흡연자
# 1) Q-Q plot 그리기
qqnorm(dfne$NICOT_ADDICT)
qqline(dfne$NICOT_ADDICT)
# 2) 히스토그램 그리기
hist(dfne$NICOT_ADDICT, prob=TRUE)
NICOT_ADDICTrange<-seq(min(dfne$NICOT_ADDICT),max(dfne$NICOT_ADDICT),length=max(max(dfne$NICOT_ADDICT)-min(dfne$NICOT_ADDICT),100))
ND<-dnorm(NICOT_ADDICTrange,mean=mean(dfne$NICOT_ADDICT),sd=sd(dfne$NICOT_ADDICT))
lines(NICOT_ADDICTrange, ND, lwd=2)
# 3) Shapiro-Wilk test 시행하기
shapiro.test(dfne$NICOT_ADDICT)

#과거 흡연자
# 1) Q-Q plot 그리기
qqnorm(dffo$NICOT_ADDICT)
qqline(dffo$NICOT_ADDICT)
# 2) 히스토그램 그리기
hist(dffo$NICOT_ADDICT, prob=TRUE)
NICOT_ADDICTrange<-seq(min(dffo$NICOT_ADDICT),max(dffo$NICOT_ADDICT),length=max(max(dffo$NICOT_ADDICT)-min(dffo$NICOT_ADDICT),100))
ND<-dnorm(NICOT_ADDICTrange,mean=mean(dffo$NICOT_ADDICT),sd=sd(dffo$NICOT_ADDICT))
lines(NICOT_ADDICTrange, ND, lwd=2)
# 3) Shapiro-Wilk test 시행하기
shapiro.test(dffo$NICOT_ADDICT)

#현재 흡연자
# 1) Q-Q plot 그리기
qqnorm(dfcu$NICOT_ADDICT)
qqline(dfcu$NICOT_ADDICT)
# 2) 히스토그램 그리기
hist(dfcu$NICOT_ADDICT, prob=TRUE)
NICOT_ADDICTrange<-seq(min(dfcu$NICOT_ADDICT),max(dfcu$NICOT_ADDICT),length=max(max(dfcu$NICOT_ADDICT)-min(dfcu$NICOT_ADDICT),100))
ND<-dnorm(NICOT_ADDICTrange,mean=mean(dfcu$NICOT_ADDICT),sd=sd(dfcu$NICOT_ADDICT))
lines(NICOT_ADDICTrange, ND, lwd=2)
# 3) Shapiro-Wilk test 시행하기
shapiro.test(dfcu$NICOT_ADDICT)


##Kruskal Wallis Test
kruskal.test(NICOT_ADDICT~SMOK, data=df) 


##Post hoc analysis
#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)

#5) Pairwise robust rank order test
rrod1<-rrod.test(dfne$NICOT_ADDICT, dffo$NICOT_ADDICT)
rrod2<-rrod.test(dfne$NICOT_ADDICT, dfcu$NICOT_ADDICT)
rrod3<-rrod.test(dffo$NICOT_ADDICT, dfcu$NICOT_ADDICT)
p1<-rrod1$p.value
p2<-rrod2$p.value
p3<-rrod3$p.value
p<-c(p1, p2, p3)
p.adjust(p, "bonferroni")
p.adjust(p, "fdr")

 

[R] 크루스칼 왈리스 검정 (비모수 일원 배치 분산 분석: Kruskal-Wallis Test) 정복 완료!

작성일: 2022.12.01.

최종 수정일: 2022.12.01.

이용 프로그램: R 4.2.2

RStudio v2022.07.2

RStudio 2022.07.2+576 "Spotted Wakerobin" Release

운영체제: Windows 10, Mac OS 12.6.1

반응형
반응형

[R] 로버스트 순위 순서 검정 (비모수 독립 표본 중앙값 검정: Robust rank order test, Flinger-Pollicello test) - rrod.test()

 

 정규분포를 따르지 않는 경우 두 분포의 중앙값이 같은지 보기 위해 맨 휘트니 U 검정 (윌콕슨 순위 합 검정, Mann Whitney U test, Wilcoxon rank sum test)을 시행할 수 있다.(2022.12.01 - [모평균 검정/R] - [R] 윌콕슨 순위 합 검정, 맨 휘트니 U 검정 (비모수 독립 표본 중앙값 검정: Wilcoxon rank sum test, Mann-Whitney U test) - wilcox.test()) 하지만 윌콕슨 순위 합 검정에는 아주 중요하고 강력한 하나의 가정이 필요한데, 두 분포의 생김새가 같다는 것이다. 만약 같지 않다면 정확도가 떨어지게 된다. 그럴 때 사용할 수 있는 로버스트 순위 순서 검정에 대해 소개하고자 한다.

 

단, 이 분석 또한 가정이 아예 필요 없는 것은 아니다. 로버스트 순위 순서 검정은 윌콕슨 순위 합 검정보다는 약한 가정 하나만 요구한다. 바로 중앙값에 대한 대칭성이다.

 

*실습용 데이터는 아래 링크를 클릭하면 다운로드할 수 있습니다.

2022.08.04 - [공지사항 및 소개] - 분석용 데이터 (update 22.12.01)

 

분석용 데이터 (update 22.12.01)

2022년 12월 01일 버전입니다. 변수는 계속하여 추가될 예정입니다. 다음 카테고리에 있는 글에서 이용된 데이터입니다. - 기술 통계 - 범주형 자료 분석 - 모평균 검정 - 반복 측정 자료 분석 - 통계

medistat.tistory.com

 

코드를 보여드리기에 앞서 워킹 디렉토리부터 지정하겠다.

워킹 디렉토리에 관한 설명은 다음 링크된 포스트에서 볼 수 있다.

2022.08.05 - [통계 프로그램 사용 방법/R] - [R] 작업 디렉토리 (Working Directory) 지정 - getwd(), setwd()

setwd("C:/Users/user/Documents/Tistory_blog")

 

데이터를 불러와 df에 객체로 저장하겠다.

데이터 불러오는 방법은 다음 링크에서 볼 수 있다.

2022.08.05 - [통계 프로그램 사용 방법/R] - [R] 데이터 불러오기 : EXCEL - read_excel(), read.xlsx()

2022.08.10 - [통계 프로그램 사용 방법/R] - [R] 데이터 저장하기 : CSV 파일 - write.csv(), write_csv()

2022.08.10 - [통계 프로그램 사용 방법/R] - [R] 데이터 불러오기 : SAS file (.sas7bdat) - read.sas7bdat(), read_sas()

install.packages("readr")
library("readr")
df<-read_csv("Data.csv")

 

목표:  고혈압 여부에 따라 신경심리검사 5 원점수 (NP5_raw)의 중앙값이 모집단 수준에서 서로 다르다고 말할 수 있는가?

만약 GGT의 분포가 정규성을 따른다면, 독립 표본 T 검정으로 이를 확인할 수 있을 것이다. 따라서 정규성 여부를 먼저 확인해보자.

 

1) 정규성 - 코드

고혈압 별로 NP5_raw의 정규성을 확인하기 위해서는 다음 두 가지의 일을 해야 한다.

1) 고혈압 여부에 따라 데이터를 나누기

2) 각각 정규성을 확인하기.

 

1)데이터 나누기

조건에 따라 데이터를 나누는 방법에 대한 설명은 다음 링크에서 볼 수 있다. 2022.11.10 - [통계 프로그램 사용 방법/R] - [R] 조건에 맞는 자료 추출하기

위 링크에서 확인할 수 있듯이 여러 가지 방법으로 나눌 수 있지만 indexing을 이용하여 나누도록 하겠다.

df_whtn<-df[df$HTN==1,]
df_wohtn<-df[df$HTN==0,]

 

2) 정규성 확인

정규성 검정에 관한 내용은 다음 링크에서 확인할 수 있다.

2022.08.11 - [기술 통계/R] - [R] 정규성 검정 (1) : Q-Q plot - qqnorm(), qqline()

2022.08.11 - [기술 통계/R] - [R] 정규성 검정 (2) : 히스토그램 - hist(), dnorm()

2022.08.11 - [기술 통계/R] - [R] 정규성 검정 (3) : 정량적 검정 (Shapiro-Wilk, Kolmogorov-Smirnov) - shapiro.test(), ks.test()

2022.08.12 - [기술 통계/R] - [R] 정규성 검정 (4) : 정량적 검정 (Lilliefors test) - lillie.test()

2022.08.16 - [기술 통계/R] - [R] 고급 Q-Q Plot - Van der Waerden, Rankit, Tukey, Blom

#고혈압 환자
# 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)

 

결과

1) 고혈압 환자

	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 test) 코드

install.packages("trend")
library("trend")

rrod.test(df_whtn$NP5_raw, df_wohtn$NP5_raw)
rrod.test(NP5_raw~HTN, data=df)

로버스트 순위 검정을 시행하는 rrod.test()함수는 trend 패키지에 들어있으므로 설치한다. 패키지 설치에 관한 내용은 다음 링크에서 확인할 수 있다. 2022.08.05 - [통계 프로그램 사용 방법/R] - [R] 패키지 설치하기 - install.packages(), library()

두 개의 코드는 같은 p-value를 제시한다.

아래 코드를 기준으로 설명하면 다음과 같다.

rrod.test(NP5_raw~HTN, data=df) : HTN에 따라 NP5_raw의 중앙값을 검정하되, 데이터는 df를 사용하라.

 

 

결과

	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) 정복 완료!

 

작성일: 2022.12.01.

최종 수정일: 2022.12.19.

이용 프로그램: R 4.2.2

RStudio v2022.07.2

RStudio 2022.07.2+576 "Spotted Wakerobin" Release

운영체제: Windows 10, Mac OS 12.6.1

반응형
반응형

[R] 윌콕슨 순위 합 검정, 맨 휘트니 U 검정 (비모수 독립 표본 중앙값 검정: Wilcoxon rank sum test, Mann-Whitney U test) - wilcox.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)에 대해 알아볼 것이다.

 

*실습용 데이터는 아래 링크를 클릭하면 다운로드할 수 있습니다.

2022.08.04 - [공지사항 및 소개] - 분석용 데이터 (update 22.11.30)

 

분석용 데이터 (update 22.11.28)

2022년 11월 28일 버전입니다. 변수는 계속하여 추가될 예정입니다. 다음 카테고리에 있는 글에서 이용된 데이터입니다. - 기술 통계 - 범주형 자료 분석 - 모평균 검정 - 반복 측정 자료 분석 - 통계

medistat.tistory.com

 

코드를 보여드리기에 앞서 워킹 디렉토리부터 지정하겠다.

워킹 디렉토리에 관한 설명은 다음 링크된 포스트에서 볼 수 있다.

2022.08.05 - [통계 프로그램 사용 방법/R] - [R] 작업 디렉토리 (Working Directory) 지정 - getwd(), setwd()

setwd("C:/Users/user/Documents/Tistory_blog")

 

데이터를 불러와 df에 객체로 저장하겠다.

데이터 불러오는 방법은 다음 링크에서 볼 수 있다.

2022.08.05 - [통계 프로그램 사용 방법/R] - [R] 데이터 불러오기 : EXCEL - read_excel(), read.xlsx()

2022.08.10 - [통계 프로그램 사용 방법/R] - [R] 데이터 저장하기 : CSV 파일 - write.csv(), write_csv()

2022.08.10 - [통계 프로그램 사용 방법/R] - [R] 데이터 불러오기 : SAS file (.sas7bdat) - read.sas7bdat(), read_sas()

install.packages("readr")
library("readr")
df<-read_csv("Data.csv")

 

목표:  음주 여부에 따라 Gamma-glutamyl transferase(GGT)의 중앙값이 모집단 수준에서 서로 다르다고 말할 수 있는가?

 

 만약 GGT의 분포가 정규성을 따른다면, 독립 표본 T 검정으로 이를 확인할 수 있을 것이다. 따라서 정규성 여부를 먼저 확인해보자.

 

1) 정규성 - 코드

음주 여부 별로 GGT의 정규성을 확인하기 위해서는 다음 두 가지의 일을 해야 한다.

1) 음주 여부에 따라 데이터를 나누기

2) 각각 정규성을 확인하기.

 

1)데이터 나누기

조건에 따라 데이터를 나누는 방법에 대한 설명은 다음 링크에서 볼 수 있다. 2022.11.10 - [통계 프로그램 사용 방법/R] - [R] 조건에 맞는 자료 추출하기

위 링크에서 확인할 수 있듯이 여러 가지 방법으로 나눌 수 있지만 indexing을 이용하여 나누도록 하겠다.

dfal<-df[df$ALCOHOL==1,]
dfno<-df[df$ALCOHOL==0,]

dfal: 음주자의 데이터

dfno: 비음주자의 데이터

 

2) 정규성 확인

정규성 검정에 관한 내용은 다음 링크에서 확인할 수 있다.

2022.08.11 - [기술 통계/R] - [R] 정규성 검정 (1) : Q-Q plot - qqnorm(), qqline()

2022.08.11 - [기술 통계/R] - [R] 정규성 검정 (2) : 히스토그램 - hist(), dnorm()

2022.08.11 - [기술 통계/R] - [R] 정규성 검정 (3) : 정량적 검정 (Shapiro-Wilk, Kolmogorov-Smirnov) - shapiro.test(), ks.test()

2022.08.12 - [기술 통계/R] - [R] 정규성 검정 (4) : 정량적 검정 (Lilliefors test) - lillie.test()

2022.08.16 - [기술 통계/R] - [R] 고급 Q-Q Plot - Van der Waerden, Rankit, Tukey, Blom

 

#음주자
# 1) Q-Q plot 그리기
qqnorm(dfal$GGT)
qqline(dfal$GGT)
# 2) 히스토그램 그리기
hist(dfal$GGT, prob=TRUE)
GGTrange<-seq(min(dfal$GGT),max(dfal$GGT),length=max(max(dfal$GGT)-min(dfal$GGT),100))
ND<-dnorm(GGTrange,mean=mean(dfal$GGT),sd=sd(dfal$GGT))
lines(GGTrange, ND, lwd=2)
# 3) Shapiro-Wilk test 시행하기
shapiro.test(dfal$GGT)

#비음주자
# 1) Q-Q plot 그리기
qqnorm(dfno$GGT)
qqline(dfno$GGT)
# 2) 히스토그램 그리기
hist(dfno$GGT, prob=TRUE)
GGTrange<-seq(min(dfno$GGT),max(dfno$GGT),length=max(max(dfno$GGT)-min(dfno$GGT),100))
ND<-dnorm(GGTrange,mean=mean(dfno$GGT),sd=sd(dfno$GGT))
lines(GGTrange, ND, lwd=2)
# 3) Shapiro-Wilk test 시행하기
shapiro.test(dfno$GGT)

 

결과

1) 음주자

 

	Shapiro-Wilk normality test

data:  dfal$GGT
W = 0.91166, p-value < 2.2e-16

 

2) 비음주자

	Shapiro-Wilk normality test

data:  dfno$GGT
W = 0.90691, p-value = 3.637e-15

 

 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) 코드

wilcox.test(dfal$GGT,dfno$GGT,correct=FALSE)

wilcox.test(dfal$GGT, dfno$GGT, correct=FALSE) : 맨 휘트니 U 검정을 시행하는데, dfal데이터의 GGT와 dfno데이터의 GGT를 비교하라. 연속성 수정은 하지 마라. (연속성 수정에 관한 내용은 다음 링크를 확인하길 바란다.2022.08.30 - [통계 이론] - [이론] 연속성을 수정한 카이 제곱 검정 (Chi-squared test with Yates's correction for continuity))

 

결과

	Wilcoxon rank sum test

data:  dfal$GGT and dfno$GGT
W = 200084, p-value < 2.2e-16
alternative hypothesis: true location shift is not equal to 0

Wilcoxon rank sum test

data:  dfal$GGT and dfno$GGT

W = 200084, p-value < 2.2e-16

alternative hypothesis: true location shift is not equal to 0

 

p-value가 매우 작으므로 유의한 결과임을 알 수 있다. 따라서 귀무가설을 기각하고 대립 가설을 택해야한다. 여기에서 귀무가설과 대립가설은 무엇일까?

 

귀무가설 $H_0=$ 두 분포는 동일하다.

대립가설 $H_1=$ 두 분포는 평균만 다를 뿐, 분포는 똑같이 생겼다.

즉, 대립 가설이 의미하는 것은 두 분포는 평행이동 관계에 있다는 것이다. 그래야 중앙값이 달라진 것인지 확인할 수 있다. 

평행이동 관계

이 말은, 맨 휘트니 U 검정에도 필요한 가정이 있다는 말이다. 두 분포는 동일하게 생겼어야 한다.

조금 어렵게 이야기하면, 두 분포의 first moment는 다르지만, 그보다 고차원의 central moment는 같아야 한다.

만약 생김새가 다르면 맨 휘트니 U 검정을 사용하면 안 되고, Robuts rank order test 등을 시행해야 한다. (물론 현실에서는 이런 내용을 무시한 채 '정규성을 따르지 않으니 맨 휘트니 U 검정을 사용한다'는 연구자가 대다수다.) 이에 관한 내용은 다음 링크를 확인하길 바란다. 2022.12.01 - [모평균 검정/R] - [R] 로버스트 순위 순서 검정 (비모수 독립 표본 중앙값 검정: Robust rank order test, Flinger-Pollicello test) - rrod.test()

 

 그런데 위 히스토그램을 보면 음주자와 비음주자의 GGT분포는 right skewed 된 분포로 어느 정도 비슷하게 생겼다. 모든 central moment를 조사하는 것은 현실적으로 힘들지만, 2nd cental moment인 분산의 차이를 검정해보면 다음과 같음을 알 수 있다. 등분산성 검정 방법은 다음 링크를 확인하길 바란다.

2022.11.12 - [모평균 검정/R] - [R] 독립 표본 T검정 (Independent samples T-test) - t.test(), var.test(), levene.test()

2022.11.20 - [모평균 검정/R] - [R] 등분산성 검정 (Homogeneity of variance) - levene.test(), bartlett.test()

 

코드

var.test(dfal$GGT,dfno$GGT)

 

결과

	F test to compare two variances

data:  dfal$GGT and dfno$GGT
F = 1.0248, num df = 589, denom df = 409, p-value = 0.7922
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
 0.8559314 1.2234543
sample estimates:
ratio of variances 
          1.024817

p-value는 0.7922로 분산에 차이가 없다고 결론 내리는 것이 합당하다.

따라서 두 분포의 모양이 같다고 생각하면 본 맨 휘트니 U 검정으로 음주자와 비음주자의 GGT값에는 차이가 있다고 결론 내릴 수 있다.

 

 

만약, 분포가 다르게 생기면 어떤 결과가 초래되길래 그렇게 강조하는 걸까?

데이터의 NP3_raw와 NP4_raw의 분포를 확인하고, 분산, 중앙값을 확인해보자

코드

hist(df$NP3_raw, breaks=20)
hist(df$NP4_raw, breaks=20)
median(df$NP3_raw)
median(df$NP4_raw)
var(df$NP3_raw)
var(df$NP4_raw)
var.test(df$NP3_raw, df$NP4_raw)
wilcox.test(df$NP3_raw, df$NP4_raw, correct=FALSE)

 

결과

> median(df$NP3_raw)
[1] 0.9968156
> median(df$NP4_raw)
[1] 0.9968156

> var(df$NP3_raw)
[1] 0.2237356
> var(df$NP4_raw)
[1] 0.2237356

> var.test(df$NP3_raw, df$NP4_raw)

	F test to compare two variances

data:  df$NP3_raw and df$NP4_raw
F = 1, num df = 999, denom df = 999, p-value = 1
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
 0.8832987 1.1321198
sample estimates:
ratio of variances 
                 1 
                 
                 
> wilcox.test(df$NP3_raw, df$NP4_raw, correct=FALSE)

	Wilcoxon rank sum test

data:  df$NP3_raw and df$NP4_raw
W = 568902, p-value = 9.515e-08
alternative hypothesis: true location shift is not equal to 0

NP3_raw는 right skewed, NP4_raw는 left skewed로 서루 매우 다르게 생겼지만, 분산, 중앙값이 완벽히 일치하며, 등분산성 검정에서도 분산에 차이가 없다고 나온다. 하지만 맨 휘트니 U 검정 결과 p-value는 매우 작게 나온다. 만약 맨 휘트니 U 검정이 "중앙값이 같은지 확인하는 검정"이라고만 알고 있다면 중앙값이 완전히 같은 이런 경우에도 "중앙값에 차이가 있다."라고 결론 내릴 것이다. 하지만 이 경우 귀무가설을 기각하고 "두 분포는 다르게 생겼다."라고 결론 내려야 옳은 결론에 다다른다. 

 

코드 정리

##워킹 디렉토리 지정
setwd("C:/Users/user/Documents/Tistory_blog")


##데이터 불러오기
install.packages("readr")
library("readr")
df<-read_csv("Data.csv")


##데이터 나누기
dfal<-df[df$ALCOHOL==1,]
dfno<-df[df$ALCOHOL==0,]


##정규성 검정
#음주자
# 1) Q-Q plot 그리기
qqnorm(dfal$GGT)
qqline(dfal$GGT)
# 2) 히스토그램 그리기
hist(dfal$GGT, prob=TRUE)
GGTrange<-seq(min(dfal$GGT),max(dfal$GGT),length=max(max(dfal$GGT)-min(dfal$GGT),100))
ND<-dnorm(GGTrange,mean=mean(dfal$GGT),sd=sd(dfal$GGT))
lines(GGTrange, ND, lwd=2)
# 3) Shapiro-Wilk test 시행하기
shapiro.test(dfal$GGT)

#비음주자
# 1) Q-Q plot 그리기
qqnorm(dfno$GGT)
qqline(dfno$GGT)
# 2) 히스토그램 그리기
hist(dfno$GGT, prob=TRUE)
GGTrange<-seq(min(dfno$GGT),max(dfno$GGT),length=max(max(dfno$GGT)-min(dfno$GGT),100))
ND<-dnorm(GGTrange,mean=mean(dfno$GGT),sd=sd(dfno$GGT))
lines(GGTrange, ND, lwd=2)
# 3) Shapiro-Wilk test 시행하기
shapiro.test(dfno$GGT)


##맨 휘트니 U 검정
wilcox.test(dfal$GGT,dfno$GGT,correct=FALSE)


##등분산성 검정
var.test(dfal$GGT,dfno$GGT)


##모양이 다른 경우
hist(df$NP3_raw, breaks=20)
hist(df$NP4_raw, breaks=20)
median(df$NP3_raw)
median(df$NP4_raw)
var(df$NP3_raw)
var(df$NP4_raw)
var.test(df$NP3_raw, df$NP4_raw)
wilcox.test(df$NP3_raw, df$NP4_raw, correct=FALSE)

 

[R] 윌콕슨 순위 합 검정, 맨 휘트니 U 검정 (비모수 독립 표본 중앙값 검정 정복 완료!

 

작성일: 2022.12.01.

최종 수정일: 2022.12.01.

이용 프로그램: R 4.2.2

RStudio v2022.07.2

RStudio 2022.07.2+576 "Spotted Wakerobin" Release

운영체제: Windows 10, Mac OS 12.6.1

반응형
반응형

[R] 등분산 가정이 성립하지 않는 일원 배치 분산 분석 (Welch's ANOVA, Brown-Forsythe ANOVA) - oneway.test(), bf.test(), tamhaneT2Test(), dunnettT3Test(), gamesHowellTest()

 

이전 포스팅에서 흡연 상태를 "1) 비흡연자", "2) 과거 흡연자", "3) 현재 흡연자"로 나누었고, 각 그룹의 폐기능 검사 중 하나인 FVC (Forced Vital Capacity)의 평균에 차이가 있는지 알아보고자 하여 일원 배치 분산 분석 (ANOVA)를 시행하였다.2022.11.23 - [모평균 검정/R] - [R] 일원 배치 분산 분석 (One-way ANOVA, ANalysis Of VAriance) - aov(), pairwise.t.test(), ScheffeTest(), scheffe.test(), duncan.test(), TukeyHSD()

 

ANOVA의 전제 조건은 두 가지였다.

1) 정규성

2) 등분산성

 

만약 '1) 정규성'은 만족하지만 '2) 등분산성'이 성립하지 않는 경우에는 어떻게 해야 할까?

당연히, 이전 포스팅의 일반적인 ANOVA는 실시할 수 없고, Welch's ANOVA 혹은 Brown-Forsythe ANOVA를 시행해야 한다.

 - 이 Brown-Forsythe ANOVA는 Brown-Forsythe의 등분산성 검정과는 다른 것이다. (아래에 각각의 코드가 등장한다.)

 

이번 포스팅에서는 이 Welch's ANOVA와 Brown-Forsythe ANOVA를 다뤄보고자 한다.

 

*실습용 데이터는 아래 링크를 클릭하면 다운로드할 수 있습니다.

2022.08.04 - [공지사항 및 소개] - 분석용 데이터 (update 22.11.28)

 

분석용 데이터 (update 22.11.28)

2022년 11월 28일 버전입니다. 변수는 계속하여 추가될 예정입니다. 다음 카테고리에 있는 글에서 이용된 데이터입니다. - 기술 통계 - 범주형 자료 분석 - 모평균 검정 - 반복 측정 자료 분석 - 통계

medistat.tistory.com

 

코드를 보여드리기에 앞서 워킹 디렉토리부터 지정하겠다.

워킹 디렉토리에 관한 설명은 다음 링크된 포스트에서 볼 수 있다.

2022.08.05 - [통계 프로그램 사용 방법/R] - [R] 작업 디렉토리 (Working Directory) 지정 - getwd(), setwd()

setwd("C:/Users/user/Documents/Tistory_blog")

 

데이터를 불러와 df에 객체로 저장하겠다.

데이터 불러오는 방법은 다음 링크에서 볼 수 있다.

2022.08.05 - [통계 프로그램 사용 방법/R] - [R] 데이터 불러오기 : EXCEL - read_excel(), read.xlsx()

2022.08.10 - [통계 프로그램 사용 방법/R] - [R] 데이터 저장하기 : CSV 파일 - write.csv(), write_csv()

2022.08.10 - [통계 프로그램 사용 방법/R] - [R] 데이터 불러오기 : SAS file (.sas7bdat) - read.sas7bdat(), read_sas()

install.packages("readr")
library("readr")
df<-read_csv("Data.csv")

 

목표:  거주 지역 단위에 따라 삶의 질의 평균이 모집단 수준에서 서로 다르다고 말할 수 있는가?

 

전제조건 (정규성)

 Welch's ANOVA 및 Brown-Forsythe ANOVA에는 하나의 전제조건이 필요하다. 본 포스팅의 예시에 맞추어 설명하면 다음과 같다.

  - 정규성: 대도시, 중소도시, 시골 거주자 별로 삶의 질의 분포를 보았을 때 각각의 분포는 정규성을 따른다. 

그리고, 분산은 같지 않아도 된다.

 

1) 정규성 - 코드

각 거주 지역 단위 별로 삶의 질 분포의 정규성을 확인하기 위해서는 다음 두 가지의 일을 해야 한다.

1) 거주 지역 단위 (대도시, 중소도시, 시골 거주자)에 따라 데이터를 나누기

2) 각각 정규성을 확인하기.

 

1)데이터 나누기

조건에 따라 데이터를 나누는 방법에 대한 설명은 다음 링크에서 볼 수 있다. 2022.11.10 - [통계 프로그램 사용 방법/R] - [R] 조건에 맞는 자료 추출하기

위 링크에서 확인할 수 있듯이 여러 가지 방법으로 나눌 수 있지만 indexing을 이용하여 나누도록 하겠다.

df_big<-df[df$RESID==0,]
df_med<-df[df$RESID==1,]
df_sma<-df[df$RESID==2,]

df_big : 대도시 지역 거주자의 데이터
df_med : 중소도시 지역 거주자의 데이터
df_sma : 시골 지역 거주자의 데이터

 

2) 정규성 확인

정규성 검정에 관한 내용은 다음 링크에서 확인할 수 있다.

2022.08.11 - [기술 통계/R] - [R] 정규성 검정 (1) : Q-Q plot - qqnorm(), qqline()

2022.08.11 - [기술 통계/R] - [R] 정규성 검정 (2) : 히스토그램 - hist(), dnorm()

2022.08.11 - [기술 통계/R] - [R] 정규성 검정 (3) : 정량적 검정 (Shapiro-Wilk, Kolmogorov-Smirnov) - shapiro.test(), ks.test()

2022.08.12 - [기술 통계/R] - [R] 정규성 검정 (4) : 정량적 검정 (Lilliefors test) - lillie.test()

2022.08.16 - [기술 통계/R] - [R] 고급 Q-Q Plot - Van der Waerden, Rankit, Tukey, Blom

#대도시 지역 거주자
# 1) Q-Q plot 그리기
qqnorm(df_big$QOL)
qqline(df_big$QOL)
# 2) 히스토그램 그리기
hist(df_big$QOL, prob=TRUE)
QOLrange<-seq(min(df_big$QOL),max(df_big$QOL),length=max(max(df_big$QOL)-min(df_big$QOL),100))
ND<-dnorm(QOLrange,mean=mean(df_big$QOL),sd=sd(df_big$QOL))
lines(QOLrange, ND, lwd=2)
# 3) Shapiro-Wilk test 시행하기
shapiro.test(df_big$QOL)

#중소도시 지역 거주자
# 1) Q-Q plot 그리기
qqnorm(df_med$QOL)
qqline(df_med$QOL)
# 2) 히스토그램 그리기
hist(df_med$QOL, prob=TRUE)
QOLrange<-seq(min(df_med$QOL),max(df_med$QOL),length=max(max(df_med$QOL)-min(df_med$QOL),100))
ND<-dnorm(QOLrange,mean=mean(df_med$QOL),sd=sd(df_med$QOL))
lines(QOLrange, ND, lwd=2)
# 3) Shapiro-Wilk test 시행하기
shapiro.test(df_med$QOL)

#시골지역 거주자
# 1) Q-Q plot 그리기
qqnorm(df_sma$QOL)
qqline(df_sma$QOL)
# 2) 히스토그램 그리기
hist(df_sma$QOL, prob=TRUE)
QOLrange<-seq(min(df_sma$QOL),max(df_sma$QOL),length=max(max(df_sma$QOL)-min(df_sma$QOL),100))
ND<-dnorm(QOLrange,mean=mean(df_sma$QOL),sd=sd(df_sma$QOL))
lines(QOLrange, ND, lwd=2)
# 3) Shapiro-Wilk test 시행하기
shapiro.test(df_sma$QOL)

 

결과

1) 대도시 거주자 (RESID=0)

	Shapiro-Wilk normality test

data:  df_big$QOL
W = 0.99468, p-value = 0.3259

 표본의 수가 2,000이 되지 않으므로 Shapiro-Wilk 통계량의 p-value를 보아야 하고, 이는 0.05보다 크다. Q-Q plot, 히스토그램에서 정규성을 띠지 않는다고 할만한 근거가 없으므로 정규성을 따른다고 결론을 짓는다.

 

 

2) 중소도시 거주자 (RESID=1)

	Shapiro-Wilk normality test

data:  df_med$QOL
W = 0.99583, p-value = 0.518

 표본의 수가 2,000이 되지 않으므로 Shapiro-Wilk 통계량의 p-value를 보아야 하고, 이는 0.05보다 크다. Q-Q plot, 히스토그램에서 정규성을 띠지 않는다고 할만한 근거가 없으므로 정규성을 따른다고 결론을 짓는다.

 

3) 시골 거주자 (RESID=2)

	Shapiro-Wilk normality test

data:  df_sma$QOL
W = 0.9981, p-value = 0.9689

 표본의 수가 2,000이 되지 않으므로 Shapiro-Wilk 통계량의 p-value를 보아야 하고, 이는 0.05보다 크다. Q-Q plot, 히스토그램에서 정규성을 띠지 않는다고 할만한 근거가 없으므로 정규성을 따른다고 결론을 짓는다.

 

따라서 정규성 전제는 따른다고 할 수 있다.

 

2) 등분산성

분산이 같은지 확인해보도록 하겠다.

등분산성 검정에 관한 내용은 다음 링크에서 확인할 수 있다.

2022.11.12 - [모평균 검정/R] - [R] 독립 표본 T검정 (Independent samples T-test) - t.test(), var.test(), levene.test()

2022.11.20 - [모평균 검정/R] - [R] 등분산성 검정 (Homogeneity of variance) - levene.test(), bartlett.test()

 

많이들 사용하는 등분산 검정 다섯 가지 중 F test는 그룹이 2개일 때에만 사용 가능하고, O'brien은 R에서 시행 불가하므로, 나머지 3가지 방법으로 검정해본다.

1) F test (그룹이 2개일 때에만 사용 가능)

2) Levene의 등분산 검정

3) O'Brien의 등분산 검정(R에서 검정 불가)

4) Brown and Forsythe의 등분산 검정

5) Bartlett의 등분산 검정

 

코드

#Levene의 등분산 검정
install.packages("lawstat")
library("lawstat")
levene.test(df$QOL, df$RESID, location="mean")

#Brown&Forsythe
install.packages("lawstat")
library("lawstat")
levene.test(df$QOL, df$RESID)

#Bartlett
bartlett.test(df$QOL, df$RESID)

 

결과

	Classical Levene's test based on the absolute deviations from the mean (
	none not applied because the location is not set to median )

data:  df$QOL
Test Statistic = 114.73, p-value < 2.2e-16



	Modified robust Brown-Forsythe Levene-type test based on the absolute
	deviations from the median

data:  df$QOL
Test Statistic = 114.64, p-value < 2.2e-16



	Bartlett test of homogeneity of variances

data:  df$QOL and df$RESID
Bartlett's K-squared = 326.64, df = 2, p-value < 2.2e-16

모두 p-value<0.05로 분산에 차이가 있다고 할 수 있다. 따라서 일반적인 ANOVA가 아니라 Welch's ANOVA 및 Brown-Forsythe ANOVA를 시행해야 한다.

 

Welch's 일원 배치 분산 분석 (Welch's One way ANOVA) 코드

oneway.test(QOL~RESID, data=df, var.equal=FALSE)

oneway.test(QOL~RESID, data=df, var.equal=FALSE) : One-way ANOVA를 시행한다 (oneway.test). RESID에 따른 QOL의 평균을 비교하라. 데이터는 df를 사용하고, 분산은 같지 않다 (var.equal=FALSE). 

 

결과

	One-way analysis of means

data:  QOL and RESID
F = 2398.8, num df = 2, denom df = 997, p-value < 2.2e-16

One-way analysis of means

data:  QOL and RESID

F = 2398.8, num df = 2, denom df = 997, p-value < 2.2e-16

 

P-value가 <0.0001이므로 귀무가설을 기각하고 대립 가설을 채택한다. 그렇다면 여기에서 귀무가설 및 대립 가설은 무엇이었는가?

 

 귀무가설:  $H_0=$세 집단의 모평균은 "모두" 동일하다.

 대립 가설:  $H_1=$ 세 집단의 모평균이 모두 동일한 것은 아니다.

 

우리는 대립 가설을 채택해야 하므로 "세 집단의 모평균이 모두 동일한 것은 아니다."라고 결론 내릴 것이다.

 

 

Brown-Forsythe's 일원 배치 분산 분석 (Brown-Forsythe's One way ANOVA) 코드

install.packages("onewaytests")
library("onewaytests")
df$RESID<-as.factor(df$RESID)
bf.test(QOL~RESID, data=df)

시행하기 위한 bf.test() 함수는 "onewaytests"라는 패키지에 있으므로 설치해야 한다. 패키지 설치 방법 안내: 2022.08.05 - [통계 프로그램 사용 방법/R] - [R] 패키지 설치하기 - install.packages(), library()

그리고 bf.test()함수는 독립변수로 문자형 혹은 범주형 변수를 요구한다. 따라서 범주형 변수로 변환한다. 변환 방법 안내: 2022.11.21 - [통계 프로그램 사용 방법/R] - [R] 변수의 유형 (타입, type) 확인 및 변경 - as.factor(), as.numeric(), as.character(), str()

bf.test(QOL~RESID, data=df) : Brown and Forsythe One way ANOVA를 시행한다 (bf.test). RESID에 따른 QOL 평균의 차이를 검정하라.

 

결과

  Brown-Forsythe Test (alpha = 0.05) 
------------------------------------------------------------- 
  data : QOL and RESID 

  statistic  : 2456.285 
  num df     : 2 
  denom df   : 683.7928 
  p.value    : 7.160917e-313 

  Result     : Difference is statistically significant. 
-------------------------------------------------------------

  Brown-Forsythe Test (alpha = 0.05) 

-------------------------------------------------------------    

data : QOL and RESID    

 

statistic  : 2456.285    

num df     : 2    

denom df   : 683.7928    

p.value    : 7.160917e-313    

 

Result     : Difference is statistically significant. 

------------------------------------------------------------- 

 

p-value가 0.05보다 작으므로 귀무가설을 기각하고 대립 가설을 채택한다. 그렇다면 여기에서 귀무가설 및 대립 가설은 무엇이었는가?

 

 귀무가설:  $H_0=$ 세 집단의 모평균은 "모두" 동일하다.

 대립 가설:  $H_1=$세 집단의 모평균이 모두 동일한 것은 아니다.

 

우리는 대립 가설을 채택해야 하므로 "세 집단의 모평균이 모두 동일한 것은 아니다."라고 결론 내릴 것이다.

 

즉, Welch, Brown-Forsythe 두 가지 방법 모두 평균에 유의미한 차이가 있다고 결론 내리고 있다.

 

그런데, 세 집단의 모평균이 모두 동일하지 않다는 말은 세 집단 중 두 개씩 골라 비교했을 때, 적어도 한 쌍에서는 차이가 난다는 것이다. 따라서 세 집단 중 두 개씩 골라 비교를 해보아야 하며, 이를 사후 분석 (post hoc analysis)이라고 한다. 세 집단에서 두 개씩 고르므로 가능한 경우의 수는$_3C_2=3$이다. 

 (1) 대도시 거주자 vs 중소도시 거주자

 (2) 대도시 거주자 vs 시골 거주자

 (3) 중소도시 거주자 vs 시골 거주자

 

사후 분석을 시행할 것이다.

이전 일반 ANOVA 포스팅(2022.11.23 - [모평균 검정/R] - [R] 일원 배치 분산 분석 (One-way ANOVA, ANalysis Of VAriance) - aov(), pairwise.t.test(), ScheffeTest(), scheffe.test(), duncan.test(), TukeyHSD())에서는 여러 사후 분석 방법이 있다고 소개했었다. 하지만 등분산이 가정되지 않는 상황에서 쓸 수 있는 사후 분석 방법은 현저히 적어진다. 많이들 쓰는 방법은 다음 네 가지다.

 

1) Tamhane's T2

2) Games-Howell

3) Dunnett's T3

4) Dunnett's C

 

하지만 마지막 Dunnett's C는 제한적으로 사용하므로 위 세 가지 방법에 대해서만 알아보도록 하겠다. 

 

사후 분석 코드

#1) Tamhane's T2
install.packages("PMCMRplus")
library("PMCMRplus")
tamhaneT2Test(QOL~RESID, data=df)

#2) Games-Howell
gamesHowellTest(QOL~RESID, data=df)

#3) Dunnett's T3
dunnettT3Test(QOL~RESID, data=df)

각 사후 분석에 쓰이는 함수들은 "PMCMRplus"라는 패키지에 있으므로 설치해야 한다. 설치 방법 링크:2022.08.05 - [통계 프로그램 사용 방법/R] - [R] 패키지 설치하기 - install.packages(), library()

tamhaneT2Test(QOL~RESID, data=df) : QOL의 평균을  RESID의 조합 가능한 경우의 수에 따라 비교(pairwise comparison)하되 방법은 Tamhane's T2 검정을 사용하라.

gamesHowellTest(QOL~RESID, data=df) : QOL의 평균을  RESID의 조합 가능한 경우의 수에 따라 비교(pairwise comparison)하되 방법은 Games-Howell 검정을 사용하라.

dunnettT3Test(QOL~RESID, data=df) : QOL의 평균을  RESID의 조합 가능한 경우의 수에 따라 비교(pairwise comparison)하되 방법은 Dunnett's T3 검정을 사용하라.

 

결과

	Pairwise comparisons using Tamhane's T2-test for unequal variances

data: QOL by RESID

  0      1     
1 <2e-16 -     
2 <2e-16 <2e-16

P value adjustment method: T2 (Sidak)
alternative hypothesis: two.sided


	Pairwise comparisons using Games-Howell test

data: QOL by RESID

  0       1      
1 1.7e-10 -      
2 < 2e-16 4.9e-10

P value adjustment method: none
alternative hypothesis: two.sided


	Pairwise comparisons using Dunnett's T3 test for multiple comparisons
		with unequal variances

data: QOL by RESID

  0      1     
1 <2e-16 -     
2 <2e-16 <2e-16

P value adjustment method: single-step
alternative hypothesis: two.sided

Pairwise comparisons using Tamhane's T2-test for unequal variances

data: QOL by RESID

  0      1     
<2e-16 -     
<2e-16 <2e-16

P value adjustment method: T2 (Sidak)
alternative hypothesis: two.sided


Pairwise comparisons using Games-Howell test

data: QOL by RESID

  0       1      
1.7e-10 -      
< 2e-16 4.9e-10

P value adjustment method: none
alternative hypothesis: two.sided


Pairwise comparisons using Dunnett's T3 test for multiple comparisons
with unequal variances

data: QOL by RESID

  0      1     
<2e-16 -     
<2e-16 <2e-16

P value adjustment method: single-step
alternative hypothesis: two.sided

 

빨간색 글씨가 p-value를 의미하는데, 이는 다음과 같이 보면 된다.

  0 1
1 대도시 거주자 (0)과 중소도시 거주자 (1) 비교  
2 대도시 거주자 (0)과 시골 거주자 (2) 비교 중소도시 거주자 (1)과 시골 거주자 (2) 비교

모든 p-value가 0.05보다 작으므로 다음과 같이 결론 내릴 수 있다.

대도시 거주자 (0)과 중소도시 거주자 (1) 사이에는 QOL 평균에 차이가 있다.

대도시 거주자 (0)과 시골 거주자 (2) 사이에는 QOL 평균에 차이가 있다.

중소도시 거주자 (1)과 시골 거주자 (2) 사이에는 QOL 평균에 차이가 있다.

 

어떤 사후 분석을 쓸 것인가

 이 논의에 대해 정답이 따로 있는 것은 아니다. 적절한 방법을 사용하여 논문에 제시하면 되고, 어떤 것이 정답이라고 콕 집어 이야기할 수는 없다. 다만, 사후 분석 방법이 여러 가지가 있다는 것은 '사후 분석 방법에 따라 산출되는 결과가 달라질 수 있다.'는 것을 의미하고, 심지어는 '어떤 사후 분석 방법을 채택하냐에 따라 유의성 여부가 달라질 수도 있다.'는 것을 의미한다. 심지어, Welch's or Brown-Forsythe's ANOVA에서는 유의한 결과가 나왔는데, 사후 분석을 해보니 유의한 차이를 보이는 경우가 없을 수도 있다. 따라서 어떤 사후 분석 방법 결과에 따른 결과인지 유의하여 해석할 필요가 있다.

 

코드 정리

##워킹 디렉토리 지정
setwd("C:/Users/user/Documents/Tistory_blog")


##데이터 불러오기
install.packages("readr")
library("readr")
df<-read_csv("Data.csv")


##데이터 나누기
df_big<-df[df$RESID==0,]
df_med<-df[df$RESID==1,]
df_sma<-df[df$RESID==2,]


##정규성 검정하기
#대도시 지역 거주자
# 1) Q-Q plot 그리기
qqnorm(df_big$QOL)
qqline(df_big$QOL)
# 2) 히스토그램 그리기
hist(df_big$QOL, prob=TRUE)
QOLrange<-seq(min(df_big$QOL),max(df_big$QOL),length=max(max(df_big$QOL)-min(df_big$QOL),100))
ND<-dnorm(QOLrange,mean=mean(df_big$QOL),sd=sd(df_big$QOL))
lines(QOLrange, ND, lwd=2)
# 3) Shapiro-Wilk test 시행하기
shapiro.test(df_big$QOL)

#중소도시 지역 거주자
# 1) Q-Q plot 그리기
qqnorm(df_med$QOL)
qqline(df_med$QOL)
# 2) 히스토그램 그리기
hist(df_med$QOL, prob=TRUE)
QOLrange<-seq(min(df_med$QOL),max(df_med$QOL),length=max(max(df_med$QOL)-min(df_med$QOL),100))
ND<-dnorm(QOLrange,mean=mean(df_med$QOL),sd=sd(df_med$QOL))
lines(QOLrange, ND, lwd=2)
# 3) Shapiro-Wilk test 시행하기
shapiro.test(df_med$QOL)

#시골지역 거주자
# 1) Q-Q plot 그리기
qqnorm(df_sma$QOL)
qqline(df_sma$QOL)
# 2) 히스토그램 그리기
hist(df_sma$QOL, prob=TRUE)
QOLrange<-seq(min(df_sma$QOL),max(df_sma$QOL),length=max(max(df_sma$QOL)-min(df_sma$QOL),100))
ND<-dnorm(QOLrange,mean=mean(df_sma$QOL),sd=sd(df_sma$QOL))
lines(QOLrange, ND, lwd=2)
# 3) Shapiro-Wilk test 시행하기
shapiro.test(df_sma$QOL)


##등분산성 검정
#Levene의 등분산 검정
install.packages("lawstat")
library("lawstat")
levene.test(df$QOL, df$RESID, location="mean")

#Brown&Forsythe
install.packages("lawstat")
library("lawstat")
levene.test(df$QOL, df$RESID)

#Bartlett
bartlett.test(df$QOL, df$RESID)


##ANOVA
#Welch's ANOVA
oneway.test(QOL~RESID, data=df, var.equal=FALSE)

#Brown and Forsythe's ANOVA
install.packages("onewaytests")
library("onewaytests")
df$RESID<-as.factor(df$RESID)
bf.test(QOL~RESID, data=df)


##사후분석
#1) Tamhane's T2
install.packages("PMCMRplus")
library("PMCMRplus")
tamhaneT2Test(QOL~RESID, data=df)

#2) Games-Howell
gamesHowellTest(QOL~RESID, data=df)

#3) Dunnett's T3
dunnettT3Test(QOL~RESID, data=df)

 

[R] 등분산 가정이 성립하지 않는 일원 배치 분산 분석 (Welch's ANOVA, Brown-Forsythe ANOVA) 정복 완료!

 

작성일: 2022.11.30.

최종 수정일: 2022.12.23.

이용 프로그램: R 4.2.2

RStudio v2022.07.2

RStudio 2022.07.2+576 "Spotted Wakerobin" Release

운영체제: Windows 10, Mac OS 12.6.1

반응형
반응형

[R] 일표본 윌콕슨 부호 순위 검정 (비모수 일표본 중앙값 검정: One-Sample Wilcoxon Signed Rank Test) - wilcox.test()

 어떤 분포의 평균이 특정 값인지 확인하는 방법을 이전에는 일표본 T검정 (One-Sample T test)로 시행했었다. (2022.11.03 - [모평균 검정/R] - [R] 일표본 T검정 (One-sample T-test) - t.test()) 하지만 여기에는 중요한 가정이 필요한데, 분포가 정규성을 따르는 것이다. 하지만 분포가 정규성을 따르지 않는다면 어떻게 해야 할까? 그럴 때 사용하는 것이 One-Sample Wilcoxon Signed Rank Test (일표본 윌콕슨 부호 순위 검정)이다. 

 

 이번 포스팅에서는 One-Sample Wilcoxon Signed Rank Test (일표본 윌콕슨 부호 순위 검정)에 대해서 알아볼 것이다.

 

*실습용 데이터는 아래 링크를 클릭하면 다운로드할 수 있습니다.

2022.08.04 - [공지사항 및 소개] - 분석용 데이터 (update 22.11.28)

 

분석용 데이터 (update 22.11.28)

2022년 11월 28일 버전입니다. 변수는 계속하여 추가될 예정입니다. 다음 카테고리에 있는 글에서 이용된 데이터입니다. - 기술 통계 - 범주형 자료 분석 - 모평균 검정 - 반복 측정 자료 분석 - 통계

medistat.tistory.com

 

코드를 보여드리기에 앞서 워킹 디렉토리부터 지정하겠다.

워킹 디렉토리에 관한 설명은 다음 링크된 포스트에서 볼 수 있다.

2022.08.05 - [통계 프로그램 사용 방법/R] - [R] 작업 디렉토리 (Working Directory) 지정 - getwd(), setwd()

 

setwd("C:/Users/user/Documents/Tistory_blog")

 

데이터를 불러와 df에 객체로 저장하겠다.

데이터 불러오는 방법은 다음 링크에서 볼 수 있다.

2022.08.05 - [통계 프로그램 사용 방법/R] - [R] 데이터 불러오기 : EXCEL - read_excel(), read.xlsx()

2022.08.10 - [통계 프로그램 사용 방법/R] - [R] 데이터 저장하기 : CSV 파일 - write.csv(), write_csv()

2022.08.10 - [통계 프로그램 사용 방법/R] - [R] 데이터 불러오기 : SAS file (.sas7bdat) - read.sas7bdat(), read_sas()

 

install.packages("readr")
library("readr")
df<-read_csv("Data.csv")

 

목표: 모집단에서 신경심리검사 1 원점수의 중앙값이 55이라고 할 수 있는가?

 이번 포스팅의 목적은 1000명의 데이터를 가지고, 이 1000명이 기원한 모집단에서 신경심리검사 1 원점수의 중앙값이 55이라고 할 수 있는지 판단하는 것이다.

 

전제조건 (정규성) - 코드

 만약 정규성을 따른다면 t-test를 하면 되므로 정규성 여부를 파악하도록 한다. 정규성 검정에 관한 분석 내용은 다음 글에서 살펴볼 수 있다. 

2022.08.11 - [기술 통계/R] - [R] 정규성 검정 (1) : Q-Q plot - qqnorm(), qqline()

2022.08.11 - [기술 통계/R] - [R] 정규성 검정 (2) : 히스토그램 - hist(), dnorm()

2022.08.11 - [기술 통계/R] - [R] 정규성 검정 (3) : 정량적 검정 (Shapiro-Wilk, Kolmogorov-Smirnov) - shapiro.test(), ks.test()

2022.08.12 - [기술 통계/R] - [R] 정규성 검정 (4) : 정량적 검정 (Lilliefors test) - lillie.test()

2022.08.16 - [기술 통계/R] - [R] 고급 Q-Q Plot - Van der Waerden, Rankit, Tukey, Blom

 

# 1) Q-Q plot 그리기
qqnorm(df$NP1_raw)
qqline(df$NP1_raw)
# 2) 히스토그램 그리기
hist(df$NP1_raw, prob=TRUE)
NP1_rawrange<-seq(min(df$NP1_raw),max(df$NP1_raw),length=max(max(df$NP1_raw)-min(df$NP1_raw),100))
ND<-dnorm(NP1_rawrange,mean=mean(df$NP1_raw),sd=sd(df$NP1_raw))
lines(NP1_rawrange, ND, lwd=2)
# 3) Shapiro-Wilk test 시행하기
shapiro.test(df$NP1_raw)

결과

	Shapiro-Wilk normality test

data:  df$NP1_raw
W = 0.97015, p-value = 1.81e-13

 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이다.

 

하지만 대립 가설은 3개가 된다.

1) $H_1=$ 모집단의 중앙값은 55보다 크다. (단측 검정)

2) $H_1=$ 모집단의 중앙값은 55보다 작다. (단측 검정)

3) $H_1=$ 모집단의 중앙값은 55보다 크거나 작다. (양측 검정)

 

대립 가설 별로 코드를 보고자 한다.

 

1)  $H_1=$ 모집단의 중앙값은 55보다 크다. (단측 검정)

wilcox.test(df$NP1_raw, mu=55, alternative="greater", correct=FALSE)

wilcox.test(df$NP1_raw, mu=55, alternative="greater", correct=FALSE) : Wilcoxon test를 시행하라. 검정할 변수는 df 데이터의 "NP1_raw"이고, 중앙값이 55와 같은지를 검정해라. 대립 가설은 모집단의 중앙값이 55보다 크다(greater)이고, 연속성 수정은 시행하지 마라. (연속성 수정은 다음 글을 참고하길 바란다.2022.08.30 - [통계 이론] - [이론] 연속성을 수정한 카이 제곱 검정 (Chi-squared test with Yates's correction for continuity))

 

결과

	Wilcoxon signed rank test

data:  df$NP1_raw
V = 156889, p-value = 1
alternative hypothesis: true location is greater than 55

Wilcoxon signed rank test

data:  df$NP1_raw

V = 156889, p-value = 1

alternative hypothesis: true location is greater than 55

 

P-value는 1.0000이다. 즉, 대립 가설 (모집단의 중앙값은 55보다 크다.)을 선택하는 것이 말이 안 된다는 것이다. 중앙값이 55보다 크다고 보는 것보다는 55이라고 보는 것이 더 합리적이라는 뜻이다. 

 

2)  $H_1=$ 모집단의 평균은 55보다 작다. (단측 검정)

wilcox.test(df$NP1_raw, mu=55, alternative="less", correct=FALSE)

wilcox.test(df$NP1_raw, mu=55, alternative="less", correct=FALSE) : Wilcoxon test를 시행하라. 검정할 변수는 df 데이터의 "NP1_raw"이고, 중앙값이 55와 같은지를 검정해라. 대립 가설은 '모집단의 중앙값이 55보다 작다(less)'이고, 연속성 수정은 시행하지 마라.

 

결과

	Wilcoxon signed rank test

data:  df$NP1_raw
V = 156889, p-value < 2.2e-16
alternative hypothesis: true location is less than 55

Wilcoxon signed rank test

data:  df$NP1_raw

V = 156889, p-value < 2.2e-16

alternative hypothesis: true location is less than 55

 

P-value2.2e-16보다 작다. 즉, 대립 가설 (모집단의 중앙값은 55보다 작다.)을 선택하는 것이 합당하다는 뜻이다. 중앙값이 55이라고 보는 것 보다는 55보다 작다고 보는 것이 더 합리적이라는 뜻이다. 

 

3)  $H_1=$ 모집단의 평균은 55보다 작거나 크다. (양측 검정)

wilcox.test(df$NP1_raw, mu=55, alternative="two.sided", correct=FALSE)

wilcox.test(df$NP1_raw, mu=55, correct=FALSE)

위 두 개의 코드는 같은 결과를 내준다. 왜냐하면 양측 검정 (alternative="two.sided")가 기본 값이어서 alternative 옵션을 지정해주지 않으면 alternative="two.sided"라고 인식하고 코드가 돌아가기 때문이다.

 

wilcox.test(df$NP1_raw, mu=55, alternative="two.sided", correct=FALSE) : Wilcoxon test를 시행하라. 검정할 변수는 df 데이터의 "NP1_raw"이고, 중앙값이 55와 같은지를 검정해라. 대립 가설은 '모집단의 중앙값은 55가 아니다(two.sided)'이고, 연속성 수정은 시행하지 마라.

 

결과

	Wilcoxon signed rank test

data:  df$NP1_raw
V = 156889, p-value < 2.2e-16
alternative hypothesis: true location is not equal to 55

Wilcoxon signed rank test

data:  df$NP1_raw

V = 156889, p-value < 2.2e-16

alternative hypothesis: true location is not equal to 55

 

P-value는 2.2e-16보다 작다. 즉, 대립 가설 (모집단의 중앙값은 55이 아니다.)을 선택하는 것이 합당하다는 뜻이다. 중앙값이 55이라고 보는 것보다는 55보다 작거나 크다고 보는 것이 더 합리적이라는 뜻이다. 

 

보통은 양측검정 결과를 사용하게 된다.

 

그러면 실제 중앙값은 얼마이길래 이런 결론이 나왔던 것일까?

median(df$NP1_raw)
[1] 49.96328

49.96328로 거의 50에 근사한다. 그렇다면 만약 50에 대해 검정하면 어떤 결과를 내보일까?

 

코드

wilcox.test(df$NP1_raw, mu=50, alternative="two.sided", correct=FALSE)

결과

	Wilcoxon signed rank test

data:  df$NP1_raw
V = 246178, p-value = 0.6558
alternative hypothesis: true location is not equal to 50

p-value가 0.6558>0.05로 귀무가설을 채택하는 것이 합리적이며, 중앙값이 50이 아니라고는 할 수 없는 상황이다.

 

 

중요한 가정: 대칭성

위 논의에서 빠진 정말 중요한 가정이 하나 있다. 신경심리검사 1 원점수 (NP1_raw)의 분포가 어떤 점수를 기점으로 좌우대칭이라는 것이다. 윌콕슨 부호 순위 검정은 좌우대칭일 때 그 검정력이 극대화되며, 좌우대칭이 아니면 정확도가 떨어진다. 신경심리검사 1 원점수 (NP1_raw)의 분포를 히스토그램을 통해 알아보자 (히스토그램 그리는 방법은 다음 링크를 참고하길 바란다.2022.08.11 - [기술 통계/R] - [R] 정규성 검정 (2) : 히스토그램 - hist(), dnorm())

#좀 더 자세히 보기 위해 breaks=40이라는 옵션을 넣었다.
hist(df$NP1_raw, breaks=40)

봉우리가 세 개인 Trimodal shape을 하고 있고, 중앙값인 50을 기준으로 대칭인 것처럼 보인다. 따라서 위 분석에는 문제가 없는 듯하다.

 

 

가정이 성립하지 않는다면?

우리가 궁금하고, 문제가 되는 상황은 가정이 성립하지 않는 상황이다. 실제로 이럴 때가 많기 때문이다. 신경심리검사 2 원점수인 NP2_raw데이터의 분포를 확인해보자.

hist(df$NP2_raw, breaks=80)

이는 심하게 오른쪽으로 치우쳐진 right skewed data다. 이때 중앙값에 대해 Wilcoxon signed rank test를 시행하면 어떻게 될까?

wilcox.test(df$NP2_raw, mu=median(df$NP2_raw), alternative="two.sided", correct=FALSE)
	Wilcoxon signed rank test

data:  df$NP2_raw
V = 295522, p-value = 7.213e-07
alternative hypothesis: true location is not equal to 1.004532

참 중앙값으로 "NP2_raw의 중앙값이 참 중앙값과 같니?"라고 물어보았는데 "아니요"라고 대답하는 참사가 벌어진다. 이는 분포에 대칭성이 없기 때문이다. 이럴 때 선택할 수 있는 옵션이 많지 않긴 한데, 그중 하나는 변환 (transformation)이다. 이런 right skewed data는 로그 변환 (log-transformation)을 하면 대칭성이 갖추어지는 경우가 많다. 로그 변환을 하고 히스토그램을 그려 분포를 확인해보자. (로그 변환에 관한 글은 다음 링크에서 확인할 수 있다.2022.11.25 - [통계 프로그램 사용 방법/R] - [R] 변수 계산 (산술 연산))

df$log_NP2_raw<-log(df$NP2_raw)
hist(df$log_NP2_raw, breaks=80)

충분히 대칭성을 띄고 있다고 할 수 있다.  "혹시 정규성을 띄지 않을까?"라는 생각으로 정규성을 확인해보면 다음과 같음을 알 수 있다. 정규성 검정 방법은 아래 링크들을 확인하길 바란다.

2022.08.11 - [기술 통계/R] - [R] 정규성 검정 (1) : Q-Q plot - qqnorm(), qqline()

2022.08.11 - [기술 통계/R] - [R] 정규성 검정 (2) : 히스토그램 - hist(), dnorm()

2022.08.11 - [기술 통계/R] - [R] 정규성 검정 (3) : 정량적 검정 (Shapiro-Wilk, Kolmogorov-Smirnov) - shapiro.test(), ks.test()

2022.08.12 - [기술 통계/R] - [R] 정규성 검정 (4) : 정량적 검정 (Lilliefors test) - lillie.test()

2022.08.16 - [기술 통계/R] - [R] 고급 Q-Q Plot - Van der Waerden, Rankit, Tukey, Blom

 

코드

# 1) Q-Q plot 그리기
qqnorm(df$log_NP2_raw)
qqline(df$log_NP2_raw)
# 2) 히스토그램 그리기
hist(df$log_NP2_raw, prob=TRUE)
log_NP2_rawrange<-seq(min(df$log_NP2_raw),max(df$log_NP2_raw),length=max(max(df$log_NP2_raw)-min(df$log_NP2_raw),100))
ND<-dnorm(log_NP2_rawrange,mean=mean(df$log_NP2_raw),sd=sd(df$log_NP2_raw))
lines(log_NP2_rawrange, ND, lwd=2)
# 3) Shapiro-Wilk test 시행하기
shapiro.test(df$log_NP2_raw)

결과

	Shapiro-Wilk normality test

data:  df$log_NP2_raw
W = 0.99284, p-value = 9.487e-05

히스토그램에서는 정규분포를 따르는 것 같아 보이지만, QQ plot에서 직선 위에 있지 않은 점들이 많고, Shapiro-wilk test에서 p-value가 매우 작게 나오므로 정규성을 띠지 않는다고 판단하는 것이 맞다. 따라서 일표본 t검정(One-sample t-test)이 아닌 일표본 윌콕슨 부호 순위 검정 (One sample Wilcoxon signed rank test)을 시행해야 한다. 위와 같이 log_NP2_raw의 중앙값에 대해 검정해보자.

 

코드

wilcox.test(df$log_NP2_raw, mu=median(df$log_NP2_raw), alternative="two.sided", correct=FALSE)

 

결과

	Wilcoxon signed rank test

data:  df$log_NP2_raw
V = 247173, p-value = 0.7363
alternative hypothesis: true location is not equal to 0.004521896

p-value가 0.7363으로 귀무가설을 기각할 수 없게 되어 중앙값은 0.004521896이 아니라고 할 수 없다고 결론을 내리게 된다.

 

 

그 외 다룰 내용 1: 정규성을 안 따르면 꼭 Wilcoxon test를 시행해야 하나?

 T-test의 전제조건은 표본 평균이 정규분포를 따르는 것이다. 그런데 중심 극한 정리에 따라 표본의 수 (우리 데이터에서는 n=1,000)가 커질수록 표본 평균의 분포는 정규성을 띠게 된다. 따라서 n수가 적당히 크기만 하면 표본이 정규분포를 따르는지와 관계없이 표본평균은 정규성을 따른다고 할 수 있다. 이런 이유로 정규성과 관계 없이 n수가 크기만 하면 t-test의 이용이 정당화되기도 한다. 합리적인 말이다.

 그런데 왜 정규분포를 따르는지 왜 다들 확인할까? 만약 모분포가 정규분포였다면 표본도 정규분포를 따를 확률이 높아진다. 물론 꼭 그런 건 아니지만 말이다. 따라서 표본이 정규분포를 따른다면, 모분포를 정규분포를 따른다고 할 수 있고, 그렇다면 n수가 작아도 표본 평균의 분포도 정규성을 따르니 t-test를 사용해도 되는 안전한 환경이 구축된다. 이런 이유로 정규성을 확인하게 된다.

 

그 외 다룰 내용 2: Wilcoxon signed rank test는 짝지어진 자료의 검정방법 아닌가?

 맞다. 하지만 one sample test에도 사용될 수 있다. 모수적 방법에서 대응 표본 T 검정 (paired T test)가 사실 일표본 T 검정 (One sample T test)와 같다는 것과 일맥상통하는 이야기다. (대응 표본 T검정은 다음 링크에서 확인할 수 있다.2022.11.25 - [반복 측정 자료 분석/R] - [R] 대응 표본 T검정 (Paired samples T-test) - t.test())

 

코드 정리

##워킹 디렉토리 지정
setwd("C:/Users/user/Documents/Tistory_blog")

##데이터 불러오기
install.packages("readr")
library("readr")
df<-read_csv("Data.csv")

###########NP1_raw###########

##정규성 검정
# 1) Q-Q plot 그리기
qqnorm(df$NP1_raw)
qqline(df$NP1_raw)
# 2) 히스토그램 그리기
hist(df$NP1_raw, prob=TRUE)
NP1_rawrange<-seq(min(df$NP1_raw),max(df$NP1_raw),length=max(max(df$NP1_raw)-min(df$NP1_raw),100))
ND<-dnorm(NP1_rawrange,mean=mean(df$NP1_raw),sd=sd(df$NP1_raw))
lines(NP1_rawrange, ND, lwd=2)
# 3) Shapiro-Wilk test 시행하기
shapiro.test(df$NP1_raw)

##Wilcoxon signed rank test
wilcox.test(df$NP1_raw, mu=55, alternative="greater", correct=FALSE)
wilcox.test(df$NP1_raw, mu=55, alternative="less", correct=FALSE)
wilcox.test(df$NP1_raw, mu=55, alternative="two.sided", correct=FALSE)
wilcox.test(df$NP1_raw, mu=55, correct=FALSE)

##중앙값 확인하기
median(df$NP1_raw)

##중앙값 근사치(50)으로 검정하기
wilcox.test(df$NP1_raw, mu=50, alternative="two.sided", correct=FALSE)

##NP1_raw 대칭성 확인
hist(df$NP1_raw)
hist(df$NP1_raw, breaks=40)

###########NP2_raw###########

##NP2_raw 대칭성 확인
hist(df$NP2_raw, breaks=80)

##NP2_raw 비대칭이지만 wilcoxon test 시행하기
wilcox.test(df$NP2_raw, mu=median(df$NP2_raw), alternative="two.sided", correct=FALSE)

##로그변환 후 히스토그램으로 대칭성 확인하기
df$log_NP2_raw<-log(df$NP2_raw)
hist(df$log_NP2_raw, breaks=80)

##로그변환된 NP2_raw의 정규성 확인하기
# 1) Q-Q plot 그리기
qqnorm(df$log_NP2_raw)
qqline(df$log_NP2_raw)
# 2) 히스토그램 그리기
hist(df$log_NP2_raw, prob=TRUE)
log_NP2_rawrange<-seq(min(df$log_NP2_raw),max(df$log_NP2_raw),length=max(max(df$log_NP2_raw)-min(df$log_NP2_raw),100))
ND<-dnorm(log_NP2_rawrange,mean=mean(df$log_NP2_raw),sd=sd(df$log_NP2_raw))
lines(log_NP2_rawrange, ND, lwd=2)
# 3) Shapiro-Wilk test 시행하기
shapiro.test(df$log_NP2_raw)

##로그변환된 변수로 Wilcoxon test 시행하기.
wilcox.test(df$log_NP2_raw, mu=median(df$log_NP2_raw), alternative="two.sided", correct=FALSE)

 

[R] 일표본 윌콕슨 부호 순위 검정 (비모수 일표본 중앙값 검정: One-Sample Wilcoxon Signed Rank Test) 정복 완료!

 

작성일: 2022.11.29.

최종 수정일: 2022.11.29.

이용 프로그램: R 4.2.2

RStudio v2022.07.2

RStudio 2022.07.2+576 "Spotted Wakerobin" Release

운영체제: Windows 10, Mac OS 12.6.1

반응형
반응형

[R] 일원 배치 분산 분석 (One-way ANOVA, ANalysis Of VAriance) - aov(), pairwise.t.test(), ScheffeTest(), scheffe.test(), duncan.test(), TukeyHSD()

 

 이전 포스팅에서 음주 여부에 따른 ALT 평균에 차이가 있는지를 알아보기 위해서는 T test를 사용한다는 내용을 알아보았다.

2022.11.12 - [모평균 검정/R] - [R] 독립 표본 T검정 (Independent samples T-test) - t.test(), var.test(), levene.test()

 음주 여부는 "1) 음주자", "2) 비음주자"로 나뉘는 이분형 변수다. 즉, 음주자의 ALT 평균과 비음주자의 ALT 평균을 비교한 것이었다. 

 

 만약, 두 개 이상의 범주로 나뉘는 변수에 대해 모평균의 검정을 하고자 한다면 어떻게 해야 할까? 예를 들어, 흡연 상태를 "1) 비흡연자", "2) 과거 흡연자", "3) 현재 흡연자"로 나누었고, 각 그룹의 폐기능 검사 중 하나인 FVC (Forced Vital Capacity)의 평균에 차이가 있는지 알아보고자 한다. 

 

 이때 사용할 수 있는 방법이 일원 배치 분산 분석인 One-way ANOVA (Analysis of Variance)이며 이번 포스팅에서 다뤄볼 주된 내용이다.

 

*실습용 데이터는 아래 링크를 클릭하면 다운로드할 수 있습니다.

2022.08.04 - [공지사항 및 소개] - 분석용 데이터 (update 22.11.21)

 

분석용 데이터 (update 22.11.21)

2022년 11월 21일 버전입니다. 변수는 계속하여 추가될 예정입니다. 다음 카테고리에 있는 글에서 이용된 데이터입니다. - 기술 통계 - 범주형 자료 분석 - 모평균 검정 - 반복 측정 자료 분석 - 통계

medistat.tistory.com

 

코드를 보여드리기에 앞서 워킹 디렉토리부터 지정하겠다.

워킹 디렉토리에 관한 설명은 다음 링크된 포스트에서 볼 수 있다.

2022.08.05 - [통계 프로그램 사용 방법/R] - [R] 작업 디렉토리 (Working Directory) 지정 - getwd(), setwd()

setwd("C:/Users/user/Documents/Tistory_blog")

 

데이터를 불러와 df에 객체로 저장하겠다.

데이터 불러오는 방법은 다음 링크에서 볼 수 있다.

2022.08.05 - [통계 프로그램 사용 방법/R] - [R] 데이터 불러오기 : EXCEL - read_excel(), read.xlsx()

2022.08.10 - [통계 프로그램 사용 방법/R] - [R] 데이터 저장하기 : CSV 파일 - write.csv(), write_csv()

2022.08.10 - [통계 프로그램 사용 방법/R] - [R] 데이터 불러오기 : SAS file (.sas7bdat) - read.sas7bdat(), read_sas()

install.packages("readr")
library("readr")
df<-read_csv("Data.csv")

 

목표: 흡연 상태에 따라 FVC 평균은 서로 다르다고 모집단 수준에서 말할 수 있는가?

 

전제조건 (정규성, 등분산성)

 ANOVA에는 두 가지 전제조건이 필요하다. 본 포스팅의 예시에 맞추어 설명하면 다음과 같다.

 1) 정규성: 현재 흡연자, 과거 흡연자, 비흡연자 별로 FVC의 분포를 보았을 때 각각의 분포는 정규성을 따른다. 

 2) 등분산성: 현재 흡연자, 과거 흡연자, 비흡연자 별로 FVC의 분포를 보았을 때 분산은 서로 같다.

 

따라서 각각을 검정해야 한다.

 

전제조건 (정규성) - 코드

 

일원 배치 분산 분석 (One-way ANOVA)의 첫 번째 전제조건은 각각의 독립 표본의 분포가 정규성을 따른다는 것이다. 즉, 여기에서는 흡연 상태 (비흡연자, 과거 흡연자, 현재 흡연자)에 따라 수축기 혈압 (SBP)의 정규성을 검정하도록 하겠다.

 

따라서 다음 두 가지의 일을 해야 한다.

1) 흡연 상태 (비흡연자, 과거 흡연자, 현재 흡연자)에 따라 데이터를 나누기

나누는 방법에 대한 설명은 다음 링크에서 볼 수 있다. 2022.11.10 - [통계 프로그램 사용 방법/R] - [R] 조건에 맞는 자료 추출하기

위 링크에서 확인할 수 있듯이 여러 가지 방법으로 나눌 수 있지만 indexing을 이용하여 나누도록 하겠다.

 

df_ns<-df[df$SMOK==0,]
df_fs<-df[df$SMOK==1,]
df_cs<-df[df$SMOK==2,]

df_ns : 비흡연자로 이루어진 데이터

df_fs : 과거 흡연자로 이루어진 데이터

df_cs : 현재 흡연자로 이루어진 데이터

 

2) 흡연 상태 (비흡연자, 과거 흡연자, 현재 흡연자)에 따라 정규성을 검정하기

정규성 검정에 관한 내용은 다음 링크에서 확인할 수 있다.

2022.08.11 - [기술 통계/R] - [R] 정규성 검정 (1) : Q-Q plot - qqnorm(), qqline()

2022.08.11 - [기술 통계/R] - [R] 정규성 검정 (2) : 히스토그램 - hist(), dnorm()

2022.08.11 - [기술 통계/R] - [R] 정규성 검정 (3) : 정량적 검정 (Shapiro-Wilk, Kolmogorov-Smirnov) - shapiro.test(), ks.test()

2022.08.12 - [기술 통계/R] - [R] 정규성 검정 (4) : 정량적 검정 (Lilliefors test) - lillie.test()

2022.08.16 - [기술 통계/R] - [R] 고급 Q-Q Plot - Van der Waerden, Rankit, Tukey, Blom

##비흡연자의 FVC 정규성 검정
# 1) Q-Q plot 그리기
qqnorm(df_ns$FVC_pPRED)
qqline(df_ns$FVC_pPRED)
# 2) 히스토그램 그리기
hist(df_ns$FVC_pPRED, prob=TRUE)
FVC_pPREDrange<-seq(min(df_ns$FVC_pPRED),max(df_ns$FVC_pPRED),length=max(max(df_ns$FVC_pPRED)-min(df_ns$FVC_pPRED),100))
ND<-dnorm(FVC_pPREDrange,mean=mean(df_ns$FVC_pPRED),sd=sd(df_ns$FVC_pPRED))
lines(FVC_pPREDrange, ND, lwd=2)
# 3) Shapiro-Wilk test 시행하기
shapiro.test(df_ns$FVC_pPRED)

##과거 흡연자의 FVC 정규성 검정
# 1) Q-Q plot 그리기
qqnorm(df_fs$FVC_pPRED)
qqline(df_fs$FVC_pPRED)
# 2) 히스토그램 그리기
hist(df_fs$FVC_pPRED, prob=TRUE)
FVC_pPREDrange<-seq(min(df_fs$FVC_pPRED),max(df_fs$FVC_pPRED),length=max(max(df_fs$FVC_pPRED)-min(df_fs$FVC_pPRED),100))
ND<-dnorm(FVC_pPREDrange,mean=mean(df_fs$FVC_pPRED),sd=sd(df_fs$FVC_pPRED))
lines(FVC_pPREDrange, ND, lwd=2)
# 3) Shapiro-Wilk test 시행하기
shapiro.test(df_fs$FVC_pPRED)

##과거 흡연자의 FVC 정규성 검정
# 1) Q-Q plot 그리기
qqnorm(df_cs$FVC_pPRED)
qqline(df_cs$FVC_pPRED)
# 2) 히스토그램 그리기
hist(df_cs$FVC_pPRED, prob=TRUE)
FVC_pPREDrange<-seq(min(df_cs$FVC_pPRED),max(df_cs$FVC_pPRED),length=max(max(df_cs$FVC_pPRED)-min(df_cs$FVC_pPRED),100))
ND<-dnorm(FVC_pPREDrange,mean=mean(df_cs$FVC_pPRED),sd=sd(df_cs$FVC_pPRED))
lines(FVC_pPREDrange, ND, lwd=2)
# 3) Shapiro-Wilk test 시행하기
shapiro.test(df_cs$FVC_pPRED)

 

정규성 검정 - 결과

1) 비흡연자

 

 

	Shapiro-Wilk normality test

data:  df_ns$FVC_pPRED
W = 0.99679, p-value = 0.4928

 N수가 2,000개 미만이므로 Shapiro-Wilk 통계량의 p-value를 보면 0.05 이상이며, Q-Q Plot은 대부분의 데이터가 선상에 있고, 히스토그램에서도 정규성을 따르는 것처럼 볼 수 있다. 

 

2) 과거 흡연자

	Shapiro-Wilk normality test

data:  df_fs$FVC_pPRED
W = 0.99727, p-value = 0.9043

 N수가 2,000개 미만이므로 Shapiro-Wilk 통계량의 p-value를 보면 0.05 이상이며, Q-Q Plot은 대부분의 데이터가 선상에 있고, 히스토그램에서도 정규성을 따르는 것처럼 볼 수 있다. 

 

3) 현재 흡연자

	Shapiro-Wilk normality test

data:  df_cs$FVC_pPRED
W = 0.98909, p-value = 0.06429

 N수가 2,000개 미만이므로 Shapiro-Wilk 통계량의 p-value를 보면 0.05 이상이며, Q-Q Plot은 대부분의 데이터가 선상에 있고, 히스토그램에서도 정규성을 따르는 것처럼 볼 수 있다. 

 

 첫 번째 전제조건이 성립한다. 즉, 흡연 상태에 따른 FVC의 분포가 정규성을 따른다고 할 수 있다. 

 

2) 등분산성 

등분산성 검정에 관한 내용은 다음 링크에서 확인할 수 있다.

2022.11.12 - [모평균 검정/R] - [R] 독립 표본 T검정 (Independent samples T-test) - t.test(), var.test(), levene.test()

2022.11.20 - [모평균 검정/R] - [R] 등분산성 검정 (Homogeneity of variance) - levene.test(), bartlett.test()

 

많이들 사용하는 등분산 검정 다섯 가지 중 F test는 그룹이 2개일 때에만 사용 가능하고, O'brien은 R에서 시행 불가하므로, 나머지 3가지 방법으로 검정해본다.

1) F test (그룹이 2개일 때에만 사용 가능)

2) Levene의 등분산 검정

3) O'Brien의 등분산 검정(R에서 검정 불가)

4) Brown and Forsythe의 등분산 검정

5) Bartlett의 등분산 검정

 

코드

#Levene의 등분산 검정
install.packages("lawstat")
library("lawstat")
levene.test(df$FVC_pPRED, df$SMOK, location="mean")

#Brown&Forsythe
install.packages("lawstat")
library("lawstat")
levene.test(df$FVC_pPRED, df$SMOK)

#Bartlett
bartlett.test(df$FVC_pPRED, df$SMOK)

 

결과

	Classical Levene's test based on the absolute deviations
	from the mean ( none not applied because the location is
	not set to median )

data:  df$FVC_pPRED
Test Statistic = 1.3065, p-value = 0.2712

	Modified robust Brown-Forsythe Levene-type test based on
	the absolute deviations from the median

data:  df$FVC_pPRED
Test Statistic = 1.2772, p-value = 0.2793

	Bartlett test of homogeneity of variances

data:  df$FVC_pPRED and df$SMOK
Bartlett's K-squared = 1.8171, df = 2, p-value = 0.4031

 Levene 검정 (p-value: 0.2712), Brown-Forsythe test (p-value: 0.2793), Bartlett test (p-value: 0.4031)로 모두 0.05보다 큰 p-value를 갖는다. 따라서 귀무가설을 기각하지 못하므로, 흡연 상태에 따라  FVC_pPRED의 분산이 동일하다는 귀무가설을 채택한다.

 

따라서 ANOVA의 두 가지 전제조건 (정규성, 등분산성)이 성립하므로 ANOVA를 실행할 수 있다.

 

범주형 변수로의 변환

ANOVA를 시행하기 앞서, df데이터의 SMOK변수는 "비흡연자", "과거 흡연자", "현재 흡연자"를 각각 0,1,2로 코딩한 범주형 변수임을 R에게 알려주어야 한다. 알려주지 않으면 0, 1, 2로 이루어진 숫자형 변수로 인식하고 있기 때문이다. 따라서 변수의 유형을 범주형으로 바꾸어준다. (상세한 내용은 다음 링크를 확인하길 바란다.2022.11.21 - [통계 프로그램 사용 방법/R] - [R] 변수의 유형 (타입, type) 확인 및 변경 - as.factor(), as.numeric(), as.character(), str())

 

df$SMOK<-as.factor(df$SMOK)

 

일원 배치 분산 분석 (One-way ANOVA) 시행

코드

a<-aov(FVC_pPRED~SMOK, data=df)
summary(a)

a<-aov(FVC_pPRED~SMOK, data=df) : df 데이터의 SMOK에 따른 FVC_pPRED의 평균이 서로 다른지 ANOVA를 시행하고, 이를 a라는 객체에 저장하라.
summary(a) : a를 요약하여 보여달라.

 

결과

             Df Sum Sq Mean Sq F value Pr(>F)    
SMOK          2  56791   28396   168.9 <2e-16 ***
Residuals   997 167579     168                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                      Df  Sum Sq  Mean Sq  F value   Pr(>F)    

SMOK             2    56791      28396     168.9 <2e-16 ***

Residuals   997  167579           168                    

 

결과는 복잡하지만 우리가 볼 곳은 빨간 글씨로 칠해놓은 단 하나의 값, p-value다. 값은 $2\times 10^{-16}$보다 작으므로 이는 0.05보다 한참 작은 값이므로 귀무가설을 기각하고 대립 가설을 채택한다. 그렇다면 여기에서 귀무가설 및 대립 가설은 무엇이었는가?

 

 귀무가설:  $H0=$ 세 집단의 모평균은 "모두" 동일하다.

 대립 가설:  $H1=$ 세 집단의 모평균이 모두 동일한 것은 아니다.

 

우리는 대립 가설을 채택해야 하므로 "세 집단의 모평균이 모두 동일한 것은 아니다."라고 결론 내릴 것이다.

이 말은, 세 집단 중 두 개씩 골라 비교했을 때, 적어도 한 쌍에서는 차이가 난다는 것이다. 따라서 세 집단 중 두 개씩 골라 비교를 해보아야 하며, 이를 사후 분석 (post hoc analysis)이라고 한다. 세 집단에서 두 개씩 고르므로 가능한 경우의 수는 $_3 C_2=3$이다.

 (1) 비흡연자 vs 과거 흡연자

 (2) 비흡연자 vs 현재 흡연자

 (3) 과거 흡연자 vs 현재 흡연자

 

사후 분석을 시행하겠다.

SAS를 이용한 ANOVA에 관한 글(2022.10.12 - [모평균 검정/SAS] - [SAS] 일원 배치 분산 분석 (One-way ANOVA, ANalysis Of VAriance) - PROC GLM, PROC ANOVA)에서 많은 사후 분석 방법이 있음을 밝혔고, 많이 사용되는 방법들에 대해서도 언급하였다. 본 포스팅에서는 많이 사용되는 다음 네 가지 사후 분석 방법을 소개한다.

 

1) Bonferroni t test

2) Scheffe's multiple comparison

3) Duncan's multiple range test

4) Tukey's studentized range test

 

 

1) Bonferroni t test - 코드

pairwise.t.test(df$FVC_pPRED, df$SMOK, p.adj="bonferroni")

pairwise.t.test(df$FVC_pPRED, df$SMOK, p.adj="bonferroni") :df 데이터의 FVC_pPRED를 SMOK변수값에 따라 각각 t-test를 시행하는데, p-value는 Bonferroni방식으로 보정하라.

 

결과

	Pairwise comparisons using t tests with pooled SD 

data:  df$FVC_pPRED and df$SMOK 

  0       1      
1 < 2e-16 -      
2 < 2e-16 8.7e-15

P value adjustment method: bonferroni

Pairwise comparisons using t tests with pooled SD 
data:  df$
FVC_pPRED and df$SMOK 

  0 1
1 <2e-16 -
2 <2e-16 8.7e-15

P value adjustment method: bonferroni 

 

결과가 복잡해 보일지라도 우리는 표만 잘 보면 된다.

SMOK에는 3가지 값 (비흡연자:0, 과거 흡연자:1, 현재 흡연자:2)이 들어있고, 이 중 둘을 뽑아 t-test를 시행할 것이다. 따라서 3개 중 2개를 뽑는 경우의 수인 3가지 p-value가 존재하게 된다. 

비흡연자(0)와 과거 흡연자(1)를 비교한 p-value는 빨간색으로 표시된 "<2e-16"이다.

비흡연자(0)와 현재 흡연자(2)를 비교한 p-value는 초록색으로 표시된 "<2e-16"이다.

과거 흡연자(1)와 현재 흡연자(2)를 비교한 p-value는 파란색으로 표시된 "8.7e-15"이다.

즉, 모든 조합에서 평균의 차이가 있다고 결론 내릴 수 있다.

 

따라서, 우리가 채택한 대립 가설 "세 집단의 모평균이 모두 동일한 것은 아니다."을 구체적으로 이야기하면 "비흡연자와 과거 흡연자, 비흡연자와 현재 흡연자, 과거 흡연자와 현재 흡연자 사이의 모평균에 차이가 있다."라고 결론내릴 수 있다.

 

 

2) Scheffe's multiple comparison - 코드 (방법 1)

install.packages("DescTools")
library("DescTools")
ScheffeTest(a)

ScheffeTest(a) :a(ANOVA 시행한 것)의 사후 분석을 Scheffe test로 시행하라.

이 함수는 "DecsTools"라는 패키지를 필요로 하므로 설치 후 사용하길 바란다.

 

결과

  Posthoc multiple comparisons of means: Scheffe Test 
    95% family-wise confidence level

$SMOK
          diff    lwr.ci     upr.ci    pval    
1-0  -9.532840 -11.90048  -7.165204 < 2e-16 ***
2-0 -18.553484 -21.07451 -16.032453 < 2e-16 ***
2-1  -9.020644 -11.77709  -6.264193 2.9e-14 ***

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

 

 

  Posthoc multiple comparisons of means: Scheffe Test 
    95% family-wise confidence level

$SMOK
                   diff         lwr.ci           upr.ci      pval    
1-0  -9.532840  -11.90048   -7.165204 < 2e-16 ***
2-0 -18.553484 -21.07451 -16.032453 < 2e-16 ***
2-1  -9.020644  -11.77709   -6.264193 2.9e-14 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

 

결과가 복잡해 보일지라도 우리는 표만 잘 보면 된다.

SMOK에는 3가지 값 (비흡연자:0, 과거 흡연자:1, 현재 흡연자:2)이 들어있고, 이 중 둘을 뽑아 사후 검정을 시행할 것이다. 따라서 3개 중 2개를 뽑는 경우의 수인 3가지 p-value가 존재하게 된다. 

비흡연자(0)와 과거 흡연자(1)를 비교한 p-value는 빨간색으로 표시된 "<2e-16"이다.

비흡연자(0)와 현재 흡연자(2)를 비교한 p-value는 초록색으로 표시된 "<2e-16"이다.

과거 흡연자(1)와 현재 흡연자(2)를 비교한 p-value는 파란색으로 표시된 "2.9e-14"이다.

즉, 모든 조합에서 평균의 차이가 있다고 결론 내릴 수 있다.

 

 

2-1) Scheffe's multiple comparison - 코드 (방법 2)

install.packages("agricolae")
library(agricolae)
scheffe.test(a,"SMOK",console=TRUE)

"agriocolae"패키지에 있는 scheffe.test() 함수를 이용할 수도 있다.

scheffe.test(a,"SMOK",console=TRUE) : a(ANOVA 시행한 것)에서 "SMOK"변수에 따라 Scheffe test를 시행하고 결과를 보여달라(console=TRUE).

 

결과

Study: a ~ "SMOK"

Scheffe Test for FVC_pPRED 

Mean Square Error  : 168.0834 

SMOK,  means

  FVC_pPRED      std   r      Min      Max
0  95.33911 13.31491 463 60.70763 130.0000
1  85.80627 12.39528 295 50.49873 124.0448
2  76.78563 12.96366 242 50.00000 116.7888

Alpha: 0.05 ; DF Error: 997 
Critical Value of F: 3.004752 

Groups according to probability of means differences and alpha level( 0.05 )

Means with the same letter are not significantly different.

  FVC_pPRED groups
0  95.33911      a
1  85.80627      b
2  76.78563      c

너무나도 복잡해 보이지만, 마지막 결과를 예쁘게 제시해주기 때문에 이 함수를 소개한 것이다.

  groups
0 a
1 b
2 c

0은 group a, 1은 group b, 2는 group c에 속한다. 즉 서로 다른 모평균을 갖는다고 결론 내릴 수 있다.

 

만약 유의미한 차이가 없는 경우 어떻게 제시될까?

df$RESID<-as.factor(df$RESID)
b<-aov(FVC_pPRED~RESID, data=df)
scheffe.test(b,"RESID",console=TRUE)

SMOK가 아니라 RESID에 대해 ANOVA를 시행하고, Scheffe test를 시행한 결과다.

Study: b ~ "RESID"

Scheffe Test for FVC_pPRED 

Mean Square Error  : 219.1342 

RESID,  means

  FVC_pPRED      std   r      Min      Max
0  84.73360 14.84186 323 50.32614 129.7551
1  88.59510 15.17535 335 50.00000 129.7026
2  90.61014 14.39179 342 50.47718 130.0000

Alpha: 0.05 ; DF Error: 997 
Critical Value of F: 3.004752 

Groups according to probability of means differences and alpha level( 0.05 )

Means with the same letter are not significantly different.

  FVC_pPRED groups
2  90.61014      a
1  88.59510      a
0  84.73360      b
  groups
0 a
1 a
2 b

위 표를 보면 비흡연자 (0)과 과거 흡연자 (1)은 모두 group a에 속해있다. 즉 비흡연자와 과거 흡연자 사이에는 유의미한 차이가 발견되지 않은 것이다. 반면 현재 흡연자는 group b에 속해있으므로 비흡연자나 과거 흡연자와는 유의미한 차이가 있었다는 것이다.

 

3) Duncan's multiple range test - 코드

 사실 Tukey나 다음에 나오는 Duncan은 각 그룹별 표본의 수가 같을 때 사용하는 것이므로 이 데이터의 사후 분석으로는 적절하지 않다. 하지만 우선 그냥 시행해보겠다.

install.packages("agricolae")
library(agricolae)
duncan.test(a,"SMOK",console=TRUE)

"agriocolae"패키지에 있는 duncan.test() 함수를 이용하므로 설치부터 한다.

duncan.test(a,"SMOK",console=TRUE) : a(ANOVA 시행한 것)에서 "SMOK"변수에 따라 Duncan test를 시행하고 결과를 보여달라(console=TRUE).

결과

Study: a ~ "SMOK"

Duncan's new multiple range test
for FVC_pPRED 

Mean Square Error:  168.0834 

SMOK,  means

  FVC_pPRED      std   r      Min      Max
0  95.33911 13.31491 463 60.70763 130.0000
1  85.80627 12.39528 295 50.49873 124.0448
2  76.78563 12.96366 242 50.00000 116.7888

Groups according to probability of means differences and alpha level( 0.05 )

Means with the same letter are not significantly different.

  FVC_pPRED groups
0  95.33911      a
1  85.80627      b
2  76.78563      c

마지막 세 줄을 보면 0은 group a, 1은 group b, 2는 group c에 속한다. 즉 서로 다른 모평균을 갖는다고 결론 내릴 수 있다.

 

4) Tukey's studentized range test 

TukeyHSD(a)

TukeyHSD(a) : a (ANOVA 시행한 것)으로 사후 분석을 시행하되 Tukey test를 사용하라.

 

결과

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

Fit: aov(formula = FVC_pPRED ~ SMOK, data = df)

$SMOK
          diff       lwr        upr p adj
1-0  -9.532840 -11.79982  -7.265858     0
2-0 -18.553484 -20.96734 -16.139628     0
2-1  -9.020644 -11.65991  -6.381377     0

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

Fit: aov(formula = FVC_pPRED ~ SMOK, data = df)

$SMOK
                 diff            lwr             upr     p adj
1-0   -9.532840 -11.79982   -7.265858     0
2-0 -18.553484 -20.96734 -16.139628     0
2-1   -9.020644 -11.65991   -6.381377     0

 

결과가 복잡하지만 "p adj"만 보면 된다. 어떤 비교든 p value가 너무 작아 0으로 나타난다. 따라서, 서로 다른 모평균을 갖는다고 결론 내릴 수 있다.

 

이 시점에서 이런 의문이 들 수 있다.

"각각의 그룹별로 평균을 비교하면 되지, 굳이 왜 ANOVA라는 방법까지 사용하는 것인가?"

 아주 논리적인 의문점이다. 하지만, 반드시 ANOVA를 사용해야 한다. 그 이유는 다음과 같다. 본 사례는 흡연 상태에 따른 조합 가능한 경우의 수가 3인데, 각각 유의성의 기준을 α=0.05로 잡아보자. 이때 세 번의 비교에서 모두 귀무가설이 학문적인 진실인데(평균에 차이가 없다.), 세 번 모두 귀무가설을 선택할 확률은 (1−0.95)3≈0.86이다. (이해가 어려우면 p-value에 대한 설명 글을 읽고 오길 바란다. 2022.09.05 - [통계 이론] - [이론] p-value에 관한 고찰)

 그런데, ANOVA의 귀무가설은 "모든 집단의 평균이 같다."이다. 따라서 모든 집단의 평균이 같은 것이 학문적 진실일 때, 적어도 한 번이라도 대립 가설을 선택하게 될 확률은 1−0.86=0.14가 된다. 즉, 유의성의 기준이 올라가게 되어, 덜 보수적인 결정을 내리게 되고, 다시 말하면 대립 가설을 잘 선택하는 쪽으로 편향되게 된다. 학문적으로는 '다중 검정 (multiple testing)을 시행하면 1종 오류가 발생할 확률이 증가하게 된다.'라고 표현한다.

 따라서, 각각을 비교해보는 것이 아니라 한꺼번에 비교하는 ANOVA를 시행해야 함이 마땅하다. 

 

여기에서 한 번 더 의문이 들 수 있다.

 "사후 분석을 할 때에는 1종 오류가 발생하지 않는가?"

그렇다. 1종 오류가 발생할 확률이 있으므로, p-value의 기준을 더 엄격하게 (0.05보다 더 작게) 잡아야 한다. P-value를 보정하는 방법은 여러 가지가 나와있는데, 그중 많이 쓰이는 네 가지를 언급한 것이다.

 

어떤 사후 분석을 쓸 것인가

 이 논의에 대해 정답이 따로 있는 것은 아니다. 적절한 방법을 사용하여 논문에 제시하면 되고, 어떤 것이 정답이라고 콕 집어 이야기할 수는 없다. 다만, 사후 분석 방법이 여러 가지가 있다는 것은 '사후 분석 방법에 따라 산출되는 결과가 달라질 수 있다.'는 것을 의미하고, 심지어는 '어떤 사후 분석 방법을 채택하냐에 따라 유의성 여부가 달라질 수도 있다.'는 것을 의미한다. 심지어, ANOVA에서는 유의한 결과가 나왔는데, 사후 분석을 해보니 유의한 차이를 보이는 경우가 없을 수도 있다. 따라서 어떤 사후 분석 방법 결과에 따른 결과인지 유의하여 해석할 필요가 있다.

 

 

코드 정리

####워킹 디렉토리 지정
setwd("C:/Users/user/Documents/Tistory_blog")


####데이터 불러오기
install.packages("readr")
library("readr")
df<-read_csv("Data.csv")


####데이터 나누기
df_ns<-df[df$SMOK==0,]
df_fs<-df[df$SMOK==1,]
df_cs<-df[df$SMOK==2,]


####정규성 검정
##비흡연자의 FVC 정규성 검정
# 1) Q-Q plot 그리기
qqnorm(df_ns$FVC_pPRED)
qqline(df_ns$FVC_pPRED)
# 2) 히스토그램 그리기
hist(df_ns$FVC_pPRED, prob=TRUE)
FVC_pPREDrange<-seq(min(df_ns$FVC_pPRED),max(df_ns$FVC_pPRED),length=max(max(df_ns$FVC_pPRED)-min(df_ns$FVC_pPRED),100))
ND<-dnorm(FVC_pPREDrange,mean=mean(df_ns$FVC_pPRED),sd=sd(df_ns$FVC_pPRED))
lines(FVC_pPREDrange, ND, lwd=2)
# 3) Shapiro-Wilk test 시행하기
shapiro.test(df_ns$FVC_pPRED)

##과거 흡연자의 FVC 정규성 검정
# 1) Q-Q plot 그리기
qqnorm(df_fs$FVC_pPRED)
qqline(df_fs$FVC_pPRED)
# 2) 히스토그램 그리기
hist(df_fs$FVC_pPRED, prob=TRUE)
FVC_pPREDrange<-seq(min(df_fs$FVC_pPRED),max(df_fs$FVC_pPRED),length=max(max(df_fs$FVC_pPRED)-min(df_fs$FVC_pPRED),100))
ND<-dnorm(FVC_pPREDrange,mean=mean(df_fs$FVC_pPRED),sd=sd(df_fs$FVC_pPRED))
lines(FVC_pPREDrange, ND, lwd=2)
# 3) Shapiro-Wilk test 시행하기
shapiro.test(df_fs$FVC_pPRED)

##과거 흡연자의 FVC 정규성 검정
# 1) Q-Q plot 그리기
qqnorm(df_cs$FVC_pPRED)
qqline(df_cs$FVC_pPRED)
# 2) 히스토그램 그리기
hist(df_cs$FVC_pPRED, prob=TRUE)
FVC_pPREDrange<-seq(min(df_cs$FVC_pPRED),max(df_cs$FVC_pPRED),length=max(max(df_cs$FVC_pPRED)-min(df_cs$FVC_pPRED),100))
ND<-dnorm(FVC_pPREDrange,mean=mean(df_cs$FVC_pPRED),sd=sd(df_cs$FVC_pPRED))
lines(FVC_pPREDrange, ND, lwd=2)
# 3) Shapiro-Wilk test 시행하기
shapiro.test(df_cs$FVC_pPRED)


####등분산 검정
#Levene의 등분산 검정
install.packages("lawstat")
library("lawstat")
levene.test(df$FVC_pPRED, df$SMOK, location="mean")

#Brown&Forsythe
install.packages("lawstat")
library("lawstat")
levene.test(df$FVC_pPRED, df$SMOK)

#Bartlett
bartlett.test(df$FVC_pPRED, df$SMOK)


####범주형 변수로의 변환
df$SMOK<-as.factor(df$SMOK)


####One-way ANOVA 시행
a<-aov(FVC_pPRED~SMOK, data=df)
summary(a)

####사후분석
#1) Bonferroni t test
pairwise.t.test(df$FVC_pPRED, df$SMOK, p.adj="bonferroni")

#2) Scheffe's multiple comparison (1)
install.packages("DescTools")
library("DescTools")
ScheffeTest(a)

#2-1) Scheffe's multiple comparison (2)
install.packages("agricolae")
library(agricolae)
scheffe.test(a,"SMOK",console=TRUE)


#3) Duncan's multiple range test 
install.packages("agricolae")
library(agricolae)
duncan.test(a,"SMOK",console=TRUE)


#4) Tukey's studentized range test 
TukeyHSD(a)

 

[R] 일원 배치 분산 분석 (One-way ANOVA, ANalysis Of VAriance) 정복 완료!

 

작성일: 2022.11.20.

최종 수정일: 2022.11.29.

이용 프로그램: R 4.2.2

RStudio v2022.07.2

RStudio 2022.07.2+576 "Spotted Wakerobin" Release

운영체제: Windows 10, Mac OS 12.6.1

반응형
반응형

 

[R] 등분산성 검정 (Homogeneity of variance) - levene.test(), bartlett.test()

 T-test, ANOVA 등 몇몇 분석에서는 분포들의 분산이 같다는 가정 (등분산성)이 필요하다. 이럴 때에는 등분산성 검정을 해야 하는데, 많이들 사용하는 방법은 크게 다섯 가지가 존재한다.

 

1) F test

2) Levene의 등분산 검정

3) O'Brien의 등분산 검정(R에서 검정 불가)

4) Brown and Forsythe의 등분산 검정

5) Bartlett의 등분산 검정

 

이 중 F test, Levene의 등분산성 검정은 T-test에 관한 글에서 다루었으므로 본 포스팅에서는 다루지 않겠다. Fooled F와 Levene 등분산성 검정에 관한 내용은 다음 링크에서 확인할 수 있다. 2022.11.12 - [모평균 검정/R] - [R] 독립 표본 T검정 (Independent samples T-test) - t.test(), var.test(), levene.test()

 

 결국 우리가 확인하고자 하는 것은, "수축기 혈압 분포의 분산이 고혈압 환자군과 일반인 사이에서 다르다고 할 수 있는가?"이다. 이에 대해 설면하도록 하겠다.

 

*실습용 데이터는 아래 링크를 클릭하면 다운로드할 수 있습니다.

2022.08.04 - [공지사항 및 소개] - 분석용 데이터 (update 22.10.11)

 

분석용 데이터 (update 22.10.11)

2022년 10월 11일 버전입니다. 변수는 계속하여 추가될 예정입니다. 다음 카테고리에 있는 글에서 이용된 데이터입니다. - 기술 통계 - 범주형 자료 분석 - 모평균 검정 - 반복 측정 자료 분석 - 통계

medistat.tistory.com

 

코드를 보여드리기에 앞서 워킹 디렉토리부터 지정하겠다.

워킹 디렉토리에 관한 설명은 다음 링크된 포스트에서 볼 수 있다.

2022.08.05 - [통계 프로그램 사용 방법/R] - [R] 작업 디렉토리 (Working Directory) 지정 - getwd(), setwd()

setwd("C:/Users/user/Documents/Tistory_blog")

 

데이터를 불러와 df에 객체로 저장하겠다.

데이터 불러오는 방법은 다음 링크에서 볼 수 있다.

2022.08.05 - [통계 프로그램 사용 방법/R] - [R] 데이터 불러오기 : EXCEL - read_excel(), read.xlsx()

2022.08.10 - [통계 프로그램 사용 방법/R] - [R] 데이터 저장하기 : CSV 파일 - write.csv(), write_csv()

2022.08.10 - [통계 프로그램 사용 방법/R] - [R] 데이터 불러오기 : SAS file (.sas7bdat) - read.sas7bdat(), read_sas()

 

install.packages("readr")
library("readr")
df<-read_csv("Data.csv")

 

코드 - 등분산성 검정

Brown&Forsythe 의 등분산 검정은 "levene.test()"라는 함수를 쓰는데, 이는 "lawstat"이라는 패키지에 있으므로 설치를 한다. Bartlett의 등분산 검정을 시행할 때 사용하는 "bartlett.test()"는 기본 내장 함수이므로 별도의 설치가 필요 없다.

설치에 관한 내용은 다음 링크에서 확인할 수 있다. 2022.08.05 - [통계 프로그램 사용 방법/R] - [R] 패키지 설치하기 - install.packages(), library()

 

install.packages("lawstat")
library("lawstat")

 

등분산성을 검정하는 코드는 다음과 같다.

#Brown&Forsythe
levene.test(df$SBP, df$HTN)

#Bartlett
bartlett.test(df$SBP, df$HTN)

#Brown&Forsythe
levene.test(df$SBP, df$HTN) : df 데이터의 SBP(수축기 혈압) 변수를 df 데이터의 HTN(고혈압 여부) 변수에 따라 Brown&Forsythe 등분산 검정을 시행하라.

#Bartlett
bartlett.test(df$SBP, df$HTN) : df 데이터의 SBP(수축기 혈압) 변수를 df 데이터의 HTN(고혈압 여부) 변수에 따라 Bartlett 등분산 검정을 시행하라.

 

결과 

	Modified robust Brown-Forsythe Levene-type test based on
	the absolute deviations from the median

data:  df$SBP
Test Statistic = 0.010831, p-value = 0.9171

	Bartlett test of homogeneity of variances

data:  df$SBP and df$HTN
Bartlett's K-squared = 3.5705e-07, df = 1, p-value =0.9995

Modified robust Brown-Forsythe Levene-type test based on
the absolute deviations from the median

data:  df$SBP
Test Statistic = 0.010831, p-value = 0.9171

Bartlett test of homogeneity of variances

data:  df$SBP and df$HTN
Bartlett's K-squared = 3.5705e-07, df = 1, p-value =0.9995

 

결과가 복잡해 보이지만 위에 색깔 처리해놓은 숫자만 봐도 충분하다.

Brown-Forsythe 검사 결과 p-value는 0.9171이고

Bartlett test 검사 결과 p-value는 0.9995이다.

 

즉 어떤 방법을 사용하더라도, 모두 p-value가 0.05 이상으로 귀무가설을 기각하지 못해 모분산이 서로 같다고 결론 내릴 수 있다. 

 

코드 정리

#워킹 디렉토리 지정
setwd("C:/Users/user/Documents/Tistory_blog")

#데이터 불러오기
install.packages("readr")
library("readr")
df<-read_csv("Data.csv")

#levene.test() 함수를 위한 패키지 설치
install.packages("lawstat")
library("lawstat")

#Brown&Forsythe 등분산 검정
levene.test(df$SBP, df$HTN)

#Bartlett 등분산 검정
bartlett.test(df$SBP, df$HTN)

 

[R] 등분산성 검정 (Homogeneity of variance: Brown&Forsythe, Bartlett test) 정복 완료!

 

작성일: 2022.11.20.

최종 수정일: 2022.11.20.

이용 프로그램: R 4.2.2

RStudio v2022.07.2

RStudio 2022.07.2+576 "Spotted Wakerobin" Release

운영체제: Windows 10, Mac OS 12.6.1

반응형
반응형

[R] 독립 표본 T검정 (Independent samples T-test) - t.test(),  var.test(),  levene.test()

  세상에 존재하는 모든 사람을 대상으로 연구를 하면 참 좋겠지만, 시간적인 이유로, 그리고 경제적인 이유로 일부를 뽑아서 연구를 진행할 수밖에 없다. 모든 사람을 모집단이라고 하고, 뽑힌 일부를 표본이라고 한다. 우리는 표본으로 시행한 연구로 모집단에 대한 결론을 도출해내고자 할 것이다. 

 1000명에게 피검사를 시행하였고, 간 기능 검사의 일환으로 ALT 수치를 모았다. 이 데이터를 기반으로 1000명이 기원한 모집단 인구에서의 ALT평균이 어떻게 될지 예측하는 것이 T-test이다. T-test는 크게 세 가지로 나눌 수 있다.

 

 1) 일표본 T검정 (One sample T-test) : 2022.11.03 - [모평균 검정/R] - [R] 일표본 T검정 (One-sample T-test) - t.test()

 : 모집단의 평균이 특정 값이라고 할 수 있는가?

 예) 모집단의 ALT 평균이 50이라고 할 수 있는가?

 

 

 2) 독립 표본 T검정 (Independent samples T-test, Two samples T-test)

 : 두 모집단의 평균이 다르다고 할 수 있는가?

 예) 고혈압 환자와 일반인의 수축기 혈압 평균이 서로 다르다고 할 수 있는가?

 

 3) 대응표본 T검정 (Paired samples T-test) : 2022.11.25 - [반복 측정 자료 분석/R] - [R] 대응 표본 T검정 (Paired samples T-test) - t.test()

 

 : 모집단의 짝지어진 변수들의 평균이 다르다고 할 수 있는가?

예) 간기능 개선제 복용 전 ALT 평균은 간기는 개선제 복용 후 ALT 평균과 다르다고 할 수 있는가?

 

이번 포스팅에서는 독립 표본 T검정 (Independent samples T-test, Two samples T-test)에 대해 알아볼 것이다.

 

*실습용 데이터는 아래 링크를 클릭하면 다운로드할 수 있습니다.

2022.08.04 - [공지사항 및 소개] - 분석용 데이터 (update 22.10.11)

 

분석용 데이터 (update 22.10.11)

2022년 10월 11일 버전입니다. 변수는 계속하여 추가될 예정입니다. 다음 카테고리에 있는 글에서 이용된 데이터입니다. - 기술 통계 - 범주형 자료 분석 - 모평균 검정 - 반복 측정 자료 분석 - 통계

medistat.tistory.com

 

코드를 보여드리기에 앞서 워킹 디렉토리부터 지정하겠다.

워킹 디렉토리에 관한 설명은 다음 링크된 포스트에서 볼 수 있다.

2022.08.05 - [통계 프로그램 사용 방법/R] - [R] 작업 디렉토리 (Working Directory) 지정 - getwd(), setwd()

setwd("C:/Users/user/Documents/Tistory_blog")

 

데이터를 불러와 df에 객체로 저장하겠다.

데이터 불러오는 방법은 다음 링크에서 볼 수 있다.

2022.08.05 - [통계 프로그램 사용 방법/R] - [R] 데이터 불러오기 : EXCEL - read_excel(), read.xlsx()

2022.08.10 - [통계 프로그램 사용 방법/R] - [R] 데이터 저장하기 : CSV 파일 - write.csv(), write_csv()

2022.08.10 - [통계 프로그램 사용 방법/R] - [R] 데이터 불러오기 : SAS file (.sas7bdat) - read.sas7bdat(), read_sas()

install.packages("readr")
library("readr")
df<-read_csv("Data.csv")

 

목표: 고혈압 환자의 수축기 혈압(SBP) 평균은 정상인의 수축기 혈압 평균과 다르다고 모집단 수준에서 말할 수 있는가?

이번 포스팅의 목적은 1000명의 데이터를 가지고, 이 1000명이 기원한 모집단에서 수축기 혈압 평균이 고혈압 유병 여부에 따라 다르다고 할 수 있는지 판단하는 것이다.

전제조건 (정규성) - 코드

 독립 표본 (Indepent samples, Two tamples) T 검정의 전제조건은 각각의 독립 표본의 분포가 정규성을 따른다는 것이다. 즉, 여기에서는 고혈압 유병 여부에 따라 수축기 혈압 (SBP)의 정규성을 검정하도록 하겠다.

 

따라서 다음 두 가지의 일을 해야 한다.

1) 고혈압 여부에 따라 데이터를 나누기

나누는 방법에 대한 설명은 다음 링크에서 볼 수 있다. 2022.11.10 - [통계 프로그램 사용 방법/R] - [R] 조건에 맞는 자료 추출하기

위 링크에서 확인할 수 있듯이 여러 가지 방법으로 나눌 수 있지만 indexing을 이용하여 나누도록 하겠다.

df_whtn<-df[df$HTN==1,]
df_wohtn<-df[df$HTN==0,]

 

2) 고혈압 여부에 따라 정규성을 검정하기

정규성 검정에 대한 설명은 다음 링크들에서 확인할 수 있다.

2022.08.11 - [기술 통계/R] - [R] 정규성 검정 (1) : Q-Q plot - qqnorm(), qqline()

2022.08.11 - [기술 통계/R] - [R] 정규성 검정 (2) : 히스토그램 - hist(), dnorm()

2022.08.11 - [기술 통계/R] - [R] 정규성 검정 (3) : 정량적 검정 (Shapiro-Wilk, Kolmogorov-Smirnov) - shapiro.test(), ks.test()

2022.08.12 - [기술 통계/R] - [R] 정규성 검정 (4) : 정량적 검정 (Lilliefors test) - lillie.test()

2022.08.16 - [기술 통계/R] - [R] 고급 Q-Q Plot - Van der Waerden, Rankit, Tukey, Blom

##고혈압 환자의 수축기 혈압 정규성 검정
# 1) Q-Q plot 그리기
qqnorm(df_whtn$SBP)
qqline(df_whtn$SBP)
# 2) 히스토그램 그리기
hist(df_whtn$SBP, prob=TRUE)
SBPrange<-seq(min(df_whtn$SBP),max(df_whtn$SBP),length=max(max(df_whtn$SBP)-min(df_whtn$SBP),100))
ND<-dnorm(SBPrange,mean=mean(df_whtn$SBP),sd=sd(df_whtn$SBP))
lines(SBPrange, ND, lwd=2)
# 3) Shapiro-Wilk test 시행하기
shapiro.test(df_whtn$SBP)

##정상인의 수축기 혈압 정규성 검정
# 1) Q-Q plot 그리기
qqnorm(df_wohtn$SBP)
qqline(df_wohtn$SBP)
# 2) 히스토그램 그리기
hist(df_wohtn$SBP, prob=TRUE)
SBPrange<-seq(min(df_wohtn$SBP),max(df_wohtn$SBP),length=max(max(df_wohtn$SBP)-min(df_wohtn$SBP),100))
ND<-dnorm(SBPrange,mean=mean(df_wohtn$SBP),sd=sd(df_wohtn$SBP))
lines(SBPrange, ND, lwd=2)
# 3) Shapiro-Wilk test 시행하기
shapiro.test(df_wohtn$SBP)

 

정규성 검정 - 결과

1) 고혈압 환자

 

	Shapiro-Wilk normality test

data:  df_whtn$SBP
W = 0.9975, p-value = 0.6677

 N수가 2,000개 미만이므로 Shapiro-Wilk 통계량의 p-value를 보면 0.05 이상이며, Q-Q Plot은 대부분의 데이터가 선상에 있고, 히스토그램에서도 정규성을 따르는 것처럼 볼 수 있다. 

 

 

2) 정상인

	Shapiro-Wilk normality test

data:  df_wohtn$SBP
W = 0.99676, p-value = 0.4119

 N수가 2,000개 미만이므로 Shapiro-Wilk 통계량의 p-value를 보면 0.05 이상이며, Q-Q Plot은 대부분의 데이터가 선상에 있고, 히스토그램에서도 정규성을 따르는 것처럼 볼 수 있다. 

 

 전제조건이 성립한다. 즉, 고혈압 유병 여부에 따른 수축기 혈압 (SBP)의 분포가 정규성을 따른다고 할 수 있다. 따라서 독립 표본 T 검정 (Independent samples T test, Two samples T test)을 시행할 수 있다.

 

등분산성 검정

 독립 표본 T 검정 (Independent samples T test, Two samples T test)은 고혈압 유병 여부에 따른 수축기 혈압의 분포의 분산이 같은지 여부에 따라 시행하는 검사 방법이 다르다. 따라서 분산이 같은지 확인하는 등분산성 검정을 먼저 시행해야 한다.

 

많이 사용하는 등분산성 검정 방법은 다음 다섯 가지가 있다.

1) F test

2) Levene의 등분산 검정

3) O'Brien의 등분산 검정(R에서 검정 불가)

4) Brown and Forsythe의 등분산 검정

5) Bartlett의 등분산 검정

 

이 글에서는 F test와 Levene의 등분산 검정까지만 해보고, 나머지는 다음 링크에서 다루기로 한다. 2022.11.20 - [모평균 검정/R] - [R] 등분산성 검정 (Homogeneity of variance) - levene.test(), bartlett.test()

 

#F test
var.test(SBP~HTN, data=df)

#Levene의 등분산 검정
install.packages("lawstat")
library("lawstat")
levene.test(df$SBP, df$HTN, location="mean")

#F test
var.test(SBP~HTN, data=df) : 등분산 검정을 하는데, HTN에 따라 SBP의 분포의 분산이 같은지 검정하라. 데이터는 df를 사용한다.

#Levene의 등분산 검정

install.packages("lawstat") : Levene의 등분산 검정을 하기 위해서는 lawstat이라는 패키지를 설치해야 한다.
library("lawstat") : lawstat패키지를 불러온다.
levene.test(df$SBP, df$HTN, location="mean") : df 데이터의 HTN에 따라 df 데이터의 SBP의 분포의 분산이 같은지 검정하라. 이때 [location="mean"]을 적어주어야 Levene의 등분산 검정을 시행한다.

 

결과

	F test to compare two variances

data:  SBP by HTN
F = 0.99995, num df = 502, denom df = 496, p-value =
0.9994
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
 0.8387988 1.1919414
sample estimates:
ratio of variances 
         0.9999465 
         
         
	Classical Levene's test based on the absolute deviations
	from the mean ( none not applied because the location is
	not set to median )

data:  df$SBP
Test Statistic = 0.010363, p-value = 0.9189

F test to compare two variances
data:  SBP by HTN
F = 0.99995, num df = 502, denom df = 496, p-value =0.9994
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
 0.8387988 1.1919414
sample estimates:
ratio of variances 
         0.9999465 
         
         
Classical Levene's test based on the absolute deviations from the mean ( none not applied because the location is not set to median )
data:  df$SBP
Test Statistic = 0.010363, p-value = 0.9189

 

결과 창에 나오는 것들이 많아 복잡해 보이지만 위에 파란색 볼드로 처리해놓은 내용만 보면 된다. F test의 p-value는 0.9994, Levene의 등분산 검정의 p-value는 0.9189이다. 모두 0.05보다 크므로 귀무가설을 기각할 수 없으며 귀무가설을 택해야 한다. 여기에서의 귀무가설은 "HTN에 따라 SBP의 분포의 분산에 차이가 없다."이므로 분산이 같다는 결론을 내린다.

 

이제 드디어 사전 점검이 모두 끝났고, t-test를 시행할 수 있다. 

코드

##T-test 
t.test(df_wohtn$SBP, df_whtn$SBP, var.equal=TRUE)


t.test(df_wohtn$SBP, df_whtn$SBP, var.equal=TRUE) : 고혈압이 없는 사람들의 수축기 혈압 (df_wohtn$SBP)과 고혈압이 있는 사람들의 수축기 혈압 (df_whtn$SBP)으로 t-test를 시행하라. 대신 두 군의 분산은 같다고 가정한다.

 

만약, 등분산성 검정 결과 분산이 다르다고 결론이 났으면 "var.equal=FALSE" 라고 적으면 된다. 혹은 t.test의 var.equal은 FALSE를 기본값으로 가지므로, "var.equal=TRUE"를 그냥 지워도 된다. 

 

결과

	Two Sample t-test

data:  df_wohtn$SBP and df_whtn$SBP
t = -46.666, df = 998, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -31.19793 -28.68000
sample estimates:
mean of x mean of y 
 120.1473  150.0862

Two Sample t-test

data:  df_wohtn$SBP and df_whtn$SBP

t = -46.666, df = 998, p-value < 2.2e-16

alternative hypothesis: true difference in means is not equal to 0

95 percent confidence interval:  -31.19793 -28.68000

sample estimates:

mean of x mean of y   

120.1473  150.0862 

 

결과가 복잡한 듯 하지만 파란색 글씨만 읽으면 된다. p-value는 $2.2 \times 10^{-16}$보다 작다. 따라서 귀무가설을 기각할 수 있으며, 본 검정에서 귀무가설은 "HTN에 따라 SBP의 평균이 같다."이므로 "HTN에 따라 SBP의 평균은 다르다."라고 결론 내릴 수 있다. 마지막 행을 보면 "mean of x"은 120.1473, "mean of y"은 150.0862라고 알려주고 있는데, x는 코드에 처음 들어간 "df_wohtn$SBP"를 의미하고, y는 두 번째로 들어간 "df_whtn$SBP"를 의미한다. 즉 고혈압이 없는 군의 수축기 혈압 평균은 120.1473, 고혈압이 없는 군의 수축기 혈압 평균은 150.0862다. 

 

코드 정리

#워킹 디렉토리 지정
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$SBP)
qqline(df_whtn$SBP)
# 2) 히스토그램 그리기
hist(df_whtn$SBP, prob=TRUE)
SBPrange<-seq(min(df_whtn$SBP),max(df_whtn$SBP),length=max(max(df_whtn$SBP)-min(df_whtn$SBP),100))
ND<-dnorm(SBPrange,mean=mean(df_whtn$SBP),sd=sd(df_whtn$SBP))
lines(SBPrange, ND, lwd=2)
# 3) Shapiro-Wilk test 시행하기
shapiro.test(df_whtn$SBP)

##정상인의 수축기 혈압 정규성 검정
# 1) Q-Q plot 그리기
qqnorm(df_wohtn$SBP)
qqline(df_wohtn$SBP)
# 2) 히스토그램 그리기
hist(df_wohtn$SBP, prob=TRUE)
SBPrange<-seq(min(df_wohtn$SBP),max(df_wohtn$SBP),length=max(max(df_wohtn$SBP)-min(df_wohtn$SBP),100))
ND<-dnorm(SBPrange,mean=mean(df_wohtn$SBP),sd=sd(df_wohtn$SBP))
lines(SBPrange, ND, lwd=2)
# 3) Shapiro-Wilk test 시행하기
shapiro.test(df_wohtn$SBP)

##등분산성 검정
#F test
var.test(SBP~HTN, data=df)

#Levene의 등분산 검정
install.packages("lawstat")
library("lawstat")
levene.test(df$SBP, df$HTN, location="mean")

##T-test 
t.test(df_wohtn$SBP, df_whtn$SBP, var.equal=TRUE)

 

R 독립 표본 T검정 (Dependent samples T-test, Two samples T-test) 정복 완료!

작성일: 2022.11.12.

최종 수정일: 2022.11.30.

이용 프로그램: R 4.2.2

RStudio v2022.07.2

RStudio 2022.07.2+576 "Spotted Wakerobin" Release

운영체제: Windows 10, Mac OS 12.6.1

반응형

+ Recent posts