diff --git a/DESCRIPTION b/DESCRIPTION index 9a620b9..1216152 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -4,22 +4,10 @@ Title: Bayesian Marginal Structural Models Version: 0.1.0 Author: Kuan Liu, Xiao Yan Maintainer: Xiao Yan -<<<<<<< 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. diff --git a/NAMESPACE b/NAMESPACE index bf05bb1..a919705 100644 --- a/NAMESPACE +++ b/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) diff --git a/R/bayesmsm.R b/R/bayesmsm.R index bfee22e..daa62ba 100644 --- a/R/bayesmsm.R +++ b/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)) { diff --git a/R/bayesweight.R b/R/bayesweight.R index b7865da..315cba0 100644 --- a/R/bayesweight.R +++ b/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) { diff --git a/R/bayesweight_cen.R b/R/bayesweight_cen.R index 33de10c..49983e3 100644 --- a/R/bayesweight_cen.R +++ b/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, diff --git a/R/plot_APO.R b/R/plot_APO.R index daeeb4b..3b07624 100644 --- a/R/plot_APO.R +++ b/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)), diff --git a/R/plot_ATE.R b/R/plot_ATE.R index 6dfc5f2..e37507e 100644 --- a/R/plot_ATE.R +++ b/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)), diff --git a/R/plot_est_box.R b/R/plot_est_box.R index f251a1c..681a39e 100644 --- a/R/plot_est_box.R +++ b/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, ...) { diff --git a/README.md b/README.md index 33ddcd9..7bd9f6d 100644 --- a/README.md +++ b/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 diff --git a/censoring_model.txt b/censoring_model.txt new file mode 100644 index 0000000..6bb46c9 --- /dev/null +++ b/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) +} diff --git a/man/bayesmsm.Rd b/man/bayesmsm.Rd index c2e61be..9dc1723 100644 --- a/man/bayesmsm.Rd +++ b/man/bayesmsm.Rd @@ -63,7 +63,9 @@ effects in Bayesian marginal structural models. It supports both continuous \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)), diff --git a/man/bayesweight.Rd b/man/bayesweight.Rd index 87c378e..6f47153 100644 --- a/man/bayesweight.Rd +++ b/man/bayesweight.Rd @@ -42,9 +42,12 @@ It uses JAGS for Bayesian inference and supports parallel computation to speed u \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, diff --git a/man/bayesweight_cen.Rd b/man/bayesweight_cen.Rd new file mode 100644 index 0000000..d9ea820 --- /dev/null +++ b/man/bayesweight_cen.Rd @@ -0,0 +1,71 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/bayesweight_cen.R +\name{bayesweight_cen} +\alias{bayesweight_cen} +\title{Bayesian Weight Estimation for Censored Data} +\usage{ +bayesweight_cen( + trtmodel.list = list(A1 ~ L11 + L21, 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), + data, + n.iter = 25000, + n.burnin = 15000, + n.thin = 5, + parallel = FALSE, + n.chains = 1, + seed = 890123 +) +} +\arguments{ +\item{trtmodel.list}{A list of formulas corresponding to each time point with the time-specific treatment variable on the left-hand side and pre-treatment covariates to be balanced on the right-hand side. Interactions and functions of covariates are allowed.} + +\item{cenmodel.list}{A list of formulas for the censoring data at each time point, with censoring indicators on the left-hand side and covariates on the right-hand side.} + +\item{data}{A data frame containing the variables in the models (treatment, censoring, and covariates).} + +\item{n.iter}{Number of iterations to run the MCMC algorithm in JAGS.} + +\item{n.burnin}{Number of iterations to discard as burn-in in the MCMC algorithm.} + +\item{n.thin}{Thinning rate for the MCMC samples.} + +\item{parallel}{Logical, indicating whether to run the MCMC sampling in parallel (default is `FALSE`).} + +\item{n.chains}{Number of MCMC chains to run. If parallel is `TRUE`, this specifies the number of chains run in parallel.} + +\item{seed}{A seed for random number generation to ensure reproducibility of the MCMC.} +} +\value{ +A vector of posterior mean weights, computed by taking the average of the weights across all MCMC iterations. +} +\description{ +This function computes posterior mean weights using Bayesian estimation for treatment models and censoring models across multiple time points. The models can be run in parallel to estimate the weights needed for causal inference with censored data. +} +\examples{ + +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), + cenmodel.list = list(C1 ~ L11 + L21, + C2 ~ L11 + L21 + A1, + C3 ~ L11 + L21 + A1 + + L12 + L22 + A2), + data = simdat_cen, + n.iter = 25000, + n.burnin = 15000, + n.thin = 5, + parallel = FALSE, + n.chains = 1, + seed = 890123) + + + +} diff --git a/man/plot_APO.Rd b/man/plot_APO.Rd index 5b3318d..e3e8d1f 100644 --- a/man/plot_APO.Rd +++ b/man/plot_APO.Rd @@ -20,7 +20,9 @@ A density plot showing the distribution of the specified average potential outco This function plots the density of APO for a specified effect type from bootstrap simulation results. } \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)), diff --git a/man/plot_ATE.Rd b/man/plot_ATE.Rd index c26e039..b2316c1 100644 --- a/man/plot_ATE.Rd +++ b/man/plot_ATE.Rd @@ -42,7 +42,9 @@ plot_ATE( Plot Average Treatment Effect Density from Bootstrap Results } \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)), diff --git a/man/plot_est_box.Rd b/man/plot_est_box.Rd index a17d769..f718a82 100644 --- a/man/plot_est_box.Rd +++ b/man/plot_est_box.Rd @@ -18,7 +18,9 @@ An error bar plot of the mean effects and their 95\% confidence intervals for co Error bar plots for treatment effects } \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)), @@ -30,7 +32,8 @@ model <- bayesmsm(ymodel = y ~ a_1+a_2, 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 } diff --git a/tests/testthat/test-bayesmsm.R b/tests/testthat/test-bayesmsm.R index fe286f3..9da32a8 100644 --- a/tests/testthat/test-bayesmsm.R +++ b/tests/testthat/test-bayesmsm.R @@ -14,10 +14,7 @@ test_that("bayesmsm works with no errors", { wmean = rep(1, 1000), nboot = 100, optim_method = "BFGS", -<<<<<<< HEAD -======= # estimand = 'RD', ->>>>>>> c02a61c92d3be0eeca4a22bb0196a74db590a6ec seed = 890123, parallel = TRUE, ncore = 2) @@ -25,8 +22,6 @@ test_that("bayesmsm works with no errors", { # Check if 'model' is a list expect_true(is.list(model)) -<<<<<<< HEAD -======= # Check if 'model' has 6 elements expect_length(model, 6) @@ -43,8 +38,6 @@ test_that("bayesmsm works with no errors", { # expect_type(model$reference, "double") # expect_type(model$comparator, "double") - ->>>>>>> c02a61c92d3be0eeca4a22bb0196a74db590a6ec testdata2 <- read.csv(system.file("extdata", "binary_outcome_data.csv", package = "bayesmsm")) model2 <- bayesmsm(ymodel = y ~ a_1+a_2, nvisit = 2, @@ -55,10 +48,7 @@ test_that("bayesmsm works with no errors", { wmean = rep(1, 1000), nboot = 100, optim_method = "BFGS", -<<<<<<< HEAD -======= # estimand = 'OR', ->>>>>>> c02a61c92d3be0eeca4a22bb0196a74db590a6ec seed = 890123, parallel = TRUE, ncore = 2) @@ -66,8 +56,6 @@ test_that("bayesmsm works with no errors", { # Check if 'model' is a list expect_true(is.list(model2)) -<<<<<<< HEAD -======= # Check if 'model' has 12 elements expect_length(model2, 12) @@ -84,6 +72,5 @@ test_that("bayesmsm works with no errors", { # expect_true(is.data.frame(model2$bootdata)) # expect_type(model2$reference, "double") # expect_type(model2$comparator, "double") ->>>>>>> c02a61c92d3be0eeca4a22bb0196a74db590a6ec }) diff --git a/tests/testthat/test-bayesweight.R b/tests/testthat/test-bayesweight.R index bd2c82a..7ec8a25 100644 --- a/tests/testthat/test-bayesweight.R +++ b/tests/testthat/test-bayesweight.R @@ -15,9 +15,6 @@ test_that("bayesweight works with no errors", { seed = 890123, parallel = TRUE) - # Check that the weights object is created and is of the expected class - expect_s3_class(weights, "mcmc.list") - # Check that the weights object has the expected length (equal to the number of observations in testdata) expect_equal(length(weights), nrow(testdata)) diff --git a/tests/testthat/test-bayesweight_cen.R b/tests/testthat/test-bayesweight_cen.R index 4a5abe2..ea67b14 100644 --- a/tests/testthat/test-bayesweight_cen.R +++ b/tests/testthat/test-bayesweight_cen.R @@ -17,9 +17,6 @@ test_that("bayesweight_cen works with no errors", { n.chains = 1, seed = 890123) - # Check that the weights_cen object is created and is of the expected class - expect_s3_class(weights_cen, "mcmc.list") - # Check that the weights_cen object has the expected dimensions (length or rows matching the data size) expect_equal(length(weights_cen), nrow(simdat_cen)) @@ -27,6 +24,6 @@ test_that("bayesweight_cen works with no errors", { expect_true(is.numeric(weights_cen)) # Check that weights are non-negative - expect_true(all(weights_cen >= 0)) + expect_true(all(na.omit(weights_cen) >= 0)) }) diff --git a/treatment_model.txt b/treatment_model.txt new file mode 100644 index 0000000..710e5df --- /dev/null +++ b/treatment_model.txt @@ -0,0 +1,44 @@ +model{ +#N = nobs +for(i in 1:N){ + +# visit 1; +# marginal treatment assignment model, visit 1; +a_1s[i] ~ dbern(pa_1s[i]) +pa_1s[i] <- ilogit(bs10) + +# conditional treatment assignment model, visit 1; +a_1[i] ~ dbern(pa_1[i]) +pa_1[i] <- ilogit(b10 + b11*w1[i] + b12*w2[i] + b13*L1_1[i] + b14*L2_1[i]) + +# visit 2; +# marginal treatment assignment model, visit 2; +a_2s[i] ~ dbern(pa_2s[i]) +pa_2s[i] <- ilogit(bs20 + bs21*a_1s[i]) + +# conditional treatment assignment model, visit 2; +a_2[i] ~ dbern(pa_2[i]) +pa_2[i] <- ilogit(b20 + b21*w1[i] + b22*w2[i] + b23*L1_1[i] + b24*L2_1[i] + b25*L1_2[i] + b26*L2_2[i] + b27*a_1[i]) + +# export quantity in full posterior specification; +w[i] <- (pa_1s[i]*pa_2s[i])/(pa_1[i]*pa_2[i]) +} + +#prior; +bs10~dnorm(0,.01) +b10~dnorm(0,.01) +b11~dnorm(0,.01) +b12~dnorm(0,.01) +b13~dnorm(0,.01) +b14~dnorm(0,.01) +bs20~dnorm(0,.01) +bs21~dnorm(0,.01) +b20~dnorm(0,.01) +b21~dnorm(0,.01) +b22~dnorm(0,.01) +b23~dnorm(0,.01) +b24~dnorm(0,.01) +b25~dnorm(0,.01) +b26~dnorm(0,.01) +b27~dnorm(0,.01) +}