########################################################## # # # 18.7 EFFECTIVENESS of MISOPROTOL IN PREVENTING # # GASTROINTESTINAL DAMAGE # # # # 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 ordinal data # Log odds ratios and Agresti's alpha based on Table 18.7 were calculated in a separate program # log.theta = log odds ratio from proportional odds model # var.log.theta = variance of log.theta # log.alpha = log of Agresti's alpha # var.log.alpha = variance of log.alpha log.theta <- c(3.5460, 4.0472, 1.9114, 3.7506, 6.5109, 1.1756, 1.1928, 1.8400, 2.9646, 2.4872, 2.5670, 0.6471, 1.1119) var.log.theta <- c(0.4344, 0.5161, 0.2819, 0.4756, 5.1947, 0.1559, 0.1518, 1.1483, 1.0743, 0.5646, 0.6358, 0.1147, 0.5035) log.alpha <- c(3.0419, 3.5750, 1.6942, 3.0368, 4.0244, 1.1172, 1.1870, 1.8369, 2.9647, 2.4872, 2.3726, 0.5977, 1.0354) var.log.alpha <- c(0.3113, 0.3871, 0.2152, 0.3473, 3.8343, 0.1418, 0.1488, 1.1338, 1.0743, 0.5646, 0.6128, 0.0969, 0.4216) # Generic inverse variance method using R package meta library(meta) # Combining log odds ratios comb.theta <- metagen(log.theta, sqrt(var.log.theta)) comb.theta improved.ci(log.theta, sqrt(var.log.theta), 0.05) # Confidence interval plot with fixed and random effects meta-analysis plot(comb.theta, comb.f = TRUE, comb.r = TRUE, main = "Misoprostol studies", xlab = "Log odds ratio", xlim = c(-1, 4)) # Funnel plot and tests for publication bias # method = rank (rank correlation) # method = linear (linear regression) funnel(comb.theta, main = "Misoprostol studies", xlab = "Log odds ratio") metabias(comb.theta, method = "rank") metabias(comb.theta, method = "linreg") # Combining log Agresti's alpha comb.alpha <- metagen(log.alpha, sqrt(var.log.alpha)) comb.alpha improved.ci(log.alpha, sqrt(var.log.alpha), 0.05) # Confidence interval plot with fixed and random effects meta-analysis plot(comb.alpha, comb.f = TRUE, comb.r = TRUE, main = "Misoprostol studies", xlab = "Log Agrestis's alpha", xlim = c(-1, 4)) # Funnel plot and tests for publication bias # method = rank (rank correlation) # method = linear (linear regression) funnel(comb.alpha, main = "Misoprostol studies", xlab = "Log Agresti's alpha") metabias(comb.alpha, method = "rank") metabias(comb.alpha, method = "linreg")