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 (IQR); 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.47680,0.54624) 0.13235 (0.11302,0.15265)
State 2 0.24973 (0.20002,0.30445) 0.13272 (0.10401,0.16425)
State 3 0.06806 (0.04021,0.10692) 0.06938 (0.04226,0.10219)
State 4 0                         0                        
        State 3                   State 4                  
State 1 0.07304 (0.05820,0.08619) 0.28293 (0.25324,0.31586)
State 2 0.14048 (0.10520,0.17748) 0.47707 (0.41579,0.54218)
State 3 0.13737 (0.08933,0.19817) 0.72519 (0.63540,0.80756)
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.47783,0.54627) 0.13235 (0.11449,0.15153)
State 2 0.24973 (0.19856,0.30540) 0.13272 (0.10487,0.16413)
State 3 0.06806 (0.04181,0.10623) 0.06938 (0.04462,0.10077)
State 4 0                         0                        
        State 3                   State 4                  
State 1 0.07304 (0.05845,0.08688) 0.28293 (0.25206,0.31657)
State 2 0.14048 (0.10418,0.17425) 0.47707 (0.41670,0.54623)
State 3 0.13737 (0.08822,0.19777) 0.72519 (0.64285,0.80373)
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.858750 1.549545 0.8251825    0
97.5% 11.441801 9.634627 5.1014553    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.449615 2.535757 2.968977    0
97.5% 10.302984 9.532278 6.592884    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.86   1   9   2  0.75   0
2   [3.04,5.97)        1 [-2.67,-1.23)  53   5   1  0.00   1   4   5  0.29   0
3  [5.97,19.46]        1 [-2.67,-1.23)  48   4   2  0.77   4   7   3  0.50   0
4      [0,3.04)        2 [-2.67,-1.23)  66  12   0  6.34   0   0   0  0.70   0
5   [3.04,5.97)        2 [-2.67,-1.23)  67   3   0  0.71   3   2   3  0.27   1
6  [5.97,19.46]        2 [-2.67,-1.23)  53   8   1  1.02   3   5   1  0.67   0
7      [0,3.04)        3 [-2.67,-1.23)  72   6   1  7.80   0   0   0  0.55   0
8   [3.04,5.97)        3 [-2.67,-1.23)  45   7   0  6.29   1   3   0  0.44   0
9  [5.97,19.46]        3 [-2.67,-1.23)  40   9   2 13.21   0   5   1  7.83   1
10     [0,3.04)        1 [-1.23,-1.18)  55   8   3  7.29   3   2   1  0.43   0
11  [3.04,5.97)        1 [-1.23,-1.18)  46   5   0  0.98   4   8   3  0.00   0
12 [5.97,19.46]        1 [-1.23,-1.18)  22   9   0  0.81   3   8   3  1.04   0
13     [0,3.04)        2 [-1.23,-1.18)  64   4   5  9.46   0   0   0  0.33   0
14  [3.04,5.97)        2 [-1.23,-1.18)  45   5   2  3.11   2   5   1  0.40   1
15 [5.97,19.46]        2 [-1.23,-1.18)  27   9   2  1.40   5   7   5  1.10   0
16     [0,3.04)        3 [-1.23,-1.18)  57   7   2  9.25   0   0   0  0.24   0
17  [3.04,5.97)        3 [-1.23,-1.18)  71  13   3  7.91   1   1   0  0.60   0
18 [5.97,19.46]        3 [-1.23,-1.18)  44  12   2 17.79   4   9   4  8.86   0
19     [0,3.04)        1 [-1.18,-1.18]  34   4   4  6.15   2   5   3  0.39   0
20  [3.04,5.97)        1 [-1.18,-1.18]  41   5   0  2.32   2  11   6  0.96   0
21 [5.97,19.46]        1 [-1.18,-1.18]  36  10   0  0.23   0  13   3  0.80   0
22     [0,3.04)        2 [-1.18,-1.18]  58   9   0  6.87   0   0   0  0.27   0
23  [3.04,5.97)        2 [-1.18,-1.18]  50   8   0  3.05   2   6   2  2.27   0
24 [5.97,19.46]        2 [-1.18,-1.18]  38   5   0  0.34   4   9   4  1.47   0
25     [0,3.04)        3 [-1.18,-1.18]  68  11   4  6.98   0   0   0  0.34   0
26  [3.04,5.97)        3 [-1.18,-1.18]  63  10   4  7.63   0   4   0  4.77   0
27 [5.97,19.46]        3 [-1.18,-1.18]  33  13   1 14.43   1  11   4 11.73   1
   3-2 3-3   3-4
1    1   0  0.22
2    0   7  1.27
3    0   1  1.41
4    0   0  0.39
5    1   3  1.57
6    0   4  1.77
7    0   0  0.39
8    0   0  2.16
9    1   2  8.82
10   1   6  0.00
11   0  11  0.96
12   2  16  0.64
13   0   0  0.00
14   0   6  1.11
15   0   8  1.06
16   0   0  0.00
17   0   0  0.93
18   0   3  9.30
19   0   0  0.39
20   1   7  1.59
21   3  11  0.77
22   0   0  0.35
23   0   9  1.95
24   2   8  0.87
25   0   0  0.26
26   0   1  3.46
27   1   4 13.36

$Expected
           Time Interval           Cov      1-1       1-2       1-3       1-4
1      [0,3.04)        1 [-2.67,-1.23) 74.09365  6.626344 1.4670353  2.672967
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.71966  3.687062 0.6781815  1.685100
4      [0,3.04)        2 [-2.67,-1.23) 67.00572  8.726522 3.0233281  5.584430
5   [3.04,5.97)        2 [-2.67,-1.23) 58.82993  6.595739 2.0018825  3.282449
6  [5.97,19.46]        2 [-2.67,-1.23) 55.70778  4.752710 1.0183301  1.541180
7      [0,3.04)        3 [-2.67,-1.23) 67.02334  9.217517 3.3171754  7.241968
8   [3.04,5.97)        3 [-2.67,-1.23) 43.09655  6.044553 2.3444202  6.804477
9  [5.97,19.46]        3 [-2.67,-1.23) 43.07550  6.392620 2.7158317 12.026053
10     [0,3.04)        1 [-1.23,-1.18) 60.44828  7.495771 1.4281502  3.917803
11  [3.04,5.97)        1 [-1.23,-1.18) 44.13860  4.810388 0.7333926  2.297617
12 [5.97,19.46]        1 [-1.23,-1.18) 27.09171  2.928338 0.4411905  1.348763
13     [0,3.04)        2 [-1.23,-1.18) 59.95063 11.256217 3.4356025  7.817555
14  [3.04,5.97)        2 [-1.23,-1.18) 42.95285  6.475975 1.5303754  4.150799
15 [5.97,19.46]        2 [-1.23,-1.18) 33.09658  3.847129 0.6313321  1.824959
16     [0,3.04)        3 [-1.23,-1.18) 52.98967 10.557638 3.3856750  8.317019
17  [3.04,5.97)        3 [-1.23,-1.18) 64.14694 13.453401 4.5373238 12.772333
18 [5.97,19.46]        3 [-1.23,-1.18) 43.51036  9.660219 3.9434897 18.675936
19     [0,3.04)        1 [-1.18,-1.18] 38.92566  5.071838 0.9298637  3.222639
20  [3.04,5.97)        1 [-1.18,-1.18] 40.12680  4.787370 0.6804085  2.725426
21 [5.97,19.46]        1 [-1.18,-1.18] 38.43489  4.542420 0.6320830  2.620608
22     [0,3.04)        2 [-1.18,-1.18] 51.75271 10.431236 2.9573631  8.728694
23  [3.04,5.97)        2 [-1.18,-1.18] 45.30714  7.843147 1.8989520  6.000757
24 [5.97,19.46]        2 [-1.18,-1.18] 35.52958  4.468962 0.6657721  2.675682
25     [0,3.04)        3 [-1.18,-1.18] 61.67547 12.940312 3.9015943 11.462620
26  [3.04,5.97)        3 [-1.18,-1.18] 57.03689 12.194893 3.8961619 11.502050
27 [5.97,19.46]        3 [-1.18,-1.18] 35.91090  7.911285 3.0057060 14.602113
         2-1        2-2       2-3        2-4        3-1       3-2        3-3
1  1.6387346  7.4566923 2.7120433 0.94252984 0.02185887 0.2615523  0.6212315
2  1.4915041  5.6120683 2.2764137 0.91001389 0.10446714 0.9782857  4.2961780
3  1.7916546  8.3762222 3.3426563 0.98946690 0.02451451 0.2031581  1.2665067
4  0.0000000  0.0000000 0.0000000 0.11895397 0.00000000 0.0000000  0.0000000
5  1.5976753  3.6859449 1.7219472 1.26443253 0.11337172 0.8660307  3.3229902
6  1.4260605  5.2496092 2.1771763 0.81715401 0.07432792 0.6562137  3.4069332
7  0.0000000  0.0000000 0.0000000 0.09963362 0.00000000 0.0000000  0.0000000
8  0.9261524  1.4740232 1.0944034 0.94542086 0.00000000 0.0000000  0.0000000
9  2.3279779  3.1317918 2.5265358 5.84369446 0.42230814 1.1734043  3.6246659
10 1.0556283  3.7721402 1.1583267 0.44390466 0.06805836 0.5079541  4.7920071
11 2.3698262  8.7268480 2.8679476 1.03537822 0.09406839 0.7573153  8.5964304
12 2.3872609  8.7267800 2.8933448 1.03261433 0.15290824 1.2221039 13.2668696
13 0.0000000  0.0000000 0.0000000 0.05345457 0.00000000 0.0000000  0.0000000
14 1.6609261  4.1481972 1.7127370 0.87813968 0.10903966 0.6559886  5.1199943
15 3.1376517 10.1141743 3.4564019 1.39177212 0.09043464 0.6436777  6.1715360
16 0.0000000  0.0000000 0.0000000 0.04258348 0.00000000 0.0000000  0.0000000
17 0.6565877  0.7638851 0.6016002 0.57792703 0.00000000 0.0000000  0.0000000
18 5.8147352  7.2955723 5.3012322 7.44846025 0.29213956 0.5648314  3.3444704
19 1.8390057  5.9826647 1.8069959 0.76133363 0.00000000 0.0000000  0.0000000
20 3.3843672 11.7936595 3.4086275 1.37334585 0.07662187 0.5811425  6.8680781
21 2.7929008  9.9820665 2.8936906 1.13134208 0.11264279 0.8602134 10.7326592
22 0.0000000  0.0000000 0.0000000 0.04400486 0.00000000 0.0000000  0.0000000
23 2.4965513  6.0551957 2.4370468 1.28120628 0.11910759 0.6946587  7.3631659
24 3.3882148 10.3176320 3.3160520 1.44810132 0.10017271 0.6907188  7.5539912
25 0.0000000  0.0000000 0.0000000 0.06275657 0.00000000 0.0000000  0.0000000
26 2.2284940  2.6238146 2.0038172 1.91387414 0.14197043 0.3349345  1.8395779
27 5.8063243  8.0677127 5.1192075 8.73675555 0.54323395 1.0089964  5.4204914
           3-4
1   0.31535736
2   2.89106919
3   0.91582066
4   0.15763468
5   2.26760746
6   1.63252519
7   0.16536285
8   1.15317599
9   7.59962166
10  1.63198051
11  2.51218594
12  3.99811820
13  0.00000000
14  2.22497750
15  2.15435164
16  0.00000000
17  0.41109626
18  8.09855868
19  0.09341149
20  2.06415753
21  3.06448455
22  0.13774001
23  2.77306780
24  2.52511730
25  0.10532310
26  2.14351716
27 12.38727823

$`Deviance*sign(O-E)`
           Time Interval           Cov         1-1          1-2           1-3
1      [0,3.04)        1 [-2.67,-1.23) -0.16019945 -1.985194945  8.5174738938
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.02104464  0.027676599  2.5779820383
4      [0,3.04)        2 [-2.67,-1.23) -0.04767243  1.236264261 -3.0233280978
5   [3.04,5.97)        2 [-2.67,-1.23)  1.14227904 -1.960336845 -2.0018825162
6  [5.97,19.46]        2 [-2.67,-1.23) -0.14054498  2.221142382 -0.0004894192
7      [0,3.04)        3 [-2.67,-1.23)  0.40919902 -1.125096012 -1.6187861717
8   [3.04,5.97)        3 [-2.67,-1.23)  0.09334872  0.152039071 -2.3444202151
9  [5.97,19.46]        3 [-2.67,-1.23) -0.23176471  1.070762237 -0.1907375693
10     [0,3.04)        1 [-1.23,-1.18) -0.54618684  0.044745112  1.7392475703
11  [3.04,5.97)        1 [-1.23,-1.18)  0.09185312  0.008831620 -0.7333926104
12 [5.97,19.46]        1 [-1.23,-1.18) -0.96438791 12.603316546 -0.4411905344
13     [0,3.04)        2 [-1.23,-1.18)  0.33321862 -4.679007694  0.7187054392
14  [3.04,5.97)        2 [-1.23,-1.18)  0.12425202 -0.339505966  0.1471574585
15 [5.97,19.46]        2 [-1.23,-1.18) -1.13663997  6.916247318  2.9709151248
16     [0,3.04)        3 [-1.23,-1.18)  0.37196649 -1.204203703 -0.5685095066
17  [3.04,5.97)        3 [-1.23,-1.18)  0.75012162 -0.018307611 -0.5215471058
18 [5.97,19.46]        3 [-1.23,-1.18)  0.02173105  0.573029039 -0.9587934289
19     [0,3.04)        1 [-1.18,-1.18] -0.67808312 -0.233397008 10.1776685908
20  [3.04,5.97)        1 [-1.18,-1.18]  0.04258332  0.012395615 -0.6804084699
21 [5.97,19.46]        1 [-1.18,-1.18] -0.15691544  6.559006337 -0.6320830085
22     [0,3.04)        2 [-1.18,-1.18]  0.80611167 -0.202641031 -2.9573630991
23  [3.04,5.97)        2 [-1.18,-1.18]  0.52286924  0.008698754 -1.8989520062
24 [5.97,19.46]        2 [-1.18,-1.18]  0.17638433  0.063806965 -0.6657720652
25     [0,3.04)        3 [-1.18,-1.18]  0.68829680 -0.295848569  0.0046041688
26  [3.04,5.97)        3 [-1.18,-1.18]  0.64121938 -0.397204454  0.0040904583
27 [5.97,19.46]        3 [-1.18,-1.18] -0.24581554  3.284449074 -1.3390680250
           1-4         2-1         2-2          2-3        2-4         3-1
1   4.78642906 -0.25065717  0.34464448 -0.192039162 -0.4349890 -0.02185887
2  -1.64074673 -0.16352265 -0.46799165  3.274909206 -0.6217696 -0.10446714
3  -0.88310111  2.73090404 -0.23617027 -0.041056001 -0.5821831 -0.02451451
4   0.63707678  0.00000000  0.00000000  0.000000000  2.8407409  0.00000000
5  -2.15928252  1.24822592 -0.77366902  0.960897649 -0.9155908  7.52563394
6  -0.61999884  1.75343343 -0.03789057 -0.640149152 -0.5687480 -0.07432792
7   0.47791915  0.00000000  0.00000000  0.000000000  2.0521212  0.00000000
8  -0.10846072  0.02401168  1.62691494 -1.094403448 -0.4545838  0.00000000
9   0.20590460 -2.32797790  1.20248628 -0.927618197  0.8584319  0.88774875
10  4.11221892  3.63293905 -0.83798910 -0.028498101 -0.4636808 -0.06805836
11 -1.06621165  1.12137613 -0.06053824  0.006080251 -1.0353782 -0.09406839
12 -0.52306592  0.16952226 -0.08337662  0.012836802  0.5815592 -0.15290824
13  0.94309105  0.00000000  0.00000000  0.000000000  1.4310438  0.00000000
14 -0.65815450  0.07794213  0.19129114 -0.299184150 -0.4980206  7.49272401
15 -0.59483154  1.12484901 -0.97060976  0.706754239 -0.5623145 -0.09043464
16  0.69220111  0.00000000  0.00000000  0.000000000  0.9347302  0.00000000
17 -1.98315017  0.23262378  0.13782332 -0.601600180  0.2961028  0.00000000
18 -0.10409315 -0.57548118  0.42897257 -0.327766905  0.4076703 -0.29213956
19  3.67257319  0.01868677 -0.16967480  0.801906380 -0.4627143  0.00000000
20 -0.49722851 -0.56833776 -0.06998678  1.989171623 -0.5644669 -0.07662187
21 -2.24658422 -2.79290079  0.93860932  0.009346107 -0.4592750 -0.11264279
22 -0.80781004  0.00000000  0.00000000  0.000000000  1.1609051  0.00000000
23 -1.82458945 -0.12296804 -0.04367538 -0.101201049  1.4785133 -0.11910759
24 -2.12037515  0.11919573 -0.18488561  0.152377362  0.4201025 -0.10017271
25 -2.06972202  0.00000000  0.00000000  0.000000000  1.2507436  0.00000000
26 -1.43323455 -2.22849399  0.87123913 -2.003817241  4.5100999 -0.14197043
27 -0.05648774 -3.97936452  1.10171341 -0.256687137  1.1935003  0.42194491
           3-2         3-3        3-4
1   2.29269350 -0.62123152 -0.2932966
2  -0.97828572  1.83735296 -1.1148390
3  -0.20315809 -0.26473005  0.5764339
4   0.00000000  0.00000000  0.3425668
5   0.05092255 -0.06886912 -0.4295771
6  -0.65621370  0.20179904  0.2271721
7   0.00000000  0.00000000  0.3072074
8   0.00000000  0.00000000  0.9947920
9  -0.04096061 -0.74948176  0.2527053
10  0.47663601  0.30451686 -1.6319805
11 -0.75731532  0.73740802 -1.2033224
12  0.49851695  0.59202380 -2.9366369
13  0.00000000  0.00000000  0.0000000
14 -0.65598855  0.21477180 -0.7816204
15 -0.64367774  0.62905551 -0.7910170
16  0.00000000  0.00000000  0.0000000
17  0.00000000  0.00000000  0.6884913
18 -0.56483139 -0.10630860  0.2379233
19  0.00000000  0.00000000  1.0359481
20  0.32350753  0.08907039 -0.4874280
21  5.35478511  0.03786605 -1.8982974
22  0.00000000  0.00000000  0.3271597
23 -0.69465868  0.45500827 -0.5113362
24  2.51230409  0.07088159 -1.2934926
25  0.00000000  0.00000000  0.2274350
26 -0.33493451 -0.43154092  0.8796763
27 -0.01033312 -0.40292685  0.1171232

$test
     stat df.lower p.lower df.upper    p.upper
 282.0627       NA      NA      237 0.02375896
  • 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.