diff --git a/.gitignore b/.gitignore index 8620c310..52d77690 100644 --- a/.gitignore +++ b/.gitignore @@ -30,3 +30,4 @@ coverage.* .vscode/ .rds node_modules +vintage_NAMESPACE diff --git a/NAMESPACE b/NAMESPACE index 5ad2b34d..1b06e70d 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -6,7 +6,6 @@ export(dbetaMix) export(dbetabinom) export(dbetabinomMix) export(dbetadiff) -export(myPlotDiff) export(oc2) export(oc3) export(ocPostprob) @@ -18,6 +17,7 @@ export(ocRctPredprobDist) export(pbetaMix) export(pbetadiff) export(plotBeta) +export(plotBetaDiff) export(plotBounds) export(plotDecision) export(plotOc) diff --git a/R/boundsPostprob.R b/R/boundsPostprob.R index a35fef8d..9cf86fe5 100644 --- a/R/boundsPostprob.R +++ b/R/boundsPostprob.R @@ -2,7 +2,7 @@ #' #' This function is used to identify the efficacy and futility #' boundaries based on the following rules: -#' Efficacy boundary: find minimum x (xU) where Pr(RR > p1 |x, n, a, b) >= tU and +#' Efficacy boundary: find minimum x (xU) where Pr(RR > p1 | x, n, a, b) >= tU and #' Futility boundary: find maximum x (xL) where Pr(RR < p0 | x, n, a, b) >= tL #' #' @inheritParams postprob @@ -48,12 +48,12 @@ boundsPostprob <- function(looks, p0, p1 = p0, tL, tU, parE = c(1, 1), weights) xL <- NA xU <- NA for (x in 0:n) { - postp_fut <- 1 - postprob(x, n, p0, parE, weights) # futility look + postp_fut <- 1 - postprob(x = x, n = n, p = p0, parE = parE, weights = weights) # futility look if (postp_fut >= tL) { # Rule is P(RR < p0) > tL postL <- postp_fut xL <- x } - postp_eff <- postprob(x, n, p1, parE, weights) # efficacy look + postp_eff <- postprob(x = x, n = n, p = p1, parE = parE, weights = weights) # efficacy look if (postp_eff >= tU) { # Rule is P(RR > p1) > tU postU <- postp_eff xU <- x diff --git a/R/plotBeta.R b/R/plotBeta.R index 98a6788c..d4d71a77 100644 --- a/R/plotBeta.R +++ b/R/plotBeta.R @@ -2,11 +2,11 @@ #' #' This function will plot the PDF of a beta distribution #' +#' @inheritParams dbetabinom #' @typed alpha : number #' first parameter of the Beta distribution #' @typed beta : number #' second parameter of the Beta distribution -#' #' @return A beta distribution density plot #' #' @importFrom graphics axis @@ -14,7 +14,7 @@ #' @example examples/plotBeta.R #' @export #' @keywords graphics -plotBeta <- function(alpha, beta) { +plotBeta <- function(alpha, beta, ...) { x_support <- seq(from = 0, to = 1, length = 1000) data <- data.frame( grid = x_support, @@ -30,141 +30,87 @@ plotBeta <- function(alpha, beta) { ggplot2::scale_x_continuous(labels = scales::percent_format()) } -#' Plot Diff Between two Beta distributions +#' Plot difference Between two Beta distributions #' #' This function will plot the PDF of a difference between two Beta distributions #' -#' @param parY non-negative parameters of the treatment Beta distribution. -#' @param parX non-negative parameters of the historical control Beta distribution -#' @param cut_B a meaningful improvement threshold -#' @param cut_W a poor improvement throshold -#' @param shade paint the two areas under the curve, default value=1 as "yes". other numbers stands for "no"; -#' @param note show values of the colored area, default value=1 as "yes". other numbers stands for "no" -#' @param \dots additional arguments to \code{plot} -#' @return nothing, only produces the plot as side effect +#' @typed parY : numeric +#' non-negative parameters of the treatment Beta distribution. +#' @typed parX : numeric +#' non-negative parameters of the historical control Beta distribution +#' @typed cut_B : number +#' a meaningful improvement threshold, the lower boundary of a meaningfully improvement in response rate +#' @typed cut_W : number +#' a poor improvement threshold, the upper boundary of a meaningfully poor improvement in response rate +#' @typed shade : flag +#' paint the two areas under the curve, default value = TRUE +#' @typed note : flag +#' show values of the colored area, default value = TRUE +#' @typed ... : +#' additional arguments to `ggplot()` +#' @return a ggplot object #' -#' @example examples/myPlotDiff.R +#' @example examples/plotBetaDiff.R #' #' @importFrom graphics par axis polygon mtext #' @importFrom stats integrate #' #' @export #' @keywords graphics -myPlotDiff <- function(parY, # parameters of phase Ib trial; - parX, # parameters of HC; - cut_B = 0.20, # a meaningful improvement threshold; - cut_W = 0.1, # a poor improvement threshold; - shade = 1, # paint the two areas under the curve, default: yes. other numbers stands for "no"; - note = 1, # show values of the colored area, default: yes. other numbers stands for "no"; - ...) { - if (note == 1) { - graphics::par(mar = c(5, 15, 1, 15) + .1) - } else { - graphics::par(mar = c(5, 5, 1, 5) + .1) - } - grid <- seq(from = -0.5, to = 0.75, length = 1000) - xticks <- seq(from = -1, to = 1, by = 0.25) - - - - graphics::plot( - x = grid, - y = dbetadiff(grid, parY = parY, parX = parX), - ylab = "", - xaxt = "n", - yaxt = "n", - type = "l", - xaxs = "i", - yaxs = "i", - ... +plotBetaDiff <- function(parY, # parameters of experimental arm + parX, # parameters of control or SOC + Go_cut = 0.20, # a meaningful improvement threshold + Stop_cut = 0.1, # a poor improvement threshold + shade = TRUE, # paint the two areas under the curve + note = TRUE) { # show values of the colored area + diff <- seq(from = -1, to = 1, length = 1000) + data <- data.frame( + grid = diff, + density = dbetadiff(z = diff, parY = parY, parX = parX) ) + data$Stop <- ifelse(diff > -1 & diff < Stop_cut, TRUE, FALSE) + data$Go <- ifelse(diff > Go_cut & diff < 1, TRUE, FALSE) - graphics::axis( - side = 1, at = xticks, - labels = - paste(ifelse(xticks >= 0, "+", ""), - xticks * 100, "%", - sep = "" - ) + Go_auc <- integrate( + f = dbetadiff, + parY = parY, + parX = parX, + lower = Go_cut, # Calculate probability of Go, if difference was at least `Go_cut`. + upper = 1 + ) + Stop_auc <- integrate( + f = dbetadiff, + parY = parY, + parX = parX, + lower = -1, + upper = Stop_cut # Calculate probability of Stop, if difference was at most `Stop_cut`. ) - ## now color the go / stop prob areas - - if (shade == 1) { - ## first stop: - stopGrid <- grid[grid <= cut_W] - nStop <- length(stopGrid) - - graphics::polygon( - x = - c( - stopGrid, - rev(stopGrid) - ), - y = - c( - rep(0, nStop), - dbetadiff(rev(stopGrid), parY = parY, parX = parX) - ), - col = "red" - ) - - A_value <- stats::integrate( - f = dbetadiff, - parY = parY, - parX = parX, - lower = -1, - upper = cut_W - ) - if (note == 1) { - graphics::mtext( - paste("Prob(diff< ", round(cut_W * 100), "%)=", - sprintf("%1.2f%%", 100 * as.numeric(A_value$value)), - sep = "" - ), - side = 2, line = 1, las = 1, cex = 1 - ) - } - - ## then go: - goGrid <- grid[grid >= cut_B] - nGo <- length(goGrid) - - graphics::polygon( - x = - c( - goGrid, - rev(goGrid) - ), - y = - c( - rep(0, nGo), - dbetadiff(rev(goGrid), parY = parY, parX = parX) - ), - col = "green" - ) - - B_value <- stats::integrate( - f = dbetadiff, - parY = parY, - parX = parX, - lower = cut_B, - upper = 1 - ) + Go_label <- paste("P(Go) is", round(Go_auc$value * 100, digits = 2), "%") + Stop_label <- paste("P(S) is", round(Stop_auc$value * 100, digits = 2), "%") + plot_title <- paste("According to Beta difference density", Go_label, "and", Stop_label) - if (note == 1) { - graphics::mtext( - paste( - sprintf("%1.2f%%", 100 * as.numeric(B_value$value)), - "=Prob(diff> ", - round(cut_B * 100), "%)", - sep = "" - ), - side = 4, - line = 1, - las = 1, - cex = 1 - ) - } + if (shade == TRUE) { + pbetadiff_plot <- ggplot2::ggplot(data = data, aes(x = grid, y = density)) + + ggplot2::geom_line(colour = "#888888") + + geom_area(data = data[data$grid < Stop_cut,], fill = "#FF0046", + mapping = aes(x = ifelse(grid < 0.2 & grid < 0.5, grid, 0))) + + geom_area(data = data[data$grid > Go_cut,], fill = "#009E73", + mapping = aes(x = ifelse(grid > 0.3, grid, 0))) + + xlab("Difference between treatment") + + ggplot2::ylab(quote(f(x))) + + ggplot2::ggtitle(plot_title) + } else { + pbetadiff_plot <- ggplot2::ggplot(data = data, aes(x = grid, y = density)) + + ggplot2::geom_line(colour = "#888888") + + xlab("Difference between treatment") + + ggplot2::ylab(quote(f(x))) + + ggplot2::ggtitle(plot_title) + } + if (note == TRUE) { + pbetadiff_plot <- pbetadiff_plot + + ggplot2::annotate("text", x = -0.5, y = 3.75, size = 5, label = Stop_label, colour = "#FF0046") + + ggplot2::annotate("text", x = -0.5, y = 3.25, size = 5, label = Go_label, colour = "#009E73") } + pbetadiff_plot } diff --git a/R/plotOc.R b/R/plotOc.R index c4765cdd..9806db9a 100644 --- a/R/plotOc.R +++ b/R/plotOc.R @@ -1,9 +1,10 @@ #' Display the operating characteristics using an oc object #' -#' Reads results from \code{\link{ocPredprob}}, \code{\link{ocPostprob}} +#' Reads results from [ocPredprob()] #' etc. and displays a bar plot of the operating characteristics #' -#' @param z returned oc value +#' @typed oc : list +#' returned oc parameters #' @return nothing, only plots as side effect #' #' @importFrom graphics barplot title @@ -11,9 +12,25 @@ #' @example examples/plotOc.R #' @export #' @keywords graphics -plotOc <- function(z) { +plotOc <- function(oc) { + + if (wiggle == FALSE) { + data <- table(oc$Decision, oc$SampleSize) / oc$params$sim + } else { + + } + +ggplot(oc, aes(x = name, y = value)) + + geom_bar(stat = "identity") + + ggtitle("Percentage of trials that Go and Stop per look") + + ylabs("Percentage %") + + xlabs("Looks and sample size") + + + ## plot function for oc.predprob or oc.postprob, or the dist versions of them - graphics::barplot(table(z$Decision, z$SampleSize) / z$params$sim, beside = TRUE) + graphics::barplot(table(oc$Decision, oc$SampleSize) / oc$params$sim, beside = TRUE) + ## get the parameter parDat <- lapply(z$params, deparse) diff --git a/examples/myPlotDiff.R b/examples/myPlotDiff.R deleted file mode 100644 index c9a4795e..00000000 --- a/examples/myPlotDiff.R +++ /dev/null @@ -1,8 +0,0 @@ -myPlotDiff( - parY = c(5, 10), - parX = c(2, 5), - cut_B = 0.2, # a meaningful improvement threshold - cut_W = 0.05, # a poor improvement threshold - shade = 1, # paint the two areas under the curve, default: yes. other numbers stands for "no"; - note = 0 -) # show values of the colored area, default: yes. other numbers stands for "no"; diff --git a/examples/plotBeta.R b/examples/plotBeta.R index 6e84087b..712b5482 100644 --- a/examples/plotBeta.R +++ b/examples/plotBeta.R @@ -1,2 +1,3 @@ +# plotBeta plotBeta(alpha = 4, beta = 5) plotBeta(alpha = 1, beta = 1) diff --git a/examples/plotBetaDiff.R b/examples/plotBetaDiff.R new file mode 100644 index 00000000..277e175e --- /dev/null +++ b/examples/plotBetaDiff.R @@ -0,0 +1,23 @@ +# The beta distribution and acceptable bounds for +# a meaningful improvement of 0.20 and worsening of 0.1 +parX <- c(1, 52) # prior parameters of experimental arm +parY <- c(5.5, 20.5) # prior parameters of control or SOC +plotBetaDiff( + parY = parY, + parX = parX, + Go_cut = 0.3, + Stop_cut = 0.1, # below a difference of 10%, is an unsuccesful trial + shade = TRUE, + note = TRUE +) + + +# a larger Go_cut with uniform prior +plotBetaDiff( + parY = c(1, 1), # prior parameters for experimental arm + parX = c(1, 1), # prior parameters for control or SOC arm + Go_cut = 0.3, + Stop_cut = 0.1, # below a difference of 10%, is an unsuccesful trial + shade = TRUE, + note = TRUE +) diff --git a/examples/plotOc.R b/examples/plotOc.R index 930ca464..53216078 100644 --- a/examples/plotOc.R +++ b/examples/plotOc.R @@ -1,9 +1,16 @@ # get operating character result from oc.postprob -res1 <- ocPostprob( +oc <- ocPostprob( nnE = c(10, 20, 30), truep = 0.4, p0 = 0.2, - p1 = 0.3, tL = 0.6, tU = 0.8, parE = c(1, 1), sim = 50000 + p1 = 0.3, tL = 0.6, tU = 0.8, parE = c(1, 1), sim = 100, wiggle = FALSE ) -res1$oc +oc$oc +# plotOc(oc) + + +oc <- ocPostprob( + nnE = c(10, 20, 30), truep = 0.4, p0 = 0.2, + p1 = 0.3, tL = 0.6, tU = 0.8, parE = c(1, 1), sim = 100, wiggle = TRUE +) +oc$oc -plotOc(res1) diff --git a/inst/WORDLIST b/inst/WORDLIST index 2faf53fe..be3cf716 100644 --- a/inst/WORDLIST +++ b/inst/WORDLIST @@ -102,6 +102,7 @@ gelman Gelman generalizable geq +ggplot grayzone grey Gsponer @@ -201,6 +202,8 @@ renewcommand reproducibility responder responders +roxygen +Roxygen Sabanes sabanes Sabanés @@ -243,6 +246,7 @@ USUBJID VAD vanillaBayes vanillaPP +vbump Vehtari WeightedBayes weightedBetaPrior diff --git a/man/boundsPostprob.Rd b/man/boundsPostprob.Rd index b21240cd..1fa8bde4 100644 --- a/man/boundsPostprob.Rd +++ b/man/boundsPostprob.Rd @@ -40,7 +40,7 @@ binomial test. \description{ This function is used to identify the efficacy and futility boundaries based on the following rules: -Efficacy boundary: find minimum x (xU) where Pr(RR > p1 |x, n, a, b) >= tU and +Efficacy boundary: find minimum x (xU) where Pr(RR > p1 | x, n, a, b) >= tU and Futility boundary: find maximum x (xL) where Pr(RR < p0 | x, n, a, b) >= tL } \examples{ diff --git a/man/myPlotDiff.Rd b/man/myPlotDiff.Rd deleted file mode 100644 index 8d949d11..00000000 --- a/man/myPlotDiff.Rd +++ /dev/null @@ -1,40 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/plotBeta.R -\name{myPlotDiff} -\alias{myPlotDiff} -\title{Plot Diff Between two Beta distributions} -\usage{ -myPlotDiff(parY, parX, cut_B = 0.2, cut_W = 0.1, shade = 1, note = 1, ...) -} -\arguments{ -\item{parY}{non-negative parameters of the treatment Beta distribution.} - -\item{parX}{non-negative parameters of the historical control Beta distribution} - -\item{cut_B}{a meaningful improvement threshold} - -\item{cut_W}{a poor improvement throshold} - -\item{shade}{paint the two areas under the curve, default value=1 as "yes". other numbers stands for "no";} - -\item{note}{show values of the colored area, default value=1 as "yes". other numbers stands for "no"} - -\item{\dots}{additional arguments to \code{plot}} -} -\value{ -nothing, only produces the plot as side effect -} -\description{ -This function will plot the PDF of a difference between two Beta distributions -} -\examples{ -myPlotDiff( - parY = c(5, 10), - parX = c(2, 5), - cut_B = 0.2, # a meaningful improvement threshold - cut_W = 0.05, # a poor improvement threshold - shade = 1, # paint the two areas under the curve, default: yes. other numbers stands for "no"; - note = 0 -) # show values of the colored area, default: yes. other numbers stands for "no"; -} -\keyword{graphics} diff --git a/man/plotBeta.Rd b/man/plotBeta.Rd index bf45b581..a17df369 100644 --- a/man/plotBeta.Rd +++ b/man/plotBeta.Rd @@ -4,7 +4,7 @@ \alias{plotBeta} \title{Plot the Beta distribution} \usage{ -plotBeta(alpha, beta) +plotBeta(alpha, beta, ...) } \arguments{ \item{alpha}{(\code{number}):\cr first parameter of the Beta distribution} @@ -18,6 +18,7 @@ A beta distribution density plot This function will plot the PDF of a beta distribution } \examples{ +# plotBeta plotBeta(alpha = 4, beta = 5) plotBeta(alpha = 1, beta = 1) } diff --git a/man/plotBetaDiff.Rd b/man/plotBetaDiff.Rd new file mode 100644 index 00000000..ac715de5 --- /dev/null +++ b/man/plotBetaDiff.Rd @@ -0,0 +1,62 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/plotBeta.R +\name{plotBetaDiff} +\alias{plotBetaDiff} +\title{Plot difference Between two Beta distributions} +\usage{ +plotBetaDiff( + parY, + parX, + Go_cut = 0.2, + Stop_cut = 0.1, + shade = TRUE, + note = TRUE +) +} +\arguments{ +\item{parY}{(\code{numeric}):\cr non-negative parameters of the treatment Beta distribution.} + +\item{parX}{(\code{numeric}):\cr non-negative parameters of the historical control Beta distribution} + +\item{shade}{(\code{flag}):\cr paint the two areas under the curve, default value = TRUE} + +\item{note}{(\code{flag}):\cr show values of the colored area, default value = TRUE} + +\item{cut_B}{(\code{number}):\cr a meaningful improvement threshold, the lower boundary of a meaningfully improvement in response rate} + +\item{cut_W}{(\code{number}):\cr a poor improvement threshold, the upper boundary of a meaningfully poor improvement in response rate} + +\item{...}{(``):\cr additional arguments to \code{ggplot()}} +} +\value{ +a ggplot object +} +\description{ +This function will plot the PDF of a difference between two Beta distributions +} +\examples{ +# The beta distribution and acceptable bounds for +# a meaningful improvement of 0.20 and worsening of 0.1 +parX <- c(1, 52) # prior parameters of experimental arm +parY <- c(5.5, 20.5) # prior parameters of control or SOC +plotBetaDiff( + parY = parY, + parX = parX, + Go_cut = 0.3, + Stop_cut = 0.1, # below a difference of 10\%, is an unsuccesful trial + shade = TRUE, + note = TRUE +) + + +# a larger Go_cut with uniform prior +plotBetaDiff( + parY = c(1, 1), # prior parameters for experimental arm + parX = c(1, 1), # prior parameters for control or SOC arm + Go_cut = 0.3, + Stop_cut = 0.1, # below a difference of 10\%, is an unsuccesful trial + shade = TRUE, + note = TRUE +) +} +\keyword{graphics} diff --git a/man/plotOc.Rd b/man/plotOc.Rd index 16389213..76965376 100644 --- a/man/plotOc.Rd +++ b/man/plotOc.Rd @@ -4,27 +4,64 @@ \alias{plotOc} \title{Display the operating characteristics using an oc object} \usage{ -plotOc(z) +plotOc(oc) } \arguments{ -\item{z}{returned oc value} +\item{oc}{(\code{list}):\cr returned oc parameters} } \value{ nothing, only plots as side effect } \description{ -Reads results from \code{\link{ocPredprob}}, \code{\link{ocPostprob}} +Reads results from \code{\link[=ocPredprob]{ocPredprob()}} etc. and displays a bar plot of the operating characteristics } \examples{ # get operating character result from oc.postprob -res1 <- ocPostprob( +oc <- ocPostprob( nnE = c(10, 20, 30), truep = 0.4, p0 = 0.2, - p1 = 0.3, tL = 0.6, tU = 0.8, parE = c(1, 1), sim = 50000 + p1 = 0.3, tL = 0.6, tU = 0.8, parE = c(1, 1), sim = 100, wiggle = FALSE ) -res1$oc +oc$oc +# plotOc(oc) -plotOc(res1) + +oc <- ocPostprob( + nnE = c(10, 20, 30), truep = 0.4, p0 = 0.2, + p1 = 0.3, tL = 0.6, tU = 0.8, parE = c(1, 1), sim = 100, wiggle = TRUE +) +oc$oc + +data <- data.frame(Decision = as.factor(oc$Decision), + Looks = oc$SampleSize) + +table(data$Decision, data$Looks) + +data \%>\% group_by(Decision, Looks) \%>\% summarise(n = n()) + +data \%>\% group_by(SampleSize) \%>\% summarise(n = n()) + +ggplot(dat, aes(x = SampleSize, y = Decision)) + + geom_bar(stat = "identity") + +data \%>\% group_by(Looks, Decision) \%>\% summarise(tot = n()) + + +table(oc$Decision, oc$SampleSize) / oc$params$sim +table(oc$Decision, oc$SampleSize) + +oc$Decision <- ifelse(oc$Decision == TRUE, "GO", + ifelse(oc$Decision == FALSE, "Stop", "Grey")) + +oc$Decision + +data.frame(looks = c(res1$wiggled_nnrE)) + +table(oc$Decision, oc$SampleSize) + +tt = data.frame(Decision = res$Decision, + look = res$SampleSize) +tt \%>\% group_by(look) \%>\% summarise(success = per()) } \keyword{graphics} diff --git a/tests/testthat/test-plotBetaDiff.R b/tests/testthat/test-plotBetaDiff.R new file mode 100644 index 00000000..e7c2a091 --- /dev/null +++ b/tests/testthat/test-plotBetaDiff.R @@ -0,0 +1,22 @@ +# plotBetaDiff +test_that("plotBetaDiff works as expected", { + result <- plotBetaDiff( + parY = c(1, 1), + parX = c(6, 10), + Go_cut = 0.3, + Stop_cut = 0.1, # below a difference of 10%, is an unsuccesful trial + shade = TRUE, + note = TRUE + ) + result <- plotBetaDiff( + parY = c(2, 3), + parX = c(4, 10), + Go_cut = 0.3, + Stop_cut = 0.1, # below a difference of 10%, is an unsuccesful trial + shade = TRUE, + note = TRUE + ) + vdiffr::expect_doppelganger("", result) + vdiffr::expect_doppelganger("", result) + +})