Werden wir Helden für einen Tag

The anatomy of meta-analysis #2

Posted on Dec 16, 2011 by Chung-hong Chan

上一篇文講了 Fixed effect methods 。其實上一次的 code 有些地方未講完。
Meta-analysis 要處理 heterogeneity 的問題。 Fixed effect methods 假設各個 trials 的 protocol 完全一樣。故此不同 trial 回報的 θ_i 不同的原因，純是因為 chance 。根據此假設， sample size 愈大的 trial 所回報的 θ_i 應該非常接近 bar θ 。
但事實卻並非如此。每個 study 的 protocol 不同之餘， study population 也全然不同。以下是一張我講書時畫過的圖。

各個 trials 之 θ_i 之不同，稱為 heterogeneity 。量
化 heterogeneity 的方法有 Q 。

從 Q 的算法可見，是 θ_i 和 bar θ_IV 的 squared difference 乘以 ω_i ，如果某 jth trial 的 sample size 超大(也即 ω_j 也超大) ，但是 θ_j 和 bar θ_IV 的分別也好大， Q 也會相應增大。 Q 太大，代表 Fixed effect methods 的「sample size 愈大的 trial 所回報的 θ_i 應該非常接近 bar θ 」假設愈不成立。 ((Q 是 follow χ² distribution with df= k-1 ，是可以進行 statistical test 測測它是否高於 0 。但是，此 test 的 power 很低，就算此 test 說 Q not significantly larger than 0 ，也不代表你的 meta-analysis 沒有 heterogeneity 。))
無論你的 meta-analysis 有否 heterogeneity 的問題，我也建議使用 random-effect model 。常用而又計算容易的方法有 DerSimonian-Laird random effect model 。
根據 random effect model ， θ_i 是假定為 normally distributed ，而其 mean 是 bar θ ， variance 是 τ² 。每個 trial 的 weight 會有所調整 (ω'_i )，變成根據 1 除以 within trial 及 between trial 的 variance 之總和。其他的算法與 fixed effect method 無異。

Source code 如下

[code]
# < - Fixed effect and random effect meta-analysis ->
# CH Chan 2011

library(metafor)
stud.label <- c("Auckland", "Block", "Dorah", "Gamsu", "Morrison", "Papageorgiou", "Tauesch")
ev.trt <- c(36, 1, 4, 14, 3, 1, 8)
ev.ctrl <- c(60, 5, 11, 20, 7, 7, 10)
n.trt <- c(532, 69, 81, 131, 67, 71, 56)
n.ctrl <- c(538, 61, 63, 137, 59, 75, 71)

# Contingency table notation
# -------------------
# | | OU+ | OU- |
# -------------------
# | EX+ | A | B |
# -------------------
# | EX- | C | D |
# -------------------

cell.a <- ev.trt
cell.b <- ev.ctrl
cell.c <- n.trt - ev.trt
cell.d <- n.ctrl - ev.ctrl

benchmark <- rma.uni(ai=cell.a, bi=cell.b, ci=cell.c, di=cell.d, measure="OR", method="FE", slab=stud.label)

## from benchmark: summary(benchmark)
## Q = 6.8597, df = 6, p=0.3340
## estimate pooled ES (log OR): -0.6003
## se: 0.1624
## z: -3.6972

## png("forestplot.png")
## forest(benchmark)
## dev.off()

# real deal

oddsratio <- (cell.a * cell.d) / (cell.b * cell.c)
ES <- log(oddsratio)
ES.var <- (1/cell.a) + (1/cell.b) + (1/cell.c) + (1/cell.d)
ES.se <- sqrt(ES.var)
oddsratio.95CI.u <- exp(ES+(qnorm(0.975)*ES.se))
oddsratio.95CI.l <- exp(ES-(qnorm(0.975)*ES.se))

## # Mantel-Haenszel method, and I don't like it. for comment only.
## total.n <- cell.a + cell.b + cell.c + cell.d
## c.OR.MH <- sum((cell.a*cell.d)/total.n) / sum((cell.b*cell.c)/total.n)
## # Just don't ask what the F these quantities are. Ask Kenneth Rothman. (Modern Epidemiology, )
## G <- (cell.a*cell.d)/total.n
## H <- (cell.b*cell.c)/total.n
## P <- (cell.a + cell.d)/total.n
## Q <- (cell.a + cell.d)/total.n
## Term1 <- sum(G*P) / (2*(sum(G))^2)
## Term2 <- (sum((G*Q )+(H*P))) / (sum(G)*sum(H)*2)
## Term3 <- sum(H*Q) / (2*(sum(H))^2)
## OR.MH.var <- Term1+Term2+Term3
## c.OR.MH.95CI.u <- exp(log(c.OR.MH) + (qnorm(0.975)*sqrt(OR.MH.var)))
## c.OR.MH.95CI.l <- exp(log(c.OR.MH) - (qnorm(0.975)*sqrt(OR.MH.var)))

# inverse variance method approximation.
# The "fixed effect" method in RevMan

ES.weight <- 1 / ES.var # inverse variance
sum.weight <- sum(ES.weight)
sum.ESweight <- sum(ES.weight * ES)
c.ES <- sum.ESweight / sum.weight
#actually, we can use the weighted.mean function
# c.ES <- weighted.mean(ES, ES.weight)
c.ES.se <- sqrt(1/sum.weight)
c.ES.95CI.u <- exp(c.ES+qnorm(0.975)*c.ES.se)
c.ES.95CI.l <- exp(c.ES-qnorm(0.975)*c.ES.se)

## useless test statistic for overall effect
## two-sided test
## Overall.p <- 2*pnorm(c.ES/c.ES.se)

Hetero.Q <- sum(ES.weight*((ES-c.ES)^2))
#test for significance of Q
k <- length(ES) #k trials
Hetero.Q.df <- k - 1
Q.pvalue <- pchisq(q=Hetero.Q, df=Hetero.Q.df, lower.tail=FALSE)

## Random effect model

benchmark.DL <- rma.uni(ai=cell.a, bi=cell.b, ci=cell.c, di=cell.d, measure="OR", method="DL", slab=stud.label)

# calculation of tausq

tausq <-(Hetero.Q-Hetero.Q.df)/((sum.weight)-((sum(ES.weight^2)/sum.weight)))
# tau sq is very small.

ES.weight.DL <- 1 / (ES.var+tausq)
c.ES.DL <- (sum(ES*ES.weight.DL))/sum(ES.weight.DL)
c.ES.DL.se <- sqrt(1/sum(ES.weight.DL))
[/code]