conduct-ccw-analysis.RmdThis lightweight package describes how to conduct clone-censor weighting (CCW) to address the problem of immortal time bias in survival analysis. This vignette will walk through the applied tutorial published by Maringe et al 2020. Refer to Hernan and Robins 2016 and Hernan et al 2016 for more technical details.
CCW is useful in the presence of immortal person-time bias in observational studies. For instance, when comparing surgery recipients vs non-recipients in non-small cell lung cancer (NSCLC), the surgery group will have a longer survival time than the non-surgery group because the non-surgery group includes patients who died before they could receive surgery. This is a form of immortal time bias.
The CCW toy dataset published by Maringe uses this exact setting as
the motivating example. Let’s explore the dataset, which comes with
survivalCCW.
library(survivalCCW)
head(toy_df)
#> id surgery timetosurgery death sex charlson perf stage emergency fup_obs
#> 1 P5958 0 NA 0 2 1 2 2 1 365.24
#> 2 P3075 0 NA 0 2 1 0 2 0 365.24
#> 3 P1410 0 NA 0 1 1 2 1 0 365.24
#> 4 P4957 0 NA 0 2 0 1 1 0 365.24
#> 5 P7340 0 NA 1 1 0 2 2 1 123.00
#> 6 P1917 0 NA 0 2 0 0 1 0 365.24
#> age deprivation
#> 1 79.57290 0.42
#> 2 82.91855 0.07
#> 3 84.03011 0.22
#> 4 78.34634 0.12
#> 5 83.63860 0.40
#> 6 82.26147 0.03Column descriptions can be found with ?toy_df:
id: patient identifierfup_obs: observed follow-up time (time to death or 1
year if censored alive)death: observed event of interest (all-cause death) 1:
dead, 0: alivetimetosurgery: time to surgery (NA if no surgery)age: age at diagnosissex: patient’s sexperf: performance status at diagnosisstage: stage at diagnosisdeprivation: deprivation scorecharlson: Charlson’s comorbidity indexemergency: route to diagnosisNote that this package addresses situations in which the covariates are all defined at baseline.
The first step is to create the clones. This can be done for any
time-to-event outcome using the survivalCCW function
create_clones. For create_clones to work, we
need to pass a one-row-per-patient data.frame with the
following columns:
id)death). Note that additional
values are not yet permitted.fup_obs)surgery).
Must be (0) or (1),timetosurgery)All other columns will be propogated for each patient. Let’s see what this looks like in practice.
# Create clones
clones <- create_clones(
df = toy_df,
id = 'id',
event = 'death',
time_to_event = 'fup_obs',
exposure = 'surgery',
time_to_exposure = 'timetosurgery',
ced_window = 182.62
)
head(clones)
#> id surgery timetosurgery death sex charlson perf stage emergency fup_obs
#> 1 P1074 0 NA 1 2 1 2 1 0 345.00
#> 2 P1074 0 NA 1 2 1 2 1 0 345.00
#> 3 P1186 1 112 0 2 0 1 1 0 365.24
#> 4 P1186 1 112 0 2 0 1 1 0 365.24
#> 5 P1189 1 41 1 1 1 1 1 0 226.00
#> 6 P1189 1 41 1 1 1 1 1 0 226.00
#> age deprivation outcome fup_outcome censor fup_censor clone
#> 1 76.67899 0.33 1 345.00 0 182.62 0
#> 2 76.67899 0.33 0 182.62 1 182.62 1
#> 3 78.73512 0.04 0 112.00 1 112.00 0
#> 4 78.73512 0.04 0 365.24 0 112.00 1
#> 5 75.92608 0.04 0 41.00 1 41.00 0
#> 6 75.92608 0.04 1 226.00 0 41.00 1Note that this object is just a data.frame with an
additional custom class which future functions will evaluate:
class(clones)
#> [1] "ccw_clones" "data.frame"Now we simply need to cast the data to long format. The
survivalCCW function cast_to_long will do this
for us. No additional arguments are needed (the clones
object is an artifact that allows you to better see and understand the
method):
clones_long <- cast_clones_to_long(clones)
head(clones_long, row.names = FALSE)
#> time_id t_start id surgery timetosurgery death sex charlson perf stage
#> 1 0 0 P1074 0 NA 1 2 1 2 1
#> 2 1 1 P1074 0 NA 1 2 1 2 1
#> 3 2 7 P1074 0 NA 1 2 1 2 1
#> 4 3 8 P1074 0 NA 1 2 1 2 1
#> 5 4 10 P1074 0 NA 1 2 1 2 1
#> 6 5 16 P1074 0 NA 1 2 1 2 1
#> emergency fup_obs age deprivation censor fup_censor clone fup_outcome
#> 1 0 345 76.67899 0.33 0 182.62 0 1
#> 2 0 345 76.67899 0.33 0 182.62 0 7
#> 3 0 345 76.67899 0.33 0 182.62 0 8
#> 4 0 345 76.67899 0.33 0 182.62 0 10
#> 5 0 345 76.67899 0.33 0 182.62 0 16
#> 6 0 345 76.67899 0.33 0 182.62 0 18
#> outcome t_stop t_event
#> 1 0 1 NA
#> 2 0 7 1
#> 3 0 8 7
#> 4 0 10 8
#> 5 0 16 10
#> 6 0 18 16Let’s pick out a single patient and look at their data:
print(clones_long[clones_long$id == "P5913", ], row.names = FALSE)
#> time_id t_start id surgery timetosurgery death sex charlson perf stage
#> 0 0 P5913 1 1 1 1 2 0 1
#> 0 0 P5913 1 1 1 1 2 0 1
#> 1 1 P5913 1 1 1 1 2 0 1
#> 2 7 P5913 1 1 1 1 2 0 1
#> 3 8 P5913 1 1 1 1 2 0 1
#> 4 10 P5913 1 1 1 1 2 0 1
#> 5 16 P5913 1 1 1 1 2 0 1
#> 6 18 P5913 1 1 1 1 2 0 1
#> 7 19 P5913 1 1 1 1 2 0 1
#> 8 20 P5913 1 1 1 1 2 0 1
#> 9 22 P5913 1 1 1 1 2 0 1
#> 10 23 P5913 1 1 1 1 2 0 1
#> 11 25 P5913 1 1 1 1 2 0 1
#> 12 26 P5913 1 1 1 1 2 0 1
#> 13 27 P5913 1 1 1 1 2 0 1
#> 14 28 P5913 1 1 1 1 2 0 1
#> 15 29 P5913 1 1 1 1 2 0 1
#> 16 30 P5913 1 1 1 1 2 0 1
#> 17 33 P5913 1 1 1 1 2 0 1
#> 18 34 P5913 1 1 1 1 2 0 1
#> 19 36 P5913 1 1 1 1 2 0 1
#> 20 37 P5913 1 1 1 1 2 0 1
#> 21 38 P5913 1 1 1 1 2 0 1
#> 22 39 P5913 1 1 1 1 2 0 1
#> 23 40 P5913 1 1 1 1 2 0 1
#> 24 41 P5913 1 1 1 1 2 0 1
#> 25 42 P5913 1 1 1 1 2 0 1
#> 26 43 P5913 1 1 1 1 2 0 1
#> 27 45 P5913 1 1 1 1 2 0 1
#> 28 46 P5913 1 1 1 1 2 0 1
#> 29 48 P5913 1 1 1 1 2 0 1
#> 30 50 P5913 1 1 1 1 2 0 1
#> 31 51 P5913 1 1 1 1 2 0 1
#> 32 57 P5913 1 1 1 1 2 0 1
#> 33 58 P5913 1 1 1 1 2 0 1
#> 34 60 P5913 1 1 1 1 2 0 1
#> 35 67 P5913 1 1 1 1 2 0 1
#> 36 76 P5913 1 1 1 1 2 0 1
#> 37 79 P5913 1 1 1 1 2 0 1
#> 38 83 P5913 1 1 1 1 2 0 1
#> 39 84 P5913 1 1 1 1 2 0 1
#> 40 88 P5913 1 1 1 1 2 0 1
#> 41 89 P5913 1 1 1 1 2 0 1
#> 42 90 P5913 1 1 1 1 2 0 1
#> 43 95 P5913 1 1 1 1 2 0 1
#> 44 97 P5913 1 1 1 1 2 0 1
#> 45 102 P5913 1 1 1 1 2 0 1
#> 46 104 P5913 1 1 1 1 2 0 1
#> 47 106 P5913 1 1 1 1 2 0 1
#> emergency fup_obs age deprivation censor fup_censor clone fup_outcome
#> 0 109 79.68515 0.37 1 1 0 1
#> 0 109 79.68515 0.37 0 1 1 1
#> 0 109 79.68515 0.37 0 1 1 7
#> 0 109 79.68515 0.37 0 1 1 8
#> 0 109 79.68515 0.37 0 1 1 10
#> 0 109 79.68515 0.37 0 1 1 16
#> 0 109 79.68515 0.37 0 1 1 18
#> 0 109 79.68515 0.37 0 1 1 19
#> 0 109 79.68515 0.37 0 1 1 20
#> 0 109 79.68515 0.37 0 1 1 22
#> 0 109 79.68515 0.37 0 1 1 23
#> 0 109 79.68515 0.37 0 1 1 25
#> 0 109 79.68515 0.37 0 1 1 26
#> 0 109 79.68515 0.37 0 1 1 27
#> 0 109 79.68515 0.37 0 1 1 28
#> 0 109 79.68515 0.37 0 1 1 29
#> 0 109 79.68515 0.37 0 1 1 30
#> 0 109 79.68515 0.37 0 1 1 33
#> 0 109 79.68515 0.37 0 1 1 34
#> 0 109 79.68515 0.37 0 1 1 36
#> 0 109 79.68515 0.37 0 1 1 37
#> 0 109 79.68515 0.37 0 1 1 38
#> 0 109 79.68515 0.37 0 1 1 39
#> 0 109 79.68515 0.37 0 1 1 40
#> 0 109 79.68515 0.37 0 1 1 41
#> 0 109 79.68515 0.37 0 1 1 42
#> 0 109 79.68515 0.37 0 1 1 43
#> 0 109 79.68515 0.37 0 1 1 45
#> 0 109 79.68515 0.37 0 1 1 46
#> 0 109 79.68515 0.37 0 1 1 48
#> 0 109 79.68515 0.37 0 1 1 50
#> 0 109 79.68515 0.37 0 1 1 51
#> 0 109 79.68515 0.37 0 1 1 57
#> 0 109 79.68515 0.37 0 1 1 58
#> 0 109 79.68515 0.37 0 1 1 60
#> 0 109 79.68515 0.37 0 1 1 67
#> 0 109 79.68515 0.37 0 1 1 76
#> 0 109 79.68515 0.37 0 1 1 79
#> 0 109 79.68515 0.37 0 1 1 83
#> 0 109 79.68515 0.37 0 1 1 84
#> 0 109 79.68515 0.37 0 1 1 88
#> 0 109 79.68515 0.37 0 1 1 89
#> 0 109 79.68515 0.37 0 1 1 90
#> 0 109 79.68515 0.37 0 1 1 95
#> 0 109 79.68515 0.37 0 1 1 97
#> 0 109 79.68515 0.37 0 1 1 102
#> 0 109 79.68515 0.37 0 1 1 104
#> 0 109 79.68515 0.37 0 1 1 106
#> 0 109 79.68515 0.37 0 1 1 109
#> outcome t_stop t_event
#> 0 1 NA
#> 0 1 NA
#> 0 7 1
#> 0 8 7
#> 0 10 8
#> 0 16 10
#> 0 18 16
#> 0 19 18
#> 0 20 19
#> 0 22 20
#> 0 23 22
#> 0 25 23
#> 0 26 25
#> 0 27 26
#> 0 28 27
#> 0 29 28
#> 0 30 29
#> 0 33 30
#> 0 34 33
#> 0 36 34
#> 0 37 36
#> 0 38 37
#> 0 39 38
#> 0 40 39
#> 0 41 40
#> 0 42 41
#> 0 43 42
#> 0 45 43
#> 0 46 45
#> 0 48 46
#> 0 50 48
#> 0 51 50
#> 0 57 51
#> 0 58 57
#> 0 60 58
#> 0 67 60
#> 0 76 67
#> 0 79 76
#> 0 83 79
#> 0 84 83
#> 0 88 84
#> 0 89 88
#> 0 90 89
#> 0 95 90
#> 0 97 95
#> 0 102 97
#> 0 104 102
#> 0 106 104
#> 1 109 106Now we simply need to generate the weights. The
survivalCCW function generate_ccw() will do
this for us.
clones_long_weights <- generate_ccw(clones_long, predvars = c("age", "sex", "perf", "stage", "deprivation", "charlson", "emergency"))
head(clones_long_weights, row.names = FALSE)
#> time_id t_start id surgery timetosurgery death sex charlson perf stage
#> 1 0 0 P1074 0 NA 1 2 1 2 1
#> 2 1 1 P1074 0 NA 1 2 1 2 1
#> 3 2 7 P1074 0 NA 1 2 1 2 1
#> 4 3 8 P1074 0 NA 1 2 1 2 1
#> 5 4 10 P1074 0 NA 1 2 1 2 1
#> 6 5 16 P1074 0 NA 1 2 1 2 1
#> emergency fup_obs age deprivation censor fup_censor clone fup_outcome
#> 1 0 345 76.67899 0.33 0 182.62 0 1
#> 2 0 345 76.67899 0.33 0 182.62 0 7
#> 3 0 345 76.67899 0.33 0 182.62 0 8
#> 4 0 345 76.67899 0.33 0 182.62 0 10
#> 5 0 345 76.67899 0.33 0 182.62 0 16
#> 6 0 345 76.67899 0.33 0 182.62 0 18
#> outcome t_stop t_event lp t hazard p_uncens weight_cox
#> 1 0 1 NA -7.420857 NA 0.0000 1.0000000 1.000000
#> 2 0 7 1 -7.420857 1 162.5508 0.9072759 1.102201
#> 3 0 8 7 -7.420857 7 162.5508 0.9072759 1.102201
#> 4 0 10 8 -7.420857 8 171.7593 0.9022882 1.108293
#> 5 0 16 10 -7.420857 10 176.4477 0.8997594 1.111408
#> 6 0 18 16 -7.420857 16 181.2410 0.8971813 1.114602Let’s pick out a single patient and look at their data:
print(clones_long_weights[clones_long_weights$id == "P5913", ], row.names = FALSE)
#> time_id t_start id surgery timetosurgery death sex charlson perf stage
#> 0 0 P5913 1 1 1 1 2 0 1
#> 0 0 P5913 1 1 1 1 2 0 1
#> 1 1 P5913 1 1 1 1 2 0 1
#> 2 7 P5913 1 1 1 1 2 0 1
#> 3 8 P5913 1 1 1 1 2 0 1
#> 4 10 P5913 1 1 1 1 2 0 1
#> 5 16 P5913 1 1 1 1 2 0 1
#> 6 18 P5913 1 1 1 1 2 0 1
#> 7 19 P5913 1 1 1 1 2 0 1
#> 8 20 P5913 1 1 1 1 2 0 1
#> 9 22 P5913 1 1 1 1 2 0 1
#> 10 23 P5913 1 1 1 1 2 0 1
#> 11 25 P5913 1 1 1 1 2 0 1
#> 12 26 P5913 1 1 1 1 2 0 1
#> 13 27 P5913 1 1 1 1 2 0 1
#> 14 28 P5913 1 1 1 1 2 0 1
#> 15 29 P5913 1 1 1 1 2 0 1
#> 16 30 P5913 1 1 1 1 2 0 1
#> 17 33 P5913 1 1 1 1 2 0 1
#> 18 34 P5913 1 1 1 1 2 0 1
#> 19 36 P5913 1 1 1 1 2 0 1
#> 20 37 P5913 1 1 1 1 2 0 1
#> 21 38 P5913 1 1 1 1 2 0 1
#> 22 39 P5913 1 1 1 1 2 0 1
#> 23 40 P5913 1 1 1 1 2 0 1
#> 24 41 P5913 1 1 1 1 2 0 1
#> 25 42 P5913 1 1 1 1 2 0 1
#> 26 43 P5913 1 1 1 1 2 0 1
#> 27 45 P5913 1 1 1 1 2 0 1
#> 28 46 P5913 1 1 1 1 2 0 1
#> 29 48 P5913 1 1 1 1 2 0 1
#> 30 50 P5913 1 1 1 1 2 0 1
#> 31 51 P5913 1 1 1 1 2 0 1
#> 32 57 P5913 1 1 1 1 2 0 1
#> 33 58 P5913 1 1 1 1 2 0 1
#> 34 60 P5913 1 1 1 1 2 0 1
#> 35 67 P5913 1 1 1 1 2 0 1
#> 36 76 P5913 1 1 1 1 2 0 1
#> 37 79 P5913 1 1 1 1 2 0 1
#> 38 83 P5913 1 1 1 1 2 0 1
#> 39 84 P5913 1 1 1 1 2 0 1
#> 40 88 P5913 1 1 1 1 2 0 1
#> 41 89 P5913 1 1 1 1 2 0 1
#> 42 90 P5913 1 1 1 1 2 0 1
#> 43 95 P5913 1 1 1 1 2 0 1
#> 44 97 P5913 1 1 1 1 2 0 1
#> 45 102 P5913 1 1 1 1 2 0 1
#> 46 104 P5913 1 1 1 1 2 0 1
#> 47 106 P5913 1 1 1 1 2 0 1
#> emergency fup_obs age deprivation censor fup_censor clone fup_outcome
#> 0 109 79.68515 0.37 1 1 0 1
#> 0 109 79.68515 0.37 0 1 1 1
#> 0 109 79.68515 0.37 0 1 1 7
#> 0 109 79.68515 0.37 0 1 1 8
#> 0 109 79.68515 0.37 0 1 1 10
#> 0 109 79.68515 0.37 0 1 1 16
#> 0 109 79.68515 0.37 0 1 1 18
#> 0 109 79.68515 0.37 0 1 1 19
#> 0 109 79.68515 0.37 0 1 1 20
#> 0 109 79.68515 0.37 0 1 1 22
#> 0 109 79.68515 0.37 0 1 1 23
#> 0 109 79.68515 0.37 0 1 1 25
#> 0 109 79.68515 0.37 0 1 1 26
#> 0 109 79.68515 0.37 0 1 1 27
#> 0 109 79.68515 0.37 0 1 1 28
#> 0 109 79.68515 0.37 0 1 1 29
#> 0 109 79.68515 0.37 0 1 1 30
#> 0 109 79.68515 0.37 0 1 1 33
#> 0 109 79.68515 0.37 0 1 1 34
#> 0 109 79.68515 0.37 0 1 1 36
#> 0 109 79.68515 0.37 0 1 1 37
#> 0 109 79.68515 0.37 0 1 1 38
#> 0 109 79.68515 0.37 0 1 1 39
#> 0 109 79.68515 0.37 0 1 1 40
#> 0 109 79.68515 0.37 0 1 1 41
#> 0 109 79.68515 0.37 0 1 1 42
#> 0 109 79.68515 0.37 0 1 1 43
#> 0 109 79.68515 0.37 0 1 1 45
#> 0 109 79.68515 0.37 0 1 1 46
#> 0 109 79.68515 0.37 0 1 1 48
#> 0 109 79.68515 0.37 0 1 1 50
#> 0 109 79.68515 0.37 0 1 1 51
#> 0 109 79.68515 0.37 0 1 1 57
#> 0 109 79.68515 0.37 0 1 1 58
#> 0 109 79.68515 0.37 0 1 1 60
#> 0 109 79.68515 0.37 0 1 1 67
#> 0 109 79.68515 0.37 0 1 1 76
#> 0 109 79.68515 0.37 0 1 1 79
#> 0 109 79.68515 0.37 0 1 1 83
#> 0 109 79.68515 0.37 0 1 1 84
#> 0 109 79.68515 0.37 0 1 1 88
#> 0 109 79.68515 0.37 0 1 1 89
#> 0 109 79.68515 0.37 0 1 1 90
#> 0 109 79.68515 0.37 0 1 1 95
#> 0 109 79.68515 0.37 0 1 1 97
#> 0 109 79.68515 0.37 0 1 1 102
#> 0 109 79.68515 0.37 0 1 1 104
#> 0 109 79.68515 0.37 0 1 1 106
#> 0 109 79.68515 0.37 0 1 1 109
#> outcome t_stop t_event lp t hazard p_uncens weight_cox
#> 0 1 NA -5.592170 NA 0 1 1
#> 0 1 NA 5.653693 NA 0 1 1
#> 0 7 1 5.653693 1 0 1 1
#> 0 8 7 5.653693 7 0 1 1
#> 0 10 8 5.653693 8 0 1 1
#> 0 16 10 5.653693 10 0 1 1
#> 0 18 16 5.653693 16 0 1 1
#> 0 19 18 5.653693 18 0 1 1
#> 0 20 19 5.653693 19 0 1 1
#> 0 22 20 5.653693 20 0 1 1
#> 0 23 22 5.653693 22 0 1 1
#> 0 25 23 5.653693 23 0 1 1
#> 0 26 25 5.653693 25 0 1 1
#> 0 27 26 5.653693 26 0 1 1
#> 0 28 27 5.653693 27 0 1 1
#> 0 29 28 5.653693 28 0 1 1
#> 0 30 29 5.653693 29 0 1 1
#> 0 33 30 5.653693 30 0 1 1
#> 0 34 33 5.653693 33 0 1 1
#> 0 36 34 5.653693 34 0 1 1
#> 0 37 36 5.653693 36 0 1 1
#> 0 38 37 5.653693 37 0 1 1
#> 0 39 38 5.653693 38 0 1 1
#> 0 40 39 5.653693 39 0 1 1
#> 0 41 40 5.653693 40 0 1 1
#> 0 42 41 5.653693 41 0 1 1
#> 0 43 42 5.653693 42 0 1 1
#> 0 45 43 5.653693 43 0 1 1
#> 0 46 45 5.653693 45 0 1 1
#> 0 48 46 5.653693 46 0 1 1
#> 0 50 48 5.653693 48 0 1 1
#> 0 51 50 5.653693 50 0 1 1
#> 0 57 51 5.653693 51 0 1 1
#> 0 58 57 5.653693 57 0 1 1
#> 0 60 58 5.653693 58 0 1 1
#> 0 67 60 5.653693 60 0 1 1
#> 0 76 67 5.653693 67 0 1 1
#> 0 79 76 5.653693 76 0 1 1
#> 0 83 79 5.653693 79 0 1 1
#> 0 84 83 5.653693 83 0 1 1
#> 0 88 84 5.653693 84 0 1 1
#> 0 89 88 5.653693 88 0 1 1
#> 0 90 89 5.653693 89 0 1 1
#> 0 95 90 5.653693 90 0 1 1
#> 0 97 95 5.653693 95 0 1 1
#> 0 102 97 5.653693 97 0 1 1
#> 0 104 102 5.653693 102 0 1 1
#> 0 106 104 5.653693 104 0 1 1
#> 1 109 106 5.653693 106 0 1 1We now have everything we need to conduct a CCW analysis. For instance, we can pipe things together to evaluate the hazard ratio for surgery vs no surgery:
library(survival)
df <- toy_df |>
create_clones(
id = 'id',
event = 'death',
time_to_event = 'fup_obs',
exposure = 'surgery',
time_to_exposure = 'timetosurgery',
ced_window = 365.25/2
) |>
cast_clones_to_long() |>
generate_ccw(c('age', 'sex', 'perf', 'stage', 'deprivation', 'charlson', 'emergency'))
coxph(Surv(t_start, t_stop, outcome) ~ clone, data = df, weights = weight_cox)
#> Call:
#> coxph(formula = Surv(t_start, t_stop, outcome) ~ clone, data = df,
#> weights = weight_cox)
#>
#> coef exp(coef) se(coef) robust se z p
#> clone -0.4988 0.6073 0.2153 0.2794 -1.785 0.0742
#>
#> Likelihood ratio test=5.51 on 1 df, p=0.01894
#> n= 22772, number of events= 62Note that we used outcome and not death in
the coxph() model. Still, there is of course a problem with
this analysis, as the cloning process renders the variance invalid. The
simplest approach to addressing this is to bootstrap the variance. I
have not made a function to do this yet, but leave the below as an
example of how to do this.
library(boot)
#>
#> Attaching package: 'boot'
#> The following object is masked from 'package:survival':
#>
#> aml
boot_cox <- function(data, indices) {
# Make long data.frame with weights
ccw_df <- data[indices, ] |>
create_clones(
id = 'id',
event = 'death',
time_to_event = 'fup_obs',
exposure = 'surgery',
time_to_exposure = 'timetosurgery',
ced_window = 182.62
) |>
cast_clones_to_long() |>
generate_ccw(c('age', 'sex', 'perf', 'stage', 'deprivation', 'charlson', 'emergency'))
# Extract HR from CoxPH
cox_ccw <- coxph(Surv(t_start, t_stop, outcome) ~ clone, data = ccw_df, weights = weight_cox)
hr <- cox_ccw |>
coef() |>
exp()
out <- c("hr" = hr)
# Create survfit objects for each of treated and untreated
surv_1 <- survfit(Surv(t_start, t_stop, outcome) ~ 1L, data = ccw_df[ccw_df$clone == 1, ], weights = weight_cox)
surv_0 <- survfit(Surv(t_start, t_stop, outcome) ~ 1L, data = ccw_df[ccw_df$clone == 0, ], weights = weight_cox)
# RMST difference
rmst_1 <- surv_1 |>
summary(rmean = 365) |>
(\(summary) summary$table)() |>
(\(table) table["rmean"])()
rmst_0 <- surv_0 |>
summary(rmean = 365) |>
(\(summary) summary$table)() |>
(\(table) table["rmean"])()
rmst_diff <- rmst_1 - rmst_0
out <- c(out, "rmst_diff" = rmst_diff)
# 1-year survival difference
# Find the index of the time point closest to 1 year
index_1yr_1 <- which.min(abs(surv_1$time - 365))
index_1yr_0 <- which.min(abs(surv_0$time - 365))
# Get the 1-year survival probabilities
surv_1_1yr <- surv_1$surv[index_1yr_1]
surv_0_1yr <- surv_0$surv[index_1yr_0]
surv_diff_1yr <- surv_1_1yr - surv_0_1yr
out <- c(out, "surv_diff_1yr" = surv_diff_1yr)
}
boot_out <- boot(data = toy_df, statistic = boot_cox, R = 10)
boot.ci(boot_out, type = "norm", index = 1)
#> BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
#> Based on 10 bootstrap replicates
#>
#> CALL :
#> boot.ci(boot.out = boot_out, type = "norm", index = 1)
#>
#> Intervals :
#> Level Normal
#> 95% ( 0.1945, 0.9168 )
#> Calculations and Intervals on Original Scale
boot.ci(boot_out, type = "norm", index = 2)
#> BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
#> Based on 10 bootstrap replicates
#>
#> CALL :
#> boot.ci(boot.out = boot_out, type = "norm", index = 2)
#>
#> Intervals :
#> Level Normal
#> 95% (-10.629, 19.976 )
#> Calculations and Intervals on Original Scale
boot.ci(boot_out, type = "norm", index = 3)
#> BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
#> Based on 10 bootstrap replicates
#>
#> CALL :
#> boot.ci(boot.out = boot_out, type = "norm", index = 3)
#>
#> Intervals :
#> Level Normal
#> 95% ( 0.0066, 0.2580 )
#> Calculations and Intervals on Original ScaleThe extension of the package functionality to issues of competing risks is trivial. We keep our outcome variable as an integer but add additional values.
toy_df$competing <- sample(0:2, nrow(toy_df), replace = TRUE)
head(toy_df)
#> id surgery timetosurgery death sex charlson perf stage emergency fup_obs
#> 1 P5958 0 NA 0 2 1 2 2 1 365.24
#> 2 P3075 0 NA 0 2 1 0 2 0 365.24
#> 3 P1410 0 NA 0 1 1 2 1 0 365.24
#> 4 P4957 0 NA 0 2 0 1 1 0 365.24
#> 5 P7340 0 NA 1 1 0 2 2 1 123.00
#> 6 P1917 0 NA 0 2 0 0 1 0 365.24
#> age deprivation competing
#> 1 79.57290 0.42 2
#> 2 82.91855 0.07 0
#> 3 84.03011 0.22 0
#> 4 78.34634 0.12 0
#> 5 83.63860 0.40 0
#> 6 82.26147 0.03 1We can then conduct the same analysis as before, but with the
competing variable as the outcome. The bootstrapped
analysis is not shown, but you can see that hte same approach can be
used.
compete_ccw <- toy_df |>
create_clones(
id = 'id',
event = 'competing',
time_to_event = 'fup_obs',
exposure = 'surgery',
time_to_exposure = 'timetosurgery',
ced_window = 365.25/2
) |>
cast_clones_to_long() |>
generate_ccw(c('age', 'sex', 'perf', 'stage', 'deprivation', 'charlson', 'emergency'))
table(compete_ccw$outcome)
#>
#> 0 1 2
#> 22625 76 71
head(compete_ccw)
#> time_id t_start id surgery timetosurgery death sex charlson perf stage
#> 1 0 0 P1074 0 NA 1 2 1 2 1
#> 2 1 1 P1074 0 NA 1 2 1 2 1
#> 3 2 7 P1074 0 NA 1 2 1 2 1
#> 4 3 8 P1074 0 NA 1 2 1 2 1
#> 5 4 10 P1074 0 NA 1 2 1 2 1
#> 6 5 16 P1074 0 NA 1 2 1 2 1
#> emergency fup_obs age deprivation competing censor fup_censor clone
#> 1 0 345 76.67899 0.33 2 0 182.625 0
#> 2 0 345 76.67899 0.33 2 0 182.625 0
#> 3 0 345 76.67899 0.33 2 0 182.625 0
#> 4 0 345 76.67899 0.33 2 0 182.625 0
#> 5 0 345 76.67899 0.33 2 0 182.625 0
#> 6 0 345 76.67899 0.33 2 0 182.625 0
#> fup_outcome outcome t_stop t_event lp t hazard p_uncens weight_cox
#> 1 1 0 1 NA -7.420857 NA 0.0000 1.0000000 1.000000
#> 2 7 0 7 1 -7.420857 1 162.5508 0.9072759 1.102201
#> 3 8 0 8 7 -7.420857 7 162.5508 0.9072759 1.102201
#> 4 10 0 10 8 -7.420857 8 171.7593 0.9022882 1.108293
#> 5 16 0 16 10 -7.420857 10 176.4477 0.8997594 1.111408
#> 6 18 0 18 16 -7.420857 16 181.2410 0.8971813 1.114602