1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
|
R Under development (unstable) (2024-06-14 r86747) -- "Unsuffered Consequences"
Copyright (C) 2024 The R Foundation for Statistical Computing
Platform: aarch64-unknown-linux-gnu
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> library(survival)
> aeq <- function(x, y, ...) all.equal(as.vector(x), as.vector(y), ...)
>
> # Check that a multi-state model, correctly set up, gives the same
> # solution as a time-dependent covariate.
> # This is a stronger test than mstrata: there the covariate which was mapped
> # into a state was constant, here it is time-dependent.
> #
> # First build the TD data set from pbcseq, with a categorical bilirubin
> pbc1 <- pbcseq
> pbc1$bili4 <- cut(pbc1$bili, c(0,1, 2,4, 100),
+ c("normal", "1-2x", "2-4x", ">4"))
> ptemp <- subset(pbc1, !duplicated(id)) # first row of each
>
> pbc2 <- tmerge(ptemp[, c("id", "age", "sex")], ptemp, id,
+ death= event(futime, status==2))
>
> pbc2 <- tmerge(pbc2, pbc1, id=id, bili = tdc(day, bili),
+ bili4 = tdc(day, bili4), bstat = event(day, as.numeric(bili4)))
> btemp <- with(pbc2, ifelse(death, 5, bstat))
>
> # a row with the same starting and ending bili4 level is not an event
> b2 <- ifelse(((as.numeric(pbc2$bili4)) == btemp), 0, btemp)
> pbc2$bstat <- factor(b2, 0:5,
+ c("censor", "normal", "1-2x", "2-4x", ">4", "death"))
> check1 <- survcheck(Surv(tstart, tstop, bstat) ~ 1, istate= bili4,
+ id = id, data=pbc2)
> check1$transitions
to
from normal 1-2x 2-4x >4 death (censored)
normal 0 81 10 3 9 77
1-2x 61 0 68 15 9 36
2-4x 2 33 0 94 12 24
>4 1 3 28 0 110 35
death 0 0 0 0 0 0
> all.equal(as.character(pbc2$bili4), as.character(check1$istate))
[1] TRUE
> # the above verifies that I created the data set correctly
>
> # Standard coxph fit with a time dependent bili4 variable.
> fit1 <- coxph(Surv(tstart, tstop, death) ~ age + bili4, pbc2, ties='breslow')
>
> # An additive multi-state fit, where bili4 is a state
> # The three forms below should all give identical models
> fit2 <- coxph(list(Surv(tstart, tstop, bstat) ~ 1,
+ c(1:4):5 ~ age / common + shared), id= id, istate=bili4,
+ data=pbc2)
> fit2b <- coxph(list(Surv(tstart, tstop, bstat) ~ 1,
+ 1:5 + 2:5 + 3:5 + 4:5 ~ age / common + shared),
+ id= id, istate=bili4, data=pbc2)
> fit2c <- coxph(list(Surv(tstart, tstop, bstat) ~ 1,
+ 0:5 ~ age / common + shared),
+ id= id, istate=bili4, data=pbc2)
>
> # Make sure the names are correct and the coefficients match
> aeq(coef(fit1), coef(fit2))
[1] TRUE
> aeq(names(coef(fit2)), c("age", "ph(2:5/1:5)", "ph(3:5/1:5)", "ph(4:5/1:5)"))
[1] TRUE
> all.equal(coef(fit2), coef(fit2b))
[1] TRUE
> all.equal(coef(fit2), coef(fit2c))
[1] TRUE
>
> # Now a model with a separate age effect for each bilirubin group
> fit3 <- coxph(Surv(tstart, tstop, death) ~ age*bili4, pbc2, ties='breslow')
> fit3b <- coxph(Surv(tstart, tstop, death) ~ bili4/age, pbc2, ties='breslow')
> fit4 <- coxph(list(Surv(tstart, tstop, bstat) ~ 1,
+ c(1:4):5 ~ age / shared), id= id, istate=bili4,
+ data=pbc2)
> all.equal(fit3$loglik, fit3b$loglik)
[1] TRUE
> all.equal(fit3$loglik, fit4$loglik)
[1] TRUE
>
> # The coefficients are quite different due to different codings for dummy vars
> # Unpack the interaction, first 4 coefs will be the age effect within each
> # bilirubin group
> temp <- c(coef(fit3)[1] + c(0, coef(fit3)[5:7]), coef(fit3)[2:4])
> names(temp)[1:4] <- c("age1", "age2", "age3", "age4")
> aeq(temp, coef(fit3b)[c(4:7, 1:3)])
[1] TRUE
> aeq(temp, coef(fit4))
[1] TRUE
>
> # Third, a model with separate baseline hazards for each bili group
> fit5 <- coxph(Surv(tstart, tstop, death) ~ strata(bili4)/age, pbc2,
+ cluster=id, ties='breslow')
> fit6 <- coxph(list(Surv(tstart, tstop, bstat) ~ 1, 0:5 ~ age),
+ id=id, istate=bili4, pbc2)
> aeq(coef(fit5), coef(fit6))
[1] TRUE
> aeq(fit5$var, fit6$var)
[1] TRUE
> aeq(fit5$naive.var, fit6$naive.var)
[1] TRUE
>
> proc.time()
user system elapsed
0.536 0.024 0.557
|
