Merge pull request #36 from Merck/prepcran

Complete improvement to calc_surv/haz calls throughout package

R/basics.R

@@ -22,6 +22,36 @@
 # basics.R
 # ==================================================================
 
+# Extract the type of a survival object - either "par" (parametric) or "spl" (splines)
+extract_type <- function(survobj) {
+  if (is.null(survobj)) {return(NA)}
+  objtype <- survobj$dlist$name
+  if (objtype=="survspline") {
+      type <- "spl"
+    } else {
+      type <- "par"
+    }
+  return(type)
+}
+
+# Extract the specification of a survival object - dependent on whether parametric or splines model
+extract_spec <- function(survobj) {
+  if (is.null(survobj)) {return(NA)}
+  objtype <- survobj$dlist$name
+  if (objtype=="survspline") {
+    spec <- list(
+      knots=survobj$knots,
+      scale=survobj$scale,
+      gamma=survobj$res[,1])
+  } else {
+    spec <- list(
+      dist=survobj$dlist$name,
+      pars=survobj$res[,1]
+    )
+  }
+  return(spec)
+}
+
 #' Calculate the distribution function for parametric functions
 #' @description Calculate the value of the distribution function for a parametric distribution, given the statistical distribution and its parameters.
 #' @param time is the time at which the distribution function should be calculated.
@@ -109,20 +139,8 @@ calc_pdist_spl <- function(time, spec) {
 calc_pdist <- function(time, type=NA, spec=NA, survobj=NULL){
   # Where a flexsurv object (survobj) is specified, pull its distribution, parameters, type 
   if (!is.null(survobj)) {
-    objtype <- survobj$dlist$name
-    if (objtype=="survspline") {
-      type <- "spl"
-      spec <- list(
-        knots=survobj$knots,
-        scale=survobj$scale,
-        gamma=survobj$res[,1])
-    } else {
-      type <- "par"
-      spec <- list(
-        dist=survobj$dlist$name,
-        pars=survobj$res[,1]
-      )
-    }
+    type <- extract_type(survobj)
+    spec <- extract_spec(survobj)
   }
   # Then call either parametric or splines function, using type and spec parameters
   if (type=="spl") {
@@ -212,20 +230,8 @@ flexsurv::hsurvspline(x = time,
 calc_haz <- function(time, type=NA, spec=NA, survobj=NULL){
   # Where a flexsurv object (survobj) is specified, pull its distribution, parameters, type 
   if (!is.null(survobj)) {
-    objtype <- survobj$dlist$name
-    if (objtype=="survspline") {
-      type <- "spl"
-      spec <- list(
-        knots=survobj$knots,
-        scale=survobj$scale,
-        gamma=survobj$res[,1])
-    } else {
-      type <- "par"
-      spec <- list(
-        dist=survobj$dlist$name,
-        pars=survobj$res[,1]
-      )
-    }
+    type <- extract_type(survobj)
+    spec <- extract_spec(survobj)
   }
   # Then call either parametric or splines function, using type and spec parameters
   if (type=="spl") {
@@ -297,20 +303,8 @@ calc_dens_spl <- function(time, spec) {
 calc_dens <- function(time, type=NA, spec=NA, survobj=NULL){
   # Where a flexsurv object (survobj) is specified, pull its distribution, parameters, type 
   if (!is.null(survobj)) {
-    objtype <- survobj$dlist$name
-    if (objtype=="survspline") {
-      type <- "spl"
-      spec <- list(
-        knots=survobj$knots,
-        scale=survobj$scale,
-        gamma=survobj$res[,1])
-    } else {
-      type <- "par"
-      spec <- list(
-        dist=survobj$dlist$name,
-        pars=survobj$res[,1]
-      )
-    }
+    type <- extract_type(survobj)
+    spec <- extract_spec(survobj)
   }
   # Then call either parametric or splines function, using type and spec parameters
   if (type=="spl") {
@@ -392,8 +386,7 @@ calc_rmd_spl <- function(Tw, spec) {
 #' - `knots` - Vector of locations of knots on the axis of log time, supplied in increasing order. Unlike in [flexsurv::flexsurvspline], these include the two boundary knots.
 #' - `scale` - Either "hazard", "odds", or "normal", as described in [flexsurv::flexsurvspline]. With the default of no knots in addition to the boundaries, this model reduces to the Weibull, log-logistic and log-normal respectively. The scale must be common to all times.
 #' @param survobj is a survival fit object from [flexsurv::flexsurvspline] or [flexsurv::flexsurvreg]
-#' @inherit calc_haz_par seealso
-#' @inherit calc_haz_par return
+#' @return the restricted mean duration, a numeric value.
 #' @export
 #' @examples
 #' calc_rmd(Tw=200,
@@ -407,20 +400,8 @@ calc_rmd_spl <- function(Tw, spec) {
 calc_rmd <- function(Tw, type=NA, spec=NA, survobj=NULL){
   # Where a flexsurv object (survobj) is specified, pull its distribution, parameters, type 
   if (!is.null(survobj)) {
-    objtype <- survobj$dlist$name
-    if (objtype=="survspline") {
-      type <- "spl"
-      spec <- list(
-        knots=survobj$knots,
-        scale=survobj$scale,
-        gamma=survobj$res[,1])
-    } else {
-      type <- "par"
-      spec <- list(
-        dist=survobj$dlist$name,
-        pars=survobj$res[,1]
-      )
-    }
+    type <- extract_type(survobj)
+    spec <- extract_spec(survobj)
   }
   # Then call either parametric or splines function, using type and spec parameters
   if (type=="spl") {

R/discrmd.R

@@ -114,23 +114,19 @@ drmd_psm <- function(ptdata, dpam, psmtype="simple", Ty=10, discrate=0, lifetabl
 # drmd_stm_cf(dpam=params, lifetable=ltable)
 drmd_stm_cf <- function(dpam, Ty=10, discrate=0, lifetable=NA, timestep=1) {
   # Declare local variables
-  Tw <- tvec <- ppd.ts <- ttp.ts <- pps.ts <- NULsppd <- sttp <- sos <- NULL
+  Tw <- tvec <- sppd <- sttp <- sos <- NULL
   adjsppd <- adjos <- vn <- pf <- os <- NULL
   # Time horizon in weeks (ceiling)
   Tw <- convert_yrs2wks(Ty)
   # Create time vector, with half-cycle addition
   tvec <- timestep*(0:floor(Tw/timestep)) + timestep/2
-  # Pull out type and spec for PPD and TTP
-  ppd.ts <- convert_fit2spec(dpam$ppd)
-  ttp.ts <- convert_fit2spec(dpam$ttp)
-  pps.ts <- convert_fit2spec(dpam$pps_cf)
   # Obtain hazard and survival functions
-  hppd <- tvec |> purrr::map_dbl(~calc_haz(.x, ppd.ts$type, ppd.ts$spec))
-  http <- tvec |> purrr::map_dbl(~calc_haz(.x, ttp.ts$type, ttp.ts$spec))
-  hpps <- tvec |> purrr::map_dbl(~calc_haz(.x, pps.ts$type, pps.ts$spec))
-  sppd <- tvec |> purrr::map_dbl(~calc_surv(.x, ppd.ts$type, ppd.ts$spec))
-  sttp <- tvec |> purrr::map_dbl(~calc_surv(.x, ttp.ts$type, ttp.ts$spec))
-  spps <- tvec |> purrr::map_dbl(~calc_surv(.x, pps.ts$type, pps.ts$spec))
+  hppd <- tvec |> purrr::map_dbl(~calc_haz(.x, survobj=dpam$ppd))
+  http <- tvec |> purrr::map_dbl(~calc_haz(.x, survobj=dpam$ttp))
+  hpps <- tvec |> purrr::map_dbl(~calc_haz(.x, survobj=dpam$pps_cf))
+  sppd <- tvec |> purrr::map_dbl(~calc_surv(.x, survobj=dpam$ppd))
+  sttp <- tvec |> purrr::map_dbl(~calc_surv(.x, survobj=dpam$ttp))
+  spps <- tvec |> purrr::map_dbl(~calc_surv(.x, survobj=dpam$pps_cf))
   # Derive the constrained S_PPD and S_PFS
   con.sppd <- constrain_survprob(sppd, lifetable=lifetable, timevec=tvec)
   con.spps <- constrain_survprob(spps, lifetable=lifetable, timevec=tvec)
@@ -171,27 +167,23 @@ drmd_stm_cf <- function(dpam, Ty=10, discrate=0, lifetable=NA, timestep=1) {
 # drmd_stm_cr(dpam=params, lifetable=ltable)
 drmd_stm_cr <- function(dpam, Ty=10, discrate=0, lifetable=NA, timestep=1) {
   # Declare local variables
-  Tw <- tvec <- ppd.ts <- ttp.ts <- sppd <- sttp <- sos <- NULL
+  Tw <- tvec <- sppd <- sttp <- sos <- NULL
   adjsppd <- adjos <- vn <- pf <- os <- NULL
   # Time horizon in weeks (ceiling)
   Tw <- convert_yrs2wks(Ty)
   # Create time vector, with half-cycle addition
   tvec <- timestep*(0:floor(Tw/timestep)) + timestep/2
-  # Pull out type and spec for PPD and TTP
-  ppd.ts <- convert_fit2spec(dpam$ppd)
-  ttp.ts <- convert_fit2spec(dpam$ttp)
-  pps.ts <- convert_fit2spec(dpam$pps_cr)
   # Obtain unconstrained survival functions
-  sppd <- tvec |> purrr::map_dbl(~calc_surv(.x, ppd.ts$type, ppd.ts$spec))
-  sttp <- tvec |> purrr::map_dbl(~calc_surv(.x, ttp.ts$type, ttp.ts$spec))
+  sppd <- tvec |> purrr::map_dbl(~calc_surv(.x, survobj=dpam$ppd))
+  sttp <- tvec |> purrr::map_dbl(~calc_surv(.x, survobj=dpam$ttp))
   # Derive the constrained S_PPD
   c.sppd <- constrain_survprob(sppd, lifetable=lifetable, timevec=tvec)
   # Integrand with constraints on S_PPD and S_PPS
   integrand <- function(u, t) {
-    i.http <- calc_haz(u, ttp.ts$type, ttp.ts$spec)
-    i.sttp <- calc_surv(u, ttp.ts$type, ttp.ts$spec)
-    i.u.sppd <- calc_surv(u, ppd.ts$type, ppd.ts$spec)
-    i.u.spps <- calc_surv(t-u, pps.ts$type, pps.ts$spec)
+    i.http <- calc_haz(u, survobj=dpam$ttp)
+    i.sttp <- calc_surv(u, survobj=dpam$ttp)
+    i.u.sppd <- calc_surv(u, survobj=dpam$ppd)
+    i.u.spps <- calc_surv(t-u, survobj=dpam$pps_cr)
     i.slxu <- calc_ltsurv(convert_wks2yrs(u), lifetable=lifetable)
     i.slxt <- calc_ltsurv(convert_wks2yrs(t), lifetable=lifetable)
     i.c.sppd <- pmin(i.u.sppd, i.slxu)

R/fitting.R

@@ -227,9 +227,11 @@ fit_ends_mods_par <- function(simdat,
 #' # Parametric modeling
 #' fits_par <- fit_ends_mods_par(bosonc)
 #' find_bestfit(fits_par$ttp, "aic")
+#' \donttest{
 #' # Splines modeling
 #' fits_spl <- fit_ends_mods_spl(bosonc)
 #' find_bestfit(fits_spl$ttp, "bic")
+#' }
 find_bestfit <- function(reglist, crit) {
   # Declare local variables
   fittable <- chosen <- NULL

R/lhoods.R

@@ -21,36 +21,6 @@
 # lhoods.R
 # ==================================================================
 
-#' Obtain the type and specification as required in other package functions from a model fit
-#' @param fitsurv Fitted model (either parametric or Royston-Parmar splines model)
-#' @return List of type and specification
-#' - `type` is "spl" for splines model or "par" for parametric model
-#' - `spec` contains distribution (`dist`) and coefficients (`coefs`) if `type=="par"`
-#' - `spec` contains gamma values (`gamma`), knot locations (log scale, `knots`) and scale (`scale`) for Royston-Parmar splines model, if `type=="spl"`
-#' @seealso [flexsurv::flexsurvspline()]
-#' @noRd
-# Examples
-# bosonc <- create_dummydata("flexbosms")
-# fits <- fit_ends_mods_spl(bosonc)
-# convert_fit2spec(fits$pfs[[3]]$result)
-convert_fit2spec <- function(fitsurv) {
-  # Declare local variables
-  par.dist <- type <- spec <- NULL
-  # Pick out distribution/splines
-  par.dist <- fitsurv$dlist$name
-  if (par.dist=="survspline") {
-    type <- "spl"
-    spec <- list(gamma = fitsurv$coefficients,
-                 knots = fitsurv$aux$knots,
-                 k = length(fitsurv$aux$knots)-2,
-                 scale = fitsurv$aux$scale)
-  } else {
-    type <- "par"
-    spec <- list(dist=par.dist, pars=fitsurv$res[,1])
-  }
-  return(list(type=type, spec=spec))
-}
-
 #' Calculate likelihood for a simple three-state partitioned survival model
 #' @description Calculate likelihood values and other summary output for a simple three-state partitioned survival model, given appropriately formatted patient-level data, a set of fitted survival regressions, and the time cut-off (if two-piece modeling is used). This function is called by [calc_likes].x three-state partitioned survival model, given appropriately formatted patient-level data, a set of fitted survival regressions, and the time cut-off (if two-piece modeling is used). This function is called by [calc_likes]. Unlike [calc_likes_psm_complex], this likelihood function assumes a progression hazard can be derived from the PFS hazard function and the ratio of progression to PFS events from PF patients.
 #' @param ptdata Dataset of patient level data. Must be a tibble with columns named:

R/ppdpps.R

@@ -62,29 +62,17 @@
 #' }
 calc_haz_psm <- function(timevar, ptdata, dpam, psmtype) {
   # Declare local variables
-  pfs.ts <- pfs.type <- pfs.spec <- NULL
-  os.ts <- os.type <- os.spec <- NULL
-  ttp.ts <- ttp.type <- ttp.spec <- NULL
   ne_pfs <- ne_ttp <- progfrac <- NULL
   http <- hpf <- hos <- sos <- spf <- NULL
   hppd_simple <- hppd_complex <- hppd <- hpps <- hdiff <- NULL
   # PFS
-  pfs.ts <- convert_fit2spec(dpam$pfs)
-  pfs.type <- pfs.ts$type
-  pfs.spec <- pfs.ts$spec
-  hpf <- calc_haz(timevar, pfs.type, pfs.spec)
-  spf <- calc_surv(timevar, pfs.type, pfs.spec)
+  hpf <- calc_haz(timevar, survobj=dpam$pfs)
+  spf <- calc_surv(timevar, survobj=dpam$pfs)
   # OS
-  os.ts <- convert_fit2spec(dpam$os)
-  os.type <- os.ts$type
-  os.spec <- os.ts$spec
-  hos <- calc_haz(timevar, os.type, os.spec)
-  sos <- calc_surv(timevar, os.type, os.spec)
+  hos <- calc_haz(timevar, survobj=dpam$os)
+  sos <- calc_surv(timevar, survobj=dpam$os)
   # TTP complex
-  ttp.ts <- convert_fit2spec(dpam$ttp)
-  ttp.type <- ttp.ts$type
-  ttp.spec <- ttp.ts$spec
-  http_complex <- calc_haz(timevar, ttp.type, ttp.spec)
+  http_complex <- calc_haz(timevar, survobj=dpam$ttp)
   # TTP simple
   ne_pfs <- sum(ptdata$pfs.flag)
   ne_ttp <- sum(ptdata$ttp.flag)
@@ -102,8 +90,6 @@ calc_haz_psm <- function(timevar, ptdata, dpam, psmtype) {
   hpps_unadj <- (sos*hos-spf*hppd)/(sos-spf)
   hpps <- pmax(0, pmin(hpps_unadj, 5000))
   hpps[timevar==0] <- 0
-  # Diff
-  # hdiff <- hos-(spf*hppd+(sos-spf)*hpps)
   # Adjusted for caps and collars
   hadj <- list(
     ttp = http,