chisq.test, t-test 및 ANOVA 통계의 자동화

[cardiomoon]

올해 들어 처음으로  출판될 예정인 제 논문제목입니다. 출판사에서 eProof를 받아서 교정후 보내서 이제 며칠 안에 출판될 예정인 논문입니다.

이 논문에 사용된 통계는 모두 R을 사용하였읍니다.(안재형선생님께 감사드립니다). 

이 논문의 표를 보면 다음과 같습니다.

즉 242명의 환자 자료를 세 군으로 나누어 범주형 변수는 Chi-square 또는 Fisher’s exact test를 했고 연속형 변수는  ANOVA test 혹은 K-W test를 했읍니다.

제가 돌린 통계를 보면 연속형변수의 경우 ANOVA를 한후(lm함수 이용)  잔차의 정규분포 가정을 만족시키는지 확인 후(shapiro.test이용)   정규분포 하지 않을 경우 K-W test를 했읍니다. 또한 표에는 median 값과 interquartile range  또는 평균±표준편차를 제시해야 하므로 tapply함수를 이용, fivenum과 평균, 표준편차를 구합니다.

물론 이 과정을 한번에 할수 있는 myanova.test 함수를 하나 만들어서 사용했지만(이 함수에 대하여는  R을 이용한 누구나하는 통계분석 게시판 70번을 참조하세요) 그래도 다음과 같은 과정을 거쳤읍니다.

myanova.test(N2BMin,Group1,pci)

myanova.test(Age,Group1,pci,1)

myanova.test(D2Bmin,Group1,pci)

myanova.test(Fluoro,Group1,pci,1)

myanova.test(Contrast,Group1,pci,1)

myanova.test(CPKpeak,Group1,pci,1)

myanova.test(CKMBPeak,Group1,pci,1)

myanova.test(EF,Group1,pci,1)

myanova.test(T1min,Group1,pci)

myanova.test(T2min,Group1,pci)

myanova.test(N2BMin-T1min-T2min, Group1,pci)

이걸 표로 만들려면 median, IQR, 평균, 표준편차를 노트에 따로 손으로 기록하고 p 값도 일일이 노트에 기록한후 다시 워드프로세서 프로그램으로 표를 기록했읍니다.

통계 돌릴 변수가 수십개라면 그리고 나누어볼 그룹이 수십개라면 이러한 작업은 아주 고된 반복작업이고 옮겨적다가 실수할 수도 있읍니다.

작년까지는 이러한 작업을 모두 일일이  R에서 수작업으로 해야 했읍니다. 수개월동안 반복하다가 통계의 자세한 과정은 생략하고 내가 원하는 결과만 나오게 하는 함수를 하나 만들자는 생각을 했읍니다. 먼저 t-test에 도전 했읍니다. 여기에 연속형 변수의 이름을 나열해주고 grouping variable을 나열해주면 내가 원하는 표에 쓸 수 있는 자료가 나오게 하고 싶었읍니다. myt.test를 개조해서 사용하였는데 다음과 같이 사용하였읍니다.(myt.test의 내용은 게시판 69번을 참조하십시요)

> dat2=read.csv("INFMI.csv")

> x1=c("Age","N2BMin","D2Bmin","FluoroTime","BMI")

> y1=c("Complete","death","MI","MACE","Entry","smoking")

> 

> myt.test(x1,y1,dat2,drawgraph=0,mode=0)

mode가 0이면 median과 IQR이 출력되게 했고 1이면 평균과 표준편차가 출력되게 했읍니다. drawgraph가 1이면 boxplot을 그렸읍니다.

> myt.test(x1,y1,dat2,drawgraph=0,mode=1)

그러다보니 grouping variable이 두그룹이면 t-test, 세그룹 이상이면 ANOVA 를 돌리고   K-W test까지 하는 함수로 extension 했읍니다.

또한 chisq.test도 이와 비슷하게 개량해서 사용하였읍니다.

> x=c("Sex","MACE","DM","HBP","lateral.wall.extension","precordial.ST.depression")

> y=c("Complete","death","MI","MACE","Entry","smoking")

> 

> mychisq.test(x,y,dat2)

마지막으로 그냥 범주형자료던지 연속형 자료던지 그냥 변수 이름을 넣어주면 알아서 범주형자료면 chisq.test를 하고 연속형자료면 t-test 또는 anova.test를 하는 함수로 개량하였읍니다. 드디어 my.test라는 함수를 만들었읍니다. 함수를 만드느라 일주일 이상 고생했지만 써보니 무지 편합니다. 

이 맛에 프로그래밍 하는 거죠. 이제 내가 보고싶은 변수를 그냥 다 넣어주면 알아서 통계를 돌려주고 표까지 그려줍니다.

> x2=c("Sex","Age","N2BMin","DM","D2Bmin","FluoroTime","BMI","HBP")

> my.test(x2,y1,dat2,mode=1) 

이 프로그램이 개인적으로 아주 맘에 듭니다. 그래서 package로 만들어 배포해볼까 생각중입니다. 원하신다면 소스를 모두 공개하겠읍니다.

[안재형]

앗, 제책에 Kruskal-Wallis를 Kurskal-Wallis라고 잘못 쓴 오타가 있었는데 여기도 그렇게 하셨네요. 3쇄에서 수정이 되었습니다. 고치셔야겠는데요^^ 죄송^^

[cardiomoon]

음,,,할수 없죠...뭐 이미 교정도 다 봐서 출판사로 넘겨서요...그냥 spelling 틀린 건 남들이 신경 안쓰기 때문에 패스할랍니다.

[강성찬]

대단하시군요. 진짜 패키지를 만들고 소개하시면 다들 쓰실 것 같네요.

[cardiomoon]

통게 전문으로 하시는 분들은 별로 맘에 안들어 하실것 같아요.. 논문 쓸때 꼭 필요한 부분만 출력해주니까요. 저처럼 통계에 문외한인 사람이 단순 반복작업을 줄이는데 필요한 겁니다.

[안재형]

패키지 만드셔도 되겠네요. LaTex 쓰는 사람은 이걸 LaTex 테이블로 출력하면 되겠네요.

[cardiomoon]

참, 그리고 제가 SPSS 쓰던 버릇 있어서 범주형자료인 경우 Yes/ No 가 아니라 0,1로 코딩하는 습관이 있읍니다.

그래서 그냥 자료를 보면 이게 0,1의 값을 갖는 연속형변수인지, 범주형자료인지 분간할수가 없읍니다. R에서는 그냥 numeric으로 인식해서요.. 그래서 제가 임의로 변수값이 5개 이내이면 그냥 범주형자료로 인식하게 했읍니다. 혹 좋은 아이디어 있으신가요?

[cardiomoon]

그리고 소스를 보고 싶어히시는 분이 아무도 없네요...흑흑...역시 통계 전문가들에게는 별 의미없는 일이겠지요

[강성찬]

의미없다니요. 사실 R의 출력결과를 문서화하기 좋은 형태로 하는 것도 굉장히 의미있는 작업입니다. 사실 코드를 올려주셔도 전 잘 이해를 못합니다.^^

[남현수]

제가 현제 만들고 있는 Evaluation 프로그램에 참고로 쓸수있도록 소스코드공개좀 부탁드립니다. 저와비슷한 분야를 패키지로 만들고 계시다니 반갑구 참고하고 싶은 욕심이 불끈합니다. 라이브러리만 가져다 C++로 짜시는건지 아님 패키지로 R스크립트로 짜시는건지 아님 저처럼 델파이나 일반프로그램으로 CMD 라인 커트롤 하시는건지도 좀 알고 싶습니다.

[변효진]

전 cardiomoon님께서 만드신 것 너무 이용하고 싶지만,,이정도도 어려워 보입니다..에효..^^ 아, 그리고 내용과 관련이 없는 줄 알지만, 제가 요즘 들은 수업에서 ANOVA에 이미 analysis가 들어가 있으므로 논문쓸 때 test를 붙이지 말라고 하시던데요. 제가 들을 것이 꼭 그런 것은 아닌가요?

[cardiomoon]

글쎄요 저는 그 문제는 잘 모르겠읍니다. 주인장님께서 답변해 주시겠죠...

[우창윤]

소스 코드 부탁 드려도 되겠는지요.. 보고 열심히 연구하겠습니다^^

[박상일]

어렵네요. 스크립트를 복사해서 그대로 실행을 해도 "Error in myt.test(x1, y1, neph, drawgraph = 0, mode = 0, 1) :
unused arguments (mode = 0, 1)" 이런 말이 나오네요.

[박상일]

우선 69번 게시글에 있는 스크립트를 그대로 옮긴다음 실행을 했고, 변수를 맞게 설정한 것 같은데 안되네요.
만들어진 함수는 스크립트로 실행하는 것이 아니라 다른 방법으로 해야 하나요?