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,
where is the reported treatment effect for the
th study,
is the reported standard error for the
th study,
is the population treatment effect of interest,
is a heterogeneity parameter, and
, and
are marginally distributed as
and
, and
are independent.
In the Copas selection model, is published (selected) if and only if the corresponding propensity score
(or the propensity to publish) is greater than zero. The propensity score
contains two parameters:
controls the overall probability of publication, and
controls how the probability of publication depends on study sample size. The reported treatment effects and propensity scores are correlated through
. If
, 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 (see Bai et al. (2020) for details). For the random effects
, 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 |
s |
An |
re.dist |
Distribution for the between-study random effects |
t.df |
Degrees of freedom for t-distribution. Only used if |
slash.shape |
Shape parameter in the slash distribution. Only used if |
init |
Optional initialization values for |
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 |
nmc |
Number of posterior samples to save. Default is |
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.samples |
MCMC samples for |
tau.hat |
Posterior mean for |
tau.samples |
MCMC samples for |
rho.hat |
Posterior median for |
rho.samples |
MCMC samples for |
gamma0.hat |
Posterior median for |
gamma0.samples |
MCMC samples for |
gamma1.hat |
Posterior median for |
gamma1.samples |
MCMC samples for |
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