Simulate patients from fitted model

DESCRIPTION

@@ -27,7 +27,7 @@ Language: en-GB
 URL: https://genentech.github.io/jmpost/
 BugReports: https://github.com/Genentech/jmpost/issues
 Roxygen: list(markdown = TRUE)
-RoxygenNote: 7.3.2
+RoxygenNote: 7.3.3
 Depends:
     R (>= 4.1.0)
 Imports:
@@ -120,13 +120,16 @@ Collate:
     'SurvivalLoglogistic.R'
     'SurvivalWeibullPH.R'
     'brier_score.R'
+    'data.R'
     'defaults.R'
     'external-exports.R'
     'jmpost-package.R'
     'link_generics.R'
     'populationHR.R'
     'settings.R'
+    'simulate.R'
     'standalone-s3-register.R'
     'zzz.R'
 VignetteBuilder: knitr
 RdMacros: Rdpack
+LazyData: true

NAMESPACE

@@ -87,6 +87,14 @@ S3method(coalesceGridTime,GridPopulation)
 S3method(coalesceGridTime,GridPrediction)
 S3method(coalesceGridTime,default)
 S3method(compileStanModel,JointModel)
+S3method(createLongitudinalSimObject,LongitudinalClaretBruno)
+S3method(createLongitudinalSimObject,LongitudinalGSF)
+S3method(createLongitudinalSimObject,LongitudinalRandomSlope)
+S3method(createLongitudinalSimObject,LongitudinalSteinFojo)
+S3method(createSurvivalSimObject,SurvivalExponential)
+S3method(createSurvivalSimObject,SurvivalGamma)
+S3method(createSurvivalSimObject,SurvivalLogLogistic)
+S3method(createSurvivalSimObject,SurvivalWeibullPH)
 S3method(dim,Quantities)
 S3method(enableGQ,JointModel)
 S3method(enableGQ,LongitudinalClaretBruno)
@@ -172,6 +180,7 @@ S3method(sampleSubjects,SimLongitudinalSteinFojo)
 S3method(sampleSubjects,SimSurvival)
 S3method(saveObject,JointModelSamples)
 S3method(set_limits,Prior)
+S3method(simulate,JointModelSamples)
 S3method(size,Parameter)
 S3method(size,ParameterList)
 S3method(subset,DataJoint)
@@ -230,13 +239,15 @@ export(SurvivalLogLogistic)
 export(SurvivalModel)
 export(SurvivalQuantities)
 export(SurvivalWeibullPH)
+export(add_pfs)
 export(as.QuantityCollapser)
 export(as.QuantityGenerator)
 export(as_formula)
 export(as_stan_list)
 export(autoplot)
 export(brierScore)
 export(compileStanModel)
+export(cut_data)
 export(enableGQ)
 export(enableLink)
 export(generateQuantities)
@@ -334,6 +345,7 @@ importFrom(stats,rlogis)
 importFrom(stats,rnorm)
 importFrom(stats,rt)
 importFrom(stats,runif)
+importFrom(stats,simulate)
 importFrom(stats,terms)
 importFrom(survival,Surv)
 importFrom(survival,coxph)

R/SimJointData.R

@@ -143,3 +143,108 @@
         return(object)
     }
 )
+
+
+#' Add PFS events at Tumour Progression to Data
+#'
+#' Adds new columns `pfs_time` and `pfs_event` based on observed changes to SLD.
+#'
+#' @param object A [SimJointData] object
+#' @param relative_threshold (`number`)\cr a multiplicative threshold for the change in SLD compared to the `min(SLD)`.
+#'  Default is 1.2 meaning a 20% increase.
+#' @param absolute_threshold (`number`)\cr an absolute threshold for the change in SLD compared to the minimum.
+#'   Default is 5.
+#' @param from_time (`number`)\cr Ignore observations before this time for determining SLD minimum.
+#' @param observed_after (`logical`)\cr If `FALSE` set longitudinal observations after the progression time to
+#'   `observed = FALSE`
+#' @details
+#' Both thresholds must be met for a progression to be declared.
+#'
+#' @export
+#' @examples
+#' data <- SimJointData(survival = SimSurvivalExponential(lambda = 1/10), longitudinal = SimLongitudinalSteinFojo())
+#' data <- add_pfs(data)
+#' data@survival # now has pfs_time and pfs_event columns
+add_pfs <- function(object, relative_threshold = 1.2, absolute_threshold = 5, from_time = 0, observed_after = FALSE) {
+    assert_class(object, "SimJointData")
+
+    pd_times <- object@longitudinal |>
+        dplyr::filter(.data$time >= from_time) |>
+        dplyr::mutate(
+            min_sld = cummin(.data$sld),
+            is_pd = .data$sld >= pmax(
+                .data$min_sld * relative_threshold,
+                .data$min_sld + absolute_threshold
+                ) & .data$observed,
+            pd_time = min(.data$time[.data$is_pd], Inf),
+            .by = subject
+        ) |>
+        dplyr::select("subject", "pd_time") |>
+        dplyr::slice_head(by = "subject")
+
+    if (isFALSE(observed_after)) {
+        object@longitudinal <- object@longitudinal |>
+            dplyr::left_join(pd_times, by = "subject") |>
+            dplyr::mutate(
+                observed = dplyr::if_else(.data$time > .data$pd_time, FALSE, .data$observed),
+                pd_time = NULL
+            )
+    }
+
+    object@survival <-
+        object@survival |>
+        dplyr::left_join(pd_times, by = "subject") |>
+        dplyr::mutate(
+            pfs_time = pmin(.data$time, .data$pd_time, na.rm = TRUE),
+            pfs_event = dplyr::if_else(.data$pfs_time < .data$time, 1, .data$event),
+            pd_time = NULL
+        )
+    object
+}
+
+
+#' Cut Study Data
+#' @param object A [SimJointData] object
+#' @param cut_time (`numeric`)\cr A vector of cut off times, either length 1 for all patients or
+#'   `nrow(object@survival)` for a time per patient.
+#' @details
+#'   All observations after this time are remove. Survival is censored at this time and any longitudinal
+#'   values are removed.
+#' @export
+#' @examples
+#' data <- SimJointData(survival = SimSurvivalExponential(lambda = 1/10), longitudinal = SimLongitudinalSteinFojo())
+#' data <- cut_data(data, 5)
+#' data@survival
+#' # Now max time is 5
+#' max(data@survival$time)
+cut_data <- function(object, cut_time) {
+    assert_class(object, "SimJointData")
+    check_len <- if (length(cut_time) > 1) nrow(object@survival) else 1
+    assert_numeric(cut_time, lower = 0, len = check_len)
+
+    object@survival <- object@survival |>
+        dplyr::mutate(
+            cut_time = cut_time,
+            event = dplyr::if_else(time < cut_time, event, 0),
+            time = pmin(cut_time, time)
+        )
+    if (all(c("pfs_event", "pfs_time") %in% colnames(object@survival))) {
+        object@survival <- object@survival |>
+            dplyr::mutate(
+                pfs_event = dplyr::if_else(pfs_time < cut_time, pfs_event, 0),
+                pfs_time = pmin(cut_time, pfs_time)
+            )
+    }
+
+    object@longitudinal <- object@longitudinal |>
+        dplyr::left_join(
+            dplyr::select(object@survival, subject, cut_time),
+            by = "subject"
+        ) |>
+        dplyr::filter(time <= cut_time) |>
+        dplyr::mutate(cut_time = NULL)
+
+    object@survival$cut_time <- NULL
+
+    object
+}

R/SimLongitudinalClaretBruno.R

@@ -163,12 +163,13 @@ as_print_string.SimLongitudinalClaretBruno <- function(object) {
 #' @export
 sampleObservations.SimLongitudinalClaretBruno <- function(object, times_df) {
     times_df |>
-        dplyr::mutate(mu_sld = clbr_sld(.data$time, .data$ind_b, .data$ind_g, .data$ind_c, .data$ind_p)) |>
-        dplyr::mutate(dsld = clbr_dsld(.data$time, .data$ind_b, .data$ind_g, .data$ind_c, .data$ind_p)) |>
-        dplyr::mutate(ttg = clbr_ttg(.data$time, .data$ind_b, .data$ind_g, .data$ind_c, .data$ind_p)) |>
-        dplyr::mutate(sld_sd = ifelse(object@scaled_variance, .data$mu_sld * object@sigma, object@sigma)) |>
-        dplyr::mutate(sld = stats::rnorm(dplyr::n(), .data$mu_sld, .data$sld_sd)) |>
         dplyr::mutate(
+            mu_sld = clbr_sld(.data$time, .data$ind_b, .data$ind_g, .data$ind_c, .data$ind_p),
+            dsld = clbr_dsld(.data$time, .data$ind_b, .data$ind_g, .data$ind_c, .data$ind_p),
+            ttg = clbr_ttg(.data$time, .data$ind_b, .data$ind_g, .data$ind_c, .data$ind_p),
+            sld_sd = ifelse(object@scaled_variance, .data$mu_sld * object@sigma, object@sigma),
+            sld = stats::rnorm(dplyr::n(), .data$mu_sld, .data$sld_sd),
+
             log_haz_link =
                 (object@link_dsld * .data$dsld) +
                 (object@link_ttg * .data$ttg) +

R/SimLongitudinalRandomSlope.R

@@ -67,10 +67,10 @@ as_print_string.SimLongitudinalRandomSlope <- function(object) {
 #' @export
 sampleObservations.SimLongitudinalRandomSlope <- function(object, times_df) {
     times_df |>
-        dplyr::mutate(err = stats::rnorm(dplyr::n(), 0, object@sigma)) |>
-        dplyr::mutate(sld_mu = .data$intercept + .data$slope_ind * .data$time) |>
-        dplyr::mutate(sld = .data$sld_mu + .data$err) |>
         dplyr::mutate(
+            err = stats::rnorm(dplyr::n(), 0, object@sigma),
+            sld_mu = .data$intercept + .data$slope_ind * .data$time,
+            sld = .data$sld_mu + .data$err,
             log_haz_link =
                 object@link_dsld * .data$slope_ind +
                 object@link_identity * .data$sld_mu
@@ -86,13 +86,13 @@ sampleSubjects.SimLongitudinalRandomSlope <- function(object, subjects_df) {
     )
 
     assert_that(
-        length(object@slope_mu) == length(unique(subjects_df[["arm"]])),
-        msg = "`length(slope_mu)` should be equal to the number of unique arms"
+        length(object@slope_mu) == nlevels(subjects_df[["arm"]]),
+        msg = "`length(slope_mu)` should be equal to the number of arms"
     )
 
     assert_that(
-        length(object@intercept) == length(unique(subjects_df[["study"]])),
-        msg = "`length(intercept)` should be equal to the number of unique studies"
+        length(object@intercept) == nlevels(subjects_df[["study"]]),
+        msg = "`length(intercept)` should be equal to the number of studies"
     )
 
     assert_that(

R/SimSurvival.R

@@ -43,6 +43,7 @@ NULL
 #' @slot beta_cat (`numeric`)\cr See arguments.
 #' @slot loghazard (`function`)\cr See arguments.
 #' @slot name (`character`)\cr See arguments.
+#' @slot beta_os_cov (`numeric`) Additional regression coefficients for survival models.
 #'
 #' @family SimSurvival
 #' @exportClass SimSurvival
@@ -55,6 +56,7 @@ NULL
         lambda_censor = "numeric",
         beta_cont = "numeric",
         beta_cat = "numeric",
+        beta_os_cov = "numeric",
         loghazard = "function",
         name = "character"
     )
@@ -186,7 +188,8 @@ sampleObservations.SimSurvival <- function(object, times_df) {
         dplyr::mutate(time = dplyr::if_else(.data$real_time <= .data$time_cen, .data$real_time, .data$time_cen)) |>
         dplyr::arrange(.data$subject)
 
-    os_dat_complete[, c("subject", "study", "arm", "time", "event", "cov_cont", "cov_cat")]
+    keep_cols <- colnames(os_dat_complete) %in% c("subject", "study", "arm", "time", "event", "cov_cont", "cov_cat")
+    os_dat_complete[, keep_cols]
 }
 
 

R/data.R

@@ -0,0 +1,64 @@
+#' OAK Study Tumor Growth Data
+#'
+#' Longitudinal tumor size measurements from the OAK clinical trial, a
+#' Phase III study comparing atezolizumab (MPDL3280A) with docetaxel in
+#' patients with previously treated non-small cell lung cancer (NSCLC).
+#' Patients were matched to those in [os_data] using baseline SLD, treatment
+#' arm, and best overall response. Only patients with at least 3 tumor size
+#' measurements are included.
+#'
+#' @format A `tibble` with one row per tumor measurement and 4 columns:
+#' \describe{
+#'   \item{id}{Patient identifier (factor), matched to [os_data]}
+#'   \item{year}{Time of measurement in years from treatment start (numeric)}
+#'   \item{sld}{Sum of longest diameters of target lesions in mm (numeric)}
+#'   \item{arm}{Treatment arm (factor): `"Docetaxel"` or `"MPDL3280A"`}
+#' }
+#' @source Tumor size measurements extracted from the S1 dataset of:
+#'   Ghaffari Laleh N, Loeffler CML, Grajek J, Staňková K, Pearson AT,
+#'   Muti HS, et al. (2022). Classical mathematical models for prediction of
+#'   response to chemotherapy and immunotherapy.
+#'   *PLoS Comput Biol* 18(2): e1009822.
+#'   <https://doi.org/10.1371/journal.pcbi.1009822>
+#'
+#'   Data preparation code adapted from:
+#'   Sabanes Bove D and Mercier F.
+#'   <https://rconis.github.io/tgi-os-training/>
+#' @seealso [os_data]
+"tumor_data"
+
+#' OAK Study Overall Survival Data
+#'
+#' Baseline covariates and survival outcomes from the OAK clinical trial, a
+#' Phase III study comparing atezolizumab (MPDL3280A) with docetaxel in
+#' patients with previously treated non-small cell lung cancer (NSCLC).
+#' Patients were matched to those in [tumor_data] using baseline SLD, treatment
+#' arm, and best overall response.
+#'
+#' @format A `tibble` with one row per patient and 12 columns:
+#' \describe{
+#'   \item{id}{Patient identifier (factor), matched to [tumor_data]}
+#'   \item{arm}{Treatment arm (factor): `"Docetaxel"` or `"MPDL3280A"`}
+#'   \item{ecog}{ECOG performance status at baseline (factor): `0` or `1`}
+#'   \item{age}{Age at baseline in years (numeric)}
+#'   \item{race}{Patient race (factor)}
+#'   \item{sex}{Patient sex (factor): `"F"` or `"M"`}
+#'   \item{sld}{Baseline sum of longest diameters of target lesions in mm (numeric)}
+#'   \item{response}{Best overall response (factor): `"CR"`, `"PR"`, `"SD"`, `"PD"`, or `"NE"`}
+#'   \item{pfs_time}{Progression-free survival time in years (numeric)}
+#'   \item{pfs_event}{Progression-free survival event indicator (logical): `TRUE` = event observed}
+#'   \item{os_time}{Overall survival time in years (numeric)}
+#'   \item{os_event}{Overall survival event indicator (logical): `TRUE` = death observed}
+#' }
+#' @source Supplementary Table 8 from:
+#'   Rittmeyer A, Barlesi F, Waterkamp D, et al. (2017). Atezolizumab versus
+#'   docetaxel in patients with previously treated non-small-cell lung cancer
+#'   (OAK): a phase 3, open-label, multicentre randomised controlled trial.
+#'   *The Lancet* 389(10066): 255–265.
+#'   <https://doi.org/10.1038/s41591-018-0134-3>
+#'
+#'   Data preparation code adapted from:
+#'   Sabanes Bove D and Mercier F.
+#'   <https://rconis.github.io/tgi-os-training/>
+#' @seealso [tumor_data]
+"os_data"