Simple and fast version for IPCW regression for just one time-point thus fitting the model
$$E( min(T, t) | X ) = exp( X^T beta) $$ or in the case of competing risks data
$$E( I(epsilon=1) (t - min(T ,t)) | X ) = exp( X^T beta) $$ thus given years lost to
cause, see binreg for the arguments.
resmeanIPCW(formula, data, outcome = c("rmst", "rmtl"), ...)When the status is binary assumes it is a survival setting and default is to consider outcome Y=min(T,t), if status has more than two levels, then computes years lost due to the specified cause, thus using the response $$ Y = (t-min(T,t)) I(status=cause) $$
Based on binomial regresion IPCW response estimating equation: $$ X ( \Delta(min(T,t)) Y /G_c(min(T,t)) - exp( X^T beta)) = 0 $$ for IPCW adjusted responses. Here $$ \Delta(min(T,t)) = I ( min(T ,t) \leq C ) $$ is indicator of being uncensored. Concretely, the uncensored observations at time t will count those with an event (of any type) before t and those with a censoring time at t or further out. One should therefore be a bit careful when data has been constructed such that some of the event times T are equivalent to t.
Can also solve the binomial regresion IPCW response estimating equation: $$ h(X) X ( \Delta(min(T,t)) Y /G_c(min(T,t)) - exp( X^T beta)) = 0 $$ for IPCW adjusted responses where $h$ is given as an argument together with iid of censoring with h.
By using appropriately the h argument we can also do the efficient IPCW estimator estimator.
Variance is based on $$ \sum w_i^2 $$ also with IPCW adjustment, and naive.var is variance under known censoring model.
When Ydirect is given it solves : $$ X ( \Delta(min(T,t)) Ydirect /G_c(min(T,t)) - exp( X^T beta)) = 0 $$ for IPCW adjusted responses.
The actual influence (type="II") function is based on augmenting with $$ X \int_0^t E(Y | T>s) /G_c(s) dM_c(s) $$ and alternatively just solved directly (type="I") without any additional terms.
Censoring model may depend on strata.
library(mets)
data(bmt); bmt$time <- bmt$time+runif(nrow(bmt))*0.001
# E( min(T;t) | X ) = exp( a+b X) with IPCW estimation
out <- resmeanIPCW(Event(time,cause!=0)~tcell+platelet+age,bmt,
time=50,cens.model=~strata(platelet),model="exp")
summary(out)
#> n events
#> 408 245
#>
#> 408 clusters
#> coeffients:
#> Estimate Std.Err 2.5% 97.5% P-value
#> (Intercept) 3.014319 0.065114 2.886698 3.141940 0.0000
#> tcell 0.106526 0.138268 -0.164474 0.377526 0.4410
#> platelet 0.247104 0.097338 0.056326 0.437882 0.0111
#> age -0.185936 0.043566 -0.271324 -0.100547 0.0000
#>
#> exp(coeffients):
#> Estimate 2.5% 97.5%
#> (Intercept) 20.37521 17.93400 23.1487
#> tcell 1.11241 0.84834 1.4587
#> platelet 1.28031 1.05794 1.5494
#> age 0.83033 0.76237 0.9043
#>
#>
## weighted GLM version RMST
out2 <- logitIPCW(Event(time,cause!=0)~tcell+platelet+age,bmt,
time=50,cens.model=~strata(platelet),model="exp",outcome="rmst")
summary(out2)
#> n events
#> 408 7467.647
#>
#> 408 clusters
#> coeffients:
#> Estimate Std.Err 2.5% 97.5% P-value
#> (Intercept) 3.020093 0.074874 2.873343 3.166844 0.0000
#> tcell 0.024791 0.186856 -0.341440 0.391021 0.8945
#> platelet 0.253188 0.126276 0.005691 0.500684 0.0450
#> age -0.180332 0.052644 -0.283513 -0.077151 0.0006
#>
#> exp(coeffients):
#> Estimate 2.5% 97.5%
#> (Intercept) 20.49320 17.69607 23.7325
#> tcell 1.02510 0.71075 1.4785
#> platelet 1.28812 1.00571 1.6498
#> age 0.83499 0.75313 0.9257
#>
#>
### time-lost
outtl <- resmeanIPCW(Event(time,cause!=0)~tcell+platelet+age,bmt,
time=50,cens.model=~strata(platelet),model="exp",outcome="rmtl")
summary(outtl)
#> n events
#> 408 245
#>
#> 408 clusters
#> coeffients:
#> Estimate Std.Err 2.5% 97.5% P-value
#> (Intercept) 3.361844 0.047699 3.268356 3.455332 0.0000
#> tcell -0.085239 0.127476 -0.335087 0.164609 0.5037
#> platelet -0.239938 0.100766 -0.437436 -0.042440 0.0173
#> age 0.175605 0.044854 0.087693 0.263518 0.0001
#>
#> exp(coeffients):
#> Estimate 2.5% 97.5%
#> (Intercept) 28.84233 26.26811 31.6688
#> tcell 0.91829 0.71528 1.1789
#> platelet 0.78668 0.64569 0.9584
#> age 1.19197 1.09165 1.3015
#>
#>
### same as Kaplan-Meier for full censoring model
bmt$int <- with(bmt,strata(tcell,platelet))
out <- resmeanIPCW(Event(time,cause!=0)~-1+int,bmt,time=30,
cens.model=~strata(platelet,tcell),model="lin")
estimate(out)
#> Estimate Std.Err 2.5% 97.5% P-value
#> inttcell=0, platelet=0 13.60 0.8316 11.97 15.23 3.830e-60
#> inttcell=0, platelet=1 18.90 1.2696 16.41 21.39 4.000e-50
#> inttcell=1, platelet=0 16.19 2.4061 11.48 20.91 1.705e-11
#> inttcell=1, platelet=1 17.77 2.4537 12.96 22.58 4.466e-13
out1 <- phreg(Surv(time,cause!=0)~strata(tcell,platelet),data=bmt)
rm1 <- resmean.phreg(out1,times=30)
summary(rm1)
#> strata times rmean se.rmean lower upper
#> tcell=0, platelet=0 0 30 13.60291 0.8315433 12.06697 15.33436
#> tcell=0, platelet=1 1 30 18.90127 1.2693294 16.57021 21.56026
#> tcell=1, platelet=0 2 30 16.19119 2.4006195 12.10803 21.65129
#> tcell=1, platelet=1 3 30 17.76601 2.4422170 13.56997 23.25953
#> years.lost
#> tcell=0, platelet=0 16.39709
#> tcell=0, platelet=1 11.09873
#> tcell=1, platelet=0 13.80881
#> tcell=1, platelet=1 12.23399
### years lost regression
outl <- resmeanIPCW(Event(time,cause!=0)~-1+int,bmt,time=30,outcome="years-lost",
cens.model=~strata(platelet,tcell),model="lin")
estimate(outl)
#> Estimate Std.Err 2.5% 97.5% P-value
#> inttcell=0, platelet=0 16.40 0.8315 14.767 18.03 1.485e-86
#> inttcell=0, platelet=1 11.10 1.2693 8.611 13.59 2.255e-18
#> inttcell=1, platelet=0 13.81 2.4006 9.104 18.51 8.810e-09
#> inttcell=1, platelet=1 12.23 2.4423 7.447 17.02 5.467e-07
## competing risks years-lost for cause 1
out <- resmeanIPCW(Event(time,cause)~-1+int,bmt,time=30,cause=1,
cens.model=~strata(platelet,tcell),model="lin")
estimate(out)
#> Estimate Std.Err 2.5% 97.5% P-value
#> inttcell=0, platelet=0 12.105 0.8508 10.438 13.773 6.162e-46
#> inttcell=0, platelet=1 6.884 1.1741 4.583 9.185 4.536e-09
#> inttcell=1, platelet=0 7.261 2.3533 2.648 11.873 2.033e-03
#> inttcell=1, platelet=1 5.780 2.0925 1.679 9.882 5.737e-03
## same as integrated cumulative incidence
rmc1 <- cif.yearslost(Event(time,cause)~strata(tcell,platelet),data=bmt,times=30)
summary(rmc1)
#> $estimate
#> strata times intF_1 intF_2 se.intF_1 se.intF_2
#> tcell=0, platelet=0 0 30 12.105148 4.291941 0.8508118 0.6161452
#> tcell=0, platelet=1 1 30 6.884173 4.214557 1.1741008 0.9057047
#> tcell=1, platelet=0 2 30 7.260728 6.548085 2.3532777 1.9703499
#> tcell=1, platelet=1 3 30 5.780427 6.453559 2.0925168 2.0815304
#> total.years.lost lower_intF_1 upper_intF_1 lower_intF_2
#> tcell=0, platelet=0 16.39709 10.547349 13.893028 3.239339
#> tcell=0, platelet=1 11.09873 4.928091 9.616672 2.765855
#> tcell=1, platelet=0 13.80881 3.846776 13.704509 3.630644
#> tcell=1, platelet=1 12.23399 2.843312 11.751551 3.429674
#> upper_intF_2
#> tcell=0, platelet=0 5.686579
#> tcell=0, platelet=1 6.422061
#> tcell=1, platelet=0 11.809864
#> tcell=1, platelet=1 12.143553
#>