From 8ffd96bd1ebd37c0731c78f9ea8ae3b763c9c64e Mon Sep 17 00:00:00 2001 From: Audrey Yeo Date: Thu, 30 Nov 2023 17:32:28 +0100 Subject: [PATCH 01/15] clean --- R/predprobDist.R | 4 +++- man/predprobDist.Rd | 7 ++++--- 2 files changed, 7 insertions(+), 4 deletions(-) diff --git a/R/predprobDist.R b/R/predprobDist.R index 270807a6..637cf2a3 100644 --- a/R/predprobDist.R +++ b/R/predprobDist.R @@ -5,9 +5,11 @@ NULL #' Compute the predictive probability that the trial will be #' successful, with a prior distribution on the SOC #' +#' @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_S + delta_0 | data) >= thetaT. Then the +#' `Pr(P_E > P_S + delta_0 | data) >= thetaT`. Then the #' predictive probability for success is: #' pp = sum over i: Pr(Y=i|x,n)*I{Pr(P_E > P_S + delta|x,Y=i)>=thetaT}, #' where Y is the number of future responses in the treatment group and x is diff --git a/man/predprobDist.Rd b/man/predprobDist.Rd index 9dffe3f2..31e0e77f 100644 --- a/man/predprobDist.Rd +++ b/man/predprobDist.Rd @@ -59,9 +59,11 @@ The predictive probability, a numeric value. In addition, a list called \code{tables} is returned as attribute of the returned number. } \description{ +\ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#experimental}{\figure{lifecycle-experimental.svg}{options: alt='[Experimental]'}}}{\strong{[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_S + delta_0 | data) >= thetaT. Then the +\code{Pr(P_E > P_S + delta_0 | data) >= thetaT}. Then the predictive probability for success is: pp = sum over i: Pr(Y=i|x,n)*I{Pr(P_E > P_S + delta|x,Y=i)>=thetaT}, where Y is the number of future responses in the treatment group and x is @@ -70,8 +72,7 @@ Prior is P_E ~ beta(a, b), default uniform which is a beta(1,1). However, also a beta mixture prior can be specified. Analogously for P_S either a classic beta prior or a beta mixture prior can be specified. -} -\details{ + Also data on the SOC might be available. Then the predictive probability is more generally defined as pp = sum over i, j: Pr(Y=i|x,n)*Pr(Z=j|xS, nS)*I{Pr(P_E > P_S + From 4518ac228542611b51869f00f6ccb9ca1eab6d3f Mon Sep 17 00:00:00 2001 From: Audrey Yeo Date: Fri, 1 Dec 2023 17:17:14 +0100 Subject: [PATCH 02/15] clean --- R/predprobDist.R | 179 ++++++++++++++++++++++++-------------------- inst/WORDLIST | 1 + man/predprobDist.Rd | 73 ++++++++---------- 3 files changed, 129 insertions(+), 124 deletions(-) diff --git a/R/predprobDist.R b/R/predprobDist.R index 637cf2a3..67a11fdd 100644 --- a/R/predprobDist.R +++ b/R/predprobDist.R @@ -11,64 +11,84 @@ NULL #' Success means that at the end of the trial the posterior probability #' `Pr(P_E > P_S + delta_0 | data) >= thetaT`. Then the #' predictive probability for success is: -#' pp = sum over i: Pr(Y=i|x,n)*I{Pr(P_E > P_S + delta|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). -#' Prior is P_E ~ beta(a, b), default uniform which is a beta(1,1). +#' `pp = sum over i: Pr(Y = i | x , n)*I{Pr(P_E > P_S + delta | 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`. +#' Prior is `P_E ~ beta(a, b)` and uniform which is a beta(1,1). #' However, also a beta mixture prior can be specified. Analogously -#' for P_S either a classic beta prior or a beta mixture prior can be +#' for `P_S` either a classic beta prior or a beta mixture prior can be #' specified. #' #' Also data on the SOC might be available. Then the predictive probability is #' more generally defined as -#' pp = sum over i, j: Pr(Y=i|x,n)*Pr(Z=j|xS, nS)*I{Pr(P_E > P_S + -#' delta|x,xS,Y=i,Z=j)>=thetaT} -#' where Z is the future number of responses in the SOC group, and xS is the +#' `pp = sum over i, j: Pr(Y = i | x, n)*Pr(Z = j | xS, nS )*I{Pr(P_E > P_S + delta | x,xS, Y=i, Z=j ) >= thetaT}` +#' where `Z` is the future number of responses in the SOC group, and `xS` is the #' current number of responses in the SOC group. #' +#' In the case where `NmaxControl = 0`, 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_S + delta | x, Y = i)` using beta posterior. +#' - `success`: indicator `I(b>thetaT)`. +#' #' A table with the following contents will be -#' included in the \code{tables} attribute of the return value +#' included in the return value #' (in the case that NmaxControl is zero): #' 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_S + delta | x, Y=i) #' bgttheta: indicator I(b>thetaT) #' -#' If NmaxControl is not zero, i.e., when data on the control treatment +#' If `NmaxControl` is not zero, i.e., when data on the control treatment #' is available in this trial, then a list with will be included with the #' following elements: #' pyz: matrix with the probabilities Pr(Y=i, Z=j | x, xS) #' b: matrix with Pr(P_E > P_S + delta | x, xS, Y=i, Z=j) #' bgttheta: matrix of indicators I(b>thetaT) #' -#' @param x number of successes (in the treatment group) -#' @param n number of patients (in the treatment group) -#' @param xS number of successes in the SOC group (default: 0) -#' @param nS number of patients in the SOC group (default: 0) -##' @param Nmax maximum number of patients at the end of the trial (in the -##' treatment group) -##' @param NmaxControl maximum number of patients at the end of the trial in the -##' SOC group (default: 0) -##' @param delta margin by which the response rate in the treatment group should -##' be better than in the SOC group (default: 0) -##' @param relativeDelta see \code{\link{postprobDist}} -##' @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. -##' @param parS beta prior parameters in the SOC group (default: uniform) -##' @param weightsS weights for the SOC group (default: uniform) -##' @param thetaT threshold on the probability to be used -##' @return The predictive probability, a numeric value. In addition, a -##' list called \code{tables} is returned as attribute of the returned number. -##' -##' @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 -##' -##' @example examples/predprobDist.R -##' @export +#' @note +#' +#' # Delta from postprobDist : +#' +#' The desired improvement is denoted as `delta`. There are two options in using `delta`. +#' The absolute case when `relativeDelta = FALSE` and relative as when `relativeDelta = TRUE`. +#' +#' 1. The absolute case is when we define an absolute delta, greater than `P_S`, +#' the response rate of the `SOC` group such that +#' the posterior is `Pr(P_E > P_S + delta | data)`. +#' +#' 2. In the relative case, we suppose that the treatment group's +#' response rate is assumed to be greater than `P_S + (1-P_S) * delta` such that +#' the posterior is `Pr(P_E > P_S + (1 - P_S) * delta | data)`. +#' +#' @typed x number of successes (in the treatment group) +#' @typed n number of patients (in the treatment group) +#' @typed xS number of successes in the SOC group (default: 0) +#' @typed nS number of patients in the SOC group (default: 0) +#' @typed Nmax maximum number of patients at the end of the trial (in the +#' treatment group) +#' @typed NmaxControl maximum number of patients at the end of the trial in the +#' SOC group (default: 0) +#' @typed delta margin by which the response rate in the treatment group should +#' be better than in the SOC group (default: 0) +#' @typed relativeDelta see \code{\link{postprobDist}} +#' @typed 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. +#' @typed weights the mixture weights of the beta mixture prior. Default are +#' uniform weights across mixture components. +#' @typed parS beta prior parameters in the SOC group (default: uniform) +#' @typed weightsS weights for the SOC group (default: uniform) +#' @typed thetaT threshold on the probability to be used +#' @return The predictive probability, a numeric value. In addition, a +#' list called \code{tables} is returned as attribute of the returned number. +#' +#' @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 +#' +#' @example examples/predprobDist.R +#' @export predprobDist <- function(x, n, xS = 0, nS = 0, Nmax, NmaxControl = 0, @@ -79,7 +99,7 @@ predprobDist <- function(x, n, parS = c(a = 1, b = 1), weightsS, thetaT) { - ## ensure reasonable numbers + # ensure reasonable numbers stopifnot( n <= Nmax, nS <= NmaxControl, @@ -87,52 +107,49 @@ predprobDist <- function(x, n, xS <= nS ) - ## remaining active patients to be seen: + # remaining active patients to be seen: mE <- Nmax - n - ## if par is a vector => situation where there is only one component + # if par is a vector => situation where there is only one component if (is.vector(parE)) { - ## check that it has exactly two entries + # check that it has exactly two entries stopifnot(identical(length(parE), 2L)) - ## and transpose to matrix with one row + # and transpose to matrix with one row parE <- t(parE) } - ## if prior weights of the beta mixture are not supplied + # 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) + # (don't need to be normalized, this is done in h_getBetamixPost) } - ## if parS is a vector => situation where there is only one component + # if parS is a vector => situation where there is only one component if (is.vector(parS)) { - ## check that it has exactly two entries + # check that it has exactly two entries stopifnot(identical(length(parS), 2L)) - ## and transpose to matrix with one row + # and transpose to matrix with one row parS <- t(parS) } - - ## if prior weights of the beta mixture are not supplied + # if prior weights of the beta mixture are not supplied if (missing(weightsS)) { weightsS <- rep(1, nrow(parS)) } - - ## now compute updated parameters for beta mixture distribution on the - ## treatment proportion + # now compute updated parameters for beta mixture distribution on the + # treatment proportion activeBetamixPost <- h_getBetamixPost(x = x, n = n, par = parE, weights = weights) - - ## now with the beta binomial mixture: - py <- with( + # now with the beta binomial mixture: + density <- with( activeBetamixPost, dbetabinomMix(x = 0:mE, m = mE, par = par, weights = weights) ) if (NmaxControl == 0) { - ## here is the only difference to predprob.R: - ## how to compute the posterior probabilities - b <- postprobDist( + # here is the only difference to predprob.R: + # how to compute the posterior probabilities + posterior <- postprobDist( x = x + c(0:mE), n = Nmax, delta = delta, @@ -142,24 +159,20 @@ predprobDist <- function(x, n, parS = parS, weightsS = weightsS ) - - ret <- structure(sum(py * (b > thetaT)), - tables = - round( - cbind( - i = c(0:mE), - py, - b, - bgttheta = (b > thetaT) - ), - 4 - ) + ret <- 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) + ) ) } else { - ## in this case also data on the SOC is available! - - ## determine remaining sample size and probabilities of response - ## counts in future SOC patients: + # in this case also data on the SOC is available! + # determine remaining sample size and probabilities of response + # counts in future SOC patients: mS <- NmaxControl - nS controlBetamixPost <- h_getBetamixPost( @@ -172,7 +185,7 @@ predprobDist <- function(x, n, dbetabinomMix(x = 0:mS, m = mS, par = par, weights = weights) ) - ## determine resulting posterior probabilities: + # determine resulting posterior probabilities: outcomesY <- x + c(0:mE) outcomesZ <- xS + c(0:mS) @@ -188,8 +201,8 @@ predprobDist <- function(x, n, for (i in seq_along(outcomesY)) { for (j in seq_along(outcomesZ)) { - ## calculate the posterior probability for this combination - ## of counts + # calculate the posterior probability for this combination + # of counts b[i, j] <- postprobDist( x = outcomesY[i], @@ -204,13 +217,13 @@ predprobDist <- function(x, n, weightsS = weightsS ) - ## what are the joint probabilities of active and control counts? - ## => because they are independent, just multiply them + # what are the joint probabilities of active and control counts? + # => because they are independent, just multiply them pyz[i, j] <- py[i] * pz[j] } } - ## should we print something? + # should we print something? ret <- structure(sum(pyz * (b > thetaT)), tables = list( @@ -221,7 +234,7 @@ predprobDist <- function(x, n, ) } - return(ret) + ret } -## todo: predprobDistFail +# todo: predprobDistFail diff --git a/inst/WORDLIST b/inst/WORDLIST index 1db7dca6..b5d955a2 100644 --- a/inst/WORDLIST +++ b/inst/WORDLIST @@ -125,6 +125,7 @@ PointMass pointmass postL postprob +postprobDist postU ppL PProbust diff --git a/man/predprobDist.Rd b/man/predprobDist.Rd index 31e0e77f..89fd25ce 100644 --- a/man/predprobDist.Rd +++ b/man/predprobDist.Rd @@ -22,37 +22,7 @@ predprobDist( ) } \arguments{ -\item{x}{number of successes (in the treatment group)} - -\item{n}{number of patients (in the treatment group)} - -\item{xS}{number of successes in the SOC group (default: 0)} - -\item{nS}{number of patients in the SOC group (default: 0)} - -\item{Nmax}{maximum number of patients at the end of the trial (in the -treatment group)} - -\item{NmaxControl}{maximum number of patients at the end of the trial in the -SOC group (default: 0)} - -\item{delta}{margin by which the response rate in the treatment group should -be better than in the SOC group (default: 0)} - -\item{relativeDelta}{see \code{\link{postprobDist}}} - -\item{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.} - -\item{weights}{the mixture weights of the beta mixture prior. Default are -uniform weights across mixture components.} - -\item{parS}{beta prior parameters in the SOC group (default: uniform)} - -\item{weightsS}{weights for the SOC group (default: uniform)} - -\item{thetaT}{threshold on the probability to be used} +\item{}{(\if{html}{\out{}}):\cr \if{html}{\out{}}} } \value{ The predictive probability, a numeric value. In addition, a @@ -65,36 +35,57 @@ Compute the predictive probability of trial success given current data. Success means that at the end of the trial the posterior probability \code{Pr(P_E > P_S + delta_0 | data) >= thetaT}. Then the predictive probability for success is: -pp = sum over i: Pr(Y=i|x,n)*I{Pr(P_E > P_S + delta|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). -Prior is P_E ~ beta(a, b), default uniform which is a beta(1,1). +\verb{pp = sum over i: Pr(Y = i | x , n)*I\{Pr(P_E > P_S + delta | x,Y=i) >= thetaT\}}, +where \code{Y} is the number of future responses in the treatment group and \code{x} is +the current number of responses in the treatment group out of \code{n}. +Prior is \code{P_E ~ beta(a, b)} and uniform which is a beta(1,1). However, also a beta mixture prior can be specified. Analogously -for P_S either a classic beta prior or a beta mixture prior can be +for \code{P_S} either a classic beta prior or a beta mixture prior can be specified. Also data on the SOC might be available. Then the predictive probability is more generally defined as -pp = sum over i, j: Pr(Y=i|x,n)*Pr(Z=j|xS, nS)*I{Pr(P_E > P_S + -delta|x,xS,Y=i,Z=j)>=thetaT} -where Z is the future number of responses in the SOC group, and xS is the +\verb{pp = sum over i, j: Pr(Y = i | x, n)*Pr(Z = j | xS, nS )*I\{Pr(P_E > P_S + delta | x,xS, Y=i, Z=j ) >= thetaT\}} +where \code{Z} is the future number of responses in the SOC group, and \code{xS} is the current number of responses in the SOC group. +In the case where \code{NmaxControl = 0}, a table with the following contents will be included in the return output : +\itemize{ +\item \code{i}: \code{Y = i}, number of future successes in \code{Nmax-n} subjects. +\item \code{density}: \code{Pr(Y = i|x)} using beta-(mixture)-binomial distribution. +\item \code{posterior}: \code{Pr(P_E > P_S + delta | x, Y = i)} using beta posterior. +\item \code{success}: indicator \code{I(b>thetaT)}. +} + A table with the following contents will be -included in the \code{tables} attribute of the return value +included in the return value (in the case that NmaxControl is zero): 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_S + delta | x, Y=i) bgttheta: indicator I(b>thetaT) -If NmaxControl is not zero, i.e., when data on the control treatment +If \code{NmaxControl} is not zero, i.e., when data on the control treatment is available in this trial, then a list with will be included with the following elements: pyz: matrix with the probabilities Pr(Y=i, Z=j | x, xS) b: matrix with Pr(P_E > P_S + delta | x, xS, Y=i, Z=j) bgttheta: matrix of indicators I(b>thetaT) } +\note{ +Delta from postprobDist : + +The desired improvement is denoted as \code{delta}. There are two options in using \code{delta}. +The absolute case when \code{relativeDelta = FALSE} and relative as when \code{relativeDelta = TRUE}. +\enumerate{ +\item The absolute case is when we define an absolute delta, greater than \code{P_S}, +the response rate of the \code{SOC} group such that +the posterior is \code{Pr(P_E > P_S + delta | data)}. +\item In the relative case, we suppose that the treatment group's +response rate is assumed to be greater than \code{P_S + (1-P_S) * delta} such that +the posterior is \code{Pr(P_E > P_S + (1 - P_S) * delta | data)}. +} +} \examples{ ## The original Lee and Liu (Table 1) example: ## Nmax=40, x=16, n=23, beta(0.6,0.4) prior distribution, From 5f9b5b8ea27cd8f52c1d75f060ffff23d97705a4 Mon Sep 17 00:00:00 2001 From: Audrey Yeo Date: Tue, 5 Dec 2023 13:29:00 +0100 Subject: [PATCH 03/15] update --- R/predprobDist.R | 102 ++++++------ design/design-doc_predprobDist.Rmd | 238 +++++++++++++++++++++++++++ examples/predprobDist.R | 66 ++++---- inst/WORDLIST | 1 + man/predprobDist.Rd | 107 +++++++----- myEnvironment.RData | Bin 0 -> 130 bytes tests/{ => testthat}/test-predprob.R | 0 7 files changed, 391 insertions(+), 123 deletions(-) create mode 100644 design/design-doc_predprobDist.Rmd create mode 100644 myEnvironment.RData rename tests/{ => testthat}/test-predprob.R (100%) diff --git a/R/predprobDist.R b/R/predprobDist.R index 67a11fdd..6daab5cb 100644 --- a/R/predprobDist.R +++ b/R/predprobDist.R @@ -61,27 +61,41 @@ NULL #' response rate is assumed to be greater than `P_S + (1-P_S) * delta` such that #' the posterior is `Pr(P_E > P_S + (1 - P_S) * delta | data)`. #' -#' @typed x number of successes (in the treatment group) -#' @typed n number of patients (in the treatment group) -#' @typed xS number of successes in the SOC group (default: 0) -#' @typed nS number of patients in the SOC group (default: 0) -#' @typed Nmax maximum number of patients at the end of the trial (in the -#' treatment group) -#' @typed NmaxControl maximum number of patients at the end of the trial in the -#' SOC group (default: 0) -#' @typed delta margin by which the response rate in the treatment group should -#' be better than in the SOC group (default: 0) -#' @typed relativeDelta see \code{\link{postprobDist}} -#' @typed parE the beta parameters matrix, with K rows and 2 columns, +#' @typed x : number +#' number of successes in the treatment group at interm +#' @typed n : number +#' number of patients in the treatment group at interim +#' @typed xS : number +#' number of successes in the SOC group at interim +#' @typed nS : number +#' number of patients in the SOC group +#' @typed Nmax : number +#' maximum number of patients in the treatment group at the end of the trial +#' @typed NmaxControl : number +#' maximum number of patients at the end of the trial in the +#' SOC group +#' @typed delta : +#' margin by which the response rate in the treatment group should +#' be better than in the SOC group +#' @typed relativeDelta : +#' see `[postprobDist()]` +#' @typed 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. -#' @typed weights the mixture weights of the beta mixture prior. Default are +#' @typed weights : +#' the mixture weights of the beta mixture prior. Default are #' uniform weights across mixture components. -#' @typed parS beta prior parameters in the SOC group (default: uniform) -#' @typed weightsS weights for the SOC group (default: uniform) -#' @typed thetaT threshold on the probability to be used -#' @return The predictive probability, a numeric value. In addition, a -#' list called \code{tables} is returned as attribute of the returned number. +#' @typed parS : +#' beta prior parameters in the SOC group +#' @typed weightsS : +#' weights for the SOC group +#' @typed thetaT : +#' threshold on the probability to be used +#' @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), @@ -90,8 +104,10 @@ NULL #' @example examples/predprobDist.R #' @export predprobDist <- function(x, n, - xS = 0, nS = 0, - Nmax, NmaxControl = 0, + xS = 0, + nS = 0, + Nmax, + NmaxControl = 0, delta = 0, relativeDelta = FALSE, parE = c(a = 1, b = 1), @@ -106,30 +122,24 @@ predprobDist <- function(x, n, x <= n, xS <= nS ) - # remaining active patients to be seen: mE <- 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 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) } - # if parS is a vector => situation where there is only one component if (is.vector(parS)) { # check that it has exactly two entries stopifnot(identical(length(parS), 2L)) - # and transpose to matrix with one row parS <- t(parS) } @@ -141,15 +151,14 @@ predprobDist <- function(x, n, # treatment proportion activeBetamixPost <- h_getBetamixPost(x = x, n = n, par = parE, weights = weights) # now with the beta binomial mixture: - density <- with( + py <- with( activeBetamixPost, dbetabinomMix(x = 0:mE, m = mE, par = par, weights = weights) ) - if (NmaxControl == 0) { # here is the only difference to predprob.R: # how to compute the posterior probabilities - posterior <- postprobDist( + b <- postprobDist( x = x + c(0:mE), n = Nmax, delta = delta, @@ -160,13 +169,13 @@ predprobDist <- function(x, n, weightsS = weightsS ) ret <- list( - result = sum(density * (posterior > thetaT)), + result = sum(py * (b > thetaT)), table = data.frame( - counts = c(0:m), - cumul_counts = n + (0:m), - density = round(density, 4), - posterior = posterior, - success = (posterior > thetaT) + counts = c(0:mE), + # cumul_counts = xS + (0:NmaxControl), + density = py, + posterior = b, + success = (b > thetaT) ) ) } else { @@ -174,21 +183,14 @@ predprobDist <- function(x, n, # determine remaining sample size and probabilities of response # counts in future SOC patients: mS <- NmaxControl - nS - - controlBetamixPost <- h_getBetamixPost( - x = xS, n = nS, par = parS, - weights = weightsS - ) - + controlBetamixPost <- h_getBetamixPost(x = xS, n = nS, par = parS, weights = weightsS) pz <- with( controlBetamixPost, dbetabinomMix(x = 0:mS, m = mS, par = par, weights = weights) ) - # determine resulting posterior probabilities: outcomesY <- x + c(0:mE) - outcomesZ <- xS + c(0:mS) - + outcomesZ <- xS + c(0:mS) # 15 more chances to get success counts pyz <- b <- matrix( nrow = 1 + mE, ncol = 1 + mS, @@ -198,8 +200,7 @@ predprobDist <- function(x, n, 0:mS ) ) - - for (i in seq_along(outcomesY)) { + for (i in seq_along(outcomesY)) { # outside? for (j in seq_along(outcomesZ)) { # calculate the posterior probability for this combination # of counts @@ -216,24 +217,21 @@ predprobDist <- function(x, n, parS = parS, weightsS = weightsS ) - # what are the joint probabilities of active and control counts? # => because they are independent, just multiply them - pyz[i, j] <- py[i] * pz[j] + pyz[i, j] <- py[i] * pz[j] # this is the matrix that Daniel was talking about } } - - # should we print something? + # should we print something? predprob part ret <- structure(sum(pyz * (b > thetaT)), tables = list( pyz = pyz, b = b, - bgttheta = (b > thetaT) + success = (b > thetaT) ) ) } - ret } diff --git a/design/design-doc_predprobDist.Rmd b/design/design-doc_predprobDist.Rmd new file mode 100644 index 00000000..b66edc79 --- /dev/null +++ b/design/design-doc_predprobDist.Rmd @@ -0,0 +1,238 @@ +--- +title: "Design Document" +output: html_document +date: "2023-12-04" +--- + +## Example inputs / Sanity check + +```{r setup, include=FALSE} +x <- 16 +n <- 23 +xS <- 5 +nS <- 10 +Nmax <- 40 +NmaxControl <- 20 +delta <- 0 +thetaT <- 0.6 +parE <- rbind(c(1, 1), c(50, 10)) +parS <- rbind(c(1, 1), c(20, 40)) +weight <- c(2, 1) +weightS <- c(2, 1) +relativeDelta <- FALSE + +# predprobDist( +# x = x, +# n = n, +# xS = xS, +# nS = nS, +# Nmax = Nmax, +# NmaxControl = NmaxControl, +# delta = delta, +# thetaT = thetaT, +# parE = parE, +# parS = parS +# ) + +# predprobDist( +# x = 16, n = 23, Nmax = 40, delta = 0, thetaT = 0.9, +# parE = c(0.6, 0.4), parS = c(7, 11) +# ) +``` + +```{r} +predprobDist( + x = 16, n = 23, Nmax = 40, delta = 0, thetaT = 0.9, + parE = c(0.6, 0.4), parS = c(7, 11) +) +``` +## R Markdown + +# helper if NmaxControl == 0, or when there is no control +```{r} +h_get_predproblist_control <- function(x, + n, + delta, + relativeDelta, + parE, + weights, + parS, + weightsS, + thetaT, + density, + mE) { + posterior <- postprobDist( + x = x + c(0:mE), + n = Nmax, + delta = delta, + relativeDelta = relativeDelta, + parE = parE, + weights = weights, + parS = parS, + weightsS = weightsS + ) + list( + result = sum(density * (posterior > thetaT)), + table = data.frame( + counts = c(0:mE), + cumul_counts = n + (0:mE), + density = round(density, 4), + posterior = posterior, + success = (posterior > thetaT) + ) + ) +} +``` + + +```{r cars} +predprobDist <- function(x, n, + xS = 0, + nS = 0, + Nmax, + NmaxControl = 0, + delta = 0, + relativeDelta = FALSE, + parE = c(a = 1, b = 1), + weights, + parS = c(a = 1, b = 1), + weightsS, + thetaT) { + # ensure reasonable numbers + stopifnot( + n <= Nmax, + nS <= NmaxControl, + x <= n, + xS <= nS + ) + # remaining active patients to be seen: + mE <- 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 + 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) + } + # if parS is a vector => situation where there is only one component + if (is.vector(parS)) { + # check that it has exactly two entries + stopifnot(identical(length(parS), 2L)) + # and transpose to matrix with one row + parS <- t(parS) + } + # if prior weights of the beta mixture are not supplied + if (missing(weightsS)) { + weightsS <- rep(1, nrow(parS)) + } + # now compute updated parameters for beta mixture distribution on the + # treatment proportion + activeBetamixPost <- h_getBetamixPost(x = x, n = n, par = parE, weights = weights) + # now with the beta binomial mixture: + py <- with( + activeBetamixPost, + dbetabinomMix(x = 0:mE, m = mE, par = par, weights = weights) + ) + if (NmaxControl == 0) { + # here is the only difference to predprob.R: + # how to compute the posterior probabilities + # posterior <- postprobDist( + # x = x + c(0:mE), + # n = Nmax, + # delta = delta, + # relativeDelta = relativeDelta, + # parE = parE, + # weights = weights, + # parS = parS, + # weightsS = weightsS + # ) + # ret <- list( + # result = sum(py * (posterior > thetaT)), + # table = data.frame( + # counts = c(0:mE), + # cumul_counts = n + (0:mE), + # density = py, + # posterior = posterior, + # success = (posterior > thetaT) + # ) + # ) + h_get_predproblist_control( + x = x, + n = n, + delta = delta, + relativeDelta = relativeDelta, + parE = parE, + weights = weights, + parS = parS, + weightsS = weightsS, + thetaT = thetaT, + density = py, + mE = mE + ) + } else { + # in this case also data on the SOC is available! + # determine remaining sample size and probabilities of response + # counts in future SOC patients: + mS <- NmaxControl - nS + controlBetamixPost <- h_getBetamixPost(x = xS, n = nS, par = parS, weights = weightsS) + pz <- with( + controlBetamixPost, + dbetabinomMix(x = 0:mS, m = mS, par = par, weights = weights) + ) + # determine resulting posterior probabilities: + outcomesY <- x + c(0:mE) + outcomesZ <- xS + c(0:mS) + pyz <- b <- matrix( + nrow = 1 + mE, + ncol = 1 + mS, + dimnames = + list( + 0:mE, + 0:mS + ) + ) + for (i in seq_along(outcomesY)) { # outside? + for (j in seq_along(outcomesZ)) { + # calculate the posterior probability for this combination + # of counts + b[i, j] <- # postprobDist part + postprobDist( + x = outcomesY[i], + n = Nmax, + xS = outcomesZ[j], + nS = NmaxControl, + delta = delta, + relativeDelta = relativeDelta, + parE = parE, + weights = weights, + parS = parS, + weightsS = weightsS + ) + # what are the joint probabilities of active and control counts? + # => because they are independent, just multiply them + pyz[i, j] <- py[i] * pz[j] + } + } + # should we print something? predprob part + ret <- structure(sum(pyz * (b > thetaT)), + tables = + list( + pyz = pyz, + b = b, + success = (b > thetaT) + ) + ) + } + ret +} +``` + + + +```{r pressure, echo=FALSE} +``` diff --git a/examples/predprobDist.R b/examples/predprobDist.R index 44b35769..90ce5ca9 100644 --- a/examples/predprobDist.R +++ b/examples/predprobDist.R @@ -1,53 +1,53 @@ -## 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. +# So the predictive probability is 56.6%. 12 responses +# in the remaining 17 patients are required for success. -## Now to the modified example, where the control response -## rate p has a beta(7, 11) distribution. -## For example if the 60% threshold stems from a trial -## with 10 patients, of which 6 were observed with response. +# Now to the modified example, where the control response +# rate p has a beta(7, 11) distribution. +# For example if the 60% threshold (control) stems from a trial +# with 10 patients, of which 6 were observed with response. predprobDist( x = 16, n = 23, Nmax = 40, delta = 0, thetaT = 0.9, parE = c(0.6, 0.4), parS = c(7, 11) ) -## Now only 11 responses are needed in the remaining 17 patients -## in order to have a successful trial. -## Also, the predictive probability for a Go decision -## is now 70.8% instead of merely 56.6%. +# Now only7 responses are needed in the remaining 17 patients +# in order to have a successful trial. +# Also, the predictive probability for a Go decision +# is now 97.78% instead of merely 56.6%. -## If the response threshold of 60% stems from an SOC estimate (e.g. 50% = 5/10) -## and a margin of e.g. 10%, then the margin can be specified with -## the delta argument: +# If the response threshold of 60% stems from an SOC estimate (e.g. 50% = 5/10) +# and a margin of e.g. 10%, then the margin can be specified with +# the delta argument: predprobDist( x = 16, n = 23, Nmax = 40, delta = 0.1, thetaT = 0.9, parE = c(1, 1), parS = c(6, 11) ) -## This again changes the result: Only 10 future responses are required, -## the predictive probability of success is now 80%. +# This again changes the result: Only 10 future responses are required, +# the predictive probability of success is now 80%. -## Next we will use a beta mixture prior, with a weight of 1/3 for an -## informative beta(50, 10) prior, on the experimental response rate: +# Next we will use a beta mixture prior, with a weight of 1/3 for an +# informative beta(50, 10) prior, on the experimental response rate: predprobDist( x = 16, n = 23, Nmax = 40, delta = 0.1, thetaT = 0.9, parE = rbind(c(1, 1), c(50, 10)), weights = c(2, 1), parS = c(6, 11) ) -## Now still 10 future responses are required for success, but the predictive -## success probability is higher (87.2%). +# Now still 10 future responses are required for success, but the predictive +# success probability is higher (87.2%). -## We can also have additional controls in our trial. For example -## assume that in the interim analysis we had 10 control patients observed -## with 5 responses out of 20 in total at the end. -## We also had historical control data (say 60 patients of -## which 20 responded), but in order to be robust against changes in the response +# We can also have additional controls in our trial. For example +# assume that in the interim analysis we had 10 control patients observed +# with 5 responses out of 20 in total at the end. +# We also had historical control data (say 60 patients of +# which 20 responded), but in order to be robust against changes in the response # rate over time, we will again use a beta mixture, putting only 1/3 of weight -## on the historical controls: +# on the historical controls: predprobDist( x = 16, n = 23, xS = 5, nS = 10, Nmax = 40, NmaxControl = 20, @@ -57,8 +57,8 @@ predprobDist( parS = rbind(c(1, 1), c(20, 40)), weightsS = c(2, 1) ) -## Now we only have 59.9% predictive probability of success. The bgttheta matrix -## lists the cases where we would have success, for all possible combinations -## of future responses in the control and experimental arms. We see that the -## success threshold for experimental responses depends on the number of -## control responses. +# Now we only have 59.9% predictive probability of success. The bgttheta matrix +# lists the cases where we would have success, for all possible combinations +# of future responses in the control and experimental arms. We see that the +# success threshold for experimental responses depends on the number of +# control responses. diff --git a/inst/WORDLIST b/inst/WORDLIST index b5d955a2..2cd5aa1f 100644 --- a/inst/WORDLIST +++ b/inst/WORDLIST @@ -91,6 +91,7 @@ increasement inefficacious Integrand integrations +interm LciU ldots leq diff --git a/man/predprobDist.Rd b/man/predprobDist.Rd index 89fd25ce..6389a81e 100644 --- a/man/predprobDist.Rd +++ b/man/predprobDist.Rd @@ -22,11 +22,42 @@ predprobDist( ) } \arguments{ -\item{}{(\if{html}{\out{}}):\cr \if{html}{\out{}}} +\item{x}{(\code{number}):\cr number of successes in the treatment group at interm} + +\item{n}{(\code{number}):\cr number of patients in the treatment group at interim} + +\item{xS}{(\code{number}):\cr number of successes in the SOC group at interim} + +\item{nS}{(\code{number}):\cr number of patients in the SOC group} + +\item{Nmax}{(\code{number}):\cr maximum number of patients in the treatment group at the end of the trial} + +\item{NmaxControl}{(\code{number}):\cr maximum number of patients at the end of the trial in the +SOC group} + +\item{delta}{(``):\cr margin by which the response rate in the treatment group should +be better than in the SOC group} + +\item{relativeDelta}{(``):\cr see \verb{[postprobDist()]}} + +\item{parE}{(``):\cr the beta parameters matrix, with K rows and 2 columns, +corresponding to the beta parameters of the K components. default is a +uniform prior.} + +\item{weights}{(``):\cr the mixture weights of the beta mixture prior. Default are +uniform weights across mixture components.} + +\item{parS}{(``):\cr beta prior parameters in the SOC group} + +\item{weightsS}{(``):\cr weights for the SOC group} + +\item{thetaT}{(``):\cr threshold on the probability to be used} } \value{ -The predictive probability, a numeric value. In addition, a -list called \code{tables} is returned as attribute of the returned number. +A \code{list} is returned with names \code{result} for predictive probability and +\code{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. } \description{ \ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#experimental}{\figure{lifecycle-experimental.svg}{options: alt='[Experimental]'}}}{\strong{[Experimental]}} @@ -87,56 +118,56 @@ the posterior is \code{Pr(P_E > P_S + (1 - P_S) * delta | data)}. } } \examples{ -## 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. +# So the predictive probability is 56.6\%. 12 responses +# in the remaining 17 patients are required for success. -## Now to the modified example, where the control response -## rate p has a beta(7, 11) distribution. -## For example if the 60\% threshold stems from a trial -## with 10 patients, of which 6 were observed with response. +# Now to the modified example, where the control response +# rate p has a beta(7, 11) distribution. +# For example if the 60\% threshold (control) stems from a trial +# with 10 patients, of which 6 were observed with response. predprobDist( x = 16, n = 23, Nmax = 40, delta = 0, thetaT = 0.9, parE = c(0.6, 0.4), parS = c(7, 11) ) -## Now only 11 responses are needed in the remaining 17 patients -## in order to have a successful trial. -## Also, the predictive probability for a Go decision -## is now 70.8\% instead of merely 56.6\%. - -## If the response threshold of 60\% stems from an SOC estimate (e.g. 50\% = 5/10) -## and a margin of e.g. 10\%, then the margin can be specified with -## the delta argument: +# Now only7 responses are needed in the remaining 17 patients +# in order to have a successful trial. +# Also, the predictive probability for a Go decision +# is now 97.78\% instead of merely 56.6\%. + +# If the response threshold of 60\% stems from an SOC estimate (e.g. 50\% = 5/10) +# and a margin of e.g. 10\%, then the margin can be specified with +# the delta argument: predprobDist( x = 16, n = 23, Nmax = 40, delta = 0.1, thetaT = 0.9, parE = c(1, 1), parS = c(6, 11) ) -## This again changes the result: Only 10 future responses are required, -## the predictive probability of success is now 80\%. +# This again changes the result: Only 10 future responses are required, +# the predictive probability of success is now 80\%. -## Next we will use a beta mixture prior, with a weight of 1/3 for an -## informative beta(50, 10) prior, on the experimental response rate: +# Next we will use a beta mixture prior, with a weight of 1/3 for an +# informative beta(50, 10) prior, on the experimental response rate: predprobDist( x = 16, n = 23, Nmax = 40, delta = 0.1, thetaT = 0.9, parE = rbind(c(1, 1), c(50, 10)), weights = c(2, 1), parS = c(6, 11) ) -## Now still 10 future responses are required for success, but the predictive -## success probability is higher (87.2\%). - -## We can also have additional controls in our trial. For example -## assume that in the interim analysis we had 10 control patients observed -## with 5 responses out of 20 in total at the end. -## We also had historical control data (say 60 patients of -## which 20 responded), but in order to be robust against changes in the response +# Now still 10 future responses are required for success, but the predictive +# success probability is higher (87.2\%). + +# We can also have additional controls in our trial. For example +# assume that in the interim analysis we had 10 control patients observed +# with 5 responses out of 20 in total at the end. +# We also had historical control data (say 60 patients of +# which 20 responded), but in order to be robust against changes in the response # rate over time, we will again use a beta mixture, putting only 1/3 of weight -## on the historical controls: +# on the historical controls: predprobDist( x = 16, n = 23, xS = 5, nS = 10, Nmax = 40, NmaxControl = 20, @@ -146,11 +177,11 @@ predprobDist( parS = rbind(c(1, 1), c(20, 40)), weightsS = c(2, 1) ) -## Now we only have 59.9\% predictive probability of success. The bgttheta matrix -## lists the cases where we would have success, for all possible combinations -## of future responses in the control and experimental arms. We see that the -## success threshold for experimental responses depends on the number of -## control responses. +# Now we only have 59.9\% predictive probability of success. The bgttheta matrix +# lists the cases where we would have success, for all possible combinations +# of future responses in the control and experimental arms. We see that the +# success threshold for experimental responses depends on the number of +# control responses. } \references{ Lee, J. J., & Liu, D. D. (2008). A predictive probability diff --git a/myEnvironment.RData b/myEnvironment.RData new file mode 100644 index 0000000000000000000000000000000000000000..167b2c222cae55323fbd0aea84d2f77a527cf746 GIT binary patch literal 130 zcmb2|=3oE=X6~X+gJ)e2k`fXU(h?HWQWDZwjU*$So$r+BN=QV^XpnJZ6Jzt6(A0E9 ziH%LoPGRy<$96^@85t>sElgA1^F3Ug?R)YJYtSBtGc7Iu7iuQ!`=9kWuc66i7s9Y34Oru;p&Z_$xmOZ60V85a2+_S*=w8UP5hFsJ|k literal 0 HcmV?d00001 diff --git a/tests/test-predprob.R b/tests/testthat/test-predprob.R similarity index 100% rename from tests/test-predprob.R rename to tests/testthat/test-predprob.R From 35a0408e353f40fc73d2a9e35f1723b1de3ac2f1 Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Tue, 5 Dec 2023 17:30:58 +0100 Subject: [PATCH 04/15] fix by adding explicit function values for boundaries --- R/dbetabinom.R | 5 ++++- tests/testthat/test-dbetabinom.R | 16 ++++++++++++++++ 2 files changed, 20 insertions(+), 1 deletion(-) diff --git a/R/dbetabinom.R b/R/dbetabinom.R index 819c99a2..89443b6a 100644 --- a/R/dbetabinom.R +++ b/R/dbetabinom.R @@ -192,11 +192,14 @@ qbetaMix <- function(p, par, weights, lower.tail = TRUE) { f <- function(pi) { pbetaMix(q = pi, par = par, weights = weights, lower.tail = lower.tail) - p } - unirootResult <- uniroot(f, lower = 0, upper = 1) + # Note: we give the lower and upper function values here in order to avoid problems for + # p = 0 or p = 1. + unirootResult <- uniroot(f, lower = 0, upper = 1, f.lower = -p, f.upper = 1 - p) if (unirootResult$iter < 0) { NA } else { assert_number(unirootResult$root) + assert_true(all.equal(f(unirootResult$root), 0)) unirootResult$root } } diff --git a/tests/testthat/test-dbetabinom.R b/tests/testthat/test-dbetabinom.R index dae311f5..e36b8159 100644 --- a/tests/testthat/test-dbetabinom.R +++ b/tests/testthat/test-dbetabinom.R @@ -129,6 +129,22 @@ test_that("qbetaMix gives a number result with beta-mixture", { expect_numeric(result) }) +test_that("qbetaMix works as expected for edge cases", { + result <- qbetaMix( + p = 0, + par = rbind(c(0.2, 0.4), c(1, 1)), + weights = c(0.6, 0.4) + ) + expect_identical(result, 0) + + result <- qbetaMix( + p = 1, + par = rbind(c(0.2, 0.4), c(1, 1)), + weights = c(0.6, 0.4) + ) + expect_identical(result, 1) +}) + # dbetaMix ---- test_that("dbetaMix gives the correct result with a 1 mixture component", { From 29df4d5cc1cff3a4a0647aaf83d206ea85fe5e9a Mon Sep 17 00:00:00 2001 From: Audrey Yeo Date: Tue, 5 Dec 2023 22:03:44 +0100 Subject: [PATCH 05/15] update --- design/design-doc_predprobDist.Rmd | 206 +++++++++++++++-------------- tests/testthat/test-predprobDist.R | 30 +++++ 2 files changed, 135 insertions(+), 101 deletions(-) create mode 100644 tests/testthat/test-predprobDist.R diff --git a/design/design-doc_predprobDist.Rmd b/design/design-doc_predprobDist.Rmd index b66edc79..4247def1 100644 --- a/design/design-doc_predprobDist.Rmd +++ b/design/design-doc_predprobDist.Rmd @@ -15,29 +15,22 @@ Nmax <- 40 NmaxControl <- 20 delta <- 0 thetaT <- 0.6 -parE <- rbind(c(1, 1), c(50, 10)) -parS <- rbind(c(1, 1), c(20, 40)) -weight <- c(2, 1) -weightS <- c(2, 1) +parE <- c(1, 1) +parS <- c(1, 1) +weights <- c(2, 1) +weightsS <- c(2, 1) relativeDelta <- FALSE -# predprobDist( -# x = x, -# n = n, -# xS = xS, -# nS = nS, -# Nmax = Nmax, -# NmaxControl = NmaxControl, -# delta = delta, -# thetaT = thetaT, -# parE = parE, -# parS = parS -# ) - -# predprobDist( -# x = 16, n = 23, Nmax = 40, delta = 0, thetaT = 0.9, -# parE = c(0.6, 0.4), parS = c(7, 11) -# ) +# helper test +mE <- Nmax - n +parE <- t(parE) +weights <- rep(1, nrow(parE)) +parS <- t(parS) +weightsS <- rep(1, nrow(parS)) + +activeBetamixPost <- h_getBetamixPost(x = x, n = n, par = parE, weights = weights) +# now with the beta binomial mixture: +py <- with(activeBetamixPost, dbetabinomMix(x = 0:mE, m = mE, par = par, weights = weights) ``` ```{r} @@ -50,18 +43,21 @@ predprobDist( # helper if NmaxControl == 0, or when there is no control ```{r} -h_get_predproblist_control <- function(x, - n, - delta, - relativeDelta, - parE, - weights, - parS, - weightsS, - thetaT, - density, - mE) { - posterior <- postprobDist( +h_get_table_singlearm <- function( + n, + x, + mE, + delta, + relativeDelta, + parE, + weights , + parS, + weightsS , + thetaT , + b, + py, + Nmax) { + b <- postprobDist( x = x + c(0:mE), n = Nmax, delta = delta, @@ -71,24 +67,43 @@ h_get_predproblist_control <- function(x, parS = parS, weightsS = weightsS ) - list( - result = sum(density * (posterior > thetaT)), - table = data.frame( - counts = c(0:mE), - cumul_counts = n + (0:mE), - density = round(density, 4), - posterior = posterior, - success = (posterior > thetaT) - ) + + ret <- structure(sum(py * (b > thetaT)), + tables = + round( + cbind( + i = c(0:mE), + py, + b, + bgttheta = (b > thetaT) + ), + 4 + ) ) + } + +``` + +```{r} +h_get_predproblist_control(x = x, + Nmax = Nmax, + delta = delta, + relativeDelta = relativeDelta, + parE = parE, + weights = weights, + parS = parS, + weightsS = weightsS, + thetaT = thetaT, + density = 0.77, + mE = 17) ``` + ```{r cars} predprobDist <- function(x, n, - xS = 0, - nS = 0, + xS = 0, nS = 0, Nmax, NmaxControl = 0, delta = 0, @@ -98,50 +113,48 @@ predprobDist <- function(x, n, parS = c(a = 1, b = 1), weightsS, thetaT) { - # ensure reasonable numbers stopifnot( n <= Nmax, nS <= NmaxControl, x <= n, xS <= nS ) - # remaining active patients to be seen: mE <- 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 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) } - # if parS is a vector => situation where there is only one component if (is.vector(parS)) { - # check that it has exactly two entries stopifnot(identical(length(parS), 2L)) - # and transpose to matrix with one row parS <- t(parS) } - # if prior weights of the beta mixture are not supplied if (missing(weightsS)) { weightsS <- rep(1, nrow(parS)) } - # now compute updated parameters for beta mixture distribution on the - # treatment proportion activeBetamixPost <- h_getBetamixPost(x = x, n = n, par = parE, weights = weights) - # now with the beta binomial mixture: py <- with( activeBetamixPost, dbetabinomMix(x = 0:mE, m = mE, par = par, weights = weights) ) if (NmaxControl == 0) { - # here is the only difference to predprob.R: - # how to compute the posterior probabilities - # posterior <- postprobDist( + h_get_table_singlearm( + x = x, + n = Nmax, + mE = mE, + delta = delta, + relativeDelta = relativeDelta, + parE = parE, + weights = weights, + parS = parS, + weightsS = weightsS, + thetaT = thetaT, + b = b, + py = py, + Nmax) + # b <- postprobDist( # x = x + c(0:mE), # n = Nmax, # delta = delta, @@ -151,40 +164,30 @@ predprobDist <- function(x, n, # parS = parS, # weightsS = weightsS # ) - # ret <- list( - # result = sum(py * (posterior > thetaT)), - # table = data.frame( - # counts = c(0:mE), - # cumul_counts = n + (0:mE), - # density = py, - # posterior = posterior, - # success = (posterior > thetaT) - # ) + # + # ret <- structure(sum(py * (b > thetaT)), + # tables = + # round( + # cbind( + # i = c(0:mE), + # py, + # b, + # bgttheta = (b > thetaT) + # ), + # 4 + # ) # ) - h_get_predproblist_control( - x = x, - n = n, - delta = delta, - relativeDelta = relativeDelta, - parE = parE, - weights = weights, - parS = parS, - weightsS = weightsS, - thetaT = thetaT, - density = py, - mE = mE - ) } else { - # in this case also data on the SOC is available! - # determine remaining sample size and probabilities of response - # counts in future SOC patients: mS <- NmaxControl - nS - controlBetamixPost <- h_getBetamixPost(x = xS, n = nS, par = parS, weights = weightsS) + controlBetamixPost <- h_getBetamixPost( + x = xS, n = nS, par = parS, + weights = weightsS + ) pz <- with( controlBetamixPost, dbetabinomMix(x = 0:mS, m = mS, par = par, weights = weights) ) - # determine resulting posterior probabilities: + ## determine resulting posterior probabilities: outcomesY <- x + c(0:mE) outcomesZ <- xS + c(0:mS) pyz <- b <- matrix( @@ -196,11 +199,9 @@ predprobDist <- function(x, n, 0:mS ) ) - for (i in seq_along(outcomesY)) { # outside? + for (i in seq_along(outcomesY)) { for (j in seq_along(outcomesZ)) { - # calculate the posterior probability for this combination - # of counts - b[i, j] <- # postprobDist part + b[i, j] <- postprobDist( x = outcomesY[i], n = Nmax, @@ -213,26 +214,29 @@ predprobDist <- function(x, n, parS = parS, weightsS = weightsS ) - # what are the joint probabilities of active and control counts? - # => because they are independent, just multiply them pyz[i, j] <- py[i] * pz[j] } } - # should we print something? predprob part ret <- structure(sum(pyz * (b > thetaT)), - tables = - list( - pyz = pyz, - b = b, - success = (b > thetaT) - ) + tables = + list( + pyz = pyz, + b = b, + bgttheta = (b > thetaT) + ) ) } ret } + ``` ```{r pressure, echo=FALSE} +predprobDist( + x = 16, n = 23, Nmax = 40, delta = 0.1, thetaT = 0.9, + parE = rbind(c(1, 1), c(50, 10)), + weights = c(2, 1), parS = c(6, 11) +) ``` diff --git a/tests/testthat/test-predprobDist.R b/tests/testthat/test-predprobDist.R new file mode 100644 index 00000000..8776c03a --- /dev/null +++ b/tests/testthat/test-predprobDist.R @@ -0,0 +1,30 @@ +x <- 16 +n <- 23 +xS <- 5 +nS <- 10 +Nmax <- 40 +NmaxControl <- 20 +delta <- 0 +thetaT <- 0.6 +parE <- c(1, 1) +parS <- c(1, 1) +weights <- c(2, 1) +weightsS <- c(2, 1) +relativeDelta <- FALSE +mE <- Nmax - n +parE <- t(parE) +weights <- rep(1, nrow(parE)) +parS <- t(parS) +weightsS <- rep(1, nrow(parS)) + +h_get_predproblist_control(x = x, + n = Nmax, + delta = delta, + relativeDelta = relativeDelta, + parE = parE, + weights = weights, + parS = parS, + weightsS = weightsS, + thetaT = thetaT, + density = 0.77, + mE = 17) From a888d482c7d39e8e9e1801ec543cce8b72828d5b Mon Sep 17 00:00:00 2001 From: Audrey Yeo Date: Wed, 6 Dec 2023 11:36:26 +0100 Subject: [PATCH 06/15] helper functions x2 --- design/design-doc_predprobDist.Rmd | 246 ++++++++++++++--------------- tests/testthat/test-predprobDist.R | 57 +++++-- 2 files changed, 168 insertions(+), 135 deletions(-) diff --git a/design/design-doc_predprobDist.Rmd b/design/design-doc_predprobDist.Rmd index 4247def1..0d7af2f5 100644 --- a/design/design-doc_predprobDist.Rmd +++ b/design/design-doc_predprobDist.Rmd @@ -30,7 +30,7 @@ weightsS <- rep(1, nrow(parS)) activeBetamixPost <- h_getBetamixPost(x = x, n = n, par = parE, weights = weights) # now with the beta binomial mixture: -py <- with(activeBetamixPost, dbetabinomMix(x = 0:mE, m = mE, par = par, weights = weights) +density_y <- with(activeBetamixPost, dbetabinomMix(x = 0:mE, m = mE, par = par, weights = weights)) ``` ```{r} @@ -41,23 +41,20 @@ predprobDist( ``` ## R Markdown -# helper if NmaxControl == 0, or when there is no control +# helper function if NmaxControl == 0, or when there is no control ```{r} -h_get_table_singlearm <- function( - n, - x, - mE, - delta, - relativeDelta, - parE, - weights , - parS, - weightsS , - thetaT , - b, - py, - Nmax) { - b <- postprobDist( +h_get_predproblist_nocontrol <- function(x, + n, + delta, + relativeDelta, + parE, + weights, + parS, + weightsS, + thetaT, + density, + mE) { + posterior_x <- postprobDist( x = x + c(0:mE), n = Nmax, delta = delta, @@ -67,40 +64,111 @@ h_get_table_singlearm <- function( parS = parS, weightsS = weightsS ) - - ret <- structure(sum(py * (b > thetaT)), - tables = - round( - cbind( - i = c(0:mE), - py, - b, - bgttheta = (b > thetaT) - ), - 4 - ) + list( + result = sum(density * (posterior_x > thetaT)), + table = data.frame( + counts = c(0:mE), + cumul_counts = n + (0:mE), + density = round(density, 4), + posterior = posterior_x, + success = (posterior_x > thetaT) + ) ) - } - ``` +# Testing helper function if NmaxControl == 0, or when there is no control ```{r} -h_get_predproblist_control(x = x, - Nmax = Nmax, - delta = delta, - relativeDelta = relativeDelta, - parE = parE, - weights = weights, - parS = parS, - weightsS = weightsS, - thetaT = thetaT, - density = 0.77, - mE = 17) +# testing +h_get_predproblist_nocontrol( + x = x, + n = Nmax, + delta = delta, + relativeDelta = relativeDelta, + parE = parE, + weights = weights, + parS = parS, + weightsS = weightsS, + thetaT = thetaT, + density = 0.77, + mE = 17 +) ``` - - - +# helper function 1 if NmaxControl != 0 +```{r} +activeBetamixPost <- h_getBetamixPost(x = 16, n = 23, par = t(c(1, 1)), weights = 1) +density_y <- with(activeBetamixPost, dbetabinomMix(x = 0:17, m = 17, par = t(c(1, 1)), weights = 1)) + +h_getlistpredprob <- function(NmaxControl, nS, xS, par, weightsS, mS, x, mE, density_y) { + mS <- NmaxControl - nS + controlBetamixPost <- h_getBetamixPost( + x = xS, + n = nS, + par = parS, + weights = weightsS + ) + pz <- with( + controlBetamixPost, + dbetabinomMix(x = 0:mS, m = mS, par = par, weights = weights) + ) + ## determine resulting posterior probabilities: + outcomesY <- x + c(0:mE) + outcomesZ <- xS + c(0:mS) + pyz <- b <- matrix( + nrow = 1 + mE, + ncol = 1 + mS, + dimnames = + list( + 0:mE, + 0:mS + ) + ) + for (i in seq_along(outcomesY)) { + for (j in seq_along(outcomesZ)) { + b[i, j] <- + postprobDist( + x = outcomesY[i], + n = Nmax, + xS = outcomesZ[j], + nS = NmaxControl, + delta = delta, + relativeDelta = relativeDelta, + parE = parE, + weights = weights, + parS = parS, + weightsS = weightsS + ) + pyz[i, j] <- density_y[i] * pz[j] + } + } + pyz + ret <- list( + result = sum(pyz * (b > thetaT)), + table = data.frame( + counts = c(0:mE), + cumul_counts = n + (0:mE), + density = pyz, + posterior = b, + success = (b > thetaT) + ) + ) + ret +} +``` +# Testing helper function if NmaxControl != 0 +```{r} +h_getlistpredprob( + NmaxControl = 20, + nS = 10, + xS = 5, + weightsS = 1, + par = c(1, 1), + x = 10, + mE = 17, + density_y = density_y +) +``` +# Main function ```{r cars} predprobDist <- function(x, n, xS = 0, nS = 0, @@ -135,15 +203,14 @@ predprobDist <- function(x, n, weightsS <- rep(1, nrow(parS)) } activeBetamixPost <- h_getBetamixPost(x = x, n = n, par = parE, weights = weights) - py <- with( + density_y <- with( activeBetamixPost, dbetabinomMix(x = 0:mE, m = mE, par = par, weights = weights) ) if (NmaxControl == 0) { - h_get_table_singlearm( + ret <- h_get_predproblist_nocontrol( x = x, n = Nmax, - mE = mE, delta = delta, relativeDelta = relativeDelta, parE = parE, @@ -151,88 +218,21 @@ predprobDist <- function(x, n, parS = parS, weightsS = weightsS, thetaT = thetaT, - b = b, - py = py, - Nmax) - # b <- postprobDist( - # x = x + c(0:mE), - # n = Nmax, - # delta = delta, - # relativeDelta = relativeDelta, - # parE = parE, - # weights = weights, - # parS = parS, - # weightsS = weightsS - # ) - # - # ret <- structure(sum(py * (b > thetaT)), - # tables = - # round( - # cbind( - # i = c(0:mE), - # py, - # b, - # bgttheta = (b > thetaT) - # ), - # 4 - # ) - # ) - } else { - mS <- NmaxControl - nS - controlBetamixPost <- h_getBetamixPost( - x = xS, n = nS, par = parS, - weights = weightsS + density = density_y, + mE = mE ) - pz <- with( - controlBetamixPost, - dbetabinomMix(x = 0:mS, m = mS, par = par, weights = weights) - ) - ## determine resulting posterior probabilities: - outcomesY <- x + c(0:mE) - outcomesZ <- xS + c(0:mS) - pyz <- b <- matrix( - nrow = 1 + mE, - ncol = 1 + mS, - dimnames = - list( - 0:mE, - 0:mS - ) - ) - for (i in seq_along(outcomesY)) { - for (j in seq_along(outcomesZ)) { - b[i, j] <- - postprobDist( - x = outcomesY[i], - n = Nmax, - xS = outcomesZ[j], - nS = NmaxControl, - delta = delta, - relativeDelta = relativeDelta, - parE = parE, - weights = weights, - parS = parS, - weightsS = weightsS - ) - pyz[i, j] <- py[i] * pz[j] - } - } - ret <- structure(sum(pyz * (b > thetaT)), - tables = - list( - pyz = pyz, - b = b, - bgttheta = (b > thetaT) - ) + } else { + ret <- h_getlistpredprob( + NmaxControl = nS, + nS = nS, xS = xS, par = parS, weightsS = weightsS, + mS = mS, x = x, mE = mE, density_y = density_y ) } ret } - ``` - - +# An example ```{r pressure, echo=FALSE} predprobDist( x = 16, n = 23, Nmax = 40, delta = 0.1, thetaT = 0.9, diff --git a/tests/testthat/test-predprobDist.R b/tests/testthat/test-predprobDist.R index 8776c03a..066f1570 100644 --- a/tests/testthat/test-predprobDist.R +++ b/tests/testthat/test-predprobDist.R @@ -16,15 +16,48 @@ parE <- t(parE) weights <- rep(1, nrow(parE)) parS <- t(parS) weightsS <- rep(1, nrow(parS)) - -h_get_predproblist_control(x = x, - n = Nmax, - delta = delta, - relativeDelta = relativeDelta, - parE = parE, - weights = weights, - parS = parS, - weightsS = weightsS, - thetaT = thetaT, - density = 0.77, - mE = 17) +h_get_predproblist_control <- function(x, + n, + delta, + relativeDelta, + parE, + weights, + parS, + weightsS, + thetaT, + density, + mE) { + posterior <- postprobDist( + x = x + c(0:mE), + n = Nmax, + delta = delta, + relativeDelta = relativeDelta, + parE = parE, + weights = weights, + parS = parS, + weightsS = weightsS + ) + list( + result = sum(density * (posterior > thetaT)), + table = data.frame( + counts = c(0:mE), + cumul_counts = n + (0:mE), + density = round(density, 4), + posterior = posterior, + success = (posterior > thetaT) + ) + ) +} +h_get_predproblist_control( + x = x, + n = Nmax, + delta = delta, + relativeDelta = relativeDelta, + parE = parE, + weights = weights, + parS = parS, + weightsS = weightsS, + thetaT = thetaT, + density = 0.77, + mE = 17 +) From dc334bb1782facef0dff635edd5c3ec617de32a6 Mon Sep 17 00:00:00 2001 From: Audrey Yeo Date: Wed, 6 Dec 2023 11:38:25 +0100 Subject: [PATCH 07/15] clean --- design/design-doc_predprobDist.Rmd | 3 --- 1 file changed, 3 deletions(-) diff --git a/design/design-doc_predprobDist.Rmd b/design/design-doc_predprobDist.Rmd index 0d7af2f5..6ac420a8 100644 --- a/design/design-doc_predprobDist.Rmd +++ b/design/design-doc_predprobDist.Rmd @@ -96,9 +96,6 @@ h_get_predproblist_nocontrol( ``` # helper function 1 if NmaxControl != 0 ```{r} -activeBetamixPost <- h_getBetamixPost(x = 16, n = 23, par = t(c(1, 1)), weights = 1) -density_y <- with(activeBetamixPost, dbetabinomMix(x = 0:17, m = 17, par = t(c(1, 1)), weights = 1)) - h_getlistpredprob <- function(NmaxControl, nS, xS, par, weightsS, mS, x, mE, density_y) { mS <- NmaxControl - nS controlBetamixPost <- h_getBetamixPost( From 7ea31819fcee477abde64664d932440395af6ab5 Mon Sep 17 00:00:00 2001 From: Audrey Yeo Date: Wed, 6 Dec 2023 11:48:19 +0100 Subject: [PATCH 08/15] rename --- design/design-doc_predprobDist.Rmd | 21 ++++++++++----------- 1 file changed, 10 insertions(+), 11 deletions(-) diff --git a/design/design-doc_predprobDist.Rmd b/design/design-doc_predprobDist.Rmd index 6ac420a8..32597365 100644 --- a/design/design-doc_predprobDist.Rmd +++ b/design/design-doc_predprobDist.Rmd @@ -69,7 +69,7 @@ h_get_predproblist_nocontrol <- function(x, table = data.frame( counts = c(0:mE), cumul_counts = n + (0:mE), - density = round(density, 4), + density = density_y, posterior = posterior_x, success = (posterior_x > thetaT) ) @@ -90,7 +90,7 @@ h_get_predproblist_nocontrol( parS = parS, weightsS = weightsS, thetaT = thetaT, - density = 0.77, + density = density_y, mE = 17 ) ``` @@ -104,14 +104,14 @@ h_getlistpredprob <- function(NmaxControl, nS, xS, par, weightsS, mS, x, mE, den par = parS, weights = weightsS ) - pz <- with( + density_z <- with( controlBetamixPost, dbetabinomMix(x = 0:mS, m = mS, par = par, weights = weights) ) ## determine resulting posterior probabilities: outcomesY <- x + c(0:mE) outcomesZ <- xS + c(0:mS) - pyz <- b <- matrix( + density_yz <- posterior_yz <- matrix( nrow = 1 + mE, ncol = 1 + mS, dimnames = @@ -122,7 +122,7 @@ h_getlistpredprob <- function(NmaxControl, nS, xS, par, weightsS, mS, x, mE, den ) for (i in seq_along(outcomesY)) { for (j in seq_along(outcomesZ)) { - b[i, j] <- + posterior_yz[i, j] <- postprobDist( x = outcomesY[i], n = Nmax, @@ -135,18 +135,17 @@ h_getlistpredprob <- function(NmaxControl, nS, xS, par, weightsS, mS, x, mE, den parS = parS, weightsS = weightsS ) - pyz[i, j] <- density_y[i] * pz[j] + density_yz[i, j] <- density_y[i] * density_z[j] } } - pyz ret <- list( - result = sum(pyz * (b > thetaT)), + result = sum(density_yz * (posterior_yz > thetaT)), table = data.frame( counts = c(0:mE), cumul_counts = n + (0:mE), - density = pyz, - posterior = b, - success = (b > thetaT) + density = density_yz, + posterior = posterior_yz, + success = (posterior_yz > thetaT) ) ) ret From 77da2a1c344f50080ff7e0ac09905f6e593879f6 Mon Sep 17 00:00:00 2001 From: Audrey Yeo Date: Wed, 6 Dec 2023 11:49:36 +0100 Subject: [PATCH 09/15] clean --- design/design-doc_predprobDist.Rmd | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/design/design-doc_predprobDist.Rmd b/design/design-doc_predprobDist.Rmd index 32597365..8ff6a910 100644 --- a/design/design-doc_predprobDist.Rmd +++ b/design/design-doc_predprobDist.Rmd @@ -33,15 +33,15 @@ activeBetamixPost <- h_getBetamixPost(x = x, n = n, par = parE, weights = weight density_y <- with(activeBetamixPost, dbetabinomMix(x = 0:mE, m = mE, par = par, weights = weights)) ``` +# Sanity Check ```{r} predprobDist( x = 16, n = 23, Nmax = 40, delta = 0, thetaT = 0.9, parE = c(0.6, 0.4), parS = c(7, 11) ) ``` -## R Markdown -# helper function if NmaxControl == 0, or when there is no control +# Helper function if NmaxControl == 0, or when there is no control ```{r} h_get_predproblist_nocontrol <- function(x, n, @@ -94,7 +94,8 @@ h_get_predproblist_nocontrol( mE = 17 ) ``` -# helper function 1 if NmaxControl != 0 + +# Helper function 1 if NmaxControl != 0 ```{r} h_getlistpredprob <- function(NmaxControl, nS, xS, par, weightsS, mS, x, mE, density_y) { mS <- NmaxControl - nS @@ -151,6 +152,7 @@ h_getlistpredprob <- function(NmaxControl, nS, xS, par, weightsS, mS, x, mE, den ret } ``` + # Testing helper function if NmaxControl != 0 ```{r} h_getlistpredprob( @@ -164,6 +166,7 @@ h_getlistpredprob( density_y = density_y ) ``` + # Main function ```{r cars} predprobDist <- function(x, n, From 59108c4e07e4e3d48e67368f65c659bbb4bf5211 Mon Sep 17 00:00:00 2001 From: Audrey Yeo Date: Wed, 6 Dec 2023 15:27:20 +0100 Subject: [PATCH 10/15] cumul counts now make sense --- design/design-doc_predprobDist.Rmd | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/design/design-doc_predprobDist.Rmd b/design/design-doc_predprobDist.Rmd index 8ff6a910..366ccfa2 100644 --- a/design/design-doc_predprobDist.Rmd +++ b/design/design-doc_predprobDist.Rmd @@ -68,7 +68,7 @@ h_get_predproblist_nocontrol <- function(x, result = sum(density * (posterior_x > thetaT)), table = data.frame( counts = c(0:mE), - cumul_counts = n + (0:mE), + cumul_counts = x + (0:mE), density = density_y, posterior = posterior_x, success = (posterior_x > thetaT) @@ -143,7 +143,7 @@ h_getlistpredprob <- function(NmaxControl, nS, xS, par, weightsS, mS, x, mE, den result = sum(density_yz * (posterior_yz > thetaT)), table = data.frame( counts = c(0:mE), - cumul_counts = n + (0:mE), + cumul_counts = x + (0:mE), density = density_yz, posterior = posterior_yz, success = (posterior_yz > thetaT) From 9da7b5d556c75954227255a73ab2824ece8970fa Mon Sep 17 00:00:00 2001 From: Audrey Yeo Date: Wed, 6 Dec 2023 17:33:03 +0100 Subject: [PATCH 11/15] Update design/design-doc_predprobDist.Rmd Co-authored-by: Daniel Sabanes Bove --- design/design-doc_predprobDist.Rmd | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/design/design-doc_predprobDist.Rmd b/design/design-doc_predprobDist.Rmd index 366ccfa2..9bb21c54 100644 --- a/design/design-doc_predprobDist.Rmd +++ b/design/design-doc_predprobDist.Rmd @@ -43,7 +43,7 @@ predprobDist( # Helper function if NmaxControl == 0, or when there is no control ```{r} -h_get_predproblist_nocontrol <- function(x, +h_predprobdist_single_arm <- function(x, n, delta, relativeDelta, From 85a04114cc539fb713d926296654e0ae8c990ccb Mon Sep 17 00:00:00 2001 From: Audrey Yeo Date: Wed, 6 Dec 2023 17:50:50 +0100 Subject: [PATCH 12/15] removed test file to external folder --- tests/testthat/test-predprobDist.R | 63 ------------------------------ 1 file changed, 63 deletions(-) delete mode 100644 tests/testthat/test-predprobDist.R diff --git a/tests/testthat/test-predprobDist.R b/tests/testthat/test-predprobDist.R deleted file mode 100644 index 066f1570..00000000 --- a/tests/testthat/test-predprobDist.R +++ /dev/null @@ -1,63 +0,0 @@ -x <- 16 -n <- 23 -xS <- 5 -nS <- 10 -Nmax <- 40 -NmaxControl <- 20 -delta <- 0 -thetaT <- 0.6 -parE <- c(1, 1) -parS <- c(1, 1) -weights <- c(2, 1) -weightsS <- c(2, 1) -relativeDelta <- FALSE -mE <- Nmax - n -parE <- t(parE) -weights <- rep(1, nrow(parE)) -parS <- t(parS) -weightsS <- rep(1, nrow(parS)) -h_get_predproblist_control <- function(x, - n, - delta, - relativeDelta, - parE, - weights, - parS, - weightsS, - thetaT, - density, - mE) { - posterior <- postprobDist( - x = x + c(0:mE), - n = Nmax, - delta = delta, - relativeDelta = relativeDelta, - parE = parE, - weights = weights, - parS = parS, - weightsS = weightsS - ) - list( - result = sum(density * (posterior > thetaT)), - table = data.frame( - counts = c(0:mE), - cumul_counts = n + (0:mE), - density = round(density, 4), - posterior = posterior, - success = (posterior > thetaT) - ) - ) -} -h_get_predproblist_control( - x = x, - n = Nmax, - delta = delta, - relativeDelta = relativeDelta, - parE = parE, - weights = weights, - parS = parS, - weightsS = weightsS, - thetaT = thetaT, - density = 0.77, - mE = 17 -) From f3611ca0606609bc18e20196c4e297312b74b00f Mon Sep 17 00:00:00 2001 From: Audrey Yeo Date: Wed, 6 Dec 2023 18:00:00 +0100 Subject: [PATCH 13/15] Update inst/WORDLIST Co-authored-by: Daniel Sabanes Bove --- inst/WORDLIST | 1 - 1 file changed, 1 deletion(-) diff --git a/inst/WORDLIST b/inst/WORDLIST index 2cd5aa1f..b5d955a2 100644 --- a/inst/WORDLIST +++ b/inst/WORDLIST @@ -91,7 +91,6 @@ increasement inefficacious Integrand integrations -interm LciU ldots leq From d78a130f49fff9d4a203b297b37dca3827ddb45e Mon Sep 17 00:00:00 2001 From: Audrey Yeo Date: Wed, 6 Dec 2023 18:00:07 +0100 Subject: [PATCH 14/15] Update R/predprobDist.R Co-authored-by: Daniel Sabanes Bove --- R/predprobDist.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/predprobDist.R b/R/predprobDist.R index 6daab5cb..3576bcfb 100644 --- a/R/predprobDist.R +++ b/R/predprobDist.R @@ -62,7 +62,7 @@ NULL #' the posterior is `Pr(P_E > P_S + (1 - P_S) * delta | data)`. #' #' @typed x : number -#' number of successes in the treatment group at interm +#' number of successes in the treatment group at interim #' @typed n : number #' number of patients in the treatment group at interim #' @typed xS : number From 26e776907a176f8e2ba476b2dd171180bee76e1d Mon Sep 17 00:00:00 2001 From: Audrey Yeo Date: Wed, 6 Dec 2023 18:03:31 +0100 Subject: [PATCH 15/15] clean --- myEnvironment.RData | Bin 130 -> 0 bytes 1 file changed, 0 insertions(+), 0 deletions(-) delete mode 100644 myEnvironment.RData diff --git a/myEnvironment.RData b/myEnvironment.RData deleted file mode 100644 index 167b2c222cae55323fbd0aea84d2f77a527cf746..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 130 zcmb2|=3oE=X6~X+gJ)e2k`fXU(h?HWQWDZwjU*$So$r+BN=QV^XpnJZ6Jzt6(A0E9 ziH%LoPGRy<$96^@85t>sElgA1^F3Ug?R)YJYtSBtGc7Iu7iuQ!`=9kWuc66i7s9Y34Oru;p&Z_$xmOZ60V85a2+_S*=w8UP5hFsJ|k