# Part I - Analysis of recurrent events

### Objective

To present four variance-corrected Cox bases models for analysing recurrent event data:

• The Andersen-Gill counting process
• The Prentice-Williams-Peterson (PWP) total time model - Conditional A
• PWP gap time model - Conditional B
• Marginal (Wei, Lin and Weissfeld) model

Those approaches do differ in how the risk set is determined but use the same COX PH model available in R or any other statistical software. The emphasis will be the creation of appropriate datasets. I will present each of the models using the psoriasis data set available at UCLA: Statistical Consulting Group.

This data set is initially described and used in the Applied Survival Analysis: Regression Modeling of Time-to-Event Data book by David W. Hosmer. All the datasets can be downloaded here.

### Psoriasis data

The data comes from (Kabat-Zinn J 1998) study of the influence of mindfullness-based stress reduction intervention in the treatment of psoriasis.

Variable name Description
id Patient ID
tape 1 - recieved instruction on the use of mindfullness-based stress reduction techniques
0 - did not receive any aspect of the mindfullness-based intervention
light 0 - phototherap
1 - photochemotherapy
yrspsor number of years the patient had psoriasis
time2 time to event
status Four different endpoints:
1: Number of days until the first response to treatment was noted by study personnel
2: Number of days until a specific major lesion showed evidence of treatment - the turning point
3: Number of days until the key lesion was half-cleared - the half-way point
4: number of days until the key lesion was completely cleared - the clearing point
> # extracting the data
> names(mydata) <- tolower(names(mydata))
>
> # first look at the data
> glimpse(mydata)
Observations: 126
Variables: 6
$id <dbl> 1, 1, 1, 1, 3, 3, 4, 4, 4, 4, 5, 6, 6, 6, 6, 8, 8, 8, 8,…$ tape    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1,…
$light <dbl> 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1,…$ yrspsor <dbl> 12, 12, 12, 12, 2, 2, 1, 1, 1, 1, 6, 22, 22, 22, 22, 4, …
$time2 <dbl> 14, 25, 44, 74, 30, 47, 5, 14, 21, 117, 21, 10, 24, 34, …$ status  <dbl> 1, 1, 1, 1, 1, 0, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1,…
> #summary(mydata)
>
> # id numbers with missing yrspsor data
> id_drop <- mydata %>%
+   filter(is.na(yrspsor)) %>%
+   distinct(id) %>%
+   pull()
>
> mydata <- mydata %>%
+   filter(!id %in% id_drop) # removing those id:s from the data set

In the psoriasis study, the follow-up times for the 32 patients is divided into 110 time intervals, of which 96 ended in an event. As an example, here is the original data structure for two patients with id = 1 and id = 3. Each patient is represented by several rows of data depending on the number of endpoints a patient had, maximum up to 4 recurrent events.

id tape light yrspsor time2 status
1 0 1 12 14 1
1 0 1 12 25 1
1 0 1 12 44 1
1 0 1 12 74 1
3 0 0 2 30 1
3 0 0 2 47 0

### Assumptions

All the four models can handle time-varying covariates. To keep things simple, we assume that the psoriases data involves only time-independent covariates. As usual, the PH assumption needs to be evaluated. Even here we will just assume that PH assumption is satisfied. The goal is not to present a correct analysis of the data, but rather to describe the usage of recurrent event models.

If you are interested, there are two posts available on my blog related to Cox PH with Time varying Covariate, post 1 and post 2.

### Counting process model

The counting process model was introduced by (Andersen 1993). It uses time defined from the beginning of the study and treats the recurrent events as being independent. The number of recurrence is not taken into account. Also, the model ignores the order of occurrence of the events. The counting process model is recommended, if it is reasonable to assume that the risk of a recurrent event remains constant, irrespective of the number of previous events.

In the psoriasis data, we need to create a time1 variable that will be the start point for each interval within each specific patient.

> # time1 are the start times
> # time2 are the stop times
> mydata_cp <- mydata %>%
+   group_by(id) %>%
+   mutate(time1 = lag(time2)) %>%
+   mutate(time1 = if_else(is.na(time1), 0, time1)) %>%
+   ungroup()
id tape light yrspsor time1 time2 status
1 0 1 12 0 14 1
1 0 1 12 14 25 1
1 0 1 12 25 44 1
1 0 1 12 44 74 1
3 0 0 2 0 30 1
3 0 0 2 30 47 0

The data for id = 1 could be described as data for four different patients. This patient contributes to the risk set at 14, 25, 44 and 74 days. The id = 3 patient has two independent intervals: an event on day 30 and is followed until day 47 and censored at that time. Thus, a patient who has experienced an event remains under risk for further events.

The no interaction Cox PH Model is: $\small \mathsf{ h(t, \mathbf{X}) = h_0(t) \exp\Big[ \beta_1 \text{ tape} + \beta_2 \text{ light} + \beta_3 \text{ yrspsor} \Big] }$ The counting process model uses the partial likelihood function expressed as the product of individual likelihoods contributed by each 96 ordered failure time.

$\mathbf{L} = \prod_{j=1}^{96} L_j$ Each individual likelihood $$L_j$$ is the conditional probability of failing at time $$t_{(j)}$$, given survival (i.e. remaining in the risk set) at $$t_{(j)}$$.

### Robust variance estimator

As mentioned, the Cox model treats each interval from the same patient as if they were independent contributions from different patients. Nevertheless, those observations are correlated, and we need to adjust for that in the analysis

(Lin and Wei 1989) proposed a robust estimator for recurrent event data. It is an extension similar to the information sandwich estimator proposed by (Zeger and Liang 1986) for generalized linear models. The approach involves adjusting the estimated variance of regression coefficients. The general formula is:

$\mathbf{\hat R (\hat \beta)} = \mathbf{\widehat{Var} (\hat \beta)} [ \mathbf{R_S' R_S} ] \mathbf{\widehat{Var} (\hat \beta)} \text{, where}$

$$\mathbf{\widehat{Var} (\hat \beta)}$$ is the information matrix of estimated variances and covariances obtained from partial ML estimation of the Cox model being fit. The $$\mathbf{R_S}$$ is the matrix of score residuals also obtained from ML estimator. The robust variance estimation applies to all the models described in this post.

> # Counting process model
> fit_cp <- coxph(Surv(time1, time2, status) ~ tape + light + yrspsor + cluster(id),
+                 data=mydata_cp,
+                 method = "breslow")
>
> # Use broom package to tidy up the output
> tidy(fit_cp, exponentiate = TRUE) %>%
+   mutate_if(is.numeric, round, 3) %>%
+   kable() %>%
+   kable_styling("striped", full_width = F, font_size = 14)
term estimate std.error robust.se statistic p.value conf.low conf.high
tape 1.477 0.216 0.124 3.135 0.002 1.157 1.885
light 1.494 0.212 0.121 3.312 0.001 1.178 1.894
yrspsor 1.009 0.012 0.006 1.612 0.107 0.998 1.020
> # Variance-covariance matrix
> vcov(fit_cp, complete=TRUE)
tape        light       yrspsor
tape     1.548908e-02 6.480624e-03 -3.843064e-05
light    6.480624e-03 1.468506e-02  9.063947e-05
yrspsor -3.843064e-05 9.063947e-05  3.113399e-05
>
> # Square-root of variance for tape = Robust SE of the estimated coefficien of the tape variable
> sqrt(coxVarCov(fit_cp)[1, 1])
[1] 0.1244551

Note that robust standard errors are smaller than the nonrobust values.

### Conditional A and B models

The both models use stratified Cox-based approaches and take into account the order of events but define the time differently. The first considers the time since study entry (total time) while the second incorporates the time since the previous event (gap time). In other words, Conditional A uses the time from the beginning of follow-up; meanwhile, Conditional B resets the time to 0 after an event is observed. Both models only include patients who have experienced the previous $$(j-1)$$ events in the risk set of the $$j$$th event.

In the total time approach, Conditional A model, the hazard for an individual $$i$$ for the $$j$$th recurrent event is:

$\mathsf{ h_{ij}(t) = h_{0j}(t) \exp \Big( \beta_j X_{ij} \Big) \text{, where} }$

$i = 1, \dots, n \quad \small j = 1, \dots, k_i \quad k_i \leq k$ In the gap time approach, Conditional B model, the hazard will be:

$\mathsf{ h_{ij}(t-t_{j-1}) = h_{0j}(t) \exp \Big( \beta_j X_{ij} \Big) \text{, where} }$

$i = 1, \dots, n \quad \small j = 1, \dots, k_i \quad k_i \leq k$ $$h_{0j}$$ is the event-specific baseline hazard for the $$j$$th event. This is achieved by stratifying on the order of the events.

> # Conditional model A
> mydata_condA <- mydata_cp %>%
+   arrange(id, time2) %>%
+   group_by(id) %>%
+   mutate(enum=row_number()) %>% #event count order
+   ungroup()
>
> # Conditional model B
> mydata_condB <- mydata_cp %>%
+   arrange(id, time2) %>%
+   group_by(id) %>%
+   mutate(enum=row_number()) %>% #event count order
+   mutate(time2=time2-time1) %>%
+   mutate(time1=0) %>%
+   ungroup()

Lets ones again look at the patient with id = 1 and id = 3 under the Conditional A and B models. The difference is in how the time1 and time2 are defined.

id tape light yrspsor time1 time2 status enum
1 0 1 12 0 14 1 1
1 0 1 12 14 25 1 2
1 0 1 12 25 44 1 3
1 0 1 12 44 74 1 4
3 0 0 2 0 30 1 1
3 0 0 2 30 47 0 2
id tape light yrspsor time1 time2 status enum
1 0 1 12 0 14 1 1
1 0 1 12 0 11 1 2
1 0 1 12 0 19 1 3
1 0 1 12 0 30 1 4
3 0 0 2 0 30 1 1
3 0 0 2 0 17 0 2

Under Conditional A model, the first patient, id = 1, enters the risk set for the second event at day 14 and leaves at day 25, while the patient with id = 3 does not enter the risk set until day 30 and is censored at day 47. Under the Conditional B model, the second event time is reset to 0. Consequently, both patients are in the risk set at the same time, from day 0 to day 11. Conditional A and B models use different entry and removal times from the risk set, but also the length of time the patients are in the risk sets together are different.

#### When use model A or model B?

Use Conditional A model if you are interested in modelling the full-time course of the recurrent event. Conditional B model can be used when the goal is to model the gap time between events. In other words, in situations when a patient may leave the risk set (e.g., lost to follow-up) but then re-enters the risk set again at some later point in time. For more details, consider reviewing the following publication (R. L. Prentice 1981). The fitted model in R is the same both for Conditional A and B models, but as already described, the risk sets do differ.

> # Conditional A model
> fit_a <- coxph(Surv(time1, time2, status) ~ tape + light + yrspsor + strata(enum) + cluster(id),
+                 data=mydata_condA, # mydata_condA is used
+                 method = "breslow")
>
> # Conditional B model
> fit_b <- coxph(Surv(time1, time2, status) ~ tape + light + yrspsor + strata(enum) + cluster(id),
+                 data=mydata_condB,  # mydata_condB is used
+                 method = "breslow") 

Outputs for the Conditional model A and B:

term estimate std.error robust.se statistic p.value conf.low conf.high
Conditional A
tape 0.838 0.247 0.250 3.355 0.001 0.349 1.328
light 1.118 0.270 0.268 4.169 0.000 0.592 1.643
yrspsor 0.021 0.013 0.014 1.539 0.124 -0.006 0.048
Conditional B
tape 1.852 0.232 0.245 2.514 0.012 1.146 2.994
light 2.251 0.232 0.224 3.630 0.000 1.452 3.488
yrspsor 1.016 0.013 0.012 1.283 0.200 0.992 1.041

### Marginal (Wei, Lin and Weissfeld) model

For each of the possible reccurent events, the marginal model considures the time from the study entry. The number of strata in the analysis will be equal to the maximum number of events in the study. Here all the time intervals will start at zero.

The hazard function for the $$j$$th recurrent event for $$i$$th individual, is:

$\mathsf{ h_{ij}(t) = h_{0j}(t) \exp \Big( \beta_j X_{ij} \Big) }$

Unlike the counting process model, this model allows a separate underlying hazard for each event. When an event is zero, a subject is no longer at risk after the last given event.

> # not the cleanest code but it works :)
>
> # keeping the last row per id, if less then 4
> last_row <- mydata_cp %>%
+   arrange(id, time2) %>%
+   group_by(id) %>%
+   mutate(enum=row_number()) %>% #event count order
+   slice(n()) %>% #keeping the last row per id
+   filter(enum!=4) %>% #remove the id if the last row is the 4th event
+   mutate(status = 0) # the status shuld be 0 in the added new rows
>
> # copying the last row 3 times - for some it will be to many since I need to only keep 4 per id
> last_row_3times <- last_row %>%
+   slice(rep(row_number(), 3))
>
> # keeping just 4 rows per id
> mydata_m  <- bind_rows(mydata_cp, last_row_3times) %>%
+   arrange(id, time2) %>%
+   group_by(id) %>%
+   mutate(enum=row_number(),
+          time1 = 0) %>%
+   filter(enum <= 4) %>%
+   ungroup()

The idea is to create 4 rows per patient since we could maximally observe 4 recurrent events. For example, consider the time to the third event for the patient 1 and 3. The first patient experienced a third event on day 44 and contributed to the risk set from 0 to 44 days. The second patient (id=3) experienced just one evet, yet under the marginal model is considered to be at risk for the second, third and fourth from time 0 to 47. The marginal model looks at each event separately and models all the available data for that specific event.

id tape light yrspsor time1 time2 status enum
1 0 1 12 0 14 1 1
1 0 1 12 0 25 1 2
1 0 1 12 0 44 1 3
1 0 1 12 0 74 1 4
3 0 0 2 0 30 1 1
3 0 0 2 0 47 0 2
3 0 0 2 0 47 0 3
3 0 0 2 0 47 0 4
> fit_m <- coxph(Surv(time1, time2, status) ~ tape + light + yrspsor + strata(enum) + cluster(id),
+                data=mydata_m,
+                method = "breslow")
>
> # Use broom package to tidy up the output
> tidy(fit_m, exponentiate = FALSE) %>%
+   mutate_if(is.numeric, round, 3) %>%
+   kable() %>%
+   kable_styling("striped", full_width = F, font_size = 14)
term estimate std.error robust.se statistic p.value conf.low conf.high
tape 1.030 0.243 0.319 3.224 0.001 0.404 1.656
light 1.540 0.254 0.307 5.019 0.000 0.939 2.142
yrspsor 0.028 0.013 0.019 1.488 0.137 -0.009 0.064

Remember that the marginal model includes all patients in each stratum; it will overestimates the covariate effects. As it is explained by (A. K. Ozga 2018)

The strata-specific effect estimators generally increase in magnitude with time as later strata in the treatment group contain more censored observations and thus the influence of a single event in the control group becomes larger resulting in an exaggerated treatment effect over time.
Furthermore, the standard deviation is higher in the Wei-Lin-Weissfeld model. This can be explained by the fact that no direct global effect but strata-specific effects and variances are estimated which are combined subsequently. As later strata contain fewer events, the strata-specific standard deviation increases with time.

#### In summary

The choice of appropriate model will be influenced by many factors:

• number of reccurent events,
• relationship among subsequence events,
• within subject correlation and varying covariates, and
• sample size.

The most appropriate model should be chosen based on the nature of the data.

#### Reproducibility

─ Session info ───────────────────────────────────────────────────────────────────────────────────────────────────────
setting  value
version  R version 3.5.3 (2019-03-11)
os       macOS Mojave 10.14.6
system   x86_64, darwin15.6.0
ui       X11
language (EN)
collate  en_US.UTF-8
ctype    en_US.UTF-8
tz       Europe/Stockholm
date     2019-11-10

─ Packages ───────────────────────────────────────────────────────────────────────────────────────────────────────────
package        * version    date       lib source
abind            1.4-5      2016-07-21 [1] CRAN (R 3.5.0)
acepack          1.4.1      2016-10-29 [1] CRAN (R 3.5.0)
assertthat       0.2.1      2019-03-21 [1] CRAN (R 3.5.2)
backports        1.1.5      2019-10-02 [1] CRAN (R 3.5.2)
base64enc        0.1-3      2015-07-28 [1] CRAN (R 3.5.0)
blogdown         0.16       2019-10-01 [1] CRAN (R 3.5.2)
bookdown         0.14       2019-10-01 [1] CRAN (R 3.5.2)
broom          * 0.5.2      2019-04-07 [1] CRAN (R 3.5.2)
callr            3.3.2      2019-09-22 [1] CRAN (R 3.5.2)
cellranger       1.1.0      2016-07-27 [1] CRAN (R 3.5.0)
checkmate        1.9.4      2019-07-04 [1] CRAN (R 3.5.2)
cli              1.1.0      2019-03-19 [1] CRAN (R 3.5.2)
cluster          2.1.0      2019-06-19 [1] CRAN (R 3.5.2)
cmprsk           2.2-9      2019-10-09 [1] CRAN (R 3.5.2)
codetools        0.2-16     2018-12-24 [1] CRAN (R 3.5.3)
colorspace       1.4-1      2019-03-18 [1] CRAN (R 3.5.2)
crayon           1.3.4      2017-09-16 [1] CRAN (R 3.5.0)
data.table     * 1.12.2     2019-04-07 [1] CRAN (R 3.5.2)
desc             1.2.0      2018-05-01 [1] CRAN (R 3.5.0)
devtools       * 2.2.1      2019-09-24 [1] CRAN (R 3.5.2)
digest           0.6.21     2019-09-20 [1] CRAN (R 3.5.2)
dplyr          * 0.8.3      2019-07-04 [1] CRAN (R 3.5.2)
ellipsis         0.3.0      2019-09-20 [1] CRAN (R 3.5.2)
evaluate         0.14       2019-05-28 [1] CRAN (R 3.5.2)
fansi            0.4.0      2018-10-05 [1] CRAN (R 3.5.0)
forcats        * 0.4.0      2019-02-17 [1] CRAN (R 3.5.2)
foreach          1.4.7      2019-07-27 [1] CRAN (R 3.5.2)
foreign          0.8-72     2019-08-02 [1] CRAN (R 3.5.2)
Formula          1.2-3      2018-05-03 [1] CRAN (R 3.5.0)
fs               1.3.1      2019-05-06 [1] CRAN (R 3.5.2)
generics         0.0.2      2018-11-29 [1] CRAN (R 3.5.0)
ggplot2        * 3.2.1      2019-08-10 [1] CRAN (R 3.5.2)
glue             1.3.1.9000 2019-10-12 [1] Github (tidyverse/glue@71eeddf)
gridExtra        2.3        2017-09-09 [1] CRAN (R 3.5.0)
gtable           0.3.0      2019-03-25 [1] CRAN (R 3.5.2)
haven            2.1.1      2019-07-04 [1] CRAN (R 3.5.2)
highr            0.8        2019-03-20 [1] CRAN (R 3.5.2)
Hmisc            4.2-0      2019-01-26 [1] CRAN (R 3.5.2)
hms              0.5.1      2019-08-23 [1] CRAN (R 3.5.2)
htmlTable        1.13.2     2019-09-22 [1] CRAN (R 3.5.2)
htmltools        0.4.0      2019-10-04 [1] CRAN (R 3.5.2)
htmlwidgets      1.5.1      2019-10-08 [1] CRAN (R 3.5.2)
httr             1.4.1      2019-08-05 [1] CRAN (R 3.5.2)
iterators        1.0.12     2019-07-26 [1] CRAN (R 3.5.2)
jsonlite         1.6        2018-12-07 [1] CRAN (R 3.5.0)
kableExtra     * 1.1.0      2019-03-16 [1] CRAN (R 3.5.2)
knitr          * 1.25       2019-09-18 [1] CRAN (R 3.5.2)
lattice          0.20-38    2018-11-04 [1] CRAN (R 3.5.3)
latticeExtra     0.6-28     2016-02-09 [1] CRAN (R 3.5.0)
lava             1.6.6      2019-08-01 [1] CRAN (R 3.5.2)
lazyeval         0.2.2      2019-03-15 [1] CRAN (R 3.5.2)
lifecycle        0.1.0      2019-08-01 [1] CRAN (R 3.5.2)
lubridate        1.7.4      2018-04-11 [1] CRAN (R 3.5.0)
magrittr         1.5        2014-11-22 [1] CRAN (R 3.5.0)
MASS             7.3-51.4   2019-03-31 [1] CRAN (R 3.5.2)
Matrix           1.2-17     2019-03-22 [1] CRAN (R 3.5.2)
MatrixModels     0.4-1      2015-08-22 [1] CRAN (R 3.5.0)
memoise          1.1.0      2017-04-21 [1] CRAN (R 3.5.0)
modelr           0.1.5      2019-08-08 [1] CRAN (R 3.5.2)
multcomp         1.4-10     2019-03-05 [1] CRAN (R 3.5.2)
munsell          0.5.0      2018-06-12 [1] CRAN (R 3.5.0)
mvtnorm          1.0-11     2019-06-19 [1] CRAN (R 3.5.2)
nlme             3.1-141    2019-08-01 [1] CRAN (R 3.5.2)
nnet             7.3-12     2016-02-02 [1] CRAN (R 3.5.3)
numDeriv         2016.8-1.1 2019-06-06 [1] CRAN (R 3.5.2)
pillar           1.4.2      2019-06-29 [1] CRAN (R 3.5.2)
pkgbuild         1.0.6      2019-10-09 [1] CRAN (R 3.5.2)
pkgconfig        2.0.3      2019-09-22 [1] CRAN (R 3.5.2)
pkgload          1.0.2      2018-10-29 [1] CRAN (R 3.5.0)
polspline        1.1.16     2019-10-04 [1] CRAN (R 3.5.2)
prettyunits      1.0.2      2015-07-13 [1] CRAN (R 3.5.0)
processx         3.4.1      2019-07-18 [1] CRAN (R 3.5.2)
prodlim        * 2018.04.18 2018-04-18 [1] CRAN (R 3.5.0)
ps               1.3.0      2018-12-21 [1] CRAN (R 3.5.0)
purrr          * 0.3.2      2019-03-15 [1] CRAN (R 3.5.2)
quantreg         5.51       2019-08-07 [1] CRAN (R 3.5.2)
R6               2.4.0      2019-02-14 [1] CRAN (R 3.5.2)
RColorBrewer     1.1-2      2014-12-07 [1] CRAN (R 3.5.0)
Rcpp             1.0.2      2019-07-25 [1] CRAN (R 3.5.2)
readr          * 1.3.1      2018-12-21 [1] CRAN (R 3.5.0)
readxl           1.3.1      2019-03-13 [1] CRAN (R 3.5.2)
remotes          2.1.0      2019-06-24 [1] CRAN (R 3.5.2)
riskRegression * 2019.01.29 2019-01-29 [1] CRAN (R 3.5.2)
rlang            0.4.0      2019-06-25 [1] CRAN (R 3.5.2)
rmarkdown        1.16       2019-10-01 [1] CRAN (R 3.5.2)
rms              5.1-3.1    2019-04-22 [1] CRAN (R 3.5.2)
rpart            4.1-15     2019-04-12 [1] CRAN (R 3.5.2)
rprojroot        1.3-2      2018-01-03 [1] CRAN (R 3.5.0)
rstudioapi       0.10       2019-03-19 [1] CRAN (R 3.5.2)
rvest            0.3.4      2019-05-15 [1] CRAN (R 3.5.2)
sandwich         2.5-1      2019-04-06 [1] CRAN (R 3.5.2)
sas7bdat       * 0.5        2014-06-04 [1] CRAN (R 3.5.0)
scales           1.0.0      2018-08-09 [1] CRAN (R 3.5.0)
sessioninfo      1.1.1      2018-11-05 [1] CRAN (R 3.5.0)
SparseM          1.77       2017-04-23 [1] CRAN (R 3.5.0)
stringi          1.4.3      2019-03-12 [1] CRAN (R 3.5.2)
stringr        * 1.4.0      2019-02-10 [1] CRAN (R 3.5.2)
survival       * 2.44-1.1   2019-04-01 [1] CRAN (R 3.5.2)
testthat         2.2.1      2019-07-25 [1] CRAN (R 3.5.2)
TH.data          1.0-10     2019-01-21 [1] CRAN (R 3.5.2)
tibble         * 2.1.3      2019-06-06 [1] CRAN (R 3.5.2)
tidyr          * 1.0.0      2019-09-11 [1] CRAN (R 3.5.2)
tidyselect       0.2.5      2018-10-11 [1] CRAN (R 3.5.0)
tidyverse      * 1.2.1      2017-11-14 [1] CRAN (R 3.5.0)
timereg          1.9.4      2019-07-30 [1] CRAN (R 3.5.2)
usethis        * 1.5.1      2019-07-04 [1] CRAN (R 3.5.2)
utf8             1.1.4      2018-05-24 [1] CRAN (R 3.5.0)
vctrs            0.2.0      2019-07-05 [1] CRAN (R 3.5.2)
viridisLite      0.3.0      2018-02-01 [1] CRAN (R 3.5.0)
webshot          0.5.1      2018-09-28 [1] CRAN (R 3.5.0)
withr            2.1.2      2018-03-15 [1] CRAN (R 3.5.0)
xfun             0.10       2019-10-01 [1] CRAN (R 3.5.2)
xml2             1.2.2      2019-08-09 [1] CRAN (R 3.5.2)
yaml             2.2.0      2018-07-25 [1] CRAN (R 3.5.0)
zeallot          0.1.0      2018-01-28 [1] CRAN (R 3.5.0)
zoo              1.8-6      2019-05-28 [1] CRAN (R 3.5.2)

[1] /Library/Frameworks/R.framework/Versions/3.5/Resources/library

#### References

A. K. Ozga, G. Rauch, M. Kieser. 2018. “A Systematic Comparison of Recurrent Event Models for Application to Composite Endpoints.” BMC Medical Research Methodology 18 (1): 2. https://doi.org/10.1186/s12874-017-0462-x.

Andersen, Borgan, P.K. 1993. Statistical Models Based on Counting Processes. Springer-Verlag.

Kabat-Zinn J, Light T, Wheeler E. 1998. “Influence of a Mindfulness Meditation-Based Stress Reduction Intervention on Rates of Skin Clearing in Patients with Moderate to Severe Psoriasis Undergoing Phototherapy (Uvb) and Photochemotherapy (Puva).” Psychosomatic Medicine 60 (5): 625–32. https://doi.org/10.1097/00006842-199809000-00020.

Lin, D. Y., and L. J. Wei. 1989. “The Robust Inference for the Cox Proportional Hazards Model.” Journal of the American Statistical Association 84 (408): 1074–8. https://doi.org/10.2307/2290085.

R. L. Prentice, A. V. Peterson, B. J. Williams. 1981. “On the Regression Analysis of Multivariate Failure Time Data.” Biometrika 68 (2): 373–79. https://doi.org/10.1093/biomet/68.2.373.

Zeger, Scott L., and Kung-Yee Liang. 1986. “Longitudinal Data Analysis for Discrete and Continuous Outcomes.” Biometrics 42 (1): 121–30. https://doi.org/10.2307/2531248.