RobustBayesianCopas (version 2.0)

BayesNonBiasCorrected: Non-bias-corrected robust Bayesian meta-analysis model

Description

This function implements the non-bias-corrected Robust Bayesian Copas selection model of Bai et al. (2020) when there is no publication bias (i.e. \(\rho=0\)). In this case, the Copas selection model reduces to the standard random effects meta-analysis model:

$$y_i = \theta + \tau u_i + s_i \epsilon_i,$$

where \(y_i\) is the reported treatment effect for the \(i\)th study, \(s_i\) is the reported standard error for the \(i\)th study, \(\theta\) is the population treatment effect of interest, \(\tau > 0\) is a heterogeneity parameter, \(\epsilon_i\) is distributed as \(N(0,1)\), and \(u_i\) and \(\epsilon_i\) are independent.

For the non-bias-corrected model, we place noninformative priors on \((\theta, \tau^2)\) (see Bai et al. (2020) for details). For the random effects \(u_i, i=1, \ldots, n\), we give the option for using normal, Student's t, Laplace, or slash distributions for the random effects. If this function is being run in order to quantify publication bias with the robust Bayesian Copas selection model, then the practitioner should use the same random effects distribution that they used for RobustBayesianCopas.

Usage

BayesNonBiasCorrected(y, s, re.dist=c("normal", "StudentsT", "Laplace", "slash"),
                      t.df = 4, slash.shape = 1, init=NULL, seed=NULL,
                      burn=10000, nmc=10000)

Arguments

y

An \(n \times 1\) vector of reported treatment effects.

s

An \(n \times 1\) vector of reported within-study standard errors.

re.dist

Distribution for the between-study random effects \(u_i, i=1, \ldots, n\). The user may specify the normal, Student's t, Laplace, or slash distribution. The default is StudentsT with 4 degrees of freedom.

t.df

Degrees of freedom for t-distribution. Only used if StudentsT is specified for re.dist. Default is 4.

slash.shape

Shape parameter in the slash distribution. Only used if slash is specified for re.dist. Default is 1.

init

Optional initialization values for \((\theta, \tau)\). If specified, they must be provided in this exact order. If they are not provided, the program estimates initial values from the data.

seed

Optional seed. This needs to be specified if you want to reproduce the exact results of your analysis.

burn

Number of burn-in samples. Default is burn=10000.

nmc

Number of posterior samples to save. Default is nmc=10000.

Value

The function returns a list containing the following components:

theta.hat

posterior mean for \(\theta\).

theta.samples

MCMC samples for \(\theta\) after burn-in. Used for plotting the posterior for \(\theta\) and performing inference of \(\theta\).

tau.hat

posterior mean for \(\tau\).

tau.samples

MCMC samples for \(\tau\) after burn-in. Used for plotting the posterior for \(\tau\) and performing inference of \(\tau\).

References

Bai, R., Lin, L., Boland, M. R., and Chen, Y. (2020). "A robust Bayesian Copas selection model for quantifying and correcting publication bias." arXiv preprint arXiv:2005.02930.

Examples

# NOT RUN {
######################################
# Example on the Barlow2014 data set #
######################################
data(Barlow2014)
attach(Barlow2014)
# Observed treatment effect
y.obs = Barlow2014[,1]
# Observed standard error
s.obs = Barlow2014[,2]

####################################
# Fit the non-bias-corrected model #
####################################

# NOTE: Use default burn-in (burn=10000) and post-burn-in samples (nmc=10000)
# Fit the model with Laplace errors.
RBCNoBias.mod = BayesNonBiasCorrected(y=y.obs, s=s.obs, re.dist="Laplace", burn=500, nmc=500)

# Point estimate for theta
theta.hat.RBCNoBias = RBCNoBias.mod$theta.hat 
# Standard error for theta
theta.se.RBCNoBias = sd(RBCNoBias.mod$theta.samples)
# 95% posterior credible interval for theta
theta.cred.int = quantile(RBCNoBias.mod$theta.samples, probs=c(0.025,0.975))

# Display results
theta.hat.RBCNoBias 
theta.se.RBCNoBias  
theta.cred.int      

# Plot the posterior for theta
hist(RBCNoBias.mod$theta.samples)

# }