# IPCW Cumulative Cost

#### Klaus Holst & Thomas Scheike

#### 2026-01-11

# Overview

We here describe how to do 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 the cost is realized and accumulated, and \(Z(t)\) is the cost (or marks) at the event times. The counting process can be a mix of random and fixed times. The data would thus be represented in counting process format with the marks/costs going along with the event times. There are many additional uses of such cumulative process, for example, when considering time-lost in a recurrent events setting, that 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 in a setting with independent right-censoring given \(X\), and based on 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 situation with a terminal this is based on a derived competing risks data that keeps track of the competing terminal event.

Regression modelling of this quantity is also possible using competing risks regression models, using for example, the cifreg function of mets.

# HF-action data

Considering the HF-action data we simulate a severity score for each event.

    library(mets)
    data(hfactioncpx12)
    hf <- hfactioncpx12
    hf$severity <- abs((5+rnorm(741)*2))[hf$id]

    proc_design <- mets:::proc_design
    ## marginal mean using formula
    outNZ <- recurrentMarginal(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.47196 0.6334775 9.301149 11.79016      0
    #>
    #> [[2]]
    #>     new.time     mean        se  CI-2.5% CI-97.5% strata
    #> 601        3 9.692896 0.7092806 8.397821 11.18769      1

    outN <- recurrentMarginal(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      1

For 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 however based on different formula 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.47196  0.63344  9.23045 11.71348       0
    #> treatment1  9.69290  0.70923  8.30283 11.08297       0
    head(iid(outNZ3))
    #>            [,1]      [,2]
    #> 1  0.0034815860 0.0000000
    #> 2  0.0068980137 0.0000000
    #> 3  0.0000000000 0.1248012
    #> 4 -0.0134963606 0.0000000
    #> 5 -0.0005487192 0.0000000
    #> 6 -0.0350784674 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.000000000

We 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
    #> coeffients:
    #>             Estimate      S.E.   dU^-1/2 P-value
    #> treatment1 -0.099914  0.090172  0.024188  0.2678
    #>
    #> exp(coeffients):
    #>            Estimate    2.5%  97.5%
    #> treatment1  0.90492 0.75832 1.0798
    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
    #> coeffients:
    #>             Estimate      S.E.   dU^-1/2 P-value
    #> treatment1 -0.110404  0.078656  0.053776  0.1604
    #>
    #> exp(coeffients):
    #>            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 number of expected events.

# Exceed threshold

Finally, we also estimate the probability of exceeding a cumulative severity at 1,5,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.06295909
    #> exceed>=1      3 0.93704091
    #> exceed>=5      3 0.90309925
    #> exceed>=10     3 0.76036644
    #> exceed>=20     3 0.44547988
    #>
    #> $`0`$se
    #>            times
    #>                3 2.99865085
    #> N<1            3 0.02704848
    #> exceed>=1      3 0.02704848
    #> exceed>=5      3 0.03045582
    #> exceed>=10     3 0.04221875
    #> exceed>=20     3 0.05065399
    #>
    #> $`0`$lower
    #>      times
    #> [1,]     3 2.9986509
    #> [2,]     3 0.1145014
    #> [3,]     3 0.8854986
    #> [4,]     3 0.8453369
    #> [5,]     3 0.6819627
    #> [6,]     3 0.3564847
    #>
    #> $`0`$upper
    #>      times
    #> [1,]     3 2.998650853
    #> [2,]     3 0.008416702
    #> [3,]     3 0.991583298
    #> [4,]     3 0.964808502
    #> [5,]     3 0.847784047
    #> [6,]     3 0.556692410
    #>
    #>
    #> $`1`
    #> $`1`$prob
    #>            times
    #>                3 2.99865085
    #> N<1            3 0.01781227
    #> exceed>=1      3 0.98218773
    #> exceed>=5      3 0.95475115
    #> exceed>=10     3 0.81657966
    #> exceed>=20     3 0.53639776
    #>
    #> $`1`$se
    #>            times
    #>                3 2.998650853
    #> N<1            3 0.009437445
    #> exceed>=1      3 0.009437445
    #> exceed>=5      3 0.014947541
    #> exceed>=10     3 0.035245302
    #> exceed>=20     3 0.057003157
    #>
    #> $`1`$lower
    #>      times
    #> [1,]     3 2.99865085
    #> [2,]     3 0.03613623
    #> [3,]     3 0.96386377
    #> [4,]     3 0.92589943
    #> [5,]     3 0.75034139
    #> [6,]     3 0.43554146
    #>
    #> $`1`$upper
    #>      times
    #> [1,]     3 2.9986509
    #> [2,]     3 0.0000000
    #> [3,]     3 1.0000000
    #> [4,]     3 0.9845019
    #> [5,]     3 0.8886653
    #> [6,]     3 0.6606089

# Cumulative time lost for recurrent events

The cumulative time lost for recurrent events has been defined as \[\begin{align*} {\cal M}(t) = E[ \int_0^\tau (\tau-s) dN(s) ] = \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.00000000

# SessionInfo

    sessionInfo()
    #> R version 4.5.2 (2025-10-31)
    #> Platform: aarch64-apple-darwin25.1.0
    #> Running under: macOS Tahoe 26.2
    #>
    #> Matrix products: default
    #> BLAS:   /Users/kkzh/.asdf/installs/r/4.5.2/lib/R/lib/libRblas.dylib
    #> LAPACK: /Users/kkzh/.asdf/installs/r/4.5.2/lib/R/lib/libRlapack.dylib;  LAPACK version 3.12.1
    #>
    #> locale:
    #> [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
    #>
    #> time zone: Europe/Copenhagen
    #> tzcode source: internal
    #>
    #> attached base packages:
    #> [1] stats     graphics  grDevices utils     datasets  methods   base
    #>
    #> other attached packages:
    #> [1] timereg_2.0.7  survival_3.8-3 mets_1.3.9
    #>
    #> loaded via a namespace (and not attached):
    #>  [1] cli_3.6.5              knitr_1.50             rlang_1.1.6
    #>  [4] xfun_0.54              jsonlite_2.0.0         future.apply_1.20.1
    #>  [7] listenv_0.10.0         lava_1.8.3             htmltools_0.5.8.1
    #> [10] sass_0.4.10            rmarkdown_2.30         grid_4.5.2
    #> [13] evaluate_1.0.5         jquerylib_0.1.4        fastmap_1.2.0
    #> [16] numDeriv_2016.8-1.1    yaml_2.3.10            mvtnorm_1.3-3
    #> [19] lifecycle_1.0.4        compiler_4.5.2         codetools_0.2-20
    #> [22] ucminf_1.2.2           Rcpp_1.1.1             future_1.68.0
    #> [25] lattice_0.22-7         digest_0.6.39          R6_2.6.1
    #> [28] parallelly_1.46.1      parallel_4.5.2         splines_4.5.2
    #> [31] Matrix_1.7-4           bslib_0.9.0            tools_4.5.2
    #> [34] RcppArmadillo_15.2.3-1 globals_0.18.0         cachem_1.1.0