最近有想為碩士班「贖罪」的感覺,一直在想,到底碩士班真的欠缺了甚麼?
為甚麼給了十二萬教出來的人都好像沒有甚麼專業知識似的。到底是我自己的問題,還是人家教育的問題。
我假定是我自己的問題。應該是我水過鴨背,也沒有自行修讀精進。別人可能講得太少,但並不代表我也要知得太少,我有責任去滿足自己的好奇心,有責任去問,有責任去讀。或者就是因為我欠缺這種責任心,人家才覺得我不是讀博士材料。 ((講笑的,我因為甚麼原因落選,我是心照的。)) ((讀碩士時,課程視為仙水級奇書的 Kenneth Rothman Modern Epidemiology 我一頁也沒有讀過。))
我的 thesis 是幹 meta-analysis 和 meta-regression ,但是我到現在都沒有試過自行做一次 meta-analysis ,只是用人家寫好的 routine 。 ((以前的碩士班,是教用 Revman 軟件,是我自己教自己用 R 來做 meta-analysis 的。)) 以前不覺得這是甚麼問題,正如你用 iPhone ,也未必需要知道 iPhone 內部怎樣運作。
但是,如果你是以成為研究人員為目的的話,不知道物件的內部運作,根本就與用家無分別。 ((人家寫好的 routine 當然是好使好用,但同時是一層的 abstraction 。)) 怎樣話比人知你是這方面的研究人員? ((當然,就算流行病學的專家,實在是沒有必要打爛沙盆問到篤各種算法如何,只需要知道何時用那種算法就可,講到尾流行病學也是應用科學。但如果你發現某算法限制太大,未必適合實制用途,想發明新的算法,那就要理解不同算法的實制步驟。))
我決定像個小孩那樣,將物品拆開去理解其運作方式。
最近才比較理解 Linear regression 、 Logistic regression 、 PCA 的實際運作方法。今次是想用教自己由零開始做一次 meta-analysis 的方法。就如上次 PCA 一樣,這是很無聊但有教育意義的活動。

有數據如下 ((Data from randomised trials before 1980 of corticosteroid therapy in premature labour and its effect on neonatal death.))

用 metafor 畫出 forest plot 會是這樣的 ((再看看此圖)) :

以下是 source code.

[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]

我是參照 Egger 本書 ((Egger M, Smith GD, Altman DG eds. Systematic reviews in health care: meta-analysis in context. 2nd Edition. BMJ books.)) Chapter 15 的 equations 來寫 code 。

要先了解何謂 meta-analysis 。 Meta-analysis 說穿了是 weighted average of effect size 。要計算 combined effect size bar θIV 只要有每個 trial 的 Effect size θi 和 weight ωi 就可。問題是 θi 和 ωi 如何計出來。
以 binary outcomes 來說,最常用的 Effect size 有 Odds ratio (OR) 、 Rate ratio (RR) 和 Rate Difference (RD) 等等。注意的是由於 OR 的分佈並非常態,在進行 Meta-analysis 時會用 log(OR) ,整合過後才用 antilog (exp) 將它變回 OR 。
而 weight ,有不同的計算方法,常用的方法有 Inverse variance method, Mantel-Haenszal method, Peto's odds ratio method 和 DerSimonian-Laird random effect models 等等。我們先由最簡單的 Inverse Variance 入手。
上圖的 bar θ 計算法是 general 的,無論用那種 effect size 或 weight 都是這樣計。但是其他的算式是適用於今次的例子。
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 。