Learn R Programming

selectMeta (version 1.0.8)

DearBegg: Compute the nonparametric weight function from Dear and Begg (1992)

Description

In Dear and Begg (1992) it was proposed to nonparametrically estimate via maximum likelihood the weight function $w$ in a selection model via pooling $p$-values in groups of 2 and assuming a piecewise constant $w$. The function DearBegg implements estimation of $w$ via a coordinate-wise Newton-Raphson algorithm as described in Dear and Begg (1992). In addition, the function DearBeggMonotone enables computation of the weight function in the same model under the constraint that it is non-increasing, see Rufibach (2011). To this end we use the differential evolution algorithm described in Ardia et al (2010a, b) and implemented in Mullen et al (2009). The functions Hij, DearBeggLoglik, and DearBeggToMinimize are not intended to be called by the user.

Usage

DearBegg(y, u, lam = 2, tolerance = 10^-10, maxiter = 1000, trace = TRUE) DearBeggMonotone(y, u, lam = 2, maxiter = 1000, CR = 0.9, NP = NA, trace = TRUE) Hij(theta, sigma, y, u, teststat) DearBeggLoglik(w, theta, sigma, y, u, hij, lam) DearBeggToMinimize(vec, y, u, lam)

Arguments

y
Normally distributed effect sizes.
u
Associated standard errors.
lam
Weight of the first entry of $w$ in the likelihood function. Dear and Begg (1992) recommend to use lam = 2.
tolerance
Stopping criterion for Newton-Raphson.
maxiter
Maximal number of iterations for Newton-Raphson.
trace
If TRUE, progress of the algorithm is shown.
CR
Parameter that is given to DEoptim. See the help file of the function DEoptim.control for details.
NP
Parameter that is given to DEoptim. See the help file of the function DEoptim.control for details.
w
Weight function, parametrized as vector of length $1 + \lfloor(n / 2)\rfloor$ where $n$ is the number of studies, i.e. the length of $y$.
theta
Effect size estimate.
sigma
Random effects variance component.
hij
Integral of density over a constant piece of $w$. See Rufibach (2011, Appendix) for details.
vec
Vector of parameters over which we maximize.
teststat
Vector of test statistics, equals $|y| / u$.

Value

w
Vector of estimated weights.
theta
Estimate of the combined effect in the Dear and Begg model.
sigma
Estimate of the random effects component variance.
p
$p$-values computed from the inputed test statistics, ordered in decreasing order.
y
Effect sizes, ordered in decreasing order of $p$-values.
u
Standard errors, ordered in decreasing order of $p$-values.
loglik
Value of the log-likelihood at the maximum.
DEoptim.res
Only available in DearBeggMonotone. Provides the object that is outputted by DEoptim.

References

Ardia, D., Boudt, K., Carl, P., Mullen, K.M., Peterson, B.G. (2010). Differential Evolution ('DEoptim') for Non-Convex Portfolio Optimization.

Ardia, D., Mullen, K.M., et.al. (2010). The 'DEoptim' Package: Differential Evolution Optimization in 'R'. Version 2.0-7.

Dear, K.B.G. and Begg, C.B. (1992). An Approach for Assessing Publication Bias Prior to Performing a Meta-Analysis. Statist. Sci., 7(2), 237--245.

Mullen, K.M., Ardia, D., Gil, D.L., Windover, D., Cline, J. (2009). 'DEoptim': An 'R' Package for Global Optimization by Differential Evolution.

Rufibach, K. (2011). Selection Models with Monotone Weight Functions in Meta-Analysis. Biom. J., 53(4), 689--704.

See Also

IyenGreen for a parametric selection model.

Examples

Run this code
## Not run: 
# ##------------------------------------------
# ## 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)
# ## End(Not run)

Run the code above in your browser using DataLab