Skip to content

Commit 110535f

Browse files
dom-mustondom-muston
authored
Merge pull request #19 from Merck/Nextrel
Revise application of liftable constraints in STM and update vignette
2 parents b9d7ac2 + 97b5ba6 commit 110535f

15 files changed

+296
-173
lines changed

‎DESCRIPTION

+1-1
Original file line numberDiff line numberDiff line change
@@ -46,7 +46,6 @@ Collate:
4646
'basics.R'
4747
'brier.R'
4848
'datasets.R'
49-
'discrmd.R'
5049
'fitting-spl.R'
5150
'fitting.R'
5251
'lhoods.R'
@@ -55,3 +54,4 @@ Collate:
5554
'probgraphs.R'
5655
'psm3mkv-package.R'
5756
'resmeans.R'
57+
'discrmd.R'

‎R/discrmd.R

+73-32
Original file line numberDiff line numberDiff line change
@@ -24,7 +24,18 @@
2424

2525
#' Discretized Restricted Mean Duration calculation for Partitioned Survival Model
2626
#' Calculate restricted mean duration (RMD) in PF, PD and OS states under a Partitioned Survival Model structure.
27+
#' @param ptdata Dataset of patient level data. Must be a tibble with columns named:
28+
#' - ptid: patient identifier
29+
#' - pfs.durn: duration of PFS from baseline
30+
#' - pfs.flag: event flag for PFS (=1 if progression or death occurred, 0 for censoring)
31+
#' - os.durn: duration of OS from baseline
32+
#' - os.flag: event flag for OS (=1 if death occurred, 0 for censoring)
33+
#' - ttp.durn: duration of TTP from baseline (usually should be equal to pfs.durn)
34+
#' - ttp.flag: event flag for TTP (=1 if progression occurred, 0 for censoring).
35+
#'
36+
#' Survival data for all other endpoints (time to progression, pre-progression death, post-progression survival) are derived from PFS and OS.
2737
#' @param dpam List of survival regressions for model endpoints. These must include time to progression (TTP) and pre-progression death (PPD).
38+
#' @param psmtype Either "simple" or "complex" PSM formulation
2839
#' @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.
2940
#' @param discrate Discount rate (%) per year (defaults to zero).
3041
#' @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 @@
4859
#' pps_cf = find_bestfit_spl(fits$pps_cf, "aic")$fit,
4960
#' pps_cr = find_bestfit_spl(fits$pps_cr, "aic")$fit
5061
#' )
51-
#' drmd_psm(dpam=params)
62+
#' drmd_psm(ptdata=bosonc, dpam=params)
5263
#' # Add a lifetable constraint
5364
#' ltable <- tibble::tibble(lttime=0:20, lx=1-lttime*0.05)
54-
#' drmd_psm(dpam=params, lifetable=ltable)
55-
drmd_psm <- function(dpam, Ty=10, discrate=0, lifetable=NA, timestep=1) {
65+
#' drmd_psm(ptdata=bosonc, dpam=params, lifetable=ltable)
66+
drmd_psm <- function(ptdata, dpam, psmtype="simple", Ty=10, discrate=0, lifetable=NA, timestep=1) {
5667
# Declare local variables
5768
Tw <- tvec <- pfprob <- osprob <- adjosprob <- adjfac <- adjprob <- vn <- NULL
5869
# Time horizon in weeks (ceiling)
5970
Tw <- convert_yrs2wks(Ty)
6071
# Create time vector, with half-cycle addition
61-
tvec <- timestep*(1:floor(Tw/timestep)) + timestep/2
62-
# Membership probabilities with lifetable constraint
72+
tvec <- timestep*(0:floor(Tw/timestep)) + timestep/2
73+
# Obtain all the hazards
74+
allh <- calc_haz_psm(timevar=tvec, ptdata=ptdata, dpam=dpam, psmtype=psmtype)$adj
75+
# PFS and OS probabilties from PSM
6376
pfprob <- prob_pf_psm(tvec, dpam)
6477
osprob <- prob_os_psm(tvec, dpam)
65-
adjosprob <- constrain_survprob(osprob, lifetable=lifetable, timevec=tvec)
66-
adjfac <- adjosprob/osprob
67-
adjfac[is.na(adjfac)] <- 1
68-
adjpfprob <- pfprob * adjfac
78+
# OS and PFS may be constrained already by definitions of PPD and PPS
79+
maxos <- exp(-cumsum(allh$os))
80+
maxpfs <- exp(-cumsum(allh$pfs))
81+
# Further constrain OS by lifetable
82+
conos <- constrain_survprob(survprob1=maxos, survprob2=osprob, lifetable=lifetable, timevec=tvec)
83+
conpfs <- constrain_survprob(survprob1=maxpfs, survprob2=pmin(pfprob, osprob), lifetable=lifetable, timevec=tvec)
6984
# Discount factor
7085
vn <- (1+discrate)^(-convert_wks2yrs(tvec+timestep/2))
7186
# Calculate RMDs
72-
pf <- sum(adjpfprob*vn) * timestep
73-
os <- sum(adjosprob*vn) * timestep
87+
pf <- sum(conpfs*vn) * timestep
88+
os <- sum(conos*vn) * timestep
7489
# Return values
7590
return(list(pf=pf, pd=os-pf, os=os))
7691
}
@@ -99,30 +114,37 @@ drmd_psm <- function(dpam, Ty=10, discrate=0, lifetable=NA, timestep=1) {
99114
#' drmd_stm_cf(dpam=params, lifetable=ltable)
100115
drmd_stm_cf <- function(dpam, Ty=10, discrate=0, lifetable=NA, timestep=1) {
101116
# Declare local variables
102-
Tw <- tvec <- ppd.ts <- ttp.ts <- sppd <- sttp <- sos <- NULL
117+
Tw <- tvec <- ppd.ts <- ttp.ts <- pps.ts <- NULsppd <- sttp <- sos <- NULL
103118
adjsppd <- adjos <- vn <- pf <- os <- NULL
104119
# Time horizon in weeks (ceiling)
105120
Tw <- convert_yrs2wks(Ty)
106121
# Create time vector, with half-cycle addition
107-
tvec <- timestep*(1:floor(Tw/timestep)) + timestep/2
122+
tvec <- timestep*(0:floor(Tw/timestep)) + timestep/2
108123
# Pull out type and spec for PPD and TTP
109124
ppd.ts <- convert_fit2spec(dpam$ppd)
110125
ttp.ts <- convert_fit2spec(dpam$ttp)
111-
# Calculate S_PPD, S_TTP and S_OS
126+
pps.ts <- convert_fit2spec(dpam$pps_cf)
127+
# Obtain hazard and survival functions
128+
hppd <- tvec |> purrr::map_dbl(~calc_haz(.x, ppd.ts$type, ppd.ts$spec))
129+
http <- tvec |> purrr::map_dbl(~calc_haz(.x, ttp.ts$type, ttp.ts$spec))
130+
hpps <- tvec |> purrr::map_dbl(~calc_haz(.x, pps.ts$type, pps.ts$spec))
112131
sppd <- tvec |> purrr::map_dbl(~calc_surv(.x, ppd.ts$type, ppd.ts$spec))
113132
sttp <- tvec |> purrr::map_dbl(~calc_surv(.x, ttp.ts$type, ttp.ts$spec))
114-
# Next line is the difference with STM-CR
115-
sos <- prob_os_stm_cf(tvec, dpam)
116-
# Apply constraints to S_PPD and S_OS
117-
adjsppd <- constrain_survprob(sppd, lifetable=lifetable, timevec=tvec)
118-
adjos <- constrain_survprob(sos, lifetable=lifetable, timevec=tvec)
133+
spps <- tvec |> purrr::map_dbl(~calc_surv(.x, pps.ts$type, pps.ts$spec))
134+
# Derive the constrained S_PPD and S_PFS
135+
con.sppd <- constrain_survprob(sppd, lifetable=lifetable, timevec=tvec)
136+
con.spps <- constrain_survprob(spps, lifetable=lifetable, timevec=tvec)
137+
# Partial prob PD
138+
con.partprobpd <- sttp*con.sppd*http/con.spps
139+
con.partprobpd[con.spps==0] <- 0
140+
con.probpd <- con.spps * cumsum(con.partprobpd)
119141
# Discount factor
120142
vn <- (1+discrate)^(-convert_wks2yrs(tvec+timestep/2))
121143
# Calculate RMDs
122-
pf <- sum(sttp*adjsppd*vn) * timestep
123-
os <- sum(adjos*vn) * timestep
144+
pf <- sum(sttp*con.sppd*vn) * timestep
145+
pd <- sum(con.probpd*vn) * timestep
124146
# Return values
125-
return(list(pf=pf, pd=os-pf, os=os))
147+
return(list(pf=pf, pd=pd, os=pf+pd))
126148
}
127149

128150
#' 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) {
154176
# Time horizon in weeks (ceiling)
155177
Tw <- convert_yrs2wks(Ty)
156178
# Create time vector, with half-cycle addition
157-
tvec <- timestep*(1:floor(Tw/timestep)) + timestep/2
179+
tvec <- timestep*(0:floor(Tw/timestep)) + timestep/2
158180
# Pull out type and spec for PPD and TTP
159181
ppd.ts <- convert_fit2spec(dpam$ppd)
160182
ttp.ts <- convert_fit2spec(dpam$ttp)
161-
# Calculate S_PPD, S_TTP and S_OS
183+
pps.ts <- convert_fit2spec(dpam$pps_cr)
184+
# Obtain unconstrained survival functions
162185
sppd <- tvec |> purrr::map_dbl(~calc_surv(.x, ppd.ts$type, ppd.ts$spec))
163186
sttp <- tvec |> purrr::map_dbl(~calc_surv(.x, ttp.ts$type, ttp.ts$spec))
164-
# Next line is the difference with STM-CF
165-
sos <- prob_os_stm_cr(tvec, dpam)
166-
# Apply constraints to S_PPD and S_OS
167-
adjsppd <- constrain_survprob(sppd, lifetable=lifetable, timevec=tvec)
168-
adjos <- constrain_survprob(sos, lifetable=lifetable, timevec=tvec)
187+
# Derive the constrained S_PPD
188+
c.sppd <- constrain_survprob(sppd, lifetable=lifetable, timevec=tvec)
189+
# Integrand with constraints on S_PPD and S_PPS
190+
integrand <- function(u, t) {
191+
i.http <- calc_haz(u, ttp.ts$type, ttp.ts$spec)
192+
i.sttp <- calc_surv(u, ttp.ts$type, ttp.ts$spec)
193+
i.u.sppd <- calc_surv(u, ppd.ts$type, ppd.ts$spec)
194+
i.u.spps <- calc_surv(t-u, pps.ts$type, pps.ts$spec)
195+
i.slxu <- calc_ltsurv(convert_wks2yrs(u), lifetable=lifetable)
196+
i.slxt <- calc_ltsurv(convert_wks2yrs(t), lifetable=lifetable)
197+
i.c.sppd <- pmin(i.u.sppd, i.slxu)
198+
i.c.spps <- pmin(i.u.spps, i.slxt/i.slxu)
199+
i.c.spps[i.slxu==0] <- 0
200+
i.c.sppd * i.sttp * i.http * i.c.spps
201+
}
202+
integrand <- Vectorize(integrand, "u")
203+
# PD membership probability is the integral
204+
probpd <- function(t) {
205+
stats::integrate(integrand, lower=0, upper=t, t=t)$value
206+
}
207+
probpd <- Vectorize(probpd, "t")
208+
# Calculate the PD membership probability for each time
209+
c.probpd <- probpd(tvec)
169210
# Discount factor
170211
vn <- (1+discrate)^(-convert_wks2yrs(tvec+timestep/2))
171212
# Calculate RMDs
172-
pf <- sum(sttp*adjsppd*vn) * timestep
173-
os <- sum(adjos*vn) * timestep
213+
pf <- sum(sttp*c.sppd*vn) * timestep
214+
pd <- sum(c.probpd*vn) * timestep
174215
# Return values
175-
return(list(pf=pf, pd=os-pf, os=os))
216+
return(list(pf=pf, pd=pd, os=pf+pd))
176217
}
177218

178219

‎R/lhoods.R

+6-6
Original file line numberDiff line numberDiff line change
@@ -134,17 +134,17 @@ calc_likes_psm_simple <- function(ptdata, dpam, cuttime=0) {
134134
hppdu = calc_haz_psm(timevar=.data$pfs.durn,
135135
ptdata=ptdata,
136136
dpam=dpam,
137-
type="simple")$adj$ppd,
137+
psmtype="simple")$adj$ppd,
138138
httpu = pmax(0, .data$hpfsu-.data$hppdu),
139139
sppstu = calc_surv_psmpps(totime=.data$os.durn,
140140
fromtime=.data$pfs.durn,
141141
ptdata=ptdata,
142142
dpam=dpam,
143-
type="simple"),
143+
psmtype="simple"),
144144
hppst = calc_haz_psm(timevar=.data$os.durn,
145145
ptdata=ptdata,
146146
dpam=dpam,
147-
type="simple")$adj$pps
147+
psmtype="simple")$adj$pps
148148
)
149149
# Replace NA for zero hppst
150150
likedata$hppst[likedata$hppst==0] <- NA
@@ -240,17 +240,17 @@ calc_likes_psm_complex <- function(ptdata, dpam, cuttime=0) {
240240
hppdu = calc_haz_psm(timevar=.data$pfs.durn,
241241
ptdata=ptdata,
242242
dpam=dpam,
243-
type="complex")$adj$ppd,
243+
psmtype="complex")$adj$ppd,
244244
httpu = pmax(0, .data$hpfsu-.data$hppdu),
245245
sppstu = calc_surv_psmpps(totime=.data$os.durn,
246246
fromtime=.data$pfs.durn,
247247
ptdata=ptdata,
248248
dpam=dpam,
249-
type="complex"),
249+
psmtype="complex"),
250250
hppst = calc_haz_psm(timevar=.data$os.durn,
251251
ptdata=ptdata,
252252
dpam=dpam,
253-
type="complex")$adj$pps
253+
psmtype="complex")$adj$pps
254254
)
255255
likedata$hppst[likedata$hppst==0] <- NA
256256
likedata <- likedata |>

‎R/ltablesurv.R

+29-15
Original file line numberDiff line numberDiff line change
@@ -97,8 +97,8 @@ vlookup <- function(indexval, indexvec, valvec, method="geom") {
9797
#' ltable <- tibble::tibble(lttime=0:10, lx=10-(0:10))
9898
#' calc_ltsurv(c(2, 2.5, 9.3), ltable)
9999
calc_ltsurv <- function(looktime, lifetable=NA){
100-
if (!is.data.frame(lifetable)) stop("Lifetable must be specified")
101-
if (lifetable$lttime[1]!=0) stop("Lifetable must run from time zero")
100+
if (!is.data.frame(lifetable)) {return(rep(1, length(looktime)))}
101+
if (lifetable$lttime[1]!=0) {stop("Lifetable must run from time zero")}
102102
vlookup(looktime, lifetable$lttime, lifetable$lx) / lifetable$lx[1]
103103
}
104104

@@ -167,8 +167,9 @@ calc_ex <- function(Ty=10, lifetable, discrate=0) {
167167

168168
#' Constrain survival probabilities according to hazards in a lifetable
169169
#' 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).
170-
#' @param survprob (Unconstrained) survival probability value or vector
171-
#' @param lifetable Lifetable
170+
#' @param survprob1 (Unconstrained) survival probability value or vector
171+
#' @param survprob2 Optional survival probability value or vector to constrain on (default = NA)
172+
#' @param lifetable Lifetable (default = NA)
172173
#' @param timevec Vector of times corresponding with survival probabilities above
173174
#' @return Vector of constrained survival probabilities
174175
#' @export
@@ -178,19 +179,32 @@ calc_ex <- function(Ty=10, lifetable, discrate=0) {
178179
#' constrain_survprob(survprob, lifetable=ltable)
179180
#' timevec <- 100*(0:4)
180181
#' constrain_survprob(survprob, lifetable=ltable, timevec=timevec)
181-
constrain_survprob <- function(survprob, lifetable, timevec=0:(length(survprob)-1)) {
182-
# Check lifetable exists or return survprob
183-
if (!is.data.frame(lifetable)) {return(survprob)}
184-
# Vector of lifetables
185-
lxprob <- calc_ltsurv(convert_wks2yrs(timevec), lifetable)
186-
# Length of survprob
187-
N <- length(survprob)
182+
#' survprob2 <- c(1,0.45,0.35,0.15,0)
183+
#' constrain_survprob(survprob, survprob2)
184+
constrain_survprob <- function(survprob1, survprob2=NA, lifetable=NA, timevec=0:(length(survprob1)-1)) {
185+
# Check what exists
186+
s2exists <- !is.na(survprob2)[1]
187+
ltexists <- !is.na(lifetable)[1]
188+
# Survprob1 and survprob2 must have equal length (when survprob2 is defined)
189+
if(s2exists) {
190+
stopifnot("Survival probability vectors must have equal length" = length(survprob1)==length(survprob1))
191+
}
192+
# Vector of lifetables (or ones, if lifetable not specified)
193+
tN <- length(timevec)
194+
lxprob <- rep(1, tN)
195+
if(ltexists) {
196+
lxprob <- calc_ltsurv(convert_wks2yrs(timevec), lifetable)
197+
}
188198
# Cycle through each element
189-
adjsurv <- slx <- sprob <- rep(NA, N)
190-
adjsurv[1] <- survprob[1]
191-
for (i in 2:N) {
199+
adjsurv <- slx <- sprob <- rep(NA, tN)
200+
adjsurv[1] <- survprob1[1]
201+
for (i in 2:tN) {
192202
slx[i] <- ifelse(lxprob[i-1]==0, 1, lxprob[i]/lxprob[i-1])
193-
sprob[i] <- ifelse(survprob[i-1]==0, 1, survprob[i]/survprob[i-1])
203+
sprob[i] <- ifelse(survprob1[i-1]==0, 1, survprob1[i]/survprob1[i-1])
204+
sprob[i] <- ifelse(s2exists,
205+
ifelse(survprob2[i-1]==0, sprob[i],
206+
pmin(sprob[i], survprob2[i]/survprob2[i-1])),
207+
sprob[i])
194208
adjsurv[i] <- adjsurv[i-1] * pmin(slx[i], sprob[i])
195209
}
196210
return(adjsurv)

0 commit comments

Comments
 (0)