<통기타> #2. 짝을 이룬 자료의 비교

들어가며

통기타 시리즈의 대부분의 내용은 

김용은 / 김지형 저자님의 <한눈에 쏙쏙 의학통계 배우기> 3판의 흐름을 참고한다.

R의 경우 코드와 샘플데이터는 김종엽 저자님의 <R 통계의 정석>을 사용한다.

 

 

 

지난 글 까지는 일반적으로 2개 이상의 그룹의 결과변수들에 차이가 있는지 확인하는 여러가지 통계 검정법들에 대해 다루었다.

예를 들면 group A에서는 (가) 수술법을, group B에서는 (나) 수술법을 시행한 결과를 비교하는 형태였고, 이때 group A와 B는 관련이 없는 군이다. 

만약 group A와 B가 모종의 연관이 있는 군이라면 어떻게 다루어야 할까?

 

• 예1) 한 사람의 양측 다리를 각기 다른 수술법으로 수술
• 예2) 한 사람에 대해 수술/intervention 전/후를 비교
• 예3) 한 사람에 대해 약물 A 투여 후, wash-out 가긴 거치고 약물 B 투여한 결과를 비교 (cross-over design)

 

장점

• 한 사람에 대해 두번 검사하는 것이라 동일 검정력에 필요한 샘플 수가 반으로 줄어든다.
• 나이, 체중 등 교란변수가 될 수 있는 demographics가 동일하기 때문에 미세한 차이를 비교하기에도 유리하다.

 

단점

• 성숙 : intervention의 효과인가? 시간이 지난 효과인가?
• 학습 : 교수법이 좋아서인가? 시험을 두번 쳐서 인가?
• 대조군 숫자 제한 : 다리는 2개 뿐인데 수술법 3종류는 비교할 수 있는가?

 

"짝을 이룬 자료"에서 유의할 점

헷갈리는 경우들이 있다. 예를 들어 남편과 아내의 월급을 비교한다고 생각해보자.

• 1) 남편 100명 그룹 vs 아내 100명 그룹의 비교 -> 짝을 이룬 자료가 아니다.
• 2) 남편 100명과 그 각각의 아내 100명을 비교 -> 짝을 이룬 자료

 

짝을 이룬 자료

지난 글에서 다룬 내용의 일부를 다시 한번 꺼내보자.

---

2개 그룹의 비교

결과변수가 연속변수, 잔차 정규성, 등분산 -> t-test

분산 다르면 -> Welch's t-test

정규성 없는 연속변수 or 순위변수 -> Mann-Whitney U test

명목변수 -> 카이제곱검정

명목변수인데 코크란법칙 위배 -> Fisher's exact test

 

3개 그룹의 비교

결과변수가 연속변수, 잔차 정규성, 등분산 -> ANOVA

분산 다르면 -> Welch's ANOVA

정규성 없는 연속변수 or 순위변수 -> Kruskal-Wallis H test

명목변수 -> M×N 카이제곱검정

명목변수인데 코크란 위배 -> Fisher

교란변수 통제 -> ANCOVA

 

꼭 해야하는 Post hoc test

---

짝을 이룬 자료도 위와 동일하게, 결과변수의 형태와 두 그룹의 등분산성 등에 따라 검정법을 선택하게 되는데,

대응되는 검정법들은 다음과 같다. 간단히 먼저 살펴보고, 세세하게 들어가보자.

 

2개 그룹의 비교

결과변수가 연속변수, 잔차 정규성 -> t-test : paired t-test

분산 다르면 -> 집단 간 비교가 아니므로 등분산 검정이 필요 없음! 

정규성 없는 연속변수 or 순위변수 -> Mann-Whitney U test : Wilcoxon signed-rank test

명목변수 -> 카이제곱검정 : McNemar test

명목변수인데 코크란법칙 위배 -> Fisher's exact test : McNemar test(exact, binorminal test, 잘 구분안함)

 

3개 그룹의 비교

결과변수가 연속변수, 잔차 정규성, 등분산 -> ANOVA : Repeated measures ANOVA (반복측정 분산분석)

분산 다르면 -> 등분산 검정이 필요 없음! 

정규성 없는 연속변수 or 순위변수 -> Kruskal-Wallis H test : Friedman test

명목변수 -> M×N 카이제곱검정 : Cochran (Q) test

명목변수인데 코크란 위배 -> Fisher

교란변수 통제 -> ANCOVA

 

 

Paired t-test - R 실습

주어진 데이터는 다음과 같다.

동물 실험의 결과로 처치 전/후의 몸무게를 stacked data로 받았다.

 

동시에 위와 같은 data를 wide type data라고 하며, 

짝을 이룬 데이터를 분석할 때에는 이러한 wide type data를 long type으로 바꿔줘야 하며, 그 모습은 다음과 같다.

 

Wide type -> Long type 변환

pD <- read.csv('./pairedData.csv')
library(tidyr)
pD2 <- pivot_longer(data = pD, cols = -ID, names_to = 'GROUP', values_to = 'RESULT')
View(pD2)

tidyr 라이브러리를 이용하여 long type의 pD2를 생성하였다.

정규성 확인

d <- pD2$RESULT[pD2$GROUP=='before'] - pD2$RESULT[pD2$GROUP=='After']
shapiro.test(d)

참고로, pD$2RESULT[pD2$GROUP=='before']의 경우에는 array나 list와 같은 형식으로 저장되는데, 정확히 무슨 type인지는 모르겠다. R 문법 공부할 것이 아니니 넘어간다.

 

p-value는 0.6141로, 귀무가설 채택 -> 정규분포를 한다.

 

 

검정

ptt <- t.test(RESULT~GROUP, data=pD2, paired=TRUE)
ptt

위처럼 test 결과를 변수에도 저장할 수 있다.

t.test 메서드의 paired=TRUE 옵션만 사용하면 간단히 검정을 시행할 수 있으며,

변수를 호출하면 다음과 같은 결과를 볼 수 있다.

 

굉장히 작은 p-value를 얻었다. 

즉, 처치 전/후의 차이가 명확하다는 결과이다. 

추가적인 정보도 제공한다.

차이의 평균값 = 194.49 과

95% 신뢰구간 = 173.4219, 215.5581 또한 제공해줌을 알 수 있다.

 

 

그래프로 나타내보기

PairedData라는 패키지를 이용한다.

library(PairedData)

before <- subset(pD2, GROUP=="before", RESULT, drop=TRUE)
After <- subset(pD2, GROUP=="After", RESULT, drop=TRUE)
pD3 <- paired(before, After)
plot(pD3, type="profile")

위와 같은 plot을 그릴 수 있으며, 기존 plot과 달리 짝지어진 데이터끼리 선으로 연결되어 있는 것을 볼 수 있다.

(이를 profile plot이라고 한다.)

 

 

Wilcoxon signed rank test - R 실습

이번에 사용할 데이터는 동일한 환자에게 두 종류 수면제를 번갈아 복용시켜 수면시간을 측정한 데이터로서, R에 내장된 sleep이라는 데이터를 사용한다. 내장 데이터를 불러올 때에는 data() 함수를 이용한다.

group에 따라 ID가 반복되므로 long type data이다.

 

정규성 검정

약간의 편의 기능을 배워보자

with 함수를 이용하여 변수 이름 앞에 변수의 소속을 매번 써주는 불편함을 줄일 수 있는 매우 유용한 기능이다.

shapiro.test(sleep$extra[sleep$group==2]
             -sleep$extra[sleep$group==1])

with(sleep, 
     shapiro.test(extra[group==2]
                  -extra[group==1]))

위의 귀찮은 코드를 아래와 같이 더 직관적으로 바꿀 수 있다.

p-value = 0.03334로 정규분포를 하지 않음을 확인하였다.

 

검정

wilcox.test 함수를 다음과 같이 이용한다.

with(sleep,
     wilcox.test(extra[group==2]-extra[group==1], 
                 exact=FALSE))

결과는 다음과 같다.

유의한 차이가 남을 알 수 있다.

 

그래프

profile plot을 그려보자.

drug1 <- subset(sleep, group==1, extra, drop=TRUE)
drug2 <- subset(sleep, group==2, extra, drop=TRUE)
sleep2 <- paired(drug1, drug2)
plot(sleep2, type='profile')

 

 

McNemar test

조금 낯선 검사이다.

짝짓지 않은 변수들에서의 카이제곱 검정과 같이, 결과변수가 명목변수일 때 사용하는 방법이며

2×2 테이블에서 역시 많이 사용하지만 그 이상에서도 가능하다.

 

다만 테이블의 뉘앙스가 카이제곱검정과는 다르다.

 

카이제곱검정에서 사용하는 테이블은 다음과 같다.

column과 row를 각각 이루는, 카이제곱검정에서 "관심"을 두는 것들은

원인-결과 이다.

 

McNemar test는 "짝을 이룬 데이터"에서 사용하는 것이기 때문에 테이블의 형태가 다음과 같다.

column과 row는 각각 paired data가 구성한다.

 

데이터를 직접 만들어서 McNemar test를 시행해보자.

이참에 R에서 자주 사용하는 데이터형인 Vector와 DataFrame도 사용해본다.

 

위와 같은 데이터프레임을 만들기 위하여 vector를 생성하자

opA <- c('만족', '만족', '불만족', '불만족')
opB <- c('만족', '불만족', '만족', '불만족')
count <- c(25, 5, 3, 30)
mydata <- data.frame(opA, opB, count)

그러면 mydata는 다음과 같은 형식이 된다.

 

이 data를 이용하여 교차표를 만들어주면

mytable <- xtabs(count~opA+opB,data=mydata)

만족 비율은 opA에서 0.4761 opB에서 0.4444로 A에서 근소하게 높다.

 

McNemar test를 시행해보면

mcnemar.test(mytable, correct = F)

이처럼 통계적으로 유의하지 않음을 알 수 있다.

 

 

Repeated measures ANOVA

앞선 흐름에서와 마찬가지로 이제 t-test가 ANOVA로 개념이 확장된 것과 같은 도구가 필요하다.

예를 들어,

"2개"의 다리에 각각 다른 수술을 하거나

A약을 먹고 시간이 지나서 B약을 먹거나

하는 것이 아닌

 

피부의 "여러" 부위에 다른 약을 발라보거나

수술 후 "3개월, 6개월, 9개월"에 차이가 있는지 확인하는 경우가 그렇다.

 

이것이 바로 repeated measures ANOVA이다.

앞에서 계속 나왔듯, 결과변수가 정규성이 있는 연속변수여야 하며

정규성 없는 연속변수 or 순위변수 -> Kruskal-Wallis H test : Friedman test

명목변수 -> M×N 카이제곱검정 : Cochran (Q) test

 

이런식으로 진행된다.

 

다음 데이터를 보자.

치료 전, 치료 후 1개월, 3개월, 6개월의 결과변수들이다. 연속변수이다.

 

정규성 검정

현재 wide type data이므로 long type으로 바꿔준 후 (tidyr library 사용)

library(tidyr)
sampleLong <- pivot_longer(data=sampleData, cols=-id, values_to = "score")

이에 대해 aov 함수를 이용한다.

out = aov(score~name, data=sampleLong)
shapiro.test(resid(out))

p value 0.954로 정규성을 확인하였다.

 

 

검정

일원 배치 반복 분산분석을 시행하자.

여기서부터는 R 코드가 좀 복잡해지는데, 그냥 통으로 갖다 쓰면 된다.

mymodel <- lm(cbind(sampleData$score0, sampleData$score1, sampleData$score3, sampleData$score6)~1)
Trials <- factor(c("score0", "score1", "score3", "score6"), ordered=F)
install.packages('car')
library(car)
model <- Anova(mymodel, idata=data.frame(Trials), idesign=~Trials, type="III")
summary(model, multivariate=F)

 

결과는 다음과 같다.

 

반복측정 ANOVA에서는 정규성 검증 외에 구형성 가정(sphericity assumption)에 대한 검증이 필요한데, 

Anova() 함수에서 자동으로 Mauchly test를 시행해주며, 위의 p-value가 0.25494이므로 구형성을 가정할 수 있다.

 

따라서 위의

"~~~ ANOVA Assuming sphericity" 항목의 p value를 선택할 수 있다.

 

 

만약 구형성 가정을 만족하지 못한다면, 

Greenhous-Geisser correction 혹은 Huynh-Feldt correction을 시행하면 되는데

이에 대한 p value가 아래에 나와 있는 2개의 p-value 들이다.

 

따라서 위의 경우에는 p-value 1.473 × 10^-14 를 최종적으로 선택할 수 있으며 반복측정에 대한 변화가 통계적으로 의미있다는 것을 알 수 있다.

 

 

사후 검정

어디서 얼마나 차이가 있는지 post hoc test가 필요하며, TukeyHSD test를 시행한다.

out = aov(score~name, data=sampleLong)
TukeyHSD(out)

따라서 시점의 변화마다 통계적으로 유의미한 차이가 있음을 알게 되었다.

 

그래프

이러한 변화를 그림으로 그려서 직관적으로 살펴보자.

gplots라는 패키지를 사용한다.

library(gplots)
means <- c(mean(sampleData$score0), mean(sampleData$score1), mean(sampleData$score3), mean(sampleData$score6))
se <- sd(sampleLong$score)/sqrt(length(sampleLong$score))
plotCI(x=means, uiw=se, type = 'l', ylab="score", xlab="month", main="One way test")

이런 식으로 그리면 된다.

뭐 지금 x 축 범위도 4개월까지밖에 안나와있고, 그리는 방법은 중요치 않고

위처럼 시간에 따른 변화와 그 오차 범위를 standard error를 계산하여 uiw 옵션을 이용하여 그려준다는 것이 중요하겠다.

 

 

Friedman test

세 가지 방법의 반복 결과를 비교하는데 결과변수가 정규분포하지 않는 경우에는 Friedman test를 시행해야 한다.

예제 데이터로는, Myles Hollander와 Douglas A. Wolfe의 1999년 책에 나온 RoundingTimes라는 데이터를 사용하며,

 

야구에서 타자가 홈플레이트에서 출발하여 1루를 거쳐 2루로 가는 3가지 방법에 대해 18명의 선수가 반복하여 얻은 훈련 데이터이다.

 

정규성 검정

wide type data를 long type으로 변환하기 위해 reshape 패키지의 melt()함수를 사용한다.

데이터의 형태가 다르기 때문에 다른 함수를 사용하는 것이며, 취사선택하면 된다.

install.packages("reshape")
library(reshape)
RT2 <- melt(RoundingTimes)
RT2
out <- aov(value~X2, data=RT2)
shapiro.test(resid(out))

p-value 0.001112로 정규성을 가정할 수 없다. 따라서 프리드먼 검정을 시행해야 한다.

 

검정

friedman.test(RoundingTimes) 
boxplot(value~X2, data=RT2)

p-value 0.003805로

세 가지 경우 중 적어도 하나는 결과에서 차이가 남을 알 수 있다.

그림으로 그려보면

위와 같다.

 

Post hoc test

패키지로 배포된 함수를 찾지 못하여, 저자가 찾은 함수를 사용하며

그 출처는 https://www.r-statistics.com이다. 코드까지 적지는 않는다.

 

해당 함수를 시행하면

friedman.test.with.post.hoc(value~X2 | X1, RT2)

위처럼 사후검정이 시행된 것을 알 수 있다.

 

(아직은) 이해할 수 없는 plot도 위와 같이 그려준다.

 

 

 

Repeated measures 2-way ANOVA

마지막으로 다룰 방법은 조금 특이한 setting에서의 검정법이다.

 

수술군 : 3개월 결과 / 6개월 결과 / 9개월 결과

대조군 : 3개월 결과 / 6개월 결과 / 9개월 결과

 

위에서 각 군 내에서는 짝을 이룬(독립이 아닌) 결과들이나, 

군과 군 사이에서는 독립인 구조가 된다. 

 

이정도만 알아두고 넘어가도록 하자.

 

 

 

 

마무리

짝을 이룬, 독립이 아닌 자료들을 다루는 방법에 대해 다뤄보았다.

Cochran(Q) test는 자세히 다루지 않아 넘어간다.

저번 글과 마찬가지로 흐름만이라도 정확히 알고 넘어갈 수 있도록 하자.

더 나아가 R 기초 문법과 자료형에 대해 알아보는 시간도 가져보도록 하자.