Compute subject-specific transition hazards with (co-)variances

This function computes subject-specific or overall cumulative transition hazards for each of the possible transitions in the multi-state model. If requested, also the variances and covariances of the estimated cumulative transition hazards are calculated.

Usage

msfit(
  object,
  newdata,
  variance = TRUE,
  vartype = c("aalen", "greenwood"),
  trans
)

Arguments

object: A coxph object describing the fit of the multi-state model
newdata: A data frame with the same variable names as those that appear in the coxph formula. Its use is somewhat different from survfit. See Details. The argument newdata may be omitted only if the right hand side of the formula in the coxph object is ~strata(trans)
variance: A logical value indicating whether the (co-)variances of the subject-specific transition hazards should be computed. Default is TRUE
vartype: A character string specifying the type of variances to be computed (so only needed if variance=TRUE). Possible values are "aalen" or "greenwood"
trans: Transition matrix describing the states and transitions in the multi-state model. See trans in msprep for more detailed information

Value

An object of class "msfit", which is a list containing

Haz: A data frame with time, Haz, trans, containing the estimated subject-specific hazards for each of the transitions in the multi-state model
varHaz: A data frame with time, Haz, trans1, trans2 containing the variances (trans1=trans2) and covariances (trans1<trans2) of the estimated hazards. This element is only returned when variance=TRUE
trans: The transition matrix used

Details

The data frame needs to have one row for each transition in the multi-state model. An additional column strata (numeric) is needed to describe for each transition to which stratum it belongs. The name has to be strata, even if in the original coxph call another variable was used. For details refer to de Wreede, Fiocco & Putter (2010). So far, the results have been checked only for the "breslow" method of dealing with ties in coxph, so this is recommended.

References

Putter H, Fiocco M, Geskus RB (2007). Tutorial in biostatistics: Competing risks and multi-state models. Statistics in Medicine 26, 2389–2430.

Therneau TM, Grambsch PM (2000). Modeling Survival Data: Extending the Cox Model. Springer, New York.

de Wreede LC, Fiocco M, and Putter H (2010). The mstate package for estimation and prediction in non- and semi-parametric multi-state and competing risks models. Computer Methods and Programs in Biomedicine 99, 261–274.

de Wreede LC, Fiocco M, and Putter H (2011). mstate: An R Package for the Analysis of Competing Risks and Multi-State Models. Journal of Statistical Software, Volume 38, Issue 7.

Author

Hein Putter H.Putter@lumc.nl

Examples


# transition matrix for illness-death model
tmat <- trans.illdeath()
# data in wide format, for transition 1 this is dataset E1 of
# Therneau & Grambsch (2000)
tg <- data.frame(illt=c(1,1,6,6,8,9),ills=c(1,0,1,1,0,1),
        dt=c(5,1,9,7,8,12),ds=c(1,1,1,1,1,1),
        x1=c(1,1,1,0,0,0),x2=c(6:1))
# data in long format using msprep
tglong <- msprep(time=c(NA,"illt","dt"),status=c(NA,"ills","ds"),
    data=tg,keep=c("x1","x2"),trans=tmat)
# events
events(tglong)
#> $Frequencies
#>          to
#> from      healthy illness death no event total entering
#>   healthy       0       4     2        0              6
#>   illness       0       0     4        0              4
#>   death         0       0     0        6              6
#> 
#> $Proportions
#>          to
#> from        healthy   illness     death  no event
#>   healthy 0.0000000 0.6666667 0.3333333 0.0000000
#>   illness 0.0000000 0.0000000 1.0000000 0.0000000
#>   death   0.0000000 0.0000000 0.0000000 1.0000000
#> 
table(tglong$status,tglong$to,tglong$from)
#> , ,  = 1
#> 
#>    
#>     2 3
#>   0 2 4
#>   1 4 2
#> 
#> , ,  = 2
#> 
#>    
#>     2 3
#>   0 0 0
#>   1 0 4
#> 
# expanded covariates
tglong <- expand.covs(tglong,c("x1","x2"))
# Cox model with different covariate
cx <- coxph(Surv(Tstart,Tstop,status)~x1.1+x2.2+strata(trans),
  data=tglong,method="breslow")
summary(cx)
#> Call:
#> coxph(formula = Surv(Tstart, Tstop, status) ~ x1.1 + x2.2 + strata(trans), 
#>     data = tglong, method = "breslow")
#> 
#>   n= 16, number of events= 10 
#> 
#>        coef exp(coef) se(coef)     z Pr(>|z|)
#> x1.1 1.4753    4.3723   1.2557 1.175    0.240
#> x2.2 0.8571    2.3563   0.8848 0.969    0.333
#> 
#>      exp(coef) exp(-coef) lower .95 upper .95
#> x1.1     4.372     0.2287    0.3731     51.24
#> x2.2     2.356     0.4244    0.4160     13.35
#> 
#> Concordance= 0.781  (se = 0.077 )
#> Likelihood ratio test= 2.93  on 2 df,   p=0.2
#> Wald test            = 2.32  on 2 df,   p=0.3
#> Score (logrank) test = 2.86  on 2 df,   p=0.2
#> 
# new data, to check whether results are the same for transition 1 as
# those in appendix E.1 of Therneau & Grambsch (2000)
newdata <- data.frame(trans=1:3,x1.1=c(0,0,0),x2.2=c(0,1,0),strata=1:3)
msfit(cx,newdata,trans=tmat)
#> $Haz
#>    time         Haz trans
#> 1     1 0.062046887     1
#> 2     5 0.062046887     1
#> 3     6 0.333333333     1
#> 4     7 0.333333333     1
#> 5     8 0.333333333     1
#> 6     9 1.333333333     1
#> 7    12 1.333333333     1
#> 8     1 0.007971491     2
#> 9     5 0.007971491     2
#> 10    6 0.007971491     2
#> 11    7 0.007971491     2
#> 12    8 0.305921143     2
#> 13    9 0.305921143     2
#> 14   12 0.305921143     2
#> 15    1 0.000000000     3
#> 16    5 1.000000000     3
#> 17    6 1.000000000     3
#> 18    7 1.500000000     3
#> 19    8 1.500000000     3
#> 20    9 2.500000000     3
#> 21   12 3.500000000     3
#> 
#> $varHaz
#>    time       varHaz trans1 trans2
#> 1     1 0.0078708181      1      1
#> 2     5 0.0078708181      1      1
#> 3     6 0.1111111111      1      1
#> 4     7 0.1111111111      1      1
#> 5     8 0.1111111111      1      1
#> 6     9 1.1111111111      1      1
#> 7    12 1.1111111111      1      1
#> 8     1 0.0000000000      1      2
#> 9     5 0.0000000000      1      2
#> 10    6 0.0000000000      1      2
#> 11    7 0.0000000000      1      2
#> 12    8 0.0000000000      1      2
#> 13    9 0.0000000000      1      2
#> 14   12 0.0000000000      1      2
#> 15    1 0.0000000000      1      3
#> 16    5 0.0000000000      1      3
#> 17    6 0.0000000000      1      3
#> 18    7 0.0000000000      1      3
#> 19    8 0.0000000000      1      3
#> 20    9 0.0000000000      1      3
#> 21   12 0.0000000000      1      3
#> 22    1 0.0009825368      2      2
#> 23    5 0.0009825368      2      2
#> 24    6 0.0009825368      2      2
#> 25    7 0.0009825368      2      2
#> 26    8 0.1352337224      2      2
#> 27    9 0.1352337224      2      2
#> 28   12 0.1352337224      2      2
#> 29    1 0.0000000000      2      3
#> 30    5 0.0000000000      2      3
#> 31    6 0.0000000000      2      3
#> 32    7 0.0000000000      2      3
#> 33    8 0.0000000000      2      3
#> 34    9 0.0000000000      2      3
#> 35   12 0.0000000000      2      3
#> 36    1 0.0000000000      3      3
#> 37    5 1.0000000000      3      3
#> 38    6 1.0000000000      3      3
#> 39    7 1.2500000000      3      3
#> 40    8 1.2500000000      3      3
#> 41    9 2.2500000000      3      3
#> 42   12 3.2500000000      3      3
#> 
#> $trans
#>          to
#> from      healthy illness death
#>   healthy      NA       1     2
#>   illness      NA      NA     3
#>   death        NA      NA    NA
#> 
#> attr(,"class")
#> [1] "msfit"