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.47806,0.54430) 0.13235 (0.11286,0.15028)
State 2 0.24973 (0.19704,0.30803) 0.13272 (0.10308,0.16378)
State 3 0.06806 (0.04187,0.10761) 0.06938 (0.04426,0.10129)
State 4 0                         0                        
        State 3                   State 4                  
State 1 0.07304 (0.05884,0.08769) 0.28293 (0.25576,0.31534)
State 2 0.14048 (0.10504,0.17385) 0.47707 (0.42134,0.54594)
State 3 0.13737 (0.08899,0.18893) 0.72519 (0.63984,0.80008)
State 4 0                         1.00000 (1.00000,1.00000)
            # ci="bootstrap", cores = 5, B=1000)

# estimating the transition probability matrix at 1 year;
pmatrix.msm(heart.msm1, t=1, 
            ci = "norm")
        State 1                    State 2                   
State 1 0.85068 (0.836468,0.86407) 0.08651 (0.076420,0.09631)
State 2 0.16324 (0.124593,0.20686) 0.56107 (0.508152,0.60537)
State 3 0.01183 (0.006648,0.01929) 0.08796 (0.055937,0.13670)
State 4 0                          0                         
        State 3                    State 4                   
State 1 0.01269 (0.010088,0.01577) 0.05012 (0.042843,0.06015)
State 2 0.17811 (0.146160,0.21372) 0.09758 (0.076336,0.13810)
State 3 0.62929 (0.558937,0.68906) 0.27093 (0.222602,0.33376)
State 4 0                          1.00000 (1.000000,1.00000)
  • 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.47709,0.54644) 0.13235 (0.11453,0.15063)
State 2 0.24973 (0.19527,0.30281) 0.13272 (0.10565,0.16457)
State 3 0.06806 (0.04069,0.10646) 0.06938 (0.04295,0.10392)
State 4 0                         0                        
        State 3                   State 4                  
State 1 0.07304 (0.05823,0.08653) 0.28293 (0.25587,0.31879)
State 2 0.14048 (0.10500,0.17487) 0.47707 (0.42621,0.54598)
State 3 0.13737 (0.08864,0.19472) 0.72519 (0.63975,0.80398)
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.913155 1.832191 0.8400134    0
97.5% 11.213078 9.724424 4.6446713    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.613802 2.611640 3.009151    0
97.5% 10.347013 9.454684 6.679838    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.61   1   9   2  0.63   0
2   [3.04,5.97)        1 [-2.67,-1.23)  53   5   1  0.00   1   4   5  0.32   0
3  [5.97,19.46]        1 [-2.67,-1.23)  48   4   2  0.84   4   7   3  0.63   0
4      [0,3.04)        2 [-2.67,-1.23)  66  12   0  6.86   0   0   0  0.75   0
5   [3.04,5.97)        2 [-2.67,-1.23)  67   3   0  0.59   3   2   3  0.34   1
6  [5.97,19.46]        2 [-2.67,-1.23)  53   8   1  0.98   3   5   1  0.56   0
7      [0,3.04)        3 [-2.67,-1.23)  72   6   1  7.53   0   0   0  0.62   0
8   [3.04,5.97)        3 [-2.67,-1.23)  45   7   0  6.41   1   3   0  0.34   0
9  [5.97,19.46]        3 [-2.67,-1.23)  40   9   2 13.18   0   5   1  7.81   1
10     [0,3.04)        1 [-1.23,-1.18)  55   8   3  7.22   3   2   1  0.34   0
11  [3.04,5.97)        1 [-1.23,-1.18)  46   5   0  1.00   4   8   3  0.00   0
12 [5.97,19.46]        1 [-1.23,-1.18)  22   9   0  0.94   3   8   3  1.05   0
13     [0,3.04)        2 [-1.23,-1.18)  64   4   5  9.88   0   0   0  0.34   0
14  [3.04,5.97)        2 [-1.23,-1.18)  45   5   2  2.93   2   5   1  0.39   1
15 [5.97,19.46]        2 [-1.23,-1.18)  27   9   2  1.43   5   7   5  1.22   0
16     [0,3.04)        3 [-1.23,-1.18)  57   7   2  8.90   0   0   0  0.32   0
17  [3.04,5.97)        3 [-1.23,-1.18)  71  13   3  8.07   1   1   0  0.61   0
18 [5.97,19.46]        3 [-1.23,-1.18)  44  12   2 17.63   4   9   4  8.73   0
19     [0,3.04)        1 [-1.18,-1.18]  34   4   4  6.26   2   5   3  0.29   0
20  [3.04,5.97)        1 [-1.18,-1.18]  41   5   0  2.33   2  11   6  0.96   0
21 [5.97,19.46]        1 [-1.18,-1.18]  36  10   0  0.27   0  13   3  0.66   0
22     [0,3.04)        2 [-1.18,-1.18]  58   9   0  7.17   0   0   0  0.33   0
23  [3.04,5.97)        2 [-1.18,-1.18]  50   8   0  3.21   2   6   2  2.14   0
24 [5.97,19.46]        2 [-1.18,-1.18]  38   5   0  0.36   4   9   4  1.41   0
25     [0,3.04)        3 [-1.18,-1.18]  68  11   4  6.57   0   0   0  0.38   0
26  [3.04,5.97)        3 [-1.18,-1.18]  63  10   4  7.46   0   4   0  4.90   0
27 [5.97,19.46]        3 [-1.18,-1.18]  33  13   1 14.37   1  11   4 11.93   1
   3-2 3-3   3-4
1    1   0  0.29
2    0   7  1.34
3    0   1  1.33
4    0   0  0.31
5    1   3  1.50
6    0   4  2.02
7    0   0  0.40
8    0   0  2.16
9    1   2  8.65
10   1   6  0.00
11   0  11  0.75
12   2  16  0.82
13   0   0  0.00
14   0   6  1.15
15   0   8  1.26
16   0   0  0.00
17   0   0  1.10
18   0   3  8.92
19   0   0  0.33
20   1   7  1.39
21   3  11  0.86
22   0   0  0.37
23   0   9  1.95
24   2   8  1.09
25   0   0  0.30
26   0   1  3.66
27   1   4 13.05

$Expected
           Time Interval           Cov      1-1       1-2       1-3       1-4
1      [0,3.04)        1 [-2.67,-1.23) 73.84465  6.611609 1.4650656  2.688680
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.78366  3.692716 0.6794497  1.684170
4      [0,3.04)        2 [-2.67,-1.23) 67.40294  8.783546 3.0416521  5.631859
5   [3.04,5.97)        2 [-2.67,-1.23) 58.72336  6.588486 1.9992342  3.278925
6  [5.97,19.46]        2 [-2.67,-1.23) 55.67237  4.748332 1.0172959  1.542000
7      [0,3.04)        3 [-2.67,-1.23) 66.84182  9.191939 3.3116859  7.184558
8   [3.04,5.97)        3 [-2.67,-1.23) 43.21207  6.058211 2.3521786  6.787544
9  [5.97,19.46]        3 [-2.67,-1.23) 42.99168  6.364300 2.6927406 12.131277
10     [0,3.04)        1 [-1.23,-1.18) 60.39190  7.489158 1.4284131  3.910525
11  [3.04,5.97)        1 [-1.23,-1.18) 44.15432  4.812956 0.7338816  2.298838
12 [5.97,19.46]        1 [-1.23,-1.18) 27.20531  2.939734 0.4431583  1.351794
13     [0,3.04)        2 [-1.23,-1.18) 60.24108 11.314736 3.4521539  7.872025
14  [3.04,5.97)        2 [-1.23,-1.18) 42.85029  6.444677 1.5198693  4.115161
15 [5.97,19.46]        2 [-1.23,-1.18) 33.12259  3.849742 0.6317488  1.825923
16     [0,3.04)        3 [-1.23,-1.18) 52.75290 10.508186 3.3694131  8.269502
17  [3.04,5.97)        3 [-1.23,-1.18) 64.26253 13.479621 4.5438649 12.783983
18 [5.97,19.46]        3 [-1.23,-1.18) 43.44934  9.659199 3.9544514 18.567008
19     [0,3.04)        1 [-1.18,-1.18] 38.99731  5.090023 0.9338956  3.238767
20  [3.04,5.97)        1 [-1.18,-1.18] 40.13725  4.787514 0.6802960  2.724939
21 [5.97,19.46]        1 [-1.18,-1.18] 38.46691  4.546709 0.6326658  2.623717
22     [0,3.04)        2 [-1.18,-1.18] 51.96588 10.473423 2.9701918  8.760508
23  [3.04,5.97)        2 [-1.18,-1.18] 45.42068  7.865370 1.9054666  6.018486
24 [5.97,19.46]        2 [-1.18,-1.18] 35.54611  4.470896 0.6659924  2.677002
25     [0,3.04)        3 [-1.18,-1.18] 61.39458 12.881524 3.8841401 11.409756
26  [3.04,5.97)        3 [-1.18,-1.18] 56.92452 12.169078 3.8849826 11.481424
27 [5.97,19.46]        3 [-1.18,-1.18] 36.04633  7.960627 3.0444353 14.318612
         2-1        2-2       2-3        2-4        3-1       3-2        3-3
1  1.6217318  7.3916597 2.6828036 0.93380496 0.02388242 0.2738355  0.6531386
2  1.4963105  5.6278956 2.2822973 0.91349653 0.10577836 0.9943227  4.3305491
3  1.8045130  8.4471814 3.3790237 0.99928190 0.02355670 0.1996134  1.2407996
4  0.0000000  0.0000000 0.0000000 0.12706398 0.00000000 0.0000000  0.0000000
5  1.6094367  3.7213255 1.7361266 1.27311114 0.10910916 0.8510740  3.3263080
6  1.4157981  5.1900820 2.1458167 0.80830325 0.07645147 0.6638986  3.5226937
7  0.0000000  0.0000000 0.0000000 0.11807201 0.00000000 0.0000000  0.0000000
8  0.9034993  1.4525869 1.0803474 0.90356647 0.00000000 0.0000000  0.0000000
9  2.3081045  3.1436586 2.5371156 5.82112133 0.41010551 1.1669965  3.6163494
10 1.0394339  3.7270435 1.1381695 0.43535310 0.06805836 0.5079541  4.7920071
11 2.3698262  8.7268480 2.8679476 1.03537822 0.09185500 0.7400871  8.4562218
12 2.3898869  8.7317095 2.8945760 1.03382770 0.15430057 1.2317961 13.3967014
13 0.0000000  0.0000000 0.0000000 0.05519433 0.00000000 0.0000000  0.0000000
14 1.6608662  4.1375350 1.7117670 0.87983169 0.10991421 0.6598304  5.1405026
15 3.1584369 10.1813246 3.4792824 1.40095609 0.09213857 0.6553626  6.3129841
16 0.0000000  0.0000000 0.0000000 0.05627295 0.00000000 0.0000000  0.0000000
17 0.6586024  0.7669832 0.6032918 0.58112267 0.00000000 0.0000000  0.0000000
18 5.7459091  7.2340379 5.2492341 7.50081896 0.28277844 0.5430587  3.2201001
19 1.8236509  5.9220295 1.7885443 0.75577525 0.00000000 0.0000000  0.0000000
20 3.3849929 11.7938854 3.4076863 1.37343534 0.07475796 0.5681223  6.7298155
21 2.7689075  9.9015865 2.8681902 1.12131577 0.11365492 0.8663711 10.7911576
22 0.0000000  0.0000000 0.0000000 0.05416649 0.00000000 0.0000000  0.0000000
23 2.4657366  5.9991023 2.4113981 1.26376295 0.12036748 0.6961683  7.3488417
24 3.3779634 10.2834977 3.3047998 1.44373907 0.10221020 0.7049259  7.7066729
25 0.0000000  0.0000000 0.0000000 0.06635919 0.00000000 0.0000000  0.0000000
26 2.2479783  2.6843499 2.0271151 1.94055672 0.14946151 0.3533013  1.9429328
27 5.8840569  8.1843325 5.1890587 8.67255196 0.54039054 0.9999186  5.3900174
           3-4
1   0.33914350
2   2.90934980
3   0.86603033
4   0.12560378
5   2.21350888
6   1.75695616
7   0.16569966
8   1.18542030
9   7.45654852
10  1.63198051
11  2.46183606
12  4.03720196
13  0.00000000
14  2.23975282
15  2.19951479
16  0.00000000
17  0.48232826
18  7.87406271
19  0.08200438
20  2.01730424
21  3.08881637
22  0.14507240
23  2.78462246
24  2.57619102
25  0.12461372
26  2.21430444
27 12.11967350

$`Deviance*sign(O-E)`
           Time Interval           Cov         1-1          1-2           1-3
1      [0,3.04)        1 [-2.67,-1.23) -0.14100261 -1.973455552  8.5376057031
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.02079696  0.026450788  2.5678855439
4      [0,3.04)        2 [-2.67,-1.23) -0.06807682  1.188451007 -3.0416520503
5   [3.04,5.97)        2 [-2.67,-1.23)  1.17314122 -1.954577321 -1.9992342334
6  [5.97,19.46]        2 [-2.67,-1.23) -0.13693093  2.229095502 -0.0004490655
7      [0,3.04)        3 [-2.67,-1.23)  0.43403088 -1.110283150 -1.6137894242
8   [3.04,5.97)        3 [-2.67,-1.23)  0.08253939  0.147302435 -2.3521785857
9  [5.97,19.46]        3 [-2.67,-1.23) -0.22164288  1.100233267 -0.1808394716
10     [0,3.04)        1 [-1.23,-1.18) -0.51914520  0.041603806  1.7353020679
11  [3.04,5.97)        1 [-1.23,-1.18)  0.09289234  0.008911365 -0.7338815869
12 [5.97,19.46]        1 [-1.23,-1.18) -1.00638357 12.512268215 -0.4431583473
13     [0,3.04)        2 [-1.23,-1.18)  0.28919977 -4.730006758  0.6998321858
14  [3.04,5.97)        2 [-1.23,-1.18)  0.13435604 -0.326648511  0.1541635264
15 [5.97,19.46]        2 [-1.23,-1.18) -1.14547217  6.905014348  2.9671907361
16     [0,3.04)        3 [-1.23,-1.18)  0.40832109 -1.176163021 -0.5578518084
17  [3.04,5.97)        3 [-1.23,-1.18)  0.72610874 -0.020079248 -0.5251789411
18 [5.97,19.46]        3 [-1.23,-1.18)  0.03274564  0.578656500 -0.9676384027
19     [0,3.04)        1 [-1.18,-1.18] -0.69739998 -0.239856214 10.1023268543
20  [3.04,5.97)        1 [-1.18,-1.18]  0.04251959  0.012332627 -0.6802960483
21 [5.97,19.46]        1 [-1.18,-1.18] -0.16114859  6.542779957 -0.6326658490
22     [0,3.04)        2 [-1.18,-1.18]  0.76007982 -0.214442635 -2.9701918250
23  [3.04,5.97)        2 [-1.18,-1.18]  0.49342348  0.006880599 -1.9054665773
24 [5.97,19.46]        2 [-1.18,-1.18]  0.17415210  0.063332380 -0.6659924207
25     [0,3.04)        3 [-1.18,-1.18]  0.75364871 -0.280309976  0.0059891104
26  [3.04,5.97)        3 [-1.18,-1.18]  0.67051988 -0.389283629  0.0049550627
27 [5.97,19.46]        3 [-1.18,-1.18] -0.26555323  3.199002195 -1.3734064472
           1-4         2-1         2-2          2-3        2-4         3-1
1   4.18689379 -0.23982674  0.37121096 -0.178191305 -0.4324353 -0.02388242
2  -1.64074673 -0.16625598 -0.47607368  3.253355682 -0.5955006 -0.10577836
3  -0.72726343  2.68145297 -0.25858723 -0.049224892 -0.5071787 -0.02355670
4   0.91714509  0.00000000  0.00000000  0.000000000  3.0577498  0.00000000
5  -2.32988135  1.21984533 -0.79934993  0.934259998 -0.8367365  7.89700076
6  -0.63439856  1.78406675 -0.02506425 -0.614392003 -0.4518935 -0.07645147
7   0.41817800  0.00000000  0.00000000  0.000000000  2.1849535  0.00000000
8  -0.08455965  0.02722631  1.69324079 -1.080347355 -0.5286321  0.00000000
9   0.18511115 -2.30810448  1.17606883 -0.935884846  0.8355861  0.91914746
10  3.63447496  3.74399328 -0.80551128 -0.022973005 -0.4560914 -0.06805836
11 -1.10170961  1.12137613 -0.06053824  0.006080251 -1.0353782 -0.09185500
12 -0.55169244  0.16786123 -0.08521052  0.012972816  0.6008827 -0.15430057
13  1.05711819  0.00000000  0.00000000  0.000000000  1.4700198  0.00000000
14 -0.72954277  0.07827925  0.19491679 -0.298594252 -0.5059438  7.37648422
15 -0.58688308  1.09522065 -1.00722041  0.684427655 -0.5850340 -0.09213857
16  0.61325827  0.00000000  0.00000000  0.000000000  1.2514329  0.00000000
17 -1.88377210  0.22915005  0.13551305 -0.603291768  0.2933288  0.00000000
18 -0.15150741 -0.54336515  0.46387055 -0.307877577  0.3603797 -0.28277844
19  3.86056375  0.02139861 -0.15034112  0.833712294 -0.5301192  0.00000000
20 -0.49961235 -0.56879238 -0.07039751  1.991458200 -0.5740712 -0.07475796
21 -2.18460705 -2.76890754  0.99439910  0.011193797 -0.5325848 -0.11365492
22 -0.75941826  0.00000000  0.00000000  0.000000000  1.4048502  0.00000000
23 -1.62752132 -0.11162276  0.04637557 -0.092865750  1.3714151 -0.12036748
24 -2.08906037  0.12512356 -0.18017499  0.159776096 -0.5060300 -0.10221020
25 -2.40079946  0.00000000  0.00000000  0.000000000  1.4952522  0.00000000
26 -1.56977353 -2.24797828  0.80031563 -2.027115062  4.8378721 -0.14946151
27  0.04717988 -4.05452227  0.99718456 -0.280206692  1.3434321  0.43741558
            3-2         3-3        3-4
1   2.168692922 -0.65313857 -0.2967815
2  -0.994322699  1.78317227 -1.0528670
3  -0.199613367 -0.31326745  0.6696241
4   0.000000000  0.00000000  0.2707308
5   0.058237597 -0.07681075 -0.4846279
6  -0.663898645  0.19452767  0.3079229
7   0.000000000  0.00000000  0.3317383
8   0.000000000  0.00000000  0.9022031
9  -0.034297243 -0.73956430  0.2360650
10  0.476636013  0.30451686 -1.6319805
11 -0.740087089  0.81201126 -1.3656669
12  0.482995098  0.53951510 -2.6989439
13  0.000000000  0.00000000  0.0000000
14 -0.659830391  0.19246164 -0.6965202
15 -0.655362564  0.53658931 -0.6330480
16  0.000000000  0.00000000  0.0000000
17  0.000000000  0.00000000  0.8202676
18 -0.543058726 -0.06630564  0.1780446
19  0.000000000  0.00000000  0.7805423
20  0.349176452  0.09408897 -0.5599713
21  5.286120009  0.03386984 -1.7826398
22  0.000000000  0.00000000  0.3488173
23 -0.696168344  0.46503298 -0.5445103
24  2.409456442  0.05543851 -1.0644602
25  0.000000000  0.00000000  0.2487815
26 -0.353301300 -0.50290319  1.0117628
27  0.009851226 -0.39086453  0.1121574

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