418-population-effects-implementation

NAMESPACE

@@ -251,6 +251,7 @@ export(linkNone)
 export(linkShrinkage)
 export(linkTTG)
 export(merge)
+export(populationHR)
 export(prior_beta)
 export(prior_cauchy)
 export(prior_gamma)

R/PopulationHR.R

@@ -1,83 +1,85 @@
 #' Calculate Population Hazard Ratios
 #'
-#' Calculates hazard ratios marginalised over subject specific random effects
+#' Calculates hazard ratios marginalised over subject specific random effects using the
+#' approach proposed by \insertCite{oudenhoven2020marginal}{jmpost}.
 #'
 #' @param object ([`JointModelSamples`]) \cr samples as drawn from a Joint Model.
-#' @param HR_terms (`character`) \cr to include in the hazard ratio calculation.
-#' By default this include the terms used in the survival model.
+#' @param HR_formula (`formula`) \cr defines the terms to include in the hazard ratio calculation.
+#' By default this uses the right side of the formula used in the survival model.
 #' Set to `NULL` not include any terms
-#' @param arm (`logical`) \cr calculate hazard ratio for `arm`.
-#' @param baseline (`string`) \cr terms to model baseline hazard using variable `"time"`.
-#' Default is a B-spline from [splines]: `"bs(time, df=10)"`
+#' @param baseline (`formula`) \cr terms to model baseline hazard using variable `time`.
+#' Default is a B-spline from [splines]: `~bs(time, df = 10)`
+#' @param quantiles (`numeric`) \cr vector of two values in (0, 1) for calculating quantiles from log hazard ratio
+#' distributions.
 #'
-#' @returns A matrix
+#' @references \insertAllCited{}
 #'
+#' @returns A matrix
+#' @export
 #' @importFrom splines bs
+populationHR <- function(
+    object,
+    HR_formula = object@data@survival@formula,
+    baseline = ~ bs(time, df = 10),
+    quantiles = c(0.025, 0.975)
+) {
+    assert_class(object, "JointModelSamples")
+    assert_formula(HR_formula)
+    assert_formula(baseline)
+    assert_numeric(quantiles, lower = 0, upper = 1, any.missing = FALSE, unique = TRUE)
+    if (!"time" %in% all.vars(baseline)) stop("baseline formula should include a time term.")
 
-PopulationHR <- function(object,
-                         HR_terms = attr(terms(object@data@survival@formula), "term.labels"),
-                         arm = TRUE,
-                         baseline = "bs(time, df=10)") {
     # Extract the variable names used in the data
     subject_var <- object@data@subject@subject
     arm_var <- object@data@subject@arm
     long_time_var <- all.vars(delete.response(terms(object@data@longitudinal@formula)))
     surv_time_var <- all.vars(object@data@survival@formula[[2]][[2]])
-    surv_covs <- all.vars(delete.response(terms(object@data@survival@formula)))
+    surv_covs <- all.vars(delete.response(terms(HR_formula)))
+    if (!all(surv_covs %in% colnames(object@data@survival@data))) {
+        stop("All variables in HR_formula must be in survival data")
+    }
 
-    population_terms <- unique(c(
-        baseline,
-        if (arm) arm_var else NULL,
-        HR_terms
-    ))
 
-    marginal_formula <- reformulate(population_terms)
+    marginal_formula <- stats::reformulate(c(
+        attr(terms(baseline), "term.labels"),
+        attr(terms(HR_formula), "term.labels")
+    ))
 
-    # Get the survival quantities at the observed longitudinal times
-    # and survival event times for each patient:
+    # Get the survival quantities at the observed longitudinal and survival times for each patient:
     times_df <- rbind(
-        object@data@longitudinal@data |>
-            select(subject = {{subject_var}}, time = {{long_time_var}}),
-        object@data@survival@data |>
-            select(subject = {{subject_var}}, time = {{surv_time_var}})
-    ) |>
-        filter(time > 0)
+        stats::setNames(object@data@longitudinal@data[, c(subject_var, long_time_var)], c("subject", "time")),
+        stats::setNames(object@data@survival@data[, c(subject_var, surv_time_var)], c("subject", "time"))
+    )
+    times_df <- times_df[order(times_df$subject, times_df$time), ]
+    times_df <- times_df[times_df$time > 0, ]
+    times_df <- times_df[!duplicated.data.frame(times_df), ]
 
     # Generate samples of the log hazard log h_i(t) for each patient i at these time points t.
     grid_spec <- split(times_df$time, times_df$subject)
-    logH <- SurvivalQuantities(
+    log_haz_samples <- SurvivalQuantities(
         object,
         grid = GridManual(grid_spec),
         type = "loghaz"
     )@quantities@quantities |> t()
 
     # Construct \tilde{X} in paper's notation with one row per patient's time
-    W_df <- times_df |>
-        select(subject, time) |>
-        left_join(
-            object@data@subject@data |>
-                select(subject = {{subject_var}}, arm = {{arm_var}}),
-            by = "subject"
-        ) |>
-        left_join(
-            object@data@survival@data |>
-                select(subject = {{subject_var}}, {{surv_covs}}),
-            by = "subject"
-        )
-
+    W_df <- dplyr::left_join(
+        times_df,
+        stats::setNames(object@data@survival@data[, c(subject_var, surv_covs)], c("subject", surv_covs)),
+        by = "subject"
+    )
     W_mat <- model.matrix(marginal_formula, W_df)
+
     # As model matrix contains baseline and covariates, we don't need the intercept term
     # but we want factor variables encoded relative to the intercept
     W_mat <- W_mat[, colnames(W_mat) != "(Intercept)", drop = FALSE]
 
-    estimates <- lm.fit(x = W_mat, y = logH)$coefficients
+    estimates <- stats::lm.fit(x = W_mat, y = log_haz_samples)$coefficients
     tidy_res <- apply(estimates, 1, function(x) {
-        c(mean = mean(x),
-          median = median(x),
-          quantile(x, c(0.025, 0.975))
-        )
-        }) |> t()
+        quantiles <- stats::quantile(x, probs = quantiles)
+        c(mean = mean(x), median = median(x), quantiles)
+    }) |>
+        t()
 
     list(tidy_res, estimates)
 }
-

_pkgdown.yml

@@ -28,14 +28,14 @@ reference:
     contents:
       - jmpost-package
       - jmpost-settings
-  
+
   - title: Data Specification
     contents:
       - DataJoint
       - DataSubject
       - DataLongitudinal
       - DataSurvival
-  
+
   - title: Data Simulation
     contents:
       - SimJointData
@@ -65,24 +65,24 @@ reference:
       - prior_std_normal
       - prior_student_t
       - prior_init_only
-  
+
   - title: Longitudinal Model Specification
     contents:
       - LongitudinalRandomSlope
       - LongitudinalGSF
       - LongitudinalSteinFojo
       - LongitudinalClaretBruno
       - LongitudinalModel
-  
+
   - title: Survival Model Specification
     contents:
       - SurvivalExponential
       - SurvivalLogLogistic
       - SurvivalModel
       - SurvivalWeibullPH
       - SurvivalGamma
-  
-  
+
+
   - title: Link Specification
     contents:
       - linkDSLD
@@ -94,16 +94,16 @@ reference:
       - Link
       - LinkComponent
       - standard-link-user
-  
+
   - title: Joint Model Specification
     contents:
       - JointModel
-  
+
   - title: Joint Model Fitting
     contents:
       - compileStanModel
       - sampleStanModel
-  
+
   - title: Postprocessing
     contents:
       - JointModelSamples
@@ -119,7 +119,8 @@ reference:
       - autoplot.SurvivalQuantities
       - brierScore
       - brierScore.SurvivalQuantities
-
+      - populationHR
+
   - title: Stan Code
     contents:
       - Parameter
@@ -133,7 +134,7 @@ reference:
       - write_stan
       - getParameters
       - getRandomEffectsNames
-  
+
   - title: Convenience Functions
     contents:
       - initialValues
@@ -185,7 +186,7 @@ reference:
       - as_formula
       - getPredictionNames
       - set_limits
-  
+
   - title: Promises
     contents:
       - Promise
@@ -194,8 +195,8 @@ reference:
       - resolvePromise
       - resolvePromise.Link
       - resolvePromise.PromiseLinkComponent
-  
-  - title: Miscellaneous 
+
+  - title: Miscellaneous
     contents:
       - STAN_BLOCKS
 

inst/REFERENCES.bib

@@ -37,3 +37,17 @@ @article{beyer2020multistate
   year={2020},
   publisher={Wiley Online Library}
 }
+
+@article{oudenhoven2020marginal,
+  author = {van Oudenhoven, Floor M. and Swinkels, Sophie H. N. and Ibrahim, Joseph G. and Rizopoulos, Dimitris},
+  title = {A marginal estimate for the overall treatment effect on a survival outcome within the joint modeling framework},
+  journal = {Statistics in Medicine},
+  volume = {39},
+  number = {28},
+  pages = {4120-4132},
+  keywords = {joint models, marginal estimates, overall treatment effect, survival outcome.},
+  doi = {https://doi.org/10.1002/sim.8713},
+  url = {https://onlinelibrary.wiley.com/doi/abs/10.1002/sim.8713},
+  eprint = {https://onlinelibrary.wiley.com/doi/pdf/10.1002/sim.8713},
+  year = {2020}
+}