########################################################## # # # 18.3 DENTIFRICE DATA # # # # 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 (standardized) mean differences from 9 studies # n.NaF = sample size in the NaF group # mean.NaF = mean in the NaF group # std.NaF = standard deviation in the NaF group # n.SMFP = sample size in the SMFP group # mean.SMFP = mean in the SMFP group # std.SMFP = standard deviation in the SMFP group study <- 1:9 n.NaF <- c( 134, 175, 137, 184, 174, 754, 209, 1151, 679) mean.NaF <- c(5.96, 4.74, 2.04, 2.70, 6.09, 4.72, 10.10, 2.82, 3.88) std.NaF <- c(4.24, 4.64, 2.59, 2.32, 4.86, 5.33, 8.10, 3.05, 4.85) n.SMFP <- c( 113, 151, 140, 179, 169, 736, 209, 1122, 673) mean.SMFP <- c(6.82, 5.07, 2.51, 3.20, 5.81, 4.76, 10.90, 3.01, 4.37) std.SMFP <- c(4.72, 5.38, 3.22, 2.46, 5.14, 5.29, 7.90, 3.32, 5.37) # Effect size: mean difference mean.diff <- mean.NaF - mean.SMFP std.meandiff <- sqrt(std.NaF^2 / n.NaF + std.SMFP^2 / n.SMFP) # Generic inverse variance method using R package meta library(meta) # Combining mean differences using functions metagen and metacont comb.meandiff <- metagen(mean.diff, std.meandiff) comb.meandiff metacont(n.NaF, mean.NaF, std.NaF, n.SMFP, mean.SMFP, std.SMFP,sm = "WMD") improved.ci(mean.diff, std.meandiff, 0.05) # Confidence interval plot with fixed and random effects meta-analysis plot(comb.meandiff, comb.f = TRUE, comb.r = TRUE, main = "Dentifrice studies", xlab = "Difference of normal means") # Funnel plot and tests for publication bias # method = rank (rank correlation) # method = linear (linear regression) funnel(comb.meandiff, main = "Dentifrice studies", xlab = "Difference of normal means") metabias(comb.meandiff, method = "rank") metabias(comb.meandiff, method = "linreg") ########################################################################## # Combining standardized mean differences using the function metacont comb.smd <- metacont(n.NaF, mean.NaF, std.NaF, n.SMFP, mean.SMFP, std.SMFP,sm = "SMD") comb.smd # Confidence interval plot with fixed and random effects meta-analysis plot(comb.smd, comb.f = TRUE, comb.r = TRUE, main = "Dentifrice studies", xlab = "Standardized difference of normal means") # Funnel plot and tests for publication bias # method = rank (rank correlation) # method = linear (linear regression) funnel(comb.smd, main = "Dentifrice studies", xlab = "Standardized difference of normal means") metabias(comb.smd, method = "rank") metabias(comb.smd, method = "linreg")