start exponential

.pre-commit-config.yaml

@@ -3,7 +3,7 @@
 # R specific hooks: https://github.com/lorenzwalthert/precommit
 repos:
   - repo: https://github.com/lorenzwalthert/precommit
-    rev: v0.4.0
+    rev: v0.4.2
     hooks:
       - id: style-files
         args: []
@@ -65,7 +65,7 @@ repos:
     hooks:
       - id: prettier
   - repo: https://github.com/pre-commit/pre-commit-hooks
-    rev: v4.5.0
+    rev: v4.6.0
     hooks:
       - id: check-added-large-files
         args: ["--maxkb=200"]
@@ -101,7 +101,7 @@ repos:
         files: '\.Rhistory|\.RData|\.Rds|\.rds$'
         # `exclude: <regex>` to allow committing specific files.
   - repo: https://github.com/igorshubovych/markdownlint-cli
-    rev: v0.39.0
+    rev: v0.40.0
     hooks:
       - id: markdownlint
         args:

DESCRIPTION

@@ -54,14 +54,14 @@ Description: Bayesian dynamic borrowing is an approach to incorporating external
    of outcomes); see Viele 2013 <doi:10.1002/pst.1589> for an overview.
    This package implements the hierarchical commensurate prior approach to dynamic borrowing
    as described in Hobbes 2011 <doi:10.1111/j.1541-0420.2011.01564.x>.
-   There are three main functionalities. First, 'psborrow2' provides a user-friendly 
+   There are three main functionalities. First, 'psborrow2' provides a user-friendly
    interface for applying dynamic borrowing on the study results handles the Markov Chain
-   Monte Carlo sampling on behalf of the user. Second, 'psborrow2' provides a 
+   Monte Carlo sampling on behalf of the user. Second, 'psborrow2' provides a
    simulation framework to compare different borrowing parameters (e.g. full borrowing, no
    borrowing, dynamic borrowing) and other trial and borrowing characteristics
    (e.g. sample size, covariates) in a unified way. Third, 'psborrow2' provides
    a set of functions to generate data for simulation studies, and also allows
-   the user to specify their own data generation process. This package is designed to 
+   the user to specify their own data generation process. This package is designed to
    use the sampling functions from 'cmdstanr' which can be installed from
    <https://mc-stan.org/r-packages/>.
 URL: https://github.com/Genentech/psborrow2, https://genentech.github.io/psborrow2/index.html
@@ -145,8 +145,10 @@ Collate:
     'make_analysis_object_list.R'
     'make_model_string_data.R'
     'make_model_string_functions.R'
+    'outcome_surv_piecewise_exponential.R'
     'make_model_string_model.R'
     'make_model_string_parameters.R'
+    'make_model_string_transf_data.R'
     'make_model_string_transf_params.R'
     'mcmc_sample.R'
     'mcmc_simulation_result.R'

NAMESPACE

@@ -50,6 +50,7 @@ export(null_event_dist)
 export(outcome_bin_logistic)
 export(outcome_cont_normal)
 export(outcome_surv_exponential)
+export(outcome_surv_piecewise_exp)
 export(outcome_surv_weibull_ph)
 export(plot)
 export(plot_pdf)

R/create_analysis_obj.R

@@ -95,6 +95,7 @@ create_analysis_obj <- function(data_matrix,
   # Model string components
   functions_str <- make_model_string_functions(analysis_obj)
   data_str <- make_model_string_data(analysis_obj)
+  trans_data_str <- make_model_string_transf_data(analysis_obj)
   param_str <- make_model_string_parameters(analysis_obj)
   trans_param_str <- make_model_string_transf_param(analysis_obj)
   model_str <- make_model_string_model(analysis_obj@borrowing, analysis_obj@outcome, analysis_obj)
@@ -106,6 +107,8 @@ create_analysis_obj <- function(data_matrix,
 
     {{data_str}}
 
+    {{trans_data_str}}
+
     {{param_str}}
 
     {{trans_param_str}}

R/make_model_string_model.R

@@ -1,3 +1,6 @@
+#' @include outcome_surv_piecewise_exponential.R
+NULL
+
 #' @rdname make_model_string_model
 setMethod(
   "make_model_string_model",
@@ -62,12 +65,17 @@ make_model_string_full_none <- function(borrowing, outcome, analysis_obj) {
   }")
 }
 
+# BorrowingFull, ANY -------
 #' @rdname make_model_string_model
 setMethod("make_model_string_model", signature("BorrowingFull", "ANY", "Analysis"), make_model_string_full_none)
 
+# BorrowingNone, ANY -------
+
 #' @rdname make_model_string_model
 setMethod("make_model_string_model", signature("BorrowingNone", "ANY", "Analysis"), make_model_string_full_none)
 
+# BorrowingHierarchicalCommensurate, ANY -------
+
 #' @rdname make_model_string_model
 setMethod(
   "make_model_string_model",
@@ -134,3 +142,125 @@ setMethod(
   }")
   }
 )
+
+# BorrowingHierarchicalCommensurate, OutcomeSurvPiecewiseExponential -----
+#' @rdname make_model_string_model
+setMethod(
+  "make_model_string_model",
+  signature("BorrowingHierarchicalCommensurate", "OutcomeSurvPiecewiseExponential", "Analysis"),
+  function(borrowing, outcome, analysis_obj) {
+    ### Treatment prior
+    beta_trt_prior <- get_prior_string(analysis_obj@treatment@trt_prior)
+
+    ### Linear predictor
+    has_covariates <- !is.null(analysis_obj@covariates)
+
+    linear_predictor <- if (has_covariates) {
+      h_glue("lp = Z*alpha_ext' + (1-Z)*alpha' + rep_matrix(trt * beta_trt + X * beta, M);")
+    } else if (!has_covariates) {
+      h_glue("lp = Z*alpha_ext' + (1-Z)*alpha' + rep_matrix(trt * beta_trt, M);")
+    }
+
+    ### Add priors for relevant parameters
+    if (NROW(analysis_obj@outcome@param_priors) > 0) {
+      names <- names(analysis_obj@outcome@param_priors)
+      values <- get_prior_string(analysis_obj@outcome@param_priors)
+      outcome_prior <- h_glue("{{names}} ~ {{values}} ;", collapse = TRUE)
+    } else {
+      outcome_prior <- ""
+    }
+
+    ### Add priors on betas
+    if (has_covariates) {
+      i <- seq_along(analysis_obj@covariates@covariates)
+      value <- get_prior_string(analysis_obj@covariates@priors)
+      index <- if (test_named(value)) get_vars(analysis_obj@covariates) else rep(1, length(i))
+      covariate_prior <- h_glue("beta[{{i}}] ~ {{value[index]}} ;", collapse = TRUE)
+    } else {
+      covariate_prior <- ""
+    }
+
+    tau_prior <- get_prior_string(analysis_obj@borrowing@tau_prior)
+    alpha_2_prior <- get_prior_string(analysis_obj@outcome@baseline_prior)
+
+    borrowing_string <- h_glue("
+      tau ~ {{tau_prior}} ;
+      real sigma;
+      sigma = 1 / tau;
+      alpha_ext ~ {{alpha_2_prior}};
+      alpha ~ normal(alpha_ext, sqrt(sigma)) ;
+      ")
+
+    ### Add in likelihood function
+    likelihood_string <- "target += sum((lp .* D) - (exp(lp) .* T));"
+
+    h_glue("
+  model {
+    matrix[N,M] lp;
+    vector[N] elp;
+    beta_trt ~ {{beta_trt_prior}};
+    {{linear_predictor}}
+    {{outcome_prior}}
+    {{covariate_prior}}
+    {{borrowing_string}}
+    {{likelihood_string}}
+  }")
+  }
+)
+
+# BorrowingNone, OutcomeSurvPiecewiseExponential -----
+#' @rdname make_model_string_model
+setMethod(
+  "make_model_string_model",
+  signature("BorrowingNone", "OutcomeSurvPiecewiseExponential", "Analysis"),
+  function(borrowing, outcome, analysis_obj) {
+    ### Treatment prior
+    beta_trt_prior <- get_prior_string(analysis_obj@treatment@trt_prior)
+
+    ### Linear predictor
+    has_covariates <- !is.null(analysis_obj@covariates)
+
+    linear_predictor <- if (has_covariates) {
+      h_glue("lp = alpha' + rep_matrix(trt * beta_trt + X * beta, M);")
+    } else if (!has_covariates) {
+      h_glue("lp = alpha' + rep_matrix(trt * beta_trt, M);")
+    }
+
+    ### Add priors for relevant parameters
+    if (NROW(analysis_obj@outcome@param_priors) > 0) {
+      names <- names(analysis_obj@outcome@param_priors)
+      values <- get_prior_string(analysis_obj@outcome@param_priors)
+      outcome_prior <- h_glue("{{names}} ~ {{values}} ;", collapse = TRUE)
+    } else {
+      outcome_prior <- ""
+    }
+
+    ### Add priors on betas
+    if (has_covariates) {
+      i <- seq_along(analysis_obj@covariates@covariates)
+      value <- get_prior_string(analysis_obj@covariates@priors)
+      index <- if (test_named(value)) get_vars(analysis_obj@covariates) else rep(1, length(i))
+      covariate_prior <- h_glue("beta[{{i}}] ~ {{value[index]}} ;", collapse = TRUE)
+    } else {
+      covariate_prior <- ""
+    }
+
+    tau_prior <- get_prior_string(analysis_obj@borrowing@tau_prior)
+    alpha_prior <- get_prior_string(analysis_obj@outcome@baseline_prior)
+    borrowing_string <- h_glue("alpha ~ {{alpha_prior}};")
+
+    ### Add in likelihood function
+    likelihood_string <- "target += sum((lp .* D) - (exp(lp) .* T));"
+
+    h_glue("
+  model {
+    matrix[N,M] lp;
+    beta_trt ~ {{beta_trt_prior}};
+    {{linear_predictor}}
+    {{outcome_prior}}
+    {{covariate_prior}}
+    {{borrowing_string}}
+    {{likelihood_string}}
+  }")
+  }
+)

R/make_model_string_transf_data.R

@@ -0,0 +1,12 @@
+make_model_string_transf_data <- function(object) {
+  if (.hasSlot(object@outcome, "transformed_data_stan_code")) {
+    str <- object@outcome@transformed_data_stan_code
+    h_glue(
+      "transformed data {
+  {{str}}
+}"
+    )
+  } else {
+    ""
+  }
+}

R/outcome_surv_piecewise_exponential.R

@@ -0,0 +1,104 @@
+#' `OutcomeSurvPiecewiseExponential` Class
+#'
+#' A class for defining a time-to-event survival analysis with a piecewise
+#' exponential survival distribution.
+#' Objects of class `OutcomeSurvPiecewiseExponential` should not be created directly
+#' but by the constructor [outcome_surv_piecewise_exp()].
+#'
+#' @slot function_stan_code character. stan function code block containing text to interpolate into stan model.
+#' Empty string for `OutcomeSurvExponential`.
+#' @slot param_stan_code character. stan parameter code block containing text to interpolate into stan model.
+#' Empty string for `OutcomeSurvExponential`.
+#' @slot likelihood_stan_code character. stan model likelihood code block containing text
+#' to interpolate into stan model.
+#' @slot transformed_data_stan_code character. stan transformed data code block for constructing piecewise data
+#' @slot n_param integer. Number of ancillary parameters for the model to estimate (0).
+#' @slot param_priors list. Named list of prior distributions on the ancillary parameters in the model.
+#' Empty for `OutcomeSurvPiecewiseExponential`.
+#' @slot time_var character. Variable used for time in `TimeToEvent` objects.
+#' @slot cens_var character. Variable used for censoring in `TimeToEvent` objects.
+#' @slot baseline_prior `Prior`. Object of class `Prior`
+#' specifying prior distribution for the baseline outcome.
+#' @slot name_beta_trt. Named vector for beta_trt.
+#' @slot name_exp_trt. Named vector for exponentiated beta_trt
+#' @slot alpha_type. How to interpret alpha.
+#' @slot name_addnl_params. Named vector for additional parameters.
+#' @include outcome_class.R
+#' @family outcome
+.outcome_surv_piecewise_exp <- setClass(
+  "OutcomeSurvPiecewiseExponential",
+  contains = "TimeToEvent",
+  slots = c(
+    start_times = "numeric",
+    transformed_data_stan_code = "character"
+  ),
+  prototype = list(
+    n_param = 0L,
+    likelihood_stan_code = ""
+  )
+)
+
+#' Piecewise Exponential survival distribution
+#'
+#' @param time_var character. Name of time variable column in model matrix
+#' @param cens_var character. Name of the censorship variable flag in model matrix
+#' @param weight_var character. Optional name of variable in model matrix for weighting the log likelihood.
+#' @param baseline_prior `Prior`. Object of class `Prior` specifying prior distribution for the baseline outcome.
+#' @param start_times numeric. A vector of times to define periods.
+#'
+#' @return Object of class [`OutcomeSurvPiecewiseExponential`][OutcomeSurvPiecewiseExponential-class].
+#' @export
+#' @family outcome models
+#'
+#' @examples
+#' es <- outcome_surv_piecewise_exp(
+#'   time_var = "time",
+#'   cens_var = "cens",
+#'   baseline_prior = prior_normal(0, 1000),
+#'   start_times = c(0, 10, 20)
+#' )
+outcome_surv_piecewise_exp <- function(time_var, cens_var, start_times, baseline_prior, weight_var = "") {
+  assert_string(time_var)
+  assert_string(cens_var)
+  assert_string(weight_var)
+  assert_numeric(start_times, lower = 0, finite = TRUE, any.missing = FALSE, min.len = 1, sorted = TRUE)
+  assert_class(baseline_prior, "Prior")
+  has_weight <- isTRUE(weight_var != "")
+  .outcome_surv_piecewise_exp(
+    time_var = time_var,
+    cens_var = cens_var,
+    baseline_prior = baseline_prior,
+    weight_var = weight_var,
+    likelihood_stan_code = "",
+    data_stan_code = h_glue(
+      "
+      vector[N] time;
+      vector[N] cens;
+      int<lower = 1> M = {{length(start_times)}};
+      vector[M] starts = [{{toString(start_times)}}]';
+      {{weight}}",
+      weight = if (has_weight) "vector[N] weight;" else ""
+    ),
+    function_stan_code = h_glue("
+    vector make_durations(vector starts, vector time) {
+      vector [rows(starts)] ends = append_row(tail(starts, rows(starts) - 1), max(time));
+      return fdim(ends, starts);
+    }"),
+    param_stan_code = "vector[M] alpha;",
+    transformed_data_stan_code = h_glue("
+      vector[M] durations;
+      matrix[N,M] T;
+      matrix[N,M] D;
+
+      durations = make_durations(starts, time);
+
+      for(j in 1:M) {
+        T[,j] = fmin(fdim(time, starts[j]), durations[j]);
+      }
+      for(j in 1:M) {
+        for(i in 1:N) {
+          D[i,j] = (starts[j] <= time[i] && time[i] < starts[j] + durations[j]) * (1 - cens[i]);
+        }
+      }")
+  )
+}