Winter Term Practicum Presentation - Xiao Yan.qmd

---
title: "Bayesian Marginal Structural Models"
subtitle: "An R Package Project"
author: "Xiao Yan | Supervisor: Kuan Liu"
format:
  revealjs: 
    theme: serif
    slide-number: true
    chalkboard: 
      buttons: false
    preview-links: auto
    css: styles.css
    footer: "Practicum Term 2 Presentation"

---

## Outline {transition="fade"}

-   Background
-   Motivation
-   Methods
-   Results
-   Summary


## Background {transition="slide"}

::: {.fragment .fade-in-then-semi-out}
-   Observational studies offer a viable, efficient, and low-cost design to readily gather evidence on exposure effects.
:::

::: {.fragment .fade-in-then-semi-out}
-   Although more practical, exposure mechanism is nonrandomized and <mark>causal inference methods</mark> are required to draw causal conclusions.
:::

::: {.fragment .fade-in-then-semi-out}
-   Popular approaches used in health research are predominantly <mark>frequentist</mark> methods.
:::

## Motivation {transition="slide"}

::: {.fragment}
-   A stream of Bayesian causal inference methods has been developed.

::: {.fragment}
-   Bayesian approaches have unique estimation features that are useful in many settings, however, there is a general lack of open-access software packages to carry out these analyses.
:::

::: {.fragment}
::: {.fragment .highlight-red}
-   Goal: build a user-friendly R package for Bayesian Marginal Structural Models (BMSMs).
:::
:::
:::

## Methodology {transition="slide"}

2 estimation steps:

::: {.fragment .fade-in-then-semi-out}
-   Step 1. Bayesian treatment effect weight (similar to PS)
:::

::: {.fragment .fade-in-then-semi-out}
-   Step 2. Bayesian non-parametric bootstrap to maximize the utility function with respect to the causal effect
:::

<!-- ## Methodology {transition="slide"} -->

<!-- Next, we introduce the methodology used in the `bayesmsm` package to perform step 2 of Bayesian Marginal Structural Models (BMSMs) analysis. -->

<!-- ::: {.fragment .highlight-red} -->
<!-- -   Main function: 'bayesmsm' -->
<!-- ::: -->
<!-- -   Used to conduct Bayesian non-parametric bootstrap to calculate causal effect. -->

## Notation {transition="fade"}

::: {.fragment .fade-in-then-semi-out}
-   A longitudinal observational study with n subjects indexed by $i$, $i = 1, \ldots, n$ and $J$ number of visits indexed by $j$, $j = 1, \ldots, J$.
:::

::: {.fragment .fade-in-then-semi-out}
-   $Y_{i}$, $X_{ij}$ and $Z_{ij}$ are random variables representing an end-of-study response, covariates and the treatment for individual $i$ at visit $j$.
:::

::: {.fragment .fade-in-then-semi-out}
-   History up to visit $j$ are denoted as $\bar{X}_{ij}$ and $\bar{Z}_{ij}$.
:::

## Methodology: BMSMs {.smaller transition="fade"}

Using Bayesian decision theory and importance sampling technique, we maximize an expected utility function (a function involving only $\theta$), $\textbf{u}_{\mathcal{E}}(\Theta, \bar{v}_{i}^*)$, via posterior predictive inference,

$$ \hat{\Theta} 
 = argmax_{\theta} \int_{\bar{v}_{i}^*}  u_{\mathcal{E}}(\Theta, \bar{v}_{i}^*)P_{\mathcal{E}}(\bar{v}_{i}^* \mid \textbf{V}_n) \ d\bar{v}_{i}^*  \nonumber \\ $$
 
$$ = argmax_{\theta}\int_{\bar{v}_{i}^*}  u_{\mathcal{E}}(\Theta, \bar{v}_{i}^*) \frac{P_{\mathcal{E}}(\bar{v}_{i}^* \mid \textbf{V}_n) }{P_{\mathcal{O}}(\bar{v}_{i}^* \mid \textbf{V}_n)}P_{\mathcal{O}}(\bar{v}_{i}^* \mid \textbf{V}_n) \ d\bar{v}_{i}^* $$

-   $u(\Theta, \bar{v}_{i}^*)= log P_{\mathcal{E}}( Y_{i}^* \mid \bar{z}_{iJ}^*; \Theta)$, <mark>utility function</mark>

-   $w_{i}^* = \frac{P_{\mathcal{E}}(\bar{v}_{i}^* \mid \textbf{v}_n)}{P_{\mathcal{O}}(\bar{v}_{i}^* \mid \textbf{v}_n)}$ expanded to <mark>treatment assignment weight</mark>

```{r, echo=FALSE}
testdata <- readr::read_csv("continuous_outcome_data.csv")

library(WeightIt)

Wmsm <- weightitMSM(
  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, 
  method = "ps",
  stabilize = TRUE)

# Wmsm$weights
```

## Weighted log-likelihood {transition="fade"}

::: panel-tabset

### Normal Y
\begin{equation}
\mathcal{l}(\theta, \sigma^2 | Y, A) = \sum_{i=1}^{n} w_i \left( -\frac{1}{2} \log(\sigma^2) - \frac{1}{2\sigma^2} (y_i - A_i \theta)^2 \right)
\end{equation}

-   $\theta$: causal parameters on the mean
-   $\sigma$: causal parameter on the sd
-   $A$: design matrix of the causal outcome model

### Binary Y

\begin{equation}
\mathcal{l}(\beta | Y, A) = \sum_{i=1}^{n} w_i \left( Y_i \eta_i - \log(1 + \exp(\eta_i)) \right)
\end{equation}

-   $\beta$: causal parameters on the log-odds scale
-   $\eta$: $A$ * $\beta$, the linear predictor
:::

::: footer
Reference: [linearML](https://www.stat.cmu.edu/~cshalizi/mreg/15/lectures/06/lecture-06.pdf) and [optim](https://www.ime.unicamp.br/~cnaber/optim_1.pdf)
:::

## Parallel computing

Parallel computing for faster bootstrap calculation.

```{.r code-line-numbers="|2|5-7"}
if (parallel == TRUE){
    numCores <- ncore
    registerDoParallel(cores = numCores)
    
    results <- foreach(i=1:nboot,
                       .combine = 'rbind',
                       .packages = 'MCMCpack') %dopar% {
                       
                       ... # Bootstrap calculation
                       
    }
}
```



## Putting it all together {transition="fade"}

The complete function 'bayesmsm', at a glance:

``` r
bayesmsm <- function(ymodel,
                     nvisit,
                     reference = c(rep(0,nvisit)), # An example of never treated
                     comparator = c(rep(1,nvisit)),
                     family = "gaussian", # "gaussian" or "binomial"
                     data,
                     wmean = rep(1, nrow(data)),
                     nboot = 1000,
                     optim_method = 'BFGS',
                     estimand = 'RD',
                     parallel = TRUE,
                     ncore = 6){

  # load all the required R packages;
  if (!require(foreach)){
    install.packages("foreach",repos="http://cran.r-project.org")
    library(foreach)
  }
  if (!require(doParallel)){
    install.packages("doParallel",repos="http://cran.r-project.org")
    library(doParallel)
  }
  if (!require(MCMCpack)){
    install.packages("MCMCpack",repos="http://cran.r-project.org")
    library(MCMCpack)
  }

  # return error message if the input weight vector has different length comparing to the outcome Y;
  if (length(wmean) != nrow(data)) {
    stop("The length of the weight vector does not match the length of Y.")
  }

  # load utility functions
  extract_variables <- function(formula) {
    # Get the terms of the formula
    formula_terms <- terms(formula)

    # Extract the response variable name (if there is one)
    response_variable <- attr(formula_terms, "response")
    response_name <- if (response_variable > 0) {
      all_vars <- all.vars(formula)
      all_vars[response_variable]
    } else {NA}

    # Extract predictor variable names
    predictor_names <- attr(formula_terms, "term.labels")

    # Return a list of response and predictor variables
    list(response = response_name, predictors = predictor_names)
  }

  variables <- extract_variables(ymodel) # Extract variable names from the formula
  Y_name <- variables$response

  Y <- data[[Y_name]]
  A_base <- data.frame(matrix(data = NA,
                              nrow = nrow(data),
                              ncol = length(variables$predictors)))
  for (i in 1:length(variables$predictors)){
    initial_vector <- variables$predictors[i]
    split_vector <- strsplit(initial_vector, ":")
    new_vector <- unlist(split_vector)
    if (length(new_vector)==1){
      A_base[,i] <-  data[, new_vector]
    } else if (length(new_vector)>1){
      A_base[,i] <-  apply(data[, new_vector],1,prod)
    }
  }

  A <- cbind(1, A_base)
  colnames(A)[2:ncol(A)]<- variables$predictors

  wloglik_normal<-function(param,
                           Y,
                           A,
                           weight){
    #number of observations;
    n <- length(Y)
    theta <- param[1:dim(A)[2]] #causal parameters on the mean
    #number of parameter is determined by number of treatment variables, plus intercept;
    sigma <- param[(dim(A)[2]+1)] # the remaining the parameter represent the standard deviation;
    mmat <- as.matrix(A) #design matrix of the causal outcome model, e.g., A = cbind(1, a_1, a_2);
    logl<- -0.5*log(sigma**2) - 0.5*((Y - mmat%*%theta)**2)/(sigma**2)
    wlogl<-sum(weight*logl)
    return(wlogl)
  }

  wloglik_binomial <- function(param,
                               Y,
                               A,
                               weight){
    # number of observations;
    n <- length(Y)
    beta <- param[1:dim(A)[2]] # causal parameters on the log-odds scale (no sigma for binomial?)
    mmat <- as.matrix(A)
    eta<-mmat %*% beta # linear predictor
    logl <- Y*eta - log(1+exp(eta))
    wlogl<-sum(weight*logl)
    return(wlogl)
  }

  expit <- function(x){exp(x) / (1+exp(x))}

  if (family == "gaussian"){
    wfn = wloglik_normal
    inits1 <- c(rep(0.1, length(A)), 4)  # Default initial values, 4 is for the SD;
  } else if (family == "binomial"){
    wfn = wloglik_binomial
    inits1 <- c(rep(0.1, length(A)))
  } else if (!family %in% c("gaussian","binomial")){
    stop("Current version only handles continuous (gaussian) and binary (binomial) outcomes.")
  }


  # Parallel computing for bootstrapping
  if (parallel == TRUE){
    numCores <- ncore
    registerDoParallel(cores = numCores)

    results <- foreach(i=1:nboot,
                       .combine = 'rbind',
                       .packages = 'MCMCpack') %dopar% {

      calculate_effect <- function(intervention_levels, variables, param_estimates) {
        # Start with the intercept term
        effect<-effect_intercept<-param_estimates[1]

        # Go through each predictor and add its contribution
        for (i in 1:length(variables$predictors)) {
          term <- variables$predictors[i]
          term_variables <- unlist(strsplit(term, ":"))
          term_index <- which(names(param_estimates) == term)

          # Calculate the product of intervention levels for the interaction term
          term_contribution <- param_estimates[term_index]
          for (term_variable in term_variables) {
            var_index <- which(variables$predictors == term_variable)
            term_contribution <- term_contribution * intervention_levels[var_index]
          }

          # Add the term contribution to the effect
          effect <- effect + term_contribution
        }

        return(effect)
      }

      results.it <- matrix(NA, 1, 3) #result matrix, three columns for bootest, effect_ref, and effect_comp;

      alpha <- as.numeric(rdirichlet(1, rep(1.0, length(Y))))

      maxim <- optim(inits1,
                     fn = wfn,
                     Y = Y,
                     A = A,
                     weight = alpha * wmean,
                     control = list(fnscale = -1),
                     method = optim_method,
                     hessian = FALSE)

      names(maxim$par) <- c("(Intercept)", variables$predictors)

      # Calculate the effects
      results.it[1,1] <- calculate_effect(reference, variables, param_estimates=maxim$par)
      results.it[1,2] <- calculate_effect(comparator, variables, param_estimates=maxim$par)

      # Calculate the ATE
      if (family == "binomial") { # binary outcomes
        if (estimand == "RD") { # Risk Difference
          results.it[1,3] <- expit(results.it[1,2]) - expit(results.it[1,1])
        } else if (estimand == "RR") { # Relative Risk
          results.it[1,3] <- expit(results.it[1,2]) / expit(results.it[1,1])
        } else if (estimand == "OR") { # Odds Ratio
          results.it[1,3] <- (expit(results.it[1,2]) / (1 - expit(results.it[1,2]))) /
            (expit(results.it[1,1]) / (1 - expit(results.it[1,1])))
        }
      } else if (family == "gaussian"){ # continuous outcomes
        if (estimand == "RD") { # Risk Difference
          results.it[1,3] <- results.it[1,2] - results.it[1,1]
        }
      }

      # combining parallel results;
      cbind(i,results.it) #end of parallel;
    }

    #saving output for the parallel setting;
    return(list(
      mean = mean(results[,4]),
      sd = sqrt(var(results[,4])),
      quantile = quantile(results[,4], probs = c(0.025, 0.975)),
      bootdata = data.frame(effect_reference = results[,2],
                             effect_comparator = results[,3],
                             ATE = results[,4]),
      reference = reference,
      comparator = comparator
    ))

  }

  else if (parallel == FALSE) {

    bootest <- numeric(nboot)
    effect_reference <- numeric(nboot)
    effect_comparator <- numeric(nboot)

    for (j in 1:nboot) {
      alpha <- as.numeric(rdirichlet(1, rep(1.0, length(Y))))

      maxim <- optim(inits1,
                     fn = wfn,
                     Y = Y,
                     A = A,
                     weight = alpha * wmean,
                     control = list(fnscale = -1),
                     method = optim_method,
                     hessian = FALSE)

      names(maxim$par) <- c("(Intercept)", variables$predictors)

      # Calculate the effects
      effect_reference[j] <- calculate_effect(reference, variables, param_estimates=maxim$par)
      effect_comparator[j] <- calculate_effect(comparator, variables, param_estimates=maxim$par)

      # Calculate the ATE
      if (family == "binomial") { # binary outcomes
        if (estimand == "RD") { # Risk Difference
          bootest[j] <- expit(effect_comparator[j]) - expit(effect_reference[j])
        } else if (estimand == "RR") { # Relative Risk
          bootest[j] <- expit(effect_comparator[j]) / expit(effect_reference[j])
        } else if (estimand == "OR") { # Odds Ratio
          bootest[j] <- (expit(effect_comparator[j]) / (1 - expit(effect_comparator[j]))) /
            (expit(effect_reference[j]) / (1 - expit(effect_reference[j])))
        }
      } else if (family == "gaussian"){ # continuous outcomes
        if (estimand == "RD") { # Risk Difference
          bootest[j] <- effect_comparator[j] - effect_reference[j]
        } else if (estimand %in% c("RR","OR")) {
          # print a warning message that say for continuous outcome, RR and OR specification are ignored. RD is the causal estimate;
          warning("For continuous outcomes, RR and OR specifications are ignored. RD is the only applicable causal estimate.")
        }
      }

    }

    #saving output for the non-parallel setting;
    return(list(
      mean = mean(bootest),
      sd = sqrt(var(bootest)),
      quantile = quantile(bootest, probs = c(0.025, 0.975)),
      bootdata = data.frame(effect_reference, effect_comparator, ATE=bootest),
      reference = reference,
      comparator = comparator
    ))

  }
}
```

```{r}
bayesmsm <- function(ymodel,
                     nvisit,
                     reference = c(rep(0,nvisit)), # An example of never treated
                     comparator = c(rep(1,nvisit)),
                     family = "gaussian", # "gaussian" or "binomial"
                     data,
                     wmean = rep(1, 1000),
                     nboot = 1000,
                     optim_method = 'BFGS',
                     estimand = 'RD',
                     parallel = TRUE,
                     ncore = 6){

  # load all the required R packages;
  if (!require(foreach)){
    install.packages("foreach",repos="http://cran.r-project.org")
    library(foreach)
  }
  if (!require(doParallel)){
    install.packages("doParallel",repos="http://cran.r-project.org")
    library(doParallel)
  }
  if (!require(MCMCpack)){
    install.packages("MCMCpack",repos="http://cran.r-project.org")
    library(MCMCpack)
  }

  # return error message if the input weight vector has different length comparing to the outcome Y;
  if (length(wmean) != nrow(data)) {
    stop("The length of the weight vector does not match the length of Y.")
  }

  # load utility functions
  extract_variables <- function(formula) {
    # Get the terms of the formula
    formula_terms <- terms(formula)

    # Extract the response variable name (if there is one)
    response_variable <- attr(formula_terms, "response")
    response_name <- if (response_variable > 0) {
      all_vars <- all.vars(formula)
      all_vars[response_variable]
    } else {NA}

    # Extract predictor variable names
    predictor_names <- attr(formula_terms, "term.labels")

    # Return a list of response and predictor variables
    list(response = response_name, predictors = predictor_names)
  }

  variables <- extract_variables(ymodel) # Extract variable names from the formula
  Y_name <- variables$response

  Y <- data[[Y_name]]
  A_base <- data.frame(matrix(data = NA,
                              nrow = nrow(data),
                              ncol = length(variables$predictors)))
  for (i in 1:length(variables$predictors)){
    initial_vector <- variables$predictors[i]
    split_vector <- strsplit(initial_vector, ":")
    new_vector <- unlist(split_vector)
    if (length(new_vector)==1){
      A_base[,i] <-  data[, new_vector]
    } else if (length(new_vector)>1){
      A_base[,i] <-  apply(data[, new_vector],1,prod)
    }
  }

  A <- cbind(1, A_base)
  colnames(A)[2:ncol(A)]<- variables$predictors

  wloglik_normal<-function(param,
                           Y,
                           A,
                           weight){
    #number of observations;
    n <- length(Y)
    theta <- param[1:dim(A)[2]] #causal parameters on the mean
    #number of parameter is determined by number of treatment variables, plus intercept;
    sigma <- param[(dim(A)[2]+1)] # the remaining the parameter represent the standard deviation;
    mmat <- as.matrix(A) #design matrix of the causal outcome model, e.g., A = cbind(1, a_1, a_2);
    logl<- -0.5*log(sigma**2) - 0.5*((Y - mmat%*%theta)**2)/(sigma**2)
    wlogl<-sum(weight*logl)
    return(wlogl)
  }

  wloglik_binomial <- function(param,
                               Y,
                               A,
                               weight){
    # number of observations;
    n <- length(Y)
    beta <- param[1:dim(A)[2]] # causal parameters on the log-odds scale (no sigma for binomial?)
    mmat <- as.matrix(A)
    eta<-mmat %*% beta # linear predictor
    logl <- Y*eta - log(1+exp(eta))
    wlogl<-sum(weight*logl)
    return(wlogl)
  }

  expit <- function(x){exp(x) / (1+exp(x))}

  if (family == "gaussian"){
    wfn = wloglik_normal
    inits1 <- c(rep(0.1, length(A)), 4)  # Default initial values, 4 is for the SD;
  } else if (family == "binomial"){
    wfn = wloglik_binomial
    inits1 <- c(rep(0.1, length(A)))
  } else if (!family %in% c("gaussian","binomial")){
    stop("Current version only handles continuous (gaussian) and binary (binomial) outcomes.")
  }


  # Parallel computing for bootstrapping
  if (parallel == TRUE){
    numCores <- ncore
    registerDoParallel(cores = numCores)

    results <- foreach(i=1:nboot,
                       .combine = 'rbind',
                       .packages = 'MCMCpack') %dopar% {

      calculate_effect <- function(intervention_levels, variables, param_estimates) {
        # Start with the intercept term
        effect<-effect_intercept<-param_estimates[1]

        # Go through each predictor and add its contribution
        for (i in 1:length(variables$predictors)) {
          term <- variables$predictors[i]
          term_variables <- unlist(strsplit(term, ":"))
          term_index <- which(names(param_estimates) == term)

          # Calculate the product of intervention levels for the interaction term
          term_contribution <- param_estimates[term_index]
          for (term_variable in term_variables) {
            var_index <- which(variables$predictors == term_variable)
            term_contribution <- term_contribution * intervention_levels[var_index]
          }

          # Add the term contribution to the effect
          effect <- effect + term_contribution
        }

        return(effect)
      }

      results.it <- matrix(NA, 1, 3) #result matrix, three columns for bootest, effect_ref, and effect_comp;

      alpha <- as.numeric(rdirichlet(1, rep(1.0, length(Y))))

      maxim <- optim(inits1,
                     fn = wfn,
                     Y = Y,
                     A = A,
                     weight = alpha * wmean,
                     control = list(fnscale = -1),
                     method = optim_method,
                     hessian = FALSE)

      names(maxim$par) <- c("(Intercept)", variables$predictors)

      # Calculate the effects
      results.it[1,1] <- calculate_effect(reference, variables, param_estimates=maxim$par)
      results.it[1,2] <- calculate_effect(comparator, variables, param_estimates=maxim$par)

      # Calculate the ATE
      if (family == "binomial") { # binary outcomes
        if (estimand == "RD") { # Risk Difference
          results.it[1,3] <- expit(results.it[1,2]) - expit(results.it[1,1])
        } else if (estimand == "RR") { # Relative Risk
          results.it[1,3] <- expit(results.it[1,2]) / expit(results.it[1,1])
        } else if (estimand == "OR") { # Odds Ratio
          results.it[1,3] <- (expit(results.it[1,2]) / (1 - expit(results.it[1,2]))) /
            (expit(results.it[1,1]) / (1 - expit(results.it[1,1])))
        }
      } else if (family == "gaussian"){ # continuous outcomes
        if (estimand == "RD") { # Risk Difference
          results.it[1,3] <- results.it[1,2] - results.it[1,1]
        }
      }

      # combining parallel results;
      cbind(i,results.it) #end of parallel;
    }

    #saving output for the parallel setting;
    return(list(
      mean = mean(results[,4]),
      sd = sqrt(var(results[,4])),
      quantile = quantile(results[,4], probs = c(0.025, 0.975)),
      bootdata = data.frame(effect_reference = results[,2],
                             effect_comparator = results[,3],
                             ATE = results[,4]),
      reference = reference,
      comparator = comparator
    ))

  }

  else if (parallel == FALSE) {

    bootest <- numeric(nboot)
    effect_reference <- numeric(nboot)
    effect_comparator <- numeric(nboot)

    for (j in 1:nboot) {
      alpha <- as.numeric(rdirichlet(1, rep(1.0, length(Y))))

      maxim <- optim(inits1,
                     fn = wfn,
                     Y = Y,
                     A = A,
                     weight = alpha * wmean,
                     control = list(fnscale = -1),
                     method = optim_method,
                     hessian = FALSE)

      names(maxim$par) <- c("(Intercept)", variables$predictors)

      # Calculate the effects
      effect_reference[j] <- calculate_effect(reference, variables, param_estimates=maxim$par)
      effect_comparator[j] <- calculate_effect(comparator, variables, param_estimates=maxim$par)

      # Calculate the ATE
      if (family == "binomial") { # binary outcomes
        if (estimand == "RD") { # Risk Difference
          bootest[j] <- expit(effect_comparator[j]) - expit(effect_reference[j])
        } else if (estimand == "RR") { # Relative Risk
          bootest[j] <- expit(effect_comparator[j]) / expit(effect_reference[j])
        } else if (estimand == "OR") { # Odds Ratio
          bootest[j] <- (expit(effect_comparator[j]) / (1 - expit(effect_comparator[j]))) /
            (expit(effect_reference[j]) / (1 - expit(effect_reference[j])))
        }
      } else if (family == "gaussian"){ # continuous outcomes
        if (estimand == "RD") { # Risk Difference
          bootest[j] <- effect_comparator[j] - effect_reference[j]
        } else if (estimand %in% c("RR","OR")) {
          # print a warning message that say for continuous outcome, RR and OR specification are ignored. RD is the causal estimate;
          warning("For continuous outcomes, RR and OR specifications are ignored. RD is the only applicable causal estimate.")
        }
      }

    }

    #saving output for the non-parallel setting;
    return(list(
      mean = mean(bootest),
      sd = sqrt(var(bootest)),
      quantile = quantile(bootest, probs = c(0.025, 0.975)),
      bootdata = data.frame(effect_reference, effect_comparator, ATE=bootest),
      reference = reference,
      comparator = comparator
    ))

  }
}
```

## Results (DAG) {transition="slide"}

Figure: Longitudinal Directed Acyclic Graph (DAG) for 2 visits

```{mermaid}
flowchart TD
    Cov1[Covariates 1] --> Cov2[Covariates 2]
    Cov1 --> Treat1[Treatment 1]
    Cov1 --> Treat2
    Cov1 --> Outcome[Outcome]
    Cov2 --> Treat2
    Cov2 --> Outcome
    Treat1 --> Treat2[Treatment 2]
    Treat1 --> Cov2
    Treat1 --> Outcome
    Treat2 --> Outcome
```

## Results {transition="fade"}

Example usage of function `bayesmsm`:

```{.r code-line-numbers="|1-2|3-4|5|7|10"}
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 = Wmsm$weights,
                  nboot = 1000,
                  optim_method = "BFGS",
                  estimand = "RD",
                  parallel = TRUE,
                  ncore = 6)
```

```{r}
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 = Wmsm$weights,
                           nboot = 1000,
                           optim_method = "BFGS",
                           estimand = "RD",
                           parallel = TRUE,
                           ncore = 6)
```

::: footer
Inspired by `gfoRmula` and `ltmle`, see: [Causal analysis with time-varying treatment](https://kuan-liu.github.io/causal_Quarto/section3.html#implementing-targeted-maximum-likelihood-estimation)
:::

## Bootstrap Results {transition="fade"}

``` r
model$bootdata
```

::: columns
::: {.column width="100%"}
```{r}
head(model$bootdata)
```
:::
:::

::: {.fragment}
This model output allows users to plot and summarize the bootstrap results.
:::

## Results {transition="fade"}

There are also other functions in this package available to visualize and interpret the results:

-   'plot_ATE'
-   'plot_APO'
-   'plot_est_box'

## Other functions in the package {.smaller .scrollable transition="fade"}

::: panel-tabset
### plot_ATE

```{r}
plot_ATE <- function(input,
                     col_density = "blue",
                     fill_density = "lightblue",
                     main = "Posterior Predictive Distribution of Average Treatment Effect (ATE）",
                     xlab = "ATE", ylab = "Posterior Predictive Distribution",
                     xlim = NULL, ylim = NULL, ...) {
  # Check if input is either a data frame or part of a model object
  if (is.list(input) && "bootdata" %in% names(input)) {
    # If input is a list and has bootdata, check for ATE column within bootdata
    if ("ATE" %in% names(input$bootdata)) {
      ate_values <- input$bootdata$ATE
    } else {
      stop("bootdata within the model object must have an 'ATE' column.")
    }
  } else if (is.data.frame(input) && "ATE" %in% names(input)) {
    ate_values <- input$ATE
  } else if (is.vector(input)) {
    ate_values <- input
  } else {
    stop("input must be a vector of ATE estimates, a data frame, or a model object containing a 'bootdata' data frame with an 'ATE' column.")
  }

  ate_density <- density(ate_values)
  ci <- quantile(ate_values, probs = c(0.025, 0.975))
  density_ci <- density(ate_values, from = ci[1], to = ci[2])

  plot(ate_density, col = col_density, main = main, xlab = xlab, ylab = ylab, xlim = xlim, ylim = ylim, ...)
  polygon(c(density_ci$x, rev(density_ci$x)), c(rep(min(ate_density$y), length(density_ci$x)), rev(density_ci$y)), col = rgb(0, 0, 1, alpha = 0.3))
  abline(v = mean(ate_values), col = "purple", lwd = 2, lty = 3)
  abline(v = ci[1], col = "darkgreen", lty = 2)
  abline(v = ci[2], col = "darkgreen", lty = 2)

  legend_text <- c("ATE Density",
                   paste("Mean:", round(mean(ate_values), 3)),
                   paste("95% CI: [", round(ci[1], 3), ",", round(ci[2], 3), "]"))

  legend("topright", legend = legend_text,
         col = c(col_density, "purple", "darkgreen"),
         lwd = 2, lty = c(1, 3, 2))
}

plot_ATE(model)
```

### plot_est_box

```{r}
plot_est_box <- function(input, ...) {
  # Extract bootdata from the model or use the data frame directly
  bootdata <- if (is.data.frame(input)) {
    input
  } else if ("bootdata" %in% names(input)) {
    input$bootdata
  } else {
    stop("Input must be a data frame or a model object containing 'bootdata'.")
  }

  # Validate bootdata
  required_columns <- c("effect_comparator", "effect_reference", "ATE")
  if (!all(required_columns %in% names(bootdata))) {
    stop("bootdata must contain 'effect_comparator', 'effect_reference', and 'ATE' columns.")
  }

  # Adjust margins if necessary
  par(mar = c(5, 4, 4, 3) + 0.1) # Adjust the last value if text is plotted outside; bottom, left, top, and right margins

  # Calculate means and standard deviations
  means <- sapply(bootdata[required_columns], mean)
  # sds <- sapply(bootdata[required_columns], sd)
  # ses <- sds / sqrt(nrow(bootdata))
  lowerbd <- sapply(bootdata[required_columns], function(x) quantile(x, probs = 0.025))
  upperbd <- sapply(bootdata[required_columns], function(x) quantile(x, probs = 0.975))

  # Define the position for each point
  position <- 1:length(means)

  # Define some offsets for text placement
  text_offset <- (max(upperbd) - min(lowerbd)) * 0.05

  # Plotting
  plot(position, means, ylim = range(lowerbd - text_offset, upperbd + text_offset), pch = 19, xaxt = "n", # round down vs round up;
       xlab = "Treatment Level", ylab = "Effect", main = "Treatment Effect Estimates", ...)
  axis(1, at = position, labels = c("Comparator Level", "Reference Level", "ATE"))

  # Error bars
  arrows(position, lowerbd, position, upperbd, angle = 90, code = 3, length = 0.1, ...)

  # Adding text for means and CIs
  # text(position, lowerbd - 0.1, labels = paste("Mean:", round(means, 2), "\n95% CI:", round(lowerbd, 2), "-", round(upperbd, 2)), adj = c(0,1))

  # for (i in seq_along(means)) {
  #   text(position[i], upperbd[i] + text_offset, labels = paste("Mean:", round(means[i], 2), "\n95% CI:", round(lowerbd[i], 2), "-", round(upperbd[i], 2)), cex = 0.8, pos = 3)
  # }

  # Adjust the y-position offset for clarity
  offset <- 0.1
  # Adding text for the first and second items with left adjustment
  text(position[1:2], lowerbd[1:2] - offset, labels = paste("Mean:", round(means[1:2], 2), "\n95% CI: [", round(lowerbd[1:2], 2), ", ", round(upperbd[1:2], 2), "]"), adj = c(0,1), ...)
  # Adding text for the third item with right adjustment
  text(position[3], upperbd[3] + offset, labels = paste("Mean:", round(means[3], 2), "\n95% CI: [", round(lowerbd[3], 2), ", ", round(upperbd[3], 2), "]"), adj = c(1,0), ...)

  # Check if the input is a model and extract treatment sequences if they exist
  has_treatment_info <- "reference" %in% names(input) && "comparator" %in% names(input)

  # Conditional treatment sequence information below x-axis labels
  if (has_treatment_info) {
    # mtext(paste("(", paste(input$reference, collapse = ", "), ")", sep = ""), side = 1, at = position[2], line = 2.5, cex = 0.7)
    # mtext(paste("(", paste(input$comparator, collapse = ", "), ")", sep = ""), side = 1, at = position[1], line = 2.5, cex = 0.7)
    mtext(paste("(", paste(input$reference, collapse = ", "), ")", sep = ""), side = 1, at = position[2], line = 2)
    mtext(paste("(", paste(input$comparator, collapse = ", "), ")", sep = ""), side = 1, at = position[1], line = 2)
  }
}

plot_est_box(model)
```
:::

## Summary {transition="slide"}

::: {.fragment .fade-in-then-semi-out}
-   Suitable for both continuous and binary Y.
:::

::: {.fragment .fade-in-then-semi-out}
-   Parallel computing option is provided in `bayesmsm` for faster calculation.
:::

::: {.fragment .fade-in-then-semi-out}
-   Simplifies complex BMSM analysis for users.
:::

::: {.fragment .fade-in-then-semi-out}
-   Our package can be downloaded on GitHub: [Kuan-Liu-Lab/bayesmsm](https://kuan-liu-lab.github.io/bayesmsm/)
:::

## Future work {transition="slide"}

::: {.fragment .fade-in-then-semi-out}
-   Bayesian parametric estimation of treatment assignment weights (step 1 of BMSM).
:::

::: {.fragment .fade-in-then-semi-out}
-   Improve computational efficiency for larger datasets.
:::

::: {.fragment .fade-in-then-semi-out}
-   Write documentation (vignette) for this package.
:::

::: {.fragment .fade-in-then-semi-out}
-   Possibly extend this package for survival outcomes.
:::

## References {transition="fade"}

-   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

## References
-   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

## Questions? {transition="fade"}

Thank you for your attention ;)

::: footer
Xiao Yan (Supervisor: Kuan Liu)
:::