Journal of Health Informatics and Statistics Brief Communication J Health Info Stat 2016;41(1): p

Transcription

1 Brief Communication J Health Info Stat 2016;41(1): pissn eissn 중재메타분석을위한 STATA 명령어 심성률 1, 신인수 2, 배종면 3 1 순천향대학교서울병원임상분자생물학연구소, 2 전주대학교사범대학교육학과, 3 제주대학교의학전문대학원예방의학교실 Intervention Meta-Analysis Using STATA Software Sung-Ryul Shim 1, In-Soo Shin 2, Jong-Myon Bae 3 1 Institute for Clinical Molecular Biology Research, Soon Chun Hyang University Hospital, Seoul; 2 Department of Education, Jeonju University, Jeonju; 3 Department of Preventive Medicine, Jeju National University School of Medicine, Jeju, Korea Objectives: The aim of this study was to demonstrate the intervention meta-analysis using STATA software. Methods: We conducted a meta-analysis using three types of data that included difference in means in continuous data, hazard ratio in time to event data, and proportion in percent data. STATA commands were metan and metaprop for overall effect size and forest plot, metareg for metaregression, and metafunnel, metabias, and metatrim for publication bias. Results: The estimated overall effect sizes, test for heterogeneity and moderator, and the publication bias were reported using STATA commands by three types of data. Especially authors stressed how to evaluate heterogeneity in intervention meta-analysis. Conclusions: Few of Korean books and articles related to the introducing meta-analysis have been published until now. However, they were not good for the Korean researchers who were non-majored in statistics due to the highly conceptual and theoretical contents. In order to encourage a conduction of meta-analysis, authors summarized commands of the STATA software related to conduct an intervention meta-analysis. Key words: Meta-analysis, Meta-regression, Forest plot, Heterogeneity, Publication bias 서론 체계적문헌고찰 (systematic review) 과메타분석 (meta-analysis) 은대상문헌의선택에있어서투명성을지니며통합된효과크기 (effect size) 를제시하여근거기반의료 (evidence based medicine, EBM) 를위한올바른의사결정을할수있게한다 [1]. 메드라인 (MEDLINE) 통계에따르면메타분석의태동기인 1990년도 272편의연구가검색되는데이수는 2000년 849편, 2010년 4,691편, 그리고 2014년에는 10,427편에이를만큼급증하는추세를나타내며의학에있어서중요한연구방법론으로자리잡고있음을알수있다 [2]. 이러한메타분석의연구방법론에대해국내에서는다수의논문들과 저서가소개하고있다 [3-6]. 그러나이들대부분이개념적이거나또는이론적인접근에만치중하고있어실제이를적용하여논문을작성하기에는구체적인적용방법상의문제점을해결하지못하는한계가있었다. Wallace et al. [7] 은메타분석소프트웨어 STATA, R, MIX, CMA, RevMan, Meta-Analyst 등에대한특성들을보고하였다. RevMan은무료이며기본적인메타분석은가능하지만조절변수의이질성파악등대부분의기능이제한적이다. CMA는상용프로그램이며그래픽사용자입력방식 (graphical user interface) 으로초보자가이용하기에적합하여 RevMan보다는더많은기능을실행하나여전히확장성은제한적이다. R과 STATA는확장성이좋아현재까지개발된중재메타분석에 Corresponding author: Jong-Myon Bae 66 Jejudaehak-ro, Jeju 63243, Korea Tel: , jmbae@jejunu.ac.kr Received: January 12, 2016 Revised: February 24, 2016 Accepted: February 25, 2016 No potential conflict of interest relevant to this article was reported. How to cite this article: Shim SR, Shin IS, Bae JM. Intervention meta-analysis using STATA software. J Health Info Stat 2016;41(1): Doi: It is identical to the Creative Commons Attribution Non-Commercial License ( whichpermit sunrestricted non-commercial use, distribution, and reproduction in any medium, provided the original work is properly cited

2 Sung-Ryul Shim, et al. Figure 1. Flow chart of a quantitative meta-analysis using STATA. *Recommend. 서부터진단검사메타분석에이르기까지다양한분석기능이가능하다. R은무료프로그램으로다양한분야에서분석방법이개발되어있으나기본적으로통계전공자를위한프로그래밍 (programming) 언어로서이를구현하기에는상당한학습이필요하다. 다만최근에는 R을이용한메타분석방법이체계적으로개발되어있어향후그쓰임새가증가할것으로예상한다. STATA는상용프로그램이며모듈 (module) 별로통계방법이개발되는데특히메타분석통계방법들은모두 STATA 저널에서검증을거치기에표준적이고신뢰할수있다. 따라서본연구는통계학이전공분야가아닌대다수의국내연구자들이보다쉽게이해하고실행할수있도록메타분석에특화된 STA- TA 프로그램을이용하여메타분석을실행하는데있어연구자가범하기쉬운오류및실수를최소화하고효과적인연구수행을위한방법및가이드라인 (guideline) 을제시하고자한다 (Figure 1). 연구방법 문헌검색및개별문헌들의질평가 (quality assessment) 는개념적인이해와연구자상호간의협업으로써수행되는것으로 Cochrane handbook [8] 을참고하기바란다. 대상문헌들을선택하였다면개별연구들에서효과크기를추출하여이를 STATA에맞게코딩 (coding) 한다 (Table 1). 예제자료본연구에서는효과크기가평균의차이 (difference in means) 로나타내어지는연속형자료 (continuous data), 위험도비 (odds ratio, risk ratio, 또는 hazard ratio) 로나타내어지는생존형자료 (time to event data), 그리고백분율로나타내어지는유병률자료 (proportion or rate data) 를가지고통합효과크기 (overall effect size) 산출, 이질성 (heterogeneity) 확인및원인파악, 그리고출판편향 (publication bias) 검토를설명토록한다 (Table 1). 그렇지만일반적으로의학에서는상관계수 (correlation coefficient, r) 로나타내어지는효과크기의활용도는떨어짐으로이에대한설명은제외하였다. 연속형자료척수손상동물모델에서줄기세포치료에따른방광기능효과를메타분석한연구에서결과지표배뇨시압력 (voiding pressure) 을예제로사용하였다 [9]. 전체연구수는 11개였고실험군 94개와대조군 93개로이루어졌다. 부집단 (subgroup) 1은 contusion model 그리고 0은 transection and hemisection model로구분하였다. 연구결과에따르면 WMD (weighted mean difference) 는 (95% CI: -8.22, -4.48) 였고 Higgins I 2 는 77.2% 그리고 Cochrane s Q statistics p-value는 0.001이었다. 생존형자료 Cetuximab and panitumumab의대장암치료시 B-Raf proto-onco

3 STATA Commands for Meta-analysis Table 1. A raw data for demonstrating the results with applying commands of STATA Continuous outcomes Time to event outcomes (same coding as dichotomous outcomes) Proportion outcomes Study T_M T_SD T_N C_M C_SD C_N G 1 Study lnhr selnhr G 2 Study F N G 3 Mitsui Bokemeyer Shim Mitsui 2005a Douillard Kim Mitsui 2005b Karapetis Lee Mitsui Seymour Park Telmeltas 2009_ Peeters Bae ,071 0 Telmeltas 2009_ Stintzing Jin WBPark 2010_ Kwon WBPark 2010_ Shim YHu Huh YingJin 2011_ Lee YingJin 2011_ Lee ,071 2 Sources in Kim, et al. [9] for continuous outcomes, Pietrantonio, et al. [10] for time to event outcomes, and sample data for proportion outcomes. T_M, treatment group mean; T_SD, treatment group standard deviation; T_N, treatment group sample size; C_M, control group mean; C_SD, control group standard deviation; C_N, control group sample size; lnhr, log hazard ratio; selnhr, standard error log hazard ratio; F, frequency; N, samplesize; G, subgroup. 1 1: contusion model, 0: transection and hemisection model. 2 1: used fluorouracil, 0: no used fluorouracil. 3 0: 1st generation, 1: 2nd generation, 2: 3rd generation. gene (BRAF) 변이에따른위험도 (Hazard Ratio, HR) 를메타분석한연구에서결과지표전체생존율 (overall survival, OS) 을예제로사용하였다 [10]. 전체연구수는 6개였고부집단 1은 fluorouracil 을사용한군그리고 0은 fluorouracil을사용하지않은군으로구분하였다. 연구결과에따르면 HR은 0.91 (95% CI: 0.62, 1.34) 였고 Higgins I 2 는 50% 그리고 Cochrane s Q statistics p-value 는 0.007이었다. 유병률자료요로감염환자에서특정항생제 (cephalosporin) 의내성률 (resistance) 을메타분석한임의의자료를예제로사용하였다. 전체연구수는 11 개였고사건 (event) 발생수는 347.3개그리고전체표본수는 3,517개로이루어졌다. 부집단은항생제의세대를나눈것으로 0은 1세대, 1은 2 세대, 그리고 2는 3세대로구분하였다. 메타분석모형선택중재법메타분석 (intervention meta-analysis) 의통합효과크기를산출하기위해일반적으로고정효과모형 (fixed effects model) 과변량효과모형 (random effects model) 이사용된다. 고정효과모형고정효과모형은개별연구들의중재 (intervention) 에따른공통효과크기 (common effect size) 가존재하며연구들간의효과차이는표본 추출오차 (sampling error) 때문이라가정한다. 따라서고정효과모형은대상집단과중재방법이동질하다고판단될경우고려할수있으며일반적으로이질성을나타내는지표인 Higgins I 2 가 50% 이하인경우또는 Cochrane Q statisticsp-value가유의하지않은경우 (p 0.1) 고려될수있다 (Cochrane handbook; Identifying and measuring heterogeneity) [8]. 고정효과모형은 inverse variance method, mantel-haenzel method, 그리고 peto method가있으며일반적으로 inverse variance method가사용된다. 개별모형에대한상세한내용은 Cochrane handbook [8] 9장을참고하기바란다. STATA에서는고정효과모형 (mantel-haenzel method) 이 default로설정되어있으며옵션에 fixed 를넣어도동일한결과를나타낸다. 변량효과모형변량효과모형은개별연구들의중재효과에대해공통의효과크기를추정하기보다는효과크기의분포 (distribution) 를추정하며정규분포를따른다고가정한다. 따라서연구들간의효과차이는표본추출오차와연구들간의실체적차이 (true difference) 를동시에고려한다 [11]. 일반적으로연구들간의이질성이높을때 (Higgins I 2 >50% or Cochrane s Q statistics p <0.1) 변량효과모형이고정효과모형보다보수적인추정치를제공하여권고된다 (Cochrane handbook; Incorporation heterogeneity into random-effects models) [8]. 연구들간의실체적차이 (tau) 가적거나 0인경우두모형은동일한값을나타낸다

4 Sung-Ryul Shim, et al. Table 2. Related commands on the STATA for intervention meta-analysis Issues Overall effect size and forest plot If the effect sizes and standard errors of each study were calculated Raw data only available of each study Continuous data Time to event data Proportion data Meta-regression for heterogeneity Publication bias Funnel plot Begg and Mazumdar s rank correlation test Egger s linear regression method test Trim and fill Examples of the related commands Metan var(effect size) var(standard error), options Metan var(sample size of tg) var(mean of tg) var(standard deviation of tg) var(sample size of cg) var(mean of cg) var(standard deviation of cg), options Metan var(effect size) var(standard error), options Metaprop var(event size) var(total sample size), options Metareg var(effect size) var(moderator), wsse (standard error) Metafunnel var(effect size) var(standard error) Metabias var(effect size) var(standard error), begg Metabias var(effect size) var(standard error), egger Metatrim var(effect size) var(standard error), funnel Options (eg. fixed, random, by, eform, and/or label etc.). var, variable; tg, treatment group; cg, control group. STATA 명령어개별연구들이여러통계값들 [eg. d(mean difference), OR/RR/HR, r(correlation)] 로나타내어져있다면이들을표준화된하나의효과크기로상호변환할수있다 (Chapter 7. Converting among effect sizes ) [1]. 이렇게직접개별연구들의효과크기를계산하였다면자료의유형에상관없이 STATA 명령어 metan에효과크기와표준오차 (standard error) 를순서대로넣어주면통합된효과크기계산이가능하다. Table 2는메타분석을위한 STATA 명령어를정리한것이다. 통합효과크기계산과숲그림 (forest plot) 을작성을위한 metan과 metaprop, 이질성의원인을파악하기위한메타회귀분석 (meta-regression) 은 metareg, 그리고출판편향 (publication bias) 을확인하기위해서는 metafunnel, metabias, 그리고 metatrim의명령어를사용한다 (Table 2). 상기옵션들은기본적으로효과크기와표준오차의변수를순서대로입력하며각자료의유형과옵션은아래에서상세히설명토록한다. 생존형자료와이분형자료 (dichotomous data) 의경우효과크기가비 (ratio) 형태인 HR 또는 OR로나타내어지기에메타분석과정이동일하다. 물론이분형자료의경우개별연구의 treatment 그룹과 control 그룹의전체표본수와사건발생수를모두추출할경우 STATA에서보다간단히처리할수있으며이에대한설명은참고문헌으로대신한다 [12]. 한다. Treatment 그룹과 control 그룹의순서는효과크기를해석하기쉬운방향으로임의로넣어도무관하다. Table 1의연속형자료를아래의명령어로이용하여나온결과가 Figure 2이며이는원저와동일하다. metan t_n t_m t_sd c_n c_m c_sd, label(namevar = study) nostandard random by(g) textsize(130) 쉼표뒤로는옵션들이위치하는데옵션의순서는상관이없으며 label(name = study) 는 study라는변수를개별연구들의이름으로설정하는것이고 nostandard는 standardized mean difference (SMD) 를사용하지말고 weighted mean difference (WMD) 를사용하라는것이고, random은이질성추정을 Mantel-Haenszel 모형에근거한변량효과모형을사용하며, by(g) 은부집단별로결과를제시하며마지막으로 textsize(130) 은 forest plot 내의글자크기를나타내는것이다. 또한 metan 명령어는 forest plot을자동으로제시한다. 이외출력된내용에대한설명은참고문헌으로대신한다 [12]. Figure 2에서전체효과크기는 (95% CI: , ) 로서통계적으로유의한값을나타낸다. Higgins I 2 는 77.2% 그리고 Cochrane s Q statistics p-value는 미만으로서이질성이존재한다는것을알수있다. 부집단분석결과 contusion model과 transection and hemisection model 에따라차이를보이고있으며 transection and 결과 hemisection model 이더낮은값을나타낸다. Table 1 의연속형자료인 Kim et al. [9] 원저는줄기세포실험데이터 통합효과크기산출연속형자료연속형자료일경우기본적으로 treatment 그룹과 control 그룹의표본수, 평균, 그리고표준편차 (standard deviation) 를순서대로넣어야 로서실험실적환경을제외하고는다른변수들이동일하다고가정하였으며동일한결과지표를사용하였기에 WMD를사용하였다. WMD 는 Cochrane에서는 mean difference 또는 difference in means으로지칭하며개별연구들에서동일한단위 (scale) 로얻어낸결괏값일경우 126

5 STATA Commands for Meta-analysis 사용이가능하며해석이용이하다는장점이있다 [Cochrane handbook; The mean difference (or difference in means)] [8]. 그렇지 만 Mean difference 는비교그룹간절댓값의차이를표본수와표준오 변량효과모형을설정하는명령어는다음과같다. metan t_n t_m t_sd c_n c_m c_sd, label(namevar = study) hedges ran dom by(g) 차의가중치로추정한것으로여기에개별연구들을각각의공통표준 편차 (pooled standard deviation) 로나누어준 SMD와의용어혼동가능성이있으므로사용을권고하지않는다. 따라서개별연구들의효과크기를종합하는메타분석방법론에충실하고연구들간의표준화된지수로서비교하기위해서는 SMD를사용하는것이바람직하다 [13]. 또한표본이작을경우이를 SMD를교정해주는 Hedges g를사용하여야한다 [14]. 따라서상기자료를일반적인 SMD (hedges g) 와 생존형자료명령어는동일하게 metan을사용하며 HR과표준오차를순서대로넣는다. 일반적으로개별연구들이 HR과 95% CI를제시하는데 metan 명령어사용을위해서는 HR과 95% CI에자연로그를취하여로그변환된 HR과표준오차를계산하여야한다. 로그변환을시키는이유는 HR 의경우지수변환된상태로서 0 이하한설정으로되어있어이를풀어 Study WMD [95% CI] % weight 1 Mitsui Mitsui 2005a Mitsui 2005b Mitsui WBPark 2010_ WBPark 2010_ YingJin 2011_ YingJin 2011_ Sub-total D+L pooled WMD Telmeltas 2009_ Telmeltas 2009_ YHu Sub-total D+L pooled WMD Overall D+L pooled WMD Study ID WMD (95% CI) % weight 1 Mitsui (-22.99, -8.61) 4.69 Mitsui 2005a (-8.67, -6.73) Mitsui 2005b (-7.07, -3.93) Mitsui (-10.36, -6.04) WBPark 2010_ (-3.56, 3.76) 9.68 WBPark 2010_ (-3.80, 8.24) 5.92 YingJin 2011_ (-8.87, -3.35) YingJin 2011_ (-8.12, -2.11) Subtotal (I-squared = 81.2%, p = 0.000) (-7.73, -3.70) Telmeltas 2009_ (-27.38, -9.42) 3.38 Telmeltas 2009_ (-18.42, 1.02) 2.98 YHu (-10.72, -3.12) 9.40 Subtotal (I-squared = 62.5%, p = 0.070) (-17.51, -3.70) Overall (I-squared = 62.5%, p = 0.000) NOTE: Weights ar from random effects analysis (-8.22, -4.48) Summary effect size A Forest plot B Figure 2. Overall effect size and forest plot for continuous outcomes. CI, confidence interval. Study ES [95% CI] % weight Bokemeyer Douillard Karapetis Seymour Peeters Stintzing D+L pooled ES Study ID ES (95% CI) % weight Bokemeyer (0.36, 1.06) Douillard (046, 1.76) Karapetis (0.20, 3.56) 5.97 Seymour (1.10, 3.08) Peeters (0.32, 1.28) Stintzing (0.47, 1.61) Overall (I-squared = 62.5%, p = 0.070) NOTE: Weights ar from random effects analysis (0.62, 1.34) Summary effect size A Forest plot B Figure 3. Overall effect size and forest plot for Time to event outcomes. CI, confidence interval

6 Sung-Ryul Shim, et al. 주기위한것이다 (Chapter 5. Effect sizes based on binary data.33-39) [1]. 표준오차는 HR의 95% CI 상 하한에자연로그를취한후상한에서하한을뺀값을 3.92로나누어주면쉽게얻을수있다 [(ln(95cihigh) - ln(95cilow))/3.92]. Table 1의생존형자료를아래의명령어로이용하여나온결과가 Figure 3이며이는원저와동일하다. metan lnhr selnhr, label(namevar = study) eform random textsize(150) lnhr은로그변환된 HR이며 selnhr은표준오차이다. 쉼표뒤로는옵션들이위치하는데나머지는연속형자료의설명과동일하며 eform은로그변환된수치로써계산을하기때문에 HR로의해석을위해서는다시이를지수변환시키라는명령어이다. 만약생존분석연구에서 HR과관련통계값을충분히제시하지않고 Kaplan-Meier 생존곡선만있는경우생존곡선에서 HR을계산할수있는방법도있으니참고하기바란다 [15]. Figure 3에서전체효과크기 HR이 (95% CI: 0.617, 1.339) 로서통계적으로유의하지않은값을나타낸다. Pietrantonio, et al. [10] 원저에서는부집단분석결과를따로제시하지않고있으나본연구에서는 fluorouracil의사용여부에따른부집단분석을옵션 by를이용하여아래와같이실시하였다. metan lnhr selnhr, label(namevar = study) eform random by(g) fluorouracil을사용한그룹의 HR은 (95% CI: 0.540, 1.005) 로위험도가약 26% 감소하였고통계적으로경계선값을 (p = 0.054) 나타내었다. 반면 fluorouracil을사용하지않은그룹의 HR은 (95% CI: 1.028, 2.748) 을나타내어위험도가통계적으로유의하게약 68% 증가하였다 (p = 0.038). 따라서 fluorouracil 의치료여부에따른 HR의차이가있음을확인할수있다. 유병률자료유병률자료는통상비율 (proportion or rate) 로표시되며질병의이환율, 사망률, 항생제내성률등을나타내는지표로자주활용되고있다. 메타분석을위해서개별연구들의효과크기가다른통계량이제시되지않고비율로만나타내어져있거나비율자체가연구의주요관심사일때는이에대한통합효과크기를계산하여야한다. 명령어는 metaprop 를사용하며사건발생빈도 (event size) 와전체표본수 (sample size) 를순서대로넣는다. Table 1의유병률자료를아래의명령어로이용하여나온결과가 Figure 4이다. metaprop f n, random by(g) label(namevar = study) 쉼표뒤로는옵션들이위치하는데 random은 inverse variance 모형 Study ES [95% CI] % weight 1st Shim Kim Lee Park Bae Sub-total Random pooled ES nd Jin Kwon Shim Sub-total Random pooled ES rd Huh Lee Lee Sub-total Random pooled ES Overall Random pooled ES Study ES (95% CI) % weight 1st Shim (0.09, 0.19) 9.07 Kim (0.07, 0.16) 9.30 Lee (0.22, 0.34) 8.37 Park (0.17, 0.28) 8.70 Bae (0.09, 0.12) Subtotal (P2=90.81%, p =0.00) 0.16 (0.10, 0.23) nd Jin (0.06, 0.27) 5.74 Kwon (0.04, 0.18) 7.89 Shim (0.03, 0.09) 9.86 Subtotal 0.07 (0.03, 0.11) (P2 = 22.39%, p = 0.28) 3rd Huh (0.00, 0.09) 9.80 Lee (0.00, 0.03) Lee (0.04, 0.07) Subtotal 0.03 (-0.01, 0.06) (P2 = 93.55%, p = 0.00) Heterogenity between groups. p = Overall (P2=95.57%, p =0.00) 0.10 (0.06, 0.14) Summary effect size A Forest plot B Figure 4. Overall effect size and forest plot for Proportion outcomes. CI, confidence interval

7 STATA Commands for Meta-analysis 에근거한변량효과모형을사용한다는뜻이다. Figure 4에서전체효과크기항생제내성률은 0.10 (95% CI: 0.06, 0.14) 이다. Higgins I 2 는 95.57% 그리고 Cochrane s Q statistics p-value는 0.01 미만으로서이질성이존재한다는것을알수있다. 항생제세대별부집단분석결과를살펴보면 1세대는 0.16 (95% CI: 0.10, 0.23), 2세대는 0.07 (95% CI: 0.03, 0.11), 그리고 3세대는 0.03 (95% CI: -0.01, 0.06) 으로서세대가높아질수록항생제내성률이감소하는경향을나타낸다. 이질성확인메타분석에서얻어낸통합효과크기를제대로해석하려면연구들간의이질성유무를확인하고만약유의한조절변수 (moderator) 가있다면이를검정하고보고하여야한다. 이러한이질성의원인은우연 (chance), 연구설계 (study design) 의차이, 연구환경, 그리고표본집단의인구사회학적요인에이르기까지매우다양하다. 정량화시킨다. 0% 에서 40% 는이질성이중요하지않을수있으며 (might not be important), 30% 에서 60% 는중간이질성 (moderate heterogeneity), 50% 에서 90% 는중대한이질성 (substantial heterogeneity), 그리고 75% 에서 100% 는심각한이질성 (considerable heterogeneity) 을나타낸다 (Cochrane handbook; Identifying and measuring heterogeneity) [8]. 연속형자료의숲그림 (Figure 2) 에서 Cochrane Q statistics p <0.001 로귀무가설을기각하고, Higgins I 2 는 77.2% 로서중대한이질성이존재한다고판단할수있다. 생존형자료의숲그림 (Figure 3) 에서 Cochrane Q statistics p = 0.073로귀무가설을기각하고, Higgins I 2 는 50.5% 로서중간이질성이존재한다고판단할수있다. 유병률자료의숲그림 (Figure 4) 에서 Cochrane Q statistics p < 0.01로귀무가설을기각하고, Higgins I 2 는 95.57% 로서심각한이질성이존재한다고판단할수있다. 시각적확인 : 숲그림 (forest plot) 개별연구들의통합효과크기와 95% CI를제시하여연구내변동과연구간변동을쉽게파악할수있도록하는것이숲그림이다 [16]. STATA에서는 metan 명령어의기본변수만넣어주면숲그림은자동으로생성된다 (Figures 2-4). 연속형자료의숲그림에서연구내변량이큰것은 Mitsui 2003, WBPark 2010_2, Telmeltas 2009_1, 그리고 Telmeltas 2009_2인것을알수있고연구간변량이큰것은 Mitsui 2003, WBPark 2010_2, 그리고 Telmeltas 2009_1인것을알수있다 (Figure 2). 생존형자료의숲그림에서연구내변량이큰것은 Karapetis 2014 이며연구간변량이큰것은 Seymour 2013인것을알수있다 (Figure 3). 유병률자료의숲그림에서연구내변량이큰것은 Jin 2003이며연구간변량이큰것은 Huh 2003, Lee 2004, 그리고 Lee 2008인것을알수있다 (Figure 4). 이처럼숲그림을통해통합효과크기각각에대해개별연구들의점추정과구간추정값을전체적으로볼수있다. 통계적검정개별연구들의이질성을통계적으로검정하는방법으로는 Higgins I 2 와 Cochrane Q statistics가있다. Cochrane Q statistics은개별연구들의중재효과가통합효과값으로부터얼마나멀리떨어져있는지를 χ 2 검정을통해실시하는것이다. 귀무가설은이질성이있다는근거가없다는것이며 p-value 0.1을통계적유의수준으로사용한다 (Cochrane handbook; Identifying and measuring heterogeneity) [8]. Higgins I 2 는 Cochrane Q statistics에서자유도 (degree of freedom) 를뺀것을다시 Cochrane Q statistics으로나누어준값으로이질성을일관성있게 이질성원인파악 (subgroup analysis and meta-regression) 숲그림을이용한시각적확인과 Cochrane Q statistics와 Higgins I 2 이용한통계적검정으로부터이질성이의심된다면이질성의원인을파악하기위하여부집단분석과메타회귀분석을실시한다. 부집단분석은앞선생존형자료에서 fluorouracil을사용한그룹의 HR (95% CI: 0.540, 1.005) 과사용하지않은그룹의 HR (95% CI: 1.028, 2.748) 이통계적으로유의한차이가있음을이미나타내었고연속형자료 (Figure 2) 와유병률자료 (Figure 4) 에서도각각의부집단에따른효과크기차이가있음을알수있었다. 따라서부집단에해당하는변수들이실제조절변수 (moderator) 로작용하는지를알아보기위하여메타회귀분석을실시한다. 메타회귀분석은일반적으로개별연구들이 10개이상일때추천된다 (Cochrane handbook; Meta-regression) [8]. 예제로사용된두연구의원저 [9,10] 는메타회귀분석까지는실시되지않았었다. 따라서본연구에서는중대한이질성요인이라고판단되는변수를 (1: contusion model, 0: transection and hemisection model [9]; 1: used fluorouracil, 0: no used fluorouracil [10]) 임의로설정하여메타회귀분석을실시하였다. 유병률자료는항생제세대를 (0: 1st, 1: 2nd, 2: 3rd generation) 부집단변수로설정하였다. STATA에서는 metareg이라는명령어를제공하는데앞에서살펴본대로 metan을사용하여통합효과크기를먼저산출한후이어서분석하여야한다. 메타회귀분석방법에는 residual maximum likelihood (REML), method of moments (MM), 그리고 empirical bayes (EB) 등이있으며 STATA에서는 REML이 default로설정되어있다. Table 1의연속형자료를아래의명령어를차례대로사용하여효과 129

8 Sung-Ryul Shim, et al. Table 3. Effects of moderators using meta-regression methods k Coef. 1 SE t p 2 95% CI Continuous outcomes Time to event outcomes Proportion outcomes Continousoutcomes (0: transection and hemisection model, 1: contusion model)[9]. Time to eventoutcomes (0: no used fluorouracil, 1: used fluorouracil) [10]. Proportion outcomes (0: 1st, 1: 2nd, 2: 3rd generation). k, number of effect sizes; SE, standard error; CI, confidence interval. 1 Regression coefficient. 2 p-values from random effect meta-regression using restricted maximum likelihood (REML) for continuous moderators. 크기와표준오차를구한다음메타회귀분석을실시한결과가 Table 3 이다. metan t_n t_m t_sd c_n c_m c_sd, label(namevar = study) hedges ran dom by(g) metareg _ES g, wsse( _sees) Metareg 명령어는효과크기 (_ES: metan에서자동으로생성됨 ) 와이질성의원인으로의심되는변수 (subgroup; 1: contusion model, 0: transection and hemisection model) 를순서대로넣는다. 쉼표뒤옵션 wsse (within study standard error) 는연구내표준오차변수 (_sees: metan에서자동으로생성됨 ) 를지정한것이다. Table 3의연속형자료에서 p값 (p = 0.727) 이통계적으로유의하지않아 contusion model과 transection and hemisection model변수는이질성요인이아닌것으로나타났다. Table 1의생존형자료를아래의명령어를차례대로사용하여효과크기와표준오차를구한다음메타회귀분석을실시한다 (Table 3). metan lnhr selnhr, label(namevar = study) eform random by(g) metareg lnhr g, wsse(selnhr) eform Metareg 명령어에효과크기 (lnhr; 로그변환된 HR) 와이질성의원인으로의심되는변수 (subgroup; 1: used fluorouracil, 0: no used fluorouracil) 를순서대로넣는다. 쉼표뒤옵션 wsse에는표준오차변수 (selnhr) 를넣는다. Table 3의생존형자료에서 p값 (p = 0.048) 이통계적으로유의하여 fluorouracil의사용여부가이질성의요인임을알수있다. 따라서부집단분석 (fluorouracil 사용한그룹, HR 0.74; 미사용한그룹, HR 1.68) 에서도알수있었듯이 fluorouracil의치료는 HR을낮추는요인으로해석할수있다. Table 1의유병률자료를아래의명령어를차례대로사용하여효과크기와표준오차를구한다음메타회귀분석을실시한다 (Table 3). metaprop f n, random by(g) label(namevar = study) metareg _ES g, wsse( _sees) Metareg 명령어에효과크기 (_ES: metaprop에서자동으로생성됨 ) 와이질성의원인으로의심되는변수 (subgroup; 0: 1st, 1: 2nd, 2: 3rd generation) 를순서대로넣는다. 쉼표뒤옵션 wsse에표준오차변수 (_sees) 를넣는다. Table 3의유병률자료에서 p값 (p = 0.006) 이통계적으로유의하여항생제세대가이질성의요인임을알수있다. 따라서부집단분석 ( 항생제내성률 1세대 0.16, 2세대 0.07, 그리고 3세대 0.03) 에서도알수있었듯이항생제세대가항생제내성률조절변수임을알수있다. 이외출력된내용에대한설명은참고문헌으로대신한다 [17]. 출판편향 (publication bias) 검토출판편향은개별연구들의특성과결과에따라연구가출판되거나출판되지않을오류이다 [18]. 일반적으로통계적유의한연구결과일경우더욱출판될가능성이높기때문에발생하는데, 이러한출판편향을고려하여해당메타연구의결과가과대또는과소추정되지는않았는지확인하여야한다. 시각적확인 (funnel plot) STATA에서는 metafunnel 명령어를제공하는데앞에서살펴본대로 metan을사용하여통합효과크기를먼저산출한후이어서분석하여야한다. Table 1의연속형자료, 생존형자료, 그리고유병률자료를아래의명령어를사용하여각각의 funnel plot을나타내었다 (Figure 5). metafunnel _ES _sees, xtitle(standard mean difference; Hedges s g) ytitle(standard error of SMD) metafunnel lnhr selnhr, xtitle(log hazard ratio) ytitle(standard error of log HR) metafunnel _ES _sees, xtitle(resistance of antibiotics) ytitle (Standard error) metafunnel 명령어에효과크기 (_ES/lnhr) 와표준오차 (_sees/selnhr) 를순서대로넣어준다. 쉼표뒤의 xtitle과 ytitle은각각의이름을지정한것이다. Funnel plot의 Y축은표본크기 ( 표준오차 ) 를 X축은효과크기를제시한다. 일반적으로작은규모의연구들은아래쪽에넓게분포되며큰규모의연구들은깔때기안상단에좁게분포된다. 따라서깔때기안상단에좌우대칭으로골고루분포되어있다면출판편향은적다고판 130

9 STATA Commands for Meta-analysis Funnel plot with pseudo 95% confidence limits Funnel plot with pseudo 95% confidence limits Funnel plot with pseudo 95% confidence limits Standard error of SMD Standard error of log HR Standard mean difference; Hedges s g Log hazard ratio Resistance of antibiotics Continous outcomes A Dichotomous outcomes B Proportion outcomes C Standard error Figure 5. Funnel plot for publication bias. 단할수있다. 연속형자료의경우 (Figure 5A) X축의중앙 SMD를기준으로깔때기바깥좌측으로 3개, 우측으로 2개의연구가, 깔때기안좌측으로 4 개, 우측으로 2개의연구가분포한다. 시각적으로판단하기에출판편향이있을것으로판단된다. 생존형자료의경우 (Figure 5B) X축의중앙 ln hazard ratio를기준으로깔때기바깥우측으로 1개의연구가, 깔때기안좌측으로 5개의연구가모두몰려있어이또한좋은분포가아니다. 유병률자료의경우 (Figure 5C) X축의중앙 Resistance of antibiotics를기준으로깔때기바깥좌측으로 1개, 우측으로 5개의연구가, 깔때기안좌측으로 1개, 우측으로 4개의연구가모두몰려있어대칭성을나타내지못하고있다. 통계적검정출판편향을통계적으로검정하는일반적인방법은 Egger s linear regression method test [19], Begg and Mazumdar s rank correlation test [20], 그리고 Orwin s failsafe N이있다. Egger s linear regression method test가 Begg and Mazumdar s rank correlation test보다효과크기의실제추정치를더정확히추정한다고보고하고있으며 Failsafe N은통계적방법에만치중하고있어실제적인임상적유의성을판단하지못하여사용을권고하지않는다 [21]. 그러나통계적검정은출판편향의영향보다는연구의수가적을경우 (small study effect) 이를정확히검정하지못하므로 Cochrane에서는권고하지않는다. Begg and Mazumdar s rank correlation test 개별연구들의표준화된효과크기와표준오차와의상관관계를보정된순위상관 (rank correlation) 으로검정한다. 순위상관검정이유의하지않다면출판편향이없음을나타낸다. STATA에서는 metabias 명령어를제공하는데앞에서살펴본대로 metan을사용하여통합효과크기를먼저산출한후이어서분석하여야한다. Table 1의연속형자료와유병률자료를아래의명령어를사용하여 Begg and Mazumdar s rank correlation test를실시하였다. metabias _ES _sees, begg 연속형자료의경우 continuity를보정한 Z는 2.96, p값이 0.003으로출판편향이있음을확인할수있으며유병률자료의경우 continuity 를보정한 Z는 1.40, p값이 로출판편향이있음을증명할수없다. Table 1의생존형자료를아래의명령어를사용하여분석한결과 continuity를보정한 Z는 0.00, p값이 1.000으로출판편향이있음을증명할수없다. metabias lnhr selnhr, begg Egger s linear regression method test 중재효과의표준오차에대한개별연구들의효과크기관계를회귀식으로나타낸것으로귀무가설은회귀식의초깃값 (intercept) 은우연에의한결과로서출판편향이있음을증명할수없다는것이다 [19]. STATA 에서는 metabias 명령어를제공하는데 metan을사용하여통합효과크기를먼저산출한후이어서분석하여야한다. Table 1의연속형자료와유병률자료를아래의명령어를사용하여 Egger s linear regression method test를실시한결과를 Table 4에나타내었다. metabias _ES _sees, egger 연속형자료의경우전체 11개의연구를사용하였고 bias항목의 Coef. 가 -9.23으로초깃값 (intercept) 을나타내며해당 p값이 0.001로서귀무가설을기각하여출판편향이있음을확인할수있다. 유병률자료의경우전체 11개의연구를사용하였고 bias항목의 Coef. 가 -5.18로초깃값 (intercept) 을나타내며해당 p값이 0.019로서귀무가설을기각하 131

10 Sung-Ryul Shim, et al. Table 4. Egger s linear regression method test for small-study effects Std_Eff Coef. SE t p >t [95% CI] Continuous Time to event Proportion Number of studies=11 Slope Bias Test of H0: no small-study effects p =0.000 Number of studies=6 Slope Bias Test of H0: no small-study effects p =0.636 Number of studies=11 Slope Bias Test of H0: no small-study effects p =0.019 SE, standard deviation; CI, confidence interval. Method Meta-analysis Pooled Est 95% CI Asymptotic No. of studies Lower Upper Z value p-value Fixed Random Filled meta-analysis Fixed Random metabias lnhr selnhr, egger 출판편향추정치보정 : trim and fill Funnel plot상에서대상연구들이비대칭성을보인다면효과크기는한쪽으로편향될가능성이있다. 따라서민감도분석 (sensitivity analysis) 의일환으로대칭성이충족된상태에서보정된효과크기를확인하는것이 trim and fill이다. Summary statistics Filled funnel plot with pseudo 95% confidence limits A Trim and fill 은메타분석에있어서출판편향을비모수적방법으로 효과크기추정치를보정한다 [22]. 방법은 Funnel plot 상에서비대칭성 을나타내는연구들을절삭 (trim) 하여남아있는대칭성을나타내는 0.4 연구들의효과크기를계산한후최초절삭된연구들과대비하여대칭 Theta, filled 을나타내는가상의연구들을채워 (fill) 넣는다. 그런다음기존연구들과가상의연구들을모두포함하여통합효과크기를계산하게되면시각적으로비대칭성을극복한상태에서보정된통합효과크기를추정할수있다. STATA에서는 metatrim 명령어를제공하는데 metan을사용하여먼저통합효과크기를산출한후이어서 metatrim 명령어에효과크기 (_ES/lnhr) 와표준오차 (_sees/selnhr) 를순서대로넣어준다 s.e. of: theta. filled Trim and fill plot Figure 6. Trim and fill. CI, confidence interval. 여출판편향이있음을확인할수있다. B 연속형자료와유병률자료의명령어는동일하며아래와같다. metatrim _ES _sees, reffect funnel 생존형자료의명령어는아래와같다. metatrim lnhr selnhr, reffect eform funnel 옵션들중 reffect는변량효과모형을사용하는것 고정효과모형이 default임 이고, funnel은채워지는가상의연구들을 funnel plot Table 1의생존형자료를아래의명령어를사용하여분석한결과전체 6개의연구를사용하였고 bias항목의 Coef. 가 으로초깃값 (intercept) 을나타내며해당 p값이 0.636로서귀무가설을기각하지않아출판편향이있음을증명할수없다. 상에표시하라는것이다. 본연구의예제로주어진 3가지유형 ( 연속형, 생존형, 그리고유병률 ) 의자료는 trim and fill에따른통합효과크기의변동이없었다. 하지만실행방법을보여주고자 Table 1의유병률자료를고정효과모형으로분 132

11 STATA Commands for Meta-analysis 석한후 trim and fill을실시하여 Figure 6에나타내었다. Figure 5C와비교해서보면전체 11개의연구들중우측으로치우친 5개를절삭한후이에대응하는가상의연구를좌측에 5개를추가하여전체 16개의연구들로보정된통합효과크기를계산하였다. 최초메타분석을실시했을때고정효과모형의값은 (95% CI: 0.040, 0.054) 이었지만 trim and fill에따른보정된값은 (95% CI: 0.030, 0.043) 로감소하였다. 요약및제한점지금까지 STATA를이용한메타분석의수행과정을자료유형별로효과크기산출, 숲그림작성, 이질성확인과메타회귀분석, 그리고출판편향등을위한프로그램명령어들을살펴보았다. 특히메타분석은개별연구들의효과크기산출이가장먼저이루어진후실행되므로효과크기에대한이해가가장중요하다고할수있다. 본연구의제한점은실질적인메타분석수행방법에집중하였기에메타분석의수행과정중통계모형선택, 통계모형내적절한옵션의선택, 그리고출판편향에대한통계적검정방법의선택등을깊이있게다루지못하였다는것이다. 특히유병률자료의경우출판편향에대한통계검정방법에따라상이한결과를나타내어혼란이있을것으로판단된다. 따라서개별연구들의특성에적합한분석방법의선택이이루어져야하며이를위해서는국문으로실시된선행연구들을 [3-6] 반드시살펴보기를권한다. 결론및제언 본연구는통계학을전공하지않은일반연구자들도쉽게수행할수있도록메타분석의실질적수행방법에집중하고통계학의이론적접근은최소화하고자노력하였다. 하지만본연구는 STATA 프로그램으로메타분석을수행하는방법만을정리한것으로서메타분석에대한심도있는이론적접근에는미치지못하였다. 따라서이를본인의연구로올바르게해석하기위해서는전문가의견청취또한소홀히하여서는안될것이다. 지금까지살펴본중재메타분석외에도진단검사메타분석 (diagnostic test accuracy meta-analysis) [23] 방법도 STATA를이용하여쉽게접근할수있으니참고하기바라며이를통하여많은국내연구자들이보다쉽게메타분석을수행할수있게되어관련연구들이활성화되기를바란다. 감사의글 본원고는메타분석연구회 ( 신인수회장 ) 연구활동의결과물임을밝힙니다. REFERENCES 1. Borensein M, Hedges LV, Higgins JPT, Rothstein HR. Introduction to Meta-Analysis.preface West Sussex, UK: John Wiley & Sons Ltd.; p National Center for Biotechnology Information, U.S. National Library of Medicine. Search Meta-Analysis [Publication Type] Filters: Publication date from each year/01/01 to each year/12/31]. Available at [accessed on 29 December 2015] 3. Lee JY. Meta-analysis. J Korean Soc Endocrinol 2008;23(6): (Korean). 4. Ahn HS. Kim HJ. An introduction to systematic review. J Korean Med Assoc 2014;57(1):49-59 (Korean). 5. Kim CB. Meta-analysis. Educational lecture supplement. Korean Soc Prev Med 2002;July: (Korean). 6. Kim SY, Park JE, Lee YJ, Son HJ, Suh HS, Shin CM. NECA s guideline for undertaking systematic reviews and meta-analysis for intervention. Seoul: National Evidence-based Healthcare Collaborating Agency; 2001 (Korean). 7. Wallace BC, Schmid CH, Lau J, Trikalinos TA. Meta-Analyst: software for meta-analysis of binary, continuous and diagnostic data. BMC Med Res Methodol. 2009;9: Higgins JPT, Green S. editors Cochrane handbook for systematic reviews of interventions. Version 5.1.0: The Cochrane Collaboration. Available at [accessed on 29 December 2015]. 9. Kim JH, Shim SR, Doo SW, Yang WJ, Yoo BW, Kim JM, et al. Bladder recovery by stem cell based cell therapy in the bladder dysfunction induced by spinal cord injury: systematic review and meta-analysis. PLoS One 2015;10(3):e Pietrantonio F, Petrelli F, Coinu A, Di Bartolomeo M, Borgonovo K, Maggi C, et al. Predictive role of BRAF mutations in patients with advanced colorectal cancer receiving cetuximab and panitumumab: a meta-analysis. Eur J Cancer 2015;51(5): DerSimonian R, Laird N. Meta-analysis in clinical trials. Control Clin Trials 1986;7: Harris R, Bradburn M, Deeks J, Harbord R, Altman D, Sterne J. Metan: fixed- and random-effects meta-analysis. Stata J 2008;8(1):

12 Sung-Ryul Shim, et al. 13. Cohen J. Statistical power analysis for the behavioral sciences. Hillside, NJ: Lawrence Erlbaum Associates; Hedges L, Olkin I. Statistical methods for meta-analysis. San Diego, CA: Academic Press; Bae JM. A practical method for using a survival curve for meta-analysis. J Health Info Stat 2015;40(2):56-64 (Korean). 16. Bae JM. An overview of systematic reviews of diagnostic tests accuracy. Epidemiol Health 2014;36:e Harbord RM, Higgins JPT. Meta-regression in Stata. Stata J 2008; 8(4): Egger M, Smith GD. Bias in location and selection of studies. BMJ 1998;316: Egger M, Smith DG, Schneider M, Minder C. Bias in meta-analysis detected by a simple, graphical test. BMJ 1997;315(7109): Begg CB, Mazumdar M. Operating charateristics of a rank correlation test for publication bias. Biometrics 1994;50: Rothstein HR, Sutton AJ, Borenstein M. Publication bias in metaanalysis: prevention, assessment and adjustments. West Sussex, UK: John Wiley & Sons Ltd.; 2006, p Duval S, Tweedi R. Trim and fill: a simple funnel-plot-based method of testing and adjusting for publication bias in meta-analysis. Biometrics 2000;56: Shim SR, Shin IS, Bae JM. Meta-analysis of diagnostic tests accuracy using STATA software. J Health Info Stat 2015;40(3): (Korean)