RobustBayesianCopas {RobustBayesianCopas}R Documentation

Robust Bayesian Copas selection model

Description

This function implements the Robust Bayesian Copas selection model of Bai et al. (2020) for the Copas selection model,

yi(zi>0)=θ+τui+siϵi,y_i | (z_i>0) = \theta + \tau u_i + s_i \epsilon_i,

zi=γ0+γ1/si+δi,z_i = \gamma_0 + \gamma_1 / s_i + \delta_i,

corr(ϵi,δi)=ρ,corr(\epsilon_i, \delta_i) = \rho,

where yiy_i is the reported treatment effect for the iith study, sis_i is the reported standard error for the iith study, θ\theta is the population treatment effect of interest, τ>0\tau > 0 is a heterogeneity parameter, and ϵi\epsilon_i, and δi\delta_i are marginally distributed as N(0,1)N(0,1) and uiu_i, and ϵi\epsilon_i are independent.

In the Copas selection model, yiy_i is published (selected) if and only if the corresponding propensity score ziz_i (or the propensity to publish) is greater than zero. The propensity score ziz_i contains two parameters: γ0\gamma_0 controls the overall probability of publication, and γ1\gamma_1 controls how the probability of publication depends on study sample size. The reported treatment effects and propensity scores are correlated through ρ\rho. If ρ=0\rho=0, then there is no publication bias and the Copas selection model reduces to the standard random effects meta-analysis model.

The RBC model places noninformative priors on (θ,τ2,ρ,γ0,γ1)(\theta, \tau^2, \rho, \gamma_0, \gamma_1) (see Bai et al. (2020) for details). For the random effects ui,i=1,,nu_i, i=1, \ldots, n, we give the option for using normal, Student's t, Laplace, or slash distributions for the random effects. The function returns the Deviance Information Criterion (DIC), which can be used to select the appropriate distribution to use for the final analysis.

Usage

RobustBayesianCopas(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×1n \times 1 vector of reported treatment effects.

s

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

re.dist

Distribution for the between-study random effects ui,i=1,,nu_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 (θ,τ,ρ,γ0,γ1)(\theta, \tau, \rho, \gamma_0, \gamma_1). 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:

DIC

Deviance Information Criterion (DIC), a measure of model fit. This can be used to compare the results for different random effects distributions. The model that gives the lowest DIC gives the best fit to the data.

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.

rho.hat

Posterior median for ρ\rho.

rho.samples

MCMC samples for ρ\rho after burn-in. Used for plotting the posterior for ρ\rho and performing inference of ρ\rho.

gamma0.hat

Posterior median for γ0\gamma_0.

gamma0.samples

MCMC samples for γ0\gamma_0 after burn-in. Used for plotting the posterior for γ0\gamma_0 and performing inference of γ0\gamma_0.

gamma1.hat

Posterior median for γ1\gamma_1.

gamma1.samples

MCMC samples for γ1\gamma_1 after burn-in. Used for plotting the posterior for γ1\gamma_1 and performing inference of γ1\gamma_1.

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

######################################
# Example on the Barlow2014 data set #
######################################
# Load data
data(Barlow2014)
attach(Barlow2014)

# Observed treatment effect
y.obs = Barlow2014[,1]
# Observed standard error
s.obs = Barlow2014[,2]

###############################################
# Fit the RBC model with slash random effects #
###############################################
# NOTE: Use default burn-in (burn=10000) and post-burn-in samples (nmc=10000)
# Fit model with slash errors
RBC.mod = RobustBayesianCopas(y=y.obs, s=s.obs, re.dist="slash", burn=500, nmc=500)

# Point estimate for rho
rho.hat.RBC = RBC.mod$rho.hat
rho.hat.RBC
# Plot posterior for rho
hist(RBC.mod$rho.samples)

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

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

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




############################################
# Example on second-hand smoking data set. #
# This is from Section 6.1 of the paper by #
# Bai et al. (2020).                       #
############################################

# Set seed, so we can reproduce the exact same result as in the paper.
seed = 1234
set.seed(seed)

# Load the full data
data(Hackshaw1997)
attach(Hackshaw1997)
# Extract the log OR
y.obs = Hackshaw1997[,2]
# Extract the observed standard error
s.obs = Hackshaw1997[,3]


###################################################
# Fit the RBC model with different random effects #
# distributions and compare them using the DIC.   #
###################################################

# Normal
RBC.mod.normal = RobustBayesianCopas(y=y.obs, s=s.obs, re.dist="normal", seed=seed) 
RBC.mod.normal$DIC # DIC=429.7854
# Student's t
RBC.mod.StudentsT = RobustBayesianCopas(y=y.obs, s=s.obs, re.dist="StudentsT", seed=seed) 
RBC.mod.StudentsT$DIC # DIC=399.1955
# Laplace
RBC.mod.Laplace = RobustBayesianCopas(y=y.obs, s=s.obs, re.dist="Laplace", seed=seed) 
RBC.mod.Laplace$DIC # DIC=410.9086
# Slash
RBC.mod.slash = RobustBayesianCopas(y=y.obs, s=s.obs, re.dist="slash", seed=seed)
RBC.mod.slash$DIC # DIC=407.431


#######################################################
# Use the model with t-distributed random errors for  #
# the final analysis since it gave the lowest DIC.    #
#######################################################

# Point estimate for rho
rho.hat.RBC = RBC.mod.StudentsT$rho.hat # rho.hat=0.459 (moderate publication bias)
# Plot posterior for rho
hist(RBC.mod.StudentsT$rho.samples)

# Point estimate for theta
theta.hat.RBC = RBC.mod.StudentsT$theta.hat # theta.hat=0.1672
# 95% posterior credible interval for theta
theta.cred.int = quantile(RBC.mod.StudentsT$theta.samples, probs=c(0.025,0.975))
# Plot the posterior for theta
hist(RBC.mod.StudentsT$theta.samples)

# Obtain odds ratio estimates
OR.samples.RBC = exp(RBC.mod.StudentsT$theta.samples) # Samples of exp(theta)
# Posterior mean OR
OR.RBC.hat = mean(OR.samples.RBC) # OR.hat=1.185
# 95% posterior credible interval for OR
OR.RBC.credint = quantile(OR.samples.RBC, probs=c(0.025,0.975)) # (1.018, 1.350)


##############################################
# Use D measure to quantify publication bias #
##############################################

# Make sure that we specify the random effects as Student's t, since that is
# the distribution that we used for our final analysis.
Bayes.nobias.mod = BayesNonBiasCorrected(y=y.obs, s=s.obs, re.dist="StudentsT", seed=seed)

# Compute D measure based on posterior samples of theta
D.RBC = D.measure(RBC.mod.StudentsT$theta.samples, Bayes.nobias.mod$theta.samples)
D.RBC # D=0.33


[Package RobustBayesianCopas version 2.0 Index]