109 my plot diff

.gitignore

@@ -30,3 +30,4 @@ coverage.*
 .vscode/
 .rds
 node_modules
+vintage_NAMESPACE

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)

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

R/plotBeta.R

@@ -2,19 +2,19 @@
 #'
 #' 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
 #'
 #' @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,91 @@ 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
 }

examples/plotBeta.R

@@ -1,2 +1,3 @@
+# plotBeta
 plotBeta(alpha = 4, beta = 5)
 plotBeta(alpha = 1, beta = 1)

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
+)

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