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

[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.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     
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

 

빨간색 글씨가 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