Survival Analysis in R Lab 2

1. Hands-on excerise fitting multistate Markov model for panel data

Panel data

  • State only observed at a finite number of times
  • The only data information that one needs is the number of all observed transitions between states and the total time at risk in each state, aggregated over all individuals
  • (intermittently-observed) don’t need to know the state between these times, e.g., managing chronic conditions at clinical visits
  • attractive modelling approach when handling large data sets.
  • msm package is developed for multistate modelling of panel data where the focus is on estimating transition probabilities between states.

1.1 Revisit Heart transplant monitoring data

  • A series of approximately yearly angiographic examinations of heart transplant recipients. The state at each time is a grade of cardiac allograft vasculopathy (cav), a deterioration of the arterial walls.

  • A data frame containing 2846 rows.

  • There are 622 patients, the rows are grouped by patient number and ordered by years after transplant, with each row representing an examination and containing additional covariates.

  • PTNUM (numeric) Patient identification number

  • age (numeric) Recipient age at examination (years)

  • years (numeric) Examination time (years after transplant)

  • dage (numeric) Age of heart donor (years)

  • sex (numeric) sex (0=male, 1=female)

  • pdiag (factor) Primary diagnosis (reason for transplant)

    • IHD=ischaemic heart disease,
    • IDC=idiopathic dilated cardiomyopathy.

  • cumrej (numeric) Cumulative number of acute rejection episodes

  • state (numeric) State at the examination.

    • State 1 represents no CAV,
    • state 2 is mild/moderate CAV
    • and state 3 is severe CAV.
    • State 4 indicates death.

  • firstobs (numeric)

    • 0 = record represents an angiogram or date of death.
    • 1 = record represents transplant (patient’s first observation)

  • statemax (numeric) Maximum observed state so far for this patient

library(msm)
library(tidyverse)
library(gtsummary)
library(DT)
# View(cav)
# str(cav)

# look at the data;
datatable(cav,
          rownames = FALSE) %>%
  formatRound(columns=c('age','years'), digits=1)
# creating baseline table
cav_baseline <- cav %>% 
  arrange(PTNUM, years) %>%
  group_by(PTNUM) %>%
  slice_head(n=1) #keeping the first visit data for each subject;

subset(cav_baseline, select = c(age, sex, dage, pdiag)) %>% 
  tbl_summary(
  missing_text = "(Missing)")
Characteristic N = 6221
age 50 (42, 55)
sex 87 (14%)
dage 29 (20, 40)
pdiag
    CVCM 19 (3.1%)
    Hyper 3 (0.5%)
    IDC 270 (44%)
    IHD 313 (51%)
    Other 4 (0.7%)
    Restr 5 (0.8%)
    (Missing) 8
1 Median (Q1, Q3); n (%)

1.2 Multistate models

  • Study objectives:
    • Estimate transition intensity (TIM) and transition probability/risk functions (TPM)
    • identify association between risk factors and CAV progression and survival
  • Summarize observed transitions
statetable.msm(state, PTNUM, cav)
    to
from    1    2    3    4
   1 1367  204   44  148
   2   46  134   54   48
   3    4   13  107   55

  • Define allowed transition as a matrix of binary indicators
    • in this example, we do not allow transition from state 1 to 3 directly
    • we do observe 44 intervals with transition 1 to 3, we assume there is an unobserved passing transition 1->2->3.
  • The transition model following the hypothesized process

\[ Q = \left[\begin{array}{cccc} -(q_{12} + q_{14}) & q_{12} & 0 & q_{14} \\ q_{21} & -(q_{21}+q_{23}+q_{24}) & q_{23} & q_{24} \\ 0 & q_{32} & -(q_{32}+q_{34}) & q_{34} \\ 0 & 0 & 0 & 0 \end{array}\right] \]

#setting initial value;
Q <- rbind (
c(0, 0.25, 0, 0.25), #i.e., from state 1 we can transit to state 2 and 4;
c(0.166, 0, 0.166, 0.166),
c(0, 0.25, 0, 0.25),
c(0, 0, 0, 0))
Q
      [,1] [,2]  [,3]  [,4]
[1,] 0.000 0.25 0.000 0.250
[2,] 0.166 0.00 0.166 0.166
[3,] 0.000 0.25 0.000 0.250
[4,] 0.000 0.00 0.000 0.000
Q.crude <- crudeinits.msm(state ~ years, PTNUM, data=cav, qmatrix=Q)
Q.crude
           [,1]        [,2]       [,3]       [,4]
[1,] -0.1173149  0.06798932  0.0000000 0.04932559
[2,]  0.1168179 -0.37584883  0.1371340 0.12189692
[3,]  0.0000000  0.04908401 -0.2567471 0.20766310
[4,]  0.0000000  0.00000000  0.0000000 0.00000000

transition matrix from state 1 can be interpreted as, we expect on average 1/(0.25+0.25) = 2 years in state 1 and we assume equal chance from state 1 to transition to state 2 and 4.

  • A more practically meaningful parameterisation of a continuous-time Markov model with transition intensities \(q_{rs}\) is in terms of the mean sojourn times \(-q_{rr}\) in each state \(r\) and the probabilities that the next move of the process when in state \(r\) is to state \(s\), \(1-q_{rs}/q_{rr}\).

1.3 Simple naive model

  • Model1: simple multistate model assuming no information on the exact time of state entry

1.2.1 TIM

# Model 1
heart.msm1 <- msm(state ~ years, subject=PTNUM, data=cav, qmatrix=Q.crude)
heart.msm1

Call:
msm(formula = state ~ years, subject = PTNUM, data = cav, qmatrix = Q.crude)

Maximum likelihood estimates

Transition intensities
                  Baseline                    
State 1 - State 1 -0.17471 (-0.19506,-0.15649)
State 1 - State 2  0.12607 ( 0.10968, 0.14491)
State 1 - State 4  0.04864 ( 0.04008, 0.05903)
State 2 - State 1  0.23789 ( 0.17790, 0.31810)
State 2 - State 2 -0.61885 (-0.72334,-0.52945)
State 2 - State 3  0.30508 ( 0.24457, 0.38056)
State 2 - State 4  0.07588 ( 0.04289, 0.13427)
State 3 - State 2  0.15067 ( 0.09222, 0.24616)
State 3 - State 3 -0.48506 (-0.61548,-0.38227)
State 3 - State 4  0.33439 ( 0.25533, 0.43794)

-2 * log-likelihood:  3986.087 
# Model estimated transition intensity matrix;
qmatrix <- qmatrix.msm(heart.msm1)
qmatrix
        State 1                      State 2                     
State 1 -0.17471 (-0.19506,-0.15649)  0.12607 ( 0.10968, 0.14491)
State 2  0.23789 ( 0.17790, 0.31810) -0.61885 (-0.72334,-0.52945)
State 3 0                             0.15067 ( 0.09222, 0.24616)
State 4 0                            0                           
        State 3                      State 4                     
State 1 0                             0.04864 ( 0.04008, 0.05903)
State 2  0.30508 ( 0.24457, 0.38056)  0.07588 ( 0.04289, 0.13427)
State 3 -0.48506 (-0.61548,-0.38227)  0.33439 ( 0.25533, 0.43794)
State 4 0                            0
  • we see patients are three times as likely to transition from no CAV to mild CAV then death.
  • Once in mild CAV state, patients are more likely to progress to sever CAV than to no CAV.
  • Once in the severe state, death is more likely to occur.

1.2.2 Ratio of transition intensities

qratio.msm(heart.msm1, ind1=c(2,3), ind2=c(2,1))
 estimate        se         L         U 
1.2824682 0.2389993 0.8900567 1.8478877
  • At state 2 (mild CAV), progress to state 3 (severe CAV) is 1.3 times more likely than recovery (no CAV).

1.2.3 Estimated mean time in each state

sojourn.msm(heart.msm1)
        estimates        SE        L        U
State 1  5.723630 0.3217075 5.126585 6.390207
State 2  1.615900 0.1286282 1.382475 1.888736
State 3  2.061599 0.2504859 1.624735 2.615929

The average amount of time a patient spends in the mild CAV state is 1.62 years (95% CI 1.38 to 1.89).

1.2.4 Transition probability matrix

# estimating the transition probability matrix at 5 years;
pmatrix.msm(heart.msm1, t=5, 
            ci = "norm")
        State 1                   State 2                  
State 1 0.51168 (0.47457,0.54511) 0.13235 (0.11358,0.15067)
State 2 0.24973 (0.19794,0.29877) 0.13272 (0.10549,0.16636)
State 3 0.06806 (0.04069,0.10487) 0.06938 (0.04308,0.10582)
State 4 0                         0                        
        State 3                   State 4                  
State 1 0.07304 (0.05832,0.08828) 0.28293 (0.25427,0.31731)
State 2 0.14048 (0.10674,0.17781) 0.47707 (0.41834,0.54033)
State 3 0.13737 (0.08676,0.19529) 0.72519 (0.63406,0.80714)
State 4 0                         1.00000 (1.00000,1.00000)
            # ci="bootstrap", cores = 5, B=1000)
  • For patients in state 1, at 5 years there is about 51% chances of staying in state 1.
  • For patients in state 3, at 5 years there is about 73% chances of being dead.

1.4 Model with covariates and exact death time

  • Assume exact death time is recorded, but not other state entry times (same covariates on all transition rates)

1.2.1 Fitting covariates

# Model 2
heart.msm2 <- msm(state ~ years, 
                  PTNUM, 
                  data=cav,
                  qmatrix=Q.crude, 
                  deathexact = 4,
                  covariates = ~ sex,
                  control=list(fnscale=5000))
# deathexact: Vector of indices of absorbing states whose time of entry is known exactly, but the individual is assumed to be in an unknown transient state ("alive") at the previous instant. This is the usual situation for times of death in chronic disease monitoring data. For example, if you specify deathexact = c(4) then states 4 is assumed to be exactly-observed death states;
# It is often worthwhile to normalize the optimisation using control=list(fnscale = a), where a is the a number of the order of magnitude of the -2 log likelihood;

heart.msm3 <- msm(state ~ years, 
                  PTNUM, 
                  data=cav,
                  qmatrix=Q.crude, 
                  deathexact = 4,
                  covariates = ~ sex + age,
                  control=list(fnscale=5000))

# example to compare covariates models;
logliks <- c(logLik(heart.msm1), logLik(heart.msm2), logLik(heart.msm3))
AICs <- AIC(heart.msm1, heart.msm2, heart.msm3)
cbind(logliks, AICs)
             logliks df      AIC
heart.msm1 -1993.044  7 4000.087
heart.msm2 -1977.649 14 3983.299
heart.msm3 -1960.507 21 3963.013
  • Model 3 is has the best fit.

1.2.2 Hazard Ratio of transition risk

hazard.msm(heart.msm3)
$sex
                         HR          L         U
State 1 - State 2 0.6599402 0.39500958  1.102558
State 1 - State 4 1.2583468 0.65425689  2.420206
State 2 - State 1 1.3602422 0.52150948  3.547891
State 2 - State 3 1.0581743 0.50005884  2.239202
State 2 - State 4 0.7770243 0.02933949 20.578640
State 3 - State 2 0.8375957 0.09833204  7.134668
State 3 - State 4 2.5174948 1.23111129  5.148016

$age
                         HR         L        U
State 1 - State 2 1.0117661 0.9985370 1.025170
State 1 - State 4 1.0618814 1.0302436 1.094491
State 2 - State 1 1.0219367 0.9902139 1.054676
State 2 - State 3 0.9786256 0.9557051 1.002096
State 2 - State 4 1.0341056 0.8700592 1.229082
State 3 - State 2 0.9543338 0.9073848 1.003712
State 3 - State 4 0.9916257 0.9659416 1.017993
  • For transition 3 to 4, female is at higher risk.
  • For patients in state 1, 1 unit increase in age is associated with 1.06 times of higher risk of death.

1.2.3 Estimating specific transition probabiliy by covariate

# estimating the transition probability matrix at 5 years;
pmatrix.msm(heart.msm1, t=5, 
            covariates = list(sex = 1, age = 65),
            ci = "norm") 
        State 1                   State 2                  
State 1 0.51168 (0.47770,0.54511) 0.13235 (0.11314,0.15062)
State 2 0.24973 (0.19882,0.30441) 0.13272 (0.10138,0.16333)
State 3 0.06806 (0.04259,0.10592) 0.06938 (0.04169,0.10172)
State 4 0                         0                        
        State 3                   State 4                  
State 1 0.07304 (0.05836,0.08738) 0.28293 (0.25618,0.31849)
State 2 0.14048 (0.10616,0.17469) 0.47707 (0.41918,0.54504)
State 3 0.13737 (0.08844,0.19452) 0.72519 (0.64078,0.80585)
State 4 0                         1.00000 (1.00000,1.00000)
# try at home, this takes a long time to run;
# pmatrix.msm(heart.msm1, t=5, 
#             covariates = list(sex = 1, age = 65),
#              ci="bootstrap", cores = 4, B=1000);
  • For a 65 year old female patient in state 1, at 5 years there is about 51% chances of staying in state 1 and about 28% chances of death.

1.2.4 Comparing expected times from state i to j in by age and sex

efpt.msm(heart.msm3, tostate = 4, 
         covariates = list(sex = 1, age = 65), ci = "norm")
           [,1]     [,2]      [,3] [,4]
est    7.147621 6.193303 1.9663185    0
2.5%   3.917184 2.278243 0.8421432    0
97.5% 11.079076 9.589995 5.0063461    0
efpt.msm(heart.msm3, tostate = 4, 
         covariates = list(sex = 0, age = 65), ci = "norm")
           [,1]     [,2]     [,3] [,4]
est    8.637568 7.791753 4.804326    0
2.5%   5.361631 2.249057 3.038707    0
97.5% 10.480166 9.676978 6.785519    0

1.2.5 Survival plots

plot(heart.msm3)
abline(h = 0.5,lwd = 1, lty = 4)

  • For patients in state 3, the median survival time is approximately 3 years.

1.5 Model diagnostics

  • (Goodness-of-fit) Checking if the model predicted transition events matches with the expected values using Pearson chi-square test.
# plot.prevalence.msm(heart.msm3, maxtime=5)
pearson.msm(heart.msm3)
Imputing sampling times after deaths...
Calculating replicates of test statistics for imputations...
$Observed
           Time Interval           Cov 1-1 1-2 1-3   1-4 2-1 2-2 2-3   2-4 3-1
1      [0,3.04)        1 [-2.67,-1.23)  71   3   5  5.79   1   9   2  0.65   0
2   [3.04,5.97)        1 [-2.67,-1.23)  53   5   1  0.00   1   4   5  0.26   0
3  [5.97,19.46]        1 [-2.67,-1.23)  48   4   2  0.81   4   7   3  0.50   0
4      [0,3.04)        2 [-2.67,-1.23)  66  12   0  6.95   0   0   0  0.62   0
5   [3.04,5.97)        2 [-2.67,-1.23)  67   3   0  0.69   3   2   3  0.38   1
6  [5.97,19.46]        2 [-2.67,-1.23)  53   8   1  0.72   3   5   1  0.64   0
7      [0,3.04)        3 [-2.67,-1.23)  72   6   1  7.26   0   0   0  0.73   0
8   [3.04,5.97)        3 [-2.67,-1.23)  45   7   0  6.31   1   3   0  0.36   0
9  [5.97,19.46]        3 [-2.67,-1.23)  40   9   2 13.47   0   5   1  7.86   1
10     [0,3.04)        1 [-1.23,-1.18)  55   8   3  6.90   3   2   1  0.33   0
11  [3.04,5.97)        1 [-1.23,-1.18)  46   5   0  0.91   4   8   3  0.00   0
12 [5.97,19.46]        1 [-1.23,-1.18)  22   9   0  0.77   3   8   3  0.96   0
13     [0,3.04)        2 [-1.23,-1.18)  64   4   5  9.47   0   0   0  0.33   0
14  [3.04,5.97)        2 [-1.23,-1.18)  45   5   2  2.93   2   5   1  0.40   1
15 [5.97,19.46]        2 [-1.23,-1.18)  27   9   2  1.29   5   7   5  1.17   0
16     [0,3.04)        3 [-1.23,-1.18)  57   7   2  9.63   0   0   0  0.34   0
17  [3.04,5.97)        3 [-1.23,-1.18)  71  13   3  8.16   1   1   0  0.60   0
18 [5.97,19.46]        3 [-1.23,-1.18)  44  12   2 17.94   4   9   4  8.87   0
19     [0,3.04)        1 [-1.18,-1.18]  34   4   4  6.15   2   5   3  0.34   0
20  [3.04,5.97)        1 [-1.18,-1.18]  41   5   0  2.39   2  11   6  0.87   0
21 [5.97,19.46]        1 [-1.18,-1.18]  36  10   0  0.21   0  13   3  0.59   0
22     [0,3.04)        2 [-1.18,-1.18]  58   9   0  6.79   0   0   0  0.31   0
23  [3.04,5.97)        2 [-1.18,-1.18]  50   8   0  3.31   2   6   2  2.34   0
24 [5.97,19.46]        2 [-1.18,-1.18]  38   5   0  0.27   4   9   4  1.39   0
25     [0,3.04)        3 [-1.18,-1.18]  68  11   4  7.06   0   0   0  0.35   0
26  [3.04,5.97)        3 [-1.18,-1.18]  63  10   4  7.30   0   4   0  4.79   0
27 [5.97,19.46]        3 [-1.18,-1.18]  33  13   1 14.52   1  11   4 12.02   1
   3-2 3-3   3-4
1    1   0  0.26
2    0   7  1.41
3    0   1  1.19
4    0   0  0.45
5    1   3  1.51
6    0   4  2.12
7    0   0  0.29
8    0   0  2.08
9    1   2  8.69
10   1   6  0.00
11   0  11  0.91
12   2  16  0.86
13   0   0  0.00
14   0   6  1.05
15   0   8  1.10
16   0   0  0.00
17   0   0  1.04
18   0   3  9.04
19   0   0  0.39
20   1   7  1.51
21   3  11  0.90
22   0   0  0.38
23   0   9  1.90
24   2   8  0.78
25   0   0  0.23
26   0   1  3.59
27   1   4 13.32

$Expected
           Time Interval           Cov      1-1       1-2       1-3       1-4
1      [0,3.04)        1 [-2.67,-1.23) 73.99266  6.629811 1.4705270  2.697001
2   [3.04,5.97)        1 [-2.67,-1.23) 52.24885  4.307707 0.8026938  1.640747
3  [5.97,19.46]        1 [-2.67,-1.23) 48.75412  3.689966 0.6786088  1.687310
4      [0,3.04)        2 [-2.67,-1.23) 67.48990  8.794203 3.0472878  5.618609
5   [3.04,5.97)        2 [-2.67,-1.23) 58.80462  6.596561 2.0012754  3.287539
6  [5.97,19.46]        2 [-2.67,-1.23) 55.44513  4.728922 1.0134760  1.532473
7      [0,3.04)        3 [-2.67,-1.23) 66.61446  9.158596 3.2953783  7.191562
8   [3.04,5.97)        3 [-2.67,-1.23) 43.13881  6.045727 2.3456883  6.779772
9  [5.97,19.46]        3 [-2.67,-1.23) 43.18464  6.396571 2.7075471 12.181244
10     [0,3.04)        1 [-1.23,-1.18) 60.13583  7.456555 1.4201011  3.887511
11  [3.04,5.97)        1 [-1.23,-1.18) 44.07886  4.804001 0.7323918  2.294744
12 [5.97,19.46]        1 [-1.23,-1.18) 27.06567  2.923046 0.4409360  1.340348
13     [0,3.04)        2 [-1.23,-1.18) 59.93098 11.260026 3.4345268  7.844471
14  [3.04,5.97)        2 [-1.23,-1.18) 42.84279  6.446570 1.5207367  4.119908
15 [5.97,19.46]        2 [-1.23,-1.18) 33.00359  3.836631 0.6297117  1.820065
16     [0,3.04)        3 [-1.23,-1.18) 53.28776 10.608262 3.4031363  8.330840
17  [3.04,5.97)        3 [-1.23,-1.18) 64.37266 13.499322 4.5487024 12.739313
18 [5.97,19.46]        3 [-1.23,-1.18) 43.68023  9.705806 3.9651973 18.588765
19     [0,3.04)        1 [-1.18,-1.18] 38.89164  5.084137 0.9359668  3.238259
20  [3.04,5.97)        1 [-1.18,-1.18] 40.19017  4.791934 0.6808205  2.727071
21 [5.97,19.46]        1 [-1.18,-1.18] 38.41851  4.540437 0.6318390  2.619218
22     [0,3.04)        2 [-1.18,-1.18] 51.70027 10.419700 2.9548379  8.715189
23  [3.04,5.97)        2 [-1.18,-1.18] 45.48929  7.879662 1.9092168  6.031835
24 [5.97,19.46]        2 [-1.18,-1.18] 35.47335  4.461563 0.6647691  2.670319
25     [0,3.04)        3 [-1.18,-1.18] 61.72207 12.952605 3.9064138 11.478907
26  [3.04,5.97)        3 [-1.18,-1.18] 56.86071 12.156588 3.8809976 11.401702
27 [5.97,19.46]        3 [-1.18,-1.18] 36.06031  7.943090 3.0163293 14.500275
         2-1        2-2       2-3        2-4        3-1       3-2        3-3
1  1.6287342  7.3952260 2.6883512 0.93768853 0.02249462 0.2679959  0.6439402
2  1.4866841  5.5962807 2.2705216 0.90651356 0.10635048 0.9989452  4.3583074
3  1.7929956  8.3729530 3.3437095 0.99034198 0.02167880 0.1768307  1.1507088
4  0.0000000  0.0000000 0.0000000 0.10548035 0.00000000 0.0000000  0.0000000
5  1.6207859  3.7271490 1.7430827 1.28898245 0.10920162 0.8478023  3.3215248
6  1.4230529  5.2332014 2.1688967 0.81484905 0.07761483 0.6720132  3.5687008
7  0.0000000  0.0000000 0.0000000 0.13638969 0.00000000 0.0000000  0.0000000
8  0.9084358  1.4575218 1.0835716 0.91047077 0.00000000 0.0000000  0.0000000
9  2.3276416  3.1486164 2.5359817 5.84776035 0.41136293 1.1831721  3.6642351
10 1.0403809  3.7137826 1.1384750 0.43736157 0.06805836 0.5079541  4.7920071
11 2.3698262  8.7268480 2.8679476 1.03537822 0.09326334 0.7518337  8.5684504
12 2.3748623  8.6799241 2.8778539 1.02735963 0.15462919 1.2338246 13.4249532
13 0.0000000  0.0000000 0.0000000 0.05365178 0.00000000 0.0000000  0.0000000
14 1.6612771  4.1471028 1.7130029 0.87861714 0.10854505 0.6516639  5.0780400
15 3.1498444 10.1533107 3.4696938 1.39715108 0.09080660 0.6463431  6.1990962
16 0.0000000  0.0000000 0.0000000 0.05836913 0.00000000 0.0000000  0.0000000
17 0.6570770  0.7614939 0.6013876 0.58004158 0.00000000 0.0000000  0.0000000
18 5.8061097  7.3288356 5.3018514 7.43320325 0.28282874 0.5513458  3.2921453
19 1.8328797  5.9476446 1.7989891 0.76048660 0.00000000 0.0000000  0.0000000
20 3.3694482 11.7419976 3.3915851 1.36696920 0.07590723 0.5753670  6.8121893
21 2.7566470  9.8619140 2.8553085 1.11613059 0.11387043 0.8685705 10.8219722
22 0.0000000  0.0000000 0.0000000 0.05078972 0.00000000 0.0000000  0.0000000
23 2.5089360  6.0880091 2.4550694 1.28798556 0.11888921 0.6931327  7.3267621
24 3.3753176 10.2730632 3.2994040 1.44221531 0.09932641 0.6851073  7.4916662
25 0.0000000  0.0000000 0.0000000 0.06440950 0.00000000 0.0000000  0.0000000
26 2.2399636  2.6617404 2.0160903 1.87220564 0.14914700 0.3447976  1.8800131
27 5.8559917  8.1734296 5.1738825 8.81669623 0.55824506 1.0141059  5.4084323
           3-4
1   0.32556926
2   2.94639692
3   0.84078166
4   0.18198399
5   2.23147133
6   1.80167126
7   0.12260934
8   1.11451870
9   7.43122995
10  1.63198051
11  2.49645260
12  4.04659302
13  0.00000000
14  2.21175111
15  2.16375416
16  0.00000000
17  0.45628608
18  7.91368012
19  0.09428993
20  2.04653639
21  3.09558689
22  0.14969858
23  2.76121607
24  2.50390007
25  0.09370359
26  2.21604238
27 12.33921678

$`Deviance*sign(O-E)`
           Time Interval           Cov         1-1          1-2           1-3
1      [0,3.04)        1 [-2.67,-1.23) -0.15475809 -1.988093890  8.4840786793
2   [3.04,5.97)        1 [-2.67,-1.23)  0.01079874  0.111258782  0.0484988822
3  [5.97,19.46]        1 [-2.67,-1.23) -0.02231177  0.027179734  2.5747083030
4      [0,3.04)        2 [-2.67,-1.23) -0.06212636  1.176884279 -3.0472877751
5   [3.04,5.97)        2 [-2.67,-1.23)  1.15015608 -1.960996546 -2.0012753935
6  [5.97,19.46]        2 [-2.67,-1.23) -0.11533405  2.264668081 -0.0003081919
7      [0,3.04)        3 [-2.67,-1.23)  0.47576329 -1.091487716 -1.5990168392
8   [3.04,5.97)        3 [-2.67,-1.23)  0.08860401  0.151459878 -2.3456883398
9  [5.97,19.46]        3 [-2.67,-1.23) -0.24626540  1.067343081 -0.1873899257
10     [0,3.04)        1 [-1.23,-1.18) -0.47766859  0.047678825  1.7647646108
11  [3.04,5.97)        1 [-1.23,-1.18)  0.09569978  0.009234123 -0.7323917772
12 [5.97,19.46]        1 [-1.23,-1.18) -0.95761769 12.651559191 -0.4409359771
13     [0,3.04)        2 [-1.23,-1.18)  0.31903271 -4.681964980  0.7180900040
14  [3.04,5.97)        2 [-1.23,-1.18)  0.13108951 -0.327026846  0.1532394619
15 [5.97,19.46]        2 [-1.23,-1.18) -1.10246598  6.960137977  2.9847988039
16     [0,3.04)        3 [-1.23,-1.18)  0.32139564 -1.232179229 -0.5797279825
17  [3.04,5.97)        3 [-1.23,-1.18)  0.69553875 -0.020399980 -0.5276429466
18 [5.97,19.46]        3 [-1.23,-1.18)  0.02518538  0.551553144 -0.9752381633
19     [0,3.04)        1 [-1.18,-1.18] -0.67762243 -0.237874538 10.0675732171
20  [3.04,5.97)        1 [-1.18,-1.18]  0.05001154  0.013149160 -0.6808205409
21 [5.97,19.46]        1 [-1.18,-1.18] -0.15474240  6.566523066 -0.6318390103
22     [0,3.04)        2 [-1.18,-1.18]  0.81547275 -0.199205820 -2.9548378636
23  [3.04,5.97)        2 [-1.18,-1.18]  0.48020689  0.006806795 -1.9092168195
24 [5.97,19.46]        2 [-1.18,-1.18]  0.18401984  0.065602168 -0.6647691056
25     [0,3.04)        3 [-1.18,-1.18]  0.67755179 -0.299305462  0.0044487275
26  [3.04,5.97)        3 [-1.18,-1.18]  0.68722093 -0.385557494  0.0054031066
27 [5.97,19.46]        3 [-1.18,-1.18] -0.27004937  3.229934181 -1.3485232261
           1-4         2-1         2-2          2-3        2-4         3-1
1   4.64120037 -0.24442680  0.37050240 -0.181059052 -0.4479827 -0.02249462
2  -1.64074673 -0.16078117 -0.45993192  3.296491137 -0.6480399 -0.10635048
3  -0.85060468  2.72654007 -0.23535179 -0.041688860 -0.5977684 -0.02167880
4   0.80548072  0.00000000  0.00000000  0.000000000  2.5137652  0.00000000
5  -2.20357125  1.19456783 -0.80327180  0.920870234 -0.7971352  8.03451343
6  -0.80450020  1.76416697 -0.03686270 -0.633708635 -0.5900204 -0.07761483
7   0.45372315  0.00000000  0.00000000  0.000000000  2.6228571  0.00000000
8  -0.09527465  0.02652953  1.67808624 -1.083571555 -0.5131690  0.00000000
9   0.22138478 -2.32764164  1.18657997 -0.935541982  0.8858265  0.94991034
10  3.19469786  3.74035090 -0.79564552 -0.023318727 -0.4465402 -0.06805836
11 -1.11446627  1.12137613 -0.06053824  0.006080251 -1.0353782 -0.09326334
12 -0.63459791  0.17609160 -0.07481357  0.013486608 -0.5515796 -0.15462919
13  0.76447482  0.00000000  0.00000000  0.000000000  1.4237426  0.00000000
14 -0.67238343  0.07784689  0.19158720 -0.299371918 -0.4979606  7.50915680
15 -0.53315030  1.10495276 -0.99033429  0.691286322 -0.5064669 -0.09080660
16  0.73871897  0.00000000  0.00000000  0.000000000  1.3624554  0.00000000
17 -1.74087686  0.23225615  0.13864088 -0.601387568  0.2936101  0.00000000
18 -0.11294832 -0.57181989  0.41415233 -0.328433722  0.4363717 -0.28282874
19  3.75051726  0.01998584 -0.15831487  0.816117991 -0.4956147  0.00000000
20 -0.66414000 -0.55864159 -0.06337903  2.025058346 -0.6184944 -0.07590723
21 -2.27771444 -2.75664697  1.02339787  0.012397582 -0.5887606 -0.11387043
22 -0.80506556  0.00000000  0.00000000  0.000000000  1.3231172  0.00000000
23 -1.56232078 -0.12530767 -0.04467015 -0.105326428  1.5470220 -0.11888921
24 -2.22933407  0.12761208 -0.18035962  0.163891887 -0.5717103 -0.09932641
25 -2.01521827  0.00000000  0.00000000  0.000000000  1.2961526  0.00000000
26 -1.66133824 -2.23996357  0.84380950 -2.016090348  4.8285559 -0.14914700
27  0.05578296 -4.02757244  1.02041106 -0.278926409  1.3474795  0.39944055
           3-2         3-3        3-4
1   2.22487061 -0.64394020 -0.3023043
2  -0.99894516  1.73940978 -1.0165347
3  -0.17683072 -0.29541459  0.5674366
4   0.00000000  0.00000000  0.3947771
5   0.06464109 -0.08581761 -0.5446536
6  -0.67201315  0.18033222  0.3093050
7   0.00000000  0.00000000  0.2308439
8   0.00000000  0.00000000  0.9352317
9  -0.04502819 -0.77419953  0.2695276
10  0.47663601  0.30451686 -1.6319805
11 -0.75183370  0.74774438 -1.2223799
12  0.48000801  0.53138177 -2.6596617
13  0.00000000  0.00000000  0.0000000
14 -0.65166388  0.22527830 -0.8076982
15 -0.64634306  0.62675041 -0.8041588
16  0.00000000  0.00000000  0.0000000
17  0.00000000  0.00000000  0.7748781
18 -0.55134579 -0.08332330  0.2062846
19  0.00000000  0.00000000  1.0088393
20  0.33082449  0.07540236 -0.4477528
21  5.26649396  0.03771463 -1.7590065
22  0.00000000  0.00000000  0.3543592
23 -0.69313266  0.46981247 -0.5299021
24  2.54920586  0.07233259 -1.3636847
25  0.00000000  0.00000000  0.1989976
26 -0.34479756 -0.46224726  0.9346160
27 -0.01286860 -0.40302849  0.1234814

$test
     stat df.lower p.lower df.upper   p.upper
 283.2702       NA      NA      237 0.0211275
  • The p-value suggest the model doesn’t fit well.

2. Hands-on excerise fitting multistate COX model for time-to-event data with complete observed transition times

  • using R packages survival and mstate and mvna, package references from the developer:
  • data contains complete information on state changes
    • know the complete process history and event times, i.e., changes of state represent events
  • can estimate transition rates under more flexible models, such as COX semi-Markov model.

2.1 Hospital acquired penumonia in ICU (illness death model)

  • This data set is a random sample drawn from the SIR-3 study that aimed at analyzing the effect of nosocomial infections on the length of ICU stay.

  • Patients were included in the study if they had stayed at least 1 day in the unit. The sample includes information to assess the effect of nosocomial pneumonia on the length of stay.

  • The endpoint is either discharge alive from the ICU or dead in the unit. These data are censoring complete as the censoring time is known for all patients.

  • Beyersmann, J., Gastmeier, P., Grundmann, H., Baerwolff, S., Geffers, C., Behnke, M., Rueden, H., and Schumacher, M. Use of multistate models to assess prolongation of intensive care unit stay due to nosocomial infection. Infection Control and Hospital Epidemiology, 27:493-499, 2006.

  • A data frame with 1421 observations on the following 8 variables.

    • id
    • start, Start of the observation time.
    • stop, Failure time.
    • status, Censoring status.
    • event, 2 is death in ICU, 3 is discharge alive
    • pneu, Nosocomial pneumonia indicator.
    • adm.cens.exit, Exit times for patients discharged alive are replaced by their administrative censoring times.
    • age, Age at inclusion
    • sex, F for female and M for male

library(survival)
library(mstate)
library(mvna) #non-parametric analysis;
library(kmi) #getting data;

data("icu.pneu")

# look at the data;
datatable(icu.pneu,
          rownames = FALSE) %>%
  formatRound(columns=c('age'), digits=1)

2.2 Non-parametric analysis

library(lattice)
# manipulate orignal data to fit the illness-death model;
my.icu.pneu <- icu.pneu %>%
           group_by(id) %>%
            mutate(leadpneu = lead(pneu)) %>%
            mutate(from = case_when(
              start ==0 ~ 0,
              start >0 & pneu ==1 ~ 1),
              to = case_when(
                from ==0 & leadpneu==1 ~ 1,
                from ==0 & is.na(leadpneu) & event ==2 ~ 2,
                from ==0 & is.na(leadpneu) & event ==3 ~ 3,
                from ==0 & leadpneu==0 & event ==2 ~ 2,
                from ==0 & leadpneu==0 & event ==3 ~ 3,
                from ==1 & event == 2 ~ 2,
                from ==1 & event == 3 ~ 3)) %>%
  select(-leadpneu, -adm.cens.exit, - pneu) %>%
  rename(entry = start, exit = stop)

# Matrix of logical indicating the possible transitions;
# state 0 admission, 1 hospital acquired penumonia,
# 2, death, 3 discharge alive;
tra <- matrix(ncol=4,nrow=4,FALSE)
tra[1,2:4] <- TRUE
tra[2,3:4] <- TRUE
tra
      [,1]  [,2]  [,3]  [,4]
[1,] FALSE  TRUE  TRUE  TRUE
[2,] FALSE FALSE  TRUE  TRUE
[3,] FALSE FALSE FALSE FALSE
[4,] FALSE FALSE FALSE FALSE
# Nelson-Aalen estimator of cumulative transition hazards
pneu <- mvna(my.icu.pneu, c("0", "1", "2","3"), tra, cens.name = NULL)
pneu
Transition 0 -> 1 
   na var.aalen var.greenwood time
 0.00      0.00          0.00    0
 0.13      0.00          0.00   17
 0.20      0.00          0.00   35
 0.22      0.00          0.00   54
 0.34      0.02          0.01   71
 0.34      0.02          0.01   96

Transition 0 -> 2 
   na var.aalen var.greenwood time
 0.00      0.00          0.00    0
 0.13      0.00          0.00   17
 0.32      0.00          0.00   35
 0.44      0.00          0.00   54
 0.79      0.03          0.03   71
 1.29      0.28          0.15   96

Transition 0 -> 3 
   na var.aalen var.greenwood time
 0.00      0.00          0.00    0
 1.25      0.00          0.00   17
 2.18      0.01          0.01   35
 3.29      0.04          0.03   54
 3.98      0.09          0.09   71
 4.93      0.34          0.26   96

Transition 1 -> 2 
   na var.aalen var.greenwood time
 0.00      0.00          0.00    4
 0.13      0.00          0.00   19
 0.25      0.01          0.01   35
 0.39      0.01          0.01   52
 0.66      0.05          0.04   63
 1.24      0.22          0.16   81

Transition 1 -> 3 
   na var.aalen var.greenwood time
 0.00      0.00          0.00    4
 0.39      0.01          0.01   19
 1.32      0.03          0.03   35
 1.99      0.07          0.06   52
 2.53      0.13          0.11   63
 2.98      0.23          0.19   81
xyplot(pneu,
        tr.choice=c("0 2","1 2","0 3","1 3"),
        aspect=1,
        strip=strip.custom(bg="white",
                           factor.levels= c("No pneumonia - Death","Pneumonia - Death","No pneumonia - Discharge","Pneumonia - Discharge"),par.strip.text=list(cex=0.9)), scales=list(alternating=1),xlab="Days", ylab="Nelson-Aalen estimates")

# estimating transition probabilities using etm package;
# Aalen-Johansen estimator;
library(etm)
pneu.etm <- etm(my.icu.pneu, c("0", "1", "2","3"), tra, cens.name = NULL, s = 0) #s is the starting state for probability calculation;

summary(pneu.etm)
Transition 0 0 
           P time          var        lower       upper n.risk n.event
 0.999238385    1 5.796151e-07 0.9977462167 1.000000000   1313       0
 0.187357197   18 1.159592e-04 0.1662514484 0.208462946    265       0
 0.055597867   36 3.998990e-05 0.0432035312 0.067992204     78       0
 0.015993907   54 1.198637e-05 0.0092082500 0.022779564     23       0
 0.004569688   71 3.464437e-06 0.0009216072 0.008217768      8       0
 0.000000000  460 0.000000e+00 0.0000000000 0.000000000      1       0

Transition 0 1 
           P time          var        lower       upper n.risk n.event
 0.000000000    1 0.000000e+00 0.000000e+00 0.000000000   1313       0
 0.051028180   18 3.688066e-05 3.912543e-02 0.062930933    265       0
 0.020563595   36 1.533948e-05 1.288727e-02 0.028239920     78       0
 0.009139375   54 6.897066e-06 3.992066e-03 0.014286685     23       0
 0.003046458   71 2.313159e-06 6.553506e-05 0.006027382      8       1
 0.000000000  460 0.000000e+00 0.000000e+00 0.000000000      1       0

Transition 0 2 
          P time          var      lower      upper n.risk n.event
 0.00000000    1 0.000000e+00 0.00000000 0.00000000   1313       0
 0.07006855   18 4.962600e-05 0.05626144 0.08387565    265       1
 0.09672506   36 6.654175e-05 0.08073702 0.11271310     78       0
 0.10357959   54 7.071657e-05 0.08709764 0.12006154     23       0
 0.10967251   71 7.436744e-05 0.09277045 0.12657456      8       0
 0.11195735  460 7.572194e-05 0.09490207 0.12901263      1       0

Transition 0 3 
            P time          var     lower       upper n.risk n.event
 0.0007616146    1 5.796151e-07 0.0000000 0.002253783   1313       1
 0.6915460777   18 1.624601e-04 0.6665644 0.716527746    265      18
 0.8271134806   36 1.089084e-04 0.8066595 0.847567508     78       5
 0.8712871287   54 8.541193e-05 0.8531734 0.889400837     23       2
 0.8827113481   71 7.885150e-05 0.8653072 0.900115509      8       1
 0.8880426504  460 7.572194e-05 0.8709874 0.905097934      1       1

Transition 1 1 
          P time          var      lower      upper n.risk n.event
 1.00000000    1 0.0000000000 1.00000000 1.00000000      0       0
 0.62231062   18 0.0034982409 0.50638673 0.73823451     72       0
 0.19637809   36 0.0014274534 0.12232741 0.27042877     27       0
 0.08416204   54 0.0005994829 0.03617362 0.13215046     12       0
 0.02104051   71 0.0001481434 0.00000000 0.04489605      4       0
 0.00000000  460 0.0000000000 0.00000000 0.00000000      0       0

Transition 1 2 
          P time         var      lower     upper n.risk n.event
 0.00000000    1 0.000000000 0.00000000 0.0000000      0       0
 0.08523632   18 0.001749435 0.00325837 0.1672143     72       1
 0.15369166   36 0.002061832 0.06469478 0.2426885     27       0
 0.17473217   54 0.002119916 0.08449044 0.2649739     12       0
 0.19577268   71 0.002179151 0.10427886 0.2872665      4       1
 0.20629293  460 0.002199978 0.11436294 0.2982229      0       0

Transition 1 3 
         P time         var     lower     upper n.risk n.event
 0.0000000    1 0.000000000 0.0000000 0.0000000      0       0
 0.2924531   18 0.002958697 0.1858430 0.3990632     72       4
 0.6499303   36 0.002648728 0.5490592 0.7508014     27       0
 0.7411058   54 0.002363766 0.6458151 0.8363964     12       0
 0.7831868   71 0.002239538 0.6904340 0.8759397      4       0
 0.7937071  460 0.002199978 0.7017771 0.8856371      0       0
par(mfrow=c(2,2))
plot(pneu.etm, tr.choice = "0 1", conf.int = TRUE,lwd = 2, legend = FALSE, ylim = c(0, 0.1),xlim = c(0, 100), xlab = "Days", main= "No pneumonia to hospital acquired penumonia", ci.fun = "cloglog")
plot(pneu.etm, tr.choice = "0 2", conf.int = TRUE,lwd = 2, legend = FALSE, ylim = c(0, 0.1),xlim = c(0, 100), xlab = "Days", main= "No pneumonia to death", ci.fun = "cloglog")
plot(pneu.etm, tr.choice = "1 2", conf.int = TRUE,lwd = 2, legend = FALSE, ylim = c(0, 0.1),xlim = c(0, 100), xlab = "Days", main= "hospital acquired penumonia to death", ci.fun = "cloglog")
plot(pneu.etm, tr.choice = "1 3", conf.int = TRUE,lwd = 2, legend = FALSE, ylim = c(0, 0.1),xlim = c(0, 100), xlab = "Days", main= "hospital acquired penumonia to discarge alive", ci.fun = "cloglog")

  • shows that the proportion of hospital infected patients increases for about the first 20 days after admission and then decreases again.

2.3 Proportional transition hazards models

my.icu.pneu <- my.icu.pneu %>%
  mutate(trans = case_when(
    from == 0 & to ==1 ~ 1,
    from == 0 & to ==2 ~ 2,
    from == 0 & to ==3 ~ 3,
    from == 1 & to ==2 ~ 4,
    from == 1 & to ==3 ~ 5,
  ),
  event2 = case_when(
    to == 3 ~ 0,
    to == 1 ~ 1,
    to == 2 ~ 2
  ))


# look at the data prior to model fit;
#  this would be the data format we want for fitting the coxph model;
datatable(my.icu.pneu[, c(-3,-4)],
          rownames = FALSE) %>%
  formatRound(columns=c('age'), digits=1)
trans1 <- coxph(Surv(entry,exit,as.factor(event2)) ~ age+sex,id=id, method="breslow", subset=(from==0),data=my.icu.pneu)

summary(trans1) # 1:2 is from state 0 to state 1 (penumonia); 1:3 is from stat 0 to state 3 (death)
Call:
coxph(formula = Surv(entry, exit, as.factor(event2)) ~ age + 
    sex, data = my.icu.pneu, subset = (from == 0), method = "breslow", 
    id = id)

  n= 1313, number of events= 234 

             coef exp(coef) se(coef) robust se     z Pr(>|z|)    
age_1:2  0.007921  1.007952 0.005632  0.005506 1.439    0.150    
sexM_1:2 0.209017  1.232466 0.201121  0.202138 1.034    0.301    
age_1:3  0.022675  1.022934 0.005732  0.005634 4.025 5.71e-05 ***
sexM_1:3 0.109237  1.115427 0.184914  0.182553 0.598    0.550    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

         exp(coef) exp(-coef) lower .95 upper .95
age_1:2      1.008     0.9921    0.9971     1.019
sexM_1:2     1.232     0.8114    0.8293     1.832
age_1:3      1.023     0.9776    1.0117     1.034
sexM_1:3     1.115     0.8965    0.7799     1.595

Concordance= 0.571  (se = 0.022 )
Likelihood ratio test= 20.32  on 4 df,   p=4e-04
Wald test            = 18.89  on 4 df,   p=8e-04
Score (logrank) test = 18.82  on 4 df,   p=9e-04,   Robust = 20.63  p=4e-04

  (Note: the likelihood ratio and score tests assume independence of
     observations within a cluster, the Wald and robust score tests do not).
trans2 <- coxph(Surv(entry,exit,as.factor(event2)) ~ age+sex,id=id, method="breslow", subset=(from==1),data=my.icu.pneu)

summary(trans2) # from stat 1 (penumonia) to state 3 (death);
Call:
coxph(formula = Surv(entry, exit, as.factor(event2)) ~ age + 
    sex, data = my.icu.pneu, subset = (from == 1), method = "breslow", 
    id = id)

  n= 108, number of events= 21 

         coef exp(coef) se(coef) robust se      z Pr(>|z|)
age  -0.01251   0.98757  0.01589   0.01408 -0.888    0.374
sexM -0.62746   0.53395  0.44016   0.41133 -1.525    0.127

     exp(coef) exp(-coef) lower .95 upper .95
age     0.9876      1.013    0.9607     1.015
sexM    0.5339      1.873    0.2384     1.196

Concordance= 0.614  (se = 0.073 )
Likelihood ratio test= 2.71  on 2 df,   p=0.3
Wald test            = 4.03  on 2 df,   p=0.1
Score (logrank) test = 2.86  on 2 df,   p=0.2,   Robust = 3.1  p=0.2

  (Note: the likelihood ratio and score tests assume independence of
     observations within a cluster, the Wald and robust score tests do not).
trans3 <- coxph(Surv(entry,exit,as.factor(event2)) ~ age+sex+strata(trans),id=id, method="breslow", data=my.icu.pneu)

summary(trans3)
Call:
coxph(formula = Surv(entry, exit, as.factor(event2)) ~ age + 
    sex + strata(trans), data = my.icu.pneu, method = "breslow", 
    id = id)

  n= 1421, number of events= 255 

              coef exp(coef)  se(coef) robust se      z Pr(>|z|)
age_1:2  -0.006546  0.993475  0.006306  0.005801 -1.128    0.259
sexM_1:2  0.107620  1.113625  0.208430  0.193554  0.556    0.578
age_1:3  -0.004240  0.995769  0.006371  0.005521 -0.768    0.443
sexM_1:3 -0.054125  0.947313  0.186267  0.180335 -0.300    0.764
age_2:3  -0.017089  0.983056  0.017713  0.013692 -1.248    0.212
sexM_2:3  0.217264  1.242672  0.474127  0.333002  0.652    0.514

         exp(coef) exp(-coef) lower .95 upper .95
age_1:2     0.9935     1.0066    0.9822     1.005
sexM_1:2    1.1136     0.8980    0.7621     1.627
age_1:3     0.9958     1.0042    0.9851     1.007
sexM_1:3    0.9473     1.0556    0.6653     1.349
age_2:3     0.9831     1.0172    0.9570     1.010
sexM_2:3    1.2427     0.8047    0.6470     2.387

Concordance= 0.507  (se = 0.023 )
Likelihood ratio test= 2.72  on 6 df,   p=0.8
Wald test            = 3.22  on 6 df,   p=0.8
Score (logrank) test = 2.72  on 6 df,   p=0.8,   Robust = 3.12  p=0.8

  (Note: the likelihood ratio and score tests assume independence of
     observations within a cluster, the Wald and robust score tests do not).

3. Quick discussion

From (Meira-Machado et al. 2009)

  • A multi-state model (MSM) is a model for a continuous time stochastic process allowing individuals to move among a finite number of states.
  • In an MSM, the transition intensities provide the hazards for movement from one state to another.

Different model assumptions can be made about the dependence of the transition rates on time.

  1. Time homogeneous models: the intensities are constant over time
  2. Markov models: the transition intensities only depend on the history of the process through the current state.
    • Because of their simplicity, Markov models have been most frequently used.
    • modelling movements of patients between units of a hospital or between stages of a disease, or in the prediction of the incidence rate of diseases
  3. Semi-Markov models (an extension of markov models): future evolution depends on the current state and the entry time into state that state. Future probability transitions depend on the sojourn time (the time spent in the current state), and the clock of each sojourn time is reset to zero after each transition into a new state.
    • some aspects of systems’ behaviour cannot be captured by Markov processes. For instance, the risk of chronic diseases such as AIDS essentially depends on the time since infection.
    • SMP has been used for predicting a disease progression and recovery progress of patients with a specific type of disease.
    • illness-death model is a three-state SMP.
  4. Standard Multistate COX Model
    • Markov Assumption: In its standard form, the multistate Cox model typically adheres to the Markov property.
    • That is, the hazard (transition intensity) for moving to another state depends only on the current state and possibly time, not on the duration already spent in the current state.
  5. Multistate COX Model extending to Semi-Markov Framework
    • Incorporating Duration Dependence: By allowing the transition hazards to depend on the time since entering the current state, the multistate Cox model can be extended to a semi-Markov framework.
    • Time-Dependent Covariates: Incorporating time-dependent covariates or specifying hazard functions that vary with the time spent in the current state can facilitate a semi-Markov structure within the Cox model.

References

Meira-Machado, Luı́s, Jacobo de Uña-Álvarez, Carmen Cadarso-Suárez, and Per K Andersen. 2009. “Multi-State Models for the Analysis of Time-to-Event Data.” Statistical Methods in Medical Research 18 (2): 195–222.