Revise application of liftable constraints in STM and update vignette

DESCRIPTION

@@ -46,7 +46,6 @@ Collate:
     'basics.R'
     'brier.R'
     'datasets.R'
-    'discrmd.R'
     'fitting-spl.R'
     'fitting.R'
     'lhoods.R'
@@ -55,3 +54,4 @@ Collate:
     'probgraphs.R'
     'psm3mkv-package.R'
     'resmeans.R'
+    'discrmd.R'

R/discrmd.R

@@ -24,7 +24,18 @@
 
 #' Discretized Restricted Mean Duration calculation for Partitioned Survival Model
 #' Calculate restricted mean duration (RMD) in PF, PD and OS states under a Partitioned Survival Model structure.
+#' @param ptdata Dataset of patient level data. Must be a tibble with columns named:
+#' - ptid: patient identifier
+#' - pfs.durn: duration of PFS from baseline
+#' - pfs.flag: event flag for PFS (=1 if progression or death occurred, 0 for censoring)
+#' - os.durn: duration of OS from baseline
+#' - os.flag: event flag for OS (=1 if death occurred, 0 for censoring)
+#' - ttp.durn: duration of TTP from baseline (usually should be equal to pfs.durn)
+#' - ttp.flag: event flag for TTP (=1 if progression occurred, 0 for censoring).
+#'
+#' Survival data for all other endpoints (time to progression, pre-progression death, post-progression survival) are derived from PFS and OS.
 #' @param dpam List of survival regressions for model endpoints. These must include time to progression (TTP) and pre-progression death (PPD).
+#' @param psmtype Either "simple" or "complex" PSM formulation
 #' @param Ty Time duration over which to calculate (defaults to 10 years). Assumes input is in years, and patient-level data is recorded in weeks.
 #' @param discrate Discount rate (%) per year (defaults to zero).
 #' @param lifetable Optional. The lifetable must be a dataframe with columns named time and lx. The first entry of the time column must be zero. Data should be sorted in ascending order by time, and all times must be unique.
@@ -48,29 +59,33 @@
 #'   pps_cf = find_bestfit_spl(fits$pps_cf, "aic")$fit,
 #'   pps_cr = find_bestfit_spl(fits$pps_cr, "aic")$fit
 #' )
-#' drmd_psm(dpam=params)
+#' drmd_psm(ptdata=bosonc, dpam=params)
 #' # Add a lifetable constraint
 #' ltable <- tibble::tibble(lttime=0:20, lx=1-lttime*0.05)
-#' drmd_psm(dpam=params, lifetable=ltable)
-drmd_psm <- function(dpam, Ty=10, discrate=0, lifetable=NA, timestep=1) {
+#' drmd_psm(ptdata=bosonc, dpam=params, lifetable=ltable)
+drmd_psm <- function(ptdata, dpam, psmtype="simple", Ty=10, discrate=0, lifetable=NA, timestep=1) {
   # Declare local variables
   Tw <- tvec <- pfprob <- osprob <- adjosprob <- adjfac <- adjprob <- vn <- NULL
   # Time horizon in weeks (ceiling)
   Tw <- convert_yrs2wks(Ty)
   # Create time vector, with half-cycle addition
-  tvec <- timestep*(1:floor(Tw/timestep)) + timestep/2
-  # Membership probabilities with lifetable constraint
+  tvec <- timestep*(0:floor(Tw/timestep)) + timestep/2
+  # Obtain all the hazards
+  allh <- calc_haz_psm(timevar=tvec, ptdata=ptdata, dpam=dpam, psmtype=psmtype)$adj
+  # PFS and OS probabilties from PSM
   pfprob <- prob_pf_psm(tvec, dpam)
   osprob <- prob_os_psm(tvec, dpam)
-  adjosprob <- constrain_survprob(osprob, lifetable=lifetable, timevec=tvec)
-  adjfac <- adjosprob/osprob
-  adjfac[is.na(adjfac)] <- 1
-  adjpfprob <- pfprob * adjfac
+  # OS and PFS may be constrained already by definitions of PPD and PPS
+  maxos <- exp(-cumsum(allh$os))
+  maxpfs <- exp(-cumsum(allh$pfs))
+  # Further constrain OS by lifetable
+  conos <- constrain_survprob(survprob1=maxos, survprob2=osprob, lifetable=lifetable, timevec=tvec)
+  conpfs <- constrain_survprob(survprob1=maxpfs, survprob2=pmin(pfprob, osprob), lifetable=lifetable, timevec=tvec)
   # Discount factor
   vn <- (1+discrate)^(-convert_wks2yrs(tvec+timestep/2))
   # Calculate RMDs
-  pf <- sum(adjpfprob*vn) * timestep
-  os <- sum(adjosprob*vn) * timestep
+  pf <- sum(conpfs*vn) * timestep
+  os <- sum(conos*vn) * timestep
   # Return values
   return(list(pf=pf, pd=os-pf, os=os))
 }
@@ -99,30 +114,37 @@ drmd_psm <- function(dpam, Ty=10, discrate=0, lifetable=NA, timestep=1) {
 #' 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 <- sppd <- sttp <- sos <- NULL
+  Tw <- tvec <- ppd.ts <- ttp.ts <- pps.ts <- NULsppd <- 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*(1:floor(Tw/timestep)) + timestep/2
+  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)
-  # Calculate S_PPD, S_TTP and S_OS
+  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))
-  # Next line is the difference with STM-CR
-  sos <- prob_os_stm_cf(tvec, dpam)
-  # Apply constraints to S_PPD and S_OS
-  adjsppd <- constrain_survprob(sppd, lifetable=lifetable, timevec=tvec)
-  adjos <- constrain_survprob(sos, lifetable=lifetable, timevec=tvec)
+  spps <- tvec |> purrr::map_dbl(~calc_surv(.x, pps.ts$type, pps.ts$spec))
+  # 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)
+  # Partial prob PD
+  con.partprobpd <- sttp*con.sppd*http/con.spps
+  con.partprobpd[con.spps==0] <- 0
+  con.probpd <- con.spps * cumsum(con.partprobpd)
   # Discount factor
   vn <- (1+discrate)^(-convert_wks2yrs(tvec+timestep/2))
   # Calculate RMDs
-  pf <- sum(sttp*adjsppd*vn) * timestep
-  os <- sum(adjos*vn) * timestep
+  pf <- sum(sttp*con.sppd*vn) * timestep
+  pd <- sum(con.probpd*vn) * timestep
   # Return values
-  return(list(pf=pf, pd=os-pf, os=os))
+  return(list(pf=pf, pd=pd, os=pf+pd))
 }
 
 #' Discretized Restricted Mean Duration calculation for State Transition Model Clock Reset structure
@@ -154,25 +176,44 @@ drmd_stm_cr <- function(dpam, Ty=10, discrate=0, lifetable=NA, timestep=1) {
   # Time horizon in weeks (ceiling)
   Tw <- convert_yrs2wks(Ty)
   # Create time vector, with half-cycle addition
-  tvec <- timestep*(1:floor(Tw/timestep)) + timestep/2
+  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)
-  # Calculate S_PPD, S_TTP and S_OS
+  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))
-  # Next line is the difference with STM-CF
-  sos <- prob_os_stm_cr(tvec, dpam)
-  # Apply constraints to S_PPD and S_OS
-  adjsppd <- constrain_survprob(sppd, lifetable=lifetable, timevec=tvec)
-  adjos <- constrain_survprob(sos, lifetable=lifetable, timevec=tvec)
+  # 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.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)
+    i.c.spps <- pmin(i.u.spps, i.slxt/i.slxu)
+    i.c.spps[i.slxu==0] <- 0
+    i.c.sppd * i.sttp * i.http * i.c.spps
+  }
+  integrand <- Vectorize(integrand, "u")
+  # PD membership probability is the integral
+  probpd <- function(t) {
+    stats::integrate(integrand, lower=0, upper=t, t=t)$value
+  }
+  probpd <- Vectorize(probpd, "t")
+  # Calculate the PD membership probability for each time
+  c.probpd <- probpd(tvec)
   # Discount factor
   vn <- (1+discrate)^(-convert_wks2yrs(tvec+timestep/2))
   # Calculate RMDs
-  pf <- sum(sttp*adjsppd*vn) * timestep
-  os <- sum(adjos*vn) * timestep
+  pf <- sum(sttp*c.sppd*vn) * timestep
+  pd <- sum(c.probpd*vn) * timestep
   # Return values
-  return(list(pf=pf, pd=os-pf, os=os))
+  return(list(pf=pf, pd=pd, os=pf+pd))
 }
 
 

R/lhoods.R

@@ -134,17 +134,17 @@ calc_likes_psm_simple <- function(ptdata, dpam, cuttime=0) {
       hppdu = calc_haz_psm(timevar=.data$pfs.durn,
                            ptdata=ptdata,
                            dpam=dpam,
-                           type="simple")$adj$ppd,
+                           psmtype="simple")$adj$ppd,
       httpu = pmax(0, .data$hpfsu-.data$hppdu),
       sppstu = calc_surv_psmpps(totime=.data$os.durn,
                                 fromtime=.data$pfs.durn,
                                 ptdata=ptdata,
                                 dpam=dpam,
-                                type="simple"),
+                                psmtype="simple"),
       hppst = calc_haz_psm(timevar=.data$os.durn,
                            ptdata=ptdata,
                            dpam=dpam,
-                           type="simple")$adj$pps
+                           psmtype="simple")$adj$pps
     )
   # Replace NA for zero hppst
   likedata$hppst[likedata$hppst==0] <- NA
@@ -240,17 +240,17 @@ calc_likes_psm_complex <- function(ptdata, dpam, cuttime=0) {
       hppdu = calc_haz_psm(timevar=.data$pfs.durn,
                            ptdata=ptdata,
                            dpam=dpam,
-                           type="complex")$adj$ppd,
+                           psmtype="complex")$adj$ppd,
       httpu = pmax(0, .data$hpfsu-.data$hppdu),
       sppstu = calc_surv_psmpps(totime=.data$os.durn,
                                 fromtime=.data$pfs.durn,
                                 ptdata=ptdata,
                                 dpam=dpam,
-                                type="complex"),
+                                psmtype="complex"),
       hppst = calc_haz_psm(timevar=.data$os.durn,
                            ptdata=ptdata,
                            dpam=dpam,
-                           type="complex")$adj$pps
+                           psmtype="complex")$adj$pps
     )
   likedata$hppst[likedata$hppst==0] <- NA
   likedata <- likedata |>

R/ltablesurv.R

@@ -97,8 +97,8 @@ vlookup <- function(indexval, indexvec, valvec, method="geom") {
 #' ltable <- tibble::tibble(lttime=0:10, lx=10-(0:10))
 #' calc_ltsurv(c(2, 2.5, 9.3), ltable)
 calc_ltsurv <- function(looktime, lifetable=NA){
-  if (!is.data.frame(lifetable)) stop("Lifetable must be specified")
-  if (lifetable$lttime[1]!=0) stop("Lifetable must run from time zero")
+  if (!is.data.frame(lifetable)) {return(rep(1, length(looktime)))}
+  if (lifetable$lttime[1]!=0) {stop("Lifetable must run from time zero")}
   vlookup(looktime, lifetable$lttime, lifetable$lx) / lifetable$lx[1]
 }
 
@@ -167,8 +167,9 @@ calc_ex <- function(Ty=10, lifetable, discrate=0) {
 
 #' Constrain survival probabilities according to hazards in a lifetable
 #' Recalculated constrained survival probabilities (by week) as the lower of the original unadjusted survival probability and the survival implied by the given lifetable (assumed indexed as years).
-#' @param survprob (Unconstrained) survival probability value or vector
-#' @param lifetable Lifetable
+#' @param survprob1 (Unconstrained) survival probability value or vector
+#' @param survprob2 Optional survival probability value or vector to constrain on (default = NA)
+#' @param lifetable Lifetable (default = NA)
 #' @param timevec Vector of times corresponding with survival probabilities above
 #' @return Vector of constrained survival probabilities
 #' @export
@@ -178,19 +179,32 @@ calc_ex <- function(Ty=10, lifetable, discrate=0) {
 #' constrain_survprob(survprob, lifetable=ltable)
 #' timevec <- 100*(0:4)
 #' constrain_survprob(survprob, lifetable=ltable, timevec=timevec)
-constrain_survprob <- function(survprob, lifetable, timevec=0:(length(survprob)-1)) {
-  # Check lifetable exists or return survprob
-  if (!is.data.frame(lifetable)) {return(survprob)}
-  # Vector of lifetables
-  lxprob <- calc_ltsurv(convert_wks2yrs(timevec), lifetable)
-  # Length of survprob
-  N <- length(survprob)
+#' survprob2 <- c(1,0.45,0.35,0.15,0)
+#' constrain_survprob(survprob, survprob2)
+constrain_survprob <- function(survprob1, survprob2=NA, lifetable=NA, timevec=0:(length(survprob1)-1)) {
+  # Check what exists
+  s2exists <- !is.na(survprob2)[1]
+  ltexists <- !is.na(lifetable)[1]
+  # Survprob1 and survprob2 must have equal length (when survprob2 is defined)
+  if(s2exists) {
+    stopifnot("Survival probability vectors must have equal length" = length(survprob1)==length(survprob1))
+    }
+  # Vector of lifetables (or ones, if lifetable not specified)
+  tN <- length(timevec)
+  lxprob <- rep(1, tN)
+  if(ltexists) {
+    lxprob <- calc_ltsurv(convert_wks2yrs(timevec), lifetable)
+  }
   # Cycle through each element
-  adjsurv <- slx <- sprob <- rep(NA, N)
-  adjsurv[1] <- survprob[1]
-  for (i in 2:N) {
+  adjsurv <- slx <- sprob <- rep(NA, tN)
+  adjsurv[1] <- survprob1[1]
+  for (i in 2:tN) {
     slx[i] <- ifelse(lxprob[i-1]==0, 1, lxprob[i]/lxprob[i-1])
-    sprob[i] <- ifelse(survprob[i-1]==0, 1, survprob[i]/survprob[i-1])
+    sprob[i] <- ifelse(survprob1[i-1]==0, 1, survprob1[i]/survprob1[i-1])
+    sprob[i] <- ifelse(s2exists,
+                       ifelse(survprob2[i-1]==0, sprob[i],
+                                 pmin(sprob[i], survprob2[i]/survprob2[i-1])),
+                        sprob[i])
     adjsurv[i] <- adjsurv[i-1] * pmin(slx[i], sprob[i])
   }
   return(adjsurv)