Two binary endpoints • BayesianQDM

This vignette describes the analysis of a proof-of-concept (PoC) study with two binary endpoints (Yamaguchi et al., submitted). vignette("overview") introduces the decision criteria for one endpoint.

Model

The responses of a subject on the two endpoints fall into one of four cells, "00", "01", "10" and "11", where the first digit is the response on endpoint 1 and the second the response on endpoint 2. In group $k$ , the counts $x_k$ of the four cells follow a multinomial distribution with cell probabilities $p_k = (p_{k,00}, p_{k,01}, p_{k,10}, p_{k,11})$ , which have a Dirichlet prior distribution $\textrm{Dir}(\alpha_k)$ (prior_dirichlet()). The posterior distribution is $\textrm{Dir}(\alpha_k + x_k)$ , and with external data $x_{e}$ borrowed with weight $a_0$ it is $\textrm{Dir}(\alpha_k + a_0 x_{e} + x_k)$ .

The response rates of the two endpoints are the margins $\pi_{k1} = p_{k,10} + p_{k,11}$ and $\pi_{k2} = p_{k,01} + p_{k,11}$ , and the treatment effects are $\theta_j = \pi_{tj} - \pi_{cj}$ ( $j = 1, 2$ ).

Regions and criteria

The target values and minimum acceptable values divide the plane of $(\theta_1, \theta_2)$ into nine regions:

qdm_regions("posterior")
#>                      theta1 > TV1 MAV1 < theta1 <= TV1 theta1 <= MAV1
#> theta2 > TV2         R1           R4                   R7            
#> MAV2 < theta2 <= TV2 R2           R5                   R8            
#> theta2 <= MAV2       R3           R6                   R9

With rule = "joint" (the default), the Go criterion is that the posterior probability of the union of go_regions is at least $\gamma_{go}$ , and the NoGo criterion that the probability of the union of nogo_regions is at least $\gamma_{nogo}$ . The defaults, region 1 for Go and region 9 for NoGo, require both effects to exceed their target values and both to be at most their minimum acceptable values. With rule = "and" or rule = "or", the criteria are applied to each endpoint separately and combined with AND or OR.

For the predictive criterion, the null values divide the plane of the future effects into four regions:

qdm_regions("predictive")
#>                 theta1 > NULL1 theta1 <= NULL1
#> theta2 > NULL2  R1             R3             
#> theta2 <= NULL2 R2             R4

The posterior probabilities of the regions are estimated with nMC draws from the Dirichlet posterior distributions. The predictive probabilities are computed exactly from the joint distribution of the numbers of future responders on the two endpoints.

Decision for observed data

The first case study of Yamaguchi et al. (submitted) has 10 subjects per group, $\textrm{Dir}(0.25, 0.25, 0.25, 0.25)$ priors, $\theta_{TV} = (0.15, 0.15)$ and $\theta_{MAV} = (0.10, 0.10)$ . Here the thresholds are $\gamma_{go} = \gamma_{nogo} = 0.6$ (the article calibrates its own thresholds), and nMC = 2000 keeps the computations quick:

prior <- prior_dirichlet(rep(0.25, 4))
des <- design_binary2(
  n_t = 10, n_c = 10, prior_t = prior, prior_c = prior,
  criterion = criterion_posterior(theta_tv = c(0.15, 0.15),
                                  theta_mav = c(0.10, 0.10),
                                  gamma_go = 0.6, gamma_nogo = 0.6),
  nMC = 2000
)
dec <- qdm_decide(des, data_t = data_binary2(c(3, 2, 2, 3)),
                  data_c = data_binary2(c(6, 2, 1, 1)), seed = 1)
dec
#> Decision: Gray
#>   Go:   P(regions {1} | data) = 0.4445 (threshold 0.6, not met)
#>   NoGo: P(regions {9} | data) = 0.0835 (threshold 0.6, not met)
#> Region probabilities: R1 = 0.4445, R2 = 0.0745, R3 = 0.227, R4 = 0.0365, R5 = 0.011, R6 = 0.0215, R7 = 0.082, R8 = 0.0195, R9 = 0.0835
#> Posterior distributions
#>   Treatment: Dirichlet(alpha = (3.25, 2.25, 2.25, 3.25))
#>   Control:   Dirichlet(alpha = (6.25, 2.25, 1.25, 1.25))

The posterior probability of region 1 is 0.445, and the decision is Gray. With rule = "and", the Go criterion requires each endpoint to meet its threshold separately:

des_and <- design_binary2(
  n_t = 10, n_c = 10, prior_t = prior, prior_c = prior,
  criterion = criterion_posterior(theta_tv = c(0.15, 0.15),
                                  theta_mav = c(0.10, 0.10),
                                  gamma_go = 0.6, gamma_nogo = 0.6,
                                  rule = "and"),
  nMC = 2000
)
qdm_decide(des_and, data_t = data_binary2(c(3, 2, 2, 3)),
           data_c = data_binary2(c(6, 2, 1, 1)), seed = 1)
#> Decision: Gray
#>   Go:   not met (rule "and")
#>     P(theta1 > 0.15 | data) = 0.746 (threshold 0.6, met)
#>     P(theta2 > 0.15 | data) = 0.563 (threshold 0.6, not met)
#>   NoGo: not met (rule "and")
#>     P(theta1 <= 0.1 | data) = 0.185 (threshold 0.6, not met)
#>     P(theta2 <= 0.1 | data) = 0.332 (threshold 0.6, not met)
#> Region probabilities: R1 = 0.4445, R2 = 0.0745, R3 = 0.227, R4 = 0.0365, R5 = 0.011, R6 = 0.0215, R7 = 0.082, R8 = 0.0195, R9 = 0.0835
#> Posterior distributions
#>   Treatment: Dirichlet(alpha = (3.25, 2.25, 2.25, 3.25))
#>   Control:   Dirichlet(alpha = (6.25, 2.25, 1.25, 1.25))

For the predictive criterion with 30 future subjects per group:

des_pred <- design_binary2(
  n_t = 10, n_c = 10, prior_t = prior, prior_c = prior,
  criterion = criterion_predictive(theta_null = c(0, 0), m_t = 30,
                                   gamma_go = 0.6, gamma_nogo = 0.6)
)
qdm_decide(des_pred, data_t = data_binary2(c(3, 2, 2, 3)),
           data_c = data_binary2(c(6, 2, 1, 1)))
#> Decision: Go
#>   Go:   P(regions {1} | data) = 0.6816 (threshold 0.6, met)
#>   NoGo: P(regions {4} | data) = 0.04534 (threshold 0.6, not met)
#> Region probabilities: R1 = 0.6816, R2 = 0.1923, R3 = 0.08079, R4 = 0.04534
#> Posterior distributions
#>   Treatment: Dirichlet(alpha = (3.25, 2.25, 2.25, 3.25))
#>   Control:   Dirichlet(alpha = (6.25, 2.25, 1.25, 1.25))

External data and hypothetical controls

External data of either group are given as cell counts and borrowed with a power prior weight:

des_ext <- design_binary2(
  n_t = 10, n_c = 10, prior_t = prior, prior_c = prior,
  criterion = criterion_posterior(theta_tv = c(0.15, 0.15),
                                  theta_mav = c(0.10, 0.10),
                                  gamma_go = 0.6, gamma_nogo = 0.6),
  external_c = borrow(data_binary2(c(12, 3, 3, 2)), a0 = 0.5),
  nMC = 2000
)
qdm_decide(des_ext, data_t = data_binary2(c(3, 2, 2, 3)),
           data_c = data_binary2(c(6, 2, 1, 1)), seed = 1)$posterior_c
#>    00    01    10    11 
#> 12.25  3.75  2.75  2.25

For a study without a concurrent control group, hypothetical() with cell counts adds the hypothetical control data to the prior of the control group:

des_hyp <- design_binary2(
  n_t = 10, n_c = 0, prior_t = prior, prior_c = prior,
  criterion = criterion_posterior(theta_tv = c(0.15, 0.15),
                                  theta_mav = c(0.10, 0.10),
                                  gamma_go = 0.6, gamma_nogo = 0.6),
  hypothetical_c = hypothetical(data = data_binary2(c(6, 2, 1, 1))),
  nMC = 2000
)
qdm_decide(des_hyp, data_t = data_binary2(c(3, 2, 2, 3)), seed = 1)
#> Decision: Gray
#>   Go:   P(regions {1} | data) = 0.4445 (threshold 0.6, not met)
#>   NoGo: P(regions {9} | data) = 0.0835 (threshold 0.6, not met)
#> Region probabilities: R1 = 0.4445, R2 = 0.0745, R3 = 0.227, R4 = 0.0365, R5 = 0.011, R6 = 0.0215, R7 = 0.082, R8 = 0.0195, R9 = 0.0835
#> Posterior distributions
#>   Treatment: Dirichlet(alpha = (3.25, 2.25, 2.25, 3.25))
#>   Control:   Dirichlet(alpha = (6.25, 2.25, 1.25, 1.25))

With fixed_c, the control response rates of the two endpoints are known values $(\pi_{c1,0}, \pi_{c2,0})$ , for example those of the standard of care reported in the literature. The effects are then $\theta_j = \pi_{tj} - \pi_{cj,0}$ , and prior_c must be NULL:

des_fix <- design_binary2(
  n_t = 10, n_c = 0, prior_t = prior,
  criterion = criterion_posterior(theta_tv = c(0.15, 0.15),
                                  theta_mav = c(0.10, 0.10),
                                  gamma_go = 0.6, gamma_nogo = 0.6),
  fixed_c = c(0.2, 0.2), nMC = 2000
)
qdm_decide(des_fix, data_t = data_binary2(c(3, 2, 2, 3)), seed = 1)
#> Decision: Go
#>   Go:   P(regions {1} | data) = 0.738 (threshold 0.6, met)
#>   NoGo: P(regions {9} | data) = 0.015 (threshold 0.6, not met)
#> Region probabilities: R1 = 0.738, R2 = 0.0545, R3 = 0.0665, R4 = 0.0505, R5 = 0.0075, R6 = 0.0055, R7 = 0.055, R8 = 0.0075, R9 = 0.015
#> Posterior distributions
#>   Treatment: Dirichlet(alpha = (3.25, 2.25, 2.25, 3.25))
#>   Control:   fixed values (0.2, 0.2)

Efficacy and safety endpoints

The favorable direction can differ between the endpoints. For an efficacy endpoint (a response, on which larger effects are favorable) and a safety endpoint (the occurrence of an adverse event, on which smaller effects are favorable), direction = c("greater", "less") reverses the inequalities of endpoint 2, so that region 1 is still the most favorable one:

qdm_regions("posterior", direction = c("greater", "less"))
#>                      theta1 > TV1 MAV1 < theta1 <= TV1 theta1 <= MAV1
#> theta2 < TV2         R1           R4                   R7            
#> TV2 <= theta2 < MAV2 R2           R5                   R8            
#> theta2 >= MAV2       R3           R6                   R9

Here the first digit of a cell is the response and the second the adverse event. The decision is Go if, with probability at least 0.5, the response rate is improved by more than 0.15 and the rate of the adverse event is increased by less than 0.1 (region 1), and NoGo if, with probability at least 0.6, the response rate is improved by at most 0.1 and the rate of the adverse event is increased by at least 0.2 (region 9):

des_safety <- design_binary2(
  n_t = 10, n_c = 10, prior_t = prior, prior_c = prior,
  criterion = criterion_posterior(theta_tv = c(0.15, 0.1),
                                  theta_mav = c(0.10, 0.2),
                                  gamma_go = 0.5, gamma_nogo = 0.6,
                                  direction = c("greater", "less")),
  nMC = 2000
)
qdm_decide(des_safety, data_t = data_binary2(c(2, 1, 5, 2)),
           data_c = data_binary2(c(6, 1, 2, 1)), seed = 1)
#> Decision: Gray
#>   Go:   P(regions {1} | data) = 0.443 (threshold 0.5, not met)
#>   NoGo: P(regions {9} | data) = 0.0185 (threshold 0.6, not met)
#> Region probabilities: R1 = 0.443, R2 = 0.186, R3 = 0.235, R4 = 0.029, R5 = 0.011, R6 = 0.0085, R7 = 0.0495, R8 = 0.0195, R9 = 0.0185
#> Posterior distributions
#>   Treatment: Dirichlet(alpha = (2.25, 1.25, 5.25, 2.25))
#>   Control:   Dirichlet(alpha = (6.25, 1.25, 2.25, 1.25))

Operating characteristics

A scenario gives the response rates of the two endpoints and their correlation in each group. cellprob_binary2() converts them into the cell probabilities; a correlation is feasible only if all four probabilities are non-negative.

cellprob_binary2(pi1 = 0.35, pi2 = 0.35, rho = 0.3)
#>      00      01      10      11 
#> 0.49075 0.15925 0.15925 0.19075

By default (enumerate = TRUE), qdm_oc() enumerates all outcomes of the study, choose(n + 3, 3) per group (286 for 10 subjects), and computes the region probabilities once for every pair of treatment and control outcomes, with the same Monte Carlo draws for all pairs. The computing time is proportional to the number of pairs times nMC; for large sample sizes, enumerate = FALSE simulates nsim studies instead.

effects <- expand.grid(theta1 = c(0, 0.15, 0.3), theta2 = c(0, 0.15, 0.3))
scen <- data.frame(pi_t1 = 0.2 + effects$theta1,
                   pi_t2 = 0.2 + effects$theta2, rho_t = 0.3,
                   pi_c1 = 0.2, pi_c2 = 0.2, rho_c = 0.3)
oc <- qdm_oc(des, scen, seed = 1)
oc
#> Operating characteristics (exact enumeration of the outcomes; posterior probabilities by Monte Carlo simulation (nMC = 2000))
#>  pi_t1 pi_t2 rho_t pi_c1 pi_c2 rho_c      Go    NoGo   Gray Miss
#>   0.20  0.20   0.3   0.2   0.2   0.3 0.01424 0.39177 0.5940    0
#>   0.35  0.20   0.3   0.2   0.2   0.3 0.03924 0.20096 0.7598    0
#>   0.50  0.20   0.3   0.2   0.2   0.3 0.06158 0.07620 0.8622    0
#>   0.20  0.35   0.3   0.2   0.2   0.3 0.03901 0.20123 0.7598    0
#>   0.35  0.35   0.3   0.2   0.2   0.3 0.12371 0.10564 0.7707    0
#>   0.50  0.35   0.3   0.2   0.2   0.3 0.21431 0.04136 0.7443    0
#>   0.20  0.50   0.3   0.2   0.2   0.3 0.06118 0.07615 0.8627    0
#>   0.35  0.50   0.3   0.2   0.2   0.3 0.21322 0.04123 0.7455    0
#>   0.50  0.50   0.3   0.2   0.2   0.3 0.39239 0.01672 0.5909    0
plot(oc)

A scenario whose correlation is outside the feasible range is not evaluable: its probabilities are NA, with a warning, and it is marked “NE” in the plot.

With overlay = TRUE, the decisions are drawn in a single heat map, as in Figure 2 of Yamaguchi et al. (submitted): each cell has the color of the decision with the highest probability, shaded in five steps of that probability (below 0.5, 0.5 to 0.6, 0.6 to 0.7, 0.7 to 0.8, and 0.8 or more). In the scenarios above, Gray decisions are the most likely in every cell; with effects up to 0.5, Go decisions become the most likely when both effects are large. cell_values = FALSE leaves out the probabilities written in the cells:

eff_wide <- expand.grid(theta1 = seq(0, 0.5, by = 0.1),
                        theta2 = seq(0, 0.5, by = 0.1))
scen_wide <- data.frame(pi_t1 = 0.2 + eff_wide$theta1,
                        pi_t2 = 0.2 + eff_wide$theta2, rho_t = 0.3,
                        pi_c1 = 0.2, pi_c2 = 0.2, rho_c = 0.3)
oc_wide <- qdm_oc(des, scen_wide, seed = 1)
plot(oc_wide, overlay = TRUE)

plot(oc_wide, overlay = TRUE, cell_values = FALSE)

The decision boundaries of qdm_boundary() (see vignette("single-binary")) are not available for two binary endpoints, because the decision depends on the four cell counts of the treatment group and cannot be shown on two axes.

Calibration of the thresholds

With rule = "joint", qdm_calibrate() chooses one threshold for the Go criterion and one for the NoGo criterion. Here the Go criterion should be met with probability at most 0.1 when neither endpoint is improved, and the NoGo criterion with probability at most 0.2 when both effects are 0.3:

cal <- qdm_calibrate(des, scenario_go = scen[1, ], scenario_nogo = scen[9, ],
                     target_go = 0.1, target_nogo = 0.2, seed = 1)
#> Warning: The Go and NoGo criteria are both met (Miss) with positive probability
#> in 2 scenario(s); see the column 'Miss'. A Miss decision is impossible if
#> gamma_go + gamma_nogo > 1.
cal
#> Calibration of the decision thresholds (exact enumeration of the outcomes; posterior probabilities by Monte Carlo simulation (nMC = 2000))
#>   gamma_go   = 0.3: P(Go criterion met | scenario_go) = 0.09792 (target 0.1)
#>   gamma_nogo = 0.17: P(NoGo criterion met | scenario_nogo) = 0.1872 (target 0.2)
#> Decision probabilities with the chosen thresholds
#>       scenario pi_t1 pi_t2 rho_t pi_c1 pi_c2 rho_c     Go   NoGo    Gray
#>    scenario_go   0.2   0.2   0.3   0.2   0.2   0.3 0.0868 0.8291 0.07297
#>  scenario_nogo   0.5   0.5   0.3   0.2   0.2   0.3 0.7192 0.1691 0.09363
#>     Miss
#>  0.01112
#>  0.01810
plot(cal)

The calibrated thresholds add up to 0.47. When the sum is 1 or less, a Go and a NoGo criterion can be met together (a Miss decision); here Miss decisions have probability up to 0.018 under the two scenarios, which the warning reports. qdm_calibrate() chooses the smallest thresholds that meet the targets, and larger thresholds also meet them, with more Gray decisions. For example, the design with $\gamma_{go} = \gamma_{nogo} = 0.6$ used above meets the Go criterion with probability 0.014 when neither endpoint is improved and the NoGo criterion with probability 0.017 when both effects are 0.3, and Miss decisions are impossible because the thresholds add up to more than 1.

Reference

Yamaguchi, Y., Hu, H., Huang, Y., Han, C., Homma, G. and Wang, A. (submitted). Bivariate Bayesian quantitative decision-making for two endpoints in rare disease drug development.