IPCW Cumulative Cost
Overview
We describe how to perform regression modelling for cumulative cost \[\begin{align*} {\cal U}(t) & = \int_0^t Z(s) dN(s) \end{align*}\] where \(N(s)\) is a counting process that registers the times at which costs are realised and accumulated, and \(Z(t)\) is the cost (or mark) at the event times. The counting process can be a mix of random and fixed times, and the data are represented in counting process format with the marks/costs attached to the event times. There are many additional uses of such cumulative processes; for example, when considering time lost in a recurrent events setting, which we return to below.
We can estimate the marginal mean of the cumulative process \[\begin{align*} \nu(t) & = E ( {\cal U}(t) ) \end{align*}\] possibly for strata with standard errors based on the derived influence function.
We provide semi-parametric regression modelling using the proportional model \[\begin{align*} E ( {\cal U}(t) | X) & = \Lambda_0(t) \exp( X^T \beta). \end{align*}\]
In addition for a fixed time-point \(t \in [0,\tau]\) we can estimate the mean given covariates \[\begin{align*} E ( {\cal U}(t) | X) & = \exp( X^T \beta) \end{align*}\] where \(\tau\) is some maximum follow-up time.
- These quantities are estimated under independent right-censoring
given \(X\), using IPCW-adjusted
estimating equations,
- similarly to the Ghosh-Lin model for recurrent events.
- A terminal event can be specified.
We also estimate the probability of exceeding thresholds over time \[\begin{align*} P ( {\cal U}(t) > k ) & = \mu_k(t), \end{align*}\] and in the setting with a terminal event this is based on a derived competing risks data structure that keeps track of the competing terminal event.
Regression modelling of this quantity is also possible using
competing risks regression models, for example via the
cifreg function in mets.
HF-action data
Using the HF-action data, we simulate a severity score for each event.
library(mets)
data(hfactioncpx12)
hf <- hfactioncpx12
set.seed(1)
hf$severity <- abs((5+rnorm(741)*2))[hf$id]
## marginal mean using formula
outNZ <- recurrent_marginal(Event(entry,time,status)~strata(treatment)+cluster(id)
+marks(severity),hf,cause=1,death.code=2)
plot(outNZ,se=TRUE)
summary(outNZ,times=3)
#> [[1]]
#> new.time mean se CI-2.5% CI-97.5% strata
#> 682 3 10.3249 0.6044746 9.2056 11.5803 0
#>
#> [[2]]
#> new.time mean se CI-2.5% CI-97.5% strata
#> 601 3 9.33613 0.6429995 8.157232 10.68541 1
outN <- recurrent_marginal(Event(entry,time,status)~strata(treatment)+cluster(id),data=hf,
cause=1,death.code=2)
plot(outN,se=TRUE,add=TRUE)
summary(outN,times=3)
#> [[1]]
#> new.time mean se CI-2.5% CI-97.5% strata
#> 682 3 2.118496 0.1138572 1.906692 2.353829 0
#>
#> [[2]]
#> new.time mean se CI-2.5% CI-97.5% strata
#> 601 3 1.924062 0.1216577 1.699801 2.177912 1For comparison we also compute the IPCW estimates with and without marks at time 3 using the linear model, and note that they are identical. Standard errors are based on different formulae that are asymptotically equivalent, and we note that they are very similar.
outNZ3 <- recregIPCW(Event(entry,time,status)~-1+treatment+cluster(id)+marks(severity),data=hf,
cause=1,death.code=2,time=3,cens.model=~strata(treatment),model="lin")
summary(outNZ3)
#> n events
#> 741 1281
#>
#> 741 clusters
#> coeffients:
#> Estimate Std.Err 2.5% 97.5% P-value
#> treatment0 10.32490 0.60444 9.14022 11.50959 0
#> treatment1 9.33613 0.64295 8.07597 10.59629 0
head(iid(outNZ3))
#> [,1] [,2]
#> 1 -0.0013037837 0.0000000
#> 2 0.0089249278 0.0000000
#> 3 0.0000000000 0.0211858
#> 4 -0.0056613816 0.0000000
#> 5 -0.0005288663 0.0000000
#> 6 -0.0342554114 0.0000000
outN3 <- recregIPCW(Event(entry,time,status)~-1+treatment+cluster(id),data=hf,cause=1,death.code=2,time=3,
cens.model=~strata(treatment),model="lin")
summary(outN3)
#> n events
#> 741 1281
#>
#> 741 clusters
#> coeffients:
#> Estimate Std.Err 2.5% 97.5% P-value
#> treatment0 2.11850 0.11385 1.89535 2.34164 0
#> treatment1 1.92406 0.12165 1.68564 2.16248 0
head(iid(outN3))
#> [,1] [,2]
#> 1 0.0004542472 0.000000000
#> 2 0.0009756994 0.000000000
#> 3 0.0000000000 0.009301496
#> 4 -0.0029668336 0.000000000
#> 5 -0.0001120764 0.000000000
#> 6 -0.0070693971 0.000000000We also apply the semiparametric proportional cost model with IPCW adjustment:
propNZ <- recreg(Event(entry,time,status)~treatment+marks(severity)+cluster(id),data=hf,cause=1,death.code=2)
summary(propNZ)
#>
#> n events
#> 2132 1391
#>
#> 741 clusters
#> coefficients:
#> Estimate S.E. dU^-1/2 P-value
#> treatment1 -0.101014 0.089005 0.024304 0.2564
#>
#> exp(coefficients):
#> Estimate 2.5% 97.5%
#> treatment1 0.90392 0.75922 1.0762
plot(propNZ,main="Baselines")
GL <- recreg(Event(entry,time,status)~treatment+cluster(id),hf,cause=1,death.code=2)
summary(GL)
#>
#> n events
#> 2132 1391
#>
#> 741 clusters
#> coefficients:
#> Estimate S.E. dU^-1/2 P-value
#> treatment1 -0.110404 0.078656 0.053776 0.1604
#>
#> exp(coefficients):
#> Estimate 2.5% 97.5%
#> treatment1 0.89547 0.76754 1.0447
plot(GL,add=TRUE,col=2)
Those treated have 14% lower cumulative severity and 11% lower expected
number of events.
Exceed threshold
Finally, we estimate the probability of exceeding cumulative severity thresholds of 1, 5, and 10:
ooNZ <- prob_exceed_recurrent(Event(entry,time,status)~strata(treatment)+cluster(id)+marks(severity),data=hf,
cause=1,death.code=2,exceed=c(1,5,10,20))
plot(ooNZ,strata=1)
plot(ooNZ,strata=2,add=TRUE)
summary(ooNZ,times=3)
#> $`0`
#> $`0`$prob
#> times
#> 3 2.99865085
#> N<1 3 0.06474819
#> exceed>=1 3 0.93525181
#> exceed>=5 3 0.90118229
#> exceed>=10 3 0.74537955
#> exceed>=20 3 0.49282394
#>
#> $`0`$se
#> times
#> 3 2.99865085
#> N<1 3 0.02657120
#> exceed>=1 3 0.02657120
#> exceed>=5 3 0.03013468
#> exceed>=10 3 0.04466101
#> exceed>=20 3 0.05149605
#>
#> $`0`$lower
#> times
#> [1,] 3 2.9986509
#> [2,] 3 0.1154033
#> [3,] 3 0.8845967
#> [4,] 3 0.8440133
#> [5,] 3 0.6627899
#> [6,] 3 0.4015580
#>
#> $`0`$upper
#> times
#> [1,] 3 2.99865085
#> [2,] 3 0.01119234
#> [3,] 3 0.98880766
#> [4,] 3 0.96222364
#> [5,] 3 0.83826056
#> [6,] 3 0.60483284
#>
#>
#> $`1`
#> $`1`$prob
#> times
#> 3 2.99865085
#> N<1 3 0.01770732
#> exceed>=1 3 0.98229268
#> exceed>=5 3 0.94934557
#> exceed>=10 3 0.79560979
#> exceed>=20 3 0.53319063
#>
#> $`1`$se
#> times
#> 3 2.99865085
#> N<1 3 0.01097663
#> exceed>=1 3 0.01097663
#> exceed>=5 3 0.01645312
#> exceed>=10 3 0.03727391
#> exceed>=20 3 0.05545584
#>
#> $`1`$lower
#> times
#> [1,] 3 2.99865085
#> [2,] 3 0.03898725
#> [3,] 3 0.96101275
#> [4,] 3 0.91763959
#> [5,] 3 0.72580801
#> [6,] 3 0.43486168
#>
#> $`1`$upper
#> times
#> [1,] 3 2.9986509
#> [2,] 3 0.0000000
#> [3,] 3 1.0000000
#> [4,] 3 0.9821470
#> [5,] 3 0.8721245
#> [6,] 3 0.6537533Cumulative time lost for recurrent events
The cumulative time lost for recurrent events is defined as \[\begin{align*} {\cal M}(t) = E\left[ \int_0^\tau (\tau-s) dN(s) \right] = \int_0^\tau \mu(s) ds \end{align*}\] where \(\mu(t) = E( N(t) )\) is the marginal mean of the recurrent events at time \(t\).
hf$lost5 <- 5-hf$time
RecLost <- recregIPCW(Event(entry,time,status)~-1+treatment+cluster(id)+marks(lost5),data=hf,
cause=1,death.code=2,time=5,cens.model=~strata(treatment),model="lin")
summary(RecLost)
#> n events
#> 741 1391
#>
#> 741 clusters
#> coeffients:
#> Estimate Std.Err 2.5% 97.5% P-value
#> treatment0 8.58300 0.42951 7.74118 9.42482 0
#> treatment1 7.66234 0.46400 6.75292 8.57177 0
head(iid(RecLost))
#> [,1] [,2]
#> 1 0.0016920221 0.00000000
#> 2 0.0073388996 0.00000000
#> 3 0.0000000000 0.02120478
#> 4 -0.0095548150 0.00000000
#> 5 -0.0005696809 0.00000000
#> 6 -0.0201750011 0.00000000SessionInfo
sessionInfo()
#> R version 4.6.0 (2026-04-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 24.04.4 LTS
#>
#> Matrix products: default
#> BLAS: /home/kkzh/.asdf/installs/r/4.6.0/lib/R/lib/libRblas.so
#> LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0 LAPACK version 3.12.0
#>
#> locale:
#> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
#> [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
#> [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
#> [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
#> [9] LC_ADDRESS=C LC_TELEPHONE=C
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
#>
#> time zone: Europe/Copenhagen
#> tzcode source: system (glibc)
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] timereg_2.0.7 survival_3.8-6 mets_1.3.10
#>
#> loaded via a namespace (and not attached):
#> [1] cli_3.6.6 knitr_1.51 rlang_1.2.0
#> [4] xfun_0.57 otel_0.2.0 future.apply_1.20.2
#> [7] listenv_0.10.1 lava_1.9.1 stats4_4.6.0
#> [10] grid_4.6.0 evaluate_1.0.5 yaml_2.3.12
#> [13] mvtnorm_1.3-7 numDeriv_2016.8-1.1 compiler_4.6.0
#> [16] codetools_0.2-20 Rcpp_1.1.1-1.1 ucminf_1.2.3
#> [19] future_1.70.0 lattice_0.22-9 digest_0.6.39
#> [22] parallelly_1.47.0 parallel_4.6.0 splines_4.6.0
#> [25] Matrix_1.7-5 tools_4.6.0 RcppArmadillo_15.2.6-1
#> [28] globals_0.19.1