Sample R Code for the BayesCACE Package
Jincheng Zhou
6/20/2021
1. Install the R package and view the datasets
The R package BayesCACE provides user-friendly functions to estimate CACE in either a single study or meta-analysis using the Baysian methods. The R package source files are available at GitHub: https://github.com/JinchengZ/BayesCACE. he BayesCACE package depends on the R packages rjags Plummer (2018), coda Plummer et al. (2006), and forestplot Gordon and Lumley (2017).
First, installing the BayeCACE R package using the following R code:
## R package BayeCACE installation
library(devtools)## Loading required package: usethisSys.setenv("TAR" = "internal")
devtools::install_github("JinchengZ/BayesCACE")## Downloading GitHub repo JinchengZ/BayesCACE@master## Warning in untar2(tarfile, files, list, exdir, restore_times): skipping pax
## global extended headers
## Warning in untar2(tarfile, files, list, exdir, restore_times): skipping pax
## global extended headers##
##
checking for file 'C:\Users\jzhou02\AppData\Local\Temp\RtmpwFIqkL\remotes2e6055029c\JinchengZ-BayesCACE-111be7f/DESCRIPTION' ...
v checking for file 'C:\Users\jzhou02\AppData\Local\Temp\RtmpwFIqkL\remotes2e6055029c\JinchengZ-BayesCACE-111be7f/DESCRIPTION' (826ms)
##
- preparing 'BayesCACE':
## checking DESCRIPTION meta-information ...
checking DESCRIPTION meta-information ...
v checking DESCRIPTION meta-information
##
- installing the package to process help pages
##
Loading required namespace: BayesCACE
##
- saving partial Rd database (2.2s)
##
- checking for LF line-endings in source and make files and shell scripts
##
- checking for empty or unneeded directories
##
- looking to see if a 'data/datalist' file should be added
##
- building 'BayesCACE_1.0.tar.gz'
##
## ## Installing package into 'C:/Users/jzhou02/Documents/R/win-library/3.6'
## (as 'lib' is unspecified)Next, open the package in your R environment and view the default datasets. The one with complete data (10 studies in total) is called epidural_c. The other one that also include incomplete compliance data (27 studies in total) is named as epidural_ic. These two data sets were obtained from Bannister-Tyrrell et al. (2015), who conducted an exploratory meta-analysis of the association between using epidural analgesia in labor and the risk of cesarean section.
library("BayesCACE")## Loading required package: rjags## Loading required package: coda## Linked to JAGS 4.3.0## Loaded modules: basemod,bugsThe dataset epidural_c contains 10 trials with full compliance information; each trial has 8 observed counts, denoted by \(N_{irto}\) and presented in columns Nirto for \(i=1, \dots, 10\) and \(r, t, o \in \{0, 1\}\). study.id contains IDs for the 10 studies, and study.name labels each study by its first author’s surname and its publication year.
data("epidural_c", package = "BayesCACE")
epidural_c## study.id study.name n000 n001 n010 n011 n100 n101 n110 n111
## 1 1 Bofill, 1997 37 2 11 1 2 0 42 5
## 2 2 Clark, 1998 72 6 68 16 7 2 134 13
## 3 3 Halpern, 2004 62 5 44 7 0 0 112 12
## 4 4 Head, 2002 51 7 2 0 3 0 43 10
## 5 5 Jain, 2003 72 11 0 0 0 2 36 7
## 6 6 Nafisi, 2006 179 19 0 0 0 0 173 24
## 7 7 Nikkola, 1997 6 0 4 0 0 0 10 0
## 8 8 Ramin, 1995 546 17 95 8 230 2 393 39
## 9 9 Sharma, 1997 336 16 5 0 114 1 231 12
## 10 10 Volmanen, 2008 23 1 3 0 1 0 23 1The other dataset, epidural_ic, represents the situation in which not all trials report complete compliance data. It contains 27 studies, only 10 out of which have full compliance information and were included in epidural_c. Each study is represented by one row in the dataset; the columns study.id and study.name have the same meanings as in the dataset epidural_c. Each study’s data are summarized in 12 numbers (columns) denoted by \(N_{irto}\) and \(N_{ir*o}\). For a particular randomization group \(r \in \{0, 1\}\), the observed counts are presented either as \(N_{irto}\) or \(N_{ir*o}\) depending on whether the compliance information is available; values for other columns are denoted by 0. The corresponding column names in the dataset are Nirto and Nirso, respectively.
data("epidural_ic", package = "BayesCACE")
head(epidural_ic)## study.id study.name n000 n001 n010 n011 n0s0 n0s1 n100 n101 n110 n111
## 1 1 Bofill, 1997 37 2 11 1 0 0 2 0 42 5
## 2 2 Clark, 1998 72 6 68 16 0 0 7 2 134 13
## 3 3 Dickinson, 2002 0 0 0 0 428 71 0 0 0 0
## 4 4 Evron, 2008 40 4 0 0 0 0 0 0 0 0
## 5 5 El Kerdawy, 2010 0 0 0 0 12 3 0 0 0 0
## 6 6 Gambling, 1998 0 0 0 0 573 34 206 10 371 29
## n1s0 n1s1
## 1 0 0
## 2 0 0
## 3 408 85
## 4 129 19
## 5 11 4
## 6 0 0Note that NA is not allowed in a dataset for the package BayesCACE, but some trials may have 0 events or 0 noncompliance rates.
Before performing the CACE analysis, one might want a visual overview of study-specific noncompliance rates in both randomization arms. The function plot.noncomp provides a forest plot of noncompliance rates in an R plot window.
plot.noncomp(data = epidural_c, overall = TRUE)
The red dot with its horizontal line shows the study-specific noncompliance rate with its 95% exact confidence interval for the patients randomized to the treatment arm, and the blue square with its horizontal line represents that rate and interval for those in the control arm. The confidence intervals are calculated by the Clopper–Pearson exact method , which is based on the cumulative distribution function of the binomial distribution. Using the default overall = TRUE, the figure also gives a summary estimate of the compliance rates per randomization group. This overall rate is estimated using a logit generalized linear mixed model. Otherwise, if the argument overall is FALSE, the plot shows only study-specific noncompliance rates.
2. CACE for a single trial with noncompliance
To estimate CACE for a single study, users need to input data with the same structure as epidural_c, containing either one row of observations for a single study, or multiple rows referring to multiple studies in a meta-analysis. This function fits a model for a single study. If the data includes more than one study, the study-specific CACEs will be estimated by retrieving data row by row.
If users do not specify their own prior distributions, the default priors are used.
set.seed(123)
out.study <- cace.study(data = epidural_c, conv.diag = TRUE, mcmc.samples =
+ TRUE, two.step = TRUE)## NA is not allowed in the input data set;
## the rows containing NA are removed.
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 2
## Unobserved stochastic nodes: 6
## Total graph size: 44
##
## Initializing model## module dic loaded## MCMC convergence diagnostic statistics are calculated and saved in conv.out
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 2
## Unobserved stochastic nodes: 6
## Total graph size: 44
##
## Initializing model
##
## MCMC convergence diagnostic statistics are calculated and saved in conv.out
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 2
## Unobserved stochastic nodes: 6
## Total graph size: 44
##
## Initializing model
##
## MCMC convergence diagnostic statistics are calculated and saved in conv.out
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 2
## Unobserved stochastic nodes: 6
## Total graph size: 44
##
## Initializing model
##
## MCMC convergence diagnostic statistics are calculated and saved in conv.out
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 2
## Unobserved stochastic nodes: 6
## Total graph size: 44
##
## Initializing model
##
## MCMC convergence diagnostic statistics are calculated and saved in conv.out
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 2
## Unobserved stochastic nodes: 6
## Total graph size: 44
##
## Initializing model
##
## MCMC convergence diagnostic statistics are calculated and saved in conv.out
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 2
## Unobserved stochastic nodes: 6
## Total graph size: 44
##
## Initializing model
##
## MCMC convergence diagnostic statistics are calculated and saved in conv.out
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 2
## Unobserved stochastic nodes: 6
## Total graph size: 44
##
## Initializing model
##
## MCMC convergence diagnostic statistics are calculated and saved in conv.out
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 2
## Unobserved stochastic nodes: 6
## Total graph size: 44
##
## Initializing model
##
## MCMC convergence diagnostic statistics are calculated and saved in conv.out
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 2
## Unobserved stochastic nodes: 6
## Total graph size: 44
##
## Initializing model
##
## MCMC convergence diagnostic statistics are calculated and saved in conv.outIf the dataset contains more than one study, e.g., the epidural_c dataset has 10 trials, then once the JAGS model compiles for the first study, it automatically continues to run on the next study’s data. The results are saved in the object out.study, a list containing the model name, posterior information for each monitored parameter, and DIC of each study.
For example, the estimates of \(\theta^\text{CACE}\) for each single study (posterior mean and standard deviation, posterior median, 95% credible interval, and time-series standard error) can be obtained by
out.study$CACE## Mean SD 2.5% 50% 97.5% Naive SE Time-series SE
## [1,] 0.04980 0.0797 -0.09510 4.46e-02 0.2180 1.45e-04 2.51e-04
## [2,] -0.02490 0.0489 -0.12200 -2.23e-02 0.0785 8.94e-05 1.48e-04
## [3,] -0.02210 0.0606 -0.12700 -2.90e-02 0.1120 1.11e-04 1.94e-04
## [4,] 0.07180 0.0758 -0.07550 7.10e-02 0.2230 1.38e-04 2.01e-04
## [5,] 0.08250 0.0768 -0.06260 8.11e-02 0.2370 1.40e-04 2.51e-04
## [6,] 0.02600 0.0319 -0.03650 2.59e-02 0.0891 5.83e-05 7.55e-05
## [7,] 0.01430 0.1580 -0.28200 2.11e-04 0.4050 2.89e-04 4.21e-04
## [8,] 0.05030 0.0248 0.00176 5.02e-02 0.0993 4.54e-05 7.34e-05
## [9,] -0.01100 0.0234 -0.05740 -1.09e-02 0.0350 4.27e-05 6.22e-05
## [10,] 0.00145 0.0655 -0.13400 -4.36e-06 0.1460 1.20e-04 1.56e-04If the argument conv.diag is specified as TRUE, the output list contains a sub-list conv.out, which outputs the point estimates of the ‘potential scale reduction factor’ (the Gelman and Rubin convergence statistic, labelled Point est.) calculated for each parameter from each single study, and their upper confidence limits (labelled Upper C.I.).
For example, the first sub-list from conv.out is
out.study$conv.out[[1]]## Point est. Upper C.I.
## CACE 1.000025 1.000046
## b1 1.000041 1.000129
## pi.a 1.000025 1.000094
## pi.c 1.000036 1.000134
## pi.n 1.000029 1.000067
## s1 1.000014 1.000018
## u1 1.000016 1.000033
## v1 1.000077 1.000185If the dataset used by the function cace.study() has more than one study, specifying the argument two.step = TRUE causes the two-step meta-analysis for \(\theta^\text{CACE}\) to be done. The outcomes are saved as a sub-list object meta. Note that users can obtain different meta-analysis estimators by changing the method argument in cace.study().
out.study$meta##
## Random-Effects Model (k = 10; tau^2 estimator: REML)
##
## tau^2 (estimated amount of total heterogeneity): 0.0002 (SE = 0.0008)
## tau (square root of estimated tau^2 value): 0.0131
## I^2 (total heterogeneity / total variability): 8.10%
## H^2 (total variability / sampling variability): 1.09
##
## Test for Heterogeneity:
## Q(df = 9) = 5.9353, p-val = 0.7464
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.0182 0.0143 1.2758 0.2020 -0.0098 0.0462
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 13. CACE for a meta-analysis with complete compliance information
The function cace.meta.c() performs the Bayesian hierarchical model method for meta-analysis when the dataset has complete compliance information for all studies.
set.seed(123)
out.meta.c <- cace.meta.c(data = epidural_c, conv.diag = TRUE, mcmc.samples= TRUE, study.specific = TRUE)## NA is not allowed in the input data set;
## the rows containing NA are removed.
## Start running MCMC...
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 20
## Unobserved stochastic nodes: 61
## Total graph size: 608
##
## Initializing model
##
## MCMC convergence diagnostic statistics are calculated and saved in conv.outIn this example, by calling the object smry from the output list out.meta.c, posterior estimates (posterior mean, standard deviation, posterior median, 95% credible interval, and time-series standard error) are displayed.
out.meta.c$smry## Mean SD 2.5% 50% 97.5% Naive SE Time-series SE
## CACE 0.020200 0.0627 -0.10200 0.018900 0.1490 1.14e-04 7.69e-04
## b1out 0.128000 0.0459 0.05970 0.121000 0.2370 8.39e-05 4.07e-04
## cacei[1] 0.043900 0.0679 -0.08130 0.040700 0.1870 1.24e-04 2.35e-04
## cacei[2] -0.023100 0.0489 -0.11500 -0.025000 0.0822 8.94e-05 1.89e-04
## cacei[3] -0.007630 0.0569 -0.11000 -0.011800 0.1130 1.04e-04 2.16e-04
## cacei[4] 0.065000 0.0678 -0.06620 0.064300 0.2010 1.24e-04 1.62e-04
## cacei[5] 0.054000 0.0686 -0.07380 0.051500 0.1960 1.25e-04 2.45e-04
## cacei[6] 0.026300 0.0309 -0.03390 0.026200 0.0875 5.64e-05 6.87e-05
## cacei[7] 0.002770 0.0931 -0.18900 -0.000142 0.2100 1.70e-04 3.71e-04
## cacei[8] 0.048300 0.0239 0.00171 0.048200 0.0956 4.36e-05 6.38e-05
## cacei[9] -0.010600 0.0224 -0.05520 -0.010500 0.0331 4.09e-05 5.49e-05
## cacei[10] 0.000228 0.0600 -0.12000 -0.001200 0.1280 1.10e-04 2.15e-04
## pia 0.114000 0.0742 0.02460 0.098200 0.3010 1.35e-04 3.92e-03
## pic 0.821000 0.0842 0.60500 0.838000 0.9350 1.54e-04 4.50e-03
## pin 0.064200 0.0397 0.01540 0.056700 0.1590 7.25e-05 2.11e-03
## s1out 0.184000 0.1040 0.04560 0.161000 0.4450 1.91e-04 9.03e-04
## u1out 0.127000 0.0473 0.05520 0.120000 0.2390 8.64e-05 6.10e-04
## v1out 0.107000 0.0406 0.04700 0.100000 0.2040 7.41e-05 4.55e-04Users can manually do model selection procedures by including different random effects and comparing DIC from the outputs. DIC and its two components are saved as an object DIC in the output list.
out.meta.c$DIC##
## D.bar 204.40801
## pD 44.74788
## DIC 249.155904. CACE for meta-analysis with incomplete compliance information
The function out.meta.ic() also estimates \(\theta^\text{CACE}\) using the Bayesian hierarchcal model but can accommodate studies with incomplete compliance data.
set.seed(123)
out.meta.ic <- cace.meta.ic(data = epidural_ic, conv.diag = TRUE,
mcmc.samples = TRUE, study.specific = TRUE)## NA is not allowed in the input data set;
## the rows containing NA are removed.
## Start running MCMC...
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 54
## Unobserved stochastic nodes: 146
## Total graph size: 1587
##
## Initializing model## Warning in jags.model(file = textConnection(modelstring), data = data.jags, :
## Adaptation incomplete## NOTE: Stopping adaptation
##
##
## MCMC convergence diagnostic statistics are calculated and saved in conv.outThe results are saved in the object out.meta.ic, a list containing posterior estimates for monitored parameters, DIC, convergence diagnostic statistics, and MCMC samples. In this example, the argument study.specific is TRUE, so the summary for each study-specific \(\theta^\text{CACE}_i\) is displayed in the object out.meta.ic$smry together with other parameters.
out.meta.ic$smry## Mean SD 2.5% 50% 97.5% Naive SE Time-series SE
## CACE 0.02750 0.0370 -0.04340 0.026400 0.1040 6.75e-05 3.58e-04
## b1out 0.15200 0.0418 0.08450 0.147000 0.2490 7.63e-05 6.32e-04
## cacei[1] 0.04170 0.0652 -0.07890 0.038200 0.1800 1.19e-04 2.14e-04
## cacei[2] -0.01850 0.0479 -0.10900 -0.020300 0.0839 8.75e-05 1.68e-04
## cacei[3] -0.00488 0.0544 -0.10300 -0.008340 0.1110 9.94e-05 1.89e-04
## cacei[4] 0.06500 0.0665 -0.06330 0.064100 0.1990 1.21e-04 1.59e-04
## cacei[5] 0.05250 0.0668 -0.07220 0.050300 0.1900 1.22e-04 2.34e-04
## cacei[6] 0.02660 0.0306 -0.03310 0.026400 0.0873 5.59e-05 6.92e-05
## cacei[7] 0.00889 0.0893 -0.17100 0.003940 0.2090 1.63e-04 2.83e-04
## cacei[8] 0.04760 0.0238 0.00147 0.047500 0.0947 4.34e-05 6.60e-05
## cacei[9] -0.00990 0.0224 -0.05430 -0.009770 0.0338 4.08e-05 5.68e-05
## cacei[10] 0.00311 0.0586 -0.11300 0.001200 0.1280 1.07e-04 1.76e-04
## cacei[11] 0.04100 0.0701 -0.08920 0.040700 0.1760 1.28e-04 2.62e-04
## cacei[12] 0.01930 0.0538 -0.06330 0.015100 0.1400 9.82e-05 2.71e-04
## cacei[13] -0.01300 0.0713 -0.15500 -0.013000 0.1340 1.30e-04 2.77e-04
## cacei[14] 0.06170 0.0983 -0.09470 0.046500 0.3110 1.79e-04 4.70e-04
## cacei[15] 0.05290 0.1510 -0.23300 0.037700 0.3890 2.75e-04 4.73e-04
## cacei[16] 0.00743 0.0705 -0.11200 0.000904 0.1680 1.29e-04 2.50e-04
## cacei[17] 0.02740 0.1120 -0.18200 0.019100 0.2810 2.04e-04 3.41e-04
## cacei[18] -0.01050 0.0849 -0.17800 -0.011900 0.1760 1.55e-04 2.95e-04
## cacei[19] -0.04510 0.0994 -0.23800 -0.045400 0.1640 1.81e-04 4.09e-04
## cacei[20] -0.00482 0.0930 -0.18800 -0.007450 0.2060 1.70e-04 3.29e-04
## cacei[21] 0.01150 0.1020 -0.18400 0.004750 0.2480 1.86e-04 3.64e-04
## cacei[22] 0.02270 0.1140 -0.19500 0.015600 0.2740 2.07e-04 3.55e-04
## cacei[23] 0.02190 0.1000 -0.16500 0.014900 0.2520 1.83e-04 3.25e-04
## cacei[24] 0.06770 0.1180 -0.14400 0.057800 0.3380 2.16e-04 4.13e-04
## cacei[25] -0.04780 0.1450 -0.34800 -0.039300 0.2500 2.64e-04 7.35e-04
## cacei[26] 0.08270 0.1730 -0.23400 0.059200 0.4690 3.15e-04 6.62e-04
## cacei[27] 0.22400 0.0966 0.06720 0.214000 0.4480 1.76e-04 4.49e-04
## pia 0.20400 0.1340 0.04290 0.170000 0.5730 2.45e-04 9.33e-03
## pic 0.69600 0.1440 0.32100 0.725000 0.8920 2.63e-04 1.02e-02
## pin 0.09970 0.0606 0.02120 0.088800 0.2410 1.11e-04 4.43e-03
## s1out 0.18700 0.0875 0.06050 0.172000 0.3980 1.60e-04 1.44e-03
## u1out 0.12100 0.0322 0.06680 0.118000 0.1920 5.87e-05 4.10e-04
## v1out 0.09320 0.0265 0.04930 0.090500 0.1530 4.84e-05 3.09e-04Note that when compiling the JAGS model, the warning “adaptation incomplete” may occasionally occur, indicating that the number of iterations for the adaptation process is not sufficient. The default value of n.adapt (the number of iterations for adaptation) is 1,000. This is an initial sampling phase during which the samplers adapt their behavior to maximize their efficiency (e.g., a Metropolis–Hastings random walk algorithm may change its step size) Plummer (2018). The “adaptation incomplete” warning indicates the MCMC algorithm may not achieve maximum efficiency, but it generally has little impact on the posterior estimates of the treatment effects. To avoid this warning, users may increase n.adapt.
5. Functions for plots
5-1. Diagnostic plots
The function plot.cacebaes() provides diagnostic plots for the MCMC, namely trace plots, auto-correlation plots and kernel density estimation plots. Both trace plots and auto-correlation plots can be used to examine whether the MCMC chains appear to be drawn from stationary distributions. A posterior density plot for a parameter visually shows the posterior distribution. Users can simply call this function on objects produced by cace.study(), cace.meta.c(), or cace.meta.ic(). In the example below we use the objects list obtained from fitting the Bayesian hierarchical model cace.meta.ic() to generate the three plots. To avoid lengthy output we just illustrate how these plots are produced for \(\theta^\text{CACE}\).
plot.cacebayes(obj = out.meta.ic)
## [[1]]
## NULL
##
## [[2]]
## NULL
##
## [[3]]
##
## Autocorrelations of series 'temp', by lag
##
## 0 1 2 3 4 5 6 7 8 9 10 11 12
## 1.000 0.836 0.762 0.703 0.652 0.609 0.569 0.534 0.501 0.471 0.444 0.419 0.394
## 13 14 15 16 17 18 19 20 21 22 23 24 25
## 0.372 0.352 0.333 0.314 0.297 0.282 0.266 0.253 0.240 0.228 0.216 0.205 0.193
## 26 27 28 29 30 31 32 33 34 35 36 37 38
## 0.183 0.173 0.165 0.156 0.147 0.139 0.131 0.124 0.117 0.111 0.105 0.101 0.095
## 39 40 41 42 43 44 45 46 47 48 49 50 51
## 0.091 0.087 0.083 0.078 0.073 0.069 0.065 0.062 0.059 0.056 0.054 0.051 0.048
## 52 53 54
## 0.046 0.044 0.042The trace plots show the parameter values sampled at each iteration versus the iteration number. Each chain is drawn as a separate trace plot to avoid overlay. Here we used the default n.chains = 3, so three trace plots are drawn. These plots show evidence that the posterior samples of \(\theta^\text{CACE}\) are drawn from the stationary distribution.
The density plot is smoothed using the R function density(). It shows that the kernel-smoothed posterior of \(\theta^\text{CACE}\) is almost symmetric. The posterior mean is not far from 0, indicating that the complier average causal effect of using epidural analgesia in labor on cesarean section is likely not significant.
The auto-correlation plot is a bar plot displaying the auto-correlation for different lags. At lag 0, the value of the chain has perfect auto-correlation with itself. As the lag becomes greater, the values become less correlated. After a lag of about 50, the auto-correlation drops below 0.1. If the plot shows high auto-correlation, users can run the chain longer or can choose a larger n.thin, e.g., n.thin = 10 would keep only 1 out of every 10 iterations, so that the thinned out chain is expected to have the auto-correlation dropping quickly.
5-2. Study-specific CACE in a forest plot
A graphical overview of the results can be obtained by creating a forest plot Lewis and Clarke (2001). The function plot.forest() draws a forest plot for \(\theta^{\text{CACE}}\) estimated from the meta-analysis. Users can call this function for the objects from cace.meta.c() or cace.meta.ic(). Here is an example using the object out.meta.ic:
plot.forest(data = epidural_ic, obj = out.meta.ic)
It is a forest plot of \(\theta^\text{CACE}_i\) for each study individually, using the Bayesian method with full random effects and default priors. The summary estimate based on the model cace.meta.ic() is automatically added to the figure, with the outer edges of the polygon indicating the confidence interval limits. The 95% credible interval of the summary \(\theta^{\text{CACE}}\) covers zero, indicating a non-significant complier average causal effect estimate for using epidural analgesia in labor on the risk of cesarean section for the meta-analysis with 27 trials. For a study with incomplete data on compliance status, a dashed horizontal line in the forest plot is used to represent the posterior 95% credible interval of \(\theta^\text{CACE}_i\) from the Bayesian hierarchical model fit. The study-specific \(\theta^{\text{CACE}}_i\) vary from negative to positive in individual studies, while most of the 95% credible intervals cover zero. As the \(\theta^\text{CACE}_i\) for a trial without complete compliance data is not estimable using only data from that single trial, dashed lines tend to have longer credible intervals than those with complete data (solid lines).
plot.forest(data = epidural_c, obj = out.study, obj2 = out.meta.c)
The function cace.study() can estimate CACE separately for an individual trial as long as it has complete compliance data. In that case, users can choose to generate a forest plot \(\theta^\text{CACE}_i\) for each individual study based on separate analyses. The pooled estimate of \(\theta^\text{CACE}\) and its 95% credible interval or confidence interval (the diamond in the plot) can be either from the Bayesian hierarchical model (the cace.meta.c() function) or from the two-step approach cace.study() function with the argument two.step = TRUE. The above code is an example of how to create such a plot. If obj contains the two-step meta-analysis result, the argument obj2 is optional and is included if users want to report the summary CACE estimate based on the Bayesian hierarchical model cace.meta.ic().
References
Bannister-Tyrrell, Melanie, Branko Miladinovic, Christine L Roberts, and Jane B Ford. 2015. “Adjustment for Compliance Behavior in Trials of Epidural Analgesia in Labor Using Instrumental Variable Meta-Analysis.” Journal Article. Journal of Clinical Epidemiology 68 (5): 525–33.
Gordon, Max, and Thomas Lumley. 2017. Forestplot: Advanced Forest Plot Using “Grid” Graphics. https://CRAN.R-project.org/package=forestplot.
Lewis, Steff, and Mike Clarke. 2001. “Forest Plots: Trying to See the Wood and the Trees.” BMJ 322 (7300): 1479–80.
Plummer, Martyn. 2018. Rjags: Bayesian Graphical Models Using Mcmc. https://CRAN.R-project.org/package=rjags.
Plummer, Martyn, Nicky Best, Kate Cowles, and Karen Vines. 2006. “CODA: Convergence Diagnosis and Output Analysis for Mcmc.” R News 6 (1): 7–11.