RobustBayesianCopas: 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,

$$y_i | (z_i>0) = \theta + \tau u_i + s_i \epsilon_i,$$ $$z_i = \gamma_0 + \gamma_1 / s_i + \delta_i,$$ $$corr(\epsilon_i, \delta_i) = \rho,$$

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, and $\epsilon_i$, and $\delta_i$ are marginally distributed as $N(0,1)$ and $u_i$, and $\epsilon_i$ are independent.

In the Copas selection model, $y_i$ is published (selected) if and only if the corresponding propensity score $z_i$ (or the propensity to publish) is greater than zero. The propensity score $z_i$ contains two parameters: $\gamma_0$ controls the overall probability of publication, and $\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 $\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 $(\theta, \tau^2, \rho, \gamma_0, \gamma_1)$ (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. 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

An $n \times 1$ vector of reported treatment effects.

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, \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 $\gamma_0$.

gamma0.samples

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

gamma1.hat

Posterior median for $\gamma_1$.

gamma1.samples

MCMC samples for $\gamma_1$ after burn-in. Used for plotting the posterior for $\gamma_1$ and performing inference of $\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

Run this code

# NOT RUN {
######################################
# 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)

# }
# NOT RUN {

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

Run the code above in your browser using DataLab