탐색적 자료분석 3-2. 두 변수의 관계 파악 - 연속 변수 vs 범주형 변수 (R)

이번 주제는 category 변수 형태의 독립변수와 연속 변수 형태의 결과변수의 관계를 분석하는 것이다. 본 주제로 들어가기 전 우선 t-test와 ANOVA에 대해 정리하고 넘어가 보자.

 

1. 두 군의 평균 비교

 

어떤 두 그룹을 비교한다고 해보자. 예를 들어 남학생과 여학생의 성적을 비교한다고 해보자. 성적을 비교한다는 것은 어떤 의미일까? 한 명 한 명의 성적을 맞대어 비교해야 하나? 통상적으로 통계분석에서 뭔가를 비교 한다는 의미는 "평균"을 비교한다는 의미이다.  평균을 비교하기 위해서는 평균이 어느 집단을 대표할 수 있는 수치가 되어야 하며, 분포가 "정규분포" 하는 것이 그 시작이다. 다시 정리해보면 우리가 사용하는 대부분의 통계기법은 평균의 비교이며, 이게 성립하기 위해서는 정규분포 해야만 한다.

 

자. 그럼 정규분포하는 두 집단이 있는데 한 집단은 평균이 0점이고, 다른 집단은 평균이 10점이라고 하자. 과연 두 집단은 차이가 있을까? 점수 차이가 있기는 한데 뭔가 좀 부족한 것을 느낄 것이다. 10점 차이가 어느 정도의 차이인지 모르기 때문일 것이다. 그 차이를 알려주는 것이 "표준편차"이다. 

library(car)
library(ggplot2)

curve(dnorm(x = x, mean = 0, sd = 2), from = -10, to = 20, xlab = "X", ylab = "Probability")
curve(dnorm(x = x , mean = 10, sd = 2), col="red", add=T)
abline(v = 0, lty=2)
abline(v= 10, lty=2, col="red")

위 그림의 검은 선의 분포는 평균 0, SD=2이며, 붉은 선은 평균 10, SD=2이다. 표준편차에 대한 정보가 있으니 비교하기 좀 더 쉬워졌을 것이다. 자 다시 정리해보자. 두 집단의 비교란 평균의 비교이며, 평균의 비교를 위해서는 정규 분포하고 분산(표준편차의 제곱)의 정보가 있어야 한다. 이를 이해하면 t-test의 이해가 쉬워진다.

df1 <- data.frame(x1=rnorm(n = 50, mean = 0, sd = 2), x2=1)
df2 <- data.frame(x1=rnorm(n = 50, mean = 10, sd = 2), x2=2)
df3 <- rbind(df1, df2)
t.test(x1~x2, data=df3, var.equal=T)

   

두 집단을 비교해 보니 group1(x2=1)의 평균은 -0.31, group2(x2=2)의 평균은 10.30이며, 두 집단의 평균 차이의 95% 신뢰구간은 -11.30~-9.91으로 유의한 차이(p-value <0.001)를 보였다.  t.test 명령어에 var.equal=TRUE 라는 문구를 넣었는데 이건 등분산(두 집단의 분산이 다르지 않다)을 가정하고 분석하라는 의미이다. 이 명령어를 넣지 않으면 R의 t.test 명령어는 등분산을 가정하지 않은 Welch test의 결과를 제시한다.

 분산이 같게 만들어졌기 때문에 등분산을 가정한 결과와 크게 다르지 않다. 그럼 등분산인지 어떻게 확인할까?

boxplot(x1~x2, data=df3)

가장 좋은 방법은 박스플롯을 그려 눈으로 확인하는 것이다. 박스플롯은 간단해 보이는 모양에 비해 많은 정보를 담고 있다. 위와 아래의 수염은 각각 4, 1 quantile, 박스의 위 와 아래는 3, 2 quantile, 박스 안의 굵은 선은 median이다. 굵은 선이 박스 안 가운데에 위치하며 위아래 수염 길이가 비슷하면 얼추 정규분포하는 것 (최소한 skewness는 없는 것)으로 본다. 거기다 양쪽 그룹의 퍼진 정도가 비슷하면 얼추 등분산 한다고 볼 수 있다. 

leveneTest(x1~factor(x2), data=df3)
bartlett.test(x1~x2, data=df3)

   

물론 R의 기본 설정 (stat package)의 bartlett test와 car package의 levene test를 시행하는 것도 좋은 방법이지만 개별적인 검사 수치가 큰 의미를 지닌다고 생각하지는 않는다. 계속 언급하는 내용이지만 논문 통계의 목적은 다변수 분석이며, 다변수 분석을 위해서는 정규성과 등분산이 기본이다. 단변수 분석에서는 정규성과 등분산 만족 못해서 다른 방법 썼는데 그 후에 다변수 분석합니다^^ 도 말이 안 되는 얘기고, 단변수 분석에서 정규성과 등분산 만족 못해서 다변수 분석 안 가고 이만 접습니다^^ 도 똑같은 얘기다. 사전 분석에서 개별 변수의 정규성 확인하고 변환해 놓으면 어지간해서는 등분산도 크게 어그러지지 않는다.  두 집단의 비교는 해석하기도 쉽다. p-value 확인해서 <0.05 면 차이가 나는 것이다. 

 

 

 

2. 세 군 이상의 평균 비교

df1 <- data.frame(x1=rpois(n = 50, lambda = 1), x2=1)
df2 <- data.frame(x1=rpois(n = 50, lambda = 1.5), x2=2)
df3 <- data.frame(x1=rpois(n = 50, lambda = 2), x2=3)
df4 <- data.frame(x1=rpois(n = 10, lambda = 2.5), x2=4)
df <- rbind(df1,df2,df3,df4)

Poisson 분포를 하는 네 변수를 묶어 하나의 자료를 만들어 보았다. 이걸 히스토그램이랑 박스플롯으로 확인하면,

hist(df$x1)
boxplot(x1~x2, data=df)

Right-shifted 되어 있는 전형적인 의학연구에서의 변수이다^^. 왠지 히스토그램만 봐도 로그변환 해야 할 느낌이 오지만, 0 이 많기 때문에 그냥 로그 변환하면 안 되고 Yeo-Johnson 변환에 의해 log(x+1)해줘야 한다 (이해 안가면). 박스플롯을 보면 위아래 수염의 길이가 제각각이다. 등분산이고 뭐고 이대로 분석하면 안 되는데...   

의외로 Levene test는 통과했다^^. 수치만 보고 그대로 믿으면 안되는 이유가 하나 더 늘었다.

로그 씌워서 다시 보면,

boxplot(log(x1+1)~x2, data=df)
leveneTest(log(x1+1)~factor(x2), data=df)
bartlett.test(log(x1+1)~x2, data=df)

 

박스플롯도 좀 더 나아 보이면서, 수치들도 좀 더 좋아졌다. (이건 누누이 말하는 것이지만, 이런 사정을 이렇게 단변수 분석했을 때 처음 알게 되는 일이 있으면, 그건 사전 조사가 망했다는 의미다. 미리 로그 변환해놓고 '이럴 줄 미리 알았어'라는 식의 접근이 필요하다.)

 

ANOVA test는 oneway.test 명령으로 시행한다. aov도 있는데 얘는 등분산이 만족하지 않는 경우는 사용 불가하다 (오히려 무조건 등분산 분석해주니 실용적일수도...)

oneway.test(log(x1+1)~x2, data=df)
oneway.test(log(x1+1)~x2, data=df, var.equal = T)

p-value 만 조금 다르고 결과는 비슷하다.

자 그럼 p-value도 유의하게 나왔고.... 차이가 있네. 하고 끝내기에는 뭔가 좀 이상하다. x2는 4집단으로 되어 있는데 뭐가 어떻게 다르다는 건가? 그걸 확인하기 위해서는 사후 분석을 해야 한다.

pairwise.t.test(log(df$x1+1), df$x2, p.adjust.method = "BH")

Benjamini-Hochberg method에 의한 사후 분석법이다. SPSS의 Hochberg 방법이다. 뒤에 보여 줄 Tukey법에 비해 결과를 눈에 보기 편하다. 1집단은 2,3,4 집단과 유의하게 차이나고, 2집단은 3집단과 차이를 보이며, 3집단과 4집단의 차이는 유의하지 않다. 결과를 보면 4집단의 결과가 이상한데... 그 이유는 n수가 적어서 이다. (다른 집단은 각각 50명인데, 4집단만 10명이다.) 집단의 n수가 적으면 통계적인 유의성이 줄어들기 때문에 이 경우 3과 4집단을 합치는 것이 하나의 고려대상일 수 있다. 아무 집단이든 무조건 합치라는 얘기는 아니다. 합쳐도 무방한 이유가 있는 경우 합치면 좀 더 좋아질 수도 있다는 얘기다. 

m1 <- aov(log(x1+1)~factor(x2), data=df)
TukeyHSD(m1)

oneway.test가 아닌 aov로 분석한 경우 Tukey 법으로 사후 분석 결과를 볼 수도 있다. 방법에 따라 별 차이는 없지만 간혹 차이가 나는 경우도 있는데 이런 경우 맘에 드는 결과를 쓰라는 유혹이 존재하기도 한다. 그런데 사실 ANOVA결과는 table 1이나 2에서 간 보는 정도로나 쓰지 진짜 결과는 어차피 다변수 분석이다. 사후 분석법이 어떤 게 더 정확하니 마니... 그런 논쟁은 큰 의미 없어 보인다. (그리고 변수 이름에 factor()를 넣는 이유는 명령어가 category 변수를 강제하는 경우 때문이다. 사실 이런 경우가 종종 있는데, 잘 안되면서 에러 메시지만 받는 경우 짜증이 날 때가 있다. str() 명령으로 변수 구조 확인해서 바꿔 주면 된다.) 

 

 

3. GLM으로 분석

자꾸 언급하는 내용이지만 t-test나 ANOVA는 본 분석 전에 하는 간 보기 같은 느낌이지, 절대로 논문의 main이 되는 분석은 될 수 없다. 어차피 GLM으로 분석한 결과가 단변수 분석의 결과이니 확인하고 넘어가자.

m2 <- lm(log(x1+1)~factor(x2), data=df)
summary(m2)
anova(m2)

사실 R에서는 선형 모형(LM, Linear Model)을 사용 할 수 없다. 위의 명령에서는 lm이라고 선형 모형을 쓰는 듯 했지만 사실 일반선형모형(GLM, General Linear Model)이다. 둘이 비슷한 거 같은데.. 결론부터 말하자면 일반선형모형 = 선형모형 + ANOVA이다. 원래 선형 모형에서는 카테고리 변수를 처리할 수 없다. 더미변수로 변환해서 넣어야 한다. 그런데 R에서는 기본이 GLM이니 더미변수 변환 없이 그냥 카테고리변수 넣어도 GLM으로 알아서 처리해 준다. 그리고 anova table을 얻을 수 있다. 먼저 summary 결과를 보면 Intercept는 0.58로 x2가 기준(reference, x2=1)일 때  log(x1+1)=0.58 이라는 것이다. 그리고 x2=2일때는 log(x1+1)=(0.58+0.24= 0.82), x2=3일때는 log(x1+1)=(0.58+0.48= 1.06)... 그리고 x2 변수 옆의 p-value는 x2가 기준 (x2=1) 일때와 비교해 유의한 차이를 보이는지이다.  그 아래 anova 결과는 x2가 변수로 투입되어서 투입되기 전에 비해 log(x1+1)에 대한 설명력이 증가하였나이다. 기존 모형은 없으니까 전체 평균만을 고려한 모형이고, x2가 도움이 되었나? 쓸모 있나?라는 것이다. 분석을 하다 보면 anova결과에서는 유의하지 않은데 glm결과에서는 하나 정도가 (예를 들면 x2=3 정도가) reference와 차이를 보이는 경우가 있다. 뭔가 유의한 것 같지만 원칙적으로는 x2 변수 자체가 도움이 안 되기 때문에 쓰면 안 된다.    

그리고 ANOVA와 anova 동일 어구를 대소문자로 쓰는 이유는 anova는 모형끼리의 비교이며 모형 내에 카테고리 변수가 하나도 없더라도 시행할 수 있다.(변수가 모형 설명에 도움이 되는지를 보는 것이기 때문에...) 하지만 ANOVA는 우리가 위에서 계속 쓰던 ANOVA이다. 

결과 제시는,

Variable		Estimates	SE	p-value
x2	1	Reference		<0.001
	2	0.241	0.109
	3	0.484	0.190
	4	0.614	0.190

정도로 한다.  (변수의 개별 항목에 속한 p-value를 제시하는 경우도 있지만 anova 결과의 p-value 제시가 더 합당하다.)   

 

그리고 R에서는 더미 변수의 작성이 그리 자주 필요하지는 않은데, 필요하다면 model.matrix 명령에 의해 만들 수 있다.

m3 <- lm(log(x1+1)~model.matrix(~factor(x2)-1), data=df)
anova(m3)
summary(m3)

(model.matrix 명령어는 변수에 따른 더미 변수를 모두 만들기 때문에 (위 같으면 중복 때문에 x2=4 분석 안 했다.) 하나는 분석이 안된 상태로 있다. 보기 지저분할 수도 있으며 reference를 직접 정하는 차원에서 따로 변수로 지정해서 하나만 빼고 (그게 reference) 투입하는 것이 정석이다. (model.matrix 명령 뒤에 -1은 intercept를 제외하는 명령이다.)

 

그리고 reference를 변경하고 싶으면,

df$x2 <- factor(df$x2)
str(df$x2)
df$x2 <- relevel(df$x2, ref = "2")
str(df$x2)

relevel 명령어에 의해 변환 가능하다 (reference 1 --> 2)

 

 

4. Take Home Message

1) 두 집단의 비교는 t-test, 3 집단 이상의 비교는 ANOVA

2) 정규성과 등분산에 너무 목숨 걸지 말자.

3) ANOVA와 anova는 다르다.