반응형

[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

반응형
반응형

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

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

 

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.04 - [통계 프로그램 사용 방법/SPSS] - [SPSS] 데이터 불러오기 및 저장하기

 

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

 

전제조건 (정규성)

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

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

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

2022.08.11 - [기술 통계/SPSS] - [SPSS] 정규성 검정

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

 

1) 정규성

1) 분석(A) > 기술통계량(E) > 데이터 탐색 (E)

 

2) 분석하고자 하는 변수인 QOL을 "종속변수"에 넣어준다. 거주 지역 단위에 따라 정규성을 검정할 것이므로 "요인(F)"에 RESID를 넣는다. 그 뒤 "도표(T)..."를 선택한다.

 

3) "히스토그램(H)", "검정과 함께 정규성 도표(O)" 체크박스를 클릭하고 "계속(C)"를 누르고, 돌아가 "확인"을 누른다.

 

결과

1) Q-Q plot

 

2) 히스토그램

 

3) Shapiro-Wilk 검정

 

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

 

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

 

Welch's & Brown-Forsythe ANOVA 실행 방법

1) 분석(A) > 평균 비교 (M) > 일원배치 분산분석(O)

 

2) 분석하고자 하는 변수 QOL를 "종속변수(E)"쪽으로 옮기고, RESID을 "요인(F)"쪽으로 옮긴다. 그리고 "옵션"을 누른다.

 

3) "분산 동질성 검정(H)" 체크박스를 선택한다. 분산이 같지 않을 경우를 대비해 "Brown-Forsythe", "Welch" 체크박스를 클릭한다.  그리고 "계속 (C)"을 누른다.

 

4)  "사후분석(H)"을 클릭한다.

분산이 같지 않을 경우에만 대비하여 "Tanhame의 T2", "Dunnett의 T3", "Games-Howell", "Dunnett의 C"를 클릭한다. 분산이 같지 않은 경우 위쪽에 있는 사후 분석 방법은 사용할 수 없고 아래 사후 분석만 사용할 수 있다.

 

 

결과

 1) 등분산성 검정

 

맨 위에 평균을 기준으로 한 검정이 일반적인 Levene의 등분산 검정이다.

그다음 중위수를 기준으로 한 검정이 Brown-Forsythe의 등분산성 검정이다.

어떤 것을 선택하는지는 개인의 선택이다. 하지만 이 경우 둘 다 0.05보다 작으므로 귀무가설을 기각하고, 분산에 차이가 있다고 결론 내린다.

 

2) 일반 ANOVA

.

분산이 다르다고 결론 내렸으므로 이 표는 보면 안 된다.

 

3) Welch, Brown-Forsythe ANOVA

두 가지 방법 모두 p-value가 0.001보다 작다고 이야기하고 있다. 따라서 귀무가설을 기각하고 대립 가설을 채택한다. 그렇다면 여기에서 귀무가설 및 대립 가설은 무엇이었는가?

 

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

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

 

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

 

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

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

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

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

 

사후 분석 결과를 볼 것이다.

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

 

1) Tamhane's T2

2) Games-Howell

3) Dunnett's T3

4) Dunnett's C

 

1) Tamhane's T2

2) Games-Howell

3) Dunnett's T3

4) Dunnett's C

 

박스 색깔과 검정방법을 매칭 시켜두었다.

위 표를 보는 법은 다음과 같다.

0 1 대도시 거주자 (0)와 중소도시 거주자 (1) 사이의 비교
2 대도시 거주자 (0)와 시골 거주자 (2) 사이의 비교
1 0 대도시 거주자 (0)와 중소도시 거주자 (1) 사이의 비교
2 중소도시 거주자 (1)와 시골 거주자 (2) 사이의 비교
2 0 대도시 거주자 (0)와 시골 거주자 (2) 사이의 비교
1 중소도시 거주자 (1)와 시골 거주자 (2) 사이의 비교

Dunnett C를 제외하고, 어떤 방법을 사용하든, 모든 비교에서 p-value는 0.001보다 작아 매우 유의하다. 그리고 Dunnett C는 신뢰구간에 0이 포함되지 않으면 유의한 것인데, 0을 포함하는 구간이 없다. 따라서 모든 비교에서 유의미한 FVC의 차이가 있다고 할 수 있다.

 

어떤 사후 분석을 쓸 것인가

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

 

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

작성일: 2022.11.30.

최종 수정일: 2022.11.30.

이용 프로그램: IBM SPSS v26

운영체제: Windows 10

반응형
반응형

[SPSS] 일원 배치 분산 분석 (One-way ANOVA, ANalysis Of VAriance)

 

이전 포스팅에서 음주 여부에 따른 ALT 평균에 차이가 있는지를 알아보기 위해서는 T test를 사용한다는 내용을 알아보았다. 2022.11.30 - [모평균 검정/SPSS] - [SPSS] 독립 표본 T검정 (Independent samples T-test)

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

 

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

 

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

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

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

 

분석용 데이터 (update 22.11.28)

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

medistat.tistory.com

 

데이터를 불러오도록 한다. 불러오는 방법은 다음 링크를 확인하도록 한다.

2022.08.04 - [통계 프로그램 사용 방법/SPSS] - [SPSS] 데이터 불러오기 및 저장하기

 

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

 

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

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

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

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

 

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

 

1) 정규성 

2022.08.11 - [기술 통계/SPSS] - [SPSS] 정규성 검정

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

 

1) 분석(A) > 기술통계량(E) > 데이터 탐색 (E)

 

2) 분석하고자 하는 변수인 FVC_pPRED을 "종속변수"에 넣어준다. 흡연 상태에 따라 분석할 것이므로 "요인(F)"에 SMOK을 넣는다. 그 뒤 "도표(T)..."를 선택한다.

 

3) "히스토그램(H)", "검정과 함께 정규성 도표(O)" 체크박스를 클릭하고 "계속(C)"를 누르고, 돌아가 "확인"을 누른다.

 

결과
1) Q-Q plot

 

2) 히스토그램

 

3) Shapiro-Wilk 검정

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

 

 따라서 정규성 전제조건은 만족한다고 할 수 있다.

 

일원 배치 분산 분석 (ANOVA) 방법

1) 분석(A) > 평균 비교 (M) > 일원배치 분산분석(O)

 

2) 분석하고자 하는 변수 FVC_pPRED를 "종속변수(E)"쪽으로 옮기고, SMOK을 "요인(F)"쪽으로 옮긴다. 그리고 "옵션"을 누른다.

 

3) "분산 동질성 검정(H)" 체크박스를 선택한다. 그리고 "계속 (C)"을 누른다.

 

4)  "사후분석(H)"을 클릭한다.

많은 사후분석 방법이 있지만 우리는 Bonferroni, Scheffe, Tukey, Duncan을 확인하도록 한다.  이 네 개를 체크한다.

 

 

결과

1) 등분산성 검정

ANOVA의 두 번째 전제조건인 등분산성에 대한 검정이다.

맨 위에 평균을 기준으로 한 검정이 일반적인 Levene의 등분산 검정이다.

그다음 중위수를 기준으로 한 검정이 Brown-Forsythe의 등분산성 검정이다.

어떤 것을 선택하는지는 개인의 선택이다. 하지만 이 경우 둘 다 0.05보다 크므로 귀무가설을 받아들여, 분산에 차이가 없다고 결론 내린다. 따라서 전제조건을 만족함을 알 수 있다.

 

2) ANOVA결과

ANOVA의 결과는 집단-간 유의 확률을 확인하면 된다. p-value<0.001으로 유의하다. 따라서 귀무가설을 기각하고 대립 가설을 선택한다. 그렇다면 여기에서 귀무가설 및 대립 가설은 무엇이었는가?

 

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

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

 

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

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

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

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

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

 

사후 분석 결과 부분을 보겠다.

 

1) 사후검정 1 - Tukey, Scheffe, Bonferroni

Duncan을 제외한 방법은 각 차이에 대한 보정된 p-value를 보여준다. 위 표를 보는 법은 다음과 같다.

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

어떤 방법을 사용하든, 모든 비교에서 p-value는 0.001보다 작아 매우 유의하므로 모든 비교에서 유의미한 FVC의 차이가 있다고 할 수 있다.

 

2) 사후검정 2 - Tukey, Duncan, Scheffe

Bonferroni는 동질적 부분집합으로 표현해주지 않으며, Duncan은 동질적 부분집합으로만 표현해준다. 

위 세 개 모두 같은 결과를 보여주고 있다. 해석하는 방법은 다음과 같다.

SMOK=2는 평균이 76.78562788이며, 부분집합 1에 속한다.

SMOK=1는 평균이 85.80627187이며, 부분집합 2에 속한다.

SMOK=0는 평균이 95.33911151이며, 부분집합 3에 속한다.

속한 부분집합이 다르면 평균에 유의미한 차이가 있다는 것이다. 부분집합 1, 2, 3은 서로 다른 평균을 가진 모집단에서 유래했다고 볼 수 있으므로 비흡연자, 과거 흡연자, 현재 흡연자는 서로 다른 평균을 갖고 있다고 할 수 있다. 

 

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

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

 아주 논리적인 의문점이다. 하지만, 반드시 ANOVA를 사용해야 한다. 그 이유는 다음과 같다. 본 사례는 흡연 상태에 따른 조합 가능한 경우의 수가 3인데, 각각 유의성의 기준을 $\alpha=0.05$로 잡아보자. 이때 세 번의 비교에서 모두 귀무가설이 학문적인 진실인데(평균에 차이가 없다.), 세 번 모두 귀무가설을 선택할 확률은 $\left(1-0.05 \right)^3 \approx 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를 보정하는 방법은 여러 가지가 나와있는데, 가장 간단하고 제일 보수적인 방법이 위에 예시에서 사용한 본페로니 방법이다. 그 외에도 굉장히 많은 방법이 있지만, Bonferroni, Duncan, Scheffe, Tukey가 가장 많이 쓰이는 듯하다. 이 중에 Duncan, Tukey는 그룹별로 표본의 수가 같을 때만 사용해야 하므로 유의하도록 한다.

 

어떤 사후 분석을 쓸 것인가

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

 

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

 

작성일: 2022.11.30.

최종 수정일: 2022.11.30.

이용 프로그램: IBM SPSS v26

운영체제: Windows 10

반응형
반응형

[SPSS] 독립 표본 T검정 (Independent samples T-test)

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

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

 

 1) 일표본 T검정 (One sample T-test) : 2022.11.29 - [모평균 검정/SPSS] - [SPSS] 일표본 T검정 (One-sample T-test)

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

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

 

 

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

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

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

 

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

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

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

 

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

 

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

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

 

분석용 데이터 (update 22.11.28)

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

medistat.tistory.com

 

데이터를 불러오도록 한다. 불러오는 방법은 다음 링크를 확인하도록 한다.

2022.08.04 - [통계 프로그램 사용 방법/SPSS] - [SPSS] 데이터 불러오기 및 저장하기

 

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

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

전제조건 (정규성) 

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

2022.08.11 - [기술 통계/SPSS] - [SPSS] 정규성 검정

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

 

1) 분석(A) > 기술통계량(E) > 데이터 탐색 (E)

 

2) 분석하고자 하는 변수인 SBP을 "종속변수"에 넣고, 고혈압에 따라 분석할 것이므로 HTN을 "요인(F)"에 넣는다. 그 뒤 "도표(T)..."를 선택한다.

 

3) "히스토그램(H)", "검정과 함께 정규성 도표(O)" 체크박스를 클릭하고 "계속(C)"를 누르고, 돌아가 "확인"을 누른다.

 

정규성 검정 - 결과 (QQ plot과 Histogram)

1) 정상인

 

2) 고혈압 환자

 

정규성 검정 - 결과 (Shapiro-Wilk 검정)

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

 

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

 

 

독립 표본 T 검정 방법

1) 분석(A) > 평균 비교 (M) > 독립표본 T 검정

 

2) 분석하고자 하는 변수 SBP를 "검정 변수(T)"쪽으로 옮기고, 아래 "집단변수(G)"에 HTN을 넣는다. 그리고 "집단 정의(D)"을 누른다.

 

3) HTN은 정상인을 0으로, 고혈압 환자를 1로 정의하고 있다. 즉 우리는 0과 1을 비교할 것이므로 집단 1과 집단 2에 각각 0과 1을 넣는다.

 

결과

 1) 기술 통계

먼저 기술통계를 보여준다. 고혈압 유무에 따른 표본 수, 평균, 표준편차, 표준오차를 보여준다.

 

 

2) T 검정 결과

표를 하나씩 뜯어보기로 한다.

 

 (1) 등분산 검정

 Levene의 등분산 검정을 시행해준다. p-value는 0.919로 0.05보다 크다. 따라서 귀무가설을 기각할 수 없다. 이 검정의 귀무가설과 대립 가설은 무엇일까?

귀무가설: 고혈압 환자와 정상인의 SBP분포는 등산이 같다.

귀무가설: 고혈압 환자와 정상인의 SBP분포는 등산이 다르다.

 

따라서 우리는 분산이 같다고 가정한다.

 

(2) 평균 차이 검정

등분산이 가정되므로 파란색 동그라미 안의 p-value를 읽는다. 만약 등분산성이 보장되지 않으면 그 아래의 p-value를 읽는다. p-value는 0.001보다 작으므로 귀무가설을 기각할 수 있으며, 본 검정에서 귀무가설은 "HTN에 따라 SBP의 평균이 같다."이므로 "HTN에 따라 SBP의 평균은 다르다."라고 결론 내릴 수 있다.

 

[SPSS] 독립 표본 T검정 (Independent samples T-test) 정복 완료!

 

작성일: 2022.11.30.

최종 수정일: 2022.11.30.

이용 프로그램: IBM SPSS v26

운영체제: Windows 10

반응형
반응형

[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

반응형
반응형

[SPSS] 일표본 T검정 (One-sample T-test)

 

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

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

 

 1) 일표본 T검정 (One sample T-test)

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

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

 

 2) 독립 표본 T검정 (Independent samples T-test): 2022.11.30 - [모평균 검정/SPSS] - [SPSS] 독립 표본 T검정 (Independent samples T-test)

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

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

 

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

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

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

 

이번 포스팅에서는 일표본 T검정 (One sample T-test)에 대해 알아볼 것이다.

 

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

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

 

분석용 데이터 (update 22.11.28)

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

medistat.tistory.com

 

데이터를 불러오도록 한다. 불러오는 방법은 다음 링크를 확인하도록 한다.

2022.08.04 - [통계 프로그램 사용 방법/SPSS] - [SPSS] 데이터 불러오기 및 저장하기

 

목표: 모집단의 ALT 평균이 50이라고 할 수 있는가?

 이번 포스팅의 목적은 1000명의 데이터를 가지고, 이 1000명이 기원한 모집단의 ALT평균이 50이라고 할 수 있는지 판단하는 것이다.

전제조건 (정규성)

일표본 (One sample) T 검정의 전제조건은 해당 분포가 정규성을 따른다는 것이다. 정규성 검정을 하도록 하겠다. 정규성 검정에 관한 분석 내용은 다음 글에서 살펴볼 수 있다. 

2022.08.11 - [기술 통계/SPSS] - [SPSS] 정규성 검정

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

 

1) 분석(A) > 기술통계량(E) > 데이터 탐색 (E)

 

2) 분석하고자 하는 변수인 ALT을 "종속변수"에 넣어준다. 그 뒤 "도표(T)..."를 선택한다.

 

3) "히스토그램(H)", "검정과 함께 정규성 도표(O)" 체크박스를 클릭하고 "계속(C)"를 누르고, 돌아가 "확인"을 누른다.

1) Q-Q plot

2) 히스토그램

3) Shapiro-Wilk 검정

기술통계도 보여주는데 평균은 35정도임을 알 수 있다.

 

 N수가 2,000개 미만이므로 Shapiro-Wilk 통계량의 p-value를 보면 0.05 이상이며, Q-Q Plot은 대부분의 데이터가 선상에 있고, 히스토그램에서도 정규성을 따르는 것처럼 볼 수 있다. 따라서 ALT로는 일표본 T검정 (One-sample T-test)를 시행할 수 있다. 시행하도록 하겠다.

 

 

 

 

1) 분석(A) > 평균 비교 (M) > 일표본 T 검정 (S)

 

2) 분석하고자 하는 변수 ALT를 "검정 변수(T)"쪽으로 옮기고, 아래 "검정값(V)"에 50을 적는다. 그리고 "확인"을 누른다.

 

결과

확인할 것은 "유의확률(양측)"의 값인 0.000이다. 이는 0.001보다 작다는 뜻으로 매우 유의하다는 말이다. 따라서 귀무가설을 기각하고 대립 가설을 선택할 것이다. 여기에서 귀무가설은 "모집단의 ALT 평균은 50이다."였고, 대립가설은 "모집단의 평균은 50보다 크거나 작다."이므로 이 집단의 모평균이 50이 아닌 것으로 결론 내리는 것이 옳다. 평균이 35였으므로 합당한 결론이다.

 

SPSS 일표본 T검정 (One-sample T-test) 정복 완료!

 

작성일: 2022.11.29.

최종 수정일: 2022.11.30.

이용 프로그램: IBM SPSS v26

운영체제: Windows 10

반응형
반응형

[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

반응형
반응형

[SAS] 등분산 가정이 성립하지 않는 일원 배치 분산 분석 (Welch's ANOVA) - PROC GLM, PROC MIXED

 

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

 

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

1) 정규성

2) 등분산성

 

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

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

 

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

 

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

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

 

분석용 데이터 (update 22.11.26)

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

medistat.tistory.com

 

시작하기 위해 라이브러리를 만들고, 파일을 불러온다.

라이브러리 만드는 방법: 2022.08.05 - [통계 프로그램 사용 방법/SAS] - [SAS] 라이브러리 만들기 - LIBNAME

파일 불러오는 방법: 2022.08.05 - [통계 프로그램 사용 방법/SAS] - [SAS] 데이터 불러오기 및 저장하기 - PROC IMPORT, PROC EXPORT

*라이브러리 지정하기;
LIBNAME hong "C:/Users/User/Documents/Tistory_blog";

*파일 불러오기;
PROC IMPORT
DATAFILE="C:\Users\user\Documents\Tistory_blog\Data.xlsx"
DBMS=EXCEL
OUT=hong.df
REPLACE;
RUN;

 

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

 

전제조건 (정규성)

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

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

 

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

 

1) 정규성 

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

2022.08.12 - [기술 통계/SAS] - [SAS] 정규성 검정 - PROC UNIVARIATE

 

PROC UNIVARIATE DATA=hong.df NORMAL PLOT;
CLASS RESID;
VAR QOL;
HISTOGRAM QOL/ NORMAL (MU=EST SIGMA=EST);
RUN;

 

결과

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

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

 

 

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

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

 

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

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

 

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

 

 

2) 등분산성

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

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

2022.10.04 - [모평균 검정/SAS] - [SAS] 독립 표본 T검정 (Independent samples T-test) - PROC TTEST

2022.10.05 - [모평균 검정/SAS] - [SAS] 등분산성 검정 (Homogeneity of variance) - PROC GLM

 

위 두 개 글에서 제시한 네 가지 방법 모두를 사용하여 등분산 검정을 시행해 보겠다.

*1) Levene;
PROC GLM DATA=hong.df;
CLASS RESID;
MODEL QOL=RESID;
MEANS RESID/ HOVTEST=LEVENE(TYPE=ABS);
RUN;

*2) O'Brien;
PROC GLM DATA=hong.df;
CLASS RESID;
MODEL QOL=RESID;
MEANS RESID/ HOVTEST=OBRIEN;
RUN;

*3) Brwon and Forsythe;
PROC GLM DATA=hong.df;
CLASS RESID;
MODEL QOL=RESID;
MEANS RESID/ HOVTEST=BF;
RUN;

*4) Bartlett;
PROC GLM DATA=hong.df;
CLASS RESID;
MODEL QOL=RESID;
MEANS RESID/ HOVTEST=BARTLETT;
RUN;
QUIT;

 

결과

1) Levene

2) O'Brien

3) Brown and Forsythe

4) Bartlett

 어떤 방법을 사용하더라도 p-value<0.0001로 등분산성은 만족하지 않는다는 것을 알 수 있다. 따라서 일원 배치 분산 분석 (One-way ANOVA)이 아닌 Welch's ANOVA를 시행해야 한다.

 

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

PROC GLM DATA=hong.df;
CLASS RESID;
MODEL QOL=RESID;
MEANS RESID/ WELCH;
RUN;
QUIT;

PROC GLM DATA=hong.df; : Welch's ANOVA를 시행할 수 있는 PROC GLM으로 코드를 시작한다. 데이터는 hong 라이브러리의 df데이터를 활용한다.
CLASS RESID; : RESID에 따른 평균을 비교할 것이다.
MODEL QOL=RESID; : RESID에 따른 QOL의 평균을 비교할 것이다.

MEANS RESID/ WELCH; : RESID의 평균을 보여주고, Welch's ANOVA를 추가로 시행하라.

결과

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

 

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

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

 

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

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

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

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

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

 

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

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

 

1) Tamhane's T2

2) Games-Howell

3) Dunnett's T3

4) Dunnett's C

 

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

 

사후 분석 코드

*1) Tamhane's T2;

PROC MIXED DATA=hong.df;
CLASS RESID;
MODEL QOL=RESID/DDFM=SATTERTH;
REPEATED/GROUP=RESID;
LSMEANS RESID/ADJUST=SIDAK ADJDFE=ROW;
RUN;

*2) Games-Howell;

PROC MIXED DATA=hong.df;
CLASS RESID;
MODEL QOL=RESID/DDFM=SATTERTH;
REPEATED/GROUP=RESID;
LSMEANS RESID/ADJUST=TUKEY ADJDFE=ROW;
RUN;

*3) Dunnett's T3;

PROC MIXED DATA=hong.df;
CLASS RESID;
MODEL QOL=RESID/DDFM=SATTERTH;
REPEATED/GROUP=RESID;
LSMEANS RESID/ADJUST=SMM ADJDFE=ROW;
RUN;

PROC MIXED DATA=hong.df; : 사후 분석을 시행하는 PROC MIXED로 시작하고, 데이터는 hong 라이브러리의 df 데이터를 사용한다.
CLASS RESID; : RESID에 따라 사후 분석을 시행한다.
MODEL QOL=RESID/DDFM=SATTERTH; : RESID에 따른 QOL 평균의 비교를 한다. (DDFM=SATTERTH는 자유도를 근사시킬 때 사용하는 옵션이다.)
REPEATED/GROUP=RESID; : 이는 R matrix 추정에 관한 내용이다.
LSMEANS RESID/ADJUST=SIDAK ADJDFE=ROW; : "ADJUST=SIDAK"과 "ADJDFE=ROW"의 조합은 Tamhane's T2 사후 검정을 하도록 지시한다.
"ADJUST="에 "TUKEY"를 넣으면 Games-Howell, "SMM"을 넣으면 Dunnett's T3검정을 시행한다. 

 

 

1) Tamhane's T2

2) Games-Howell

3) Dunnett's T3

봐야 할 것은 "Adj P", 즉 보정된 p-value이다. 방법에 관계없이 p-value는 모두 0.0001보다 작다. 따라서 다음과 같이 결론 내릴 수 있다.

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

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

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

 

어떤 사후 분석을 쓸 것인가

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

 

코드 정리

*라이브러리 지정하기;
LIBNAME hong "C:/Users/User/Documents/Tistory_blog";

*파일 불러오기;
PROC IMPORT
DATAFILE="C:\Users\user\Documents\Tistory_blog\Data.xlsx"
DBMS=EXCEL
OUT=hong.df
REPLACE;
RUN;

*정규성 검정하기;
PROC UNIVARIATE DATA=hong.df NORMAL PLOT;
CLASS RESID;
VAR QOL;
HISTOGRAM QOL/ NORMAL (MU=EST SIGMA=EST);
RUN;

*등분산성 검정하기;
*1) Levene;
PROC GLM DATA=hong.df;
CLASS RESID;
MODEL QOL=RESID;
MEANS RESID/ HOVTEST=LEVENE(TYPE=ABS);
RUN;

*2) O'Brien;
PROC GLM DATA=hong.df;
CLASS RESID;
MODEL QOL=RESID;
MEANS RESID/ HOVTEST=OBRIEN;
RUN;

*3) Brwon and Forsythe;
PROC GLM DATA=hong.df;
CLASS RESID;
MODEL QOL=RESID;
MEANS RESID/ HOVTEST=BF;
RUN;

*4) Bartlett;
PROC GLM DATA=hong.df;
CLASS RESID;
MODEL QOL=RESID;
MEANS RESID/ HOVTEST=BARTLETT;
RUN;
QUIT;

*Welch's ANOVA 시행하기;
PROC GLM DATA=hong.df;
CLASS RESID;
MODEL QOL=RESID;
MEANS RESID/ WELCH;
RUN;
QUIT;

*사후분석 시행하기;
*1) Tamhane's T2;

PROC MIXED DATA=hong.df;
CLASS RESID;
MODEL QOL=RESID/DDFM=SATTERTH;
REPEATED/GROUP=RESID;
LSMEANS RESID/ADJUST=SIDAK ADJDFE=ROW;
RUN;

*2) Games-Howell;

PROC MIXED DATA=hong.df;
CLASS RESID;
MODEL QOL=RESID/DDFM=SATTERTH;
REPEATED/GROUP=RESID;
LSMEANS RESID/ADJUST=TUKEY ADJDFE=ROW;
RUN;

*3) Dunnett's T3;

PROC MIXED DATA=hong.df;
CLASS RESID;
MODEL QOL=RESID/DDFM=SATTERTH;
REPEATED/GROUP=RESID;
LSMEANS RESID/ADJUST=SMM ADJDFE=ROW;
RUN;

 

사실 등분산성이 만족하지 않을 때 시행할 수 있는 ANOVA로 Welch's ANOVA 외에도 Brown and Forsythe ANOVA가 있다. 이것과 Dunnett's C 사후 검정을 SPSS에서는 시행할 수 있으므로 필요한 경우 SPSS에서 시행해볼 수도 있겠다.

 

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

 

작성일: 2022.11.26.

최종 수정일: 2022.11.26.

이용 프로그램: SAS v9.4

운영체제: Windows 10

반응형

+ Recent posts