DearBegg: Compute the nonparametric weight function from Dear and Begg...
In selectMeta: Estimation of Weight Functions in Meta Analysis

Description Usage Arguments Value Author(s) References See Also Examples

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.

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)

`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.

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`.

Kaspar Rufibach (maintainer), kaspar.rufibach@gmail.com,
http://www.kasparrufibach.ch

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.

IyenGreen for a parametric selection model.

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