Skip to contents

Computes the posterior or posterior predictive probability that the bivariate treatment effect \(\theta = \mu_t - \mu_c\) falls in each of the 9 (posterior) or 4 (predictive) rectangular regions defined by the decision thresholds.

Usage

pbayespostpred2cont(
  prob,
  design,
  prior,
  theta_TV1 = NULL,
  theta_MAV1 = NULL,
  theta_TV2 = NULL,
  theta_MAV2 = NULL,
  theta_NULL1 = NULL,
  theta_NULL2 = NULL,
  n_t,
  n_c = NULL,
  ybar_t,
  S_t,
  ybar_c = NULL,
  S_c = NULL,
  m_t = NULL,
  m_c = NULL,
  kappa0_t = NULL,
  nu0_t = NULL,
  mu0_t = NULL,
  Lambda0_t = NULL,
  kappa0_c = NULL,
  nu0_c = NULL,
  mu0_c = NULL,
  Lambda0_c = NULL,
  r = NULL,
  ne_t = NULL,
  ne_c = NULL,
  alpha0e_t = NULL,
  alpha0e_c = NULL,
  bar_ye_t = NULL,
  bar_ye_c = NULL,
  se_t = NULL,
  se_c = NULL,
  nMC = 10000L,
  CalcMethod = "MC"
)

Arguments

prob

A character string specifying the probability type. Must be 'posterior' or 'predictive'.

design

A character string specifying the trial design. Must be 'controlled', 'uncontrolled', or 'external'.

prior

A character string specifying the prior distribution. Must be 'vague' or 'N-Inv-Wishart'.

theta_TV1

A numeric scalar giving the target value (TV) threshold for Endpoint 1. Required when prob = 'posterior'; must satisfy theta_TV1 > theta_MAV1. Set to NULL when prob = 'predictive'.

theta_MAV1

A numeric scalar giving the minimum acceptable value (MAV) threshold for Endpoint 1. Required when prob = 'posterior'; must satisfy theta_TV1 > theta_MAV1. Set to NULL when prob = 'predictive'.

theta_TV2

A numeric scalar giving the target value (TV) threshold for Endpoint 2. Required when prob = 'posterior'; must satisfy theta_TV2 > theta_MAV2. Set to NULL when prob = 'predictive'.

theta_MAV2

A numeric scalar giving the minimum acceptable value (MAV) threshold for Endpoint 2. Required when prob = 'posterior'; must satisfy theta_TV2 > theta_MAV2. Set to NULL when prob = 'predictive'.

theta_NULL1

A numeric scalar giving the null hypothesis threshold for Endpoint 1. Required when prob = 'predictive'; set to NULL when prob = 'posterior'.

theta_NULL2

A numeric scalar giving the null hypothesis threshold for Endpoint 2. Required when prob = 'predictive'; set to NULL when prob = 'posterior'.

n_t

A positive integer giving the number of patients in the treatment group in the proof-of-concept (PoC) trial.

n_c

A positive integer giving the number of patients in the control group in the PoC trial. For design = 'uncontrolled', this is the hypothetical control sample size (required for consistency with other designs).

ybar_t

A numeric vector of length 2 or a numeric matrix with 2 columns giving the sample mean(s) for the treatment group. When a matrix with \(N\) rows is supplied, the function computes probabilities for all \(N\) observations simultaneously and returns an \(N \times n_{\rm regions}\) matrix.

S_t

A 2x2 numeric matrix giving the sum-of-squares matrix for the treatment group (single observation), or a list of \(N\) such matrices (vectorised call). Must be consistent with ybar_t: if ybar_t is a matrix, S_t must be a list of the same length.

ybar_c

A numeric vector of length 2, a numeric matrix with 2 columns, or NULL giving the sample mean(s) for the control group. Not used when design = 'uncontrolled'.

S_c

A 2x2 numeric matrix, a list of 2x2 matrices, or NULL giving the sum-of-squares matrix/matrices for the control group. Not used when design = 'uncontrolled'.

m_t

A positive integer giving the number of patients in the treatment group for the future trial. Required when prob = 'predictive'; otherwise set to NULL.

m_c

A positive integer giving the number of patients in the control group for the future trial. Required when prob = 'predictive'; otherwise set to NULL.

kappa0_t

A positive numeric scalar giving the NIW prior concentration parameter for the treatment group. Required when prior = 'N-Inv-Wishart'; otherwise set to NULL.

nu0_t

A numeric scalar giving the NIW prior degrees of freedom for the treatment group. Must be greater than 3. Required when prior = 'N-Inv-Wishart'; otherwise set to NULL.

mu0_t

A numeric vector of length 2 giving the NIW prior mean for the treatment group. Required when prior = 'N-Inv-Wishart'; otherwise set to NULL.

Lambda0_t

A 2x2 positive-definite numeric matrix giving the NIW prior scale matrix for the treatment group. Required when prior = 'N-Inv-Wishart'; otherwise set to NULL.

kappa0_c

A positive numeric scalar giving the NIW prior concentration parameter for the control group. Required when prior = 'N-Inv-Wishart' and design != 'uncontrolled'; otherwise set to NULL.

nu0_c

A numeric scalar giving the NIW prior degrees of freedom for the control group. Must be greater than 3. Required when prior = 'N-Inv-Wishart' and design != 'uncontrolled'; otherwise set to NULL.

mu0_c

A numeric vector of length 2 giving the NIW prior mean for the control group, or the hypothetical control location when design = 'uncontrolled'. Required when prior = 'N-Inv-Wishart'; otherwise set to NULL.

Lambda0_c

A 2x2 positive-definite numeric matrix giving the NIW prior scale matrix for the control group. Required when prior = 'N-Inv-Wishart' and design != 'uncontrolled'; otherwise set to NULL.

r

A positive numeric scalar giving the variance scaling factor for the hypothetical control distribution. Required when design = 'uncontrolled'; otherwise set to NULL.

ne_t

A positive integer giving the external treatment group sample size. Required when design = 'external' and external treatment data are used; otherwise set to NULL.

ne_c

A positive integer giving the external control group sample size. Required when design = 'external' and external control data are used; otherwise set to NULL.

alpha0e_t

A numeric scalar in (0, 1] giving the power prior weight for the external treatment data. Required when external treatment data are used; otherwise set to NULL.

alpha0e_c

A numeric scalar in (0, 1] giving the power prior weight for the external control data. Required when external control data are used; otherwise set to NULL.

bar_ye_t

A numeric vector of length 2 giving the external treatment group sample mean. Required when external treatment data are used; otherwise set to NULL.

bar_ye_c

A numeric vector of length 2 giving the external control group sample mean. Required when external control data are used; otherwise set to NULL.

se_t

A 2x2 numeric matrix giving the external treatment group sum-of-squares matrix. Required when external treatment data are used; otherwise set to NULL.

se_c

A 2x2 numeric matrix giving the external control group sum-of-squares matrix. Required when external control data are used; otherwise set to NULL.

nMC

A positive integer giving the number of Monte Carlo draws used to estimate region probabilities. Default is 10000. Required when CalcMethod = 'MC'. May be set to NULL when CalcMethod = 'MM' and \(\nu_k > 4\) (the MM method uses mvtnorm::pmvt analytically); if CalcMethod = 'MM' but \(\nu_k \le 4\) causes a fallback to MC, nMC must be a positive integer.

CalcMethod

A character string specifying the computation method. Must be 'MC' (Monte Carlo, default) or 'MM' (Moment-Matching via mvtnorm::pmvt). When CalcMethod = 'MM' and \(\nu_k \le 4\), a warning is issued and the function falls back to CalcMethod = 'MC'.

Value

When ybar_t is a length-2 vector (single observation): a named numeric vector of length 9 (R1R9) for prob = 'posterior' or length 4 (R1R4) for prob = 'predictive'. All elements are non-negative and sum to 1.

    When \code{ybar_t} is an \eqn{N \times 2} matrix (vectorised call):
    a numeric matrix with \eqn{N} rows and 9 (or 4) columns, with
    column names \code{R1}--\code{R9} (or \code{R1}--\code{R4}). Each
    row sums to 1.

Details

Model. Both endpoints follow a bivariate Normal distribution \(y_{k,j} \sim N_2(\mu_k, \Sigma_k)\) for group \(k \in \{t, c\}\). The treatment effect is \(\theta = \mu_t - \mu_c\).

Posterior distribution (vague prior). $$\mu_k | Y_k \sim t_{n_k - 2}\!\left(\bar{y}_k,\; \frac{S_k}{n_k(n_k - 2)}\right)$$

Posterior distribution (NIW prior). $$\mu_k | Y_k \sim t_{\nu_{nk} - 1}\!\left(\mu_{nk},\; \frac{\Lambda_{nk}}{\kappa_{nk}(\nu_{nk} - 1)}\right)$$ with updated hyperparameters \(\kappa_{nk} = \kappa_{0k} + n_k\), \(\nu_{nk} = \nu_{0k} + n_k\), \(\mu_{nk} = (\kappa_{0k}\mu_{0k} + n_k\bar{y}_k)/\kappa_{nk}\), and \(\Lambda_{nk} = \Lambda_{0k} + S_k + \kappa_{0k}n_k(\bar{y}_k - \mu_{0k})(\bar{y}_k - \mu_{0k})^T / \kappa_{nk}\).

Predictive distribution. The scale matrix of a single future observation is inflated by \((1 + n_k)/n_k\) (vague) or \((1 + \kappa_{nk})/\kappa_{nk}\) (NIW) relative to the posterior. The mean of \(m_k\) future observations has scale divided by \(m_k\).

Posterior probability regions (prob = 'posterior'). Row-major 3x3 grid; Endpoint 1 varies slowest:

  • R1: \(\theta_1 > TV_1\) AND \(\theta_2 > TV_2\)

  • R2: \(\theta_1 > TV_1\) AND \(TV_2 \ge \theta_2 > MAV_2\)

  • R3: \(\theta_1 > TV_1\) AND \(\theta_2 \le MAV_2\)

  • R4: \(TV_1 \ge \theta_1 > MAV_1\) AND \(\theta_2 > TV_2\)

  • R5: \(TV_1 \ge \theta_1 > MAV_1\) AND \(TV_2 \ge \theta_2 > MAV_2\)

  • R6: \(TV_1 \ge \theta_1 > MAV_1\) AND \(\theta_2 \le MAV_2\)

  • R7: \(\theta_1 \le MAV_1\) AND \(\theta_2 > TV_2\)

  • R8: \(\theta_1 \le MAV_1\) AND \(TV_2 \ge \theta_2 > MAV_2\)

  • R9: \(\theta_1 \le MAV_1\) AND \(\theta_2 \le MAV_2\)

Predictive probability regions (prob = 'predictive'). Row-major 2x2 grid; Endpoint 1 varies slowest:

  • R1: \(\tilde\theta_1 > \theta_{\rm NULL1}\) AND \(\tilde\theta_2 > \theta_{\rm NULL2}\)

  • R2: \(\tilde\theta_1 > \theta_{\rm NULL1}\) AND \(\tilde\theta_2 \le \theta_{\rm NULL2}\)

  • R3: \(\tilde\theta_1 \le \theta_{\rm NULL1}\) AND \(\tilde\theta_2 > \theta_{\rm NULL2}\)

  • R4: \(\tilde\theta_1 \le \theta_{\rm NULL1}\) AND \(\tilde\theta_2 \le \theta_{\rm NULL2}\)

Uncontrolled design. Under NIW prior: \(\mu_c \sim t_{\nu_{nt} - 1}(\mu_{0c},\; r\Lambda_{nt} / [\kappa_{nt}(\nu_{nt} - 1)])\). The parameter r allows the variance scale of the hypothetical control to differ from the treatment group.

External design. Incorporated via the power prior with NIW conjugate representation. The effective posterior hyperparameters are obtained by constructing the power prior from external data with weight \(a_0\), then updating with current PoC data (see Conjugate_power_prior.pdf, Theorem 5).

Vectorised usage. When ybar_t is supplied as an \(N \times 2\) matrix and S_t as a list of \(N\) scatter matrices, the function computes region probabilities for all \(N\) observations in a single call, returning an \(N \times n_{\rm regions}\) matrix. For CalcMethod = 'MC', standard normal and chi-squared variates are pre-generated once (size nMC) and reused across all observations, with only the Cholesky factor of the replicate-specific scale matrix recomputed per observation. For CalcMethod = 'MM', the moment-matching parameters are computed per observation (since they depend on the replicate-specific \(V_k\)) and mvtnorm::pmvt is called once per region per observation.

Examples

# Example 1: Controlled design - posterior probability, vague prior
S_t <- matrix(c(18.0, 3.6, 3.6, 9.0), 2, 2)
S_c <- matrix(c(16.0, 2.8, 2.8, 8.5), 2, 2)
pbayespostpred2cont(
  prob = 'posterior', design = 'controlled', prior = 'vague',
  theta_TV1 = 1.5, theta_MAV1 = 0.5,
  theta_TV2 = 1.0, theta_MAV2 = 0.3,
  theta_NULL1 = NULL, theta_NULL2 = NULL,
  n_t = 12L, n_c = 12L,
  ybar_t = c(3.5, 2.1), S_t = S_t,
  ybar_c = c(1.8, 1.0), S_c = S_c,
  m_t = NULL, m_c = NULL,
  kappa0_t = NULL, nu0_t = NULL, mu0_t = NULL, Lambda0_t = NULL,
  kappa0_c = NULL, nu0_c = NULL, mu0_c = NULL, Lambda0_c = NULL,
  r = NULL,
  ne_t = NULL, ne_c = NULL, alpha0e_t = NULL, alpha0e_c = NULL,
  bar_ye_t = NULL, bar_ye_c = NULL, se_t = NULL, se_c = NULL,
  nMC = 1000L
)
#>    R1    R2    R3    R4    R5    R6    R7    R8    R9 
#> 0.421 0.211 0.019 0.139 0.169 0.013 0.007 0.020 0.001 

# Example 2: Uncontrolled design - posterior probability, NIW prior
S_t <- matrix(c(18.0, 3.6, 3.6, 9.0), 2, 2)
L0  <- matrix(c(20.0, 0.0, 0.0, 10.0), 2, 2)
pbayespostpred2cont(
  prob = 'posterior', design = 'uncontrolled', prior = 'N-Inv-Wishart',
  theta_TV1 = 1.5, theta_MAV1 = 0.5,
  theta_TV2 = 1.0, theta_MAV2 = 0.3,
  theta_NULL1 = NULL, theta_NULL2 = NULL,
  n_t = 12L, n_c = NULL,
  ybar_t = c(3.5, 2.1), S_t = S_t,
  ybar_c = NULL, S_c = NULL,
  m_t = NULL, m_c = NULL,
  kappa0_t = 2.0, nu0_t = 5.0, mu0_t = c(2.0, 1.0), Lambda0_t = L0,
  kappa0_c = NULL, nu0_c = NULL, mu0_c = c(1.0, 0.5), Lambda0_c = NULL,
  r = 1.0,
  ne_t = NULL, ne_c = NULL, alpha0e_t = NULL, alpha0e_c = NULL,
  bar_ye_t = NULL, bar_ye_c = NULL, se_t = NULL, se_c = NULL,
  nMC = 1000L
)
#>    R1    R2    R3    R4    R5    R6    R7    R8    R9 
#> 0.750 0.134 0.007 0.076 0.024 0.004 0.003 0.002 0.000 

# Example 3: External design - posterior probability, NIW prior
S_t  <- matrix(c(18.0, 3.6, 3.6, 9.0), 2, 2)
S_c  <- matrix(c(16.0, 2.8, 2.8, 8.5), 2, 2)
L0   <- matrix(c(20.0, 0.0, 0.0, 10.0), 2, 2)
se_c <- matrix(c(15.0, 2.5, 2.5, 7.5), 2, 2)
pbayespostpred2cont(
  prob = 'posterior', design = 'external', prior = 'N-Inv-Wishart',
  theta_TV1 = 1.5, theta_MAV1 = 0.5,
  theta_TV2 = 1.0, theta_MAV2 = 0.3,
  theta_NULL1 = NULL, theta_NULL2 = NULL,
  n_t = 12L, n_c = 12L,
  ybar_t = c(3.5, 2.1), S_t = S_t,
  ybar_c = c(1.8, 1.0), S_c = S_c,
  m_t = NULL, m_c = NULL,
  kappa0_t = 2.0, nu0_t = 5.0, mu0_t = c(2.0, 1.0), Lambda0_t = L0,
  kappa0_c = 2.0, nu0_c = 5.0, mu0_c = c(1.0, 0.5), Lambda0_c = L0,
  r = NULL,
  ne_t = NULL, ne_c = 10L, alpha0e_t = NULL, alpha0e_c = 0.5,
  bar_ye_t = NULL, bar_ye_c = c(1.5, 0.8), se_t = NULL, se_c = se_c,
  nMC = 1000L
)
#>    R1    R2    R3    R4    R5    R6    R7    R8    R9 
#> 0.352 0.246 0.013 0.178 0.160 0.020 0.012 0.016 0.003 

# Example 4: Controlled design - posterior predictive probability, vague prior
S_t <- matrix(c(18.0, 3.6, 3.6, 9.0), 2, 2)
S_c <- matrix(c(16.0, 2.8, 2.8, 8.5), 2, 2)
pbayespostpred2cont(
  prob = 'predictive', design = 'controlled', prior = 'vague',
  theta_TV1 = NULL, theta_MAV1 = NULL,
  theta_TV2 = NULL, theta_MAV2 = NULL,
  theta_NULL1 = 0.5, theta_NULL2 = 0.3,
  n_t = 12L, n_c = 12L,
  ybar_t = c(3.5, 2.1), S_t = S_t,
  ybar_c = c(1.8, 1.0), S_c = S_c,
  m_t = 30L, m_c = 30L,
  kappa0_t = NULL, nu0_t = NULL, mu0_t = NULL, Lambda0_t = NULL,
  kappa0_c = NULL, nu0_c = NULL, mu0_c = NULL, Lambda0_c = NULL,
  r = NULL,
  ne_t = NULL, ne_c = NULL, alpha0e_t = NULL, alpha0e_c = NULL,
  bar_ye_t = NULL, bar_ye_c = NULL, se_t = NULL, se_c = NULL,
  nMC = 1000L
)
#>    R1    R2    R3    R4 
#> 0.995 0.002 0.003 0.000 

# Example 5: Uncontrolled design - posterior predictive probability, NIW prior
S_t <- matrix(c(18.0, 3.6, 3.6, 9.0), 2, 2)
L0  <- matrix(c(20.0, 0.0, 0.0, 10.0), 2, 2)
pbayespostpred2cont(
  prob = 'predictive', design = 'uncontrolled', prior = 'N-Inv-Wishart',
  theta_TV1 = NULL, theta_MAV1 = NULL,
  theta_TV2 = NULL, theta_MAV2 = NULL,
  theta_NULL1 = 0.5, theta_NULL2 = 0.3,
  n_t = 12L, n_c = NULL,
  ybar_t = c(3.5, 2.1), S_t = S_t,
  ybar_c = NULL, S_c = NULL,
  m_t = 30L, m_c = 30L,
  kappa0_t = 2.0, nu0_t = 5.0, mu0_t = c(2.0, 1.0), Lambda0_t = L0,
  kappa0_c = NULL, nu0_c = NULL, mu0_c = c(1.0, 0.5), Lambda0_c = NULL,
  r = 1.0,
  ne_t = NULL, ne_c = NULL, alpha0e_t = NULL, alpha0e_c = NULL,
  bar_ye_t = NULL, bar_ye_c = NULL, se_t = NULL, se_c = NULL,
  nMC = 1000L
)
#> R1 R2 R3 R4 
#>  1  0  0  0 

# Example 6: External design - posterior predictive probability, NIW prior
S_t  <- matrix(c(18.0, 3.6, 3.6, 9.0), 2, 2)
S_c  <- matrix(c(16.0, 2.8, 2.8, 8.5), 2, 2)
L0   <- matrix(c(20.0, 0.0, 0.0, 10.0), 2, 2)
se_c <- matrix(c(15.0, 2.5, 2.5, 7.5), 2, 2)
pbayespostpred2cont(
  prob = 'predictive', design = 'external', prior = 'N-Inv-Wishart',
  theta_TV1 = NULL, theta_MAV1 = NULL,
  theta_TV2 = NULL, theta_MAV2 = NULL,
  theta_NULL1 = 0.5, theta_NULL2 = 0.3,
  n_t = 12L, n_c = 12L,
  ybar_t = c(3.5, 2.1), S_t = S_t,
  ybar_c = c(1.8, 1.0), S_c = S_c,
  m_t = 30L, m_c = 30L,
  kappa0_t = 2.0, nu0_t = 5.0, mu0_t = c(2.0, 1.0), Lambda0_t = L0,
  kappa0_c = 2.0, nu0_c = 5.0, mu0_c = c(1.0, 0.5), Lambda0_c = L0,
  r = NULL,
  ne_t = NULL, ne_c = 10L, alpha0e_t = NULL, alpha0e_c = 0.5,
  bar_ye_t = NULL, bar_ye_c = c(1.5, 0.8), se_t = NULL, se_c = se_c,
  nMC = 1000L
)
#>    R1    R2    R3    R4 
#> 0.992 0.004 0.004 0.000