---
title: "Quickstart"
author: Miriam Pedrera Gomez, Isaac Gravestock, and Marcel Wolbers
output: 
  bookdown::html_document2:
    toc: true
    toc_depth: 4
    number_sections: true
    citation_package: natbib
    base_format: rmarkdown::html_vignette
bibliography: "references.bib"
link-citations: true   
linkcolor: blue
pkgdown:
  as_is: true
vignette: >
  %\VignetteIndexEntry{Quickstart}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
    collapse = TRUE,
    comment = "#>"
)

if (isFALSE(requireNamespace("cmdstanr", quietly = TRUE)) &&  packageVersion("rstan") == "2.32.7") {
  knitr::opts_chunk$set(eval = FALSE)
}
```

# Introduction

The `bonsaiforest2` package consists of 3 core functions which are typically called in sequence:

1.  `run_brms_analysis()` - Prepares the model formula and fits the Bayesian shrinkage model using `brms`
2.  `summary_subgroup_effects()` - Derives standardized treatment effect estimates from the fitted model and displays them
3.  `plot()` - Creates a basic forest plot from the summary object

The package supports both **one-way shrinkage models**, which fit a separate Bayesian shrinkage model for each subgrouping variable, and **global shrinkage models**, which fit a single model including all subgrouping variables at once. For both types of models, treatment effect estimates in subgroups are derived via standardization (G-computation) .
The package supports continuous, binary, time-to-event (Cox), and count outcomes. It leverages `brms` and `Stan` for model fitting, allowing for a flexible choice of priors, including normal, regularized horseshoe, and R2-D2.

A comprehensive description of the methodology implemented in this R package can be found in @wolbersUnifiedShrinkage. 

This vignette demonstrates how to use the package to fit and compare these different modeling approaches You'll learn how to:

-   Fit one-way models 
-   Fit global models 
-   Derive standardized treatment effect estimates from the fitted model and displays them
-   Visualize and compare results from different model specifications

It assumes that [`bonsaiforest2`](https://github.com/openpharma/bonsaiforest2), [`brms`](https://paulbuerkner.com/brms/), and a R interface to Stan (e.g.,  [`cmdstanr`](https://mc-stan.org/cmdstanr/)) have been installed.

# The Data

We will use the simulated `shrink_data` dataset which is included in the package. 
For this vignette, we will use the continuous outcome variable `y`, the treatment assignment `trt`, and three categorical subgrouping variables `x1`, `x2`, and `x3`. 

Let's load the libraries and the data.

```{r load-packages, message=FALSE, warning=FALSE}
# Load the main package
library(bonsaiforest2)
library(brms)

# Use cmdstanr by default if available
if (requireNamespace("cmdstanr", quietly = TRUE)) {
  prev_backend <- options(brms.backend = "cmdstanr")
} 
```

```{r ex-data}
# Load the shrink_data package dataset
shrink_data <- bonsaiforest2::shrink_data

print(head(shrink_data))
```

# Choice of hyperprior parameters 

In our example, we anchor parameters of the hyperprior in trial assumptions about the anticipated size of the overall treatment effect and the standard deviation as proposed by @wolbersUnifiedShrinkage. We assume that the data is from a randomized trial which targeted a treatment effect of $\delta_{plan}=0.3$ and assumed a standard deviation of $\sigma_{plan}=1$. 

Specifically, we use the following priors for the shrunken predictive effects in the one-way and global models below:

- *One-way shrinkage model*: a normal prior with a half-normal hyperprior and heterogeneity parameter $\phi=\delta_{plan}=0.3$ for the standard deviation
- *Global shrinkage models*: a regularized horseshoe priors with parameters $\tau_0 = \delta_{plan}=0.3$ (`scale_global`), $s = 2\sigma_{plan} = 2$ (`scale_slab`) and $\nu = 4$ (`df_slab`)

# One-way shrinkage models

A simple one-way shrinkage model for a categorical subgrouping variable `x` is a Bayesian regression model of the form `y ~ 1 + trt + x + trt:x` with a normal shrinkage prior applied to the treatment-by-subgroup interactions. For efficient computation, it is preferable, to reformulate the model as a Bayesian random effects model of the form `y ~ 1 + trt + x + (0+trt || x)` where the term `(0+trt || x)` refers to random treatment effects per subgroups defined by levels of `x`, i.e. treatment-by-subgroup interactions. As discussed, we specify a normal distribution for the random effects representing treatment-by-subgroup interactions with a half-normal hyperprior with heterogeneity parameter $\phi=\delta_{plan}=0.3$ for the standard deviation.  

In the code below, we fit three separate models for `x1`, `x2`, and `x3`, respectively.  We extend the simple model described above by including all three subgrouping variables as unshrunken prognostic terms. While this extension is not mandatory, adjustment for prognostic variables can improve the precision of treatment effect estimates as illustrated in the simulation study for a continuous endpoint reported in @wolbersUnifiedShrinkage. 

## One-way model for `x1`

```{r oneway-x1, message=FALSE, warning=FALSE}
# Fit one-way model using only x1 as a predictive subgrouping variable
# Random effects notation (0 + trt || x1) estimates varying treatment slopes by levels of x1
oneway_x1 <- run_brms_analysis(
  data = shrink_data,
  response_type = "continuous",
  response_formula = y ~ trt,
  unshrunk_terms_formula = ~ x1 + x2 +x3, 
  shrunk_predictive_formula = ~ (0 + trt || x1),
  shrunk_predictive_prior = set_prior("normal(0, 0.3)", class = "sd"),
  chains = 2, iter = 2000, warmup = 1000, cores = 2,
  refresh = 0
)

summary_oneway_x1 <- summary_subgroup_effects(brms_fit = oneway_x1)
print(summary_oneway_x1)
```

## One-way model for `x2`

```{r oneway-x2, message=FALSE, warning=FALSE}
oneway_x2 <- run_brms_analysis(
  data = shrink_data,
  response_type = "continuous",
  response_formula = y ~ trt,
  unshrunk_terms_formula = ~ x1 + x2 +x3,
  shrunk_predictive_formula = ~ (0 + trt || x2),
  shrunk_predictive_prior = set_prior("normal(0, 0.3)", class = "sd"),
  chains = 2, iter = 2000, warmup = 1000, cores = 2,
  refresh = 0
)

summary_oneway_x2 <- summary_subgroup_effects(brms_fit = oneway_x2)
```

## One-way model for `x3`

```{r oneway-x3, message=FALSE, warning=FALSE}
oneway_x3 <- run_brms_analysis(
  data = shrink_data,
  response_type = "continuous",
  response_formula = y ~ trt,
  unshrunk_terms_formula = ~ x1 + x2 +x3,
  shrunk_predictive_formula = ~ (0 + trt || x3),
  shrunk_predictive_prior = set_prior("normal(0, 0.3)", class = "sd"),
  chains = 2, iter = 2000, warmup = 1000, cores = 2,
  refresh = 0
)

summary_oneway_x3 <- summary_subgroup_effects(brms_fit = oneway_x3)
```

## Forest plot of results from all one-way models

You can combine and visualize results from multiple models using `combine_summaries()`:

```{r compare-all-oneway, fig.width=9, fig.height=8}
# Combine all one-way models
combined_oneway <- combine_summaries(list(
  "x1" = summary_oneway_x1,
  "x2" = summary_oneway_x2,
  "x3" = summary_oneway_x3
))

plot(combined_oneway, title = "One-way Models: All Subgrouping Variables")
```

# Global shrinkage model

A simple global shrinkage model for our setting has the form ` y~ 1 + trt + x1 + x2 + x3 + trt:x1 + trt:x2 + trt:x3` where a shrinkage prior is applied to all treatment-by-subgroup interaction terms. As discussed, we use a regularized horseshoe prior with parameters $\tau_0 = \delta_{plan}=0.3$ (`scale_global`), $s = 2\sigma_{plan} = 2$ (`scale_slab`) and $\nu = 4$ (`df_slab`) in the example below. Standardized treatment effects in subgroup defined by the levels of a single subgrouping variable at a time are subsequently derived from this model via G-computation. 

## Model fitting

```{r global-fit, message=FALSE, warning=FALSE}
# Fit a single unified model with ALL subgrouping variables simultaneously using global approach
global_shrinkage_model <- run_brms_analysis(
  data = shrink_data,
  response_type = "continuous",
  response_formula = y ~ trt,
  unshrunk_terms_formula  = ~ 1 + x1 + x2 + x3,
  shrunk_predictive_formula = ~ 0 + trt:x1 + trt:x2 + trt:x3,
  shrunk_predictive_prior = "horseshoe(scale_global=0.3, scale_slab = 2, df_slab = 4, autoscale = FALSE)",
  chains = 2, iter = 2000, warmup = 1000, cores = 2,
  refresh = 0
)
```

## Summary of subgroup effects

The function `summary_subgroup_effects()` derives standardized treatment effects in subgroups from the global model. Internally, subgrouping variables are identified as all terms describing treatment-by-covariate interactions provided in the formulas. 

```{r global-summary}
global_summary <- summary_subgroup_effects(brms_fit = global_shrinkage_model)

# Print the summary of subgroup-specific treatment effects
print(global_summary)
```

## Global Model: Visualization

The `plot()` function creates  a basic forest plot from the summary object:

```{r global-plot, fig.width=8, fig.height=10}
plot(global_summary, title = "Global Model: All subgrouping variables")
```

# Comparing Multiple Models in One Plot

The `plot()` function supports comparing multiple models side-by-side by passing a in a forest plot by passing a named list of `subgroup_summary` objects. 

```{r compare-models, fig.width=9, fig.height=7}
# Combine summaries for comparison
combined <- combine_summaries(list(
  "One-way" = combined_oneway,
  "Global" = global_summary
))

# Plot the comparison
plot(combined, title = "Comparing One-way vs Global Models")
```

```{r}
# reset stan backend
if (requireNamespace("cmdstanr", quietly = TRUE)) {
  options(prev_backend)
}
```

# References