0.1 주의 사항 및 강의 목적

메타분석은 체계적 문헌 고찰을 전제로 하고 있습니다. 본 발표에서는 체계적 문헌 고찰 이후 데이터가 모였을 경우 메타분석을 어떻게 시행하는지에 대한 내용입니다. 또한 개념 파악을 중심으로 진행하며, 여러 오류에 대한 내용은 다루지 않습니다. 강의 목적은 향후 CDM 으로 데이터가 산출 되었을 때 어떻게 내용을 병합할지 계획하기 위해 메타분석의 기본 개념을 익히는 데 있습니다.

0.2 메타분석 기본 개념

4 + 6의 평균은 얼마인가요? 당연히 (4+6)/2 = 5 입니다. 여기서 4가 하나, 6이 하나 있기 때문에 나누기 2를 한 것입니다. 만약 4가 9개고 6이 1개이면 어떤가요? (4 * 9 + 6 * 1)/(9 + 1) =4.2 입니다. 즉 갯수가 많은 쪽으로 평균이 바뀌게 됩니다. 당연하지요. 여기서 첫번째 개념은 무언가의 대표값을 나타낼때 (그게 평균이라면)) 갯수가 많은 쪽으로 가중치를 주어서 표현할 수 있다는 것입니다. 비교위험도(Odds Ratio, OR)의 크기가 어떤 A 연구에서는 4이고 어떤 B 연구에서는 6라면 두 연구의 결과를 합하면 어떨까요? (4 +6)/2로 할 수 있을 까요? 이때도 가중치를 생각할 필요가 있습니다. 이번에는 가중치로 각 연구에 참여했던 연구 대상자 수를 생각해 보겠습니다. 만약 숫자적으로 A가 9명으로 B가 1명으로 연구한 것이라면 위의 4+6의 예처럼 4에 가까운 대표값이 나올 것입니다. 즉 대상자가 많은 A 연구 결과에 가까운 대표값이 나오게 됩니다. 이렇게 각 연구의 값을 연구들의 가중치를 고려하여 계산하여 대표값을 산출하는 것이 메타분석이라고 이야기할 수 있습니다. 그리고 가중치로 가장 많이 사용되는 것은 대상자 수 (sample size)입니다. Standard Error (SE) 가 sample size 의 제곱근의 반비례 하는데 따라서 1/SE를 자주 사용합니다.

0.3 실습 데이터 불러오기

나눠 드린 압출파일을 풀고, meta.Rproj를 더블 클릭하십시오. 그러면 아래와 같은 화면이 표시될 것 입니다. R 화면 stomach.csv라는 파일이 보일 것 입니다. 이 파일을 불러오겠습니다.

mm<-read.csv("stomach.csv")
mm <-mm[, 1:9]

데이터의 구조를 보겠습니다.

library(kableExtra)
kable(mm[1:5,]) %>% kable_styling()
Author nation Year jf jt esT esR ll ul
Partanen Finland 1994 ceramic 1 OR 4.91 0.48 49.80
Coggon D UK 1990 ceramic 1 RR 1.50 0.80 3.00
N Aragones sweden 2002 ceramic 1 RR 1.26 0.91 1.74
Miguel Santibanez spain 2011 ceramic 1 OR 0.57 0.08 3.77
Miguel Santibanez spain 2011 construction 2 OR 1.20 0.65 2.22

저자, 나라, 출판년도, 직종, 직종번호, OR|RR, lower limit, upper limit 이 있습니다. 이를 이용해서 meta analysis를 해보겠습니다. 본 연구의 목적은 분진 노출자와 위암 발생률의 메타분석입니다. 그리고 여기서는 OR|RR을 구분하지 않고 RR로 생각해서 분석해 보겠습니다.

0.4 Data step

필요한 packages 를 불러옵니다.

library(dslabs)
library(tidyverse)
library(metafor)

위 데이터에서 SE 를 구해서 가중치로 사용할 것 입니다. Upper limit 은 OR + 1.96 * SE 이니까, 이를 이용해서 SE를 구합니다. 다만 Exp(Beta) = OR 이므로 log를 이용해 SE를 구합니다. log(UL) =log(OR)+1.96 *SE --> SE = (log(UL) - log (OR))/1.96 (mm데이터에 esT가 factor이네요 numberic으로 바꾸겠습니다.)

mm <- mm %>% 
  mutate(lnesR=log(esR)) %>%
  mutate(se=(log(ul)-log(esR))/1.96)

0.5 기초 분석

metafor packages에서 effect Size (OR|RR)과 SE 그리고 데이터를 지정하고, fixed effect|random effect를 고려합니다. fixed effect 는 가중치를 1/se 로 random effect는 1/(se^2 + v) 를 이요합니다. 간단히 설명하면 이질성이 크면 random effect로 가중치를 v 만큰 보정해 준다는 것입니다. 아래 결과에서 이질성 검정의 P value를 보고 판단하면 됩니다.[참고: df=31일때 Chisq 값은 44.9입니다. 그런데 Q 값이 85로 크게 나왔네요 이러면 이질성이 있는 것입니다.]

res <-rma(yi=lnesR, sei=se, measure="RR", digits=1, data=mm)
summary(res)
## 
## Random-Effects Model (k = 32; tau^2 estimator: REML)
## 
##   logLik  deviance       AIC       BIC      AICc 
##    -10.8      21.5      25.5      28.4      26.0   
## 
## tau^2 (estimated amount of total heterogeneity): 0.0 (SE = 0.0)
## tau (square root of estimated tau^2 value):      0.1
## I^2 (total heterogeneity / total variability):   59.77%
## H^2 (total variability / sampling variability):  2.49
## 
## Test for Heterogeneity:
## Q(df = 31) = 85.3, p-val < .1
## 
## Model Results:
## 
## estimate   se  zval  pval  ci.lb  ci.ub 
##      0.3  0.0   5.7   <.1    0.2    0.4  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

각 스터디의 요약값과 가중치를 고려해서 합한 값을 그림으로 봅니다. Randome Effect model로 1.31 (1.19-1.44) 로 나왔네요.

forest(res, digits=2, atransf=exp)

0.6 오류 발견 및 수정

퍼블리케이션 바이어스는 출판 논문의 대칭성을 봅니다. 우선 아래 그림을 보면

funnel(res)

이를 다시 그려보겠습니다. 이질성이 증가했습니다. 그리고 빈 동그라미가 생겼습니다. 만약 대칭성이 있었다면 빈 동그라미에 해댱하는 논문이 있었어야 하는데 없다는 뜻입니다.

rtf <-trimfill(res)
rtf
## 
## Estimated number of missing studies on the left side: 4 (SE = 3.7252)
## 
## Random-Effects Model (k = 36; tau^2 estimator: REML)
## 
## tau^2 (estimated amount of total heterogeneity): 0.0214 (SE = 0.0141)
## tau (square root of estimated tau^2 value):      0.1463
## I^2 (total heterogeneity / total variability):   56.62%
## H^2 (total variability / sampling variability):  2.31
## 
## Test for Heterogeneity:
## Q(df = 35) = 89.4341, p-val < .0001
## 
## Model Results:
## 
## estimate      se    zval    pval   ci.lb   ci.ub 
##   0.2517  0.0467  5.3916  <.0001  0.1602  0.3432  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
funnel(rtf)

통계적으로 대칭적인지도 살펴 봅니다. 크게 이상이 없네요.

regtest(res, model="lm")
## 
## Regression Test for Funnel Plot Asymmetry
## 
## model:     weighted regression with multiplicative dispersion
## predictor: standard error
## 
## test for funnel plot asymmetry: t = 2.7, df = 30, p < .1

0.6.1 직종 보정해 보기

이질성이 있는 이유를 생각해보면 직종이 달라서 일 수 있습니다.

res2 <-rma(yi=lnesR, sei=se, measure="RR",mods=cbind(jt),  data=mm)
res2
## 
## Mixed-Effects Model (k = 32; tau^2 estimator: REML)
## 
## tau^2 (estimated amount of residual heterogeneity):     0.0285 (SE = 0.0180)
## tau (square root of estimated tau^2 value):             0.1688
## I^2 (residual heterogeneity / unaccounted variability): 57.51%
## H^2 (unaccounted variability / sampling variability):   2.35
## R^2 (amount of heterogeneity accounted for):            0.00%
## 
## Test for Residual Heterogeneity:
## QE(df = 30) = 66.3071, p-val = 0.0001
## 
## Test of Moderators (coefficient 2):
## QM(df = 1) = 1.3130, p-val = 0.2519
## 
## Model Results:
## 
##          estimate      se    zval    pval    ci.lb   ci.ub 
## intrcpt    0.1309  0.1346  0.9725  0.3308  -0.1329  0.3947    
## jt         0.0444  0.0388  1.1458  0.2519  -0.0316  0.1204    
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
forest(res2, digits=2, atransf=exp)

각 직종별로 보면 어떤 직종이 가장 위험도가 높은가요? 일반적으로 높은 농도에 노출되는 직종에서 위암 위험도가 높다고 할 수 있을까요?

tt<-predict(res2, newmods=cbind(seq(from=1, to=5, by=1)), transf=exp, addx=TRUE)
tt %>% kable() %>% kable_styling()
pred ci.lb ci.ub cr.lb cr.ub X.intrcpt X.jt
1.191647 0.9798748 1.449189 0.8113796 1.750134 1 1
1.245794 1.0872472 1.427460 0.8711327 1.781591 1 2
1.302401 1.1765435 1.441721 0.9213899 1.840966 1 3
1.361579 1.2112803 1.530528 0.9586333 1.933898 1 4
1.423447 1.2020145 1.685672 0.9817231 2.063925 1 5

이를 그림으로 그려 보겠습니다.

preds <- predict(res2, newmods = cbind(1:5), transf = exp)
jtp<-c(1:5)
pp<-cbind(jtp, data.frame(preds))
wi<-1/(mm$se^2)
size <- 5 + 10 * (wi - min(wi))/(max(wi) - min(wi))

plot(mm$jt, exp(mm$lnesR), pch = 19 ,col = alpha(blues9, 0.4), cex = size/2, 
     xlab = "Occupation", ylab = "Relative Risk",
     las = 1, bty = "l", log = "y", xlim=c(1,5.5))
lines(1:5, preds$pred, lwd=2, col='grey')
lines(1:5, preds$ci.lb, lty = "dashed")
lines(1:5, preds$ci.ub, lty = "dashed")
abline(h = 1, lty = "dotted")