다변수분석법 (R)

normal.Rdata

0.01MB

탐색적 자료 분석이 끝난 후 본 분석인 다변수 분석을 시행해야 한다. 탐색적 자료 분석을 성실히 시행했으면 다변수 분석 결과는 이미 머릿속에 대충 그려질 것이다. 그 내용을 논문 양식에 맞춰서 제시하는 것이 이번 포스트의 중심 내용이다. 

 

 

1. 다변수 분석의 필요성

 

 사실 예전에 개인용 컴퓨터가 일반화되지 않았을 때는 통계작업을 전부 수기로 진행했었다고 한다. 그러니 다변수 분석 같은 것은 상상도 하기 힘들고 전문 수학자나 통계학자 같은 경우만 선형 회귀, 로지스틱 회귀 같은 분석을 시행하고, 나머지 대부분의 학자들은 카이 스퀘어 검정 하나에만 목매달고 있었다고 한다. 세월이 지나서 더 이상 수기로 통계 계산을 할 필요가 없어지고 개인용 컴퓨터가 보급되기 시작하면서 학문 영역에서 쓰이는 통계기법도 점점 발달하게 되었다. SPSS(Statistical Software for Social Science)가 1970년 생이며, SAS(Statistical Analysis System)는 1976년생,  R이 1995년생 (stable beta는 2000년)이니 그 전 사람들은 어쩔 수 없었을 것이다.

 근데 일이 쉬워진다고 꼭 좋은것만은 아니다. 이 바닥이 어차피 경쟁사회라서 남들이 지겹게 쓴 내용 그대로 써서는 논문에 실리기 힘들다. 뭔가 새로운 것이 있어야 하니 자꾸 이전에 안 쓰던 기법을 이용할 수밖에 없다. 그러다 보니 특히 NEJM 같은 경우에는 정상적인 교육을 받은 전문의들이 도저히 이해할 수 없는 기법을 이용한 연구가 대부분이다. 심지어는 RCT에도 뭔가 새로운 기법들이 자꾸 도입되어서 보는 사람의 눈을 어지럽게 만든다.

 원래 전문의 수련 과정에 논문 작성이 포함되어 있는 이유는, 최신 의학지식을 받아들여서 그 결과를 임상현장에 적용하라는 의미인데, 최근의 의학연구는 너무 복잡한 통계기법이 사용되기 때문에 정상적인 교육과정 과정 (선형회귀, 로지스틱 회귀, 생존분석)을 거친 의사들은 최신 의학연구를 이해할 수도, 연구에 대한 아이디어가 있다고 쳐도 본인이 스스로 연구를 진행해서 해외 학술지에 논문을 실을 수 없게 되었다.  이와 같은 괴리를 벗어나는 방법은 최신 지견이 이해될 때까지 통계적 지식을 늘리거나, 아니면 최신 연구결과를 참조하는 것을 완전히 포기하고 연수교육이나 교과서 등의 방법을 이용해 지식을 습득하는 방법 (제대로 교육받은 의사 = 의학자라는 대전제를 포기) 밖에는 남지 않은 것 같다. 

 쓸데없는 쪽으로 말이 길어졌는데, 제 길로 다시 돌아와 말하자면, 지금 세상은 개인용 컴퓨터가 너무도 흔한 세상이니 최소한 다변수 분석은 필수라는 것이다. 그게 싫으면 1970년대 이전으로 돌아가는 수밖에 없다.

 

 첨부된 파일 내의 다이아몬드 데이터를 꺼내보자.

library(car)
library(psych)
library(dplyr)
library(ggplot2)

ggplot2에 내장된 diamonds파일의 축약본이다.

다이아몬드 가격을 결정짓는 요소로 4C라는게 존재한다고 한다 (Carat, Cut, Color, Clarity). 결혼반지 살때는 이런 게 있는 줄도 몰랐고, 크기(carat)만 생각했었는데...

cut에 따른 다이아의 가격을 본다고 해보자.

with(diamond,boxplot(price~cut))
with(diamond, leveneTest(price~cut))
with(diamond, pairwise.t.test(price,cut, p.adjust.method = "BH"))

  

등분산 검정은 르벤 검정에서 등분산 만족 못하는 것(p=0.049)으로 나왔다. 그래서 oneway.test로 결과를 보니 최소한 한 그룹은 다른 그룹들과 차이가 있었다 (p=0.004). 그래서 pairwise comparisone 결과를 보니 Good vs  Premium, Premium vs Ideal 이 유의한 차이를 보였다. 

	subgroup
Fair	A
Good	A
Very Good	A
Premium	A, B
Ideal	A, C

어? Ideal 보다 Premium cut이 비싸네? 당연히 Ideal cut이 비싼줄 알았는데... 세계에서 내가 첫 번째로 알아낸 사실이다^^. 세계 최초 너무 좋아하면 안된다. 세계 최초라는 것은 두 가지 의미가 있다. 내가 천재 거나 또라이거나. 근데 내가 천재일 가능성은 그리 많지 않으니 당연히... 세상은 절대로 만만치 않다. 그리고 이런 결과는 밖에 유출되면 절대로 안된다. 

 

 얼마전부터 인터넷 뉴스들을 보면 건강 부분에 의학연구결과들이 심심치 않게 실리는 것을 볼 수 있다. 알고 보니 연구결과가 언론을 타면 지도교수(지도전문의)에게 인센티브를 주는 곳도 있단다. 그러니 아무리 허접한 논문이라도 어디든 실리면 무조건 보도자료 뿌린단다. 결과가 기존 상식과 큰 차이 없는 결과면 다행인데, 문제는 기존 상식과 다른 경우 발생한다. 위 결과와 같이... 기존 상식과 다르니 기자들은 얼씨구나 하고 보도자료대로 뉴스 낼 거고... 다치는 사람은 결국 주저 자이다.  

 

 다변수 분석하면 괜찮다고? 정답이다. 근데... 다변수 분석을 잘 해야 한다. 제대로 안 하면 제대로 보정되지 않아 단변수결과와 다름없이 나온다. 그럼 위의 결과와 비슷해진다. 

 

 연구에 대해 어설프게 아는 사람의 최고 문제점 중 하나가, 실제(real world)의 법칙이 이론과 어긋나는 점을 보면 실제 세상의 문제점을 찾으려 드는 점이다. 본인의 이론이 틀렸다는 생각은 잘하지 않는다. 겸손해야 한다. 내 지식과 실제 세상의 모습이 다르면 내게 무슨 문제가 있는지 찾고 또 찾아봐야 한다. 개인적으로 학자 인생 중 하나라도 뭔가를 최초로 찾아냈다면 그 사람의 학자로서의 인생은 그리 나쁘지 않았다 생각한다. 그런 일생에 한두 번 있는 이벤트이니 만큼 기존 상식과 어긋나는 결과는 고찰하고 또 고찰해라.   

diamond %>% ggplot(aes(cut, log(price)))+geom_boxplot()+facet_wrap(~cut_number(carat,10), scales = "free_y")

다이아몬드 가격을 결정하는 가장 중요한 요소는 무게(carat)이다. 이런 경우, 특히 관찰연구에서, 중요 요소를 먼저 보정하지 않으면 나머지 결과는 당연히 엉망으로 나온다. 이게 다변수 분석을 해야 하는 이유이다. 

 

 

2. 다변수분석의 실제

 

 앞의 포스트에서 carat과 price는 로그 씌워 분석해야 한다는 얘기 까진 했는데... 범주형 변수를 다루는 포스트에서는 다른 데이터를 사용한 이유로 cut, color, clarity에 대한 탐색적 자료 분석에 대한 결과는... 

 다변수 분석에 들어갈 자료는 단변수 분석에서 유의했던 변수 + 이전 연구들에서 유의했던 변수 + 내가 중요하다고 생각하는 변수 이다. 대개 method 마지막에 statistical analysis 부분이 있는데 거기 다변수 분석에 투입한 변수를 결정한 기준을 넣으면 된다. 단변수 분석기준 p<0.1 보다 낮은 경우 투입한다는 사람도 있는데... 그건 연구자의 맘인 것 같다. 중요한 변수를 다변수 분석에 넣지 않아서 문제가 되는 경우는 흔하지만 중요하지 않은 변수를 넣어서 문제가 되는 경우는 드물다. 그리고 상호작용의 경우 유의한 상호작용만 본 분석에 투입하고 그 내용을 method 부분에 언급하면 문제 없을 것 같다.

 하여간 앞의 포스트에서 설명한 방법대로 분석했고 문제가 없었다는 가정하에 다변수 분석을 해보자.

m1 <- lm(log(price)~log(carat)*(cut+clarity)+color*clarity, data=diamond)
anova(m1)
summary(m1)

다변수 분석의 anova table이다. 아 근데 혹시 m1 회귀식 이해 못하려나? 그리 어렵지 않다. 일반 괄호, 더하기, 곱하기와 완전 동일하다. 다만 겹치는 항목이 있을 경우 하나만 적용된다. 그거 이해하면 그리 어렵지 않을 것이다. 그리고 범주형변수*볌주형변수의 상호작용은 보지 않는다 했는데... 예들이 자꾸 설명한 내용과 다른 것 같다 TT. 하여간 유의하긴 하다. 결과를 제시하는 게 만만치 않아 그렇지.

 결과가 컴퓨터 화면 하나를 넘는다TT. color와 clarity가 여러 개의 범주를 가진 변수이다 보니 상호작용의 개수도 늘어나서 그런다. 만약 실제 연구였으면 나 같으면 color*clarity 상호작용은 제외하고 분석하고 발표했을 거다 (결과가 달라지지 않는다는 전제 하에). 어쨌든 그 결정은 연구자가 하는 거다. 

 

이제 잔차검정 시작해야 한다. 사실 단변수 분석 및 상호작용 제대로 봤으면 잔차분석 그리 필요 없을 수도 있는데... 그렇게 방심하고 있으면 꼭 놓치는 게 있다. 확인해 주자.

residualPlots(m1)

  

    car package의 residualPlots 명령어로 잔차검정을 할 수 있다. 위의 숫자 결과는 2차항 필요성에 대한 Tukey 검정 결과이다. 단변수 분석에서는 분명 직선 관계로 보았는데... log(carat)의 Tukey 검정 결과 (p<0.001) 및 잔차도표 상 웃는 모양의 파란 선... 둘 다 2차항의 필요성을 보여주고 있다. 여기서 2차항을 투입해도 되지만 난 안 할 거다. 선형 모형으로도 충분히 잘 설명되고 있다 생각하기 때문이다. 

m2 <- lm(formula = log(price) ~ log(carat^2)+log(carat) * (cut + clarity) + color *clarity, data = diamond)
anova(m1,m2)

    

하지만 2차항 투입한 모형이 통계적으로 더 우월한 것은 사실이다. 

Breush-Pagan test(잔차 등분산 검사) 하고 싶으면 lmtest package의 bptest 명령어 사용하면 되고, 잔차 정규성 보고 싶으면 shapiro.test(resid(m1)) 쓰면 된다.

 

outlierTest(m1)
influencePlot(m1)

 

이상치 검정의 경우 이전 포스트에서 의학연구에서는 이상치라고 막 제거하는 등의 폭거를 시행하면 안 된다 했다고 했었다.  우리가 살펴봐야 할 부분은 이상치로 언급되는 관측치 (340번) 및 회귀식에 많은 영향을 끼치는 관측치 (다른 애들과 엇나가는 애들, 22, 56,84, 98, 156, 340번)을 확인하면 된다. 이상치에 대해 잘 봐야 하는 이유는 이상치가 있으면 그냥 자기 혼자 이상한 걸로 끝나는 게 아니라 최종 결과를 꽤 많이 건드린다. 

diamond[c(22,56,84,98,156,340),]

22번은 0.3캐럿짜리인데 1069$이다. 사실 다이아 가격은 잘 모르지만 3부 짜리 다이아를 1000만원 넘게 주고 사라면... 그 외에도 언급된 애들 보면 가격이 좀 이상한 것 같기는 하다. 그럼 여기서 할 내용은 이 데이터들의 완결성을 다시 확인하는 거다. 누가 입력 실수한 건지... 원 자료를 찾아서 다시 확인해야 한다. 계속 언급했지만 이상치라고 그냥 삭제하는 폭거는 저지르지 말라고 했다. 비록 이번 데이터는 다이아몬드 가격이라서 문제없지만 자료가 사람인 경우는, 이상한 결과가 나온 이유를 찾아서 그 부분을 수정해야 한다.  

m3 <- update(m1, subset=-c(22,56,84,98,156,340))
compareCoefs(m1,m3)

기울기 차이가 꽤 난다 (1.22 vs 2.43). 이상치에 대해 확인해 볼 필요는 충분하다.

 

df <- expand.grid(carat=seq(from=0.1, to = 2, length.out = 100), cut=factor(unique(diamond$cut)), color=factor(unique(diamond$color)), clarity=factor(unique(diamond$clarity)))
df$pred <- exp(predict(m1, newdata = df))
df %>% ggplot(aes(carat, group=cut, col=cut))+stat_summary(aes(y=pred),fun.y = mean, geom="line")+theme_classic(base_size = 12)+theme(legend.position = c(0.2,0.8))

 여태까지 했던 결과의 재현이다. log를 씌워 분석했기 때문에 결과는 exponential growth model이 되었다.  그리고 마지막에 사족을 달자면... 상호작용 항을 없애도 충분히 괜찮은 모형이다. 결과를 자세히 보면 알겠지만 특정 등급이 최하인 경우 (cut=Fair, color=J, clarity=SI1 )인 경우 carat의 증가에도 불구하고 가격을 제대로 못받는 현상이 보인다(위의 Fair도 혼자만 따로 논다). 그런 최하등급품들의 문제가 회귀식을 더 복잡하게 만든 이유다.   

 

 지금까지 다변수 선형회귀 분석을 시행하였다. 탐색적 자료 분석부터 시작하면 꽤 많은 과정이 필요함을 알 수 있다. 거기에다 실제 연구에서는 통상적으로 이번 예제에 사용된 변수들보다 많은 변수를 다루고 있어 걸리는 시간은 더욱 늘어날 것으로 생각된다. 하지만 어차피 반복 과정이며 자주 하다 보면 점점 빨라진다.

 

 

 

Take Home Message

 

1. 제대로 된 다변수 분석이 중요하다.

2. 잔차검정, 이상치 분석

3. 연구를 하다 보면 어느 정도 적당히 합의할 필요도 있다.