if (FALSE) {
##------------------------------------------
## Analysis of Hedges & Olkin dataset
## re-analyzed in Iyengar & Greenhouse, Dear & Begg
##------------------------------------------
data(education)
t <- education$t
q <- education$q
N <- education$N
y <- education$theta
u <- sqrt(2 / N)
n <- length(y)
k <- 1 + floor(n / 2)
lam1 <- 2
## compute p-values
p <- 2 * pnorm(-abs(t))
##------------------------------------------
## compute all weight functions available
## in this package
##------------------------------------------
## weight functions from Iyengar & Greenhouse (1988)
res1 <- IyenGreenMLE(t, q, N, type = 1)
res2 <- IyenGreenMLE(t, q, N, type = 2)
## weight function from Dear & Begg (1992)
res3 <- DearBegg(y, u, lam = lam1)
## monotone version of Dear & Begg, as introduced in Rufibach (2011)
set.seed(1977)
res4 <- DearBeggMonotone(y, u, lam = lam1, maxiter = 1000, CR = 1)
## plot
plot(0, 0, type = "n", xlim = c(0, 1), ylim = c(0, 1), xlab = "p-values",
ylab = "estimated weight function")
ps <- seq(0, 1, by = 0.01)
rug(p, lwd = 3)
lines(ps, IyenGreenWeight(-qnorm(ps / 2), b = res1$beta, q = 50,
type = 1, alpha = 0.05), lwd = 3, col = 2)
lines(ps, IyenGreenWeight(-qnorm(ps / 2), b = res2$beta, q = 50,
type = 2, alpha = 0.05), lwd = 3, col = 4)
weightLine(p, w = res3$w, col0 = 3, lwd0 = 3, lty0 = 2)
weightLine(p, w = res4$w, col0 = 6, lwd0 = 2, lty0 = 1)
legend("topright", c(expression("Iyengar & Greenhouse (1988) w"[1]),
expression("Iyengar & Greenhouse (1988) w"[2]), "Dear and Begg (1992)",
"Rufibach (2011)"), col = c(2, 4, 3, 6), lty = c(1, 1, 2, 1),
lwd = c(3, 3, 3, 2), bty = "n")
## compute selection bias
eta <- sqrt(res4$sigma ^ 2 + res4$u ^ 2)
bias <- effectBias(res4$y, res4$u, res4$w, res4$theta, eta)
bias
##------------------------------------------
## Compute p-value to assess null hypothesis of no selection,
## as described in Rufibach (2011, Section 6)
## We use the package 'meta' to compute initial estimates for
## theta and sigma
##------------------------------------------
library(meta)
## compute null parameters
meta.edu <- metagen(TE = y, seTE = u, sm = "MD", level = 0.95,
comb.fixed = TRUE, comb.random = TRUE)
theta0 <- meta.edu$TE.random
sigma0 <- meta.edu$tau
M <- 1000
res <- DearBeggMonotonePvalSelection(y, u, theta0, sigma0, lam = lam1,
M = M, maxiter = 1000)
## plot all the computed monotone functions
plot(0, 0, xlim = c(0, 1), ylim = c(0, 1), type = "n", xlab = "p-values",
ylab = expression(w(p)))
abline(v = 0.05, lty = 3)
for (i in 1:M){weightLine(p, w = res$res.mono[1:k, i], col0 = grey(0.8),
lwd0 = 1, lty0 = 1)}
rug(p, lwd = 2)
weightLine(p, w = res$mono0, col0 = 2, lwd0 = 1, lty0 = 1)
## =======================================================================
##------------------------------------------
## Analysis second-hand tobacco smoke dataset
## Rothstein et al (2005), Publication Bias in Meta-Analysis, Appendix A
##------------------------------------------
data(passive_smoking)
u <- passive_smoking$selnRR
y <- passive_smoking$lnRR
n <- length(y)
k <- 1 + floor(n / 2)
lam1 <- 2
res2 <- DearBegg(y, u, lam = lam1)
set.seed(1)
res3 <- DearBeggMonotone(y = y, u = u, lam = lam1, maxiter = 2000, CR = 1)
plot(0, 0, type = "n", xlim = c(0, 1), ylim = c(0, 1), pch = 19, col = 1,
xlab = "p-values", ylab = "estimated weight function")
weightLine(rev(sort(res2$p)), w = res2$w, col0 = 2, lwd0 = 3, lty0 = 2)
weightLine(rev(sort(res3$p)), w = res3$w, col0 = 4, lwd0 = 2, lty0 = 1)
legend("bottomright", c("Dear and Begg (1992)", "Rufibach (2011)"), col =
c(2, 4), lty = c(2, 1), lwd = c(3, 2), bty = "n")
## compute selection bias
eta <- sqrt(res3$sigma ^ 2 + res3$u ^ 2)
bias <- effectBias(res3$y, res3$u, res3$w, res3$theta, eta)
bias
##------------------------------------------
## Compute p-value to assess null hypothesis of no selection
##------------------------------------------
## compute null parameters
meta.toba <- metagen(TE = y, seTE = u, sm = "MD", level = 0.95,
comb.fixed = TRUE, comb.random = TRUE)
theta0 <- meta.toba$TE.random
sigma0 <- meta.toba$tau
M <- 1000
res <- DearBeggMonotonePvalSelection(y, u, theta0, sigma0, lam = lam1,
M = M, maxiter = 2000)
## plot all the computed monotone functions
plot(0, 0, xlim = c(0, 1), ylim = c(0, 1), type = "n", xlab = "p-values",
ylab = expression(w(p)))
abline(v = 0.05, lty = 3)
for (i in 1:M){weightLine(p, w = res$res.mono[1:k, i], col0 = grey(0.8),
lwd0 = 1, lty0 = 1)}
rug(p, lwd = 2)
weightLine(p, w = res$mono0, col0 = 2, lwd0 = 1, lty0 = 1)
}
Run the code above in your browser using DataLab