# beanz: Bayesian Analysis of Heterogeneous Treatment Effect

#### Chenguang Wang and Ravi Varadhan

#### 2016-11-13

    ## Loading required package: beanz

    ## Loading required package: Rcpp

# 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 r.get.subgrp.raw for estimating subgroup effect for each subgroup. The return value from r.get.subgrp.raw is a data frame with the format of Summary treatment effect data.

The example is as follows:

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

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

    ##   Subgroup sodium lvef 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.06776454 0.04629163 223
    ## 4        4      0    1                   1 -0.23655764 0.02400353 341
    ## 5        5      1    0                   0 -0.79235451 0.03939983 237
    ## 6        6      1    0                   1 -0.39334304 0.02969421 250
    ## 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 call.stan 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.cov    <- c("sodium", "lvef", "any.vasodilator.use");
    var.estvar <- c("Estimate", "Variance");

    rst.nse <- call.stan("nse", dat.sub=subgrp.effect,
                         var.estvar = var.estvar, var.cov = var.cov,
                         lst.par.pri = list(vtau=1000, vrange=c(0,0)),
                         chains=1, iter=4000,
                         warmup=2000, thin=2, seed=1000);

    ##
    ## SAMPLING FOR MODEL 'nse' NOW (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)
    ##  Elapsed Time: 0.107619 seconds (Warm-up)
    ##                0.068626 seconds (Sampling)
    ##                0.176245 seconds (Total)

    rst.sr  <- call.stan("sr", dat.sub=subgrp.effect,
                         var.estvar = var.estvar, var.cov = var.cov,
                         lst.par.pri = list(vtau=1000, vgamma=1000, vrange=c(0,0)),
                         chains=1, iter=4000,
                         warmup=2000, thin=2, seed=1000);

    ##
    ## SAMPLING FOR MODEL 'sr' NOW (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)
    ##  Elapsed Time: 0.209661 seconds (Warm-up)
    ##                0.115246 seconds (Sampling)
    ##                0.324907 seconds (Total)

    rst.bs  <- call.stan("bs", dat.sub=subgrp.effect,
                         var.estvar = var.estvar, var.cov = var.cov,
                         lst.par.pri = list(vtau=1000, vw=100, vrange=c(-0.1,0.1)),
                         chains=1, iter=4000, warmup=2000, thin=2, seed=1000);

    ##
    ## SAMPLING FOR MODEL 'bs' NOW (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)
    ##  Elapsed Time: 0.186847 seconds (Warm-up)
    ##                0.175148 seconds (Sampling)
    ##                0.361995 seconds (Total)

# Results presentation

## Posterior subgroup treatment effect summary

Posterior subgroup treatment effect can be summarized and presented by functions r.summary.stan, r.plot.stan and r.forest.stan. 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 <- r.summary.stan(rst.sr, ref.stan.rst=rst.nse, ref.sel.grps=1);
    print(tbl.sub);

    ##                       Subgroup                Mean     SD      2.5%
    ## Subgroup 1            "Subgroup 1"            "-0.398" "0.091" "-0.571"
    ## Subgroup 2            "Subgroup 2"            "-0.38"  "0.082" "-0.537"
    ## Subgroup 3            "Subgroup 3"            "-0.062" "0.128" "-0.313"
    ## Subgroup 4            "Subgroup 4"            "-0.044" "0.117" "-0.278"
    ## Subgroup 5            "Subgroup 5"            "-0.482" "0.129" "-0.73"
    ## Subgroup 6            "Subgroup 6"            "-0.464" "0.122" "-0.698"
    ## Subgroup 7            "Subgroup 7"            "-0.146" "0.156" "-0.432"
    ## Subgroup 8            "Subgroup 8"            "-0.128" "0.145" "-0.408"
    ## No subgroup effect(1) "No subgroup effect(1)" "-0.322" "0.054" "-0.433"
    ##                       25%      Median   75%      97.5%    Prob < 0
    ## Subgroup 1            "-0.461" "-0.397" "-0.34"  "-0.214" "1"
    ## Subgroup 2            "-0.435" "-0.383" "-0.325" "-0.213" "1"
    ## Subgroup 3            "-0.15"  "-0.059" "0.024"  "0.187"  "0.683"
    ## Subgroup 4            "-0.125" "-0.04"  "0.033"  "0.189"  "0.639"
    ## Subgroup 5            "-0.571" "-0.484" "-0.402" "-0.218" "1"
    ## Subgroup 6            "-0.548" "-0.467" "-0.383" "-0.221" "1"
    ## Subgroup 7            "-0.252" "-0.149" "-0.038" "0.163"  "0.822"
    ## Subgroup 8            "-0.225" "-0.127" "-0.03"  "0.159"  "0.811"
    ## No subgroup effect(1) "-0.358" "-0.32"  "-0.287" "-0.214" "1"

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

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

### Basic shrinkage model

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

    ##                       Subgroup                Mean     SD      2.5%
    ## Subgroup 1            "Subgroup 1"            "-0.366" "0.101" "-0.563"
    ## Subgroup 2            "Subgroup 2"            "-0.333" "0.092" "-0.509"
    ## Subgroup 3            "Subgroup 3"            "-0.061" "0.198" "-0.403"
    ## Subgroup 4            "Subgroup 4"            "-0.252" "0.138" "-0.5"
    ## Subgroup 5            "Subgroup 5"            "-0.61"  "0.194" "-1.006"
    ## Subgroup 6            "Subgroup 6"            "-0.363" "0.153" "-0.664"
    ## Subgroup 7            "Subgroup 7"            "-0.079" "0.254" "-0.533"
    ## Subgroup 8            "Subgroup 8"            "-0.097" "0.228" "-0.485"
    ## No subgroup effect(1) "No subgroup effect(1)" "-0.322" "0.054" "-0.433"
    ##                       25%      Median   75%      97.5%    Prob < 0
    ## Subgroup 1            "-0.435" "-0.366" "-0.298" "-0.171" "0.999"
    ## Subgroup 2            "-0.399" "-0.331" "-0.269" "-0.161" "1"
    ## Subgroup 3            "-0.205" "-0.073" "0.07"   "0.355"  "0.639"
    ## Subgroup 4            "-0.351" "-0.253" "-0.165" "0.052"  "0.961"
    ## Subgroup 5            "-0.738" "-0.605" "-0.471" "-0.266" "1"
    ## Subgroup 6            "-0.466" "-0.358" "-0.266" "-0.063" "0.99"
    ## Subgroup 7            "-0.259" "-0.093" "0.079"  "0.477"  "0.647"
    ## Subgroup 8            "-0.268" "-0.117" "0.052"  "0.372"  "0.68"
    ## No subgroup effect(1) "-0.358" "-0.32"  "-0.287" "-0.214" "1"

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

    r.forest.stan(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 r.summary.comp, r.plot.comp and r.forest.comp.

### Simple regression model

    tbl.sub <- r.summary.comp(rst.sr, sel.grps=sel.grps);
    print(tbl.sub);

    ##              Comparison     Mean     SD      2.5%     25%      Median
    ## Subgroup 4-1 "Subgroup 4-1" "0.356"  "0.148" "0.077"  "0.256"  "0.354"
    ## Subgroup 5-1 "Subgroup 5-1" "-0.08"  "0.162" "-0.406" "-0.183" "-0.085"
    ## Subgroup 5-4 "Subgroup 5-4" "-0.436" "0.171" "-0.76"  "-0.548" "-0.437"
    ##              75%      97.5%    Prob <0
    ## Subgroup 4-1 "0.455"  "0.648"  "0.01"
    ## Subgroup 5-1 "0.027"  "0.235"  "0.707"
    ## Subgroup 5-4 "-0.324" "-0.075" "0.992"

    r.plot.stan(rst.sr, sel.grps = sel.grps);

    r.forest.stan(rst.sr, sel.grps = sel.grps);

### Basic shrinkage model

    tbl.sub <- r.summary.comp(rst.bs, sel.grps=sel.grps);
    print(tbl.sub);

    ##              Comparison     Mean     SD      2.5%     25%      Median
    ## Subgroup 4-1 "Subgroup 4-1" "0.115"  "0.173" "-0.225" "0.005"  "0.111"
    ## Subgroup 5-1 "Subgroup 5-1" "-0.236" "0.223" "-0.71"  "-0.386" "-0.229"
    ## Subgroup 5-4 "Subgroup 5-4" "-0.355" "0.231" "-0.823" "-0.514" "-0.35"
    ##              75%      97.5%   Prob <0
    ## Subgroup 4-1 "0.229"  "0.474" "0.242"
    ## Subgroup 5-1 "-0.066" "0.165" "0.859"
    ## Subgroup 5-4 "-0.194" "0.076" "0.937"

    r.plot.comp(rst.bs, sel.grps = sel.grps);

    r.forest.comp(rst.bs, sel.grps = sel.grps);

## Overall summary

beanz provides function r.rpt.tbl 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 <- r.rpt.tbl(lst.rst, dat.sub = subgrp.effect, var.cov = var.cov);
    print(tbl.summary);

    ##            Subgroup sodium lvef any.vasodilator.use   Mean    SD Prob < 0
    ## Subgroup 1        1      0    0                   0 -0.366 0.101    0.999
    ## Subgroup 2        2      0    0                   1 -0.333 0.092    1.000
    ## Subgroup 3        3      0    1                   0 -0.061 0.198    0.639
    ## Subgroup 4        4      0    1                   1 -0.252 0.138    0.961
    ## Subgroup 5        5      1    0                   0 -0.610 0.194    1.000
    ## Subgroup 6        6      1    0                   1 -0.363 0.153    0.990
    ## Subgroup 7        7      1    1                   0 -0.079 0.254    0.647
    ## Subgroup 8        8      1    1                   1 -0.097 0.228    0.680

# Predictive distribution

Function r.pred.subgrp.effect generates the predictive distribution of the subgrooup treatment effects.

    pred.dist <- r.pred.subgrp.effect(rst.sr,
                                      dat.sub=subgrp.effect,
                                      var.estvar = var.estvar,
                                      vrange = c(0,0));
    head(pred.dist);

    ##
    ## iterations       [,1]       [,2]        [,3]         [,4]       [,5]
    ##       [1,] -0.2176025 -0.3856852  0.06607013 -0.139958161 -0.4460334
    ##       [2,] -0.4170360 -0.3146029 -0.17863836 -0.011249011 -0.7085243
    ##       [3,] -0.6200378 -0.5654614  0.06944788  0.014510103 -0.7175829
    ##       [4,] -0.4422932 -0.2139824 -0.27109531 -0.007374267 -0.2818747
    ##       [5,] -0.4339711 -0.2571437 -0.24024606 -0.040295363  0.1101854
    ##       [6,] -0.5366364 -0.2758277 -0.10846071 -0.253520355 -0.4519036
    ##
    ## iterations        [,6]        [,7]        [,8]
    ##       [1,] -0.85429190 -0.01060591 -0.72668000
    ##       [2,] -0.41420544  0.19974309 -0.02774356
    ##       [3,] -0.52636155 -0.01108116 -0.34838555
    ##       [4,] -0.02455673 -0.35312284  0.36870280
    ##       [5,] -0.25490106 -0.21428155  0.11761253
    ##       [6,] -0.35239622 -0.56282302 -0.08960953

# 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

    run.beanz();

# Toolbox

Package beanz provides function r.gailsimon that implements the Gail-Simon test for qualitative interactions:

    gs.pval <- r.gailsimon(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.