CoxPH with one ID contributing many comparator observations

Question

UPDATE-----------------

I have data where one ID can serve as a comparator several times.

This is done by creating a variable to group the treated in variable GROUP. Assume variable GROUP, as an example, containing the number of visits to a physician until a therapy is initiated. Thus, ID 4 initiated the therapy after 3 visits. A possible comparator is selected if it also has 3 physician visits but does not initiate the therapy. Therefore, both belong to GROUP == 3. Note that this could happen on different points in time specified in variable TIME_POINT. As ID 1 never starts therapy it contributes also a possible comparator observation for someone who initiated therapy after 5 physician visits in GROUP == 5. Next, within every value of variable GROUP, the covariates are collected in a time window before the TIME_POINT on which an ID belongs to a GROUP.
In this minimum example, assume covariates to be balanced using a weight specified in variable WEIGHT. In addition, the data outside the minimum example has several comparators per GROUP.

Kept short:
ID 1 is a comparator as it is untreated (TREAT == 0) and has one event in total.
ID 2, 3, 4, 5 and 6 are different from each other and were treated (TREAT == 1).

ID 1 therefore contributes to different GROUPs as a comparator for ID 2, 3, 4, 5 and 6.

ID 1 has one event in total.
Accordingly, the duration of the follow-up time (Variable FOLLOW_UP_TIME) naturally shortens as value of GROUP and the associated variable TIME_POINT increases.

Data:

data <- data.frame(
  ID = c(1, 1, 1, 1, 1, 2, 3, 4, 5, 6),
  TREAT = c(0, 0, 0, 0, 0, 1, 1, 1, 1, 1),
  GROUP = c(1, 2, 3, 4, 5, 1, 2, 3, 4, 5),
  TIME_POINT = as.Date(c("2019-05-13", "2019-06-01", "2020-03-01", "2021-02-01", "2022-09-01", "2021-02-01", "2021-04-01", "2021-07-01", "2021-01-13", "2021-01-01")),
  EVENT_DATE = as.Date(c("2022-10-01", "2022-10-01", "2022-10-01", "2022-10-01", "2022-10-01", "2022-01-03", "2021-10-03", NA, NA, NA)),
  END_DATE = as.Date(c("2022-10-01", "2022-10-01", "2022-10-01", "2022-10-01", "2022-10-01", "2022-01-03", "2022-10-03", "2021-10-04", "2021-09-03", "2021-11-02")),
  EVENT = c(1, 1, 1, 1, 1, 1, 1, 0, 0, 0),
  FOLLOW_UP_TIME = c(1237, 1218, 944, 607, 30, 336, 185, 95, 233, 305),
  WEIGHT = c(0.834, 0.543, 0.234, 0.743, 0.323, 1, 1, 1, 1, 1)
)

ID TREAT GROUP TIME_POINT EVENT_DATE   END_DATE EVENT FOLLOW_UP_TIME WEIGHT

1     0     1 2019-05-13 2022-10-01 2022-10-01  1     1237            0.834
1     0     2 2019-06-01 2022-10-01 2022-10-01  1     1218            0.543
1     0     3 2020-03-01 2022-10-01 2022-10-01  1     944             0.234
1     0     4 2021-02-01 2022-10-01 2022-10-01  1     607             0.743
1     0     5 2022-09-01 2022-10-01 2022-10-01  1     30              0.323


2     1     1 2021-02-01 2022-01-03 2022-01-03  1     336                1
3     1     2 2021-04-01 2021-10-03 2021-10-03  1     185                1
4     1     3 2021-07-01       <NA> 2021-10-04  0     95                 1
5     1     4 2021-01-13       <NA> 2021-09-03  0     233                1
6     1     5 2021-01-01       <NA> 2021-11-02  0     305                1

Example of GROUP == 5:

ID TREAT GROUP TIME_POINT EVENT_DATE   END_DATE EVENT FOLLOW_UP_TIME WEIGHT
1     0     5 2022-09-01 2022-10-01 2022-10-01  1     30              0.323
6     1     5 2021-01-01       <NA> 2021-11-02  0     305                1

Since an ID can occur several times as a comparator, I use bootstrapped standard errors for the cox model and tried to specify the cluster argument to account for the different time points of comparability, specified in variable GROUP:

library(survival)
library(boot)

r <- 1000

cox_fit <- function(data, indices) {
  subset_data <- data[indices, ]
  cox_model <- coxph(Surv(FOLLOW_UP_TIME, EVENT) ~ TREAT, data = subset_data, robust = TRUE, weights = WEIGHT, cluster = GROUP)
  return(coef(cox_model))
}

set.seed(123)
bootstrap_results <- boot(data = data, statistic = cox_fit, R = r)
print(bootstrap_results, digits = 2)

conf_level <- 0.95
boot_ci <- boot.ci(bootstrap_results, type = c("bca"), conf = conf_level)
boot_ci

Model:

ph <- coxph(Surv(FOLLOW_UP_TIME, EVENT) ~ TREAT, data = data, robust = TRUE, weights = WEIGHT, cluster = GROUP)

> ph
Call:
coxph(formula = Surv(FOLLOW_UP_TIME, EVENT) ~ TREAT, data = data, 
    weights = WEIGHT, robust = TRUE, cluster = GROUP)

       coef exp(coef) se(coef) robust se    z     p
TREAT 2.053     7.792    1.950     1.396 1.47 0.141

Likelihood ratio test=1.58  on 1 df, p=0.2088
n= 10, number of events= 7 

Is the cluster argument sufficient? How to handle such situation where one contributes many times as a comparator since the model assumes 7 events?

Answer 1

Your design of the GROUP variable troubles me.

Assume variable GROUP, as an example, containing the number of visits to a physician until a therapy is initiated.

In that formulation, your ID1 case has a right-censored value for GROUP. It would have been at least 5 visits, and might have been more if the event hadn't already happened.

That also leads to the uncomfortable situation of having 5 different FOLLOW_UP_TIME values for the same individual. The Surv(FOLLOW_UP_TIME, EVENT) formula depends on a consistent choice of time = 0. In this situation that might be the time from study entry or the time when therapy began, but there should only be a single time = 0 reference for each individual.

The setup in the OP thus seems to be inconsistent. It seems to use the time when therapy began as the time = 0 reference for those who started therapy, but the date of study entry (first visit) for ID1 who never started therapy, at least in the first data row. That could lead to substantial bias.

A modification of the counting-process Cox model with Surv(start_time, stop_time, status) proposed by Lukas Lohse would be a way to handle this, if you wanted to use the date of study entry as the time = 0 reference. You would then define a row with (start_time, stop_time, status) for each individual for each interval between visits. The status would then be the status at the end of that time interval. One way to handle the initiation of therapy would be to add an indicator of whether therapy had begun at the start of the interval to each row, as a time-varying covariate. For ID1 that indicator of therapy would always be 0; for others, it would be zero for all intervals until the one at which therapy began. Other ways to represent therapy over time might be considered.

If you want to use the date of starting therapy as a reference, you could use the GROUP indicator variable to represent the number of visits before therapy was begun and the Surv(FOLLOW_UP_TIME, EVENT) Cox model form. But then you should not include individuals like ID1 who never started therapy. Your model would then be conditional upon having started therapy.

I sense additional problems with the design, as the choice of when to start therapy presumably represents an estimate of current disease status. Then the number of visits prior to therapy is perhaps just a measure of how early in the disease process an individual was first seen. You should try to get some experienced local statistical advice to help work through the possible sources of bias in this study and how best to formulate and text your hypothesis.

Answer 2

@EdM has good points in his comments and I would also question, that if TIME_POINT is study entry, them it can't really change for ID-1. But if I accept for a moment, then GROUP seems a straight forward example of a time dependent covariate, see here: https://cran.r-project.org/web/packages/survival/vignettes/timedep.pdf

The key inside here is that ID 1 was in group 1 not for 365 days but only for January 2021 and it didn't have the event while in group 5 but when it was in group 1. To implement this we need to use the start, stop Surv-notation

library(survival)
library(tidyverse)

dat_surv <- data %>%  
  group_by(ID) %>% 
  arrange(ID, TIME_POINT) %>% 
  mutate(start_time = TIME_POINT - min(TIME_POINT)) %>% # min stays in group defined by group_by
  mutate(stop_time = pmin(END_DATE - min(TIME_POINT), lead(start_time), na.rm = T)) %>% # lead stays in group, pmin is rowwise
  mutate(status = (EVENT == 1) & (TIME_POINT == max(TIME_POINT))) # only timeframe with event is marked



# unsuprisingly the minimal example with 3 events does not work well       
ph <- coxph(Surv(start_time, stop_time, status) ~ TREAT, data = dat_surv, robust = TRUE, weights = WEIGHT, cluster = GROUP)
summary(ph)

transformed data:

> dat_surv %>%  select(ID, TIME_POINT,start_time, stop_time, status)
# A tibble: 10 × 5
# Groups:   ID [6]
      ID TIME_POINT start_time stop_time status
   <dbl> <date>     <drtn>     <drtn>    <lgl> 
 1     1 2021-01-01   0 days    31 days  FALSE 
 2     1 2021-02-01  31 days    59 days  FALSE 
 3     1 2021-03-01  59 days    90 days  FALSE 
 4     1 2021-04-01  90 days   120 days  FALSE 
 5     1 2021-05-01 120 days   365 days  TRUE  
 6     2 2021-05-01   0 days   247 days  TRUE  
 7     3 2021-04-01   0 days   246 days  TRUE  
 8     4 2021-03-01   0 days   217 days  FALSE 
 9     5 2021-02-01   0 days   214 days  FALSE 
10     6 2021-01-01   0 days   305 days  FALSE