diff --git a/.Rbuildignore b/.Rbuildignore index f9a0624..a2b9b19 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -21,5 +21,4 @@ ^library$ ^data-raw$ ^\.merck\.yaml$ -^\.github$ -^src/stan/\*\.rds$ \ No newline at end of file +^\.github$ \ No newline at end of file diff --git a/.gitignore b/.gitignore index fc08723..050b0af 100644 --- a/.gitignore +++ b/.gitignore @@ -12,4 +12,5 @@ Meta inst/library/ library src/stan/*.rds -docs/ \ No newline at end of file +docs/ +.DS_Store \ No newline at end of file diff --git a/DESCRIPTION b/DESCRIPTION index 2027e31..2aafd19 100644 --- a/DESCRIPTION +++ b/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 diff --git a/NAMESPACE b/NAMESPACE index a96c5b7..f46ebff 100644 --- a/NAMESPACE +++ b/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) diff --git a/NEWS.md b/NEWS.md index 3ae9392..0cf11a7 100644 --- a/NEWS.md +++ b/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. diff --git a/R/gen_pos.R b/R/gen_pos.R index 43d84d8..f1eb04d 100644 --- a/R/gen_pos.R +++ b/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')) + diff --git a/R/get_beta_params.R b/R/get_beta_params.R new file mode 100644 index 0000000..29b6f5e --- /dev/null +++ b/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) +} diff --git a/R/globals.R b/R/globals.R new file mode 100644 index 0000000..df46cb4 --- /dev/null +++ b/R/globals.R @@ -0,0 +1 @@ +utils::globalVariables(c("omega", "pos", "reject", "theta_P3_hat")) diff --git a/R/resp2oddsratio_estimate_ctrl.R b/R/resp2oddsratio_estimate_ctrl.R new file mode 100644 index 0000000..322f961 --- /dev/null +++ b/R/resp2oddsratio_estimate_ctrl.R @@ -0,0 +1,136 @@ +#' @title Estimate Log Odds Ratio Using Observed Response Data and Prior on +#' Control Arm +#' +#' @description +#' Estimates the log odds ratio (log OR) between treatment and control arms +#' in a single-arm trial by combining observed responses in treatment arm with +#' prior knowledge of the control response rate. The prior distribution on the +#' control arm's response rate is specified through a credible interval. +#' +#' @section Specification: +#' \if{latex}{ +#' \itemize{ +#' \item Compute the observed response rate in the treatment arm. +#' \item Translate the control response rate bounds into a normal prior on +#' the logit scale. +#' \item Fit a Bayesian model to estimate the posterior distribution of the +#' log OR. +#' \item Return the posterior mean and standard deviation of the log OR. +#' } +#' } +#' \if{html}{ +#' The contents of this section are shown in the PDF user manual only. +#' } +#' +#' @param n_resp_trt number of responses in the treatment arm +#' @param n_trt sample size in the treatment arm +#' @param low_soc_rr lower bound of the control arm's response rate +#' @param upp_soc_rr upper bound of the control arm's response rate +#' @param ci_rr confidence level (e.g., 0.80) for the control arm response rate +#' interval. +#' @param niter number of iterations to be used in stan run, Default: 1000 +#' @param nchains number of chains to be used in stan run, Default: 4 +#' @param ncores number of cores to be used in stan run, Default: 4 +#' @param seed seed to be used in stan run +#' @param refresh integer, progress indicator, Default: 0 (turned off) +#' @param ... params to pass to stan run +#' @return A list containing: +#' \describe{ +#' \item{\code{est}}{Posterior mean of the log odds ratio.} +#' \item{\code{se}}{Posterior standard deviation (i.e., standard error) of the +#' log odds ratio.} +#' } +#' +#' @details +#' This function is useful in single-arm trials where the control arm is not +#' directly observed. A prior on the control arm's response rate provided as a +#' credible interval, is transformed into a normal prior on the logit scale. +#' The posterior distribution of the log odds ratio is then estimated using +#' Bayesian inference via Stan. +#' +#' The Stan model (`estimate_ctrl.stan`) is loaded from the installed package directory. +#' +#' @examples +#' \dontrun{ +#' resp2oddsratio_estimate_ctrl( +#' n_resp_trt = 40, +#' n_trt = 100, +#' low_soc_rr = 0.05, +#' upp_soc_rr = 0.30, +#' ci_rr = 0.80 +#' ) +#' } +#' +#' @rdname resp2oddsratio_estimate_ctrl +#' @export +#' @importFrom stats sd qnorm +resp2oddsratio_estimate_ctrl <- function( + n_resp_trt, + n_trt, + low_soc_rr, + upp_soc_rr, + ci_rr = 0.8, + niter = 1000, + nchains = 1, + ncores = 1, + seed = 123, + refresh = 0, + ... + ) { + + # Compute observed treatment ORR + p_resp_trt <- n_resp_trt / n_trt + + # Compute prior on logit(p_ctrl) + logit_low <- log(low_soc_rr / (1 - low_soc_rr)) + logit_upp <- log(upp_soc_rr / (1 - upp_soc_rr)) + mu_logit_ctrl <- (logit_low + logit_upp) / 2 + sigma_logit_ctrl <- (logit_upp - logit_low) / (2 * qnorm(1 - (1 - ci_rr)/2)) + + # Prepare data list for Stan + stan_data <- list( + n_resp_trt = n_resp_trt, + n_trt = n_trt, + mu_logit_ctrl = mu_logit_ctrl, + sigma_logit_ctrl = sigma_logit_ctrl + ) + + rstan::rstan_options(auto_write = TRUE) + + stan_file_name <- "estimate_ctrl.stan" + stan_mod <- .load_stan_model(stan_file_name) + # Run Stan model + if (!is.null(stan_mod)) { + fit <- rstan::sampling( + object = stan_mod, + data = stan_data, + iter = niter, + chains = nchains, + cores = ncores, + seed = seed, + refresh = refresh, + ... + ) + } else { + file_path <- system.file('stan', stan_file_name, package = 'oncoPoS') + fit <- rstan::stan( + file = file_path, + data = stan_data, + iter = niter, + chains = nchains, + cores = ncores, + seed = seed, + refresh = refresh, + ... + ) + .save_stan_model(fit@stanmodel, stan_file_name) + } + + # Extract log(OR) samples + posterior <- rstan::extract(fit) + log_or_samples <- posterior$log_or + + # Return posterior mean and sd + return(list(est = mean(log_or_samples), + se = sd(log_or_samples))) +} diff --git a/R/run_stan.R b/R/run_stan.R index 475be35..456b461 100644 --- a/R/run_stan.R +++ b/R/run_stan.R @@ -1,4 +1,4 @@ -#' @title Run stan to generate a PoS prediction +#' @title Run Stan to Predict Probability of Success (PoS) for Phase III Trial #' @description A wrapper function to run stan in order to generate a PoS #' prediction #' @@ -9,11 +9,17 @@ #' \item When `use_PFS` is TRUE, calculate log hazard ratio PFS estimate and #' its standard error from a phase 2 study #' \item When `use_ORR` is TRUE, calculate log odds ratio for ORR estimate -#' and its standard error estimate from a phase 2 study +#' and its standard error estimate using either two-arm or single-arm data +#' from a phase 2 study +#' \item Construct a Beta prior for \code{omega} using user-specified +#' mean and variance. #' \item Derive sigma_P1 and sigma_P2, i.e., the standard deviations in the #' mixture prior, based on the input threshold value #' \item Calculate the unit standard deviation and the covariance matrix of #' the observed treatment effect using `gen_sigma()` +#' \item Select indication-specific regression parameters for surrogate +#' -primary endpoint relationship when the primary endpoint (e.g., PFS) is +#' not directly available from the earlier-phase study. #' \item Derive the hyperparameter for the half normal distribution #' \item Create a list for `rstan::stan()` input #' \item Update the list based on `use_PFS` and `use_ORR` @@ -30,8 +36,10 @@ #' @param nevents3 numeric vector of target number of events in phase 3 study #' @param hr_bound numeric vector of hazard ratio bounds for analyses in phase #' 3 study -#' @param omega probability that a treatment effect comes from a enthusiastic -#' prior component, i.e., initial benchmarking probability for the study success +#' @param omega_mean Mean of the Beta prior for \code{omega}. \code{omega} is +#' the probability that a treatment effect comes from a enthusiastic prior +#' component, i.e., initial benchmarking probability for the study success +#' @param omega_var Variance of the Beta prior for \code{omega}. #' @param est_obs_pfs estimated PFS hazard ratio based on prior/earlier study #' @param low_obs_pfs lower bound of estimated PFS hazard ratio based on #' prior/earlier study @@ -48,8 +56,15 @@ #' study #' @param n_trt2 sample size in treatment arm from a prior/earlier study #' @param n_ctrl2 sample size in control arm from a prior/earlier study +#' @param low_soc_rr Lower bound of historical control response rate for +#' single-arm estimation. +#' @param upp_soc_rr Upper bound of historical control response rate for +#' single-arm estimation. +#' @param ci_rr Confidence level for control response rate bounds, Default: 0.8. #' @param use_orr whether response data from a prior/earlier study should be #' used, Default: FALSE +#' @param single_arm whether ORR data is from a single-arm trial, Default: +#' FALSE. #' @param use_pfs whether PFS data from a prior/earlier study should be used, #' Default: FALSE #' @param het_degree_p2 indicates the degree of heterogeneity of a study level @@ -58,19 +73,33 @@ #' @param het_degree_p3 indicates the degree of heterogeneity of a study level #' parameter for a phase 3 study and must be one of "large", #' "substantial", "moderate", "small", "very small", Default: 'very small' +#' @param indication Integer from 1 to 6 for selecting indication-specific +#' ORR-PFS regression parameters, as follows: +#' 1 = hematologic malignancies +#' 2 = gynecologic cancers +#' 3 = thoracic cancers +#' 4 = other solid tumors +#' 5 = breast cancer +#' 6 = any tumor type +#' Default: 6 #' @param m_0 intercept for linear regression of log treatment effect of PFS on -#' log treatment effect on response. A value is expected only when `use_orr = TRUE` -#' and `use_pfs = TRUE`, Default: NA +#' log treatment effect on response. A value is expected only when +#' `use_orr = TRUE` and `use_pfs = TRUE`. Auto-filled based on \code{indication} +#' if not supplied. #' @param m_1 slope for linear regression of log treatment effect of PFS on -#' log treatment effect on response. A value is expected only when `use_orr = TRUE` -#' and `use_pfs = TRUE`, Default: NA +#' log treatment effect on response. A value is expected only when +#' `use_orr = TRUE` and `use_pfs = TRUE`. Auto-filled based on \code{indication} +#' if not supplied. #' @param nu_0 standard error of `m_0`. A value is expected only when -#' `use_orr = TRUE` and `use_pfs = TRUE`, Default: NA +#' `use_orr = TRUE` and `use_pfs = TRUE`. Auto-filled based on \code{indication} +#' if not supplied. #' @param nu_1 standard error of `m_1`. A value is expected only when -#' `use_orr = TRUE` and `use_pfs = TRUE`, Default: NA +#' `use_orr = TRUE` and `use_pfs = TRUE`. Auto-filled based on \code{indication} +#' if not supplied. #' @param lm_sd linear regression residual variance of log treatment effect of #' PFS on log treatment effect on response. A value is expected only when -#' `use_orr = TRUE` and `use_pfs = TRUE`, Default: NA +#' `use_orr = TRUE` and `use_pfs = TRUE`. Auto-filled based on \code{indication} +#' if not supplied. #' @param niter number of iterations to be used in stan run, Default: 1000 #' @param nchains number of chains to be used in stan run, Default: 4 #' @param ncores number of cores to be used in stan run, Default: 4 @@ -80,20 +109,22 @@ #' @return list which includes the generated stan object, list of data that #' was supplied to `rstan::stan()`, and the name of the stan file which was run #' @examples -#' # use PFS data from a prior study +#' # use single arm ORR data from a prior study #' run_stan( -#' 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) +#' target_hr = 0.70, +#' J = 2, +#' nevents3 = c(370, 468), +#' hr_bound = c(0.7790, 0.8204), +#' thres = 0.01, +#' 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=111) #' @seealso #' \code{\link[stats]{Normal}} #' \code{\link[pracma]{erf}} @@ -107,7 +138,8 @@ run_stan <- function( target_hr, J, nevents3, - hr_bound, omega, + hr_bound, + omega_mean = 0.52, omega_var = 0.02, est_obs_pfs, low_obs_pfs, upp_obs_pfs, @@ -117,11 +149,16 @@ run_stan <- 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, @@ -132,12 +169,13 @@ run_stan <- function( ncores = 4, seed, ... - ) { +) { # check the 'three dots' input - #add an error if the length of the nevenet s, hr_bound equal to J - + #add an error if the length of the nevent s, hr_bound equal to J + delta_P <- log(target_hr) loghr_bound <- log(hr_bound) + params <- get_beta_params(omega_mean, omega_var) #log HR PFS estimate and its SE from phase II: if (use_pfs) { @@ -147,13 +185,28 @@ run_stan <- function( #log Odds Ratio for ORR estimate and its SE from phase II: if (use_orr) { - #orr % estimates in each arm: - log_odds_est <- resp2oddsratio( - n_resp_trt = n_resp_trt2, - n_resp_ctrl = n_resp_ctrl2, - n_trt = n_trt2, - n_ctrl = n_ctrl2) - } + if (!single_arm){ + log_odds_est <- resp2oddsratio( + n_resp_trt = n_resp_trt2, + n_trt = n_trt2, + n_resp_ctrl = n_resp_ctrl2, + n_ctrl = n_ctrl2) + n_patients <- n_trt2 + n_ctrl2 # Use two-arm data + } + if (single_arm){ + log_odds_est <- resp2oddsratio_estimate_ctrl( + n_resp_trt = n_resp_trt2, + n_trt = n_trt2, + low_soc_rr = low_soc_rr, + upp_soc_rr = upp_soc_rr, + ci_rr = ci_rr, + niter = niter, + nchains = nchains, + seed = seed, + ...) + n_patients <- n_trt2 # Use single-arm data + } + } sigma_P1 <- sigma_P2 <- delta_P/stats::qnorm(thres) @@ -165,9 +218,24 @@ run_stan <- function( z2 <- sigma_unit*get_denominator(het_degree_p2)/pracma::erfinv(0.5)/sqrt(2) z3 <- sigma_unit*get_denominator(het_degree_p3)/pracma::erfinv(0.5)/sqrt(2) # (very small degree of heterogeneity for phase 3) + + # regression parameters getting from Indication-Specific Bayesian Hierarchical Model + m_0_vec <- c(-0.005942885, -0.189576348, -0.154820463, -0.023896981, 0.209979369, -0.03014321) + nu_0_vec <- c(0.1648034, 0.2021437, 0.1637840, 0.1596490, 0.1559987, 0.1789729) + m_1_vec <- c(1.450295, 2.118675, 1.502140, 1.698569, 2.061274, 1.771319) + nu_1_vec <- c(0.4741983, 0.5534602, 0.3682789, 0.3440543, 0.2589501, 0.3598188) + if (is.na(lm_sd)) lm_sd <- 2.02218 #mean(posterior_draws$sigma_wls) + + if (is.na(m_0)) m_0 <- m_0_vec[indication] + if (is.na(m_1)) m_1 <- m_1_vec[indication] + if (is.na(nu_0)) nu_0 <- nu_0_vec[indication] + if (is.na(nu_1)) nu_1 <- nu_1_vec[indication] + #define the data list for the stan model based on the available phase II data stan_list <- list( - omega = omega, + #omega = omega, + omega_alpha = params$alpha, + omega_beta = params$beta, delta_P = delta_P, sigma_P1 = sigma_P1, sigma_P2 = sigma_P2, @@ -184,7 +252,7 @@ run_stan <- function( list( theta_hat = lest_obs_pfs, theta_hat_sd = lse_obs_pfs - ) + ) ) stan_file <- "phase23_interim_pfs.stan" } @@ -196,8 +264,8 @@ run_stan <- function( orr_hat = log_odds_est$est, orr_hat_sd = log_odds_est$se, m_0 = m_0, m_1 = m_1, nu_0 = nu_0, nu_1 = nu_1, - wls_sd = lm_sd/sqrt(n_trt2 + n_ctrl2) - ) + wls_sd = lm_sd/sqrt(n_patients) + ) ) stan_file <- "phase23_interim_orr.stan" @@ -209,19 +277,31 @@ run_stan <- function( } rstan::rstan_options(auto_write = TRUE) - - package_path <- find.package('oncoPoS',lib.loc = .libPaths()) - file_path <- file.path(package_path,'bin', 'stan', stan_file) - - fit_rstan <- rstan::stan( - file = file_path, - data = stan_list, - iter = niter, - chains = nchains, - cores = ncores, - seed = seed, - ... + + stan_mod <- .load_stan_model(stan_file) + if (!is.null(stan_mod)) { + fit_rstan <- rstan::sampling( + object = stan_mod, + data = stan_list, + iter = niter, + chains = nchains, + cores = ncores, + seed = seed, + ... ) - + } else { + file_path <- system.file('stan', stan_file, package = 'oncoPoS') + fit_rstan <- rstan::stan( + file = file_path, + data = stan_list, + iter = niter, + chains = nchains, + cores = ncores, + seed = seed, + ... + ) + .save_stan_model(fit_rstan@stanmodel, stan_file) + } + return(list(fit_rstan = fit_rstan, stan_list = stan_list, stan_file = stan_file)) } diff --git a/R/stan_platform.R b/R/stan_platform.R new file mode 100644 index 0000000..74e1d5e --- /dev/null +++ b/R/stan_platform.R @@ -0,0 +1,246 @@ +#' @importFrom rstan sampling stan_model +NULL + +# --------------------------------------------------------------------------- +# Internal helpers for platform-aware Stan model caching +# +# Two-tier cache lookup (checked in order): +# +# 1. Package bundle -- inst/stan/precompiled// +# Committed to version control. Covers platforms that the maintainer has +# pre-compiled for. Read-only after installation. +# +# 2. User cache -- tools::R_user_dir("oncoPoS", "cache")// +# Always writable (Linux system installs, read-only HPC mounts, etc.). +# Populated automatically on first use for any platform not in the bundle. +# Survives package upgrades until the platform key changes. +# +# The platform key encodes OS, CPU architecture, R version, rstan version, and +# StanHeaders version -- the full set of factors that determine compiled model +# compatibility. +# --------------------------------------------------------------------------- + +#' Build a platform-specific cache key for Stan models +#' +#' Combines OS name, CPU architecture, R version, rstan version, and +#' StanHeaders version into a single string. Two sessions that produce the +#' same key can share compiled \code{.rds} files without recompilation. +#' +#' @return A character scalar, e.g. +#' \code{"Windows_x86-64_4.5.1_2.21.9_2.26.28"}. +#' @keywords internal +.stan_platform_key <- function() { + info <- Sys.info() + paste( + info[["sysname"]], + info[["machine"]], + paste(R.version$major, R.version$minor, sep = "."), + as.character(utils::packageVersion("rstan")), + as.character(utils::packageVersion("StanHeaders")), + sep = "_" + ) +} + +#' Locate the bundled precompiled directory for the current platform +#' +#' Searches \code{inst/stan/precompiled/} inside the installed package. +#' +#' @param platform_key Character scalar from \code{.stan_platform_key()}. +#' @return Absolute path to the directory, or \code{""} if it does not exist. +#' @keywords internal +.stan_precompiled_dir <- function(platform_key) { + system.file("stan", "precompiled", platform_key, package = "oncoPoS") +} + +#' Locate the user-writable cache directory for the current platform +#' +#' Uses \code{tools::R_user_dir("oncoPoS", "cache")} so it is always writable +#' regardless of where the package is installed (system library, HPC, etc.). +#' +#' @param platform_key Character scalar from \code{.stan_platform_key()}. +#' @return Absolute path to the directory (created if necessary). +#' @keywords internal +.stan_user_cache_dir <- function(platform_key) { + base <- tools::R_user_dir("oncoPoS", "cache") + file.path(base, platform_key) +} + +#' Load a pre-compiled Stan model for the current platform +#' +#' Checks the bundled \code{inst/stan/precompiled/} directory first, then +#' falls back to the user-writable cache (\code{tools::R_user_dir()}). +#' +#' @param stan_file Filename of the Stan model (e.g. +#' \code{"phase23_interim_none.stan"}). +#' @return A \code{stanmodel} object if a matching pre-compiled \code{.rds} +#' exists in either cache tier, otherwise \code{NULL}. +#' @keywords internal +.load_stan_model <- function(stan_file) { + key <- .stan_platform_key() + rds_name <- sub("\\.stan$", ".rds", stan_file) + + # Tier 1: bundled precompiled (committed to version control) + bundle_dir <- .stan_precompiled_dir(key) + if (nzchar(bundle_dir)) { + rds_path <- file.path(bundle_dir, rds_name) + if (file.exists(rds_path)) { + mod <- tryCatch(readRDS(rds_path), error = function(e) NULL) + if (!is.null(mod)) return(mod) + } + } + + # Tier 2: user cache (always writable) + cache_dir <- .stan_user_cache_dir(key) + rds_path <- file.path(cache_dir, rds_name) + if (file.exists(rds_path)) { + tryCatch(readRDS(rds_path), error = function(e) NULL) + } else { + NULL + } +} + +#' Save a compiled Stan model to the user-writable cache +#' +#' Writes the \code{stanmodel} to \code{tools::R_user_dir("oncoPoS", "cache")} +#' so it persists across sessions on any platform, including read-only system +#' library installations. Creates a \code{manifest.dcf} on the first call. +#' Silently no-ops if the cache directory cannot be created or written. +#' +#' @param stan_model A \code{stanmodel} object (e.g. \code{fit@@stanmodel}). +#' @param stan_file Filename of the Stan source (e.g. +#' \code{"phase23_interim_none.stan"}). +#' @keywords internal +.save_stan_model <- function(stan_model, stan_file) { + key <- .stan_platform_key() + out_dir <- .stan_user_cache_dir(key) + dir.create(out_dir, recursive = TRUE, showWarnings = FALSE) + + rds_path <- file.path(out_dir, sub("\\.stan$", ".rds", stan_file)) + tryCatch( + saveRDS(stan_model, rds_path), + error = function(e) invisible(NULL) + ) + + # Write manifest once per platform directory + manifest_path <- file.path(out_dir, "manifest.dcf") + if (!file.exists(manifest_path)) { + tryCatch({ + info <- Sys.info() + write.dcf( + data.frame( + sysname = info[["sysname"]], + machine = info[["machine"]], + r_version = paste(R.version$major, R.version$minor, sep = "."), + rstan_version = as.character(utils::packageVersion("rstan")), + stanheaders_version = as.character(utils::packageVersion("StanHeaders")), + compiled_at = format(Sys.time(), "%Y-%m-%dT%H:%M:%S"), + stringsAsFactors = FALSE + ), + manifest_path + ) + }, error = function(e) invisible(NULL)) + } + + invisible(out_dir) +} + +# --------------------------------------------------------------------------- +# Exported developer utility +# --------------------------------------------------------------------------- + +#' @title Compile all Stan models and cache them for the current platform +#' +#' @description +#' Run this function once from the package root directory when setting up a +#' new computing environment, or after modifying Stan source files. +#' +#' Each \code{.stan} model in \code{stan_src_dir} is compiled with +#' \code{\link[rstan]{stan_model}} and saved to +#' \code{/precompiled//}. A human-readable +#' \code{manifest.dcf} records the build environment for traceability. +#' +#' \strong{Committing the \code{precompiled/} tree to version control} lets +#' collaborators on identical platforms (same OS, CPU arch, R version, rstan +#' version, and StanHeaders version) skip recompilation entirely -- analogous to +#' how \pkg{renv} stores per-platform package binaries in +#' \code{renv/library///}. +#' +#' @section Platform key: +#' The subdirectory name is constructed by \code{.stan_platform_key()} and +#' has the form \cr +#' \code{____} +#' \cr e.g. \code{Windows_x86-64_4.5.1_2.21.9_2.26.28}. +#' +#' @section Traceability: +#' Each platform directory contains a \code{manifest.dcf} file with fields: +#' \describe{ +#' \item{sysname}{Operating system (from \code{Sys.info()}).} +#' \item{machine}{CPU architecture.} +#' \item{r_version}{R major.minor version.} +#' \item{rstan_version}{rstan package version.} +#' \item{stanheaders_version}{StanHeaders package version.} +#' \item{compiled_at}{ISO-8601 timestamp of compilation.} +#' \item{stan_files}{Comma-separated list of compiled model filenames.} +#' } +#' +#' @param stan_src_dir Character. Path to the \code{inst/stan/} directory. +#' Defaults to \code{"inst/stan"}, which is correct when running from the +#' package root. +#' @param verbose Logical; if \code{TRUE} (default), print progress messages. +#' +#' @return Invisibly, the path to the platform-specific precompiled directory. +#' +#' @examples +#' \dontrun{ +#' # Run once from the package root to cache Stan models for this machine +#' compile_stan_models() +#' +#' # Or point at an arbitrary Stan source directory +#' compile_stan_models(stan_src_dir = "path/to/stan") +#' } +#' +#' @export +compile_stan_models <- function(stan_src_dir = file.path("inst", "stan"), + verbose = TRUE) { + key <- .stan_platform_key() + out_dir <- file.path(stan_src_dir, "precompiled", key) + dir.create(out_dir, recursive = TRUE, showWarnings = FALSE) + + stan_files <- list.files(stan_src_dir, pattern = "\\.stan$", full.names = TRUE) + if (length(stan_files) == 0L) { + stop("No .stan files found in: ", stan_src_dir) + } + + if (verbose) message("Platform key: ", key) + + for (stan_path in stan_files) { + stan_name <- basename(stan_path) + if (verbose) message("Compiling ", stan_name, " ...") + mod <- rstan::stan_model(file = stan_path, verbose = FALSE) + rds_name <- sub("\\.stan$", ".rds", stan_name) + saveRDS(mod, file.path(out_dir, rds_name)) + if (verbose) message(" -> saved ", rds_name) + } + + info <- Sys.info() + write.dcf( + data.frame( + sysname = info[["sysname"]], + machine = info[["machine"]], + r_version = paste(R.version$major, R.version$minor, sep = "."), + rstan_version = as.character(utils::packageVersion("rstan")), + stanheaders_version = as.character(utils::packageVersion("StanHeaders")), + compiled_at = format(Sys.time(), "%Y-%m-%dT%H:%M:%S"), + stan_files = paste(basename(stan_files), collapse = ", "), + stringsAsFactors = FALSE + ), + file.path(out_dir, "manifest.dcf") + ) + + if (verbose) { + message("All models compiled.") + message("Manifest written to: ", file.path(out_dir, "manifest.dcf")) + } + + invisible(out_dir) +} diff --git a/README.md b/README.md index a2297e5..2297745 100644 --- a/README.md +++ b/README.md @@ -10,6 +10,43 @@ oncology study using a Bayesian Hierarchical model. The prior in this model is based on the data observed in an earlier study, design features for the phase 3 study that is being assessed and the industry benchmark for success. +The benchmark PoS is specified via a Beta distribution using its +**mean and variance**, which can be derived from historical success rates or +predictive models (e.g., random forest). +The framework supports both **single-arm and two-arm** early-phase designs and +allows for **indication-specific ORR–PFS regression** when surrogate endpoints +are used. + +## Indication-specific ORR-PFS regression +When early endpoint objective response rate (ORR) is used to predict phase 3 outcome +progression-free survival (PFS), the relationship varies by cancer type. To account +for this heterogeneity, `oncoPoS` groups cancer indications into five categories, +each associated with a distinct set of ORR–PFS regression parameters derived from +prior Bayesian hierarchical modeling. + +**Indication Groups** + +- **Group 1: Hematologic malignancies** + Includes classical Hodgkin lymphoma (CHL), diffuse large B-cell lymphoma (DLBCL), + follicular lymphoma (FL), multiple myeloma (MM), non-Hodgkin lymphoma (NHL), + and peripheral T-cell lymphoma (PTCL). + +- **Group 2: Gynecologic cancers** + Includes cervical, endometrial, and ovarian cancers. + +- **Group 3: Thoracic malignancies** + Includes non-small cell lung cancer (NSCLC), small cell lung cancer (SCLC), + and mesothelioma. + +- **Group 4: Urologic and gastrointestinal solid tumors** + Includes bladder cancer, gastric cancer, and renal cell carcinoma (RCC). + +- **Group 5: Breast cancer** + Includes breast cancer. + +If no indication is specified by the user, `oncoPoS` defaults to using the +average ORR–PFS relationship across all indication groups. + ## Installation @@ -19,7 +56,7 @@ You can install development version of `oncoPoS` from GitHub with: if (!requireNamespace("remotes")) { install.packages("remotes") } -remotes::install_github("MSDLLCPapers/oncoPoS") +remotes::install_github("MSDLLCPapers/oncoPoS", ref = "develop") ``` ## Example @@ -27,6 +64,10 @@ remotes::install_github("MSDLLCPapers/oncoPoS") Below is a simple example which assumes the following information is available for the an oncology phase 3 trial PoS calculation: +- Disease area: + - This trial is planned for breast cancer, which falls under Indication Group + 5 in the model’s ORR–PFS regression framework. + - Design features of a phase 3 trial: - Progression-free survival (PFS) is the primary endpoint with the target hazard ratio (HR) of 0.7; @@ -38,27 +79,40 @@ for the an oncology phase 3 trial PoS calculation: - The above phase 3 trial is planned following promising results in an earlier phase 2 study, which reported PFS HR (95% CI) of 0.53 (0.31, 0.91). -- The industry benchmark for success is 0.57 for this type of phase 3 trial. +- In addition to PFS, objective response rate (ORR) data is available from a + single-arm phase 2 trial, with 40 responders out of 100 patients and + historical control response rate assumed to lie between 0.05 and 0.20. + +- The benchmark probability of success (`omega`) is modeled using a Beta prior + with mean 0.3 and variance 0.03, reflecting prior belief about the historical + success rate in similar studies. All the above information is synthesized in `oncoPoS::gen_pos()` using Bayesian Hierarchical model to generate a PoS estimate: ```r gen_pos( - target_hr = 0.70, - J = 2, - nevents3 = c(370, 468), - hr_bound = c(0.7790, 0.8204), - est_obs_pfs = 0.53, - low_obs_pfs = 0.31, - upp_obs_pfs = 0.91, - omega = 0.57, - seed = 574, - ratio = 2, - use_pfs = TRUE, - nchains = 4, - niter = 1000 -) + target_hr = 0.70, + J = 2, + nevents3 = c(370, 468), + hr_bound = c(0.7790, 0.8204), + thres = 0.01, + ratio = 2, + omega_mean = 0.3, + omega_var = 0.03, + est_obs_pfs = 0.53, + low_obs_pfs = 0.31, + upp_obs_pfs = 0.91, + use_pfs = TRUE, + indication = 5, + n_trt2 = 100, + n_resp_trt2 = 40, + low_soc_rr = 0.05, + upp_soc_rr = 0.2, + use_orr = TRUE, + single_arm = TRUE, + seed = 222 + ) ``` diff --git a/cleanup b/cleanup deleted file mode 100755 index 1e9fd15..0000000 --- a/cleanup +++ /dev/null @@ -1,30 +0,0 @@ -#!/usr/bin/env sh -# -# This script was copied and modified from the source code of the {configure} -# R package by Kevin Ushey. License and copyright are in the comment below. -# -# MIT License -# ===================== -# -# Copyright 2021 Kevin Ushey -# -# Permission is hereby granted, free of charge, to any person obtaining a copy -# of this software and associated documentation files (the “Software”), to deal -# in the Software without restriction, including without limitation the rights -# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -# copies of the Software, and to permit persons to whom the Software is -# furnished to do so, subject to the following conditions: -# -# The above copyright notice and this permission notice shall be included -# in all copies or substantial portions of the Software. -# -# THE SOFTWARE IS PROVIDED “AS IS”, WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL -# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN -# THE SOFTWARE. -if [ -d "bin" ]; then - rm -rf bin -fi diff --git a/cleanup.win b/cleanup.win deleted file mode 100755 index db713b0..0000000 --- a/cleanup.win +++ /dev/null @@ -1,31 +0,0 @@ -#!/usr/bin/env sh -# -# This script was copied and modified from the source code of the {configure} -# R package by Kevin Ushey. License and copyright are in the comment below. -# -# MIT License -# ===================== -# -# Copyright 2021 Kevin Ushey -# -# Permission is hereby granted, free of charge, to any person obtaining a copy -# of this software and associated documentation files (the “Software”), to deal -# in the Software without restriction, including without limitation the rights -# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -# copies of the Software, and to permit persons to whom the Software is -# furnished to do so, subject to the following conditions: -# -# The above copyright notice and this permission notice shall be included -# in all copies or substantial portions of the Software. -# -# THE SOFTWARE IS PROVIDED “AS IS”, WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL -# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN -# THE SOFTWARE. -if [ -d "bin" ]; then - rm -rf bin -fi - diff --git a/inst/stan/estimate_ctrl.stan b/inst/stan/estimate_ctrl.stan new file mode 100644 index 0000000..75307a4 --- /dev/null +++ b/inst/stan/estimate_ctrl.stan @@ -0,0 +1,23 @@ +data { + int n_resp_trt; // responders in treatment arm + int n_trt; // total in treatment arm + real mu_logit_ctrl; // mean of logit(p_resp_soc) + real sigma_logit_ctrl; // SD of logit(p_resp_soc) +} + +parameters { + real p_resp_trt; // treatment response probability + real logit_p_resp_soc; // latent logit(p_resp_soc) +} + +transformed parameters { + real p_resp_soc = inv_logit(logit_p_resp_soc); + real log_odds_orr_trt = log(p_resp_trt / (1 - p_resp_trt)); + real log_odds_orr_soc = log(p_resp_soc / (1 - p_resp_soc)); + real log_or = log_odds_orr_soc - log_odds_orr_trt; +} + +model { + logit_p_resp_soc ~ normal(mu_logit_ctrl, sigma_logit_ctrl); // prior on control arm ORR + n_resp_trt ~ binomial(n_trt, p_resp_trt); // Binomial likelihood for treatment arm +} diff --git a/src/stan/phase23_interim_both.stan b/inst/stan/phase23_interim_both.stan similarity index 86% rename from src/stan/phase23_interim_both.stan rename to inst/stan/phase23_interim_both.stan index bcf6687..91cacb9 100644 --- a/src/stan/phase23_interim_both.stan +++ b/inst/stan/phase23_interim_both.stan @@ -7,7 +7,9 @@ data { // hyperparameters: real tau_sd2; // half normal prior variance real tau_sd3; // half normal prior variance - real omega; // mixture weight + //real omega; // mixture weight + real omega_alpha; + real omega_beta; real delta_P; // mean for pessimistic scenario real sigma_P1; // stdev for optimistic scenario real sigma_P2; // stdev for pessimistic scenario @@ -38,6 +40,8 @@ parameters { // phase III parameters real theta_P2_raw; // real phase 3 tmt effects real theta_P3_raw; // real phase 3 tmt effects + + real omega; } transformed parameters { @@ -47,12 +51,14 @@ transformed parameters { real beta_1; real orr_P; + //orr_P = mu_P + tau_P2 * orr_P_raw; theta_P2 = mu_P + tau_P2*theta_P2_raw; theta_P3 = mu_P + tau_P3*theta_P3_raw; beta_0 = m_0 + nu_0*beta_0_raw; beta_1 = m_1 + nu_1*beta_1_raw; + //theta_P2 = beta_0 + beta_1*orr_P + wls_sd*theta_P2_raw; //orr_P as the independent variable orr_P = beta_0 + beta_1*theta_P2 + wls_sd*orr_P_raw; } @@ -60,6 +66,8 @@ model { // stan will try to find posterior of them here // linear relationship // population level + //omega ~ beta(2, 2); // weakly informative prior + omega ~ beta(omega_alpha, omega_beta); // informative prior from step 1 target += log_mix(omega, normal_lpdf(mu_P|delta_P,sigma_P1), normal_lpdf(mu_P|0,sigma_P2)); // regression parameter; diff --git a/src/stan/phase23_interim_none.stan b/inst/stan/phase23_interim_none.stan similarity index 87% rename from src/stan/phase23_interim_none.stan rename to inst/stan/phase23_interim_none.stan index e66fcbd..d0d7356 100644 --- a/src/stan/phase23_interim_none.stan +++ b/inst/stan/phase23_interim_none.stan @@ -3,7 +3,9 @@ data { // hyperparameters: real tau_sd2; real tau_sd3; // half normal prior variance - real omega; // mixture weight + //real omega; // mixture weight + real omega_alpha; + real omega_beta; real delta_P; // mean for pessimistic scenario real sigma_P1; // stdev for optimistic scenario real sigma_P2; // stdev for pessimistic scenario @@ -20,13 +22,14 @@ transformed data { } parameters { + real omega; + real mu_P; real tau_P2; real tau_P3; // study-level treatment variance at phase 3 // phase III parameters real theta_P2_raw; // real phase 3 tmt effects real theta_P3_raw; // real phase 3 tmt effects - } transformed parameters { @@ -41,6 +44,8 @@ model { // stan will try to find posterior of them here // linear relationship // population level + //omega ~ beta(2, 2); // weakly informative prior + omega ~ beta(omega_alpha, omega_beta); // informative prior from step 1 target += log_mix(omega, normal_lpdf(mu_P|delta_P,sigma_P1), normal_lpdf(mu_P|0,sigma_P2)); // study level diff --git a/src/stan/phase23_interim_orr.stan b/inst/stan/phase23_interim_orr.stan similarity index 86% rename from src/stan/phase23_interim_orr.stan rename to inst/stan/phase23_interim_orr.stan index 9094e9b..40eee5e 100644 --- a/src/stan/phase23_interim_orr.stan +++ b/inst/stan/phase23_interim_orr.stan @@ -5,7 +5,9 @@ data { // hyperparameters: real tau_sd2; // half normal prior variance real tau_sd3; // half normal prior variance - real omega; // mixture weight + //real omega; // mixture weight + real omega_alpha; + real omega_beta; real delta_P; // mean for pessimistic scenario real sigma_P1; // stdev for optimistic scenario real sigma_P2; // stdev for pessimistic scenario @@ -37,7 +39,7 @@ parameters { real theta_P2_raw; real theta_P3_raw; - + real omega; } @@ -48,18 +50,23 @@ transformed parameters { real beta_1; real orr_P; + //orr_P = mu_P + tau_P2 * orr_P_raw; theta_P2 = mu_P + tau_P2*theta_P2_raw; theta_P3 = mu_P + tau_P3*theta_P3_raw; beta_0 = m_0 + nu_0*beta_0_raw; beta_1 = m_1 + nu_1*beta_1_raw; + //theta_P2 = beta_0 + beta_1*orr_P + wls_sd*theta_P2_raw; //orr_P as the independent variable orr_P = beta_0 + beta_1*theta_P2 + wls_sd*orr_P_raw; } model { + // population level + omega ~ beta(omega_alpha, omega_beta); // informative prior from step 1 + target += log_mix(omega, normal_lpdf(mu_P|delta_P,sigma_P1), normal_lpdf(mu_P|0,sigma_P2)); // regression parameter; diff --git a/src/stan/phase23_interim_pfs.stan b/inst/stan/phase23_interim_pfs.stan similarity index 87% rename from src/stan/phase23_interim_pfs.stan rename to inst/stan/phase23_interim_pfs.stan index d894c3d..2a807e4 100644 --- a/src/stan/phase23_interim_pfs.stan +++ b/inst/stan/phase23_interim_pfs.stan @@ -5,7 +5,9 @@ data { real theta_hat_sd; // sampling error variance for theta_hat real tau_sd2; // half normal prior variance real tau_sd3; // half normal prior variance - real omega; // mixture weight + //real omega; // mixture weight + real omega_alpha; + real omega_beta; real delta_P; // mean for pessimistic scenario real sigma_P1; // stdev for optimistic scenario real sigma_P2; // stdev for pessimistic scenario @@ -25,6 +27,8 @@ parameters { // phase III parameters real theta_P2_raw; // real phase 3 tmt effects real theta_P3_raw; // real phase 3 tmt effects + + real omega; } transformed parameters { @@ -37,6 +41,8 @@ transformed parameters { model { // population level + //omega ~ beta(2, 2); // weakly informative prior + omega ~ beta(omega_alpha, omega_beta); // informative prior from step 1 target += log_mix(omega, normal_lpdf(mu_P|delta_P,sigma_P1), normal_lpdf(mu_P|0,sigma_P2)); // study level diff --git a/inst/stan/precompiled/Windows_x86-64_4.5.1_2.32.7_2.32.10/estimate_ctrl.rds b/inst/stan/precompiled/Windows_x86-64_4.5.1_2.32.7_2.32.10/estimate_ctrl.rds new file mode 100644 index 0000000..d38d98f Binary files /dev/null and b/inst/stan/precompiled/Windows_x86-64_4.5.1_2.32.7_2.32.10/estimate_ctrl.rds differ diff --git a/inst/stan/precompiled/Windows_x86-64_4.5.1_2.32.7_2.32.10/manifest.dcf b/inst/stan/precompiled/Windows_x86-64_4.5.1_2.32.7_2.32.10/manifest.dcf new file mode 100644 index 0000000..03096a2 --- /dev/null +++ b/inst/stan/precompiled/Windows_x86-64_4.5.1_2.32.7_2.32.10/manifest.dcf @@ -0,0 +1,9 @@ +sysname: Windows +machine: x86-64 +r_version: 4.5.1 +rstan_version: 2.32.7 +stanheaders_version: 2.32.10 +compiled_at: 2026-06-03T14:21:39 +stan_files: estimate_ctrl.stan, phase23_interim_both.stan, + phase23_interim_none.stan, phase23_interim_orr.stan, + phase23_interim_pfs.stan diff --git a/inst/stan/precompiled/Windows_x86-64_4.5.1_2.32.7_2.32.10/phase23_interim_both.rds b/inst/stan/precompiled/Windows_x86-64_4.5.1_2.32.7_2.32.10/phase23_interim_both.rds new file mode 100644 index 0000000..3a41810 Binary files /dev/null and b/inst/stan/precompiled/Windows_x86-64_4.5.1_2.32.7_2.32.10/phase23_interim_both.rds differ diff --git a/inst/stan/precompiled/Windows_x86-64_4.5.1_2.32.7_2.32.10/phase23_interim_none.rds b/inst/stan/precompiled/Windows_x86-64_4.5.1_2.32.7_2.32.10/phase23_interim_none.rds new file mode 100644 index 0000000..fd2c4d0 Binary files /dev/null and b/inst/stan/precompiled/Windows_x86-64_4.5.1_2.32.7_2.32.10/phase23_interim_none.rds differ diff --git a/inst/stan/precompiled/Windows_x86-64_4.5.1_2.32.7_2.32.10/phase23_interim_orr.rds b/inst/stan/precompiled/Windows_x86-64_4.5.1_2.32.7_2.32.10/phase23_interim_orr.rds new file mode 100644 index 0000000..3fb5475 Binary files /dev/null and b/inst/stan/precompiled/Windows_x86-64_4.5.1_2.32.7_2.32.10/phase23_interim_orr.rds differ diff --git a/inst/stan/precompiled/Windows_x86-64_4.5.1_2.32.7_2.32.10/phase23_interim_pfs.rds b/inst/stan/precompiled/Windows_x86-64_4.5.1_2.32.7_2.32.10/phase23_interim_pfs.rds new file mode 100644 index 0000000..a887edd Binary files /dev/null and b/inst/stan/precompiled/Windows_x86-64_4.5.1_2.32.7_2.32.10/phase23_interim_pfs.rds differ diff --git a/man/compile_stan_models.Rd b/man/compile_stan_models.Rd new file mode 100644 index 0000000..ba3ebb7 --- /dev/null +++ b/man/compile_stan_models.Rd @@ -0,0 +1,65 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/stan_platform.R +\name{compile_stan_models} +\alias{compile_stan_models} +\title{Compile all Stan models and cache them for the current platform} +\usage{ +compile_stan_models(stan_src_dir = file.path("inst", "stan"), verbose = TRUE) +} +\arguments{ +\item{stan_src_dir}{Character. Path to the \code{inst/stan/} directory. +Defaults to \code{"inst/stan"}, which is correct when running from the +package root.} + +\item{verbose}{Logical; if \code{TRUE} (default), print progress messages.} +} +\value{ +Invisibly, the path to the platform-specific precompiled directory. +} +\description{ +Run this function once from the package root directory when setting up a +new computing environment, or after modifying Stan source files. + +Each \code{.stan} model in \code{stan_src_dir} is compiled with +\code{\link[rstan]{stan_model}} and saved to +\code{/precompiled//}. A human-readable +\code{manifest.dcf} records the build environment for traceability. + +\strong{Committing the \code{precompiled/} tree to version control} lets +collaborators on identical platforms (same OS, CPU arch, R version, rstan +version, and StanHeaders version) skip recompilation entirely -- analogous to +how \pkg{renv} stores per-platform package binaries in +\code{renv/library///}. +} +\section{Platform key}{ + +The subdirectory name is constructed by \code{.stan_platform_key()} and +has the form \cr +\code{____} +\cr e.g. \code{Windows_x86-64_4.5.1_2.21.9_2.26.28}. +} + +\section{Traceability}{ + +Each platform directory contains a \code{manifest.dcf} file with fields: +\describe{ +\item{sysname}{Operating system (from \code{Sys.info()}).} +\item{machine}{CPU architecture.} +\item{r_version}{R major.minor version.} +\item{rstan_version}{rstan package version.} +\item{stanheaders_version}{StanHeaders package version.} +\item{compiled_at}{ISO-8601 timestamp of compilation.} +\item{stan_files}{Comma-separated list of compiled model filenames.} +} +} + +\examples{ +\dontrun{ +# Run once from the package root to cache Stan models for this machine +compile_stan_models() + +# Or point at an arbitrary Stan source directory +compile_stan_models(stan_src_dir = "path/to/stan") +} + +} diff --git a/man/dot-load_stan_model.Rd b/man/dot-load_stan_model.Rd new file mode 100644 index 0000000..46f5e46 --- /dev/null +++ b/man/dot-load_stan_model.Rd @@ -0,0 +1,20 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/stan_platform.R +\name{.load_stan_model} +\alias{.load_stan_model} +\title{Load a pre-compiled Stan model for the current platform} +\usage{ +.load_stan_model(stan_file) +} +\arguments{ +\item{stan_file}{Filename of the Stan model (e.g. +\code{"phase23_interim_none.stan"}).} +} +\value{ +A \code{stanmodel} object if a matching pre-compiled \code{.rds} +exists, otherwise \code{NULL}. +} +\description{ +Load a pre-compiled Stan model for the current platform +} +\keyword{internal} diff --git a/man/dot-save_stan_model.Rd b/man/dot-save_stan_model.Rd new file mode 100644 index 0000000..6fa18da --- /dev/null +++ b/man/dot-save_stan_model.Rd @@ -0,0 +1,21 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/stan_platform.R +\name{.save_stan_model} +\alias{.save_stan_model} +\title{Save a compiled Stan model to the platform-specific cache} +\usage{ +.save_stan_model(stan_model, stan_file) +} +\arguments{ +\item{stan_model}{A \code{stanmodel} object (e.g. \code{fit@stanmodel}).} + +\item{stan_file}{Filename of the Stan source (e.g. +\code{"phase23_interim_none.stan"}).} +} +\description{ +Writes the \code{stanmodel} to the appropriate +\code{precompiled//} subdirectory and creates a +\code{manifest.dcf} on the first call. Silently no-ops if the target +directory is not writable. +} +\keyword{internal} diff --git a/man/dot-stan_platform_key.Rd b/man/dot-stan_platform_key.Rd new file mode 100644 index 0000000..52b5aa5 --- /dev/null +++ b/man/dot-stan_platform_key.Rd @@ -0,0 +1,18 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/stan_platform.R +\name{.stan_platform_key} +\alias{.stan_platform_key} +\title{Build a platform-specific cache key for Stan models} +\usage{ +.stan_platform_key() +} +\value{ +A character scalar, e.g. +\code{"Windows_x86-64_4.5.1_2.21.9_2.26.28"}. +} +\description{ +Combines OS name, CPU architecture, R version, rstan version, and +StanHeaders version into a single string. Two sessions that produce the +same key can share compiled \code{.rds} files without recompilation. +} +\keyword{internal} diff --git a/man/dot-stan_precompiled_dir.Rd b/man/dot-stan_precompiled_dir.Rd new file mode 100644 index 0000000..5d8ac3e --- /dev/null +++ b/man/dot-stan_precompiled_dir.Rd @@ -0,0 +1,20 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/stan_platform.R +\name{.stan_precompiled_dir} +\alias{.stan_precompiled_dir} +\title{Locate the precompiled directory for the current platform} +\usage{ +.stan_precompiled_dir(platform_key) +} +\arguments{ +\item{platform_key}{Character scalar from \code{.stan_platform_key()}.} +} +\value{ +Absolute path to the directory, or \code{""} if it does not exist. +} +\description{ +Searches \code{stan/precompiled/} inside the installed or source +package (works for both \code{devtools::load_all()} and installed contexts +because \code{system.file()} is overridden by pkgload in both cases). +} +\keyword{internal} diff --git a/man/gen_pos.Rd b/man/gen_pos.Rd index b2adc73..28ca858 100644 --- a/man/gen_pos.Rd +++ b/man/gen_pos.Rd @@ -9,7 +9,8 @@ gen_pos( J, nevents3, hr_bound, - omega, + omega_mean = 0.52, + omega_var = 0.02, est_obs_pfs, low_obs_pfs, upp_obs_pfs, @@ -19,19 +20,24 @@ gen_pos( 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, + 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, ... @@ -48,8 +54,11 @@ alternative hypothesis)} \item{hr_bound}{numeric vector of hazard ratio bounds for analyses in phase 3 study} -\item{omega}{probability that a treatment effect comes from a enthusiastic -prior component, i.e., initial benchmarking probability for the study success} +\item{omega_mean}{Mean of the Beta prior for \code{omega}. \code{omega} is +the probability that a treatment effect comes from a enthusiastic prior +component, i.e., initial benchmarking probability for the study success} + +\item{omega_var}{Variance of the Beta prior for \code{omega}.} \item{est_obs_pfs}{estimated PFS hazard ratio based on prior/earlier study} @@ -76,9 +85,20 @@ study} \item{n_resp_ctrl2}{number of responses in control arm from a prior/earlier study} +\item{low_soc_rr}{Lower bound of historical control response rate for +single-arm estimation.} + +\item{upp_soc_rr}{Upper bound of historical control response rate for +single-arm estimation.} + +\item{ci_rr}{Confidence level for control response rate bounds, Default: 0.8.} + \item{use_orr}{whether response data from a prior/earlier study should be used, Default: FALSE} +\item{single_arm}{whether ORR data is from a single-arm trial, Default: +FALSE.} + \item{use_pfs}{whether PFS data from a prior/earlier study should be used, Default: FALSE} @@ -92,23 +112,38 @@ parameter for a phase 3 study and must be one of "large", \item{ratio}{randomization ratio of experimental arm compared to control} +\item{indication}{Integer from 1 to 6 for selecting indication-specific +ORR-PFS regression parameters, as follows: +1 = hematologic malignancies +2 = gynecologic cancers +3 = thoracic cancers +4 = other solid tumors +5 = breast cancer +6 = any tumor type +Default: 6} + \item{m_0}{intercept for linear regression of log treatment effect of PFS on -log treatment effect on response. A value is expected only when \code{use_orr = TRUE} -and \code{use_pfs = TRUE}, Default: NA} +log treatment effect on response. A value is expected only when +\code{use_orr = TRUE} and \code{use_pfs = TRUE}. Auto-filled based on \code{indication} +if not supplied.} \item{m_1}{slope for linear regression of log treatment effect of PFS on -log treatment effect on response. A value is expected only when \code{use_orr = TRUE} -and \code{use_pfs = TRUE}, Default: NA} +log treatment effect on response. A value is expected only when +\code{use_orr = TRUE} and \code{use_pfs = TRUE}. Auto-filled based on \code{indication} +if not supplied.} \item{nu_0}{standard error of \code{m_0}. A value is expected only when -\code{use_orr = TRUE} and \code{use_pfs = TRUE}, Default: NA} +\code{use_orr = TRUE} and \code{use_pfs = TRUE}. Auto-filled based on \code{indication} +if not supplied.} \item{nu_1}{standard error of \code{m_1}. A value is expected only when -\code{use_orr = TRUE} and \code{use_pfs = TRUE}, Default: NA} +\code{use_orr = TRUE} and \code{use_pfs = TRUE}. Auto-filled based on \code{indication} +if not supplied.} \item{lm_sd}{linear regression residual variance of log treatment effect of PFS on log treatment effect on response. A value is expected only when -\code{use_orr = TRUE} and \code{use_pfs = TRUE}, Default: NA} +\code{use_orr = TRUE} and \code{use_pfs = TRUE}. Auto-filled based on \code{indication} +if not supplied.} \item{niter}{number of iterations to be used in stan run, Default: 1000} @@ -124,12 +159,19 @@ printed, Default: FALSE} \item{...}{params to pass to stan run} } \value{ -tibble of PoS estimates and the corresponding standard errors for -each analysis. If \code{plots_out} is turned on, then the MCMC chains mixing and +tibble of PoS estimates, the corresponding standard errors for +each analysis, and the posterior mean and variance for omega. +If \code{plots_out} is turned on, then the MCMC chains mixing and autocorrelation plots are provided as well. } \description{ -PoS estimation +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}{ @@ -140,6 +182,7 @@ PoS estimation \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 } @@ -150,25 +193,37 @@ The contents of this section are shown in PDF user manual only. } \examples{ -# use PFS data from a prior study - - 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) +# Using both ORR and PFS from a prior single-arm study with a Beta prior on +# omega + +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}} } diff --git a/man/get_beta_params.Rd b/man/get_beta_params.Rd new file mode 100644 index 0000000..1f4eb26 --- /dev/null +++ b/man/get_beta_params.Rd @@ -0,0 +1,66 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/get_beta_params.R +\name{get_beta_params} +\alias{get_beta_params} +\title{Generate Beta Distribution Parameters for Omega} +\usage{ +get_beta_params(mean, var) +} +\arguments{ +\item{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.} + +\item{var}{Numeric value, specifying the prior variance of the Beta +distribution. This may reflect uncertainty in the model-based prediction.} +} +\value{ +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.} +} +} +\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. +} +\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). +} +\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. +} +} + +\examples{ +# Example using prior mean and variance estimated from a Benchmark +# probability model +get_beta_params(mean = 0.52, var = 0.02) + +} diff --git a/man/resp2oddsratio_estimate_ctrl.Rd b/man/resp2oddsratio_estimate_ctrl.Rd new file mode 100644 index 0000000..b8f9d1c --- /dev/null +++ b/man/resp2oddsratio_estimate_ctrl.Rd @@ -0,0 +1,97 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/resp2oddsratio_estimate_ctrl.R +\name{resp2oddsratio_estimate_ctrl} +\alias{resp2oddsratio_estimate_ctrl} +\title{Estimate Log Odds Ratio Using Observed Response Data and Prior on +Control Arm} +\usage{ +resp2oddsratio_estimate_ctrl( + n_resp_trt, + n_trt, + low_soc_rr, + upp_soc_rr, + ci_rr = 0.8, + niter = 1000, + nchains = 1, + ncores = 1, + seed = 123, + refresh = 0, + ... +) +} +\arguments{ +\item{n_resp_trt}{number of responses in the treatment arm} + +\item{n_trt}{sample size in the treatment arm} + +\item{low_soc_rr}{lower bound of the control arm's response rate} + +\item{upp_soc_rr}{upper bound of the control arm's response rate} + +\item{ci_rr}{confidence level (e.g., 0.80) for the control arm response rate +interval.} + +\item{niter}{number of iterations to be used in stan run, Default: 1000} + +\item{nchains}{number of chains to be used in stan run, Default: 4} + +\item{ncores}{number of cores to be used in stan run, Default: 4} + +\item{seed}{seed to be used in stan run} + +\item{refresh}{integer, progress indicator, Default: 0 (turned off)} + +\item{...}{params to pass to stan run} +} +\value{ +A list containing: +\describe{ +\item{\code{est}}{Posterior mean of the log odds ratio.} +\item{\code{se}}{Posterior standard deviation (i.e., standard error) of the +log odds ratio.} +} +} +\description{ +Estimates the log odds ratio (log OR) between treatment and control arms +in a single-arm trial by combining observed responses in treatment arm with +prior knowledge of the control response rate. The prior distribution on the +control arm's response rate is specified through a credible interval. +} +\details{ +This function is useful in single-arm trials where the control arm is not +directly observed. A prior on the control arm's response rate provided as a +credible interval, is transformed into a normal prior on the logit scale. +The posterior distribution of the log odds ratio is then estimated using +Bayesian inference via Stan. + +The Stan model (\code{estimate_ctrl.stan}) is loaded from the installed package directory. +} +\section{Specification}{ + +\if{latex}{ + \itemize{ + \item Compute the observed response rate in the treatment arm. + \item Translate the control response rate bounds into a normal prior on + the logit scale. + \item Fit a Bayesian model to estimate the posterior distribution of the + log OR. + \item Return the posterior mean and standard deviation of the log OR. + } +} +\if{html}{ + The contents of this section are shown in the PDF user manual only. +} +} + +\examples{ +\dontrun{ +resp2oddsratio_estimate_ctrl( + n_resp_trt = 40, + n_trt = 100, + low_soc_rr = 0.05, + upp_soc_rr = 0.30, + ci_rr = 0.80 +) +} + +} diff --git a/man/run_stan.Rd b/man/run_stan.Rd index 5ec8bb8..64136ce 100644 --- a/man/run_stan.Rd +++ b/man/run_stan.Rd @@ -2,14 +2,15 @@ % Please edit documentation in R/run_stan.R \name{run_stan} \alias{run_stan} -\title{Run stan to generate a PoS prediction} +\title{Run Stan to Predict Probability of Success (PoS) for Phase III Trial} \usage{ run_stan( target_hr, J, nevents3, hr_bound, - omega, + omega_mean = 0.52, + omega_var = 0.02, est_obs_pfs, low_obs_pfs, upp_obs_pfs, @@ -19,11 +20,16 @@ run_stan( 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, + indication = 6, m_0 = NA, m_1 = NA, nu_0 = NA, @@ -47,8 +53,11 @@ alternative hypothesis)} \item{hr_bound}{numeric vector of hazard ratio bounds for analyses in phase 3 study} -\item{omega}{probability that a treatment effect comes from a enthusiastic -prior component, i.e., initial benchmarking probability for the study success} +\item{omega_mean}{Mean of the Beta prior for \code{omega}. \code{omega} is +the probability that a treatment effect comes from a enthusiastic prior +component, i.e., initial benchmarking probability for the study success} + +\item{omega_var}{Variance of the Beta prior for \code{omega}.} \item{est_obs_pfs}{estimated PFS hazard ratio based on prior/earlier study} @@ -75,9 +84,20 @@ study} \item{n_resp_ctrl2}{number of responses in control arm from a prior/earlier study} +\item{low_soc_rr}{Lower bound of historical control response rate for +single-arm estimation.} + +\item{upp_soc_rr}{Upper bound of historical control response rate for +single-arm estimation.} + +\item{ci_rr}{Confidence level for control response rate bounds, Default: 0.8.} + \item{use_orr}{whether response data from a prior/earlier study should be used, Default: FALSE} +\item{single_arm}{whether ORR data is from a single-arm trial, Default: +FALSE.} + \item{use_pfs}{whether PFS data from a prior/earlier study should be used, Default: FALSE} @@ -91,23 +111,38 @@ parameter for a phase 3 study and must be one of "large", \item{ratio}{randomization ratio of experimental arm compared to control} +\item{indication}{Integer from 1 to 6 for selecting indication-specific +ORR-PFS regression parameters, as follows: +1 = hematologic malignancies +2 = gynecologic cancers +3 = thoracic cancers +4 = other solid tumors +5 = breast cancer +6 = any tumor type +Default: 6} + \item{m_0}{intercept for linear regression of log treatment effect of PFS on -log treatment effect on response. A value is expected only when \code{use_orr = TRUE} -and \code{use_pfs = TRUE}, Default: NA} +log treatment effect on response. A value is expected only when +\code{use_orr = TRUE} and \code{use_pfs = TRUE}. Auto-filled based on \code{indication} +if not supplied.} \item{m_1}{slope for linear regression of log treatment effect of PFS on -log treatment effect on response. A value is expected only when \code{use_orr = TRUE} -and \code{use_pfs = TRUE}, Default: NA} +log treatment effect on response. A value is expected only when +\code{use_orr = TRUE} and \code{use_pfs = TRUE}. Auto-filled based on \code{indication} +if not supplied.} \item{nu_0}{standard error of \code{m_0}. A value is expected only when -\code{use_orr = TRUE} and \code{use_pfs = TRUE}, Default: NA} +\code{use_orr = TRUE} and \code{use_pfs = TRUE}. Auto-filled based on \code{indication} +if not supplied.} \item{nu_1}{standard error of \code{m_1}. A value is expected only when -\code{use_orr = TRUE} and \code{use_pfs = TRUE}, Default: NA} +\code{use_orr = TRUE} and \code{use_pfs = TRUE}. Auto-filled based on \code{indication} +if not supplied.} \item{lm_sd}{linear regression residual variance of log treatment effect of PFS on log treatment effect on response. A value is expected only when -\code{use_orr = TRUE} and \code{use_pfs = TRUE}, Default: NA} +\code{use_orr = TRUE} and \code{use_pfs = TRUE}. Auto-filled based on \code{indication} +if not supplied.} \item{niter}{number of iterations to be used in stan run, Default: 1000} @@ -135,11 +170,17 @@ prediction \item When `use_PFS` is TRUE, calculate log hazard ratio PFS estimate and its standard error from a phase 2 study \item When `use_ORR` is TRUE, calculate log odds ratio for ORR estimate - and its standard error estimate from a phase 2 study + and its standard error estimate using either two-arm or single-arm data + from a phase 2 study + \item Construct a Beta prior for \code{omega} using user-specified + mean and variance. \item Derive sigma_P1 and sigma_P2, i.e., the standard deviations in the mixture prior, based on the input threshold value \item Calculate the unit standard deviation and the covariance matrix of the observed treatment effect using `gen_sigma()` + \item Select indication-specific regression parameters for surrogate + -primary endpoint relationship when the primary endpoint (e.g., PFS) is + not directly available from the earlier-phase study. \item Derive the hyperparameter for the half normal distribution \item Create a list for `rstan::stan()` input \item Update the list based on `use_PFS` and `use_ORR` @@ -152,20 +193,22 @@ The contents of this section are shown in PDF user manual only. } \examples{ -# use PFS data from a prior study +# use single arm ORR data from a prior study run_stan( - 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) + target_hr = 0.70, + J = 2, + nevents3 = c(370, 468), + hr_bound = c(0.7790, 0.8204), + thres = 0.01, + 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=111) } \seealso{ \code{\link[stats]{Normal}} diff --git a/src/Makevars b/src/Makevars deleted file mode 100644 index 8b13789..0000000 --- a/src/Makevars +++ /dev/null @@ -1 +0,0 @@ - diff --git a/src/Makevars.win b/src/Makevars.win deleted file mode 100644 index 8b13789..0000000 --- a/src/Makevars.win +++ /dev/null @@ -1 +0,0 @@ - diff --git a/src/install.libs.R b/src/install.libs.R deleted file mode 100755 index bad5940..0000000 --- a/src/install.libs.R +++ /dev/null @@ -1,34 +0,0 @@ -libs <- file.path(R_PACKAGE_DIR, "libs", R_ARCH) -dir.create(libs, recursive = TRUE, showWarnings = FALSE) -for (file in c("symbols.rds", Sys.glob(paste0("*", SHLIB_EXT)))) { - if (file.exists(file)) { - file.copy(file, file.path(libs, file)) - } -} -inst_stan <- file.path("..", "inst", "stan") -if (dir.exists(inst_stan)) { - warning( - "Stan models in inst/stan/ are deprecated in {instantiate} ", - ">= 0.0.4.9001 (2024-01-03). Please put them in src/stan/ instead." - ) - if (file.exists("stan")) { - warning("src/stan/ already exists. Not copying models from inst/stan/.") - } else { - message("Copying inst/stan/ to src/stan/.") - fs::dir_copy(path = inst_stan, new_path = "stan") - } -} -bin <- file.path(R_PACKAGE_DIR, "bin") -if (!file.exists(bin)) { - dir.create(bin, recursive = TRUE, showWarnings = FALSE) -} -bin_stan <- file.path(bin, "stan") -fs::dir_copy(path = "stan", new_path = bin_stan) -callr::r( - func = function(bin_stan) { - for(i in list.files(bin_stan, pattern = '\\.stan$', full.names = TRUE)){ - rstan::stan_model(i, auto_write = TRUE, verbose = TRUE) - } - }, - args = list(bin_stan = bin_stan) -) diff --git a/tests/testthat/fixtures/generate_fixtures.R b/tests/testthat/fixtures/generate_fixtures.R new file mode 100644 index 0000000..3fcf3ef --- /dev/null +++ b/tests/testthat/fixtures/generate_fixtures.R @@ -0,0 +1,101 @@ +# Run this script once from the package root to regenerate Stan fixtures. +# Only needed after modifying Stan models or on a new machine/rstan version. +# Commit the generated .rds files to version control alongside the precompiled +# platform cache so CI environments can run tests without Stan installed. +# +# Usage (from package root): +# Rscript tests/testthat/fixtures/generate_fixtures.R + +devtools::load_all(quiet = TRUE) + +fixture_dir <- file.path("tests", "testthat", "fixtures") +dir.create(fixture_dir, showWarnings = FALSE, recursive = TRUE) + +rstan::rstan_options(auto_write = TRUE) + +cat("Generating Stan sampling fixtures (niter = 200)...\n") + +# --------------------------------------------------------------------------- +# phase23 model fixtures — use the same inputs as the unit tests +# --------------------------------------------------------------------------- +common_args <- list( + target_hr = 0.7, + J = 2, + nevents3 = c(370, 468), + hr_bound = c(0.779, 0.8204), + omega_mean = 0.7, + omega_var = 0.03, + est_obs_pfs = 0.53, + low_obs_pfs = 0.31, + upp_obs_pfs = 0.91, + obs_pfs_conf_level = 0.95, + thres = 0.01, + n_trt2 = 60, + n_ctrl2 = 63, + n_resp_trt2 = 33, + n_resp_ctrl2 = 18, + het_degree_p2 = "small", + het_degree_p3 = "very small", + ratio = 1, + m_0 = 0.05, + m_1 = 0.4, + nu_0 = 0.05, + nu_1 = 0.2, + lm_sd = 5, + niter = 200, + nchains = 1, + ncores = 1, + seed = 123 +) + +for (cfg in list( + list(name = "none", use_pfs = FALSE, use_orr = FALSE), + list(name = "pfs", use_pfs = TRUE, use_orr = FALSE), + list(name = "orr", use_pfs = FALSE, use_orr = TRUE), + list(name = "both", use_pfs = TRUE, use_orr = TRUE) +)) { + cat(" phase23_interim_", cfg$name, "...\n", sep = "") + suppressWarnings({ + fit <- do.call(run_stan, c(common_args, + list(use_pfs = cfg$use_pfs, + use_orr = cfg$use_orr)))$fit_rstan + }) + saveRDS(fit, file.path(fixture_dir, paste0("stanfit_", cfg$name, ".rds"))) +} + +# --------------------------------------------------------------------------- +# estimate_ctrl fixtures — two scenarios for the directional test +# --------------------------------------------------------------------------- +mod <- oncoPoS:::.load_stan_model("estimate_ctrl.stan") + +logit_low <- log(0.05 / 0.95) +logit_upp <- log(0.30 / 0.70) +mu_ctrl <- (logit_low + logit_upp) / 2 +sigma_ctrl <- (logit_upp - logit_low) / (2 * qnorm(1 - (1 - 0.8) / 2)) + +for (cfg in list( + list(name = "estimate_ctrl_high", n_resp_trt = 80L), + list(name = "estimate_ctrl_low", n_resp_trt = 40L), + list(name = "estimate_ctrl_struct", n_resp_trt = 33L) +)) { + cat(" estimate_ctrl (n_resp_trt =", cfg$n_resp_trt, ")...\n") + stan_data <- list( + n_resp_trt = cfg$n_resp_trt, + n_trt = 100L, + mu_logit_ctrl = mu_ctrl, + sigma_logit_ctrl = sigma_ctrl + ) + suppressWarnings({ + fit <- rstan::sampling(mod, data = stan_data, + iter = 200, chains = 1, seed = 123, refresh = 0) + }) + saveRDS(fit, file.path(fixture_dir, paste0("stanfit_", cfg$name, ".rds"))) +} + +cat("Done. Fixtures written to:", fixture_dir, "\n") +cat("Platform:", paste( + Sys.info()[["sysname"]], Sys.info()[["machine"]], + paste(R.version$major, R.version$minor, sep = "."), + as.character(utils::packageVersion("rstan")), + sep = "_" +), "\n") diff --git a/tests/testthat/fixtures/stanfit_both.rds b/tests/testthat/fixtures/stanfit_both.rds new file mode 100644 index 0000000..c35e790 Binary files /dev/null and b/tests/testthat/fixtures/stanfit_both.rds differ diff --git a/tests/testthat/fixtures/stanfit_estimate_ctrl_high.rds b/tests/testthat/fixtures/stanfit_estimate_ctrl_high.rds new file mode 100644 index 0000000..ccb6307 Binary files /dev/null and b/tests/testthat/fixtures/stanfit_estimate_ctrl_high.rds differ diff --git a/tests/testthat/fixtures/stanfit_estimate_ctrl_low.rds b/tests/testthat/fixtures/stanfit_estimate_ctrl_low.rds new file mode 100644 index 0000000..694f0ce Binary files /dev/null and b/tests/testthat/fixtures/stanfit_estimate_ctrl_low.rds differ diff --git a/tests/testthat/fixtures/stanfit_estimate_ctrl_struct.rds b/tests/testthat/fixtures/stanfit_estimate_ctrl_struct.rds new file mode 100644 index 0000000..90213cd Binary files /dev/null and b/tests/testthat/fixtures/stanfit_estimate_ctrl_struct.rds differ diff --git a/tests/testthat/fixtures/stanfit_none.rds b/tests/testthat/fixtures/stanfit_none.rds new file mode 100644 index 0000000..88dcb4e Binary files /dev/null and b/tests/testthat/fixtures/stanfit_none.rds differ diff --git a/tests/testthat/fixtures/stanfit_orr.rds b/tests/testthat/fixtures/stanfit_orr.rds new file mode 100644 index 0000000..df66075 Binary files /dev/null and b/tests/testthat/fixtures/stanfit_orr.rds differ diff --git a/tests/testthat/fixtures/stanfit_pfs.rds b/tests/testthat/fixtures/stanfit_pfs.rds new file mode 100644 index 0000000..64d84f7 Binary files /dev/null and b/tests/testthat/fixtures/stanfit_pfs.rds differ diff --git a/tests/testthat/helper-mock-stan.R b/tests/testthat/helper-mock-stan.R new file mode 100644 index 0000000..aba29a6 --- /dev/null +++ b/tests/testthat/helper-mock-stan.R @@ -0,0 +1,58 @@ +# --------------------------------------------------------------------------- +# Stan sampling mock — analogous to VCR cassettes in API testing. +# +# Fixtures are pre-generated stanfit objects (200 iter each) stored in +# tests/testthat/fixtures/. They contain real MCMC draws so all downstream +# consumers (rstan::extract, tidybayes::spread_draws, bayesplot) work without +# change. +# +# Routing: the mock inspects the `data` list passed to rstan::sampling / +# rstan::stan to determine which cassette to load: +# +# n_resp_trt key present → estimate_ctrl model (keyed by n_resp_trt value) +# theta_hat + orr_hat → phase23_interim_both +# theta_hat only → phase23_interim_pfs +# orr_hat only → phase23_interim_orr +# (none of the above) → phase23_interim_none +# +# Usage inside a test_that() block: +# +# local_mocked_bindings( +# sampling = mock_rstan_sampling, +# stan = mock_rstan_stan, +# .package = "rstan" +# ) +# --------------------------------------------------------------------------- + +.stan_fixture_path <- function(name) { + testthat::test_path("fixtures", paste0("stanfit_", name, ".rds")) +} + +.load_stan_fixture <- function(name) { + readRDS(.stan_fixture_path(name)) +} + +.route_stan_fixture <- function(data) { + if ("n_resp_trt" %in% names(data)) { + if (data$n_resp_trt >= 60L) return(.load_stan_fixture("estimate_ctrl_high")) + if (data$n_resp_trt >= 40L) return(.load_stan_fixture("estimate_ctrl_low")) + return(.load_stan_fixture("estimate_ctrl_struct")) + } + if ("theta_hat" %in% names(data) && "orr_hat" %in% names(data)) + return(.load_stan_fixture("both")) + if ("theta_hat" %in% names(data)) + return(.load_stan_fixture("pfs")) + if ("orr_hat" %in% names(data)) + return(.load_stan_fixture("orr")) + .load_stan_fixture("none") +} + +# Drop-in replacement for rstan::sampling (used when precompiled model exists) +mock_rstan_sampling <- function(object, data, iter, chains, cores, seed, ...) { + .route_stan_fixture(data) +} + +# Drop-in replacement for rstan::stan (used as fallback when no precompiled model) +mock_rstan_stan <- function(file, data, iter, chains, cores, seed, ...) { + .route_stan_fixture(data) +} diff --git a/tests/testthat/test-independent-testing-gen_pos.R b/tests/testthat/test-independent-testing-gen_pos.R index 01a4393..7d60d88 100644 --- a/tests/testthat/test-independent-testing-gen_pos.R +++ b/tests/testthat/test-independent-testing-gen_pos.R @@ -1,9 +1,15 @@ test_that("check results of gen_pos", { + local_mocked_bindings( + sampling = mock_rstan_sampling, + stan = mock_rstan_stan, + .package = "rstan" + ) target_hr = 0.7 J = 2 nevents3 = c(370, 468) hr_bound = c(0.779, 0.8204) - omega = 0.7 + omega_mean = 0.7 + omega_var = 0.03 est_obs_pfs = 0.53 low_obs_pfs = 0.31 upp_obs_pfs = 0.91 @@ -21,136 +27,140 @@ test_that("check results of gen_pos", { nu_0 = 0.05 nu_1 = 0.2 lm_sd = 5 - niter = 2000 + niter = 3000 nchains = 1 ncores = 1 seed = 123 re_NoPFS_NoORR = gen_pos( - target_hr, - J, - nevents3, - hr_bound, - omega, - est_obs_pfs, - low_obs_pfs, - upp_obs_pfs, - obs_pfs_conf_level, - thres, - n_trt2, - n_ctrl2, - n_resp_trt2, - n_resp_ctrl2, + target_hr = target_hr, + J = J, + nevents3 = nevents3, + hr_bound = hr_bound, + 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, + obs_pfs_conf_level = obs_pfs_conf_level, + thres = thres, + n_trt2 = n_trt2, + n_ctrl2 = n_ctrl2, + n_resp_trt2 = n_resp_trt2, + n_resp_ctrl2 = n_resp_ctrl2, use_orr = FALSE, use_pfs = FALSE, - het_degree_p2 = "small", - het_degree_p3 = "very small", - ratio, - m_0, - m_1, - nu_0, - nu_1, - lm_sd, - niter, - nchains, - ncores, - seed, + het_degree_p2 = het_degree_p2, + het_degree_p3 = het_degree_p3, + ratio = ratio, + m_0 = m_0, + m_1 = m_1, + nu_0 = nu_0, + nu_1 = nu_1, + lm_sd = lm_sd, + niter = niter, + nchains = nchains, + ncores = ncores, + seed = seed, plots_out = TRUE ) re_PFS_NoORR = gen_pos( - target_hr, - J, - nevents3, - hr_bound, - omega, - est_obs_pfs, - low_obs_pfs, - upp_obs_pfs, - obs_pfs_conf_level, - thres, - n_trt2, - n_ctrl2, - n_resp_trt2, - n_resp_ctrl2, + target_hr = target_hr, + J = J, + nevents3 = nevents3, + hr_bound = hr_bound, + 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, + obs_pfs_conf_level = obs_pfs_conf_level, + thres = thres, + n_trt2 = n_trt2, + n_ctrl2 = n_ctrl2, + n_resp_trt2 = n_resp_trt2, + n_resp_ctrl2 = n_resp_ctrl2, use_orr = FALSE, use_pfs = TRUE, - het_degree_p2 = "small", - het_degree_p3 = "very small", - ratio, - m_0, - m_1, - nu_0, - nu_1, - lm_sd, - niter, - nchains, - ncores, - seed, + het_degree_p2 = het_degree_p2, + het_degree_p3 = het_degree_p3, + ratio = ratio, + m_0 = m_0, + m_1 = m_1, + nu_0 = nu_0, + nu_1 = nu_1, + lm_sd = lm_sd, + niter = niter, + nchains = nchains, + ncores = ncores, + seed = seed, plots_out = TRUE ) re_NoPFS_ORR = gen_pos( - target_hr, - J, - nevents3, - hr_bound, - omega, - est_obs_pfs, - low_obs_pfs, - upp_obs_pfs, - obs_pfs_conf_level, - thres, - n_trt2, - n_ctrl2, - n_resp_trt2, - n_resp_ctrl2, + target_hr = target_hr, + J = J, + nevents3 = nevents3, + hr_bound = hr_bound, + 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, + obs_pfs_conf_level = obs_pfs_conf_level, + thres = thres, + n_trt2 = n_trt2, + n_ctrl2 = n_ctrl2, + n_resp_trt2 = n_resp_trt2, + n_resp_ctrl2 = n_resp_ctrl2, use_orr = TRUE, use_pfs = FALSE, - het_degree_p2 = "small", - het_degree_p3 = "very small", - ratio, - m_0, - m_1, - nu_0, - nu_1, - lm_sd, - niter, - nchains, - ncores, - seed, + het_degree_p2 = het_degree_p2, + het_degree_p3 = het_degree_p3, + ratio = ratio, + m_0 = m_0, + m_1 = m_1, + nu_0 = nu_0, + nu_1 = nu_1, + lm_sd = lm_sd, + niter = niter, + nchains = nchains, + ncores = ncores, + seed = seed, plots_out = TRUE ) re_PFS_ORR = gen_pos( - target_hr, - J, - nevents3, - hr_bound, - omega, - est_obs_pfs, - low_obs_pfs, - upp_obs_pfs, - obs_pfs_conf_level, - thres, - n_trt2, - n_ctrl2, - n_resp_trt2, - n_resp_ctrl2, + target_hr = target_hr, + J = J, + nevents3 = nevents3, + hr_bound = hr_bound, + 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, + obs_pfs_conf_level = obs_pfs_conf_level, + thres = thres, + n_trt2 = n_trt2, + n_ctrl2 = n_ctrl2, + n_resp_trt2 = n_resp_trt2, + n_resp_ctrl2 = n_resp_ctrl2, use_orr = TRUE, use_pfs = TRUE, - het_degree_p2 = "small", - het_degree_p3 = "very small", - ratio, - m_0, - m_1, - nu_0, - nu_1, - lm_sd, - niter, - nchains, - ncores, - seed, + het_degree_p2 = het_degree_p2, + het_degree_p3 = het_degree_p3, + ratio = ratio, + m_0 = m_0, + m_1 = m_1, + nu_0 = nu_0, + nu_1 = nu_1, + lm_sd = lm_sd, + niter = niter, + nchains = nchains, + ncores = ncores, + seed = seed, plots_out = TRUE ) diff --git a/tests/testthat/test-independent-testing-get_beta_params.R b/tests/testthat/test-independent-testing-get_beta_params.R new file mode 100644 index 0000000..03122e4 --- /dev/null +++ b/tests/testthat/test-independent-testing-get_beta_params.R @@ -0,0 +1,71 @@ +test_that("check the output equals to expected values given different mean and var", { + mean1 = 0.52 + var1 = 0.02 + + mean2 = 0.3 + var2 = 0.05 + + re1 = get_beta_params(mean1, var1) + re2 = get_beta_params(mean2, var2) + + # compare to manually computed values + alpha1 = mean1 * ((mean1 * (1 - mean1)) / var1 - 1) + beta1 = (1 - mean1) * ((mean1 * (1 - mean1)) / var1 - 1) + + alpha2 = mean2 * ((mean2 * (1 - mean2)) / var2 - 1) + beta2 = (1 - mean2) * ((mean2 * (1 - mean2)) / var2 - 1) + + expect_equal(re1$alpha, alpha1) + expect_equal(re1$beta, beta1) + + expect_equal(re2$alpha, alpha2) + expect_equal(re2$beta, beta2) +}) + +test_that("check the output list contains alpha and beta", { + re = get_beta_params(mean = 0.52, var = 0.02) + + expect_true(all(c("alpha", "beta") %in% names(re))) + expect_equal(length(re), 2) +}) + +test_that("check alpha and beta are both positive for valid inputs", { + re1 = get_beta_params(mean = 0.5, var = 0.01) + re2 = get_beta_params(mean = 0.2, var = 0.05) + re3 = get_beta_params(mean = 0.8, var = 0.03) + + expect_gt(re1$alpha, 0) + expect_gt(re1$beta, 0) + + expect_gt(re2$alpha, 0) + expect_gt(re2$beta, 0) + + expect_gt(re3$alpha, 0) + expect_gt(re3$beta, 0) +}) + +test_that("check the relationship between alpha and beta reflects the mean", { + # when mean > 0.5, alpha should be greater than beta + re_high = get_beta_params(mean = 0.7, var = 0.02) + expect_gt(re_high$alpha, re_high$beta) + + # when mean < 0.5, alpha should be less than beta + re_low = get_beta_params(mean = 0.3, var = 0.02) + expect_lt(re_low$alpha, re_low$beta) + + # when mean = 0.5, alpha should equal beta + re_mid = get_beta_params(mean = 0.5, var = 0.02) + expect_equal(re_mid$alpha, re_mid$beta) +}) + + +test_that("check that alpha / (alpha + beta) recovers the original mean", { + mean1 = 0.52; var1 = 0.02 + mean2 = 0.35; var2 = 0.04 + + re1 = get_beta_params(mean1, var1) + re2 = get_beta_params(mean2, var2) + + expect_equal(re1$alpha / (re1$alpha + re1$beta), mean1) + expect_equal(re2$alpha / (re2$alpha + re2$beta), mean2) +}) diff --git a/tests/testthat/test-independent-testing-resp2oddsratio_estimate_ctrl.R b/tests/testthat/test-independent-testing-resp2oddsratio_estimate_ctrl.R new file mode 100644 index 0000000..153fef2 --- /dev/null +++ b/tests/testthat/test-independent-testing-resp2oddsratio_estimate_ctrl.R @@ -0,0 +1,96 @@ +test_that("resp2oddsratio_estimate_ctrl() produces a smaller mean odds ratio when the observed RR in trt is higher", { + local_mocked_bindings( + sampling = mock_rstan_sampling, + stan = mock_rstan_stan, + .package = "rstan" + ) + n_resp_trt1 <- 80 + n_resp_trt2 <- 40 + n_trt <- 100 + low_soc_rr <- 0.05 + upp_soc_rr <- 0.30 + + log_or_high <- resp2oddsratio_estimate_ctrl( + n_resp_trt = n_resp_trt1, + n_trt = n_trt, + low_soc_rr = low_soc_rr, + upp_soc_rr = upp_soc_rr, + ci_rr = 0.80, + ncores = 1 + ) + + log_or_low <- resp2oddsratio_estimate_ctrl( + n_resp_trt = n_resp_trt2, + n_trt = n_trt, + low_soc_rr = low_soc_rr, + upp_soc_rr = upp_soc_rr, + ci_rr = 0.80, + ncores = 1 + ) + + expect_lt(exp(log_or_high$est), exp(log_or_low$est)) + +}) + +test_that("check the prior computation of mu and sigma on logit scale", { + low_soc_rr = 0.15 + upp_soc_rr = 0.35 + ci_rr = 0.8 + + logit_low = log(low_soc_rr / (1 - low_soc_rr)) + logit_upp = log(upp_soc_rr / (1 - upp_soc_rr)) + + mu_expected = (logit_low + logit_upp) / 2 + sigma_expected = (logit_upp - logit_low) / (2 * qnorm(1 - (1 - ci_rr) / 2)) + + # mu should be the midpoint of logit-transformed bounds + expect_equal(mu_expected, (logit_low + logit_upp) / 2) + + # sigma should be positive + expect_gt(sigma_expected, 0) + + # higher ci_rr => narrower normal quantile => larger sigma + sigma_narrow_ci = (logit_upp - logit_low) / (2 * qnorm(1 - (1 - 0.5) / 2)) + sigma_wide_ci = (logit_upp - logit_low) / (2 * qnorm(1 - (1 - 0.95) / 2)) + expect_gt(sigma_narrow_ci, sigma_wide_ci) +}) + +test_that("check the output structure of resp2oddsratio_estimate_ctrl", { + local_mocked_bindings( + sampling = mock_rstan_sampling, + stan = mock_rstan_stan, + .package = "rstan" + ) + n_resp_trt = 33 + n_trt = 60 + low_soc_rr = 0.15 + upp_soc_rr = 0.35 + ci_rr = 0.8 + niter = 2000 + nchains = 1 + seed = 123 + + re = resp2oddsratio_estimate_ctrl( + n_resp_trt = n_resp_trt, + n_trt = n_trt, + low_soc_rr = low_soc_rr, + upp_soc_rr = upp_soc_rr, + ci_rr = ci_rr, + niter = niter, + nchains = nchains, + seed = seed + ) + + # check output contains est and se + expect_true(all(c("est", "se") %in% names(re))) + expect_equal(length(re), 2) + + # check est and se are numeric scalars + expect_true(is.numeric(re$est)) + expect_true(is.numeric(re$se)) + expect_equal(length(re$est), 1) + expect_equal(length(re$se), 1) + + # se must be positive + expect_gt(re$se, 0) +}) diff --git a/tests/testthat/test-independent-testing-run_stan.R b/tests/testthat/test-independent-testing-run_stan.R index 1304140..c3d0adc 100644 --- a/tests/testthat/test-independent-testing-run_stan.R +++ b/tests/testthat/test-independent-testing-run_stan.R @@ -1,10 +1,15 @@ test_that("Test the stan_list and stan_file when different combinations of endpoints are used", { - # use KN189 (phase 3) and KN021-G as test example: + local_mocked_bindings( + sampling = mock_rstan_sampling, + stan = mock_rstan_stan, + .package = "rstan" + ) target_hr = 0.7 J = 2 nevents3 = c(370, 468) hr_bound = c(0.779, 0.8204) - omega = 0.7 + omega_mean = 0.7 + omega_var = 0.03 est_obs_pfs = 0.53 low_obs_pfs = 0.31 upp_obs_pfs = 0.91 @@ -22,126 +27,138 @@ test_that("Test the stan_list and stan_file when different combinations of endpo nu_0 = 0.05 nu_1 = 0.2 lm_sd = 5 - niter = 2000 + niter = 3000 nchains = 1 ncores = 1 seed = 123 NoPFS_NoORR <- run_stan( - target_hr, - J, - nevents3, - hr_bound, omega, - est_obs_pfs, - low_obs_pfs, - upp_obs_pfs, - obs_pfs_conf_level, - thres, - n_trt2, - n_ctrl2, - n_resp_trt2, - n_resp_ctrl2, + target_hr = target_hr, + J = J, + nevents3 = nevents3, + hr_bound = hr_bound, + 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, + obs_pfs_conf_level = obs_pfs_conf_level, + thres = thres, + n_trt2 = n_trt2, + n_ctrl2 = n_ctrl2, + n_resp_trt2 = n_resp_trt2, + n_resp_ctrl2 = n_resp_ctrl2, use_orr = FALSE, use_pfs = FALSE, - het_degree_p2, - het_degree_p3, - ratio, - m_0, - m_1, - nu_0, - nu_1, - lm_sd, - niter, - nchains, - ncores, - seed) + het_degree_p2 = het_degree_p2, + het_degree_p3 = het_degree_p3, + ratio = ratio, + m_0 = m_0, + m_1 = m_1, + nu_0 = nu_0, + nu_1 = nu_1, + lm_sd = lm_sd, + niter = niter, + nchains = nchains, + ncores = ncores, + seed = seed + ) PFS_NoORR <- run_stan( - target_hr, - J, - nevents3, - hr_bound, omega, - est_obs_pfs, - low_obs_pfs, - upp_obs_pfs, - obs_pfs_conf_level, - thres, - n_trt2, - n_ctrl2, - n_resp_trt2, - n_resp_ctrl2, + target_hr = target_hr, + J = J, + nevents3 = nevents3, + hr_bound = hr_bound, + 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, + obs_pfs_conf_level = obs_pfs_conf_level, + thres = thres, + n_trt2 = n_trt2, + n_ctrl2 = n_ctrl2, + n_resp_trt2 = n_resp_trt2, + n_resp_ctrl2 = n_resp_ctrl2, use_orr = FALSE, use_pfs = TRUE, - het_degree_p2, - het_degree_p3, - ratio, - m_0, - m_1, - nu_0, - nu_1, - lm_sd, - niter, - nchains, - ncores, - seed) + het_degree_p2 = het_degree_p2, + het_degree_p3 = het_degree_p3, + ratio = ratio, + m_0 = m_0, + m_1 = m_1, + nu_0 = nu_0, + nu_1 = nu_1, + lm_sd = lm_sd, + niter = niter, + nchains = nchains, + ncores = ncores, + seed = seed + ) NoPFS_ORR <- run_stan( - target_hr, - J, - nevents3, - hr_bound, omega, - est_obs_pfs, - low_obs_pfs, - upp_obs_pfs, - obs_pfs_conf_level, - thres, - n_trt2, - n_ctrl2, - n_resp_trt2, - n_resp_ctrl2, + target_hr = target_hr, + J = J, + nevents3 = nevents3, + hr_bound = hr_bound, + 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, + obs_pfs_conf_level = obs_pfs_conf_level, + thres = thres, + n_trt2 = n_trt2, + n_ctrl2 = n_ctrl2, + n_resp_trt2 = n_resp_trt2, + n_resp_ctrl2 = n_resp_ctrl2, use_orr = TRUE, use_pfs = FALSE, - het_degree_p2, - het_degree_p3, - ratio, - m_0, - m_1, - nu_0, - nu_1, - lm_sd, - niter, - nchains, - ncores, - seed) + het_degree_p2 = het_degree_p2, + het_degree_p3 = het_degree_p3, + ratio = ratio, + m_0 = m_0, + m_1 = m_1, + nu_0 = nu_0, + nu_1 = nu_1, + lm_sd = lm_sd, + niter = niter, + nchains = nchains, + ncores = ncores, + seed = seed + ) PFS_ORR <- run_stan( - target_hr, - J, - nevents3, - hr_bound, omega, - est_obs_pfs, - low_obs_pfs, - upp_obs_pfs, - obs_pfs_conf_level, - thres, - n_trt2, - n_ctrl2, - n_resp_trt2, - n_resp_ctrl2, + target_hr = target_hr, + J = J, + nevents3 = nevents3, + hr_bound = hr_bound, + 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, + obs_pfs_conf_level = obs_pfs_conf_level, + thres = thres, + n_trt2 = n_trt2, + n_ctrl2 = n_ctrl2, + n_resp_trt2 = n_resp_trt2, + n_resp_ctrl2 = n_resp_ctrl2, use_orr = TRUE, use_pfs = TRUE, - het_degree_p2, - het_degree_p3, - ratio, - m_0, - m_1, - nu_0, - nu_1, - lm_sd, - niter, - nchains, - ncores, - seed) + het_degree_p2 = het_degree_p2, + het_degree_p3 = het_degree_p3, + ratio = ratio, + m_0 = m_0, + m_1 = m_1, + nu_0 = nu_0, + nu_1 = nu_1, + lm_sd = lm_sd, + niter = niter, + nchains = nchains, + ncores = ncores, + seed = seed + ) expect_equal(NoPFS_NoORR$stan_file, "phase23_interim_none.stan") expect_equal(PFS_NoORR$stan_file, "phase23_interim_pfs.stan") @@ -163,14 +180,14 @@ test_that("Test the stan_list and stan_file when different combinations of endpo expect_equal(NoPFS_ORR$stan_list$wls_sd, PFS_ORR$stan_list$wls_sd) # check stan_list items - expect_true(all(c("omega","delta_P","sigma_P1","sigma_P2","tau_sd2","tau_sd3","J","Sigma") + expect_true(all(c("omega_alpha","omega_beta","delta_P","sigma_P1","sigma_P2","tau_sd2","tau_sd3","J","Sigma") %in% names(NoPFS_NoORR$stan_list))) - expect_true(all(c("omega","delta_P","sigma_P1","sigma_P2","tau_sd2","tau_sd3","J","Sigma", + expect_true(all(c("omega_alpha","omega_beta","delta_P","sigma_P1","sigma_P2","tau_sd2","tau_sd3","J","Sigma", "theta_hat","theta_hat_sd") %in% names(PFS_NoORR$stan_list))) - expect_true(all(c("omega","delta_P","sigma_P1","sigma_P2","tau_sd2","tau_sd3","J","Sigma", + expect_true(all(c("omega_alpha","omega_beta","delta_P","sigma_P1","sigma_P2","tau_sd2","tau_sd3","J","Sigma", "orr_hat","orr_hat_sd","m_0","m_1","nu_0","nu_1","wls_sd") %in% names(NoPFS_ORR$stan_list))) - expect_true(all(c("omega","delta_P","sigma_P1","sigma_P2","tau_sd2","tau_sd3","J","Sigma", + expect_true(all(c("omega_alpha","omega_beta","delta_P","sigma_P1","sigma_P2","tau_sd2","tau_sd3","J","Sigma", "theta_hat","theta_hat_sd", "orr_hat","orr_hat_sd","m_0","m_1","nu_0", "nu_1","wls_sd") %in% names(PFS_ORR$stan_list))) }) diff --git a/vignettes/oncoPoS.Rmd b/vignettes/oncoPoS.Rmd index c888a35..af8e879 100644 --- a/vignettes/oncoPoS.Rmd +++ b/vignettes/oncoPoS.Rmd @@ -3,8 +3,6 @@ title: "Introduction to oncoPoS" output: rmarkdown::html_vignette bibliography: "oncoPoS.bib" -resource_files: - - package_manual.pdf vignette: > %\VignetteIndexEntry{Introduction to oncoPoS} %\VignetteEngine{knitr::rmarkdown} @@ -16,6 +14,8 @@ editor_options: + + # Overview `oncoPoS` performs Probability of Success (PoS) calculations for a phase 3 @@ -29,9 +29,9 @@ vignette for details. The main functions of the `oncoPoS` are `run_stan` and `gen_pos`, which generate the PoS estimate and its standard error (SE) at each planned analysis by -running `rstan`. Note that the SE is +running `rstan`. The SE is calculated as $\sqrt{\hat{p}(1-\hat{p})/N_{\text{draws}}}$, where $\hat{p}$ is the estimated PoS and $N_{\text{draws}}$ is the total number of post-warmup MCMC draws. -The PoS estimate could be generated without considering an observed treatment effect in an earlier study, +The PoS estimate can be generated without considering an observed treatment effect from an earlier study, or it can incorporate prior ORR and/or PFS data as described in the sections below. # PoS estimate example @@ -48,9 +48,13 @@ While the overall alpha level of 2.5% and power 95.5% could be seen as design features in general, these are not considered as direct design features for the PoS estimation and therefore are not inputs in `oncoPoS` relevant functions. -To proceed with the PoS estimation, an industry benchmark is required as well. -We'll assume that, without considering specific design features of this trial, -there is a 50% chance that such trial will be successful. +To proceed with the PoS estimation, an industry benchmark is required as +well. We treat the benchmark as a random variable $\omega$, reflecting +uncertainty in the prior probability of trial success. Without +considering the specific design features of this trial, we assume there +is a 50% chance that such a trial will be successful, i.e., +$\omega_{\text{mean}} = 0.5$, with a default variance of +$\omega_{\text{var}} = 0.02$. ## No earlier study data @@ -61,24 +65,20 @@ each analysis is calculated as following: ``` r -library(oncoPoS) - gen_pos( - target_hr = 0.7, - J = 2, - nevents3 = c(370, 468), - hr_bound = c(0.779, 0.8204), - omega = 0.5, - seed = 245 - ) -``` - -``` -## # A tibble: 2 × 3 -## J pos pos_se -## -## 1 1 0.414 0.0110 -## 2 2 0.492 0.0112 + target_hr = 0.7, + J = 2, + nevents3 = c(370, 468), + hr_bound = c(0.779, 0.8204), + omega_mean = 0.5, + seed = 245, + refresh = 0 +) +#> # A tibble: 2 × 5 +#> J pos pos_se omega_mean omega_var +#> +#> 1 1 0.416 0.0220 0.500 0.0208 +#> 2 2 0.474 0.0223 0.500 0.0208 ``` ## Prior PFS data @@ -94,25 +94,23 @@ the PoS estimates as follows: ``` r gen_pos( - target_hr = 0.7, - J = 2, - nevents3 = c(370, 468), - hr_bound = c(0.779, 0.8204), - omega = 0.5, - use_pfs = TRUE, - est_obs_pfs = 0.53, - low_obs_pfs = 0.31, + target_hr = 0.7, + J = 2, + nevents3 = c(370, 468), + hr_bound = c(0.779, 0.8204), + omega_mean = 0.5, + use_pfs = TRUE, + est_obs_pfs = 0.53, + low_obs_pfs = 0.31, upp_obs_pfs = 0.91, - seed = 245 - ) -``` - -``` -## # A tibble: 2 × 3 -## J pos pos_se -## -## 1 1 0.73 0.00993 -## 2 2 0.812 0.00874 + seed = 245, + refresh = 0 +) +#> # A tibble: 2 × 5 +#> J pos pos_se omega_mean omega_var +#> +#> 1 1 0.732 0.0198 0.515 0.0173 +#> 2 2 0.816 0.0173 0.515 0.0173 ``` ## Prior ORR data @@ -122,36 +120,89 @@ available, it can be used for PoS estimation as well. We will assume that, based on the earlier study 33 out of 60 and 18 out of 63 participants in the experimental and control arms had responses respectively. -In this case, when the ORR data is used, coefficients for the linear -relationship between the log treatment effect on ORR and PFS must be specified: +In this case, when the ORR data is used, the coefficients for the linear +relationship between the log treatment effect on ORR and PFS are estimated +indication-specific using prior meta-analytic data through a Bayesian hierarchical model. +If no indication is specified, the default coefficients (calculated as the mean across +all indication groups) will be used. + +### Indication Groups + +- **Group 1: Hematologic malignancies** + Includes classical Hodgkin lymphoma (cHL), diffuse large B-cell lymphoma (DLBCL), + follicular lymphoma (FL), multiple myeloma (MM), non-Hodgkin lymphoma (NHL), + and peripheral T-cell lymphoma (PTCL). +- **Group 2: Gynecologic cancers** + Includes cervical, endometrial, and ovarian cancers. +- **Group 3: Thoracic malignancies** + Includes non-small cell lung cancer (NSCLC), small cell lung cancer (SCLC), + and mesothelioma. +- **Group 4: Urologic and gastrointestinal solid tumors** + Includes bladder cancer, gastric cancer, and renal cell carcinoma (RCC). +- **Group 5: Breast cancer** + Includes breast cancer. +Here we specify the indication group as breast cancer (Group 5): ``` r gen_pos( - target_hr = 0.7, - J = 2, - nevents3 = c(370, 468), - hr_bound = c(0.779, 0.8204), - omega = 0.5, - use_orr = TRUE, - n_resp_trt2 = 33, - n_trt2 = 60, - n_resp_ctrl2 = 18, - n_ctrl2 = 63, - m_0 = -0.20, - m_1 = 2.2, - nu_0 = 0.0541, - nu_1 = 0.2145, - lm_sd = 5.1158, - seed = 245 - ) + target_hr = 0.7, + J = 2, + nevents3 = c(370, 468), + hr_bound = c(0.779, 0.8204), + omega_mean = 0.5, + use_orr = TRUE, + n_resp_trt2 = 33, + n_trt2 = 60, + n_resp_ctrl2 = 18, + n_ctrl2 = 63, + indication = 5, + seed = 245, + refresh = 0 +) +#> # A tibble: 2 × 5 +#> J pos pos_se omega_mean omega_var +#> +#> 1 1 0.86 0.0155 0.533 0.0187 +#> 2 2 0.898 0.0135 0.533 0.0187 ``` -``` -## # A tibble: 2 × 3 -## J pos pos_se -## -## 1 1 0.572 0.0111 -## 2 2 0.656 0.0106 +### Single-arm setting + +When no concurrent control arm is available in the earlier study, the SOC +response rate $p_{\text{SOC}}$ cannot be observed directly. Instead, a +distribution is placed over $p_{\text{SOC}}$ using user-specified lower +and upper bounds on the control ORR: + +- `low_soc_rr`: lower bound for the control ORR +- `upp_soc_rr`: upper bound for the control ORR + +In this setting, `n_resp_ctrl2` and `n_ctrl2` are replaced by `low_soc_rr` +and `upp_soc_rr`. For example, assuming the SOC ORR is bounded between 5% +and 20%: + + +``` r +gen_pos( + target_hr = 0.7, + J = 2, + nevents3 = c(370, 468), + hr_bound = c(0.779, 0.8204), + omega_mean = 0.5, + use_orr = TRUE, + n_resp_trt2 = 33, + n_trt2 = 60, + single_arm = TRUE, + low_soc_rr = 0.05, + upp_soc_rr = 0.20, + indication = 5, + seed = 245, + refresh = 0 +) +#> # A tibble: 2 × 5 +#> J pos pos_se omega_mean omega_var +#> +#> 1 1 0.864 0.0153 0.539 0.0206 +#> 2 2 0.91 0.0128 0.539 0.0206 ``` diff --git a/vignettes/oncoPoS.Rmd.orig b/vignettes/oncoPoS.Rmd.orig index 4871e27..0f47c48 100644 --- a/vignettes/oncoPoS.Rmd.orig +++ b/vignettes/oncoPoS.Rmd.orig @@ -3,8 +3,6 @@ title: "Introduction to oncoPoS" output: rmarkdown::html_vignette bibliography: "oncoPoS.bib" -resource_files: - - package_manual.pdf vignette: > %\VignetteIndexEntry{Introduction to oncoPoS} %\VignetteEngine{knitr::rmarkdown} @@ -12,6 +10,16 @@ editor_options: chunk_output_type: console --- +```{r setup, include=FALSE} +knitr::opts_chunk$set( + collapse = TRUE, + comment = "#>", + eval = TRUE, + message = FALSE, + warning = FALSE +) +``` + @@ -29,9 +37,9 @@ vignette for details. The main functions of the `oncoPoS` are `run_stan` and `gen_pos`, which generate the PoS estimate and its standard error (SE) at each planned analysis by -running `rstan`. Note that the SE is +running `rstan`. The SE is calculated as $\sqrt{\hat{p}(1-\hat{p})/N_{\text{draws}}}$, where $\hat{p}$ is the estimated PoS and $N_{\text{draws}}$ is the total number of post-warmup MCMC draws. -The PoS estimate could be generated without considering an observed treatment effect in an earlier study, +The PoS estimate can be generated without considering an observed treatment effect from an earlier study, or it can incorporate prior ORR and/or PFS data as described in the sections below. # PoS estimate example @@ -48,9 +56,13 @@ While the overall alpha level of 2.5% and power 95.5% could be seen as design features in general, these are not considered as direct design features for the PoS estimation and therefore are not inputs in `oncoPoS` relevant functions. -To proceed with the PoS estimation, an industry benchmark is required as well. -We'll assume that, without considering specific design features of this trial, -there is a 50% chance that such trial will be successful. +To proceed with the PoS estimation, an industry benchmark is required as +well. We treat the benchmark as a random variable $\omega$, reflecting +uncertainty in the prior probability of trial success. Without +considering the specific design features of this trial, we assume there +is a 50% chance that such a trial will be successful, i.e., +$\omega_{\text{mean}} = 0.5$, with a default variance of +$\omega_{\text{var}} = 0.02$. ## No earlier study data @@ -59,17 +71,16 @@ phase 3 through an objective response rate (ORR) or PFS, the PoS estimate at each analysis is calculated as following: -```{r} -library(oncoPoS) - +```{r no-earlier-data} gen_pos( - target_hr = 0.7, - J = 2, - nevents3 = c(370, 468), - hr_bound = c(0.779, 0.8204), - omega = 0.5, - seed = 245 - ) + target_hr = 0.7, + J = 2, + nevents3 = c(370, 468), + hr_bound = c(0.779, 0.8204), + omega_mean = 0.5, + seed = 245, + refresh = 0 +) ``` ## Prior PFS data @@ -82,19 +93,20 @@ thus design a phase 3 trial. This additional information is incorporated in the PoS estimates as follows: -```{r} +```{r prior-pfs} gen_pos( - target_hr = 0.7, - J = 2, - nevents3 = c(370, 468), - hr_bound = c(0.779, 0.8204), - omega = 0.5, - use_pfs = TRUE, - est_obs_pfs = 0.53, - low_obs_pfs = 0.31, + target_hr = 0.7, + J = 2, + nevents3 = c(370, 468), + hr_bound = c(0.779, 0.8204), + omega_mean = 0.5, + use_pfs = TRUE, + est_obs_pfs = 0.53, + low_obs_pfs = 0.31, upp_obs_pfs = 0.91, - seed = 245 - ) + seed = 245, + refresh = 0 +) ``` ## Prior ORR data @@ -104,27 +116,77 @@ available, it can be used for PoS estimation as well. We will assume that, based on the earlier study 33 out of 60 and 18 out of 63 participants in the experimental and control arms had responses respectively. -In this case, when the ORR data is used, coefficients for the linear -relationship between the log treatment effect on ORR and PFS must be specified: +In this case, when the ORR data is used, the coefficients for the linear +relationship between the log treatment effect on ORR and PFS are estimated +indication-specific using prior meta-analytic data through a Bayesian hierarchical model. +If no indication is specified, the default coefficients (calculated as the mean across +all indication groups) will be used. + +### Indication Groups + +- **Group 1: Hematologic malignancies** + Includes classical Hodgkin lymphoma (cHL), diffuse large B-cell lymphoma (DLBCL), + follicular lymphoma (FL), multiple myeloma (MM), non-Hodgkin lymphoma (NHL), + and peripheral T-cell lymphoma (PTCL). +- **Group 2: Gynecologic cancers** + Includes cervical, endometrial, and ovarian cancers. +- **Group 3: Thoracic malignancies** + Includes non-small cell lung cancer (NSCLC), small cell lung cancer (SCLC), + and mesothelioma. +- **Group 4: Urologic and gastrointestinal solid tumors** + Includes bladder cancer, gastric cancer, and renal cell carcinoma (RCC). +- **Group 5: Breast cancer** + Includes breast cancer. + +Here we specify the indication group as breast cancer (Group 5): + +```{r prior-orr} +gen_pos( + target_hr = 0.7, + J = 2, + nevents3 = c(370, 468), + hr_bound = c(0.779, 0.8204), + omega_mean = 0.5, + use_orr = TRUE, + n_resp_trt2 = 33, + n_trt2 = 60, + n_resp_ctrl2 = 18, + n_ctrl2 = 63, + indication = 5, + seed = 245, + refresh = 0 +) +``` + +### Single-arm setting + +When no concurrent control arm is available in the earlier study, the SOC +response rate $p_{\text{SOC}}$ cannot be observed directly. Instead, a +distribution is placed over $p_{\text{SOC}}$ using user-specified lower +and upper bounds on the control ORR: + +- `low_soc_rr`: lower bound for the control ORR +- `upp_soc_rr`: upper bound for the control ORR +In this setting, `n_resp_ctrl2` and `n_ctrl2` are replaced by `low_soc_rr` +and `upp_soc_rr`. For example, assuming the SOC ORR is bounded between 5% +and 20%: -```{r} +```{r prior-orr-single-arm} gen_pos( - target_hr = 0.7, - J = 2, - nevents3 = c(370, 468), - hr_bound = c(0.779, 0.8204), - omega = 0.5, - use_orr = TRUE, + target_hr = 0.7, + J = 2, + nevents3 = c(370, 468), + hr_bound = c(0.779, 0.8204), + omega_mean = 0.5, + use_orr = TRUE, n_resp_trt2 = 33, - n_trt2 = 60, - n_resp_ctrl2 = 18, - n_ctrl2 = 63, - m_0 = -0.20, - m_1 = 2.2, - nu_0 = 0.0541, - nu_1 = 0.2145, - lm_sd = 5.1158, - seed = 245 - ) + n_trt2 = 60, + single_arm = TRUE, + low_soc_rr = 0.05, + upp_soc_rr = 0.20, + indication = 5, + seed = 245, + refresh = 0 +) ``` diff --git a/vignettes/oncoPoS.bib b/vignettes/oncoPoS.bib index 828c6c3..8b52844 100644 --- a/vignettes/oncoPoS.bib +++ b/vignettes/oncoPoS.bib @@ -18,4 +18,15 @@ @article{blumenthal2015overall pages={1008}, year={2015}, publisher={American Society of Clinical Oncology} -} \ No newline at end of file +} + +@article{weber2021applying, + title={Applying Meta-Analytic-Predictive Priors with the {R} {B}ayesian Evidence Synthesis Tools}, + author={Weber, Sebastian and Li, Yue and Seaman III, John W and Kakizume, Tomoyuki and Schmidli, Heinz}, + journal={Journal of Statistical Software}, + volume={100}, + number={19}, + pages={1--39}, + year={2021}, + doi={10.18637/jss.v100.i19} +} diff --git a/vignettes/pos_bayes_framework.Rmd b/vignettes/pos_bayes_framework.Rmd index 4d8e27a..9357f2b 100644 --- a/vignettes/pos_bayes_framework.Rmd +++ b/vignettes/pos_bayes_framework.Rmd @@ -3,8 +3,6 @@ title: "PoS Bayesian Framework for Pivotal Oncology Trials" output: rmarkdown::html_vignette bibliography: "oncoPoS.bib" -resource_files: - - package_manual.pdf vignette: > %\VignetteEngine{knitr::rmarkdown} %\VignetteIndexEntry{PoS Bayesian Framework for Pivotal Oncology Trials} @@ -18,12 +16,6 @@ knitr::opts_chunk$set( ) ``` -```{r, include=FALSE, eval = FALSE} -manual <-devtools::build_manual(path = ".") -manual <- basename(grep("--output=", manual$command, value = TRUE)) -file.rename(from = manual, to = "package_manual.pdf") -``` - # Introduction @@ -75,16 +67,21 @@ will follow Supplementary Materials E in @hampson2022improving. # Population level model \label{sec:gen_pop_level} The population level treatment effect, $\mu_P$, is assumed to come from a -mixture prior: +mixture prior with a random mixing weight: \begin{align} -\mu_{P} \sim \omega Normal (\delta_{P}, \sigma_{P1}^2) + (1-\omega) Normal(0, \sigma_{P2}^2) , +\mu_{P} & \sim \omega Normal (\delta_{P}, \sigma_{P1}^2) + (1-\omega) Normal(0, \sigma_{P2}^2),\\ +\omega &\sim \text{Beta}(\alpha, \beta), \end{align} with the following components: -- $w$ is a probability that $\mu_P$ comes from the enthusiastic prior component. -The value of $w$ is determined by the industry benchmark. +- $\omega$ is the probability that $\mu_P$ comes from the enthusiastic prior component. +$\omega$ is treated as a random variable to incorporate uncertainty and +variability in the benchmark probability of success (PoS). The Beta prior +allows integration of historical information (e.g., industry Phase III +success rates or machine-learning-based predictions) while permitting +data-driven updating. - $Normal(\delta_P, \sigma^2_{P1})$ is the enthusiastic component, i.e., a distribution which is centered at the target treatment effect $\delta_P$ (i.e., @@ -189,3 +186,155 @@ Based on the approximated correlation between $\theta_{ORR, 2}$ and $\theta_{PFS, 2}$, a distribution for $\theta_{PFS, 2}$ can be obtained and, therefore, a predicated efficacy PoS can be estimated as using models that are outlined above. + +## Indication-specific surrogate-primary endpoint relationships (e.g., ORR $\rightarrow$ PFS) {#sec:indication_specific} + +When early endpoint objective response rate (ORR) is used to predict +phase III progression-free survival (PFS), the strength and direction +of the association may vary substantially across cancer types. +To account for this heterogeneity, we group cancer indications into five categories, +each associated with a distinct set of ORR–PFS regression parameters derived from +prior Bayesian hierarchical modeling. + +Trial indexed by: $j = 1, \dots, J$. Indication indexed by: $k = z_j \in \{1, \dots, K\}$. + +$\theta^P_j$: treatment effect on endpoint $P$ at trial $j$. +$\sigma^P_j$: standard error of treatment effect on endpoint $P$ at trial $j$. + +### Level 1: Observed trial-level model + +Each trial reports an observed log-odds ratio for ORR, $\hat\theta^{\text{ORR}}_j$, +and an observed log-hazard ratio for PFS, $\hat\theta^{\text{PFS}}_j$, modeled as: + +$$ +\hat\theta^{\text{ORR}}_j \sim \mathcal{N}\!\left(\theta^{\text{ORR}}_j,\ \left(\sigma^{\text{ORR}}_j\right)^2\right) +$$ + +$$ +\hat\theta^{\text{PFS}}_j \sim \mathcal{N}\!\left(\theta^{\text{PFS}}_j,\ \left(\sigma^{\text{PFS}}_j\right)^2\right) +$$ + +Conditional on the (latent) true PFS effect, the true ORR effect follows an +indication-specific regression: + +$$ +\theta^{\text{ORR}}_j \;\sim\; \mathcal{N}\!\left(\alpha_{z_j} + \beta_{z_j}\,\theta^{\text{PFS}}_j,\; +\frac{\sigma_{\text{WLS}}^{2}}{N_{\text{patients},j}}\right) +$$ + +### Level 2: Indication-specific regression parameters + +Information is shared across indications via hierarchical modeling. +Each pair $(\alpha_k, \beta_k)$ is learned adaptively across indications: + +$$ +\alpha_k \sim \mathcal{N}(a_0,\ s_0^2), \qquad +\beta_k \sim \mathcal{N}(b_0,\ s_1^2), \qquad +k = 1, \dots, 5 +$$ + +### Level 3: Population-level hyperpriors + +$$ +a_0 \sim \mathcal{N}(0,\ 5^2), \quad +b_0 \sim \mathcal{N}(2,\ 5^2), \quad +s_0 \sim \mathcal{N}(0,\ 5^2), \quad +s_1 \sim \mathcal{N}(0,\ 5^2), \quad +\sigma_{\text{WLS}} \sim \mathcal{N}(0,\ 5^2) +$$ + +### Indication Groups + +**Group 1: Hematologic malignancies** + +Includes classical Hodgkin lymphoma (cHL), diffuse large B-cell lymphoma (DLBCL), +follicular lymphoma (FL), multiple myeloma (MM), non-Hodgkin lymphoma (NHL), and peripheral T-cell lymphoma (PTCL). + +**Group 2: Gynecologic cancers** + +Includes cervical, endometrial, and ovarian cancers. + +**Group 3: Thoracic malignancies** + +Includes non-small cell lung cancer (NSCLC), small cell lung cancer (SCLC), +and mesothelioma. + +**Group 4: Urologic and gastrointestinal solid tumors** + +Includes bladder cancer, gastric cancer, and renal cell carcinoma (RCC). + +**Group 5: Breast cancer** + +Includes breast cancer. + +> If no indication is specified, the average ORR–PFS relationship across all +> indication groups will be used. + + +## Prior ORR data {#sec:orr_data} + +The method for computing the observed log-odds ratio +$\hat\theta^{\text{ORR}}$ and its standard error differs depending on whether +the earlier study was a two-arm or single-arm trial. + +### Two-arm setting + +In the two-arm setting, $\hat\theta^{\text{ORR}}$ and its standard error are +estimated directly from the observed response counts in each arm using +frequentist estimation: + +$$ +\hat\theta^{\text{ORR}} = \log\!\left(\frac{x_{\text{SOC}}}{n_{\text{SOC}} - x_{\text{SOC}}}\right) - \log\!\left(\frac{x_{\text{trt}}}{n_{\text{trt}} - x_{\text{trt}}}\right) +$$ + +$$ +SE\!\left(\hat\theta^{\text{ORR}}\right) = \sqrt{\frac{1}{x_{\text{trt}}} + \frac{1}{n_{\text{trt}} - x_{\text{trt}}} + \frac{1}{x_{\text{SOC}}} + \frac{1}{n_{\text{SOC}} - x_{\text{SOC}}}} +$$ + +where $x_{\text{SOC}}$, $x_{\text{trt}}$ are the number of responders and +$n_{\text{SOC}}$, $n_{\text{trt}}$ are the total number of patients in the SOC +and treatment groups, respectively. + +### Single-arm setting + +When only single-arm data are available (i.e., no concurrent control +arm), the observed ORR treatment effect cannot be computed directly +from two-arm counts. Instead, following @weber2021applying, uncertainty +in the control ORR is incorporated by placing a distribution over the +Standard of Care (SOC) response rate, $p_{\text{SOC}}$, using +user-specified lower and upper bounds. + +Let `low_soc_rr` and `upp_soc_rr` denote the lower and upper bounds for +the control ORR. The corresponding logit-scale bounds are: + +$$ +\texttt{logit_low} = \log\!\left(\frac{\texttt{low_soc_rr}}{1 - \texttt{low_soc_rr}}\right), \qquad +\texttt{logit_upp} = \log\!\left(\frac{\texttt{upp_soc_rr}}{1 - \texttt{upp_soc_rr}}\right) +$$ + +Rather than fixing $p_{\text{SOC}}$ to a point estimate, a normal +distribution is placed on the logit scale: + +$$ +\text{logit}(p_{\text{SOC}}) \sim \mathcal{N}\!\left(\mu_{\text{SOC}},\, \sigma_{\text{SOC}}^2\right) +\;\Rightarrow\; p_{\text{SOC}} = \frac{e^x}{1 + e^x}, \quad x \sim \mathcal{N}\!\left(\mu_{\text{SOC}},\, \sigma_{\text{SOC}}^2\right) +$$ + +where the mean and standard deviation are derived from the specified +bounds: + +$$ +\mu_{\text{SOC}} = \frac{\texttt{logit_low} + \texttt{logit_upp}}{2} +$$ + +$$ +\sigma_{\text{SOC}} = \frac{\texttt{logit_upp} - \texttt{logit_low}}{2\cdot\Phi^{-1}\!\left(1 - \dfrac{1 - \texttt{ci_rr}}{2}\right)} +$$ + +$\Phi^{-1}(\cdot)$ denotes the quantile function of the standard normal +distribution. `ci_rr` is the user-specified confidence level for the SOC +response bounds (default = 80%). Lower values may be used for more +conservative assumptions or when the SOC data are uncertain. In the single-arm setting, +$x_{\text{SOC}}$ and $n_{\text{SOC}}$ entering the expressions above are drawn +from the distribution over $p_{\text{SOC}}$ rather than observed directly. +