Bug fix and function testing updated

1. Fixed various bugs, including deleting the text accompanied with compiling the package website. For example:
<<<<<<< HEAD
=======
>>>>>>> c02a61c,
2. Changed the DESCRIPTION, NAMESPACE files to reflect correct imports
3. Fixed a warning message in function bayesweight when running parallel computation
4. Fixed function testing for functions bayesweight, bayesweight_cen
5. Trimmed example code margin in each function to avoid warning in function check
6. Updated bayesweight_cen and correctly documented the function
7. Remaining problem: should have used package::function() instead of library() or require() (this will lead to a warning message)

DESCRIPTION

@@ -4,22 +4,10 @@ Title: Bayesian Marginal Structural Models
 Version: 0.1.0
 Author: Kuan Liu, Xiao Yan
 Maintainer: Xiao Yan <[email protected]>
-<<<<<<< HEAD
 Depends:
     R (>= 4.2.0)
-Imports:
-	MCMCpack,
-	doParallel,
-	foreach,
-	parallel,
-	R2jags,
-	coda (>= 0.19-4)
-License: MIT + file LICENSE
-Encoding: UTF-8
-LazyData: true
-RoxygenNote: 7.3.2
 Suggests: 
-=======
+    testthat (>= 3.0.0)
 License: MIT + file LICENSE
 Encoding: UTF-8
 LazyData: true
@@ -31,9 +19,7 @@ Imports:
 	  parallel,
 	  R2jags,
 	  knitr,
-	  coda (>= 0.19-4),
->>>>>>> c02a61c92d3be0eeca4a22bb0196a74db590a6ec
-    testthat (>= 3.0.0)
+	  coda (>= 0.19-4)
 NeedsCompilation: no
 Config/testthat/edition: 3
 Description: This package provides tools for estimating causal effects using Bayesian Marginal Structural Models. It includes functions for estimating Bayesian weights using JAGS and for Bayesian non-parametric bootstrap to calculate causal effects.

NAMESPACE

@@ -2,15 +2,15 @@
 
 export(bayesmsm)
 export(bayesweight)
+export(bayesweight_cen)
 export(plot_APO)
 export(plot_ATE)
 export(plot_est_box)
-<<<<<<< HEAD
 export(summary_bayesmsm)
 import(MCMCpack)
-=======
->>>>>>> c02a61c92d3be0eeca4a22bb0196a74db590a6ec
 import(doParallel)
+import(foreach)
+import(parallel)
 importFrom(MCMCpack,rdirichlet)
 importFrom(R2jags,jags)
 importFrom(coda,mcmc)

R/bayesmsm.R

@@ -29,19 +29,18 @@
 #'
 #' @importFrom foreach "%dopar%"
 #' @import doParallel
-<<<<<<< HEAD
 #' @import parallel
 #' @import MCMCpack
-=======
 #' @importFrom MCMCpack rdirichlet
->>>>>>> c02a61c92d3be0eeca4a22bb0196a74db590a6ec
 #'
 #' @export
 #'
 #' @examples
 #'
 #' # Continuous outcome
-#' testdata <- read.csv(system.file("extdata", "continuous_outcome_data.csv", package = "bayesmsm"))
+#' testdata <- read.csv(system.file("extdata",
+#'                                  "continuous_outcome_data.csv",
+#'                                  package = "bayesmsm"))
 #' model <- bayesmsm(ymodel = y ~ a_1+a_2,
 #'                            nvisit = 2,
 #'                            reference = c(rep(0,2)),
@@ -71,12 +70,10 @@ bayesmsm <- function(ymodel,
                      ncore = 6){
 
   # load all the required R packages;
-<<<<<<< HEAD
   # require(foreach)
   # require(doParallel)
   # require(MCMCpack)
   # require(parallel)
-=======
   # if (!require(foreach)){
   #   install.packages("foreach",repos="http://cran.r-project.org")
   #   library(foreach)
@@ -93,7 +90,6 @@ bayesmsm <- function(ymodel,
   require(foreach)
   require(doParallel)
   require(MCMCpack)
->>>>>>> c02a61c92d3be0eeca4a22bb0196a74db590a6ec
 
   # return error message if the input weight vector has different length comparing to the outcome Y;
   if (length(wmean) != nrow(data)) {

R/bayesweight.R

@@ -16,23 +16,22 @@
 #'
 #' @importFrom R2jags jags
 #' @importFrom coda mcmc
-<<<<<<< HEAD
 #' @import parallel
 #' @import doParallel
 #' @import foreach
-=======
-#' @import doParallel
 #' @importFrom foreach "%dopar%"
 #'
->>>>>>> c02a61c92d3be0eeca4a22bb0196a74db590a6ec
 #' @export
 #'
 #' @examples
 #'
 #' # Continuous outcome
-#' testdata <- read.csv(system.file("extdata", "continuous_outcome_data.csv", package = "bayesmsm"))
+#' testdata <- read.csv(system.file("extdata",
+#'                                  "continuous_outcome_data.csv",
+#'                                  package = "bayesmsm"))
 #' weights <- bayesweight(trtmodel.list = list(a_1 ~ w1 + w2 + L1_1 + L2_1,
-#'                                             a_2 ~ w1 + w2 + L1_1 + L2_1 + L1_2 + L2_2 + a_1),
+#'                                             a_2 ~ w1 + w2 + L1_1 + L2_1 +
+#'                                                   L1_2 + L2_2 + a_1),
 #'                        data = testdata,
 #'                        n.iter = 2500,
 #'                        n.burnin = 1500,
@@ -53,14 +52,12 @@ bayesweight <- function(trtmodel.list,
                         parallel = TRUE){
 
   # Load all the required R packages;
-<<<<<<< HEAD
   # require(foreach)
   # require(doParallel)
   # require(MCMCpack)
   # require(parallel)
   # require(R2jags)
   # require(coda)
-=======
   # if (!require(R2jags)){
   #   install.packages("R2jags",repos="http://cran.r-project.org")
   #   library(R2jags)
@@ -78,7 +75,6 @@ bayesweight <- function(trtmodel.list,
   require(coda)
   require(doParallel)
   require(foreach)
->>>>>>> c02a61c92d3be0eeca4a22bb0196a74db590a6ec
 
   create_marginal_treatment_models <- function(trtmodel.list) {
     # Initialize the list for the marginal treatment models
@@ -292,15 +288,21 @@ bayesweight <- function(trtmodel.list,
     if (n.chains >= available_cores) {
       stop(paste("Parallel MCMC requires 1 core per chain. You have", available_cores, "cores. We recommend using", available_cores - 2, "cores."))
     }
-<<<<<<< HEAD
     # Run JAGS model in parallel
     cl <- parallel::makeCluster(n.chains)
     doParallel::registerDoParallel(cl)
-=======
     # Run JAGS model in parallel;
-    cl <- makeCluster(n.chains)
-    registerDoParallel(cl)
->>>>>>> c02a61c92d3be0eeca4a22bb0196a74db590a6ec
+    # cl <- makeCluster(n.chains)
+    # registerDoParallel(cl)
+
+    # Ensure the cluster is stopped when the function exits, even in case of error
+    on.exit({
+      if (!is.null(cl)) {
+        parallel::stopCluster(cl)
+        foreach::registerDoSEQ()  # Reset to sequential processing correctly
+      }
+    }, add = TRUE)
+
     jags.model.wd <- paste(getwd(), '/treatment_model.txt',sep='')
 
     posterior <- foreach::foreach(i=1:n.chains, .packages=c('R2jags'),
@@ -319,7 +321,8 @@ bayesweight <- function(trtmodel.list,
                            return(do.call(rbind, lapply(out.mcmc, as.matrix)))
 
                          }
-    parallel::stopCluster(cl)
+    # parallel::stopCluster(cl)
+    # doParallel::registerDoSEQ()
 
   } else if (parallel == FALSE) {
 

R/bayesweight_cen.R

@@ -17,13 +17,19 @@
 #'
 #' @examples
 #'
-#' simdat_cen <- read.csv(system.file("extdata", "sim_causal.csv", package = "bayesmsm"))
+#' simdat_cen <- read.csv(system.file("extdata",
+#'                                    "sim_causal.csv",
+#'                                    package = "bayesmsm"))
 #' weights_cen <- bayesweight_cen(trtmodel.list = list(A1 ~ L11 + L21,
-#'                                                     A2 ~ L11 + L21 + L12 + L22 + A1,
-#'                                                     A3 ~ L11 + L21 + L12 + L22 + A1 + L13 + L23 + A2),
+#'                                                     A2 ~ L11 + L21 + L12 +
+#'                                                          L22 + A1,
+#'                                                     A3 ~ L11 + L21 + L12 +
+#'                                                          L22 + A1 + L13 +
+#'                                                          L23 + A2),
 #'                                cenmodel.list = list(C1 ~ L11 + L21,
 #'                                                     C2 ~ L11 + L21 + A1,
-#'                                                     C3 ~ L11 + L21 + A1 + L12 + L22 + A2),
+#'                                                     C3 ~ L11 + L21 + A1 +
+#'                                                          L12 + L22 + A2),
 #'                                data = simdat_cen,
 #'                                n.iter = 25000,
 #'                                n.burnin = 15000,

R/plot_APO.R

@@ -10,7 +10,9 @@
 #' @export
 #'
 #' @examples
-#' testdata <- read.csv(system.file("extdata", "continuous_outcome_data.csv", package = "bayesmsm"))
+#' testdata <- read.csv(system.file("extdata",
+#'                                  "continuous_outcome_data.csv",
+#'                                  package = "bayesmsm"))
 #' model <- bayesmsm(ymodel = y ~ a_1+a_2,
 #'                            nvisit = 2,
 #'                            reference = c(rep(0,2)),

R/plot_ATE.R

@@ -14,7 +14,9 @@
 #' @export
 #'
 #' @examples
-#' testdata <- read.csv(system.file("extdata", "continuous_outcome_data.csv", package = "bayesmsm"))
+#' testdata <- read.csv(system.file("extdata",
+#'                                  "continuous_outcome_data.csv",
+#'                                  package = "bayesmsm"))
 #' model <- bayesmsm(ymodel = y ~ a_1+a_2,
 #'                            nvisit = 2,
 #'                            reference = c(rep(0,2)),

R/plot_est_box.R

@@ -7,7 +7,9 @@
 #' @export
 #'
 #' @examples
-#' testdata <- read.csv(system.file("extdata", "continuous_outcome_data.csv", package = "bayesmsm"))
+#' testdata <- read.csv(system.file("extdata",
+#'                                  "continuous_outcome_data.csv",
+#'                                  package = "bayesmsm"))
 #' model <- bayesmsm(ymodel = y ~ a_1+a_2,
 #'                            nvisit = 2,
 #'                            reference = c(rep(0,2)),
@@ -19,7 +21,8 @@
 #'                            optim_method = "BFGS",
 #'                            parallel = FALSE,
 #'                            ncore = 2)
-#' plot_est_box(model$bootdata) # without reference & comparator information below labels
+#' plot_est_box(model$bootdata) # without reference & comparator information
+#'                              # below labels
 #' plot_est_box(model) # with reference & comparator information below labels
 #'
 plot_est_box <- function(input, ...) {

README.md

@@ -94,5 +94,5 @@ If you use the `bayesmsm` package in your research, please cite the reference pa
 
 # Developers
 
--  Kuan Liu
 -  Xiao Yan
+-  Kuan Liu

censoring_model.txt

@@ -0,0 +1,93 @@
+model{
+
+for (i in 1:N1) {
+
+# conditional model;
+A1[i] ~ dbern(p1[i])
+logit(p1[i]) <- b10 + b11*L11[i] + b12*L21[i]
+C1[i] ~ dbern(cp1[i])
+logit(cp1[i]) <- s10 + s11*L11[i] + s12*L21[i]
+
+# marginal model;
+A1s[i] ~ dbern(p1s[i])
+logit(p1s[i]) <- bs10
+C1s[i] ~ dbern(cp1s[i])
+logit(cp1s[i]) <- ts10
+}
+
+for (i in 1:N2) {
+
+# conditional model;
+A2[i] ~ dbern(p2[i])
+logit(p2[i]) <- b20 + b21*L11[i] + b22*L21[i] + b23*L12[i] + b24*L22[i] + b25*A1[i]
+C2[i] ~ dbern(cp2[i])
+logit(cp2[i]) <- s20 + s21*L11[i] + s22*L21[i] + s23*A1[i]
+
+# marginal model;
+A2s[i] ~ dbern(p2s[i])
+logit(p2s[i]) <- bs20 + bs21*A1s[i]
+C2s[i] ~ dbern(cp2s[i])
+logit(cp2s[i]) <- ts20 + ts21*A1s[i]
+}
+
+for (i in 1:N3) {
+
+# conditional model;
+A3[i] ~ dbern(p3[i])
+logit(p3[i]) <- b30 + b31*L11[i] + b32*L21[i] + b33*L12[i] + b34*L22[i] + b35*A1[i] + b36*L13[i] + b37*L23[i] + b38*A2[i]
+C3[i] ~ dbern(cp3[i])
+logit(cp3[i]) <- s30 + s31*L11[i] + s32*L21[i] + s33*A1[i] + s34*L12[i] + s35*L22[i] + s36*A2[i]
+
+# marginal model;
+A3s[i] ~ dbern(p3s[i])
+logit(p3s[i]) <- bs30 + bs31*A1s[i] + bs32*A2s[i]
+C3s[i] ~ dbern(cp3s[i])
+logit(cp3s[i]) <- ts30 + ts31*A1s[i] + ts32*A2s[i]
+}
+
+# Priors
+b10 ~ dunif(-10, 10)
+b11 ~ dunif(-10, 10)
+b12 ~ dunif(-10, 10)
+s10 ~ dunif(-10, 10)
+s11 ~ dunif(-10, 10)
+s12 ~ dunif(-10, 10)
+bs10 ~ dunif(-10, 10)
+ts10 ~ dunif(-10, 10)
+b20 ~ dunif(-10, 10)
+b21 ~ dunif(-10, 10)
+b22 ~ dunif(-10, 10)
+b23 ~ dunif(-10, 10)
+b24 ~ dunif(-10, 10)
+b25 ~ dunif(-10, 10)
+s20 ~ dunif(-10, 10)
+s21 ~ dunif(-10, 10)
+s22 ~ dunif(-10, 10)
+s23 ~ dunif(-10, 10)
+bs20 ~ dunif(-10, 10)
+bs21 ~ dunif(-10, 10)
+ts20 ~ dunif(-10, 10)
+ts21 ~ dunif(-10, 10)
+b30 ~ dunif(-10, 10)
+b31 ~ dunif(-10, 10)
+b32 ~ dunif(-10, 10)
+b33 ~ dunif(-10, 10)
+b34 ~ dunif(-10, 10)
+b35 ~ dunif(-10, 10)
+b36 ~ dunif(-10, 10)
+b37 ~ dunif(-10, 10)
+b38 ~ dunif(-10, 10)
+s30 ~ dunif(-10, 10)
+s31 ~ dunif(-10, 10)
+s32 ~ dunif(-10, 10)
+s33 ~ dunif(-10, 10)
+s34 ~ dunif(-10, 10)
+s35 ~ dunif(-10, 10)
+s36 ~ dunif(-10, 10)
+bs30 ~ dunif(-10, 10)
+bs31 ~ dunif(-10, 10)
+bs32 ~ dunif(-10, 10)
+ts30 ~ dunif(-10, 10)
+ts31 ~ dunif(-10, 10)
+ts32 ~ dunif(-10, 10)
+}