########################################################## # # # 18.6 SECOND-HAND SMOKING # # # # DATE: 08/24/2008 # # # ########################################################## 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 the results of 19 studies # where effect estimates plus 95% confidence intervals are known # rr = estimated relative risk # low.rr = lower bound of the 95% CI on rr # upp.rr = upper bound of the 95% CI on rr study <- 1:19 rr <- c(1.52, 1.52, 0.81, 0.75, 2.07, 1.19, 1.31, 2.16, 2.34, 2.55, 0.79, 1.55, 1.65, 2.01, 1.03, 1.28, 1.26, 2.13, 1.41) low.rr <- c(0.88, 0.39, 0.34, 0.43, 0.82, 0.82, 0.87, 1.08, 0.81, 0.74, 0.25, 0.90, 1.16, 1.09, 0.41, 0.76, 0.57, 1.19, 0.54) upp.rr <- c(2.63, 5.99, 1.90, 1.30, 5.25, 1.73, 1.98, 4.29, 6.75, 8.78, 2.45, 2.67, 2.35, 3.71, 2.55, 2.15, 2.82, 3.83, 3.67) # combining the results on the log scale # first we have to extract the standard error according Section 17.1 log.rr <- log(rr) se.logrr <- 0.5 * (log(upp.rr) - log(low.rr)) / qnorm(0.975) # Generic inverse variance method using R package meta library(meta) # Combining d's comb.logrr <- metagen(log.rr, se.logrr) comb.logrr improved.ci(log.rr, se.logrr, 0.05) # Confidence interval plot with fixed and random effects meta-analysis plot(comb.logrr, comb.f = TRUE, comb.r = TRUE, main = "Second-hand smoking", xlab = "Log relative risk") # Funnel plot and tests for publication bias # method = rank (rank correlation) # method = linear (linear regression) funnel(comb.logrr, main = "Second-hand smoking", xlab = "Log relative risk") metabias(comb.logrr, method = "rank") metabias(comb.logrr, method = "linreg")