E.3 Reproducing the Output

이제 위의 결과가 어떻게 나왔는지 중간 과정을 재현해 보면 다음과 같다.

E.3.1 Relative Risk and Its Confidence Limit

먼저 RR과 그것의 신뢰구간을 구하는 함수는 8장의 함수를 그대로 사용한다. 이제 RR 함수로 자료를 일차적으로 분석(각 study의 RR과 SE, 신뢰구간)하면 다음과 같다.

r9 = sasLM::RR(e.t, n.t, e.c, n.c) ; r9

      p1    p2    RR SElog lower upper
1 0.0797 0.107 0.742 0.179 0.522  1.05
2 0.0580 0.083 0.699 0.189 0.483  1.01
3 0.0852 0.104 0.822 0.249 0.505  1.34
4 0.1226 0.148 0.827 0.124 0.649  1.05
5 0.1049 0.128 0.819 0.165 0.593  1.13
6 0.1085 0.097 1.118 0.088 0.941  1.33

RR들을 병합하려면 log변환하여야 한다. (SE도 log 변환 상태에서의 SE이다.) 그것들을 thi라 하면 다음과 같다.

thi = log(r9$RR) ; thi

[1] -0.298 -0.358 -0.195 -0.190 -0.199  0.112

E.3.2 Fixed Effect Model

Mantel-Haenszel method에 의한 fixed effect들의 weight는 다음과 같이 계산된다.

wi0 = e.c*n.t/(n.t + n.c)    # Mantel-Haenszel method
pwi = wi0/sum(wi0)*100 ; pwi # percent w_i. compare with the metabin output

[1] 11.55 11.02  5.63 21.65 12.03 38.12

하지만, 우리는 다음과 같이 inverse variance로 weight를 줄 것이다.

Vi = r9$SElog^2
wi = 1/Vi ; wi              # inverse variance method

[1]  31.1  28.0  16.2  65.1  36.6 129.1

다음은 계산의 효율성을 위해 sum(wi)를 미리 계산한다.

sumwi = sum(wi) ; sumwi

[1] 306

이제 fixed effect model로 병합(합동)추정한 RR의 점추정치(PE)와 표준오차(SE)는 다음과 같다.

th.hat = sum(wi*thi)/sumwi ; th.hat # TE.fixed in str(r6)

[1] -0.0904

seth.hat = sqrt(1/sumwi) ; seth.hat # seTE.fixed in str(r6)

[1] 0.0572

이것을 다음과 같이 원척도(original scale)로 변환하고 신뢰구간으로 표시해야 알아보기가 쉽다.

exp(th.hat)

[1] 0.914

exp(th.hat + c(-1, 1)*1.96*seth.hat) 

[1] 0.817 1.022

Inverse variance method를 썼기 때문에 metabin과는 약간 차이가 있다. 하지만, 이 방법도 잘못된 것은 아니다.

다음은 연구(study, 시험)간의 이질성(heterogeneity)을 보기위한 Q를 계산하면 다음과 같다.

Q = sum(wi*(thi - th.hat)^2) ; Q # 9.88

[1] 9.88

이것의 자유도와 카이제곱분포를 이용하여 다음과 같이 p-value를 구할 수 있다.

k = length(e.t) ; k  # 6

[1] 6

1 - pchisq(Q, k - 1) # p-value

[1] 0.0786

E.3.3 Random Effect Model

Inter-study variability (\(\tau^2\))를 추정하는 방법은 여러가지 있는데, metabin함수는 기본적(default)으로 REML방법을 사용한다.

우선 가장 단순한 method of moment로 구한 \(\tau^2\) 에 대한 추정값은 다음과 같다.

(Q - (k - 1))/(sumwi - sum(wi^2)/sumwi)  # method of moment, not by REML nor ML

[1] 0.0215

여기에는 위에서 구한 Q가 사용되었다.

metabin의 결과인 r6에서 tau2를 찾아보면 0.021로 비슷하다. 이제 이것으로 계속 계산하면 다음과 같다.

tau2 = r6$tau2 # tau2 by REML

Random effect method에서 각 study의 weight는 summary(r6)의 첫번째 부분 마지막 열에 나와 있다. 이것을 수작업으로 구해보면 다음과 같다. sum(wsi)는 계산의 효율성을 위해 sumwsi로 미리 계산해 두었다.

wsi = 1/(Vi + tau2) ; wsi

[1] 18.9 17.7 12.1 27.6 20.7 34.9

sumwsi = sum(wsi)
pwsi = wsi/sumwsi*100 ; pwsi # percent w star i

[1] 14.31 13.39  9.16 20.92 15.74 26.48

이를 이용하여 구한 log scale에서의 점추정치와 표준오차는 다음과 같다.

th.hat.ran = sum(wsi*thi)/sumwsi ; th.hat.ran # TE.ran in str(r6)

[1] -0.15

seth.hat.ran = sqrt(1/sumwsi) ; seth.hat.ran  # seTE.ran in str(r6)

[1] 0.0871

원래 scale에서의 RR과 그 신뢰구간은 다음과 같다.

exp(th.hat.ran)

[1] 0.861

exp(th.hat.ran + c(-1, 1)*1.96*seth.hat.ran)

[1] 0.726 1.021

이로써 meta-anlaysis의 한가지 예로서 그 계산 과정을 살펴 보았다.

다음은 위의 과정을 하나의 함수로 만든 예시이다. Risk Difference와 Odds Ratio의 경우에도 제시하였다.

require(sasLM)
metaRR = function(e.t, n.t, e.c, n.c, conf.level=0.95)
{
  r1 = RR(e.t, n.t, e.c, n.c, conf.level=conf.level)
  thi = log(r1$RR)
  wi = 1/(r1$SElog)^2
  wi0 = e.c*n.t/(n.t + n.c)
  pwi = wi0/sum(wi0)*100
  r1$pwi = pwi
  sumwi = sum(wi) ; sumwi

  th.hat = sum(wi*thi)/sumwi
  seth.hat = sqrt(1/sumwi)
  z.crit = qnorm(0.5 + conf.level/2)
  ci.lower = exp(th.hat - z.crit*seth.hat) 
  ci.upper = exp(th.hat + z.crit*seth.hat) 
  r2 = data.frame(RR = exp(th.hat), CI.lower = ci.lower, CI.upper = ci.upper)

  Q = sum(wi*(thi - th.hat)^2)
  k = length(e.t)
  pQ = 1 - pchisq(Q, k - 1)
  r3 = data.frame(Q=Q, prob=pQ)

  tau2 = (Q - (k - 1))/(sumwi - sum(wi^2)/sumwi) # method of moment
  wsi = 1/(1/wi + tau2)
  sumwsi = sum(wsi)
  pwsi = wsi/sumwsi*100
  r1$pwsi = pwsi

  th.hat.ran = sum(wsi*thi)/sumwsi
  se2 = sqrt(1/sumwsi)
  ci.lower2 = exp(th.hat.ran - z.crit*se2)
  ci.upper2 = exp(th.hat.ran + z.crit*se2)
  r4 = data.frame(RR = exp(th.hat.ran), CI.lower = ci.lower2, CI.upper = ci.upper2)
  Res = list(RelRisk = r1, Heterogeneity = r3, tau2 = tau2, Fixed = r2, Random = r4)
  return(Res)
}
metaRR(e.t, n.t, e.c, n.c)

$RelRisk
      p1    p2    RR SElog lower upper   pwi  pwsi
1 0.0797 0.107 0.742 0.179 0.522  1.05 11.55 14.36
2 0.0580 0.083 0.699 0.189 0.483  1.01 11.02 13.45
3 0.0852 0.104 0.822 0.249 0.505  1.34  5.63  9.23
4 0.1226 0.148 0.827 0.124 0.649  1.05 21.65 20.88
5 0.1049 0.128 0.819 0.165 0.593  1.13 12.03 15.77
6 0.1085 0.097 1.118 0.088 0.941  1.33 38.12 26.31

$Heterogeneity
     Q   prob
1 9.88 0.0786

$tau2
[1] 0.0215

$Fixed
     RR CI.lower CI.upper
1 0.914    0.817     1.02

$Random
    RR CI.lower CI.upper
1 0.86    0.724     1.02

metaRD = function(e.t, n.t, e.c, n.c, conf.level=0.95)
{
  r1 = RD(e.t, n.t, e.c, n.c, conf.level=conf.level)
  wi = 1/(r1$SE)^2
  sumwi = sum(wi)
  pe = sum(wi*r1$RD)/sumwi
  se = sqrt(1/sumwi)
  z.crit = qnorm(0.5 + conf.level/2)
  ci.lower = pe - z.crit*se
  ci.upper = pe + z.crit*se
  r2 = data.frame(PE=pe, SE=se, CI.lower=ci.lower, CI.upper=ci.upper)
  Q = sum(wi*(r1$RD - pe)^2)
  
  k = length(wi)
  pQ = 1 - pchisq(Q, k - 1)
  r3 = data.frame(Q=Q, prob=pQ)
  
  tau2 = max(0, (Q - (k - 1))/(sumwi - sum(wi^2)/sumwi))
  wsi = 1/(1/wi + tau2)
  sumwsi = sum(wsi)
  pe2 = sum(wsi*r1$RD)/sumwsi
  se2 = sqrt(1/sumwsi)
  ci.lower2 = pe2 - z.crit*se2
  ci.upper2 = pe2 + z.crit*se2
  r4 = data.frame(PE=pe2, SE=se2, CI.lower=ci.lower2, CI.upper=ci.upper2)
  
  Res = list(RiskDiff = r1, Heterogeneity = r3, tau2 = tau2, Fixed = r2, Random = r4)
  return(Res)
}
metaRD(e.t, n.t, e.c, n.c)

$RiskDiff
      p1    p2      RD      SE    lower    upper
1 0.0797 0.107 -0.0277 0.01652 -0.06007 0.004677
2 0.0580 0.083 -0.0250 0.01307 -0.05058 0.000658
3 0.0852 0.104 -0.0184 0.02337 -0.06419 0.027421
4 0.1226 0.148 -0.0256 0.01667 -0.05831 0.007030
5 0.1049 0.128 -0.0231 0.01977 -0.06190 0.015616
6 0.1085 0.097  0.0115 0.00903 -0.00621 0.029175

$Heterogeneity
     Q   prob
1 9.57 0.0884

$tau2
[1] 0.000204

$Fixed
       PE     SE CI.lower CI.upper
1 -0.0098 0.0058  -0.0212  0.00156

$Random
       PE      SE CI.lower CI.upper
1 -0.0151 0.00866   -0.032  0.00191

metaOR = function(e.t, n.t, e.c, n.c, method="MH", conf.level=0.95)
{
  Res1 = OR(e.t, n.t, e.c, n.c, conf.level)

  ai = e.t
  bi = n.t - e.t
  ci = e.c
  di = n.c - e.c
  ni = ai + bi + ci + di
   
  if (method == "MH") {
    pe = sum(ai*di/ni)/sum(bi*ci/ni)
    Ri = ai*di/ni
    Si = bi*ci/ni
    Pi = (ai + di)/ni
    Qi = (bi + ci)/ni
    Rsum = sum(Ri)
    Ssum = sum(Si)
    selog = sqrt(sum(Pi*Ri/(2*Rsum^2) + (Pi*Si + Qi*Ri)/(2*Rsum*Ssum) + Qi*Si/(2*Ssum^2)))
  } else {
    method = "Inverse Variance"
    thi = log(Res1$OR)
    vi = Res1$SElog^2
    wi = 1/vi
    sumwi = sum(wi)
    pe = exp(sum(wi*thi)/sumwi)
    selog = sqrt(1/sumwi)
  }
  
  z.crit = qnorm(0.5 + conf.level/2)
  ci.lower = exp(log(pe) - z.crit*selog)
  ci.upper = exp(log(pe) + z.crit*selog)
  Res2 = data.frame(OR=pe, SElog=selog, CI.lower=ci.lower, CI.upper=ci.upper)
  return(list(OR=Res1, Method=method, Fixed=Res2))
}
metaOR(e.t, n.t, e.c, n.c)

$OR
    odd1   odd2    OR  SElog lower upper
1 0.0866 0.1203 0.720 0.1972 0.489  1.06
2 0.0616 0.0905 0.681 0.2029 0.457  1.01
3 0.0931 0.1155 0.806 0.2745 0.471  1.38
4 0.1397 0.1740 0.803 0.1431 0.606  1.06
5 0.1172 0.1469 0.798 0.1876 0.553  1.15
6 0.1217 0.1075 1.133 0.0981 0.935  1.37

$Method
[1] "MH"

$Fixed
     OR  SElog CI.lower CI.upper
1 0.903 0.0637    0.797     1.02

단순한 2x2 표로 나타나는 자료의 경우에, meta-analysis나 strata로 구분된 경우의 분석법이 개념적으로는 동일하다. 따라서, 2x2 표로 결과가 제시된 여러 연구를 병합분석할 때는 sasLM의 RDmn, RRmn, ORmn, ORcmh 등의 함수를 활용할 수 있다.