Introduction

In patient-centered outcomes research, it is vital to assess the heterogeneity of treatment effects (HTE) when making health care decisions for an individual patient or a group of patients. Nevertheless, it remains challenging to evaluate HTE based on information collected from clinical studies that are often designed and conducted to evaluate the efficacy of a treatment for the overall population. The Bayesian framework offers a principled and flexible approach to estimate and compare treatment effects across subgroups of patients defined by their characteristics.

R package beanz provides functions to facilitate the conduct of Bayesian analysis of HTE and a web-based graphical user interface for users to conduct such Bayesian analysis in an interactive and user-friendly manner.

Data accepted by beanz

There are two types of data structures that beanz recognizes:

Summary treatment effect data: Each row should represent a subgroup with covariates that define the subgroup, estimated treatment effect in the subgroup and variance for the estimation.
Patient level raw data: Each row should represent a patient with covariates that define the subgroup in which the patient belongs to, treatment indicator and outcome. The outcome can be binary, continuous, or time to event.

The beanz package provides dataset solvd.sub from the SOLVD trial as an example Patient level raw data dataset.

Estimate subgroup effect

If Patient level raw data is provided, the package provides function bzGetSubgrpRaw for estimating subgroup effect for each subgroup. The return value from bzGetSubgrpRaw is a data frame with the format of Summary treatment effect data.

The example is as follows:

var.cov    <- c("lvef", "sodium", "any.vasodilator.use");
var.resp   <- "y";
var.trt    <- "trt";
var.censor <- "censor";
resptype   <- "survival";

subgrp.effect <- bzGetSubgrpRaw(solvd.sub,
                                  var.resp   = var.resp,
                                  var.trt    = var.trt,
                                  var.cov    = var.cov,
                                  var.censor = var.censor,
                                  resptype   = resptype);
print(subgrp.effect);

##   Subgroup lvef sodium any.vasodilator.use    Estimate   Variance   N
## 1        1    0      0                   0 -0.37783038 0.01212786 562
## 2        2    0      0                   1 -0.34655336 0.01004499 695
## 3        3    0      1                   0 -0.79235451 0.03939983 237
## 4        4    0      1                   1 -0.39334304 0.02969421 250
## 5        5    1      0                   0  0.06776454 0.04629163 223
## 6        6    1      0                   1 -0.23655764 0.02400353 341
## 7        7    1      1                   0  0.15435495 0.10365396 104
## 8        8    1      1                   1  0.05947290 0.07761840 123

Bayesian HTE models

The function bzCallStan calls rstan::sampling to draw samples for different Bayesian models. The following models are available in the current version of beanz:

nse: No subgroup effect model
fs: Full stratification model
sr: Simple regression model
bs: Basic shrinkage model
srs: Simple regression with shrinkage model
ds: Dixon-Simon model
eds: Extended Dixon-Simon model.

The following examples show how No subgroup effect model (nse), Simple regression model* (sr) and Basic shrinkage model (bs) are called:

var.estvar <- c("Estimate", "Variance");

rst.nse <- bzCallStan("nse", dat.sub=subgrp.effect,
                     var.estvar = var.estvar, var.cov = var.cov,
                     par.pri = c(B=1000),
                     chains=4, iter=4000,
                     warmup=2000, seed=1000, cores=1);

## 
## SAMPLING FOR MODEL 'nse' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 3.7e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.37 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 1: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 1: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 1: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 1: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 1: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 1: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 1: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 1: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 1: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 1: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 1: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 0.137673 seconds (Warm-up)
## Chain 1:                0.089251 seconds (Sampling)
## Chain 1:                0.226924 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL 'nse' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 1e-05 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.1 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 2: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 2: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 2: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 2: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 2: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 2: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 2: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 2: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 2: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 2: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 2: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 0.136251 seconds (Warm-up)
## Chain 2:                0.089365 seconds (Sampling)
## Chain 2:                0.225616 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL 'nse' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 1.1e-05 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.11 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 3: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 3: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 3: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 3: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 3: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 3: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 3: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 3: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 3: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 3: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 3: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 0.142632 seconds (Warm-up)
## Chain 3:                0.089981 seconds (Sampling)
## Chain 3:                0.232613 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL 'nse' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 1e-05 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.1 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 4: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 4: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 4: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 4: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 4: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 4: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 4: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 4: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 4: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 4: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 4: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 0.141991 seconds (Warm-up)
## Chain 4:                0.091796 seconds (Sampling)
## Chain 4:                0.233787 seconds (Total)
## Chain 4:

## Warning: Relative effective sample sizes ('r_eff' argument) not specified.
## For models fit with MCMC, the reported PSIS effective sample sizes and 
## MCSE estimates will be over-optimistic.

rst.sr  <- bzCallStan("sr", dat.sub=subgrp.effect,
                     var.estvar = var.estvar, var.cov = var.cov,
                     par.pri = c(B=1000, C=1000),
                     chains=4, iter=4000,
                     warmup=2000,  seed=1000, cores=1);

## 
## SAMPLING FOR MODEL 'sr' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 4.3e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.43 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 1: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 1: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 1: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 1: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 1: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 1: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 1: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 1: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 1: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 1: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 1: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 0.212377 seconds (Warm-up)
## Chain 1:                0.18307 seconds (Sampling)
## Chain 1:                0.395447 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL 'sr' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 1.5e-05 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.15 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 2: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 2: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 2: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 2: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 2: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 2: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 2: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 2: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 2: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 2: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 2: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 0.217548 seconds (Warm-up)
## Chain 2:                0.130562 seconds (Sampling)
## Chain 2:                0.34811 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL 'sr' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 1.2e-05 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.12 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 3: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 3: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 3: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 3: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 3: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 3: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 3: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 3: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 3: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 3: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 3: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 0.206407 seconds (Warm-up)
## Chain 3:                0.158516 seconds (Sampling)
## Chain 3:                0.364923 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL 'sr' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 1.6e-05 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.16 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 4: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 4: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 4: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 4: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 4: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 4: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 4: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 4: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 4: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 4: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 4: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 0.209421 seconds (Warm-up)
## Chain 4:                0.176035 seconds (Sampling)
## Chain 4:                0.385456 seconds (Total)
## Chain 4:

## Warning: Relative effective sample sizes ('r_eff' argument) not specified.
## For models fit with MCMC, the reported PSIS effective sample sizes and 
## MCSE estimates will be over-optimistic.

## Warning: Some Pareto k diagnostic values are too high. See help('pareto-k-diagnostic') for details.

rst.bs  <- bzCallStan("bs", dat.sub=subgrp.effect,
                     var.estvar = var.estvar, var.cov = var.cov,
                     par.pri = c(B=1000, D=1),
                     chains=4, iter=4000, warmup=2000,  seed=1000, cores=1);

## 
## SAMPLING FOR MODEL 'bs' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 2.8e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.28 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 1: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 1: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 1: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 1: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 1: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 1: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 1: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 1: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 1: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 1: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 1: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 0.316251 seconds (Warm-up)
## Chain 1:                0.268294 seconds (Sampling)
## Chain 1:                0.584545 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL 'bs' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 1.1e-05 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.11 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 2: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 2: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 2: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 2: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 2: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 2: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 2: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 2: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 2: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 2: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 2: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 0.312 seconds (Warm-up)
## Chain 2:                0.220448 seconds (Sampling)
## Chain 2:                0.532448 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL 'bs' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 1.2e-05 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.12 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 3: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 3: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 3: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 3: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 3: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 3: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 3: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 3: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 3: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 3: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 3: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 0.340872 seconds (Warm-up)
## Chain 3:                0.220584 seconds (Sampling)
## Chain 3:                0.561456 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL 'bs' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 1.2e-05 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.12 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:    1 / 4000 [  0%]  (Warmup)
## Chain 4: Iteration:  400 / 4000 [ 10%]  (Warmup)
## Chain 4: Iteration:  800 / 4000 [ 20%]  (Warmup)
## Chain 4: Iteration: 1200 / 4000 [ 30%]  (Warmup)
## Chain 4: Iteration: 1600 / 4000 [ 40%]  (Warmup)
## Chain 4: Iteration: 2000 / 4000 [ 50%]  (Warmup)
## Chain 4: Iteration: 2001 / 4000 [ 50%]  (Sampling)
## Chain 4: Iteration: 2400 / 4000 [ 60%]  (Sampling)
## Chain 4: Iteration: 2800 / 4000 [ 70%]  (Sampling)
## Chain 4: Iteration: 3200 / 4000 [ 80%]  (Sampling)
## Chain 4: Iteration: 3600 / 4000 [ 90%]  (Sampling)
## Chain 4: Iteration: 4000 / 4000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 0.305328 seconds (Warm-up)
## Chain 4:                0.30374 seconds (Sampling)
## Chain 4:                0.609068 seconds (Total)
## Chain 4:

## Warning: There were 8 divergent transitions after warmup. See
## https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## to find out why this is a problem and how to eliminate them.

## Warning: Examine the pairs() plot to diagnose sampling problems

## Warning: Relative effective sample sizes ('r_eff' argument) not specified.
## For models fit with MCMC, the reported PSIS effective sample sizes and 
## MCSE estimates will be over-optimistic.

## Warning: Some Pareto k diagnostic values are too high. See help('pareto-k-diagnostic') for details.

Results presentation

Posterior subgroup treatment effect summary

Posterior subgroup treatment effect can be summarized and presented by functions bzSummary, bzPlot and bzForest. These functions allows to include a subgroup from another model (i.e. No subgroup effect model) as a reference in the results.

Simple regression model

sel.grps <- c(1,4,5);
tbl.sub <- bzSummary(rst.sr, ref.stan.rst=rst.nse, ref.sel.grps=1);
print(tbl.sub);

##                Subgroup   Mean    SD   Q025    Q25 Median    Q75   Q975 ProbLT0
## 1            Subgroup 1 -0.401 0.094 -0.589 -0.463 -0.401 -0.339 -0.219       1
## 2            Subgroup 2  -0.38 0.088 -0.549 -0.439  -0.38 -0.321 -0.205       1
## 3            Subgroup 3 -0.489 0.134 -0.757 -0.579 -0.488 -0.399 -0.223       1
## 4            Subgroup 4 -0.467 0.127 -0.716 -0.553 -0.467 -0.382 -0.222       1
## 5            Subgroup 5 -0.063 0.134  -0.32 -0.154 -0.062  0.029  0.197    0.68
## 6            Subgroup 6 -0.041 0.121 -0.277 -0.122 -0.039   0.04  0.192   0.632
## 7            Subgroup 7  -0.15 0.161 -0.465 -0.259 -0.152 -0.042  0.163   0.828
## 8            Subgroup 8 -0.129 0.147 -0.415 -0.227 -0.131 -0.032  0.167   0.813
## 9 No subgroup effect(1) -0.322 0.055 -0.432 -0.359 -0.321 -0.285 -0.215       1

bzPlot(rst.sr, sel.grps = sel.grps, ref.stan.rst=rst.nse, ref.sel.grps=1);

bzForest(rst.sr, sel.grps = sel.grps, ref.stan.rst=rst.nse, ref.sel.grps=1);

Basic shrinkage model

tbl.sub <- bzSummary(rst.bs, ref.stan.rst=rst.nse, ref.sel.grps=1);
print(tbl.sub);

##                Subgroup   Mean    SD   Q025    Q25 Median    Q75   Q975 ProbLT0
## 1            Subgroup 1 -0.348 0.094 -0.544 -0.409 -0.345 -0.284 -0.171       1
## 2            Subgroup 2 -0.333 0.087 -0.509 -0.389 -0.332 -0.277 -0.165       1
## 3            Subgroup 3 -0.514 0.187 -0.929 -0.639  -0.49 -0.368 -0.225       1
## 4            Subgroup 4 -0.346 0.129 -0.621 -0.421 -0.341 -0.264 -0.094   0.996
## 5            Subgroup 5  -0.15 0.181 -0.427 -0.287 -0.177 -0.036  0.251   0.795
## 6            Subgroup 6 -0.266 0.122 -0.489 -0.348 -0.275  -0.19 -0.004   0.977
## 7            Subgroup 7 -0.172 0.212   -0.5 -0.318 -0.212 -0.056  0.328    0.81
## 8            Subgroup 8 -0.182 0.192 -0.483 -0.318 -0.215  -0.07  0.265   0.832
## 9 No subgroup effect(1) -0.322 0.055 -0.432 -0.359 -0.321 -0.285 -0.215       1

bzPlot(rst.bs, sel.grps = sel.grps, ref.stan.rst=rst.nse, ref.sel.grps=1);

bzForest(rst.bs, sel.grps = sel.grps, ref.stan.rst=rst.nse, ref.sel.grps=1);

Posterior subgroup treatment effect comparison

Posterior subgroup treatment effect can be compared between subgroups by functions bzSummaryComp, bzPlotComp and bzForestComp.

Simple regression model

tbl.sub <- bzSummaryComp(rst.sr, sel.grps=sel.grps);
print(tbl.sub);

##     Comparison   Mean    SD   Q025    Q25 Median   Q75  Q975 ProbLT0
## 1 Subgroup 4-1 -0.067 0.159 -0.374 -0.173 -0.067  0.04 0.249   0.668
## 2 Subgroup 5-1  0.337 0.166  0.014  0.224  0.337 0.449 0.665    0.02
## 3 Subgroup 5-4  0.403 0.184  0.053   0.28  0.403 0.526 0.768   0.011

bzPlot(rst.sr, sel.grps = sel.grps);

bzForest(rst.sr, sel.grps = sel.grps);

Basic shrinkage model

tbl.sub <- bzSummaryComp(rst.bs, sel.grps=sel.grps);
print(tbl.sub);

##     Comparison  Mean    SD   Q025    Q25 Median   Q75  Q975 ProbLT0
## 1 Subgroup 4-1     0 0.162 -0.327 -0.105  0.002 0.105 0.317   0.496
## 2 Subgroup 5-1 0.201 0.205 -0.148  0.049   0.18 0.334 0.642   0.163
## 3 Subgroup 5-4 0.197 0.224 -0.195  0.037  0.182 0.347 0.672   0.192

bzPlotComp(rst.bs, sel.grps = sel.grps);

bzForestComp(rst.bs, sel.grps = sel.grps);

Overall summary

beanz provides function bzRptTbl to generate the summary posterior subgroup treatment effect table from the model selected by DIC (i.e. the model with the smallest DIC):

lst.rst     <- list(nse=rst.nse, sr=rst.sr, bs=rst.bs);
tbl.summary <- bzRptTbl(lst.rst, dat.sub = subgrp.effect, var.cov = var.cov);

## Warning: Accessing looic using '$' is deprecated and will be removed in a
## future release. Please extract the looic estimate from the 'estimates'
## component instead.

## Warning: Accessing looic using '$' is deprecated and will be removed in a
## future release. Please extract the looic estimate from the 'estimates'
## component instead.

## Warning: Accessing looic using '$' is deprecated and will be removed in a
## future release. Please extract the looic estimate from the 'estimates'
## component instead.

print(tbl.summary);

##                         Model Subgroup lvef sodium any.vasodilator.use   Mean
## Subgroup 1 No subgroup effect        1    0      0                   0 -0.322
## Subgroup 2 No subgroup effect        2    0      0                   1 -0.322
## Subgroup 3 No subgroup effect        3    0      1                   0 -0.322
## Subgroup 4 No subgroup effect        4    0      1                   1 -0.322
## Subgroup 5 No subgroup effect        5    1      0                   0 -0.322
## Subgroup 6 No subgroup effect        6    1      0                   1 -0.322
## Subgroup 7 No subgroup effect        7    1      1                   0 -0.322
## Subgroup 8 No subgroup effect        8    1      1                   1 -0.322
##               SD Prob < 0
## Subgroup 1 0.055        1
## Subgroup 2 0.055        1
## Subgroup 3 0.055        1
## Subgroup 4 0.055        1
## Subgroup 5 0.055        1
## Subgroup 6 0.055        1
## Subgroup 7 0.055        1
## Subgroup 8 0.055        1

Predictive distribution

Function bzPredSubgrp generates the predictive distribution of the subgrooup treatment effects.

pred.dist <- bzPredSubgrp(rst.sr,
                                  dat.sub=subgrp.effect,
                                  var.estvar = var.estvar);
head(pred.dist);

##            [,1]       [,2]       [,3]        [,4]        [,5]        [,6]
## [1,] -0.5679901 -0.4997621 -0.4687660 -0.76574757  0.26725847  0.01626289
## [2,] -0.3420504 -0.2521794 -0.7345765 -0.47152359  0.36613344 -0.07645110
## [3,] -0.4338145 -0.2525592 -0.7803938 -0.54795933 -0.01929658 -0.03668858
## [4,] -0.4814998 -0.5537714 -0.6324113 -0.39108840 -0.15506442 -0.12222381
## [5,] -0.4676442 -0.4078786 -0.6636981 -0.51484721 -0.30893251 -0.54646626
## [6,] -0.3273695 -0.2048737 -0.5903188 -0.09467884 -0.34219340 -0.12112388
##             [,7]         [,8]
## [1,]  0.12626214  0.002695502
## [2,]  0.46222121  0.081932323
## [3,] -0.74327267 -0.875854377
## [4,]  0.03560672  0.246599939
## [5,] -0.79005015 -0.039686095
## [6,] -0.12416157 -0.076585617

Graphical User Interface

With package shiny installed, beaz provides a web-based graphical user interface (GUI) for conducting the HTE analysis in an user-friendly interactive manner. The GUI can be started by

bzShiny();

Toolbox

Package beanz provides function bzGailSimon that implements the Gail-Simon test for qualitative interactions:

gs.pval <- bzGailSimon(subgrp.effect$Estimate,
                       sqrt(subgrp.effect$Variance));
print(gs.pval);

## [1] 0.9191656

The result show that there is no significant qualitative interactions according to the Gail-Simon test.

beanz: Bayesian Analysis of Heterogeneous Treatment Effect

Chenguang Wang and Ravi Varadhan

2023-08-07

Introduction

Data accepted by beanz

Estimate subgroup effect

Bayesian HTE models

Results presentation

Posterior subgroup treatment effect summary

Simple regression model

Basic shrinkage model

Posterior subgroup treatment effect comparison

Simple regression model

Basic shrinkage model

Overall summary

Predictive distribution

Graphical User Interface

Toolbox