From 0b100246934845ee205a38b9f3423eb1d0d1bc71 Mon Sep 17 00:00:00 2001 From: Isaac Gravestock Date: Thu, 6 Feb 2025 11:51:26 +0100 Subject: [PATCH 1/6] add stan models --- inst/stan/bin/logit_fpp.stan | 45 ++++++++++++++++++++ inst/stan/cont/gauss_fpp.stan | 40 ++++++++++++++++++ inst/stan/surv/exp_fpp.stan | 52 +++++++++++++++++++++++ inst/stan/surv/pem_fpp.stan | 58 ++++++++++++++++++++++++++ inst/stan/surv/weib_ph_fpp.stan | 73 +++++++++++++++++++++++++++++++++ 5 files changed, 268 insertions(+) create mode 100644 inst/stan/bin/logit_fpp.stan create mode 100644 inst/stan/cont/gauss_fpp.stan create mode 100644 inst/stan/surv/exp_fpp.stan create mode 100644 inst/stan/surv/pem_fpp.stan create mode 100644 inst/stan/surv/weib_ph_fpp.stan diff --git a/inst/stan/bin/logit_fpp.stan b/inst/stan/bin/logit_fpp.stan new file mode 100644 index 00000000..70dce24c --- /dev/null +++ b/inst/stan/bin/logit_fpp.stan @@ -0,0 +1,45 @@ +// Logistic regression +// Fixed power prior + +data { + + int N; // number of observations + vector[N] trt; // treatment indicator + array[N] int y; // outcome + vector[N] power; // power parameter + + {{ weights.data }} + {{ cov.data }} + +} + +parameters { + + real beta_trt; // treatment effect + real alpha; // baseline log-odds + + {{ cov.parameters }} + +} + +transformed parameters { + + real OR_trt = exp(beta_trt); + +} + +model { + + vector[N] lp; + + {{ trt.prior }} + {{ cov.priors }} + {{ baseline.prior }} + + lp = alpha + trt * beta_trt {{ cov.linpred }} ; + + for (i in 1:N) { + target += bernoulli_logit_lupmf(y[i] | lp[i]) * power[i] {{ weights.likelihood }}; + } + +} diff --git a/inst/stan/cont/gauss_fpp.stan b/inst/stan/cont/gauss_fpp.stan new file mode 100644 index 00000000..648e7d84 --- /dev/null +++ b/inst/stan/cont/gauss_fpp.stan @@ -0,0 +1,40 @@ +// Continuous normal +// Fixed power prior + +data { + + int N; // number of observations + vector[N] trt; // treatment indicator + array[N] real y; // outcome + vector[N] power; // power parameter + + {{ weights.data }} + {{ cov.data }} + +} + +parameters { + + real beta_trt; // treatment effect + real alpha; // baseline mean + real std_dev_outcome; // standard deviation of outcome + + {{ cov.parameters }} + +} + +model { + vector[N] lp; + + {{ trt.prior }} + {{ cov.priors }} + {{ baseline.prior }} + {{ stdev.prior }} + + lp = alpha + trt * beta_trt {{ cov.linpred }} ; + + for (i in 1:N) { + target += normal_lupdf(y[i] | lp[i], std_dev_outcome) * power[i] {{ weights.likelihood }}; + } + +} diff --git a/inst/stan/surv/exp_fpp.stan b/inst/stan/surv/exp_fpp.stan new file mode 100644 index 00000000..320ec6ef --- /dev/null +++ b/inst/stan/surv/exp_fpp.stan @@ -0,0 +1,52 @@ +// Exponential survival model +// Fixed power prior + +data { + + int N; // number of observations + vector[N] trt; // treatment indicator + vector[N] time; // survival time + vector[N] cens; // censoring indicator + vector[N] power; // power parameter + + {{ weights.data }} + {{ cov.data }} + +} + +parameters { + + real beta_trt; // treatment effect + real alpha; // baseline hazard + + {{ cov.parameters }} + +} + +transformed parameters { + + real HR_trt = exp(beta_trt); + +} + +model { + + vector[N] lp; + vector[N] elp; + + {{ trt.prior }} + {{ cov.priors }} + {{ baseline.prior }} + + lp = alpha + trt * beta_trt {{ cov.linpred }} ; + elp = exp(lp); + + for (i in 1:N) { + if (cens[i] == 1) { + target += exponential_lccdf(time[i] | elp[i]) * power[i] {{ weights.likelihood }}; + } else { + target += exponential_lpdf(time[i] | elp[i]) * power[i] {{ weights.likelihood }}; + } + } + +} diff --git a/inst/stan/surv/pem_fpp.stan b/inst/stan/surv/pem_fpp.stan new file mode 100644 index 00000000..9fd41129 --- /dev/null +++ b/inst/stan/surv/pem_fpp.stan @@ -0,0 +1,58 @@ +// Piecewise exponential survival model +// Fixed power prior + +data { + + int N; // number of observation-periods + vector[N] trt; // treatment indicator + vector[N] time; // survival time + vector[N] cens; // censoring indicator + int N_periods; // number of periods + matrix[N, N_periods] Z; // period indicators + vector[N] power; // power parameter + + {{ weights.data }} + {{ cov.data }} + +} + +parameters { + + real beta_trt; // treatment effect + vector[N_periods] alpha; // baseline hazard + + {{ cov.parameters }} + +} + +transformed parameters { + + real HR_trt = exp(beta_trt); + +} + +model { + + vector[N] lp; + vector[N] elp; + real sigma; + + {{ trt.prior }} + {{ cov.priors }} + + for (i in 1:N_periods) { + alpha[i] ~ {{ baseline.prior }}; + } + + lp = Z * alpha + trt * beta_trt {{ cov.linpred }} ; + elp = exp(lp); + + for (i in 1:N) { + if (cens[i] == 1) { + target += exponential_lccdf(time[i] | elp[i]) * power[i] {{ weights.likelihood }}; + } else { + target += exponential_lpdf(time[i] | elp[i]) * power[i] {{ weights.likelihood }}; + } + } + +} diff --git a/inst/stan/surv/weib_ph_fpp.stan b/inst/stan/surv/weib_ph_fpp.stan new file mode 100644 index 00000000..0e96cb40 --- /dev/null +++ b/inst/stan/surv/weib_ph_fpp.stan @@ -0,0 +1,73 @@ +// Weibull proportional hazards survival model +// Fixed power prior + +functions { + + real weibull_ph_lpdf(real y, real alpha, real lambda) { + real lprob = log(alpha) + log(lambda) + (alpha - 1) * log(y) - lambda * (y^alpha); + return lprob; + } + + real weibull_ph_lcdf(real y, real alpha, real lambda) { + real lprob = log(1 - exp(-lambda * y^alpha)); + return lprob; + } + + real weibull_ph_lccdf(real y, real alpha, real lambda) { + real lprob = -lambda * y^alpha; + return lprob; + } + +} + +data { + + int N; // number of observations + vector[N] trt; // treatment indicator + vector[N] time; // survival time + vector[N] cens; // censoring indicator + vector[N] power; // power parameter + + {{ weights.data }} + {{ cov.data }} + +} + +parameters { + + real beta_trt; // treatment effect + real alpha; // baseline hazard + real shape_weibull; // weibull shape parameter + + {{ cov.parameters }} + +} + +transformed parameters { + + real HR_trt = exp(beta_trt); + +} + +model { + + vector[N] lp; + vector[N] elp; + + {{ trt.prior }} + {{ cov.priors }} + {{ baseline.prior }} + {{ shape.prior }} + + lp = alpha + trt * beta_trt {{ cov.linpred }} ; + elp = exp(lp); + + for (i in 1:N) { + if (cens[i] == 1) { + target += weibull_ph_lccdf(time[i] | shape_weibull, elp[i]) * power[i] {{ weights.likelihood }}; + } else { + target += weibull_ph_lpdf(time[i] | shape_weibull, elp[i]) * power[i] {{ weights.likelihood }}; + } + } + +} From aa6aa1bc320e4e13ddea08312c49ca3e7adcf3ce Mon Sep 17 00:00:00 2001 From: Isaac Gravestock Date: Thu, 6 Feb 2025 11:53:15 +0100 Subject: [PATCH 2/6] add borrowing class --- DESCRIPTION | 1 + NAMESPACE | 1 + R/borrowing_fixed_power_prior.R | 56 +++++++++++++++++++ inst/WORDLIST | 14 +++++ man/Borrowing-class.Rd | 3 +- man/BorrowingFixedPowerPrior-class.Rd | 28 ++++++++++ man/BorrowingFull-class.Rd | 3 +- ...BorrowingHierarchicalCommensurate-class.Rd | 3 +- man/BorrowingNone-class.Rd | 3 +- man/borrowing_fixed_power_prior.Rd | 25 +++++++++ man/trim_cols.Rd | 5 +- 11 files changed, 137 insertions(+), 5 deletions(-) create mode 100644 R/borrowing_fixed_power_prior.R create mode 100644 man/BorrowingFixedPowerPrior-class.Rd create mode 100644 man/borrowing_fixed_power_prior.Rd diff --git a/DESCRIPTION b/DESCRIPTION index 9d8fb07a..70519b08 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -122,6 +122,7 @@ Collate: 'outcome_class.R' 'analysis_class.R' 'borrowing_details.R' + 'borrowing_fixed_power_prior.R' 'borrowing_full.R' 'borrowing_hierarchical_commensurate.R' 'borrowing_none.R' diff --git a/NAMESPACE b/NAMESPACE index 92e1526c..c4f98c60 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -9,6 +9,7 @@ export(beta_prior) export(bin_var) export(binary_cutoff) export(borrowing_details) +export(borrowing_fixed_power_prior) export(borrowing_full) export(borrowing_hierarchical_commensurate) export(borrowing_none) diff --git a/R/borrowing_fixed_power_prior.R b/R/borrowing_fixed_power_prior.R new file mode 100644 index 00000000..fb8fe8a7 --- /dev/null +++ b/R/borrowing_fixed_power_prior.R @@ -0,0 +1,56 @@ +#' `BorrowingFixedPowerPrior` class +#' +#' @slot method_name string. The name of the method. +#' @slot ext_flag_col character. Name of the external flag column in the matrix. +#' @slot powers numeric. The fixed values to be used as the powers in the power prior. +#' +#' @include borrowing_class.R +#' @family borrowing classes +.borrowing_fixed_power_prior <- setClass( + "BorrowingFixedPowerPrior", + slots = c( + power_col = "character" + ), + prototype = list( + method_name = "Borrowing with fixed power prior" + ), + contains = "Borrowing" +) + +#' Fixed Power Prior Borrowing +#' +#' @param ext_flag_col character. Name of the external flag column in the matrix. +#' @param power_col character. Name of the column containing values to be used as the power parameters. +#' +#' @return Object of class [`BorrowingFixedPowerPrior`][BorrowingFixedPowerPrior-class]. +#' @export +#' @examples +#' borrowing_fixed_power_prior( +#' ext_flag_col = "ext", +#' power_col = "power") +#' ) +borrowing_fixed_power_prior <- function(ext_flag_col, power_col) { + assert_string(ext_flag_col) + assert_string(power_col) + .borrowing_fixed_power_prior(ext_flag_col = ext_flag_col, power_col = power_col) +} + +# show ---- +setMethod( + f = "show", + signature = "BorrowingFixedPowerPrior", + definition = function(object) { + callNextMethod() + } +) + +# trim cols ---- +#' @rdname trim_cols +#' @include generics.R +setMethod( + f = "trim_cols", + signature = "BorrowingFixedPowerPrior", + definition = function(borrowing_object, analysis_object) { + c(ext_flag_col = object@ext_flag_col, power_col = power_col) + } +) diff --git a/inst/WORDLIST b/inst/WORDLIST index 3dde3d00..c05ff956 100644 --- a/inst/WORDLIST +++ b/inst/WORDLIST @@ -15,6 +15,7 @@ biopharmaceutical bootswatch Bove cauchy +cens ci CJ CMD @@ -30,8 +31,10 @@ costructor cov dataLayer desc +dev doi eg +elp Enas estimands et @@ -70,8 +73,16 @@ Kexin Khanal Kinnersley knitr +lccdf +lcdf Lifecycle Lindborg +linpred +lp +lpdf +lprob +lupdf +lupmf Mandrekar Manoj MatchIt @@ -87,6 +98,7 @@ operatorname ORCID parametrized patsubst +ph Pharmaceut poisson propto @@ -113,6 +125,7 @@ sourcenotes src stan Statist +stdev survial tabset tbody @@ -126,6 +139,7 @@ UA vgwzvpznzo Viele Walkthrough +weibull wpfbyktkdu www Yichen diff --git a/man/Borrowing-class.Rd b/man/Borrowing-class.Rd index 3b1ec53c..59537dc1 100644 --- a/man/Borrowing-class.Rd +++ b/man/Borrowing-class.Rd @@ -20,7 +20,8 @@ A class for defining borrowing details. Objects of class \seealso{ Prior constructor functions: \code{\link[=borrowing_full]{borrowing_full()}}, \code{\link[=borrowing_hierarchical_commensurate]{borrowing_hierarchical_commensurate()}}, \code{\link[=borrowing_none]{borrowing_none()}} -Other borrowing classes: +Other borrowing classes: +\code{\link{BorrowingFixedPowerPrior-class}}, \code{\link{BorrowingFull-class}}, \code{\link{BorrowingHierarchicalCommensurate-class}}, \code{\link{BorrowingNone-class}} diff --git a/man/BorrowingFixedPowerPrior-class.Rd b/man/BorrowingFixedPowerPrior-class.Rd new file mode 100644 index 00000000..dba3c605 --- /dev/null +++ b/man/BorrowingFixedPowerPrior-class.Rd @@ -0,0 +1,28 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/borrowing_fixed_power_prior.R +\docType{class} +\name{BorrowingFixedPowerPrior-class} +\alias{BorrowingFixedPowerPrior-class} +\alias{.borrowing_fixed_power_prior} +\title{\code{BorrowingFixedPowerPrior} class} +\description{ +\code{BorrowingFixedPowerPrior} class +} +\section{Slots}{ + +\describe{ +\item{\code{method_name}}{string. The name of the method.} + +\item{\code{ext_flag_col}}{character. Name of the external flag column in the matrix.} + +\item{\code{powers}}{numeric. The fixed values to be used as the powers in the power prior.} +}} + +\seealso{ +Other borrowing classes: +\code{\link{Borrowing-class}}, +\code{\link{BorrowingFull-class}}, +\code{\link{BorrowingHierarchicalCommensurate-class}}, +\code{\link{BorrowingNone-class}} +} +\concept{borrowing classes} diff --git a/man/BorrowingFull-class.Rd b/man/BorrowingFull-class.Rd index 309ab99c..3056d881 100644 --- a/man/BorrowingFull-class.Rd +++ b/man/BorrowingFull-class.Rd @@ -22,8 +22,9 @@ should not be created directly but by the constructor }} \seealso{ -Other borrowing classes: +Other borrowing classes: \code{\link{Borrowing-class}}, +\code{\link{BorrowingFixedPowerPrior-class}}, \code{\link{BorrowingHierarchicalCommensurate-class}}, \code{\link{BorrowingNone-class}} } diff --git a/man/BorrowingHierarchicalCommensurate-class.Rd b/man/BorrowingHierarchicalCommensurate-class.Rd index 5187dba1..352d0048 100644 --- a/man/BorrowingHierarchicalCommensurate-class.Rd +++ b/man/BorrowingHierarchicalCommensurate-class.Rd @@ -23,8 +23,9 @@ should not be created directly but by the constructor }} \seealso{ -Other borrowing classes: +Other borrowing classes: \code{\link{Borrowing-class}}, +\code{\link{BorrowingFixedPowerPrior-class}}, \code{\link{BorrowingFull-class}}, \code{\link{BorrowingNone-class}} } diff --git a/man/BorrowingNone-class.Rd b/man/BorrowingNone-class.Rd index 6ca73db4..477b3205 100644 --- a/man/BorrowingNone-class.Rd +++ b/man/BorrowingNone-class.Rd @@ -20,8 +20,9 @@ should not be created directly but by the constructor }} \seealso{ -Other borrowing classes: +Other borrowing classes: \code{\link{Borrowing-class}}, +\code{\link{BorrowingFixedPowerPrior-class}}, \code{\link{BorrowingFull-class}}, \code{\link{BorrowingHierarchicalCommensurate-class}} } diff --git a/man/borrowing_fixed_power_prior.Rd b/man/borrowing_fixed_power_prior.Rd new file mode 100644 index 00000000..4fb14444 --- /dev/null +++ b/man/borrowing_fixed_power_prior.Rd @@ -0,0 +1,25 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/borrowing_fixed_power_prior.R +\name{borrowing_fixed_power_prior} +\alias{borrowing_fixed_power_prior} +\title{Fixed Power Prior Borrowing} +\usage{ +borrowing_fixed_power_prior(ext_flag_col, power_col) +} +\arguments{ +\item{ext_flag_col}{character. Name of the external flag column in the matrix.} + +\item{power_col}{character. Name of the column containing values to be used as the power parameters.} +} +\value{ +Object of class \code{\link[=BorrowingFixedPowerPrior-class]{BorrowingFixedPowerPrior}}. +} +\description{ +Fixed Power Prior Borrowing +} +\examples{ +borrowing_fixed_power_prior( + ext_flag_col = "ext", + power_col = "power") +) +} diff --git a/man/trim_cols.Rd b/man/trim_cols.Rd index 5345369c..660bcf3b 100644 --- a/man/trim_cols.Rd +++ b/man/trim_cols.Rd @@ -1,9 +1,10 @@ % Generated by roxygen2: do not edit by hand % Please edit documentation in R/generics.R, R/borrowing_class.R, -% R/borrowing_hierarchical_commensurate.R +% R/borrowing_fixed_power_prior.R, R/borrowing_hierarchical_commensurate.R \name{trim_cols} \alias{trim_cols} \alias{trim_cols,Borrowing-method} +\alias{trim_cols,BorrowingFixedPowerPrior-method} \alias{trim_cols,BorrowingHierarchicalCommensurate-method} \title{Trim columns from Data Matrix Based on Borrowing object type} \usage{ @@ -11,6 +12,8 @@ trim_cols(borrowing_object, analysis_object) \S4method{trim_cols}{Borrowing}(borrowing_object, analysis_object) +\S4method{trim_cols}{BorrowingFixedPowerPrior}(borrowing_object, analysis_object) + \S4method{trim_cols}{BorrowingHierarchicalCommensurate}(borrowing_object, analysis_object) } \arguments{ From 0a9968b7290df00c420a877dae7933fce37a8cf4 Mon Sep 17 00:00:00 2001 From: Isaac Gravestock Date: Thu, 6 Feb 2025 12:03:41 +0100 Subject: [PATCH 3/6] data preparation --- R/prepare_stan_data_inputs.R | 96 ++++++++++++++++++++++++++++++++++-- 1 file changed, 93 insertions(+), 3 deletions(-) diff --git a/R/prepare_stan_data_inputs.R b/R/prepare_stan_data_inputs.R index 0eab972a..76243cb3 100644 --- a/R/prepare_stan_data_inputs.R +++ b/R/prepare_stan_data_inputs.R @@ -41,7 +41,6 @@ setGeneric( #' @return Named list of data inputs with covariates and weights added. #' @noRd add_covariates_and_weights <- function(data_in, analysis_obj, data_matrix) { - ## Covariate additions if (!is.null(analysis_obj@covariates)) { data_in[["K"]] <- NROW(analysis_obj@covariates@covariates) @@ -84,6 +83,27 @@ setMethod( } ) +### Fixed power prior borrowing ---- +setMethod( + f = "prepare_stan_data_inputs", + signature = c("OutcomeSurvExponentialWeibull", "BorrowingFixedPowerPrior", "ANY"), + definition = function(outcome, borrowing, analysis_obj) { + data_matrix <- trim_data_matrix(analysis_obj) + data_in <- list( + N = nrow(data_matrix), + trt = data_matrix[, analysis_obj@treatment@trt_flag_col], + time = data_matrix[, outcome@time_var], + cens = data_matrix[, outcome@cens_var], + power = data_matrix[, borrowing@power_col] + ) + + # Add covariates and weights + data_in <- add_covariates_and_weights(data_in, analysis_obj, data_matrix) + + return(data_in) + } +) + ### Hierarchical commensurate prior borrowing ---- setMethod( f = "prepare_stan_data_inputs", @@ -115,7 +135,6 @@ setMethod( f = "prepare_stan_data_inputs", signature = c("OutcomeSurvPEM", "BorrowingNoneFull", "ANY"), definition = function(outcome, borrowing, analysis_obj) { - analysis_obj@data_matrix <- trim_data_matrix(analysis_obj) analysis_obj <- cast_mat_to_long_pem(analysis_obj) data_matrix <- analysis_obj@data_matrix @@ -142,12 +161,43 @@ setMethod( } ) +### Fixed power prior borrowing ---- +setMethod( + f = "prepare_stan_data_inputs", + signature = c("OutcomeSurvPEM", "BorrowingFixedPowerPrior", "ANY"), + definition = function(outcome, borrowing, analysis_obj) { + analysis_obj@data_matrix <- trim_data_matrix(analysis_obj) + analysis_obj <- cast_mat_to_long_pem(analysis_obj) + data_matrix <- analysis_obj@data_matrix + + n_periods <- analysis_obj@outcome@n_periods + Z <- matrix(0, nrow = nrow(data_matrix), ncol = n_periods) + for (i in seq_len(nrow(data_matrix))) { + period <- data_matrix[i, "__period__"] + Z[i, period] <- 1 + } + data_in <- list( + N = nrow(data_matrix), + trt = data_matrix[, analysis_obj@treatment@trt_flag_col], + time = data_matrix[, outcome@time_var], + cens = data_matrix[, outcome@cens_var], + N_periods = max(data_matrix[, "__period__"]), + Z = Z, + power = data_matrix[, borrowing@power_col] + ) + + # Add covariates and weights + data_in <- add_covariates_and_weights(data_in, analysis_obj, data_matrix) + + return(data_in) + } +) + ### Hierarchical Commensurate Borrowing ---- setMethod( f = "prepare_stan_data_inputs", signature = c("OutcomeSurvPEM", "BorrowingHierarchicalCommensurate", "ANY"), definition = function(outcome, borrowing, analysis_obj) { - analysis_obj <- cast_mat_to_long_pem(analysis_obj) data_matrix <- analysis_obj@data_matrix @@ -197,6 +247,26 @@ setMethod( } ) +### Fixed power prior borrowing ---- +setMethod( + f = "prepare_stan_data_inputs", + signature = c("OutcomeBinaryLogistic", "BorrowingFixedPowerPrior", "ANY"), + definition = function(outcome, borrowing, analysis_obj) { + data_matrix <- trim_data_matrix(analysis_obj) + data_in <- list( + N = nrow(data_matrix), + trt = data_matrix[, analysis_obj@treatment@trt_flag_col], + y = data_matrix[, outcome@binary_var], + power = data_matrix[, borrowing@power_col] + ) + + # Add covariates and weights + data_in <- add_covariates_and_weights(data_in, analysis_obj, data_matrix) + + return(data_in) + } +) + ### Hierarchical commensurate prior borrowing ---- setMethod( f = "prepare_stan_data_inputs", @@ -241,6 +311,26 @@ setMethod( } ) +### Fixed power prior borrowing ---- +setMethod( + f = "prepare_stan_data_inputs", + signature = c("OutcomeContinuousNormal", "BorrowingFixedPowerPrior", "ANY"), + definition = function(outcome, borrowing, analysis_obj) { + data_matrix <- trim_data_matrix(analysis_obj) + data_in <- list( + N = nrow(data_matrix), + trt = data_matrix[, analysis_obj@treatment@trt_flag_col], + y = data_matrix[, outcome@continuous_var], + power = data_matrix[, borrowing@power_col] + ) + + # Add covariates and weights + data_in <- add_covariates_and_weights(data_in, analysis_obj, data_matrix) + + return(data_in) + } +) + ### Hierarchical commensurate prior borrowing ---- setMethod( f = "prepare_stan_data_inputs", From 482730595de624b72b7c001f5728a44c35362a85 Mon Sep 17 00:00:00 2001 From: Isaac Gravestock Date: Thu, 6 Feb 2025 12:10:45 +0100 Subject: [PATCH 4/6] load models --- R/load_and_interpolate_stan_model.R | 90 ++++++++++++++++++++++++++++- 1 file changed, 89 insertions(+), 1 deletion(-) diff --git a/R/load_and_interpolate_stan_model.R b/R/load_and_interpolate_stan_model.R index d07a4680..40020c8b 100644 --- a/R/load_and_interpolate_stan_model.R +++ b/R/load_and_interpolate_stan_model.R @@ -8,7 +8,7 @@ #' @return String containing the interpolated Stan model #' @include outcome_surv_pem.R build_model_string <- function(template_domain, template_filename, outcome, borrowing, analysis_obj, ...) { - + # Load the Stan template template <- load_stan_file(template_domain, template_filename) @@ -72,6 +72,23 @@ setMethod( } ) +### Fixed power prior borrowing ---- +setMethod( + f = "load_and_interpolate_stan_model", + signature = c("OutcomeSurvExponential", "BorrowingFixedPowerPrior", "ANY"), + definition = function(outcome, borrowing, analysis_obj) { + model_string <- build_model_string( + template_domain = "surv", + template_filename = "exp_fpp.stan", + outcome = outcome, + borrowing = borrowing, + analysis_obj = analysis_obj + ) + + return(model_string) + } +) + ### Hierarchical commensurate prior borrowing ---- setMethod( f = "load_and_interpolate_stan_model", @@ -110,6 +127,24 @@ setMethod( } ) +### Fixed power prior borrowing ---- +setMethod( + f = "load_and_interpolate_stan_model", + signature = c("OutcomeSurvWeibullPH", "BorrowingFixedPowerPrior", "ANY"), + definition = function(outcome, borrowing, analysis_obj) { + model_string <- build_model_string( + template_domain = "surv", + template_filename = "weib_ph_fpp.stan", + outcome = outcome, + borrowing = borrowing, + analysis_obj = analysis_obj, + shape.prior = h_glue("shape_weibull ~ {{get_prior_string(outcome@param_priors$shape_weibull)}} ;") + ) + + return(model_string) + } +) + ### Hierarchical commensurate prior borrowing ---- setMethod( f = "load_and_interpolate_stan_model", @@ -149,6 +184,24 @@ setMethod( } ) +### Fixed power prior borrowing ---- +setMethod( + f = "load_and_interpolate_stan_model", + signature = c("OutcomeSurvPEM", "BorrowingFixedPowerPrior", "ANY"), + definition = function(outcome, borrowing, analysis_obj) { + model_string <- build_model_string( + template_domain = "surv", + template_filename = "pem_fpp.stan", + outcome = outcome, + borrowing = borrowing, + analysis_obj = analysis_obj, + baseline.prior = get_prior_string(outcome@baseline_prior) + ) + + return(model_string) + } +) + ### Hierarchical commensurate prior borrowing ---- setMethod( f = "load_and_interpolate_stan_model", @@ -188,6 +241,23 @@ setMethod( } ) +### Fixed power prior borrowing ---- +setMethod( + f = "load_and_interpolate_stan_model", + signature = c("OutcomeBinaryLogistic", "BorrowingFixedPowerPrior", "ANY"), + definition = function(outcome, borrowing, analysis_obj) { + model_string <- build_model_string( + template_domain = "bin", + template_filename = "logit_fpp.stan", + outcome = outcome, + borrowing = borrowing, + analysis_obj = analysis_obj + ) + + return(model_string) + } +) + ### Hierarchical commensurate prior borrowing ---- setMethod( f = "load_and_interpolate_stan_model", @@ -228,6 +298,24 @@ setMethod( } ) +### Fixed power prior borrowing ---- +setMethod( + f = "load_and_interpolate_stan_model", + signature = c("OutcomeContinuousNormal", "BorrowingFixedPowerPrior", "ANY"), + definition = function(outcome, borrowing, analysis_obj) { + model_string <- build_model_string( + template_domain = "cont", + template_filename = "gauss_fpp.stan", + outcome = outcome, + borrowing = borrowing, + analysis_obj = analysis_obj, + stdev.prior = h_glue("std_dev_outcome ~ {{get_prior_string(outcome@param_priors$std_dev_outcome)}} ;") + ) + + return(model_string) + } +) + ### Hierarchical commensurate prior borrowing ---- setMethod( f = "load_and_interpolate_stan_model", From 3cf822895a1cb2f2546f5e71131841e9dec5ddeb Mon Sep 17 00:00:00 2001 From: Isaac Gravestock Date: Fri, 7 Feb 2025 15:02:20 +0100 Subject: [PATCH 5/6] add tests --- .Rbuildignore | 3 +- R/borrowing_fixed_power_prior.R | 14 +- R/check_data_matrix_has_columns.R | 4 +- man/borrowing_fixed_power_prior.Rd | 2 +- man/get_vars.Rd | 9 +- .../test-borrowing_fixed_power_prior.R | 14 ++ tests/testthat/test-create_analysis_obj.R | 14 ++ tests/testthat/test-mcmc_sample-analysis.R | 212 ++++++++++++++++++ tests/testthat/test-trim_data_matrix.R | 20 ++ 9 files changed, 283 insertions(+), 9 deletions(-) create mode 100644 tests/testthat/test-borrowing_fixed_power_prior.R diff --git a/.Rbuildignore b/.Rbuildignore index fafba974..7f222f31 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -29,4 +29,5 @@ develop ^dev$ ^\.pre-commit-config\.yaml$ ^vignettes/original$ -^vignettes/Makefile$ \ No newline at end of file +^vignettes/Makefile$ +^node_modules$ diff --git a/R/borrowing_fixed_power_prior.R b/R/borrowing_fixed_power_prior.R index fb8fe8a7..21fe80eb 100644 --- a/R/borrowing_fixed_power_prior.R +++ b/R/borrowing_fixed_power_prior.R @@ -27,7 +27,7 @@ #' @examples #' borrowing_fixed_power_prior( #' ext_flag_col = "ext", -#' power_col = "power") +#' power_col = "power" #' ) borrowing_fixed_power_prior <- function(ext_flag_col, power_col) { assert_string(ext_flag_col) @@ -51,6 +51,16 @@ setMethod( f = "trim_cols", signature = "BorrowingFixedPowerPrior", definition = function(borrowing_object, analysis_object) { - c(ext_flag_col = object@ext_flag_col, power_col = power_col) + get_vars(analysis_object) + } +) +# get vars ----- +#' @rdname get_vars +#' @include generics.R +setMethod( + f = "get_vars", + signature = "BorrowingFixedPowerPrior", + definition = function(object) { + c(ext_flag_col = object@ext_flag_col, power_col = object@power_col) } ) diff --git a/R/check_data_matrix_has_columns.R b/R/check_data_matrix_has_columns.R index 37ce78a1..a33b6e34 100644 --- a/R/check_data_matrix_has_columns.R +++ b/R/check_data_matrix_has_columns.R @@ -39,8 +39,8 @@ check_data_matrix_has_columns <- function(object) { error_col <- c(error_col, get_vars(object@treatment)) } - if (!get_vars(object@borrowing) %in% data_cols) { - error_col <- c(error_col, get_vars(object@borrowing)) + if (any(missing_outcomes <- !get_vars(object@borrowing) %in% data_cols)) { + error_col <- c(error_col, get_vars(object@borrowing)[missing_outcomes]) } missing_covs <- setdiff(get_vars(object@covariates), data_cols) diff --git a/man/borrowing_fixed_power_prior.Rd b/man/borrowing_fixed_power_prior.Rd index 4fb14444..04124efa 100644 --- a/man/borrowing_fixed_power_prior.Rd +++ b/man/borrowing_fixed_power_prior.Rd @@ -20,6 +20,6 @@ Fixed Power Prior Borrowing \examples{ borrowing_fixed_power_prior( ext_flag_col = "ext", - power_col = "power") + power_col = "power" ) } diff --git a/man/get_vars.Rd b/man/get_vars.Rd index aeaf4373..bc1aa9fc 100644 --- a/man/get_vars.Rd +++ b/man/get_vars.Rd @@ -1,9 +1,9 @@ % Generated by roxygen2: do not edit by hand % Please edit documentation in R/generics.R, R/covariate_class.R, % R/treatment_class.R, R/borrowing_class.R, R/outcome_class.R, -% R/analysis_class.R, R/sim_treatment_list.R, R/sim_outcome_list.R, -% R/sim_borrowing_list.R, R/sim_covariate_list.R, R/simulation_class.R, -% R/simulate_data_baseline.R +% R/analysis_class.R, R/borrowing_fixed_power_prior.R, R/sim_treatment_list.R, +% R/sim_outcome_list.R, R/sim_borrowing_list.R, R/sim_covariate_list.R, +% R/simulation_class.R, R/simulate_data_baseline.R \name{get_vars} \alias{get_vars} \alias{get_vars,Covariates-method} @@ -14,6 +14,7 @@ \alias{get_vars,ContinuousOutcome-method} \alias{get_vars,Analysis-method} \alias{get_vars,NULL-method} +\alias{get_vars,BorrowingFixedPowerPrior-method} \alias{get_vars,SimTreatmentList-method} \alias{get_vars,SimOutcomeList-method} \alias{get_vars,SimBorrowingList-method} @@ -40,6 +41,8 @@ get_vars(object) \S4method{get_vars}{NULL}(object) +\S4method{get_vars}{BorrowingFixedPowerPrior}(object) + \S4method{get_vars}{SimTreatmentList}(object) \S4method{get_vars}{SimOutcomeList}(object) diff --git a/tests/testthat/test-borrowing_fixed_power_prior.R b/tests/testthat/test-borrowing_fixed_power_prior.R new file mode 100644 index 00000000..6ad83d4e --- /dev/null +++ b/tests/testthat/test-borrowing_fixed_power_prior.R @@ -0,0 +1,14 @@ +test_that("borrowing_fixed_power_prior works as expected", { + fpp <- borrowing_fixed_power_prior("ext", "pp") + expect_class(fpp, "Borrowing") + expect_class(fpp, "BorrowingFixedPowerPrior") + expect_equal(fpp@ext_flag_col, "ext") + expect_equal(fpp@power_col, "pp") +}) + +test_that("get_vars works for BorrowingFixedPowerPrior", { + expect_identical( + get_vars(borrowing_fixed_power_prior("ext_fl", "power")), + c(ext_flag_col = "ext_fl", power_col = "power") + ) +}) diff --git a/tests/testthat/test-create_analysis_obj.R b/tests/testthat/test-create_analysis_obj.R index 19a003c8..9704caac 100644 --- a/tests/testthat/test-create_analysis_obj.R +++ b/tests/testthat/test-create_analysis_obj.R @@ -196,6 +196,20 @@ test_that("Columns in analysis_obj should be in matrix", { ), "The following specified variables were not found in `data_matrix`:\n ext_flag_col: tira" ) + + expect_error( + create_analysis_obj( + data_matrix = example_matrix, + covariates = ac, + outcome = esd, + treatment = td, + borrowing = borrowing_fixed_power_prior( + ext_flag_col = "ext", + power_col = "power_par" + ) + ), + "The following specified variables were not found in `data_matrix`:\n power_col: power_par" + ) }) diff --git a/tests/testthat/test-mcmc_sample-analysis.R b/tests/testthat/test-mcmc_sample-analysis.R index ee33a29c..403bb0c6 100644 --- a/tests/testthat/test-mcmc_sample-analysis.R +++ b/tests/testthat/test-mcmc_sample-analysis.R @@ -713,3 +713,215 @@ test_that("mcmc_sample for Analysis works for BDB, piecewise exponential dist", expect_true(tau_incommens_aggr > tau_incommens_conserv) }) + +# Fixed Power Prior ------ +test_that("mcmc_sample for Analysis works for fixed power prior borrowing, exponential dist", { + skip_on_cran() + skip_on_ci() + power <- 1 - 0.3 * example_matrix[, "ext"] + + exp_c1 <- exp(coef(flexsurv::flexsurvreg( + survival::Surv(time, 1 - cnsr) ~ trt + cov1, + data = as.data.frame(example_matrix), + dist = "exponential", + weights = power + ))[["trt"]]) + + + fpp_exp_bayes_ao <- create_analysis_obj( + data_matrix = cbind(example_matrix, power = power), + outcome = outcome_surv_exponential("time", "cnsr", prior_normal(0, 100000)), + borrowing = borrowing_fixed_power_prior( + ext_flag_col = "ext", + power_col = "power" + ), + treatment = treatment_details("trt", prior_normal(0, 100000)) + ) + + fpp_exp_bayes <- mcmc_sample( + fpp_exp_bayes_ao, + iter_warmup = 2000, + iter_sampling = 2000, + chains = 1 + ) + expect_r6(fpp_exp_bayes, "CmdStanMCMC") + expect_equal( + fpp_exp_bayes$summary("HR_trt", "median")[[2]], + exp_c1, + tolerance = .05 + ) +}) + +test_that("mcmc_sample for Analysis works for fixed power prior, exponential dist, one covariate", { + skip_on_cran() + skip_on_ci() + power <- 1 - 0.3 * example_matrix[, "ext"] + fpp_exp_c1 <- exp(coef(flexsurv::flexsurvreg( + survival::Surv(time, 1 - cnsr) ~ trt + cov1, + data = as.data.frame(example_matrix), + dist = "exponential", + weights = power + ))[["trt"]]) + + fpp_exp_bayes_c1_ao <- create_analysis_obj( + data_matrix = cbind(example_matrix, power = power), + covariates = add_covariates("cov1", prior_normal(0, 100000)), + outcome = outcome_surv_exponential("time", "cnsr", prior_normal(0, 100000)), + borrowing = borrowing_fixed_power_prior( + ext_flag_col = "ext", + power_col = "power" + ), + treatment = treatment_details("trt", prior_normal(0, 100000)) + ) + + fpp_exp_bayes_c1 <- mcmc_sample( + fpp_exp_bayes_c1_ao, + iter_warmup = 2000, + iter_sampling = 2000, + chains = 1 + ) + expect_r6(fpp_exp_bayes_c1, "CmdStanMCMC") + expect_equal( + fpp_exp_bayes_c1$summary("HR_trt", "median")[[2]], + fpp_exp_c1, + tolerance = .05 + ) +}) + +test_that("mcmc_sample for Analysis works for fixed power prior borrowing, normal dist", { + set.seed(123) + power <- 1 - 0.5 * example_matrix[, "ext"] + outcome_col <- 5 + example_matrix[, "trt"] + example_matrix[, "cov1"] + 2 * example_matrix[, "cov2"] + + 0.5 * example_matrix[, "cov3"] - 1.5 * example_matrix[, "cov4"] + rnorm(500, 0, 1) + + outcome <- outcome_cont_normal( + continuous_var = "outcome", + baseline_prior = prior_normal(0, 100), + std_dev_prior = prior_half_cauchy(1, 5) + ) + borrowing <- borrowing_fixed_power_prior( + ext_flag_col = "ext", + power_col = "power" + ) + treatment <- treatment_details( + trt_flag_col = "trt", + trt_prior = prior_normal(0, 1000) + ) + anls_obj <- create_analysis_obj( + data_matrix = cbind(example_matrix, outcome = outcome_col, power = power), + outcome = outcome, + borrowing = borrowing, + treatment = treatment, + quiet = FALSE + ) + result <- mcmc_sample( + anls_obj, + iter_warmup = 2000, + iter_sampling = 2000, + chains = 1 + ) + + result_summary <- result$summary(c("alpha", "beta_trt")) + expect_equal(result_summary[["median"]], c(6.44, 0.72), tolerance = .05) + expect_equal(result_summary[["q5"]], c(6.25, 0.38), tolerance = .05) + expect_equal(result_summary[["q95"]], c(6.62, 1.06), tolerance = .05) +}) + + +test_that("mcmc_sample for Analysis works for fixed power prior borrowing, logistic dist", { + skip_on_cran() + skip_on_ci() + set.seed(123) + power <- 1 - 0.5 * example_matrix[, "ext"] + fpp_bin_bayes_ao <- create_analysis_obj( + data_matrix = cbind(example_matrix, power = power), + outcome = outcome_bin_logistic("resp", prior_normal(0, 100000)), + borrowing = borrowing_fixed_power_prior("ext", power_col = "power"), + treatment = treatment_details("trt", prior_normal(0, 100000)) + ) + + fpp_bin_bayes <- mcmc_sample( + fpp_bin_bayes_ao, + iter_warmup = 2000, + iter_sampling = 2000, + chains = 1 + ) + expect_r6(fpp_bin_bayes, "CmdStanMCMC") + expect_equal( + fpp_bin_bayes$summary("OR_trt", "median")[[2]], + 1.59, + tolerance = .05 + ) +}) + +test_that("mcmc_sample for Analysis works for fixed power prior borrowing, PEM dist", { + skip_on_cran() + skip_on_ci() + cuts <- c(1, 5, 10) + set.seed(1234) + internal_as_external <- example_matrix[example_matrix[, "ext"] == 0 & example_matrix[, "trt"] == 0, ] + internal_as_external[, "ext"] <- 1 + internal_as_external[, "id"] <- seq(10000, 10000 + nrow(internal_as_external) - 1) + data_matrix <- rbind( + example_matrix[example_matrix[, "ext"] == 0, ], + internal_as_external + ) + power <- 1 - 0.5 * data_matrix[, "ext"] + data_matrix <- cbind(data_matrix, power = power) + + ## Conservative commensurate + fpp_pem_ao <- create_analysis_obj( + data_matrix = data_matrix, + outcome = outcome_surv_pem("time", "cnsr", prior_normal(0, 100000), cut_points = cuts), + borrowing = borrowing_fixed_power_prior("ext", "power"), + treatment = treatment_details("trt", prior_normal(0, 100000)) + ) + + fpp_pem <- mcmc_sample( + fpp_pem_ao, + iter_warmup = 2000, + iter_sampling = 5000, + chains = 2 + ) + beta_trt <- fpp_pem$summary("beta_trt")[["median"]] + expect_equal(beta_trt, -0.11, tolerance = 0.05) +}) + +test_that("mcmc_sample for Analysis works for fixed power prior borrowing, Weibull dist", { + skip_on_cran() + skip_on_ci() + + data_matrix <- cbind(example_matrix, power = 1 - 0.5 * example_matrix[, "ext"]) + + fpp_weib <- exp(coef(flexsurv::flexsurvreg( + survival::Surv(time, 1 - cnsr) ~ trt, + data = as.data.frame(data_matrix), + dist = "weibullPH", + weights = power + ))[["trt"]]) + + fpp_weib_bayes_ao <- create_analysis_obj( + data_matrix = data_matrix, + outcome = outcome_surv_weibull_ph( + "time", + "cnsr", + prior_normal(0, 100000), + prior_normal(0, 100000) + ), + borrowing = borrowing_fixed_power_prior("ext", "power"), + treatment = treatment_details("trt", prior_normal(0, 100000)) + ) + + fpp_weib_bayes <- mcmc_sample( + fpp_weib_bayes_ao, + iter_warmup = 2000, + iter_sampling = 2000, + chains = 1 + ) + expect_r6(fpp_weib_bayes, "CmdStanMCMC") + expect_equal( + fpp_weib_bayes$summary("HR_trt", "median")[[2]], + fpp_weib, + tolerance = .05 + ) +}) diff --git a/tests/testthat/test-trim_data_matrix.R b/tests/testthat/test-trim_data_matrix.R index 7770dadc..9e1e9b88 100644 --- a/tests/testthat/test-trim_data_matrix.R +++ b/tests/testthat/test-trim_data_matrix.R @@ -49,3 +49,23 @@ test_that("data matrix trimming works with No Borrowing", { colnames(result), c("time", "cnsr", "trt", "cov1") ) }) + + +test_that("data matrix trimming works with Fixed Power Prior Borrowing", { + object <- psborrow2:::.analysis_obj( + data_matrix = cbind(example_matrix, power = runif(500)), + outcome = outcome_surv_exponential("time", "cnsr", prior_normal(0, 1000)), + treatment = treatment_details("trt", trt_prior = prior_normal(0, 1000)), + covariates = add_covariates("cov1", prior_normal(0, 1000)), + borrowing = borrowing_fixed_power_prior( + ext_flag_col = "ext", + power_col = "power" + ) + ) + + result <- psborrow2:::trim_data_matrix(object) + expect_matrix(result, mode = "numeric", nrows = 500, ncols = 6) + expect_set_equal( + colnames(result), c("time", "cnsr", "ext", "trt", "power", "cov1") + ) +}) From 14698dc5562c72a95de5a888c604dfab845963e4 Mon Sep 17 00:00:00 2001 From: Isaac Gravestock Date: Wed, 26 Feb 2025 15:39:37 +0100 Subject: [PATCH 6/6] check powers are on [0,1] --- R/prepare_stan_data_inputs.R | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/R/prepare_stan_data_inputs.R b/R/prepare_stan_data_inputs.R index 76243cb3..b3ef55a1 100644 --- a/R/prepare_stan_data_inputs.R +++ b/R/prepare_stan_data_inputs.R @@ -89,6 +89,9 @@ setMethod( signature = c("OutcomeSurvExponentialWeibull", "BorrowingFixedPowerPrior", "ANY"), definition = function(outcome, borrowing, analysis_obj) { data_matrix <- trim_data_matrix(analysis_obj) + + assert_numeric(data_matrix[, borrowing@power_col], lower = 0, upper = 1, any.missing = FALSE) + data_in <- list( N = nrow(data_matrix), trt = data_matrix[, analysis_obj@treatment@trt_flag_col], @@ -170,6 +173,8 @@ setMethod( analysis_obj <- cast_mat_to_long_pem(analysis_obj) data_matrix <- analysis_obj@data_matrix + assert_numeric(data_matrix[, borrowing@power_col], lower = 0, upper = 1, any.missing = FALSE) + n_periods <- analysis_obj@outcome@n_periods Z <- matrix(0, nrow = nrow(data_matrix), ncol = n_periods) for (i in seq_len(nrow(data_matrix))) { @@ -253,6 +258,9 @@ setMethod( signature = c("OutcomeBinaryLogistic", "BorrowingFixedPowerPrior", "ANY"), definition = function(outcome, borrowing, analysis_obj) { data_matrix <- trim_data_matrix(analysis_obj) + + assert_numeric(data_matrix[, borrowing@power_col], lower = 0, upper = 1, any.missing = FALSE) + data_in <- list( N = nrow(data_matrix), trt = data_matrix[, analysis_obj@treatment@trt_flag_col], @@ -317,6 +325,9 @@ setMethod( signature = c("OutcomeContinuousNormal", "BorrowingFixedPowerPrior", "ANY"), definition = function(outcome, borrowing, analysis_obj) { data_matrix <- trim_data_matrix(analysis_obj) + + assert_numeric(data_matrix[, borrowing@power_col], lower = 0, upper = 1, any.missing = FALSE) + data_in <- list( N = nrow(data_matrix), trt = data_matrix[, analysis_obj@treatment@trt_flag_col],