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