inst/doc/bayesmsm-censoring.R

In bayesmsm: Fitting Bayesian Marginal Structural Models for Longitudinal Observational Data

## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.align = 'center',
  fig.width = 9,
  fig.height = 6,
  warning = F,
  message = F
)

## ----warning = FALSE, results=FALSE, force=TRUE, include=FALSE----------------
if (!requireNamespace("bayesmsm", quietly = TRUE)) {
  stop("The package 'bayesmsm' is required to run this vignette. Please install it manually using:
       devtools::install_github('Kuan-Liu-Lab/bayesmsm')")
} else {
  library(bayesmsm)
}

## -----------------------------------------------------------------------------
# 1) Define coefficient lists for 2 visits
amodel <- list(
  # Visit 1: logit P(A1=1) = -0.3 + 0.4*L1_1 - 0.2*L2_1
  c("(Intercept)" = -0.3, "L1_1" = 0.4, "L2_1" = -0.2),
  # Visit 2: logit P(A2=1) = -0.1 + 0.3*L1_2 - 0.1*L2_2 + 0.5*A_prev
  c("(Intercept)" = -0.1, "L1_2" = 0.3, "L2_2" = -0.1, "A_prev" = 0.5)
)

# 2) Define outcome model: logistic on treatments and last covariates
ymodel <- c(
  "(Intercept)" = -0.8,
  "A1"          = 0.2,
  "A2"          = 0.4,
  "L1_2"        = 0.3,
  "L2_2"        = -0.3
)

# 3) Define right-censoring models at each visit
cmodel <- list(
  # Censor at visit 1 based on baseline covariates and A1
  c("(Intercept)" = -1.5, "L1_1" = 0.2, "L2_1" = -0.2, "A" = 0.2),
  # Censor at visit 2 based on visit-2 covariates and A2
  c("(Intercept)" = -1.5, "L1_2" = 0.1, "L2_2" = -0.1, "A" = 0.3)
)

# 4) Load package and simulate data
simdat_cen <- simData(
  n                = 100,
  n_visits         = 2,
  covariate_counts = c(2, 2),
  amodel           = amodel,
  ymodel           = ymodel,
  y_type           = "binary",
  right_censor     = TRUE,
  cmodel           = cmodel,
  seed             = 123
)

# 5) Inspect first rows
head(simdat_cen)

## -----------------------------------------------------------------------------
weights_cen <- bayesweight_cen(
  trtmodel.list = list(
    A1 ~ L1_1 + L2_1,
    A2 ~ L1_2 + L2_2 + A1),
  cenmodel.list = list(
    C1 ~ L1_1 + L2_1 + A1,
    C2 ~ L1_2 + L2_2 + A2),
  data = simdat_cen,
  n.chains = 1,
  n.iter = 200,
  n.burnin = 100,
  n.thin = 1,
  seed = 890123,
  parallel = FALSE)

summary(weights_cen$weights)

## -----------------------------------------------------------------------------
cat(weights_cen$model_string)

## -----------------------------------------------------------------------------
# Remove all NAs (censored observations) from the original dataset and weights
simdat_cen$weights <- weights_cen$weights
simdat_cen2 <- na.omit(simdat_cen)

model <- bayesmsm(ymodel = Y ~ A1 + A2,
  nvisit = 2,
  reference = c(rep(0,2)),
  comparator = c(rep(1,2)),
  family = "binomial",
  data = simdat_cen2,
  wmean = simdat_cen2$weights,
  nboot = 50,
  optim_method = "BFGS",
  parallel = TRUE,
  seed = 890123,
  ncore = 2)
str(model)

## -----------------------------------------------------------------------------
summary_bayesmsm(model)

## -----------------------------------------------------------------------------
plot_ATE(model)

## -----------------------------------------------------------------------------
plot_APO(model, effect_type = "effect_comparator")
plot_APO(model, effect_type = "effect_reference")

## -----------------------------------------------------------------------------
plot_est_box(model)