DearBegg {selectMeta}

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

A list consisting of the following elements:

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.

Author(s)

Kaspar Rufibach (maintainer), [e-mail],
http://www.kasparrufibach.ch

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

## 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)