Simulate patients from fitted model

.lintr

@@ -3,7 +3,6 @@ linters: linters_with_defaults(
     line_length_linter(120),
     object_name_linter = NULL,     # Stops complaints about functions not being in snake_case
     object_usage_linter = NULL,    # Stops spurious warnings over imported functions
-    cyclocomp_linter = NULL,
     indentation_linter = indentation_linter(4),
     object_length_linter = object_length_linter(40),
     return_linter = NULL

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/Quantities.R

@@ -146,7 +146,7 @@ collapse_quantities <- function(quantities_raw, collapser) {
         ncol = length(collapser)
     )
 
-    for (idx in seq_len(length(collapser))) {
+    for (idx in seq_along(collapser)) {
         quantities[, idx] <- quantities_raw[
             ,
             collapser@indexes[[idx]],

R/SimJointData.R

@@ -143,3 +143,114 @@ setMethod(
         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(.data$time < .data$cut_time, .data$event, 0),
+            time = pmin(.data$cut_time, .data$time)
+        )
+    if (all(c("pfs_event", "pfs_time") %in% colnames(object@survival))) {
+        object@survival <- object@survival |>
+            dplyr::mutate(
+                pfs_event = dplyr::if_else(.data$pfs_time < .data$cut_time, .data$pfs_event, 0),
+                pfs_time = pmin(.data$cut_time, .data$pfs_time)
+            )
+    }
+
+    object@longitudinal <- object@longitudinal |>
+        dplyr::left_join(
+            dplyr::select(object@survival, "subject", "cut_time"),
+            by = "subject"
+        ) |>
+        dplyr::filter(.data$time <= .data$cut_time) |>
+        dplyr::mutate(cut_time = NULL)
+
+    object@survival$cut_time <- NULL
+
+    object
+}

R/SimLongitudinal.R

@@ -36,6 +36,6 @@ setMethod(
 )
 
 #' @rdname as_print_string
-as_print_string.SimLongitudinal <- function(object) {
+as_print_string.SimLongitudinal <- function(object, ...) {
     return("SimLongitudinal")
 }

R/SimLongitudinalClaretBruno.R

@@ -155,20 +155,21 @@ setValidity(
 )
 
 #' @rdname as_print_string
-as_print_string.SimLongitudinalClaretBruno <- function(object) {
+as_print_string.SimLongitudinalClaretBruno <- function(object, ...) {
     return("SimLongitudinalClaretBruno")
 }
 
 #' @rdname sampleObservations
 #' @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/SimLongitudinalGSF.R

@@ -156,20 +156,20 @@ setValidity(
 )
 
 #' @rdname as_print_string
-as_print_string.SimLongitudinalGSF <- function(object) {
+as_print_string.SimLongitudinalGSF <- function(object, ...) {
     return("SimLongitudinalGSF")
 }
 
 #' @rdname sampleObservations
 #' @export
 sampleObservations.SimLongitudinalGSF <- function(object, times_df) {
     times_df |>
-        dplyr::mutate(mu_sld = gsf_sld(.data$time, .data$psi_b, .data$psi_s, .data$psi_g, .data$psi_phi)) |>
-        dplyr::mutate(dsld = gsf_dsld(.data$time, .data$psi_b, .data$psi_s, .data$psi_g, .data$psi_phi)) |>
-        dplyr::mutate(ttg = gsf_ttg(.data$time, .data$psi_b, .data$psi_s, .data$psi_g, .data$psi_phi)) |>
-        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 = gsf_sld(.data$time, .data$psi_b, .data$psi_s, .data$psi_g, .data$psi_phi),
+            dsld = gsf_dsld(.data$time, .data$psi_b, .data$psi_s, .data$psi_g, .data$psi_phi),
+            ttg = gsf_ttg(.data$time, .data$psi_b, .data$psi_s, .data$psi_g, .data$psi_phi),
+            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) +
@@ -194,30 +194,31 @@ sampleSubjects.SimLongitudinalGSF <- function(object, subjects_df) {
 
     res <- subjects_df |>
         dplyr::distinct(.data$subject, .data$arm, .data$study) |>
-        dplyr::mutate(study_idx = as.numeric(.data$study)) |>
-        dplyr::mutate(arm_idx = as.numeric(.data$arm)) |>
-        dplyr::mutate(psi_b = stats::rlnorm(
-            dplyr::n(),
-            object@mu_b[.data$study_idx],
-            object@omega_b[.data$study_idx]
-        )) |>
-        dplyr::mutate(psi_s = stats::rlnorm(
-            dplyr::n(),
-            object@mu_s[.data$arm_idx],
-            object@omega_s[.data$arm_idx]
-        )) |>
-        dplyr::mutate(psi_g = stats::rlnorm(
-            dplyr::n(),
-            object@mu_g[.data$arm_idx],
-            object@omega_g[.data$arm_idx]
-        )) |>
-        dplyr::mutate(psi_phi_logit = stats::rnorm(
-            dplyr::n(),
-            object@mu_phi[.data$arm_idx],
-            object@omega_phi[.data$arm_idx]
-        )) |>
-        dplyr::mutate(psi_phi = stats::plogis(.data$psi_phi_logit))
-
+        dplyr::mutate(
+            study_idx = as.numeric(.data$study),
+            arm_idx = as.numeric(.data$arm),
+            psi_b = stats::rlnorm(
+                dplyr::n(),
+                object@mu_b[.data$study_idx],
+                object@omega_b[.data$study_idx]
+            ),
+            psi_s = stats::rlnorm(
+                dplyr::n(),
+                object@mu_s[.data$arm_idx],
+                object@omega_s[.data$arm_idx]
+            ),
+            psi_g = stats::rlnorm(
+                dplyr::n(),
+                object@mu_g[.data$arm_idx],
+                object@omega_g[.data$arm_idx]
+            ),
+            psi_phi_logit = stats::rnorm(
+                dplyr::n(),
+                object@mu_phi[.data$arm_idx],
+                object@omega_phi[.data$arm_idx]
+            ),
+            psi_phi = stats::plogis(.data$psi_phi_logit)
+        )
     res[, c("subject", "arm", "study", "psi_b", "psi_s", "psi_g", "psi_phi")]
 }
 

R/SimLongitudinalRandomSlope.R

@@ -59,18 +59,18 @@ SimLongitudinalRandomSlope <- function(
 }
 
 #' @rdname as_print_string
-as_print_string.SimLongitudinalRandomSlope <- function(object) {
+as_print_string.SimLongitudinalRandomSlope <- function(object, ...) {
     return("SimLongitudinalRandomSlope")
 }
 
 #' @rdname sampleObservations
 #' @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(