Two continuous endpoints • BayesianQDM

This vignette describes the analysis of a proof-of-concept (PoC) study with two continuous endpoints (Yamaguchi et al., submitted). vignette("two-binary") describes the regions and the criteria for two endpoints, which are the same here.

Model

In group $k$ ( $k = t, c$ ), the responses of a subject on the two endpoints have a bivariate normal distribution $N_2(\mu_k, \Sigma_k)$ , and the data are summarized by the number of subjects, the vector of sample means and the sample covariance matrix (data_continuous2()). With the Normal-Inverse-Wishart prior (prior_niw()), $\mu_k \mid \Sigma_k \sim N_2(\mu_0, \Sigma_k / \kappa_0), \qquad \Sigma_k \sim \textrm{Inv-Wishart}(\nu_0, \Lambda_0),$ the posterior distribution has the same form with parameters $(\mu_n, \kappa_n, \nu_n, \Lambda_n)$ , and the mean has the bivariate t distribution with $\nu_n - 1$ degrees of freedom, location $\mu_n$ and scale matrix $\Lambda_n / \{\kappa_n (\nu_n - 1)\}$ . With the vague prior $p(\mu, \Sigma) \propto |\Sigma|^{-3/2}$ (prior_vague()), the posterior distribution of the mean is $t_{n - 2}(\bar y, S / \{n (n - 2)\})$ , where $S$ is $n - 1$ times the sample covariance matrix.

For the predictive criterion, the mean of $m$ future subjects has the same distribution with the scale matrix $(1 / \kappa_n + 1 / m) \Lambda_n / (\nu_n - 1)$ . Version 0.1.0 of the package and Yamaguchi et al. (submitted) used $(1 + \kappa_n) \Lambda_n / \{\kappa_n (\nu_n - 1) m\}$ , which treats the future subjects as independent draws from the predictive distribution of one subject.

The treatment effects are $\theta_j = \mu_{tj} - \mu_{cj}$ ( $j = 1, 2$ ), and the probabilities of the regions are computed by one of two methods (argument method of design_continuous2()):

"MM" (the default): the difference of the two bivariate t distributions is approximated by a bivariate t distribution with matched moments (Yamaguchi, Homma, Takeda and Maruo, submitted, Theorem 3), whose rectangle probabilities are computed by numerical integration for real degrees of freedom; it requires more than 4 degrees of freedom in both groups;
"MC": Monte Carlo simulation with nMC draws from each distribution.

Decision for observed data

The second case study of Yamaguchi et al. (submitted) has 30 subjects per group, Normal-Inverse-Wishart priors, $\theta_{TV} = (10, 2)$ and $\theta_{MAV} = (5, 1)$ . Here the thresholds are $\gamma_{go} = \gamma_{nogo} = 0.6$ (the article calibrates its own thresholds):

L0 <- matrix(c(80^2, 15, 15, 6^2), 2, 2)
crit <- criterion_posterior(theta_tv = c(10, 2), theta_mav = c(5, 1),
                            gamma_go = 0.6, gamma_nogo = 0.6)
des <- design_continuous2(
  n_t = 30, n_c = 30,
  prior_t = prior_niw(mu0 = c(0, 1), kappa0 = 0.1, nu0 = 4, Lambda0 = L0),
  prior_c = prior_niw(mu0 = c(-10, -1), kappa0 = 0.1, nu0 = 4,
                      Lambda0 = L0),
  criterion = crit
)
data_t <- data_continuous2(30, mean = c(12, 3),
                           cov = matrix(c(6400, 190, 190, 36), 2, 2))
data_c <- data_continuous2(30, mean = c(-10, -1),
                           cov = matrix(c(5600, 150, 150, 42), 2, 2))
dec <- qdm_decide(des, data_t = data_t, data_c = data_c)
dec
#> Decision: Go
#>   Go:   P(regions {1} | data) = 0.6788 (threshold 0.6, met)
#>   NoGo: P(regions {9} | data) = 0.01376 (threshold 0.6, not met)
#> Region probabilities: R1 = 0.6788, R2 = 0.03923, R3 = 0.01266, R4 = 0.06662, R5 = 0.007634, R6 = 0.003074, R7 = 0.1521, R8 = 0.02611, R9 = 0.01376
#> Posterior distributions
#>   Treatment: mean ~ t(df = 33, location = (11.960, 2.993), scale = (193.310, 5.565; 5.565, 1.088))
#>   Control:   mean ~ t(df = 33, location = (-10, -1), scale = (169.939, 4.394; 4.394, 1.262))

The decision is Go. Monte Carlo simulation gives similar probabilities of the regions:

des_mc <- design_continuous2(
  n_t = 30, n_c = 30,
  prior_t = prior_niw(mu0 = c(0, 1), kappa0 = 0.1, nu0 = 4, Lambda0 = L0),
  prior_c = prior_niw(mu0 = c(-10, -1), kappa0 = 0.1, nu0 = 4,
                      Lambda0 = L0),
  criterion = crit, method = "MC", nMC = 10000
)
dec_mc <- qdm_decide(des_mc, data_t = data_t, data_c = data_c, seed = 1)
rbind(MM = dec$prob_regions, MC = dec_mc$prob_regions)
#>           R1         R2         R3         R4          R5          R6        R7
#> MM 0.6787969 0.03922545 0.01265575 0.06661762 0.007633996 0.003073625 0.1521253
#> MC 0.6803000 0.03950000 0.01230000 0.06970000 0.007800000 0.002800000 0.1507000
#>            R8         R9
#> MM 0.02611161 0.01375978
#> MC 0.02520000 0.01170000

Hypothetical control

For a study without a concurrent control group, hypothetical(mean, r) gives the control mean the distribution of the treatment mean with location mean and the scale matrix multiplied by r. The prior of the control group is then not used and is left NULL:

des_hyp <- design_continuous2(
  n_t = 30, n_c = 0,
  prior_t = prior_niw(mu0 = c(0, 1), kappa0 = 0.1, nu0 = 4, Lambda0 = L0),
  criterion = crit,
  hypothetical_c = hypothetical(mean = c(-10, -1), r = 1)
)
dec_hyp <- qdm_decide(des_hyp, data_t = data_t)
dec_hyp
#> Decision: Go
#>   Go:   P(regions {1} | data) = 0.6812 (threshold 0.6, met)
#>   NoGo: P(regions {9} | data) = 0.0132 (threshold 0.6, not met)
#> Region probabilities: R1 = 0.6812, R2 = 0.03398, R3 = 0.009255, R4 = 0.06648, R5 = 0.007249, R6 = 0.00251, R7 = 0.1584, R8 = 0.02774, R9 = 0.0132
#> Posterior distributions
#>   Treatment: mean ~ t(df = 33, location = (11.960, 2.993), scale = (193.310, 5.565; 5.565, 1.088))
#>   Control:   mean ~ t(df = 33, location = (-10, -1), scale = (193.310, 5.565; 5.565, 1.088))

Fixed control

With fixed_c, the control means are known values, for example the means of the standard of care reported in the literature. The probabilities of the regions are then those of the bivariate t distribution of the treatment mean, which are computed by numerical integration for any degrees of freedom, so method is not used:

des_fix <- design_continuous2(
  n_t = 30, n_c = 0,
  prior_t = prior_niw(mu0 = c(0, 1), kappa0 = 0.1, nu0 = 4, Lambda0 = L0),
  criterion = crit, fixed_c = c(-10, -1)
)
dec_fix <- qdm_decide(des_fix, data_t = data_t)
dec_fix
#> Decision: Go
#>   Go:   P(regions {1} | data) = 0.7864 (threshold 0.6, met)
#>   NoGo: P(regions {9} | data) = 0.001807 (threshold 0.6, not met)
#> Region probabilities: R1 = 0.7864, R2 = 0.01441, R3 = 0.001277, R4 = 0.0778, R5 = 0.004079, R6 = 0.0004697, R7 = 0.1035, R8 = 0.0103, R9 = 0.001807
#> Posterior distributions
#>   Treatment: mean ~ t(df = 33, location = (11.960, 2.993), scale = (193.310, 5.565; 5.565, 1.088))
#>   Control:   fixed values (-10, -1)

The control means are not uncertain, so the probability of the Go region is 0.786, compared with 0.681 for the hypothetical control with the same location.

Operating characteristics

qdm_oc() simulates the sample means and the matrix of sums of squares and cross-products of nsim studies under each scenario, with the same random numbers for all scenarios. A scenario gives the means, standard deviations and correlation of the two endpoints in each group. To keep this vignette quick, 1000 studies are simulated per scenario (the default is 10000):

effects <- expand.grid(theta1 = c(0, 10, 20), theta2 = c(0, 2, 4))
scen <- data.frame(mu_t1 = -10 + effects$theta1,
                   mu_t2 = -1 + effects$theta2, sd_t1 = 80, sd_t2 = 6,
                   rho_t = 0.3, mu_c1 = -10, mu_c2 = -1, sd_c1 = 80,
                   sd_c2 = 6, rho_c = 0.3)
oc <- qdm_oc(des, scen, nsim = 1000, seed = 1)
oc
#> Operating characteristics (simulation of 1000 studies; probabilities by moment-matching approximation)
#>  mu_t1 mu_t2 sd_t1 sd_t2 rho_t mu_c1 mu_c2 sd_c1 sd_c2 rho_c    Go  NoGo  Gray
#>    -10    -1    80     6   0.3   -10    -1    80     6   0.3 0.019 0.300 0.681
#>      0    -1    80     6   0.3   -10    -1    80     6   0.3 0.032 0.188 0.780
#>     10    -1    80     6   0.3   -10    -1    80     6   0.3 0.047 0.108 0.845
#>    -10     1    80     6   0.3   -10    -1    80     6   0.3 0.095 0.090 0.815
#>      0     1    80     6   0.3   -10    -1    80     6   0.3 0.153 0.057 0.790
#>     10     1    80     6   0.3   -10    -1    80     6   0.3 0.212 0.034 0.754
#>    -10     3    80     6   0.3   -10    -1    80     6   0.3 0.190 0.012 0.798
#>      0     3    80     6   0.3   -10    -1    80     6   0.3 0.319 0.008 0.673
#>     10     3    80     6   0.3   -10    -1    80     6   0.3 0.468 0.005 0.527
#>  Miss
#>     0
#>     0
#>     0
#>     0
#>     0
#>     0
#>     0
#>     0
#>     0
plot(oc)

The simulation standard error of the probabilities is at most $0.5 / \sqrt{nsim}$ , 0.016 for 1000 studies. The moment-matching approximation is much faster than Monte Carlo estimation of the region probabilities, whose computing time is proportional to nMC for every simulated study.

overlay = TRUE draws a single heat map of the decision with the highest probability in each cell (see vignette("two-binary")), here over a wider range of effects:

eff_wide <- expand.grid(theta1 = seq(0, 40, by = 10),
                        theta2 = seq(0, 8, by = 2))
scen_wide <- data.frame(mu_t1 = -10 + eff_wide$theta1,
                        mu_t2 = -1 + eff_wide$theta2, sd_t1 = 80,
                        sd_t2 = 6, rho_t = 0.3, mu_c1 = -10, mu_c2 = -1,
                        sd_c1 = 80, sd_c2 = 6, rho_c = 0.3)
plot(qdm_oc(des, scen_wide, nsim = 1000, seed = 1), overlay = TRUE)

Decision boundary

For each observed treatment mean of endpoint 1, qdm_boundary() gives the intervals of the observed treatment mean of endpoint 2 that lead to Go, NoGo and Gray decisions, for a given sample covariance matrix, together with the observed differences in means. The ends of the intervals are located by bisection between the values of the grid mean_t[[2]]. With a concurrent control group, the observed control data are given by data_c:

qdm_boundary(des, mean_t = list(c(0, 10, 20, 30), seq(-3, 7, by = 0.25)),
             cov_t = data_t$cov, data_c = data_c)
#> Decision boundary of a design with two continuous endpoints
#>   n_t = 30
#>   Control data: two continuous endpoints, n = 30, mean = (-10, -1), cov = (5600, 150; 150, 42)
#> 
#> Observed treatment mean of endpoint 2
#>   effect1   Gray       Go
#>   10        all        -
#>   20        <= 2.39    >= 2.39
#>   30        <= 1.696   >= 1.696
#>   40        <= 1.489   >= 1.489
#> 
#> Observed difference in means of endpoint 2 (control mean -1)
#>   effect1   Gray       Go
#>   10        all        -
#>   20        <= 3.39    >= 3.39
#>   30        <= 2.696   >= 2.696
#>   40        <= 2.489   >= 2.489

With more values of endpoint 1, the plot shows the regions of the decisions over the observed differences in means of the two endpoints:

b <- qdm_boundary(des, mean_t = list(seq(-30, 50, by = 2),
                                     seq(-3, 7, by = 0.25)),
                  cov_t = data_t$cov, data_c = data_c)
plot(b)

Calibration of the thresholds

The thresholds are calibrated so that the Go criterion is met with probability at most 0.1 when neither endpoint is improved, and the NoGo criterion with probability at most 0.2 when the effects are 20 and 4:

cal <- qdm_calibrate(des, scenario_go = scen[1, ], scenario_nogo = scen[9, ],
                     target_go = 0.1, target_nogo = 0.2, nsim = 1000,
                     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 (simulation of 1000 studies; probabilities by moment-matching approximation)
#>   gamma_go   = 0.3: P(Go criterion met | scenario_go) = 0.1 (target 0.1)
#>   gamma_nogo = 0.06: P(NoGo criterion met | scenario_nogo) = 0.195 (target 0.2)
#> Decision probabilities with the chosen thresholds
#>       scenario mu_t1 mu_t2 sd_t1 sd_t2 rho_t mu_c1 mu_c2 sd_c1 sd_c2 rho_c
#>    scenario_go   -10    -1    80     6   0.3   -10    -1    80     6   0.3
#>  scenario_nogo    10     3    80     6   0.3   -10    -1    80     6   0.3
#>     Go  NoGo  Gray  Miss
#>  0.065 0.861 0.039 0.035
#>  0.711 0.142 0.094 0.053
plot(cal)

The calibrated thresholds add up to 0.36. 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.053 under the two scenarios, which the warning reports. Larger thresholds also meet the targets, with more Gray decisions: the design with $\gamma_{go} = \gamma_{nogo} = 0.6$ used above meets the Go criterion with probability 0.019 when neither endpoint is improved and the NoGo criterion with probability 0.005 when the effects are 20 and 4, and Miss decisions are impossible because the thresholds add up to more than 1.

References

Yamaguchi, Y., Homma, G., Takeda, K. and Maruo, K. (submitted). Moment-matching approximation for t-differences in clinical trials with continuous endpoints.

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.