Add vignettes, test functions, and indication-specific ORR-PFS updates

.Rbuildignore

@@ -21,5 +21,4 @@
 ^library$
 ^data-raw$
 ^\.merck\.yaml$
-^\.github$
-^src/stan/\*\.rds$
+^\.github$

.gitignore

@@ -12,4 +12,5 @@ Meta
 inst/library/
 library
 src/stan/*.rds
-docs/
+docs/
+.DS_Store

DESCRIPTION

@@ -3,14 +3,15 @@ Type: Package
 Title: Probability of Success (PoS) Calculation for Oncology Pivotal Studies
 Version: 0.1.0
 Authors@R: c(
-    person("Yulia", "Sidi", email = "yulia.sidi@msd.com", role = c("aut", "cre")),
-    person("Xiang", "Peng", email = "xiang.peng@msd.com", role = c("aut")),
+    person("Yulia", "Sidi", email = "yulia.sidi@merck.com", role = c("aut", "cre")),
+    person("Xiang", "Peng", email = "xiang.peng@merck.com", role = c("aut")),
+    person("Shun", "Rao", email = "srao7@mdanderson.org", role = c("aut")),
     person("Alex", "Ziyu Jiang", email = "jiang14@uw.edu", role = c("aut")),
     person("BARDS, MSD", role = c("cph")) )
 Description: Performs Bayesian PoS calculations for phase 3 oncology studies.
-Depends: R (>= 3.6.0)
+Depends: R (>= 4.1.0)
 License: file LICENSE
-URL: https://github.com/MSDLLCPapers/onsoPoS, https://msdllcpapers.github.io/oncoPoS/
+URL: https://github.com/MSDLLCPapers/oncoPoS, https://msdllcpapers.github.io/oncoPoS/
 Encoding: UTF-8
 LazyData: true
 Imports:
@@ -34,4 +35,4 @@ Suggests:
 VignetteBuilder: knitr
 Config/testthat/edition: 3
 Roxygen: list(markdown = TRUE)
-RoxygenNote: 7.3.2
+Config/roxygen2/version: 8.0.0

NAMESPACE

@@ -1,8 +1,11 @@
 # Generated by roxygen2: do not edit by hand
 
+export(compile_stan_models)
 export(gen_pos)
 export(gen_sigma)
+export(get_beta_params)
 export(resp2oddsratio)
+export(resp2oddsratio_estimate_ctrl)
 export(run_stan)
 importFrom(bayesplot,mcmc_acf)
 importFrom(bayesplot,mcmc_trace)
@@ -12,7 +15,10 @@ importFrom(dplyr,summarise)
 importFrom(pracma,erfinv)
 importFrom(purrr,map2_dbl)
 importFrom(purrr,set_names)
+importFrom(rstan,sampling)
 importFrom(rstan,stan)
+importFrom(rstan,stan_model)
 importFrom(stats,qnorm)
+importFrom(stats,sd)
+importFrom(stats,var)
 importFrom(tidybayes,spread_draws)
-importFrom(utils,globalVariables)

NEWS.md

@@ -1,4 +1,4 @@
-# oncoPoS 0.1.0
+# mkstep2pos 0.1.0
 
 - Initial version.
 - Added a `NEWS.md` file to track changes to the package.

R/gen_pos.R

@@ -1,5 +1,12 @@
 #' @title PoS estimation
-#' @description PoS estimation
+#' @description
+#' Computes the predicted probability of success (PoS) for a phase III clinical 
+#' trial by integrating early-phase efficacy data -- objective response rate (ORR)
+#' and/or progression-free survival (PFS) -- with prior beliefs about study 
+#' success. The prior distribution for the benchmark PoS is specified as a Beta 
+#' distribution, parameterized by its mean and variance. Treatment effect 
+#' estimation for ORR supports both two-arm and single-arm designs, with the 
+#' latter incorporating historical control information.
 #' 
 #' @section Specification:
 #' \if{latex}{
@@ -9,6 +16,7 @@
 #'   \item Transfer the hazard ratio bound to log-scale
 #'   \item Extract simulated phase 3 trials treatment effects estimates and
 #'   derive their standard errors
+#'   \item Calculates posterior mean and variance of \code{omega}.
 #'   \item If `plots_out` is TRUE, generate plots for MCMC chains and 
 #'   autocorrelation
 #'   }
@@ -20,42 +28,56 @@
 #' @param plots_out whether plots for MCMC chains and autocorrelation should be 
 #' printed, Default: FALSE
 #' @inheritParams run_stan
-#' @return tibble of PoS estimates and the corresponding standard errors for
-#' each analysis. If `plots_out` is turned on, then the MCMC chains mixing and
+#' @return tibble of PoS estimates, the corresponding standard errors for
+#' each analysis, and the posterior mean and variance for omega. 
+#' If `plots_out` is turned on, then the MCMC chains mixing and
 #' autocorrelation plots are provided as well.
 #' @examples 
-#' # use PFS data from a prior study
+#' # Using both ORR and PFS from a prior single-arm study with a Beta prior on 
+#' # omega
 #' 
-#'   gen_pos(
-#'     target_hr = 0.7, 
-#'     J = 2,  
-#'     nevents3 = c(370, 468), 
-#'     hr_bound = c(0.779, 0.8204), 
-#'     omega = 0.5, 
-#'     est_obs_pfs = 0.88, 
-#'     low_obs_pfs = 0.74, 
-#'     upp_obs_pfs = 1.05, 
-#'     use_pfs = TRUE, 
-#'     seed = 325,
-#'     ncores = 1,
-#'     nchains = 1) 
+#' gen_pos(
+#'   target_hr = 0.70,
+#'   J = 2,
+#'   nevents3 = c(370, 468),
+#'   hr_bound = c(0.7790, 0.8204),
+#'   thres = 0.01,
+#'   omega_mean = 0.3,
+#'   omega_var = 0.03,
+#'   est_obs_pfs = 0.73,
+#'   low_obs_pfs = 0.61,
+#'   upp_obs_pfs = 0.91,
+#'   use_pfs = TRUE,
+#'   n_trt2 = 100,
+#'   n_resp_trt2 = 40,
+#'   low_soc_rr = 0.05,
+#'   upp_soc_rr = 0.2,
+#'   use_orr = TRUE,
+#'   single_arm = TRUE,
+#'   ncores = 1,
+#'   nchains = 1,
+#'   seed = 222
+#' )
 #' @seealso 
 #'  \code{\link[tidybayes]{gather_draws}}
-#'  \code{\link[dplyr]{mutate}}, \code{\link[dplyr]{group_by}}, \code{\link[dplyr]{summarise}}
+#'  \code{\link[dplyr]{mutate}}, \code{\link[dplyr]{summarise}},
+#'  \code{\link[dplyr]{group_by}}
 #'  \code{\link[purrr]{map2}}, \code{\link[purrr]{reexports}}
-#'  \code{\link[bayesplot]{MCMC-traces}}, \code{\link[bayesplot]{MCMC-diagnostics}}
+#'  \code{\link[bayesplot]{MCMC-traces}}, 
+#'  \code{\link[bayesplot]{MCMC-diagnostics}}
 #' @rdname gen_pos
 #' @export 
 #' @importFrom tidybayes spread_draws
-#' @importFrom dplyr mutate group_by summarise
+#' @importFrom dplyr mutate summarise group_by
 #' @importFrom purrr map2_dbl set_names
 #' @importFrom bayesplot mcmc_trace mcmc_acf
+#' @importFrom stats var
 gen_pos <- function(
     target_hr,
     J,
     nevents3,
     hr_bound, 
-    omega,
+    omega_mean = 0.52, omega_var = 0.02,
     est_obs_pfs,
     low_obs_pfs,
     upp_obs_pfs,
@@ -65,30 +87,35 @@ gen_pos <- function(
     n_ctrl2,
     n_resp_trt2, 
     n_resp_ctrl2,
+    low_soc_rr,
+    upp_soc_rr,
+    ci_rr = 0.8,
     use_orr = FALSE,
+    single_arm = FALSE,
     use_pfs = FALSE,
     het_degree_p2 = "small",
     het_degree_p3 = "very small",
     ratio = 1, # ratio of experimental over control
+    indication = 6,
     m_0 = NA,
     m_1 = NA,
     nu_0 = NA,
     nu_1 = NA,
     lm_sd = NA,
     niter = 1000,
-    nchains = 4,
-    ncores = 4,
+    nchains = 1,
+    ncores = 1,
     seed,
     plots_out = FALSE,
     ...
-    ) {
-
+) {
+  
   run_stan_res <- run_stan(
     target_hr = target_hr,
     J = J,
     nevents3 = nevents3,
     hr_bound = hr_bound, 
-    omega = omega,
+    omega_mean = omega_mean, omega_var = omega_var,   
     est_obs_pfs = est_obs_pfs,
     low_obs_pfs = low_obs_pfs,
     upp_obs_pfs = upp_obs_pfs,
@@ -98,11 +125,13 @@ gen_pos <- function(
     n_ctrl2 = n_ctrl2,
     n_resp_trt2 = n_resp_trt2, 
     n_resp_ctrl2 = n_resp_ctrl2,
-    use_orr = use_orr,
+    low_soc_rr = low_soc_rr, upp_soc_rr = upp_soc_rr, ci_rr = ci_rr,
+    use_orr = use_orr, single_arm = single_arm,
     use_pfs = use_pfs,
     het_degree_p2 = het_degree_p2,
     het_degree_p3 = het_degree_p3,
     ratio = ratio, 
+    indication = indication,
     m_0 = m_0,
     m_1 = m_1,
     nu_0 = nu_0,
@@ -113,41 +142,50 @@ gen_pos <- function(
     ncores = ncores,
     seed = seed,
     ...
-    )
+  )
 
   fit_rstan <- run_stan_res$fit_rstan
 
   loghr_bound <- log(hr_bound) # transform hr bound to log scale
 
+  omega_summary <- fit_rstan |>
+    tidybayes::spread_draws(omega) |>
+    dplyr::summarise(
+      omega_mean = mean(omega),
+      omega_var = stats::var(omega)
+      )
+
   # extract simulated treatment effects for phase III experiments
-  out <- tidybayes::spread_draws(model = fit_rstan, theta_P3_hat[J]) |>
+  out <- fit_rstan |>
+    tidybayes::spread_draws(theta_P3_hat[J]) |>
     dplyr::mutate(
-      reject = purrr::map2_dbl(J, theta_P3_hat,  function(x, y){
+      reject = purrr::map2_dbl(J, theta_P3_hat, function(x, y) {
         as.numeric(y < loghr_bound[[x]])
       })
     ) |>
     dplyr::group_by(J) |>
     dplyr::summarise(pos = mean(reject)) |>
     dplyr::mutate(
-      pos_se = sqrt(pos * (1 - pos)/(nchains*niter/2))
+      pos_se = sqrt(pos * (1 - pos)/(nchains * niter / 2)),
+      omega_mean = omega_summary$omega_mean,
+      omega_var = omega_summary$omega_var
     )
 
   if (plots_out) {
     suppressWarnings({
       p1 <- bayesplot::mcmc_trace(
         fit_rstan,
-        pars = c("mu_P", "tau_P2", "tau_P3", "theta_P2", "theta_P3")
+        pars = c("mu_P", "tau_P2", "tau_P3", "theta_P2", "theta_P3", "omega")
       )
       p2 <- bayesplot::mcmc_acf(
         fit_rstan,
-        pars = c("mu_P", "tau_P2", "tau_P3", "theta_P2", "theta_P3")
+        pars = c("mu_P", "tau_P2", "tau_P3", "theta_P2", "theta_P3", "omega")
       )      
     })
-
+    
     out <- list(p1, p2, out) |> purrr::set_names(c("trace", "auto", "pos_est"))
   }
   return(out)
 }
 
-#' @importFrom utils globalVariables
-utils::globalVariables(c('pos', 'reject', 'theta_P3_hat'))
+

R/get_beta_params.R

@@ -0,0 +1,59 @@
+#' @title Generate Beta Distribution Parameters for Omega
+#' @description Computes the shape parameters of a Beta distribution based on 
+#' the specified mean and variance. This is used to parameterize a Beta prior 
+#' distribution for the benchmark probability for study success, \code{omega} 
+#' in Bayesian modeling of phase III trial success.
+#' 
+#' @section Specification:
+#' \if{latex}{
+#'   \itemize{
+#'     \item Accepts prior mean and variance for a Beta-distributed quantity.
+#'     \item These can be obtained from predictive modeling (e.g., Benchmark
+#'     probability model using random forest) or empirical prior estimates.
+#'     \item Converts the mean and variance into Beta shape parameters 
+#'     \eqn{\alpha} and \eqn{\beta}.
+#'     \item Returns the resulting shape parameters as a named list.
+#'   }
+#' }
+#' \if{html}{
+#'   The contents of this section are shown in the PDF user manual only.
+#' }
+#'
+#' @param mean Numeric value between 0 and 1, specifying the prior mean of the 
+#' Beta distribution. This may be estimated from a Benchmark probability model 
+#' such as random forest.
+#' @param var Numeric value, specifying the prior variance of the Beta 
+#' distribution. This may reflect uncertainty in the model-based prediction.
+#'            
+#' @return A named list with elements:
+#' \describe{
+#'   \item{\code{alpha}}{First shape parameter of the Beta distribution.}
+#'   \item{\code{beta}}{Second shape parameter of the Beta distribution.}
+#' }
+#' 
+#' @details
+#' The Beta distribution is parameterized by two positive shape parameters, 
+#' \eqn{\alpha} and \eqn{\beta}, 
+#' which can be derived from a given mean \eqn{\mu} and variance \eqn{\sigma^2} 
+#' using:
+#' \deqn{
+#'   \alpha = \mu \left( \frac{\mu(1 - \mu)}{\sigma^2} - 1 \right), \quad
+#'   \beta = (1 - \mu) \left( \frac{\mu(1 - \mu)}{\sigma^2} - 1 \right)
+#' }
+#' These parameters allow for a flexible specification of prior distributions, 
+#' and are particularly useful when the prior belief is derived from a 
+#' predictive model in an earlier step (e.g., machine learning model 
+#' estimating historical success probabilities).
+#'
+#' @examples
+#' # Example using prior mean and variance estimated from a Benchmark 
+#' # probability model
+#' get_beta_params(mean = 0.52, var = 0.02)
+#'
+#' @rdname get_beta_params
+#' @export
+get_beta_params <- function(mean, var) {
+  alpha <- mean * ((mean * (1 - mean)) / var - 1)
+  beta <- (1 - mean) * ((mean * (1 - mean)) / var - 1)
+  list(alpha = alpha, beta = beta)
+}

R/globals.R

@@ -0,0 +1 @@
+utils::globalVariables(c("omega", "pos", "reject", "theta_P3_hat"))