RobustBayesianCopas: Robust Bayesian Copas selection model
In RobustBayesianCopas: Robust Bayesian Copas Selection Model
Description
Usage
Arguments
Value
References
Examples
View source: R/RobustBayesianCopas.R
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) = θ + τ u_i + s_i ε_i,
z_i = γ_0 + γ_1 / s_i + δ_i,
corr(ε_i, δ_i) = ρ,
where y_i is the reported treatment effect for the ith study, s_i is the reported standard error for the ith study, θ is the population treatment effect of interest, τ > 0 is a heterogeneity parameter, and ε_i, and δ_i are marginally distributed as N(0,1) and u_i, and ε_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: γ_0 controls the overall probability of publication, and γ_1 controls how the probability of publication depends on study sample size. The reported treatment effects and propensity scores are correlated through ρ. If ρ=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) (see Bai et al. (2020) for details). For the random effects u_i, i=1, …, 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
1
2
3
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 \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, …, 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). 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.samples
MCMC samples for θ after burn-in. Used for plotting the posterior for θ and performing inference of θ.
tau.hat
Posterior mean for τ.
tau.samples
MCMC samples for τ after burn-in. Used for plotting the posterior for τ and performing inference of τ.
rho.hat
Posterior median for ρ.
rho.samples
MCMC samples for ρ after burn-in. Used for plotting the posterior for ρ and performing inference of ρ.
gamma0.hat
Posterior median for γ_0.
gamma0.samples
MCMC samples for γ_0 after burn-in. Used for plotting the posterior for γ_0 and performing inference of γ_0.
gamma1.hat
Posterior median for γ_1.
gamma1.samples
MCMC samples for γ_1 after burn-in. Used for plotting the posterior for γ_1 and performing inference of γ_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
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
######################################
# 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
RobustBayesianCopas documentation built on Jan. 13, 2021, 12:50 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Description Usage Arguments Value References Examples
View source: R/RobustBayesianCopas.R
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) = θ + τ u_i + s_i ε_i,
z_i = γ_0 + γ_1 / s_i + δ_i,
corr(ε_i, δ_i) = ρ,
where y_i is the reported treatment effect for the ith study, s_i is the reported standard error for the ith study, θ is the population treatment effect of interest, τ > 0 is a heterogeneity parameter, and ε_i, and δ_i are marginally distributed as N(0,1) and u_i, and ε_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: γ_0 controls the overall probability of publication, and γ_1 controls how the probability of publication depends on study sample size. The reported treatment effects and propensity scores are correlated through ρ. If ρ=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) (see Bai et al. (2020) for details). For the random effects u_i, i=1, …, 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
1 2 3 | 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 \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, …, n. The user may specify the normal, Student's T, Laplace, or slash distribution. The default is |
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 (θ, τ, ρ, γ_0, γ_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 |
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 θ after burn-in. Used for plotting the posterior for θ and performing inference of θ. |
tau.hat |
Posterior mean for τ. |
tau.samples |
MCMC samples for τ after burn-in. Used for plotting the posterior for τ and performing inference of τ. |
rho.hat |
Posterior median for ρ. |
rho.samples |
MCMC samples for ρ after burn-in. Used for plotting the posterior for ρ and performing inference of ρ. |
gamma0.hat |
Posterior median for γ_0. |
gamma0.samples |
MCMC samples for γ_0 after burn-in. Used for plotting the posterior for γ_0 and performing inference of γ_0. |
gamma1.hat |
Posterior median for γ_1. |
gamma1.samples |
MCMC samples for γ_1 after burn-in. Used for plotting the posterior for γ_1 and performing inference of γ_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
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 | ######################################
# 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
|
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.