-=======
-
-
-
-
-
-
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.
-
-
Reference
-
-
Reference paper on Bayesian marginal structural models:
-
-Saarela, O., Stephens, D. A., Moodie, E. E., & Klein, M. B. (2015). On Bayesian estimation of marginal structural models. Biometrics, 71(2), 279-288.
Liu, K., Saarela, O., Feldman, B. M., & Pullenayegum, E. (2020). Estimation of causal effects with repeatedly measured outcomes in a Bayesian framework. Statistical methods in medical research, 29(9), 2507-2519.
-
$ About
-
-
Installation
-
-
Install using devtools package:
-
-
-
Dependency
-
-
This package depends on the following R packages: - MCMCpack - doParallel - foreach - parallel - R2jags - coda
-
-
Examples
-
-
Here are some examples demonstrating how to use the bayesmsm package:
-
-# Load example data
-testdata <- read.csv(system.file("extdata", "continuous_outcome_data.csv", package = "bayesmsm"))
-
-# Calculate Bayesian weights
-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
- ),
- data = testdata,
- n.iter = 250,
- n.burnin = 150,
- n.thin = 5,
- n.chains = 2,
- seed = 890123,
- parallel = TRUE
-)
-
-# Perform Bayesian non-parametric bootstrap
-model <- bayesmsm(
- ymodel = y ~ a_1 + a_2,
- nvisit = 2,
- reference = c(rep(0, 2)),
- comparator = c(rep(1, 2)),
- family = "gaussian",
- data = testdata,
- wmean = weights,
- nboot = 1000,
- optim_method = "BFGS",
- seed = 890123,
- parallel = TRUE,
- ncore = 2
-)
-
-# View model summary
-summary.bayesmsm(model)
-
-
Documentation
-
-
For comprehensive documentation and examples, visit the pkgdown website.
-
-
License
-
-
This package is licensed under the MIT License. See the LICENSE file for details.
-
-
Citation
-
-
If you use the bayesmsm package in your research, please cite the reference papers listed above.
-
-
Developers
-
-
-- Kuan Liu
-- Xiao Yan (Maintainer)
-
-
-
+
1. Introduction
+
+The bayesmsm package enables easy implementation of
+the Bayesian marginal structural models (BMSMs) for longitudinal data.
+The methodology of BMSMs can be divided into 2 estimation steps:
+
+- Step 1. Bayesian treatment effect weight estimation
+- Step 2. Bayesian non-parametric bootstrap to maximize the utility
+function with respect to the causal effect
+
For Step 1, we estimate treatment weights \(w_{ij}\) using posterior samples of the
+\(\alpha\) and \(\beta\) via fitting a series of logistic
+regressions in a Bayesian framework. The package incorporates both
+Inverse Probability of Treatment Weighting (IPTW) and Inverse
+Probability of Censoring Weighting (IPCW) to handle longitudinal data
+without and with right-censoring. For Step 2, \(P_n(v_{ij})\) is estimated via
+non-parametric Bayesian bootstrap with \(Dir(1,\ldots,1)\) sampling
+weights.
The main functions in this package include:
+
+bayesweight: Calculates Bayesian weights for
+subject-specific treatment effects.bayesmsm: Estimates marginal structural models using
+the calculated Bayesian weights.plot_ATE: Plots the estimated Average Treatment Effect
+(ATE).plot_APO: Plots the estimated Average Potential Outcome
+(APO).plot_est_box: Plots the distribution of estimated
+treatment effects.
Installation
+
+- To install the bayesmsm package, you can use the
+
devtools package to install it directly from GitHub:
+
+
devtools::install_github("Kuan-Liu-Lab/bayesmsm")
+library(bayesmsm)
+
2. Simulated observational data with a time-dependent
+treatment
+
+- The simulated dataset (continuous outcome)
+
+- 1000 patients and 3 visits (2 of which patients were assigned a
+treatment)
+- y, an end-of-study continuous outcome
+- z, a binary treatment
+- w1 and w2 are two baseline covariates (one continuous and one
+binary, mimicking age and sex)
+- L1 and L2 are two time-dependent covariates (also one continuous and
+one binary)
+- no missing data
+
+- The simulated DAG
+
+
library(DiagrammeR)
+grViz("
+ digraph causal {
+ # Nodes
+ node [shape=plaintext]
+ W [label = 'w1, w2']
+ L1 [label = 'L11, L21']
+ Z1 [label = 'Z1']
+ L2 [label = 'L12, L22']
+ Z2 [label = 'Z2']
+ Y [label = 'Y']
+
+ # Edges
+ edge [color=black, arrowhead=vee]
+ rankdir = LR
+ W->L1
+ W->Z1
+ W->L2
+ W->Z2
+ W->Y
+ L1->Z1
+ L1->L2
+ L1->Z2
+ L1->Y
+ Z1->L2
+ Z1->Z2
+ Z1->Y
+ L2->Z2
+ L2->Y
+ Z2->Y
+
+ # Graph
+ graph [overlap=true, fontsize=14]
+ }")
+
+
+
library(DT)
+options(scipen = 999)
+
+testdata <- read.csv(system.file("extdata", "continuous_outcome_data.csv", package = "bayesmsm"))
+
+# look at the data;
+datatable(testdata,
+ rownames = FALSE,
+ options = list(dom = 't')) %>%
+ formatRound(columns=c('w2', 'L2_1', 'L2_2', 'y'), digits=2)
+
+
+
+- Frequency Counts by Treatment Combinations
+
+
# frequency counts by treatment combinations;
+table(testdata$a_1, testdata$a_2)
+#>
+#> 0 1
+#> 0 520 201
+#> 1 111 168
+
+- Suppose the causal parameter of interest is the average treatment
+effect between always treated and never treated,
+
+
\[
+ATE = E(Y \mid Z_1 = 1, Z_2 = 1) - E(Y \mid Z_1 = 0, Z_2 = 0)
+\]
+
3. Bayesian treatment effect weight estimation using
+bayesweight
+
+- The following code calls the function
bayesweight to
+run JAGS and calculate the weights.
+
+- Non-parallel computing requires that
n.chains = 1.
+Parallel MCMC requires at least 2 chains because computing is running on
+1 core per chain, and we recommend using at most 2 chains less than the
+number of available cores on your computer.
+- Running this function automatically saves a JAGS model file in the
+working directory, which the user can check to review the model
+specifications.
+
+- Parameters Description:
+
+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.data: The dataset containing all the variables
+specified in trtmodel.list.n.iter: Total number of iterations for each chain
+(including burn-in).n.burnin: Number of iterations to discard at the
+beginning of the simulation (burn-in).n.thin: Thinning rate for the MCMC sampler.n.chains: Number of MCMC chains to run. For
+non-parallel execution, this should be set to 1. For parallel execution,
+it requires at least 2 chains.seed: Seed to ensure reproducibility.parallel: Logical flag indicating whether to run the
+MCMC chains in parallel. Default is TRUE.
+
+
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),
+ data = testdata,
+ n.iter = 10000,
+ n.burnin = 5000,
+ n.thin = 5,
+ n.chains = 4,
+ seed = 890123,
+ parallel = TRUE)
+str(weights)
+#> num [1:1000] 1.233 1.138 1.016 0.83 0.794 ...
+
+- It returns a list containing:
+
+weights: The calculated weights for subject-specific
+treatment effects.
+
+
4. Bayesian non-parametric bootstrap to maximize the utility
+function with respect to the causal effect using
+bayesmsm
+
The function bayesmsm estimates causal effect of
+time-varying treatments. It uses subject-specific treatment assignmennt
+weights weights calculated using bayesweight, and
+performs Bayesian non-parametric bootstrap to estimate the causal
+parameters.
+
+- Parameters Description:
+
+ymodel: A formula representing the outcome model, which
+can include interactions and functions of covariates.nvisit: Specifies the number of visits or time points
+considered in the model.reference: The baseline or reference intervention
+across all visits, typically represented by a vector of zeros indicating
+no treatment (default is a vector of all zeros).comparator: The comparison intervention across all
+visits, typically represented by a vector of ones indicating full
+treatment (default is a vector of all ones).family: Specifies the outcome distribution family; use
+“gaussian” for continuous outcomes or “binomial” for binary outcomes
+(default is “gaussian”).data: The dataset containing all variables required for
+the model.wmean: A vector of treatment assignment weights.
+Default is a vector of ones, implying equal weighting.nboot: The number of bootstrap iterations to perform
+for estimating the uncertainty around the causal estimates.optim_method: The optimization method used to find the
+best parameters in the model (default is ‘BFGS’).seed: A seed value to ensure reproducibility of
+results.parallel: A logical flag indicating whether to perform
+computations in parallel (default is TRUE).ncore: The number of cores to use for parallel
+computation (default is 4).
+
+
model <- bayesmsm(ymodel = y ~ a_1+a_2,
+ nvisit = 2,
+ reference = c(rep(0,2)),
+ comparator = c(rep(1,2)),
+ family = "gaussian",
+ data = testdata,
+ wmean = rep(1, 1000),
+ nboot = 1000,
+ optim_method = "BFGS",
+ parallel = TRUE,
+ seed = 890123,
+ ncore = 4)
+str(model)
+#> List of 6
+#> $ RD_mean : num -3.16
+#> $ RD_sd : num 0.0947
+#> $ RD_quantile: Named num [1:2] -3.34 -2.98
+#> ..- attr(*, "names")= chr [1:2] "2.5%" "97.5%"
+#> $ bootdata :'data.frame': 1000 obs. of 3 variables:
+#> ..$ effect_reference : num [1:1000] 2.36 2.29 2.24 2.33 2.3 ...
+#> ..$ effect_comparator: num [1:1000] -0.951 -0.825 -0.77 -0.811 -0.75 ...
+#> ..$ RD : num [1:1000] -3.31 -3.11 -3.01 -3.14 -3.05 ...
+#> $ reference : num [1:2] 0 0
+#> $ comparator : num [1:2] 1 1
+
+- It returns a model object which contains:
+
+mean, sd, quantile: the mean,
+standard deviation and 95% credible interval of the estimated causal
+effect (ATE). From the above results, the mean of ATE is approximately
+-3.161, which indicates that the expected outcome for always treated
+patients is, on average, 3.161 units less than that for never treated
+patients.bootdata: a data frame containing the bootstrap samples
+for the reference effect, comparator effect, and average treatment
+effect (ATE).reference, comparator: the reference level
+and comparator level the user chooses to compare. Here the reference
+level is never treated (0,0), and the comparator level is always treated
+(1,1).
+
+
5. Visualization functions: plot_ATE,
+plot_APO, plot_est_box
+
The bayesmsm package also provides several other
+functions for visualizing the above results: plot_ATE,
+plot_APO, and plot_est_box. These functions
+help the user better interpret the estimated causal effects.
+
+- Plotting the Average Treatment Effect (ATE)
+
+- The
plot_ATE function generates a plot of the estimated
+ATE with its 95% credible interval.
+
+
+
+
+
+- Plotting the Average Potential Outcome (APO)
+
+- Similarly, the
plot_APO function plots the estimated
+APO for both the reference and comparator level effects.
+
+
+
plot_APO(model, "effect_reference")
+
+
plot_APO(model, "effect_comparator")
+
+
+- Plotting the Distribution of Estimated Treatment Effects
+
+- The
plot_est_box function generates an error bar plot
+of the estimated treatment effects (APO and ATE) from the bootstrap
+samples.
+
+
+
+
+
Reference
+
+- Liu, K. (2021). Bayesian causal inference with longitudinal data.
+Tspace.library.utoronto.ca. https://tspace.library.utoronto.ca/handle/1807/109330
+- Saarela, O., Stephens, D. A., Moodie, E. E. M., & Klein, M. B.
+(2015). On Bayesian estimation of marginal structural models.
+Biometrics, 71(2), 279–288. https://doi.org/10.1111/biom.12269
+- Robins, J. M., Hernán, M. A., & Brumback, B. (2000). Marginal
+structural models and causal inference in epidemiology. Epidemiology,
+11(5), 550–560. https://doi.org/10.1097/00001648-200009000-00011
+- Liu, K., Saarela, O., Feldman, B. M., & Pullenayegum, E. (2020).
+Estimation of causal effects with repeatedly measured outcomes in a
+Bayesian framework. Statistical Methods in Medical Research, 29(9),
+2507–2519. https://doi.org/10.1177/0962280219900362
+
-
-
+
-
-
+
+
-
+
diff --git a/docs/pkgdown.js b/docs/pkgdown.js
deleted file mode 100644
index 9bd6621..0000000
--- a/docs/pkgdown.js
+++ /dev/null
@@ -1,156 +0,0 @@
-/* http://gregfranko.com/blog/jquery-best-practices/ */
-(function($) {
- $(function() {
-
- $('nav.navbar').headroom();
-
- Toc.init({
- $nav: $("#toc"),
- $scope: $("main h2, main h3, main h4, main h5, main h6")
- });
-
- if ($('#toc').length) {
- $('body').scrollspy({
- target: '#toc',
- offset: $("nav.navbar").outerHeight() + 1
- });
- }
-
- // Activate popovers
- $('[data-bs-toggle="popover"]').popover({
- container: 'body',
- html: true,
- trigger: 'focus',
- placement: "top",
- sanitize: false,
- });
-
- $('[data-bs-toggle="tooltip"]').tooltip();
-
- /* Clipboard --------------------------*/
-
- function changeTooltipMessage(element, msg) {
- var tooltipOriginalTitle=element.getAttribute('data-bs-original-title');
- element.setAttribute('data-bs-original-title', msg);
- $(element).tooltip('show');
- element.setAttribute('data-bs-original-title', tooltipOriginalTitle);
- }
-
- if(ClipboardJS.isSupported()) {
- $(document).ready(function() {
- var copyButton = "
";
-
- $("div.sourceCode").addClass("hasCopyButton");
-
- // Insert copy buttons:
- $(copyButton).prependTo(".hasCopyButton");
-
- // Initialize tooltips:
- $('.btn-copy-ex').tooltip({container: 'body'});
-
- // Initialize clipboard:
- var clipboard = new ClipboardJS('[data-clipboard-copy]', {
- text: function(trigger) {
- return trigger.parentNode.textContent.replace(/\n#>[^\n]*/g, "");
- }
- });
-
- clipboard.on('success', function(e) {
- changeTooltipMessage(e.trigger, 'Copied!');
- e.clearSelection();
- });
-
- clipboard.on('error', function(e) {
- changeTooltipMessage(e.trigger,'Press Ctrl+C or Command+C to copy');
- });
-
- });
- }
-
- /* Search marking --------------------------*/
- var url = new URL(window.location.href);
- var toMark = url.searchParams.get("q");
- var mark = new Mark("main#main");
- if (toMark) {
- mark.mark(toMark, {
- accuracy: {
- value: "complementary",
- limiters: [",", ".", ":", "/"],
- }
- });
- }
-
- /* Search --------------------------*/
- /* Adapted from https://github.com/rstudio/bookdown/blob/2d692ba4b61f1e466c92e78fd712b0ab08c11d31/inst/resources/bs4_book/bs4_book.js#L25 */
- // Initialise search index on focus
- var fuse;
- $("#search-input").focus(async function(e) {
- if (fuse) {
- return;
- }
-
- $(e.target).addClass("loading");
- var response = await fetch($("#search-input").data("search-index"));
- var data = await response.json();
-
- var options = {
- keys: ["what", "text", "code"],
- ignoreLocation: true,
- threshold: 0.1,
- includeMatches: true,
- includeScore: true,
- };
- fuse = new Fuse(data, options);
-
- $(e.target).removeClass("loading");
- });
-
- // Use algolia autocomplete
- var options = {
- autoselect: true,
- debug: true,
- hint: false,
- minLength: 2,
- };
- var q;
-async function searchFuse(query, callback) {
- await fuse;
-
- var items;
- if (!fuse) {
- items = [];
- } else {
- q = query;
- var results = fuse.search(query, { limit: 20 });
- items = results
- .filter((x) => x.score <= 0.75)
- .map((x) => x.item);
- if (items.length === 0) {
- items = [{dir:"Sorry 😿",previous_headings:"",title:"No results found.",what:"No results found.",path:window.location.href}];
- }
- }
- callback(items);
-}
- $("#search-input").autocomplete(options, [
- {
- name: "content",
- source: searchFuse,
- templates: {
- suggestion: (s) => {
- if (s.title == s.what) {
- return `${s.dir} >
${s.title}
`;
- } else if (s.previous_headings == "") {
- return `${s.dir} >
${s.title}
> ${s.what}`;
- } else {
- return `${s.dir} >
${s.title}
> ${s.previous_headings} > ${s.what}`;
- }
- },
- },
- },
- ]).on('autocomplete:selected', function(event, s) {
- window.location.href = s.path + "?q=" + q + "#" + s.id;
- });
- });
-})(window.jQuery || window.$)
-
-
diff --git a/docs/pkgdown.yml b/docs/pkgdown.yml
deleted file mode 100644
index 7e0cd95..0000000
--- a/docs/pkgdown.yml
+++ /dev/null
@@ -1,15 +0,0 @@
-pandoc: 3.1.1
-pkgdown: 2.0.7
-pkgdown_sha: ~
-<<<<<<< HEAD
-articles:
- vignette: vignette.html
-last_built: 2024-06-27T10:06Z
-urls:
- reference: https://kuan-liu-lab.github.io/bayesmsm/reference
- article: https://kuan-liu-lab.github.io/bayesmsm/articles
-=======
-articles: {}
-last_built: 2024-01-25T17:31Z
->>>>>>> parent of b799ead (Updated functions bayesmsm, plot_ATE, plot_APO, and vignette)
-
diff --git a/docs/reference/index.html b/docs/reference/index.html
deleted file mode 100644
index a8c11a1..0000000
--- a/docs/reference/index.html
+++ /dev/null
@@ -1,119 +0,0 @@
-
-
Function reference • bayesmsm
-
Skip to contents
-
-
-<<<<<<< HEAD
-