########################################################## # # # 18.8 PREVENTION OF TUBERCULOSIS # # # # DATE: 08/24/2008 # # # ########################################################## # simple function for computing the improved confidence interval # in the generic intervse variance method according to # Hartung and Knapp (2001a,b) using the DerSimonian-Laird # estimator of the heterogeneity parameter # theta = vector of study-specific estimates # xi = vector of study-specific standard deviations # alpha = 1 - confidence coefficients improved.ci <- function(theta, xi, alpha) { k <- length(theta) # number of studies v <- 1 / xi^2 # fixed effects weights theta.fix <- sum(v * theta) / sum(v) # fixed effects estimator of theta cochran.q <- sum(v * (theta - theta.fix)^2) # Cochran's Q tau.dsl <- (cochran.q - k + 1) / (sum(v) - sum(v^2) / sum(v)) # DerSimonian-Laird (DSL) estimator tau.est <- max(0, tau.dsl) # truncated DSL estimator w <- 1 / (xi^2 + tau.est) # random effects weights theta.ran <- sum(w * theta) / sum(w) # random effects estimator of theta hartung.q <- sum(w * (theta - theta.ran)^2) / ((k - 1) * sum(w)) # improved variance estimator ci.lower <- theta.ran - sqrt(hartung.q) * qt(1 - alpha/2, k - 1) # lower bound of CI on theta ci.upper <- theta.ran + sqrt(hartung.q) * qt(1 - alpha/2, k - 1) # upper bound of CI on theta test <- theta.ran / sqrt(hartung.q) # test statistic p.value <- 2 * (1 - pt(abs(test), k - 1)) # p-value return(cbind(tau.est, theta.ran, ci.lower, ci.upper, p.value)) } # Combining binary data # vac.dis = number of diseases in the vaccinated group # vac.nodis = number of no diseases in the vaccinated group # novac.dis = number of diseases in the non-vaccinated group # novac.nodis = number of no diseases in the non-vaccinated group # vac = number of indivduals in the vaccinated group # novac = number of individuals in the non-vaccinated group vac.dis <- c( 4, 6, 3, 62, 33, 180, 8, 505, 29, 17, 186, 5, 27) vac.nodis <- c(119, 300, 228, 13536, 5036, 1361, 2537, 87886, 7470, 1699, 50448, 2493, 16886) novac.dis <- c( 11, 29, 11, 248, 47, 372, 10, 499 , 45, 65, 414, 3, 29) novac.nodis <- c(128, 274, 209, 12619, 5761, 1079, 619, 87892, 7232, 1600, 27197, 2338, 17825) vac <- vac.dis + vac.nodis novac <- novac.dis + novac.nodis ############################################################## # # # Effect size: log relative risk and standard error # # # ############################################################## p.vac <- vac.dis / vac p.novac <- novac.dis / novac log.rr <- log(p.vac / p.novac) std.logrr <- sqrt(1 / vac.dis - 1 / vac + 1 / novac.dis - 1 / novac) # Generic inverse variance method using R package meta library(meta) # Combining log relative risks using functions metagen and metabin # Note: results are backtransformed to the relative risk scale comb.logrr <- metagen(log.rr, std.logrr) comb.logrr metabin(vac.dis, vac, novac.dis, novac, method = "Inverse", sm = "RR") improved.ci(log.rr, std.logrr, 0.05) # Confidence interval plot with fixed and random effects meta-analysis plot(comb.logrr, comb.f = TRUE, comb.r = TRUE, main = "Tuberculosis studies", xlab = "Log relative risk", xlim = c(-2.5,1)) # Funnel plot and tests for publication bias # method = rank (rank correlation) # method = linear (linear regression) funnel(comb.logrr, main = "Tuberculosis studies", xlab = "Log relative risk") metabias(comb.logrr, method = "rank") metabias(comb.logrr, method = "linreg") # B E Y O N D T H E B O O K # Mantel-Haenszel weights in the fixed effects model of meta-analyis # Analysis in random effects model is generic inverse variance method metabin(vac.dis, vac, novac.dis, novac, method = "MH", sm = "RR")