---
title: "5. Incorporating propensity scores analysis in psborrow2"
author: "Isaac Gravestock and Matt Secrest"
output:rmarkdown::html_vignette:
vignette: >
  %\VignetteIndexEntry{5. Incorporating propensity scores analysis in psborrow2}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
bibliography: references.bib
link-citations: yes
---




``` r
library(psborrow2)
library(cmdstanr)
```

Propensity scores (PS) methods offer various ways to adjust analyses for differences in groups of patients.
@austin2013 discusses various approaches for using PS with survival analyses to obtain effect measures similar to
randomized controlled trials. @wang2021 discuss using PS for IPTW, matching and stratification in combination with a Bayesian analysis.
These methods allow for the separation of the design and analysis into two stages, which may be attractive in a regulatory setting.
Another approach is the direct inclusion of the PS as a covariate in the outcome model.

## Alternative PS Weights with `WeightIt`

The [`WeightIt`](https://cran.r-project.org/package=WeightIt) package can calculate PS and other balancing
weights with a number of different methods, such as generalized boosted modeling (`method = "gbm"`).
In addition, weights can be calculated differently for different estimands. Here, we specifying
`estimand = "ATT"`, to calculate weights for estimating the average treatment effect among the treated (ATT).


``` r
library(WeightIt)
#> Error in `library()`:
#> ! there is no package called 'WeightIt'

example_dataframe <- as.data.frame(example_matrix)
example_dataframe$int <- 1 - example_dataframe$ext

weightit_model <- weightit(
  int ~ cov1 + cov2 + cov3 + cov4,
  data = example_dataframe,
  method = "gbm",
  estimand = "ATT"
)
#> Error in `weightit()`:
#> ! could not find function "weightit"
summary(weightit_model)
#> Error:
#> ! object 'weightit_model' not found
```

Another useful package is [cobalt](https://cran.r-project.org/package=cobalt), which provides tools for
assessing balance between groups after weighting or matching. It is compatible with many matching and weighting packages.
See the [vignette](https://CRAN.R-project.org/package=cobalt/vignettes/cobalt.html) for more details.
We can use the `cobalt` package to assess balance with `bal.plot()`.


``` r
library(cobalt)
#> Error in `library()`:
#> ! there is no package called 'cobalt'
bal.plot(weightit_model)
#> Error in `bal.plot()`:
#> ! could not find function "bal.plot"
```

"Love plots" can also be generated that compare the populations before and after weighting:


``` r
love.plot(weightit_model, stars = "std")
#> Error in `love.plot()`:
#> ! could not find function "love.plot"
```

The PS values can be copied into the data set and the analysis object can be constructed as before.




``` r
example_dataframe$att <- weightit_model$weights
#> Error:
#> ! object 'weightit_model' not found

example_matrix_att <- create_data_matrix(
  example_dataframe,
  ext_flag_col = "ext",
  outcome = c("time", "cnsr"),
  trt_flag_col = "trt",
  weight_var = "att"
)
#> Error in `create_data_matrix()`:
#> ! Assertion on 'weight_var' failed: Must be a subset of {'id','ext','trt','cov4','cov3','cov2','cov1','time','status','cnsr','resp','int'}, but has additional elements {'att'}.

analysis_att <- create_analysis_obj(
  data_matrix = example_matrix_att,
  outcome = outcome_surv_exponential("time", "cnsr", prior_normal(0, 10000), weight_var = "att"),
  borrowing = borrowing_full("ext"),
  treatment = treatment_details("trt", prior_normal(0, 10000)),
  quiet = TRUE
)
#> Error:
#> ! object 'example_matrix_att' not found

result_att <- mcmc_sample(analysis_att, seed = 123)
#> Error in `h()`:
#> ! error in evaluating the argument 'x' in selecting a method for function 'mcmc_sample': object 'analysis_att' not found
```



## Matching with `MatchIt`

A variety of matching methods, including PS matching are implemented in the
[MatchIt](https://cran.r-project.org/package=MatchIt) package.

As described in the [Getting Started vignette](https://CRAN.R-project.org/package=MatchIt/vignettes/MatchIt.html),
it can be useful to check the imbalance before matching.


``` r
library(MatchIt)
#> Error in `library()`:
#> ! there is no package called 'MatchIt'
# No matching; constructing a pre-match matchit object
no_match <- matchit(trt ~ cov1 + cov2 + cov3 + cov4,
  data = example_dataframe,
  method = NULL, distance = "glm"
)
#> Error in `matchit()`:
#> ! could not find function "matchit"
summary(no_match)
#> Error:
#> ! object 'no_match' not found
```

Here we are matching treated to untreated to select the most comparable control group, regardless of whether they are
internal or external. For simplicity let's try a 1:1 nearest matching approach.


``` r
match_11 <- matchit(trt ~ cov1 + cov2 + cov3 + cov4,
  data = example_dataframe,
  method = "nearest", distance = "glm"
)
#> Error in `matchit()`:
#> ! could not find function "matchit"
summary(match_11)
#> Error:
#> ! object 'match_11' not found
```


``` r
set.seed(123)
plot(match_11, type = "jitter", interactive = FALSE)
#> Error in `h()`:
#> ! error in evaluating the argument 'x' in selecting a method for function 'plot': object 'match_11' not found
```

Determining whether the balance after matching is appropriate is beyond the scope of this vignette. You can read more in
the [MatchIt Assessing Balance](https://CRAN.R-project.org/package=MatchIt/vignettes/assessing-balance.html)
vignette. Again the `cobalt` package can be useful here.

However, if you are happy with the results of the matching procedure, you can extract the data for use in `psborrow2`.


``` r
example_matrix_match <- create_data_matrix(
  data = example_dataframe[match_11$weights == 1, ],
  ext_flag_col = "ext",
  outcome = c("time", "cnsr"),
  trt_flag_col = "trt"
)
#> Error:
#> ! object 'match_11' not found
```




``` r
analysis_match <- create_analysis_obj(
  data_matrix = example_matrix_match,
  outcome = outcome_surv_exponential("time", "cnsr", prior_normal(0, 10000)),
  borrowing = borrowing_full("ext"),
  treatment = treatment_details("trt", prior_normal(0, 10000)),
  quiet = TRUE
)
#> Error:
#> ! object 'example_matrix_match' not found

result_match <- mcmc_sample(analysis_match, seed = 123)
#> Error in `h()`:
#> ! error in evaluating the argument 'x' in selecting a method for function 'mcmc_sample': object 'analysis_match' not found
```

## Combined Weighting and Dynamic Borrowing

The models also support fixed weights on the likelihood contributions from each observation. This
is equivalent to fixed power prior weights. This allows for the combination of models, such as an
IPTW + commensurate prior approach.


``` r
analysis_iptw_bdb <- create_analysis_obj(
  data_matrix = example_matrix_att,
  outcome = outcome_surv_exponential("time", "cnsr", prior_normal(0, 10000), weight_var = "att"),
  borrowing = borrowing_hierarchical_commensurate("ext", prior_gamma(0.01, 0.01)),
  treatment = treatment_details("trt", prior_normal(0, 10000)),
  quiet = TRUE
)
#> Error:
#> ! object 'example_matrix_att' not found

result_iptw_bdb <- mcmc_sample(analysis_iptw_bdb, seed = 123)
#> Error in `h()`:
#> ! error in evaluating the argument 'x' in selecting a method for function 'mcmc_sample': object 'analysis_iptw_bdb' not found
```

## Fixed Weights

We can also use weights to specify a fixed power prior model. Here we set the power parameter $\alpha = 0.1$
for the external controls.


``` r
example_matrix_pp01 <- cbind(example_matrix, pp_alpha = ifelse(example_matrix[, "ext"] == 1, 0.1, 1))

analysis_pp01 <- create_analysis_obj(
  data_matrix = example_matrix_pp01,
  outcome = outcome_surv_exponential("time", "cnsr", prior_normal(0, 10000), weight_var = "pp_alpha"),
  borrowing = borrowing_full("ext"),
  treatment = treatment_details("trt", prior_normal(0, 10000)),
  quiet = TRUE
)

result_pp01 <- mcmc_sample(analysis_pp01, seed = 123)
```

## Reference Models

For comparison, we also fit a full borrowing, a no borrowing, and a BDB model without weights.


``` r
result_full <- mcmc_sample(
  create_analysis_obj(
    data_matrix = example_matrix,
    outcome = outcome_surv_exponential("time", "cnsr", prior_normal(0, 10000)),
    borrowing = borrowing_full("ext"),
    treatment = treatment_details("trt", prior_normal(0, 10000)),
    quiet = TRUE
  ),
  seed = 123
)

result_none <- mcmc_sample(
  create_analysis_obj(
    data_matrix = example_matrix,
    outcome = outcome_surv_exponential("time", "cnsr", prior_normal(0, 10000)),
    borrowing = borrowing_none("ext"),
    treatment = treatment_details("trt", prior_normal(0, 10000)),
    quiet = TRUE
  ),
  seed = 123
)

result_bdb <- mcmc_sample(
  create_analysis_obj(
    data_matrix = example_matrix_att,
    outcome = outcome_surv_exponential("time", "cnsr", prior_normal(0, 10000)),
    borrowing = borrowing_hierarchical_commensurate("ext", prior_gamma(0.01, 0.01)),
    treatment = treatment_details("trt", prior_normal(0, 10000)),
    quiet = TRUE
  ),
  seed = 123
)
#> Error in `h()`:
#> ! error in evaluating the argument 'x' in selecting a method for function 'mcmc_sample': object 'example_matrix_att' not found
```




# Comparison of Results


``` r
models <- list(
  "No borrowing" = result_none,
  "Full borrowing" = result_full,
  "Power Prior 01" = result_pp01,
  "ATT Weights" = result_att,
  "IPTW + BDB" = result_iptw_bdb,
  "Matching 1:1" = result_match,
  "BDB" = result_bdb
)
#> Error:
#> ! object 'result_att' not found
```

We can use `summary()` to extract the variable of interest and specify summary statistics.

``` r
results_table <- do.call(rbind, lapply(
  models,
  function(i) i$summary("HR_trt", c("mean", "median", "sd", "quantile2"))
))
#> Error:
#> ! object 'models' not found
knitr::kable(cbind(models = names(models), results_table), digits = 3)
#> Error:
#> ! object 'models' not found
```

We can extract a `draws` object from each model and plot the posterior distribution of the
treatment hazard ratio.

``` r
plot(density(models[[1]]$draws("HR_trt")),
  col = 1, xlim = c(0, 2), ylim = c(0, 9), lwd = 2,
  xlab = "HR_trt",
  main = "Posterior Hazard Ratio"
)
#> Error in `h()`:
#> ! error in evaluating the argument 'x' in selecting a method for function 'plot': object 'models' not found
for (i in 2:7) {
  lines(density(models[[i]]$draws("HR_trt")), col = i, lty = i, lwd = 2)
}
#> Error:
#> ! object 'models' not found
legend("topright", col = seq_along(models), lty = seq_along(models), legend = names(models))
#> Error:
#> ! object 'models' not found
```

Here we see no borrowing and full borrowing at the extremes and the other methods in between.

## References
