R을 이용해서 검정력분석(power analysis)을 해보자. -깜신의 통계 이야기-

들어가는 글

---

검색을 통해 이 글을 읽으시는 분들이라면 검정력 분석(power analysis)이 무엇이고, 왜 필요한지에 대해 이미 알고 계시리라 생각됩니다. 간단히 요약하자면, 실험의 결과에 대한 P 값이 0.05보다 크게 나와서 통계적으로 의미가 없어 보일 때, 이에 정말 의미가 없는(실험군과 대조군 간에 차이가 없는)건지 아니면, 대상자(n)수가 작아서 의미가 없어 보이는 건지 확인하기 위한 절차이지요. 그래서 열등성검정(실험군이 대조군에 비해 열등하지 않다는 가설을 세우고 진행하는 실험)에서는 특히나 꼭 필요한 것이 바로 검정력 분석입니다. 하지만 검정력 분석을 SAS나 SPSS 등 다수 유료프로그램에서 조차 지원하고 있지 않아서 임상연구 초보자 분들이 애를 먹지요. 

 그래서 이번 글에서는 통계프로그램 ‘R’을 활용해서 검정력을 분석하는 과정을 정리해보았습니다. R은 통계전문가 그룹이 사용하는 고급소프트웨어이지만, GUI(윈도우 시스템처럼 그래픽 아이콘을 마우스로 클릭해서 프로그램을 실행하는 인터페이스)가 아닌 coding(예전 MS-DOS 시절처럼 명령어를 직접 입력해서 프로그램을 실행하는) 방식을 취하고 있어서 어려워 보이는 게 사실입니다. 하지만 조금만 배워보면 MS-DOS 시절 배치파일을 만들어가며 ‘페르시아의 왕자’나 ‘바바리안’을 실행했던 것처럼, 여러분도 t-test나 Wilcoxon Signed-Rank Test를 하실 수 있을 거라 확신합니다. 그럼, 밝고 푸른 R의 세계로 들어가 보도록 하죠.

묻지 말고 따라 하기! R을 이용한 검정력 분석

---

검정력 분석을 위한 준비물

1) 일단, R이 필요하겠죠. R은 앞서 설명해 드린 대로 무료소프트웨어라서 인터넷에서 자신의 컴퓨터 환경에 맞는 버전을 다운받아 설치하시면 됩니다. 

http://www.r-project.org 로 접속하시고요.

왼쪽 메뉴바의 중단에 있는 CRAN을 클릭하세요.

Mirrors 중에서 Korea를 찾으시면 됩니다. (우리나라에는 두 곳이 있는데요. 둘 중 마음 가는 곳에서 다운받으세요.)

R은 총 세 가지 버전으로 개발되어 있습니다. 보시는 것처럼, 리눅스와 맥 그리고 윈도우 버전이 있으니까 자신의 환경에 맞는 프로그램을 받아 설치하도록 하세요. (참고로 저는 mac 환경에서 작업합니다.)

2) 그리고, 연습해볼 데이터가 필요하겠죠. 저도 첫 임상연구를 시작할 때 통계를 잘 아는 친구에게 n 수는 몇 명 정도가 적당할지 물어본 적이 있었습니다. 지금 생각해보면 참 부족한 질문이죠. 왜냐하면, 검정력 분석을 통해 적정 n 수를 계산하기 위해서는 선행연구가 필수이기 때문입니다. 선행연구를 통해서 실험군과 대조군의 실험값에 대한 평균과 표준편차가 구해져야 검정력 분석을 할 수 있거든요. 이번 글에서는 제가 선행연구 결과에 해당하는 데이터를 미리 드리도록 하겠습니다. 

R에서도 엑셀과 비슷한 화면을 제공하고 그 창을 통해 데이터를 직접 입력할 수도 있지만, 데이터 입력은 아무래도 평소 익숙한 엑셀이 더 편하더군요. 그래서 저는 행과 열을 만들고 데이터를 입력하는 과정까지는 엑셀이나 Number(맥에서 사용하는 엑셀과 같은 스프레드시트 프로그램입니다.)를 이용합니다. 그리고 확장자 .CSV 파일로 저장한 후 R에서 데이터를 불러들입니다. 

시간이 조금 여유로우시다면, 이참에 엑셀에서 위의 데이터들을 수작업으로 입력한 후 blog.csv 파일을 직접 만들어 보시는 것도 좋을 것 같습니다. 하지만 바쁜 현대인이시라면 링크를 걸어놓은 CSV 파일을 바탕화면에 저장해놓고, 다음 단계로 저장하셔도 무방합니다.

blog.csv

 

R의 워크스페이스(workspace:작업공간)로 데이터 불러오기

---

자, 이제 데이터를 R로 불러들일 차례입니다. 

R을 실행해보세요. (처음이라면 무척 당황스러우실 겁니다. ^^; 마음 단단히 먹고 클릭하세요.)

깜빡이는 프롬프트가 보이시나요?

이제 R의 명령어를 하나하나 배워보도록 하겠습니다.

첫 번째 명령어는 csv 파일을 불러오는 명령어 read.csv()입니다.

먼저, 앞으로 쓸 ‘mydata’라는 오브젝트에 우리가 미리 만들어놓은 blog.csv 데이터를 옮겨볼게요. 

프롬프트에 다음처럼 입력해보세요.

핑크 사각형 안의 공간은 예상하시는 것처럼, 여러분의 blog.csv 파일이 위치해있는 경로로 수정해서 입력해주셔야 합니다.

(주의: 윈도우 사용자라면, 작은따옴표 안에 ‘C:/blog.csv’라고 경로를 지정해주시면 됩니다. 이때 주의하셔야 할 게, ‘C:\’가 아니라, ‘C:/‘로 쓰셔야 한다는 거죠. R은 기본적으로 Unix 체제를 따르기 때문에 그렇습니다.)

작업 폴더를 미리 지정해두는 방법도 있는데, 그건 다음에 알아보도록 하죠.

자, 해보셨나요?

에러 메시지가 없이, 그냥 또 프롬프트만 깜빡이고 있다면, 성공하신 겁니다. 뭔가에 성공하면, 확인 버튼을 눌러달라고 꼬리를 흔드는 윈도우하고는 달라도 너무 다르죠, R은 잘 된 일에 대해서는 무심한 듯 한편으론 시크하게 그냥 넘어갑니다. ^^; 익숙해지면 이게 또 매력이에요. 그렇다면, 이제 ‘mydata’라고 입력해보세요.

이런 결과가 나왔다면, 여러분이 방금 mydata 라는 오브젝트를 생성하고 그 안에 데이터를 잘 안착시킨 겁니다.

그다음 배울 함수는 attach()입니다. 이 함수는 오브젝트를 메모리에 상주시키는 일을 하는데요. ‘이 오브젝트는 내가 계속 쓸 거니까, 네가 잘 기억하고 있어.’라고 R에게 명령을 내리는 함수입니다. 이렇게 attach 시켜놓으면 그때부터는 detach() 함수를 써서, 잊으라고 할 때까지 R은 우리가 지정한 오브젝트를 늘 기억하고 있죠. 그래서 오브젝트 안의 변수들을 지정할 때 매번 ‘mydata$age’나 ‘mydata$group’이라 길게 표현하지 않고, 그냥 ‘age’나 ‘group’이라고만 하면 됩니다. R은 이 변수들을 attach(mydata) 시켜놓은 오브젝트 안에서 찾아 불러오거든요. 워낙 많이 쓰는 ‘윈도우 탐색기’ 같은 녀석이니까, 꼭 기억해두시길 바랄게요.

attach(mydata)라고 입력하면 또, 마찬가지로 시크하게 웃고 넘깁니다. 

본격적인 검정력 분석(POWER ANALYSIS)

---

자, 이제 준비운동은 끝났습니다. 본격적으로 검정력 분석을 시작해보죠. R이 검정력을 분석하기 위해서는 추가로 프로그램 패키지를 설치해주어야 하는데요.

윈도우처럼, 패키지 설치가 복잡하지도 않습니다. 그냥 프롬프트에 ‘install.packages(“pwr”)'이라고 입력하세요.

처음에는 실수를 많이 하는데요. 스펠링이 package 아닙니다. packages 고요. 함수 안의 패키지 이름은 큰따옴표로 포장해주셔야 에러가 나지 않습니다. 제대로 명령어를 입력하면 다음과 같습니다.

그런데 이렇게 인스톨만 해주어서는 해당 패키지를 바로 사용할 수 없습니다. 프롬프트 창에 library(pwr)이라고 입력해주세요. 직전에 다운받은 pwr 패키지를 R의 라이브러리에 장착하는 과정입니다.

Two sample t test에서 검정력을 분석해보자.

---

이번 글에서는 Two sample t test에서 적정 n 수를 계산하는 방법을 알아보겠습니다. 임상 연구에서 가장 많이 활용하는 통계가 동질 두 집단 간의 비교니까요. 이 과정을 제대로 이해하고 나면, ANOVA나 비율 검정에서도 pwr 패키지를 이용하시는 데 큰 어려움은 없을 겁니다.

Two sample t test에서 pwr 패키지를 이용해 검정력을 분석하는 공식은 다음과 같습니다.

pwr.t.test(d=‘effect size’,power=.8,sig.level=.05,type="two.sample",alternative="two.sided")

임상연구에서 가장 많이 쓰이는 power level 0.8, 유의수준 0.05, 양측 검정을 기준으로 한 공식입니다. 문제는 effect size는 따로 계산해주어야 pwr.t.test()함수를 이용할 수 있다는 거죠.

Effect size를 구하는 Effective한 방법

---

저도 처음 검정력을 계산할 때 effect size를 구하는 과정에서 가장 애를 먹었는데요. Effect size를 구하는 공식은 연구 설계에 따라 여러 방법이 있습니다만, Two sample t test에서 가장 보편적으로 사용되는 공식은 Cohen’ d 입니다. 그래서 pwr 패키지에서도 effect size의 약자를 ‘d’로 표기하지요.

Effect size, d는 다음 공식으로 구합니다.

이때, Mean1은 실험군의 평균, Mean2는 대조군의 평균을 의미합니다. 그렇다면, 또, pooled standard deviation(Sp)은 어떻게 구할까요? 공식은 다음과 같습니다.

선행연구(Pilot study)에서 실험군과 대조군의 n 수가 같은 경우:

(S1은 실험군의 표준편차, S2는 대조군의 표준편차)

하지만, 선행연구에서 n 수가 똑같이 딱 떨어지는 경우는 많지 않죠. n 수가 다르다면, n 수에 대한 가중치를 계산해주어야 합니다. 다음처럼요.

(n1은 실험군의 n 수, n2는 대조군의 n수)

그럼, 이제 제가 나눠드린 예제를 가지고 effect size를 계산해보죠. 위의 공식들을 꼼꼼히 살펴보시면, effect size(d)를 구하기 위해서는 두 그룹의 평균과 각각의 표준편차 그리고 각각의 n 수가 필요함을 알 수 있습니다.

엑셀이나 Numbers에서도 평균과 표준편차, n 수 등은 쉽게 계산할 수 있으니, 편하신 대로 구하시면 됩니다. 물론, R에서도 mean()함수와 sd()함수 등을 이용해서 각각의 값들을 구할 수 있습니다.

실험군과 대조군의 n 수, 평균값, 표준편차가 다음과 같이 구해지셨나요?!

실험군(group=1)

n1=20

Mean1=14.5

S1=18.69844

대조군(group=2)

n2=21

Mean2=35.90476

S2=26.54224

이제 이 값들을 공식에 넣어서 d 값을 구해보죠.

이제, pwr.t.test()함수 공식에 사용할 변수들이 완성되었네요.

프롬프트 창에 

pwr.t.test(d=0.9283539,power=.8,sig.level=.05,type="two.sample",alternative="two.sided")

라고 띄어 쓰거나 줄 바꾸지 말고, 입력해보세요.

이런 결과가 나오셨나요?!

축하합니다. 지금까지 잘 따라오셨네요. 핑크 사각형 안에 들어 있는 19.2XXXX가 바로 적정 n 수입니다. 물론, 이 수는 실험군과 대조군 각각에서 필요한 n 수고 사람은 소수점이 있을 수 없으니, 각각 20명씩 총 40명의 대상자가 있어야 한다는 결론을 얻으신 겁니다.

처음 R을 접하신 분들도 R을 가지고 검정력을 계산할 수 있도록 쉽게 쓴다고 썼는데, 어땠는지 모르겠네요. 작은 도움이라도 되셨길 바라면서 이 글을 마칩니다.

씨앗글

en.m.wikipedia.org—Statistical_power

cran.r-project.org—pwr.pdf

www.ats.ucla.edu—t_test_power2.htm

en.wikipedia.org—Effect_size

en.wikiversity.org—Cohen's_d

www.sfu.ca—Wk07_1_Full.pdf