Skip to content

Commit

Permalink
24_predprob.R (#25)
Browse files Browse the repository at this point in the history
* clean

* clean

* clean

* clean lintr errors

* .R: Auto stash before merge of "24_predprob.R" and "main"

* clean

* clean

* clean

* Update R/predprob.R

Co-authored-by: Daniel Sabanes Bove <[email protected]>

* Update R/predprob.R

Co-authored-by: Daniel Sabanes Bove <[email protected]>

* clean

---------

Co-authored-by: Daniel Sabanes Bove <[email protected]>
  • Loading branch information
audreyyeoCH and danielinteractive authored Nov 30, 2023
1 parent 6af999c commit 06ae0a2
Show file tree
Hide file tree
Showing 4 changed files with 242 additions and 129 deletions.
114 changes: 53 additions & 61 deletions R/predprob.R
Original file line number Diff line number Diff line change
Expand Up @@ -2,91 +2,83 @@
#' @include postprob.R
NULL

#' Compute the predictive probability that the trial will be
#' successful, with a fixed response rate threshold
#' Predictive probability of trial success
#'
#' @description `r lifecycle::badge("experimental")`
#'
#' Compute the predictive probability of trial success given current data.
#' Success means that at the end of the trial the posterior probability
#' Pr(P_E > p) >= thetaT,
#' Success means that at the end of the trial the posterior probability is
#' `Pr(P_E > p) >= thetaT`,
#' where p is the fixed response rate threshold.
#' Then the predictive probability for success is:
#' pp = sum over i: Pr(Y=i|x,n)*I{Pr(P_E > p|x,Y=i)>=thetaT}
#' where Y is the number of future responses in the treatment group and x is
#' the current number of responses in the treatment group (out of n).
#' The prior is P_E ~ beta(a, b), default uniform which is a beta(1,1).
#' However, also a beta mixture prior can be specified.
#' `pp = sum over i: Pr(Y=i|x,n)*I{Pr(P_E > p|x,Y=i)>=thetaT}`
#' where `Y` is the number of future responses in the treatment group and `x` is
#' the current number of responses in the treatment group out of n.
#' The prior for the response rate in the experimental arm is `P_E ~ beta(a, b)`.
#'
#' A table with the following contents will be included in the \code{table} attribute of
#' the return value:
#' i: Y=i (number of future successes in Nmax-n subjects)
#' py: Pr(Y=i|x) using beta-(mixture)-binomial distribution
#' b: Pr(P_E > p | x, Y=i) using beta (mixture) posterior
#' bgttheta: indicator I(b>thetaT)
#' A table with the following contents will be included in the return output :
#' - `i`: `Y = i`, number of future successes in `Nmax-n` subjects.
#' - `density`: `Pr(Y = i|x)` using beta-(mixture)-binomial distribution.
#' - `posterior`: `Pr(P_E > p | x, Y=i)` using beta posterior.
#' - `success`: indicator `I(b>thetaT)`.
#'
#' @param x number of successes
#' @param n number of patients
#' @param Nmax maximum number of patients at the end of the trial
#' @param p threshold on the response rate
#' @param thetaT threshold on the probability to be above p
#' @param parE the beta parameters matrix, with K rows and 2 columns,
#' corresponding to the beta parameters of the K components. default is a
#' uniform prior.
#' @param weights the mixture weights of the beta mixture prior. Default are
#' uniform weights across mixture components.
#' @return The predictive probability, a numeric value. In addition, a
#' list called \code{table} is returned as attribute of the returned number.
#' @typed x : number
#' number of successes.
#' @typed n : number
#' number of patients.
#' @typed Nmax : number
#' maximum number of patients at the end of the trial.
#' @typed p : number
#' threshold on the response rate.
#' @typed thetaT : number
#' threshold on the probability to be above p.
#' @typed parE : numeric
#' the beta parameters matrix, with K rows and 2 columns,
#' corresponding to the beta parameters of the K components.
#' @typed weights : numeric
#' the mixture weights of the beta mixture prior.
#' @return A `list` is returned with names `result` for predictive probability and
#' `table` of numeric values with counts of responses in the remaining patients,
#' probabilities of these counts, corresponding probabilities to be above threshold,
#' and trial success indicators.
#'
#' @references Lee, J. J., & Liu, D. D. (2008). A predictive probability
#' design for phase II cancer clinical trials. Clinical Trials, 5(2),
#' 93-106. doi:10.1177/1740774508089279
#' design for phase II cancer clinical trials. Clinical Trials, 5(2),
#' 93-106. doi:10.1177/1740774508089279
#'
#' @example examples/predprob.R
#' @export
predprob <- function(x, n, Nmax, p, thetaT, parE = c(1, 1),
weights) {
## m = Nmax - n future observations
# m = Nmax - n future observations
m <- Nmax - n

## if par is a vector => situation where there is only one component
if (is.vector(parE)) {
## check that it has exactly two entries
stopifnot(identical(length(parE), 2L))

## and transpose to matrix with one row
assert_true(identical(length(parE), 2L))
# and transpose to matrix with one row
parE <- t(parE)
}

## if prior weights of the beta mixture are not supplied
if (missing(weights)) {
weights <- rep(1, nrow(parE))
## (don't need to be normalized, this is done in h_getBetamixPost)
}

## now compute updated parameters
betamixPost <- h_getBetamixPost(x = x, n = n, par = parE, weights = weights)

py <- with(
density <- with(
betamixPost,
dbetabinomMix(x = 0:m, m = m, par = par, weights = weights)
)

## now with the beta binomial mixture:
## posterior probabilities to be above threshold p
b <- postprob(x = x + c(0:m), n = Nmax, p = p, parE = parE, weights = weights)

## prepare the return value
ret <- structure(sum(py * (b > thetaT)),
table =
round(
cbind(
i = c(0:m),
py,
b,
bgttheta = (b > thetaT)
),
4
)
assert_numeric(density, lower = 0, upper = 1, finite = TRUE, any.missing = FALSE)
assert_numeric(posterior, lower = 0, upper = 1, finite = TRUE, any.missing = FALSE)
assert_number(thetaT, lower = 0, upper = 1, finite = TRUE)
# posterior probabilities to be above threshold p
posterior <- postprob(x = x + c(0:m), n = Nmax, p = p, parE = parE, weights = weights)
list(
result = sum(density * (posterior > thetaT)),
table = data.frame(
counts = c(0:m),
cumul_counts = n + (0:m),
density = round(density, 4),
posterior = posterior,
success = (posterior > thetaT)
)
)

return(ret)
}
27 changes: 7 additions & 20 deletions examples/predprob.R
Original file line number Diff line number Diff line change
@@ -1,39 +1,26 @@
## The original Lee and Liu (Table 1) example:
## Nmax=40, x=16, n=23, beta(0.6,0.4) prior distribution,
## thetaT=0.9. The control response rate is 60%:
# The original Lee and Liu (Table 1) example:
# Nmax = 40, x = 16, n = 23, beta(0.6,0.4) prior distribution,
# thetaT = 0.9. The control response rate is 60%:
predprob(
x = 16, n = 23, Nmax = 40, p = 0.6, thetaT = 0.9,
parE = c(0.6, 0.4)
)
## So the predictive probability is 56.6%. 12 responses
## in the remaining 17 patients are required for success.

## Lowering/Increasing the probability threshold thetaT of course increases
## /decreases the predictive probability of success, respectively:
# Lowering/Increasing the probability threshold thetaT of course increases
# /decreases the predictive probability of success, respectively:
predprob(
x = 16, n = 23, Nmax = 40, p = 0.6, thetaT = 0.8,
parE = c(0.6, 0.4)
)
## 70.8%

predprob(
x = 16, n = 23, Nmax = 40, p = 0.6, thetaT = 0.95,
parE = c(0.6, 0.4)
)
## 40.7%

## Instead of a fixed beta prior on the response rate, dynamic borrowing
## from a previous data set is possible by using a beta-mixture prior.
## For example, assume we have a previous trial where we saw 25 responses
## in 40 patients. We would like to use information worth 10 patients in our
## trial, but be robust against deviations from that response rate (maybe
## because it was conducted in slightly different disease or patients).
## Then we can use a beta mixture prior, with the informative component getting
## weight of 1/4:
# Mixed beta prior
predprob(
x = 20, n = 23, Nmax = 40, p = 0.6, thetaT = 0.9,
parE = rbind(c(1, 1), c(25, 15)),
weights = c(3, 1)
)
## Since the response rate in the historical dataset is lower (62.5% < 69.6%)
## than in the trial, but similar, we notice that the predictive probability
## of success is now lower.
83 changes: 35 additions & 48 deletions man/predprob.Rd

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

Loading

0 comments on commit 06ae0a2

Please sign in to comment.