Extensions for the cmprsk package.
See intro vignette
To install:
# install.packages('devtools')
devtools::install_github('raredd/cmprsk2', build_vignettes = TRUE)## model deaths with ltx and withdraw as competing events
cr1 <- crr2(Surv(futime, event(censored) == death) ~ age + sex + abo,
data = transplant)
## include an all-cause death overall survival model
cr2 <- crr2(Surv(futime, event(censored) == death) ~ age + sex + abo,
data = transplant,
cox = Surv(futime, event == 'death') ~ age + sex + abo)summary(cr1)## $`CRR: death`
## HR L95 U95 p
## age 1.02 1.00 1.04 0.09
## sexf 0.65 0.40 1.07 0.09
## aboB 1.48 0.70 3.12 0.31
## aboAB 1.14 0.34 3.87 0.83
## aboO 1.48 0.85 2.56 0.17
##
## $`CRR: ltx`
## HR L95 U95 p
## age 0.99 0.99 1.00 0.15
## sexf 1.09 0.93 1.29 0.28
## aboB 0.70 0.54 0.91 0.01
## aboAB 0.91 0.59 1.39 0.65
## aboO 0.59 0.49 0.70 0.00
##
## $`CRR: withdraw`
## HR L95 U95 p
## age 0.98 0.95 1.02 0.28
## sexf 1.42 0.75 2.68 0.28
## aboB 2.35 0.82 6.74 0.11
## aboAB 3.04 0.81 11.43 0.10
## aboO 2.37 1.04 5.41 0.04
library('htmlTable')
summary(
cr2,
html = TRUE, n = TRUE, ref = TRUE,
htmlArgs = list(
caption = 'CRR models.',
rgroup = c('Age', 'Sex', 'Blood type'),
rnames = c('+1 year change', 'Female', 'B', 'AB', 'O'),
css.cell = 'white-space: nowrap; padding: 0px 5px 0px;'
)
)| CRR models. | |||||||||||||||||
| Cox PH | CRR: death | CRR: ltx | CRR: withdraw | ||||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Total n = 797 (%) |
Events n = 66 (8) |
HR (95% CI) | p | Events n = 66 (8) |
HR (95% CI) | p | Events n = 618 (78) |
HR (95% CI) | p | Events n = 37 (5) |
HR (95% CI) | p | |||||
| Age | |||||||||||||||||
| +1 year change | 797 (100) | 66 (100) | 1.02 (0.99, 1.04) | 0.16 | 66 (100) | 1.02 (1.00, 1.04) | 0.090 | 618 (100) | 0.99 (0.99, 1.00) | 0.15 | 37 (100) | 0.98 (0.95, 1.02) | 0.28 | ||||
| Sex | |||||||||||||||||
| Female | 438 (55) | 42 (64) | 42 (64) | 332 (54) | 17 (46) | ||||||||||||
| B | 359 (45) | 24 (36) | 0.69 (0.42, 1.14) | 0.15 | 24 (36) | 0.65 (0.40, 1.07) | 0.093 | 286 (46) | 1.09 (0.93, 1.29) | 0.28 | 20 (54) | 1.42 (0.75, 2.68) | 0.28 | ||||
| Blood type | |||||||||||||||||
| AB | 317 (40) | 21 (32) | 21 (32) | 261 (42) | 8 (22) | ||||||||||||
| O | 102 (13) | 10 (15) | 1.15 (0.54, 2.45) | 0.72 | 10 (15) | 1.48 (0.70, 3.12) | 0.31 | 77 (12) | 0.70 (0.54, 0.91) | 0.008 | 6 (16) | 2.35 (0.82, 6.74) | 0.11 | ||||
| NA | 40 (5) | 3 (5) | 1.18 (0.35, 3.96) | 0.79 | 3 (5) | 1.14 (0.34, 3.87) | 0.83 | 32 (5) | 0.91 (0.59, 1.39) | 0.65 | 3 (8) | 3.04 (0.81, 11.43) | 0.10 | ||||
| NA | 338 (42) | 32 (48) | 0.99 (0.57, 1.73) | 0.98 | 32 (48) | 1.48 (0.85, 2.56) | 0.17 | 248 (40) | 0.59 (0.49, 0.70) | < 0.001 | 20 (54) | 2.37 (1.04, 5.41) | 0.040 | ||||
## can use same formula as crr2
ci1 <- cuminc2(Surv(futime, event(censored) == death) ~ abo,
data = transplant)
## but event indicator is not required
ci2 <- cuminc2(Surv(futime, event(censored)) ~ sex,
data = transplant)summary(ci1)## $est
## 0 500 1000 1500 2000
## A death 0.00000 0.06191 0.06191 NA NA
## B death 0.00000 0.08799 0.08799 0.12717 0.12717
## AB death 0.00000 0.07317 NA NA NA
## O death 0.00289 0.09106 0.09598 0.09598 NA
## A ltx 0.00000 0.81022 0.84174 NA NA
## B ltx 0.00000 0.71905 0.77478 0.77478 0.77478
## AB ltx 0.02439 0.80488 NA NA NA
## O ltx 0.00000 0.69130 0.76563 0.78129 NA
## A withdraw 0.00000 0.02481 0.02481 NA NA
## B withdraw 0.00000 0.05886 0.05886 0.05886 0.05886
## AB withdraw 0.00000 0.07317 NA NA NA
## O withdraw 0.00000 0.04413 0.06470 0.06470 NA
##
## $var
## 0 500 1000 1500 2000
## A death 0.00000 0.00018 0.00018 NA NA
## B death 0.00000 0.00081 0.00081 0.00259 0.00259
## AB death 0.00000 0.00192 NA NA NA
## O death 0.00001 0.00024 0.00027 0.00027 NA
## A ltx 0.00000 0.00048 0.00046 NA NA
## B ltx 0.00000 0.00202 0.00233 0.00233 0.00233
## AB ltx 0.00059 0.00426 NA NA NA
## O ltx 0.00000 0.00064 0.00058 0.00062 NA
## A withdraw 0.00000 0.00008 0.00008 NA NA
## B withdraw 0.00000 0.00056 0.00056 0.00056 0.00056
## AB withdraw 0.00000 0.00177 NA NA NA
## O withdraw 0.00000 0.00012 0.00020 0.00020 NA
##
## $events
## 0 500 1000 1500 2000
## A death 0 20 20 21 21
## B death 0 9 9 10 10
## AB death 0 3 3 3 3
## O death 1 31 32 32 32
## A ltx 0 262 269 269 269
## B ltx 0 74 77 77 77
## AB ltx 1 33 33 33 33
## O ltx 0 234 254 256 256
## A withdraw 0 8 8 8 8
## B withdraw 0 6 6 6 6
## AB withdraw 0 3 3 3 3
## O withdraw 0 15 20 20 20
##
## $total_events
## A death B death AB death O death A ltx B ltx AB ltx O ltx
## 21 10 3 32 269 78 33 256
## A withdraw B withdraw AB withdraw O withdraw
## 8 6 3 20
##
## $total_groups
## A B AB O
## 325 103 41 346
##
## $total_atrisk
## 0 500 1000 1500 2000
## 811 94 20 2 1
summary(ci1, times = 0:10 * 100)$est## 0 100 200 300 400 500 600 700 800 900 1000
## A death 0.00000 0.03390 0.04317 0.05562 0.06191 0.06191 0.06191 0.06191 0.06191 0.06191 0.06191
## B death 0.00000 0.04854 0.06796 0.07767 0.08799 0.08799 0.08799 0.08799 0.08799 0.08799 0.08799
## AB death 0.00000 0.04878 0.04878 0.04878 0.04878 0.07317 NA NA NA NA NA
## O death 0.00289 0.05250 0.07597 0.08185 0.08779 0.09106 0.09106 0.09106 0.09106 0.09106 0.09598
## A ltx 0.00000 0.48743 0.72831 0.77824 0.80346 0.81022 0.83623 0.84174 0.84174 0.84174 0.84174
## B ltx 0.00000 0.33010 0.50485 0.69903 0.71905 0.71905 0.74866 0.74866 0.74866 0.77478 0.77478
## AB ltx 0.02439 0.56098 0.60976 0.73171 0.80488 0.80488 NA NA NA NA NA
## O ltx 0.00000 0.23354 0.43015 0.54777 0.64341 0.69130 0.73601 0.75159 0.75580 0.76071 0.76563
## A withdraw 0.00000 0.00926 0.01853 0.02167 0.02481 0.02481 0.02481 0.02481 0.02481 0.02481 0.02481
## B withdraw 0.00000 0.02913 0.03883 0.04854 0.05886 0.05886 0.05886 0.05886 0.05886 0.05886 0.05886
## AB withdraw 0.00000 0.02439 0.07317 0.07317 0.07317 0.07317 NA NA NA NA NA
## O withdraw 0.00000 0.01167 0.02927 0.03810 0.04413 0.04413 0.04413 0.06470 0.06470 0.06470 0.06470
par(mfrow = c(2, 2))
# ciplot(ci2)
plot(ci2, add = TRUE) ## equivalently
plot(ci2, split = 'event', add = TRUE, wh.events = 'est')
## convenience wrapper
par(mfrow = c(2, 2))
ciplot_by(
rhs = 'sex', time = 'futime', event = 'event',
data = transplant, by = 'abo', xlim = c(0, 1500),
events = FALSE, single = FALSE, events.total = 2100
)
## pairwise gray tests
cuminc_pairs(ci1)$p.value## $death
## A B AB O
## A NA 1.000 1.000 1
## B 0.391 NA 1.000 1
## AB 0.788 0.776 NA 1
## O 0.201 0.877 0.712 NA
##
## $ltx
## A B AB O
## A NA 0.035 0.805 0.000
## B 0.007 NA 0.442 0.442
## AB 0.805 0.162 NA 0.051
## O 0.000 0.147 0.013 NA
##
## $withdraw
## A B AB O
## A NA 0.446 0.446 0.253
## B 0.097 NA 1.000 1.000
## AB 0.089 0.753 NA 1.000
## O 0.042 0.885 0.402 NA
timepoints2(
ci2, html = TRUE,
htmlArgs = list(
caption = 'cuminc estimates at specific time points (<code>cuminc::timepoints</code>).'
)
)
cuminc estimates at specific time points (cuminc::timepoints). | |||||
| 0 | 500 | 1000 | 1500 | 2000 | |
|---|---|---|---|---|---|
| m death | 0.002 | 0.091 | 0.097 | 0.104 | 0.104 |
| f death | 0.000 | 0.063 | 0.063 | - | - |
| m ltx | 0.000 | 0.738 | 0.786 | 0.800 | 0.800 |
| f ltx | 0.003 | 0.761 | 0.822 | - | - |
| m withdraw | 0.000 | 0.030 | 0.045 | 0.045 | 0.045 |
| f withdraw | 0.000 | 0.052 | 0.056 | - | - |
AIC(cr1$`CRR: death`)## AIC
## 867.6161
sapply(cr1, BIC)## CRR: death.BIC CRR: ltx.BIC CRR: withdraw.BIC
## 878.5644 7494.9446 496.6589
logLik(cr1$`CRR: death`)## 'log Lik.' -428.8081 (df=5)
deviance(cr1$`CRR: death`)## Call:
## crr(transplant[, "futime"], transplant[, "event"], cov1 = model.matrix(~age +
## sex + abo, transplant)[, -1L, drop = FALSE], cencode = "censored",
## failcode = "death", variance = TRUE, cengroup = rep(1L, nrow(transplant)),
## gtol = 1e-06, maxiter = 10, init = c(0L, 0L, 0L, 0L, 0L))
##
## Deviance = 7.15 on 5 df, 0.20984
crrFits(cr1$`CRR: death`)## Model selection table
##
## 0: Null Model
##
## 1: Model 1 call:
## crr(transplant[, "futime"], transplant[, "event"], cov1 = model.matrix(~age +
## sex + abo, transplant)[, -1L, drop = FALSE], cencode = "censored",
## failcode = "death", variance = TRUE, cengroup = rep(1L, nrow(transplant)),
## gtol = 1e-06, maxiter = 10, init = c(0L, 0L, 0L, 0L, 0L))
##
##
## n loglik df k -2logLik -2logLik diff AIC AIC diff BIC BIC diff
## 0 797 -432.38 0 1 864.76 7.1484 864.76 0.0000 864.76 0.0
## 1 797 -428.81 5 6 857.62 0.0000 867.62 2.8516 878.56 13.8
crrwald.test(cr1$`CRR: death`)## chi2 df P
## Overall 8.11151956 5 0.15019572
## age 2.86948310 1 0.09027386
## sexf 2.82220518 1 0.09296860
## aboB 1.04828004 1 0.30590354
## aboAB 0.04546327 1 0.83115446
## aboO 1.91025511 1 0.16693492