--- title: "While Alive estimands for Recurrent Events" author: Klaus Holst & Thomas Scheike date: "2026-07-29" output: rmarkdown::html_vignette: fig_caption: yes vignette: > %\VignetteIndexEntry{While Alive estimands for Recurrent Events} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- While Alive estimands for Recurrent Events ========================================== We consider two while-alive estimands for recurrent events data \begin{align*} \frac{E(N(D \wedge t))}{E(D \wedge t)} \end{align*} and the mean of the subject specific events per time-unit \begin{align*} E( \frac{N(D \wedge t)}{D \wedge t} ) \end{align*} for two treatment-groups in the case of an RCT. For the mean of events per time-unit it has been seen that when the sample size is small one can improve the finite sample properties by employing a transformation such as square or cube-root, and thus consider \begin{align*} E( (\frac{N(D \wedge t)}{D \wedge t})^.33 ) \end{align*} ``` r data(hfactioncpx12) dtable(hfactioncpx12,~status) #> #> status #> 0 1 2 #> 617 1391 124 dd <- WA_recurrent(Event(entry,time,status)~treatment+cluster(id),hfactioncpx12,time=2,death.code=2) summary(dd) #> While-Alive summaries: #> #> RMST, E(min(D,t)) #> Estimate Std.Err 2.5% 97.5% P-value #> treatment0 1.859 0.02108 1.817 1.900 0 #> treatment1 1.924 0.01502 1.894 1.953 0 #> #> Estimate Std.Err 2.5% 97.5% P-value #> [treatment0] - [treat.... -0.06517 0.02588 -0.1159 -0.01444 0.0118 #> mean events, E(N(min(D,t))): #> Estimate Std.Err 2.5% 97.5% P-value #> treatment0 1.572 0.09573 1.384 1.759 1.375e-60 #> treatment1 1.453 0.10315 1.251 1.656 4.376e-45 #> #> Estimate Std.Err 2.5% 97.5% P-value #> [treatment0] - [treat.... 0.1185 0.1407 -0.1574 0.3943 0.4 #> _______________________________________________________ #> Ratio of means E(N(min(D,t)))/E(min(D,t)) #> Estimate Std.Err 2.5% 97.5% P-value #> p1 0.8457 0.05264 0.7425 0.9488 4.411e-58 #> p2 0.7555 0.05433 0.6490 0.8619 5.963e-44 #> #> Estimate Std.Err 2.5% 97.5% P-value #> [p1] - [p2] 0.09022 0.07565 -0.05805 0.2385 0.233 #> _______________________________________________________ #> Mean of Events per time-unit E(N(min(D,t))/min(D,t)) #> Estimate Std.Err 2.5% 97.5% P-value #> treat0 1.0725 0.1222 0.8331 1.3119 1.645e-18 #> treat1 0.7552 0.0643 0.6291 0.8812 7.508e-32 #> #> Estimate Std.Err 2.5% 97.5% P-value #> [treat0] - [treat1] 0.3173 0.1381 0.04675 0.5879 0.02153 dd <- WA_recurrent(Event(entry,time,status)~treatment+cluster(id),hfactioncpx12,time=2, death.code=2,trans=.333) summary(dd,type="log") #> While-Alive summaries, log-scale: #> #> RMST, E(min(D,t)) #> Estimate Std.Err 2.5% 97.5% P-value #> treatment0 0.6199 0.011340 0.5977 0.6421 0 #> treatment1 0.6543 0.007807 0.6390 0.6696 0 #> #> Estimate Std.Err 2.5% 97.5% P-value #> [treatment0] - [treat.... -0.03446 0.01377 -0.06145 -0.007478 0.01231 #> mean events, E(N(min(D,t))): #> Estimate Std.Err 2.5% 97.5% P-value #> treatment0 0.4523 0.06090 0.3329 0.5716 1.119e-13 #> treatment1 0.3739 0.07097 0.2348 0.5130 1.376e-07 #> #> Estimate Std.Err 2.5% 97.5% P-value #> [treatment0] - [treat.... 0.07835 0.09352 -0.1049 0.2616 0.4022 #> _______________________________________________________ #> Ratio of means E(N(min(D,t)))/E(min(D,t)) #> Estimate Std.Err 2.5% 97.5% P-value #> treatment0 -0.1676 0.06224 -0.2896 -0.04563 7.081e-03 #> treatment1 -0.2804 0.07192 -0.4214 -0.13947 9.651e-05 #> #> Estimate Std.Err 2.5% 97.5% P-value #> [treatment0] - [treat.... 0.1128 0.09511 -0.07361 0.2992 0.2356 #> _______________________________________________________ #> Mean of Events per time-unit E(N(min(D,t))/min(D,t)) #> Estimate Std.Err 2.5% 97.5% P-value #> treat0 -0.3833 0.04939 -0.4801 -0.2865 8.487e-15 #> treat1 -0.5380 0.05666 -0.6491 -0.4270 2.191e-21 #> #> Estimate Std.Err 2.5% 97.5% P-value #> [treat0] - [treat1] 0.1548 0.07517 0.007459 0.3021 0.03948 ``` We see that the ratio of means is not very different between groups, but that the subject-specific mean of events per time-unit shows that those on the active treatment have fewer events per time-unit on average. We can also fit a regression model for the mean of the subject-specific events per time-unit, here using the exponential link. ``` r set.seed(10) hfactioncpx12$age <- rnorm(741)[hfactioncpx12$id] hfactioncpx12$sex <- rbinom(741,1,0.5)[hfactioncpx12$id] dd <- WA_reg(Event(entry,time,status)~treatment+age+sex+cluster(id),hfactioncpx12,time=2,death.code=2) summary(dd) #> n events #> 741 86 #> #> 741 clusters #> coeffients: #> Estimate Std.Err 2.5% 97.5% P-value #> (Intercept) -1.6366e-02 1.1187e-01 -2.3564e-01 2.0290e-01 0.8837 #> treatment1 -3.5688e-01 1.4242e-01 -6.3602e-01 -7.7743e-02 0.0122 #> age -1.0527e-01 5.3727e-02 -2.1058e-01 2.9638e-05 0.0501 #> sex 1.5723e-01 1.4784e-01 -1.3253e-01 4.4700e-01 0.2876 #> #> exp(coeffients): #> Estimate 2.5% 97.5% #> (Intercept) 0.98377 0.79007 1.2250 #> treatment1 0.69986 0.52940 0.9252 #> age 0.90008 0.81012 1.0000 #> sex 1.17027 0.87587 1.5636 ``` Composite outcomes involving death and marks ============================================ The event count can be generalised in various ways by using outcomes other than $N(D \wedge t)$, for example, \begin{align*} \tilde N(D \wedge t) = \int_0^t I(D \geq s) M(s) dN(s) + \sum_j M_j I(D \leq t,\epsilon=j) ) \end{align*} where $M(s)$ are the marks associated with $N(s)$ and $M_j$ are marks associated with the different causes of the terminal event. This provides an extension of the weighted composite outcomes measure of Mao \& Lin (2022). The marks (or weights) can be stochastic, for example when counting hospital expenses, and are stored as a vector in the data frame. The marks for the event times (defined through the causes) are then used. Here we weight death by 2 and otherwise count recurrent events as before (with weight 1): ``` r hfactioncpx12$marks <- runif(nrow(hfactioncpx12)) ##ddmg <- WA_recurrent(Event(entry,time,status)~treatment+cluster(id),hfactioncpx12,time=2, ##cause=1:2,death.code=2,marks=hfactioncpx12$marks) ##summary(ddmg) ddm <- WA_recurrent(Event(entry,time,status)~treatment+cluster(id),hfactioncpx12,time=2, cause=1:2,death.code=2,marks=hfactioncpx12$status) ``` SessionInfo ============ ``` r 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] splines stats graphics grDevices utils datasets methods #> [8] base #> #> other attached packages: #> [1] prodlim_2026.03.11 timereg_2.0.7 survival_3.8-6 mets_1.3.12 #> #> loaded via a namespace (and not attached): #> [1] cli_3.6.6 knitr_1.51 rlang_1.2.0 #> [4] xfun_0.57 KernSmooth_2.23-26 otel_0.2.0 #> [7] data.table_1.18.4 future.apply_1.20.2 listenv_0.10.1 #> [10] lava_1.9.2.1 stats4_4.6.0 grid_4.6.0 #> [13] evaluate_1.0.5 yaml_2.3.12 mvtnorm_1.3-7 #> [16] numDeriv_2016.8-1.1 compiler_4.6.0 codetools_0.2-20 #> [19] Rcpp_1.1.1-1.1 ucminf_1.2.3 future_1.70.0 #> [22] lattice_0.22-9 digest_0.6.39 parallelly_1.47.0 #> [25] parallel_4.6.0 Matrix_1.7-5 tools_4.6.0 #> [28] RcppArmadillo_15.2.6-1 globals_0.19.1 ```