# Werden wir Helden für einen Tag

Home | About | Archive

# The anatomy of meta-analysis #1

Posted on Dec 15, 2011 by Chung-hong Chan

[code]
# < - Fixed 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
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.df <- sum(ES.weight*((ES-c.ES)^2))
#test for significance of Q
Hetero.Q.df <- length(ES) - 1
Q.pvalue <- pchisq(q=Hetero.Q, df=Hetero.Q.df, lower.tail=FALSE)
[/code]

Inverse variance 就是將 weight 設定為 1 除 effect size 的 variance 。故此， sample size 愈大的 trial ， variance 愈細，故此 weight 愈大。本例計算的 effect size 是 log OR ，而 log OR 的 Variance 如上圖算法。

Source code 是 self-explanatory 的。算出來每個值都與 metafor 的 benchmark case 相同。
inverse variance 屬於 fixed-effect methods ，下一篇文會講講怎樣用 DerSimonian-Laird random effect models 。