########################################################## # # # 18.1 VALIDITY STUDIES # # # # 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 correlation coefficents from 20 studies # n = sample size # r = product moment correlation coefficent # std.r = standard deviation of r study <- 1:20 n <- c( 10, 20, 13, 22, 28, 12, 12, 36, 19, 12, 36, 75, 33, 121, 37, 14, 40, 16, 14, 20) r <- c(0.68, 0.56, 0.23, 0.64, 0.49, -0.04, 0.49, 0.33, 0.58, 0.18, -0.11, 0.27, 0.26, 0.40, 0.49, 0.51, 0.40, 0.34, 0.42, 0.16) std.r <- sqrt((1 - r^2)^2 / (n - 1)) # Fisher's z-transformation and standard devation z <- 0.5 * log ((1 + r) / (1 - r)) std.z <- sqrt(1 / (n - 3)) # Generic inverse variance method using R package meta library(meta) # Combining r's comb.r <- metagen(r, std.r) comb.r improved.ci(r, std.r, 0.05) # Confidence interval plot with fixed and random effects meta-analysis plot(comb.r, comb.f = TRUE, comb.r = TRUE, main = "Validity studies", xlab = "Correlation coefficient") # Funnel plot and tests for publication bias # method = rank (rank correlation) # method = linear (linear regression) funnel(comb.r, main = "Validity studies", xlab = "Correlation coefficient") metabias(comb.r, method = "rank") metabias(comb.r, method = "linreg") # Combining z's comb.z <- metagen(z, std.z) comb.z improved.ci(z, std.z, 0.05) # Confidence interval plot with fixed and random effects meta-analysis plot(comb.z, comb.f = TRUE, comb.r = TRUE, main = "Validity studies", xlab = "Fisher's z-transform") # Funnel plot and tests for publication bias # method = rank (rank correlation) # method = linear (linear regression) funnel(comb.z, main = "Validity studies", xlab = "Fisher's z-transform") metabias(comb.z, method = "rank") metabias(comb.z, method = "linreg")