Thall et al. (2003) described a method for analysing treatment effects of a common intervention in several sub-types of a single disease. The treatment effects are assumed to be different but exchangeable and correlated. Observing good efficacy in one cohort, for example, increases one’s expectations of efficacy in other cohorts.

They demonstrate the hierarchical model in a trial with binary response outcomes and in another with time-to-event outcomes. This vignette describes the sarcoma example with binary response outcomes. The authors provide WinBUGS code in the appendix of their paper (Thall et al. 2003). We port their model to Stan and illustrate usage with the example given in their paper.

`trialr`

Statistically, the authors assume that in a trial of \(k\) disease subtypes, the treatment effects are drawn from the same common normal distribution

\(\rho_1, ..., \rho_k \sim N(\mu, \sigma^2)\)

As is the convention in BUGS, the authors define normal distributions by a precision parameter \(\tau\) as opposed to the standard deviation parameter \(\sigma\) used here. We have re-specified the model to comply with the Stan convention of using standard deviation. The authors use a normal hyperprior on \(\mu\), and a gamma hyperprior on \(\tau\), equivalent to an inverse gamma hyperprior on \(\sigma^2\).

We observe \(x_i\) responses in
\(n_i\) patients in disease sub-type
\(i\). The rate of response in subtype
\(i\) is modelled as \(p_i = \text{logit}^{-1}(\rho_i)\). Each
\(x_i\) is assumed to be binomially
distributed with success parameter \(p_i\). In Stan, that relationship is
described as `x ~ binomial_logit(n, rho);`

The treatment is judged to be worthy of further investigation in cohort \(i\) if

\(\text{Pr}\left\{ p_i > \theta | \mathcal{D} \right\} > q\)

where \(\theta\) is the minimum required response rate, and \(q\) is the required certainty to approve. In their provided example, Thall et al. (2003) use \(\theta = 0.3\).

`library(trialr)`

We reproduce Thall *etl al.*’s example. We have outcomes in 10
disease subgroups. The number of responses is stored in
`group_responses`

and the number of patients in
`group_sizes`

. There have been 3 / 7 responses in subgroup 4,
for example, but 0 / 2 responses in subgroup 2, and zero patients
treated at all in subgroups 1, 6 and 10.

```
<- stan_hierarchical_response_thall(
fit group_responses = c(0, 0, 1, 3, 5, 0, 1, 2, 0, 0),
group_sizes = c(0, 2 ,1, 7, 5, 0, 2, 3, 1, 0),
mu_mean = -1.3863,
mu_sd = sqrt(1 / 0.1),
tau_alpha = 2,
tau_beta = 20)
```

` fit`

```
## Inference for Stan model: ThallHierarchicalBinary.
## 4 chains, each with iter=2000; warmup=1000; thin=1;
## post-warmup draws per chain=1000, total post-warmup draws=4000.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff
## mu -0.03 0.03 1.37 -2.77 -0.91 -0.02 0.83 2.66 2004
## sigma2 12.42 0.38 11.50 3.33 6.15 9.32 14.60 38.69 909
## rho[1] -0.06 0.07 3.80 -8.09 -2.28 0.03 2.14 7.58 2920
## rho[2] -2.97 0.06 2.50 -9.18 -4.21 -2.51 -1.23 0.66 1645
## rho[3] 2.38 0.06 2.88 -1.93 0.56 1.93 3.67 9.23 1979
## rho[4] -0.31 0.01 0.81 -1.98 -0.82 -0.29 0.21 1.24 3735
## rho[5] 3.72 0.06 2.33 0.58 2.17 3.25 4.74 9.73 1570
## rho[6] 0.08 0.07 3.80 -7.53 -2.20 0.14 2.40 7.42 2857
## rho[7] -0.02 0.02 1.46 -2.97 -0.95 -0.02 0.91 2.88 3596
## rho[8] 0.75 0.02 1.28 -1.68 -0.08 0.72 1.53 3.44 3719
## rho[9] -2.30 0.06 2.61 -8.45 -3.69 -1.98 -0.56 1.96 1991
## rho[10] -0.02 0.07 3.74 -7.77 -2.09 -0.03 2.19 7.21 3193
## sigma 3.31 0.04 1.20 1.83 2.48 3.05 3.82 6.22 880
## prob_response[1] 0.50 0.01 0.38 0.00 0.09 0.51 0.89 1.00 3583
## prob_response[2] 0.15 0.00 0.19 0.00 0.01 0.07 0.23 0.66 3699
## prob_response[3] 0.77 0.00 0.25 0.13 0.64 0.87 0.98 1.00 4504
## prob_response[4] 0.43 0.00 0.17 0.12 0.31 0.43 0.55 0.78 4029
## prob_response[5] 0.92 0.00 0.10 0.64 0.90 0.96 0.99 1.00 4067
## prob_response[6] 0.51 0.01 0.38 0.00 0.10 0.54 0.92 1.00 3542
## prob_response[7] 0.50 0.00 0.26 0.05 0.28 0.49 0.71 0.95 3452
## prob_response[8] 0.64 0.00 0.23 0.16 0.48 0.67 0.82 0.97 4180
## prob_response[9] 0.23 0.00 0.25 0.00 0.02 0.12 0.36 0.88 3673
## prob_response[10] 0.50 0.01 0.38 0.00 0.11 0.49 0.90 1.00 3419
## lp__ -34.21 0.14 3.62 -42.70 -36.17 -33.69 -31.61 -28.89 700
## Rhat
## mu 1
## sigma2 1
## rho[1] 1
## rho[2] 1
## rho[3] 1
## rho[4] 1
## rho[5] 1
## rho[6] 1
## rho[7] 1
## rho[8] 1
## rho[9] 1
## rho[10] 1
## sigma 1
## prob_response[1] 1
## prob_response[2] 1
## prob_response[3] 1
## prob_response[4] 1
## prob_response[5] 1
## prob_response[6] 1
## prob_response[7] 1
## prob_response[8] 1
## prob_response[9] 1
## prob_response[10] 1
## lp__ 1
##
## Samples were drawn using NUTS(diag_e) at Sat Mar 11 15:48:06 2023.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).
```

The probability of response in each subgroup is calculated under
`prob_response`

. We can use the generic summary function from
rstan to get a nice summary:

`::kable(rstan::summary(fit, par = 'prob_response')$summary, digits = 3) knitr`

mean | se_mean | sd | 2.5% | 25% | 50% | 75% | 97.5% | n_eff | Rhat | |
---|---|---|---|---|---|---|---|---|---|---|

prob_response[1] | 0.498 | 0.006 | 0.378 | 0.000 | 0.093 | 0.507 | 0.895 | 0.999 | 3582.526 | 1.000 |

prob_response[2] | 0.151 | 0.003 | 0.185 | 0.000 | 0.015 | 0.075 | 0.226 | 0.659 | 3699.373 | 1.001 |

prob_response[3] | 0.771 | 0.004 | 0.253 | 0.126 | 0.637 | 0.873 | 0.975 | 1.000 | 4503.532 | 0.999 |

prob_response[4] | 0.433 | 0.003 | 0.173 | 0.121 | 0.306 | 0.428 | 0.552 | 0.775 | 4029.447 | 1.000 |

prob_response[5] | 0.925 | 0.002 | 0.097 | 0.640 | 0.898 | 0.962 | 0.991 | 1.000 | 4066.571 | 1.000 |

prob_response[6] | 0.512 | 0.006 | 0.384 | 0.001 | 0.100 | 0.536 | 0.917 | 0.999 | 3542.444 | 1.001 |

prob_response[7] | 0.496 | 0.004 | 0.263 | 0.049 | 0.278 | 0.494 | 0.714 | 0.947 | 3451.631 | 1.002 |

prob_response[8] | 0.637 | 0.003 | 0.225 | 0.157 | 0.480 | 0.673 | 0.821 | 0.969 | 4180.176 | 1.000 |

prob_response[9] | 0.228 | 0.004 | 0.253 | 0.000 | 0.024 | 0.121 | 0.363 | 0.877 | 3673.320 | 0.999 |

prob_response[10] | 0.500 | 0.006 | 0.376 | 0.000 | 0.110 | 0.492 | 0.899 | 0.999 | 3419.327 | 0.999 |

Let us say that we are willing to approve the treatment for further study in a subgroup if we have at least \(q = 70\)% certainty that the probability of efficacy exceeds the target response rate of 30%.

`colMeans(as.data.frame(fit, pars = 'prob_response') > 0.3)`

```
## prob_response[1] prob_response[2] prob_response[3] prob_response[4]
## 0.60175 0.18425 0.92250 0.75900
## prob_response[5] prob_response[6] prob_response[7] prob_response[8]
## 1.00000 0.60850 0.72425 0.90250
## prob_response[9] prob_response[10]
## 0.30525 0.60875
```

On that basis, at this interim stage, we would be eager to approve the treatment in subgroups 3, 5 and 8. Subgroups 4 and 7 are close to the boundary.

Some distribution plots of the probabilities of efficacy in some subgroups may be informative.

```
library(ggplot2)
library(rstan)
library(dplyr)
plot(fit, pars = 'prob_response') +
geom_vline(xintercept = 0.3, col = 'orange', linetype = 'dashed') +
labs(title = 'Partially-pooled analysis of response rate in 10 sarcoma subtypes')
```

We see that the inferred efficacy in subgroup 5 is very high. In contrast, efficacy is subject to more uncertainty in subgroup 4, but the majority of the mass clearly lies to the right of 30%.

`trialr`

is available at https://github.com/brockk/trialr and https://CRAN.R-project.org/package=trialr

Thall, Peter F., J. Kyle Wathen, B. Nebiyou Bekele, Richard E. Champlin,
Laurence H. Baker, and Robert S. Benjamin. 2003. “Hierarchical
Bayesian Approaches to Phase II Trials in
Diseases with Multiple Subtypes.” *Statistics in Medicine*
22 (5): 763–80. https://doi.org/10.1002/sim.1399.