# Time-varying coefficients based on binary variables

I am doing a survival analysis, and for each individual i have one event occuring (the one of interest), then i have baseline events (inc1 & inc2), so no problem for those ones. But then i have events that can occur during the follow-up (fu1 & fu2) at anytime (or before or after the event of interest). And i am testing if the group they are part of has an impact. For each event occuring i have the time of when it occured. And my purpose is to use the other events as time-varying but i don't know how to incorporate them in the cox model properly. Here is an example of the data:

d <- structure(list(id = 1:10, event = c(0, 0, 0, 0, 0, 0, 0, 1, 1, 0), group=c(1, 2, 1, 2, 1, 2, 1, 2, 1, 2),
time = c(4855, 3693, 4542, 3978, 3809, 4216, 3858, 1068, 214, 4659),
inc1 = c(0, 1, 0, 0, 0, 0, 0, 1, 0, 0), inc2 = c(0, 0, 1, 0, 0, 0, 0, 0, 0, 0),
fu1_event = c(0, 0, 1, 1, 1, 0, 1, 1, 1, 0),
fu1_time = c(4855, 3693, 4542, 3978, 3809, 4216, 3858, 4207, 1827, 4659),
fu2_event = c(1, 1, 1, 0, 0, 0, 0, 0, 1, 0),
fu2_time = c(4855, 3693, 4542, 3978, 3809, 4216, 3858, 4207, 3997, 4659)),
row.names = c(NA, 10L), class = "data.frame")
d
id event group time inc1 inc2 fu1_event fu1_time fu2_event fu2_time
1   1     0     1 4855    0    0         0     4855         1     4855
2   2     0     2 3693    1    0         0     3693         1     3693
3   3     0     1 4542    0    1         1     4542         1     4542
4   4     0     2 3978    0    0         1     3978         0     3978
5   5     0     1 3809    0    0         1     3809         0     3809
6   6     0     2 4216    0    0         0     4216         0     4216
7   7     0     1 3858    0    0         1     3858         0     3858
8   8     1     2 1068    1    0         1     4207         0     4207
9   9     1     1  214    0    0         1     1827         1     3997
10 10     0     2 4659    0    0         0     4659         0     4659


With this data my first model without adjustment would be :

library(survival)
coxph(Surv(time, event)~group, data=d)


And with baseline adjustments it would be :

coxph(Surv(time, event)~group+inc1+inc2, data=d)


My question is how to add the fu1_event and the fu2_event properly ? First of all i surely got to put them to 0 if they appear after the event of interest (or the same day), but that is kind of an example of the raw data. I saw something like taking every unique time of event and then split the database with those as timepoints, but i don't know if it's a really valid way to do