Metaanalýzy
Lenka Filova
Meranie veľkosti efektov
Za efekt môžeme považovať akúkoľvej vhodnú sumárnu štatistiku.
Príklad 1: BCG vakcína
Uvažujme 13 štúdií vakcíny Bacillus Calmette-Guerin (BCG), ktoré porovnávali zaočkovaných a nezaočkovaných jedincov s cieľom prevencie tuberkulózy.
Dáta obsahujú číslo experimentu, autora, rok publikovania, informácie o počte vakcinovaných subjektov, ktorí mali pozitívny a negatívny test na tuberkulózu (tpos a tneg), a podobne pre kontrolnú, nevakcinovanú, skupinu (cpos a cneg). Naviac máme k dispozícii geografickú polohu oblasti, kde sa čtúdia vykonávala a spôsob alokovanie pacientov.
library(metafor) ## Warning: package 'metafor' was built under R version 3.6.3## Loading required package: Matrix## Loading 'metafor' package (version 2.4-0). For an overview
## and introduction to the package please type: help(metafor).data(dat.bcg)
str(dat.bcg) ## 'data.frame': 13 obs. of 9 variables:
## $ trial : int 1 2 3 4 5 6 7 8 9 10 ...
## $ author: chr "Aronson" "Ferguson & Simes" "Rosenthal et al" "Hart & Sutherland" ...
## $ year : int 1948 1949 1960 1977 1973 1953 1973 1980 1968 1961 ...
## $ tpos : int 4 6 3 62 33 180 8 505 29 17 ...
## $ tneg : int 119 300 228 13536 5036 1361 2537 87886 7470 1699 ...
## $ cpos : int 11 29 11 248 47 372 10 499 45 65 ...
## $ cneg : int 128 274 209 12619 5761 1079 619 87892 7232 1600 ...
## $ ablat : int 44 55 42 52 13 44 19 13 27 42 ...
## $ alloc : chr "random" "random" "random" "random" ...Označme a,b počty pozitívnych/negatívnych testov vo vakcinovanej skupine a c,d počty pozitívnych/negatívnych testov v kontrolnej skupine.
Vypočítame tzv. log-odds ratio \[ LOR=\log\frac{ad}{bc} \] a varianciu \(V=\frac{1}{a}+\frac{1}{b}+\frac{1}{c}+\frac{1}{d}\).
Y <- with(dat.bcg, log(tpos * cneg/(tneg * cpos))) #log odds ratio
V <- with(dat.bcg, 1/tpos + 1/cneg + 1/tneg + 1/cpos) #variancia
cbind(Y, V)## Y V
## [1,] -0.93869414 0.357124952
## [2,] -1.66619073 0.208132394
## [3,] -1.38629436 0.433413078
## [4,] -1.45644355 0.020314413
## [5,] -0.21914109 0.051951777
## [6,] -0.95812204 0.009905266
## [7,] -1.63377584 0.227009675
## [8,] 0.01202060 0.004006962
## [9,] -0.47174604 0.056977124
## [10,] -1.40121014 0.075421726
## [11,] -0.34084965 0.012525134
## [12,] 0.44663468 0.534162172
## [13,] -0.01734187 0.071635117V balíku metafor máme na výpočet veľkosti efektov k dispozícii funkciu escalc. To, aký efekt chceme počítať, špecifikujeme pomocou argumentu measure:
ES <- escalc(ai = tpos, bi = tneg, ci = cpos, di = cneg,data = dat.bcg,measure = "OR")
ES## trial author year tpos tneg cpos cneg ablat alloc
## 1 1 Aronson 1948 4 119 11 128 44 random
## 2 2 Ferguson & Simes 1949 6 300 29 274 55 random
## 3 3 Rosenthal et al 1960 3 228 11 209 42 random
## 4 4 Hart & Sutherland 1977 62 13536 248 12619 52 random
## 5 5 Frimodt-Moller et al 1973 33 5036 47 5761 13 alternate
## 6 6 Stein & Aronson 1953 180 1361 372 1079 44 alternate
## 7 7 Vandiviere et al 1973 8 2537 10 619 19 random
## 8 8 TPT Madras 1980 505 87886 499 87892 13 random
## 9 9 Coetzee & Berjak 1968 29 7470 45 7232 27 random
## 10 10 Rosenthal et al 1961 17 1699 65 1600 42 systematic
## 11 11 Comstock et al 1974 186 50448 141 27197 18 systematic
## 12 12 Comstock & Webster 1969 5 2493 3 2338 33 systematic
## 13 13 Comstock et al 1976 27 16886 29 17825 33 systematic
## yi vi
## 1 -0.9387 0.3571
## 2 -1.6662 0.2081
## 3 -1.3863 0.4334
## 4 -1.4564 0.0203
## 5 -0.2191 0.0520
## 6 -0.9581 0.0099
## 7 -1.6338 0.2270
## 8 0.0120 0.0040
## 9 -0.4717 0.0570
## 10 -1.4012 0.0754
## 11 -0.3408 0.0125
## 12 0.4466 0.5342
## 13 -0.0173 0.0716cbind(ES$yi, ES$vi)## [,1] [,2]
## [1,] -0.93869414 0.357124952
## [2,] -1.66619073 0.208132394
## [3,] -1.38629436 0.433413078
## [4,] -1.45644355 0.020314413
## [5,] -0.21914109 0.051951777
## [6,] -0.95812204 0.009905266
## [7,] -1.63377584 0.227009675
## [8,] 0.01202060 0.004006962
## [9,] -0.47174604 0.056977124
## [10,] -1.40121014 0.075421726
## [11,] -0.34084965 0.012525134
## [12,] 0.44663468 0.534162172
## [13,] -0.01734187 0.071635117Modely s fixnými a s náhodnými efektmi pomocou príkazu rma
Model s fixnými efektmi špecifikujeme v rma pomocou argumentu method="FE":
result.or <- rma(yi = Y, vi = V, method = "FE") # Log Odds Ratio
summary(result.or)##
## Fixed-Effects Model (k = 13)
##
## logLik deviance AIC BIC AICc
## -76.0290 163.1649 154.0580 154.6229 154.4216
##
## I^2 (total heterogeneity / total variability): 92.65%
## H^2 (total variability / sampling variability): 13.60
##
## Test for Heterogeneity:
## Q(df = 12) = 163.1649, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## -0.4361 0.0423 -10.3190 <.0001 -0.5190 -0.3533 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Model s náhodnými efektmi:
result.or <- rma(yi=Y, vi=V)rma vracia objekt s nasledovnými zložkami:
names(result.or) ## [1] "b" "beta" "se" "zval" "pval"
## [6] "ci.lb" "ci.ub" "vb" "tau2" "se.tau2"
## [11] "tau2.fix" "tau2.f" "k" "k.f" "k.eff"
## [16] "k.all" "p" "p.eff" "parms" "m"
## [21] "QE" "QEp" "QM" "QMp" "I2"
## [26] "H2" "R2" "vt" "int.only" "int.incl"
## [31] "allvipos" "coef.na" "yi" "vi" "X"
## [36] "weights" "yi.f" "vi.f" "X.f" "weights.f"
## [41] "M" "ai.f" "bi.f" "ci.f" "di.f"
## [46] "x1i.f" "x2i.f" "t1i.f" "t2i.f" "ni"
## [51] "ni.f" "ids" "not.na" "subset" "slab"
## [56] "slab.null" "measure" "method" "weighted" "test"
## [61] "dfs" "s2w" "btt" "intercept" "digits"
## [66] "level" "control" "verbose" "add" "to"
## [71] "drop00" "fit.stats" "formula.yi" "formula.mods" "version"
## [76] "model" "call" "time"Najdôležitejšie zložky:
b: výsledný efektci.lb, ci.ub: dolná, resp. horná hranica intervalu spoľahlivostiyi: vektor veľkostí efektovvi: vektor variancií efektov
Príklad 2: Amlodipín
Máme 8 randomizovaných štúdií amlodipine vs placebo. Cieľom bolo zmierniť následky srdcovej angíny.
library(meta)## Warning: package 'meta' was built under R version 3.6.3## Loading 'meta' package (version 4.11-0).
## Type 'help(meta)' for a brief overview.##
## Attaching package: 'meta'## The following objects are masked from 'package:metafor':
##
## baujat, forest, funnel, funnel.default, labbe, radial, trimfilldata(amlodipine)
str(amlodipine) ## 'data.frame': 8 obs. of 7 variables:
## $ study : chr "Protocol 154" "Protocol 156" "Protocol 157" "Protocol 162A" ...
## $ n.amlo : int 46 30 75 12 32 31 27 46
## $ mean.amlo: num 0.232 0.281 0.189 0.093 0.162 ...
## $ var.amlo : num 0.2254 0.1441 0.1981 0.1389 0.0961 ...
## $ n.plac : int 48 26 72 12 34 31 27 47
## $ mean.plac: num -0.0027 0.027 0.0443 0.2277 0.0056 ...
## $ var.plac : num 0.0007 0.1139 0.4972 0.0488 0.0955 ...Za efekt teraz vezmime priemerné rozdiely:
result.md <- rma(m1 = mean.amlo, m2 = mean.plac,sd1 = sqrt(var.amlo), sd2 = sqrt(var.plac),n1 = n.amlo, n2 = n.plac,method = "FE", measure = "MD",data = amlodipine)Vypíšeme výsledný efekt a interval spoľahlivosti:
result.md$b## [,1]
## intrcpt 0.161895c(result.md$ci.lb, result.md$ci.ub)## [1] 0.09860271 0.22518735Akým percentuálnym podielom sa každá štúdia podieľa na výsledku?
contributions <- 1/result.md$vi/sum(1/result.md$vi) * 100
cbind(contributions)## contributions
## [1,] 21.218718
## [2,] 11.354529
## [3,] 10.923051
## [4,] 6.666883
## [5,] 17.942547
## [6,] 10.848040
## [7,] 1.661018
## [8,] 19.385215Najviac výsledky ovplyvňuje prvá štúdia:
max(contributions)## [1] 21.21872amlodipine$study[which(contributions == max(contributions))]## [1] "Protocol 154"contributions <- 1/result.md$vi/sum(1/result.md$vi) * 100
par(mar = c(5, 10, 5, 5))
barplot(contributions, names = amlodipine$study,xlim = c(0, 50), las = 2, horiz = TRUE,col = "blue")
Model s náhodnými efektmi
Pozrime sa teraz bližšie na model s náhodnými efektmi.
Tu potrebujeme zvolit metodu na odhad \(\tau\), napr. REML: restricted maximum likelihood
result.md <- rma(m1 = mean.amlo, m2 = mean.plac,sd1 = sqrt(var.amlo), sd2 = sqrt(var.plac),n1 = n.amlo, n2 = n.plac,method = "REML", measure = "MD",data = amlodipine)
summary(result.md)##
## Random-Effects Model (k = 8; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## 3.3094 -6.6188 -2.6188 -2.7270 0.3812
##
## tau^2 (estimated amount of total heterogeneity): 0.0001 (SE = 0.0042)
## tau (square root of estimated tau^2 value): 0.0116
## I^2 (total heterogeneity / total variability): 1.54%
## H^2 (total variability / sampling variability): 1.02
##
## Test for Heterogeneity:
## Q(df = 7) = 12.3311, p-val = 0.0902
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.1617 0.0326 4.9584 <.0001 0.0978 0.2257 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Vypíšeme odhad, smerodajnú odchýlku a interval spoľahlivosti pre odhad \(\tau\):
result.md$tau2## [1] 0.0001352638result.md$se.tau2## [1] 0.00423431result.md$tau2 + 1.96 * c(-1, 1) * result.md$se.tau2 #95% CI## [1] -0.008163984 0.008434512Rôzne metódy odhadov môžu často dávať veľmi odlišné výsledky. Vykreslime si, ako to vyzerá v našom príklade:
estimators <- c("DL", "REML", "HE", "HS", "SJ", "ML", "EB")
taus <- sapply(estimators, function(method) {
rma(m1 = mean.amlo, m2 = mean.plac, sd1 = sqrt(var.amlo), sd2 = sqrt(var.plac), n1 = n.amlo, n2 = n.plac,
method = method, measure = "MD", data = amlodipine)$tau2})
plot(y = taus, x = 1:length(taus), type = "h", pch = 19, axes = FALSE, xlab = "Estimators")
axis(2, las = 1)
axis(1, at = 1:length(taus), lab = estimators)
Testovanie heterogenity
Na testovanie heterogenity môžeme použiť Q test pre priemerné rozdiely:
MD <- with(amlodipine, mean.amlo - mean.plac)
W <- 1/with(amlodipine, var.amlo/n.amlo + var.plac/n.plac)
Q <- sum(W * (MD - sum(W * MD)/sum(W))^2)
Q## [1] 12.33106S akou pravdepodobnosťou je Q menšie ako 0?
df <- length(MD) - 1
pchisq(Q, df = df, lower = FALSE) ## [1] 0.09018396curve(1 - pchisq(x, df = df), 0, 20)
abline(v = Q, col = "red", lwd = 2)
result.md$QE## [1] 12.33106result.md$QEp## [1] 0.09018396confint(result.md)##
## estimate ci.lb ci.ub
## tau^2 0.0001 0.0000 0.1667
## tau 0.0116 0.0000 0.4082
## I^2(%) 1.5397 0.0000 95.0658
## H^2 1.0156 1.0000 20.2666Heterogenita sa často vizualizuje pomocou forest plotu, čo je graf veľkostí efektov a ich presností. Ide o najčastejší spôsob znázorňovania výsledkov metaanalýzy (dostal sa napíklad aj do loga spoločnosti Cochrane, ktorá sa venuje práve meatanalýzam a systematizáci výsledkov klinických štúdií).
forest(result.md)
Častým problémom pri metaanalýzach je publication bias, teda fakt, že signifikantné výsledky sú publikované s väčšou pravdepodobnosťou. Či je v našom prípade publication bias, vieme vizuálne posúdiť pomocou funnel plotu:
funnel(result.md)
Posúdiť asymetriu z tohto obrázka ale nemusí byť vždy jednoduché. Pomôcť si vieme trim-and-fill metódou, ktorá sa pokúša uhádnuť chýbajúce (nepublikované) štúdie. Výstupom je obázok, v ktorom sú doplnené prázdne bodky symetricky k plným (teda k štúdiám, ktoré máme k dispozícii).
result.rd <- rma(ai = tpos, bi = tneg, ci = cpos, di = cneg,data = dat.bcg,measure = "RD",method = "DL")
trimfill(result.rd) ##
## Estimated number of missing studies on the right side: 4 (SE = 2.4132)
##
## Random-Effects Model (k = 17; tau^2 estimator: DL)
##
## tau^2 (estimated amount of total heterogeneity): 0.0000 (SE = 0.0000)
## tau (square root of estimated tau^2 value): 0.0051
## I^2 (total heterogeneity / total variability): 95.83%
## H^2 (total variability / sampling variability): 23.98
##
## Test for Heterogeneity:
## Q(df = 16) = 383.6062, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## -0.0049 0.0018 -2.7858 0.0053 -0.0084 -0.0015 **
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1funnel(trimfill(result.rd))
Pre príklad s BCG vakcínami vyzerajú tieto dve vizualizácie nasledovne:
forest(result.or)
funnel(result.or)