markov_msm  R Documentation 
A numerically efficient algorithm to calculate predictions from a continuous time, nonhomogeneous Markov multistate model. The main inputs are the models for the transition intensities, the initial values, the transition matrix and the covariate patterns. The predictions include state occupancy probabilities (possibly with discounting and utilities), length of stay and costs. Standard errors are calculated using the delta method. Includes, differences, ratios and standardisation.
markov_msm(x, trans, t = c(0,1), newdata = NULL, init=NULL, tmvar = NULL, sing.inf = 1e+10, method="adams", rtol=1e10, atol=1e10, slow=FALSE, min.tm=1e8, utility=function(t) rep(1, nrow(trans)), utility.sd=rep(0,nrow(trans)), use.costs=FALSE, transition.costs=function(t) rep(0, sum(!is.na(trans))), # per transition transition.costs.sd=rep(0,sum(!is.na(trans))), state.costs=function(t) rep(0,nrow(trans)), # per unit time state.costs.sd=rep(0,nrow(trans)), discount.rate = 0, block.size=500, spline.interpolation=FALSE, debug=FALSE, ...) ## S3 method for class 'markov_msm' vcov(object, ...) ## S3 method for class 'markov_msm' as.data.frame(x, row.names=NULL, optional=FALSE, ci=TRUE, P.conf.type="logit", L.conf.type="log", C.conf.type="log", P.range=c(0,1), L.range=c(0,Inf), C.range=c(0,Inf), state.weights=NULL, obs.weights=NULL, ...) ## S3 method for class 'markov_msm_diff' as.data.frame(x, row.names=NULL, optional=FALSE, P.conf.type="plain", L.conf.type="plain", C.conf.type="plain", P.range=c(Inf,Inf), L.range=c(Inf,Inf), C.range=c(Inf,Inf), ...) ## S3 method for class 'markov_msm_ratio' as.data.frame(x, row.names=NULL, optional=FALSE, ...) standardise(x, ...) ## S3 method for class 'markov_msm' standardise(x, weights = rep(1,nrow(x$newdata)), normalise = TRUE, ...) ## S3 method for class 'markov_msm' plot(x, y, stacked=TRUE, which=c('P','L'), xlab="Time", ylab=NULL, col=2:6, border=col, ggplot2=FALSE, lattice=FALSE, alpha=0.2, strata=NULL, ...) ## S3 method for class 'markov_msm' subset(x, subset, ...) ## S3 method for class 'markov_msm' diff(x, y, ...) ratio_markov_msm(x, y, ...) ## S3 method for class 'markov_msm' rbind(..., deparse.level=1) ## S3 method for class 'markov_msm' transform(`_data`, ...) collapse_markov_msm(object, which=NULL, sep="; ") zeroModel(object) hrModel(object,hr=1,ci=NULL,seloghr=NULL) aftModel(object,af=1,ci=NULL,selogaf=NULL) addModel(...) hazFun(f, tmvar="t", ...) splineFun(time,rate,method="natural",scale=1,...)
For markov_msm
:
x 
list of functions or parametric or penalised survival models. Currently
the models include
combinations of 
trans 
Transition matrix describing the states and transitions
in the multistate model. If S is the number of states in the
multistate model, 
t 
numerical vector for the times to evaluation the predictions. Includes the start time 
newdata 

init 
vector of the initial values with the same length as the number of states. Defaults to the first state having an
initial value of 1 (i.e. 
tmvar 
specifies the name of the time variable. This should be set for
regression models that do not specify this (e.g. 
sing.inf 
If there is a singularity in the observed hazard,
for example a Weibull distribution with 
method 
For For 
rtol 
relative error tolerance, either a
scalar or an array as long as the number of states. Passed to 
atol 
absolute error tolerance, either a scalar or an array as
long as the number of states. Passed to 
slow 
logical to show whether to use the slow 
min.tm 
Minimum time used for evaluations. Avoids log(0) for some models. 
utility 
a function of the form 
utility.sd 
a function of the form 
use.costs 
logical for whether to use costs. Default: FALSE 
transition.costs 
a function of the form 
transition.costs.sd 
a function of the form 
state.costs 
a function of the form 
state.costs.sd 
a function of the form 
discount.rate 
numerical value for the proportional reduction (per unit time) in the length of stay and costs 
block.size 
divide 
spline.interpolation 
logical for whether to use spline interpolation for the transition hazards rather than the model predictions directly (default=TRUE). 
debug 
logical flag for whether to keep the full output from the ordinary differential equation in the 
... 
other arguments. For 
For as.data.frame.markov_msm
:
row.names 
add in row names to the output dataframe 
optional 
(not currently used) 
ci 
logical for whether to include confidence intervals. Default: TRUE 
P.conf.type 
type of transformation for the confidence interval
calculation for the state occupancy probabilities. Default: loglog transformation. This is changed for

L.conf.type 
type of transformation for the confidence interval
calculation for the length of stay calculation. Default: log transformation. This is changed for

C.conf.type 
type of transformation for the confidence interval
calculation for the length of stay calculation. Default: log transformation. This is changed for

P.range 
valid values for the state occupancy probabilities. Default: (0,1). This is changed for

L.range 
valid values for the state occupancy probabilities. Default: (0,Inf). This is changed for

C.range 
valid values for the state occupancy probabilities. Default: (0,Inf). This is changed for

state.weights 
Not currently documented 
obs.weights 
Not currently documented 
For standardise.markov_msm
:
weights 
numerical vector to use in standardising the state occupancy probabilities, length of stay and costs. Default: 1 for each observation. 
normalise 
logical for whether to normalise the weights to 1. Default: TRUE 
For plot.markov_msm
:
y 
(currently ignored) 
stacked 
logical for whether to stack the plots. Default: TRUE 
xlab 
xaxis label 
ylab 
xaxis label 
col 
colours (ignored if 
border 
border colours for the 
ggplot2 
use 
alpha 
alpha value for confidence bands (ggplot) 
lattice 
use 
strata 
formula for the stratification factors for the plot 
For subset.markov_msm
:
subset 
expression that is evaluated on the 
For transform.markov_msm
:
_data 
an object of class 
For rbind.markov_msm
:
deparse.level 
not currently used 
For collapse.states
:
which 
either an index of the states to collapse or a character vector of the state names to collapse 
sep 
separator to use for the collapsed state names 
For zeroModel
to predict zero rates:
object 
survival regression object to be wrapped 
For hrModel
to predict rates times a hazard ratio:
hr 
hazard ratio 
seloghr 
alternative specification for the se of the log(hazard ratio); see also 
For aftModel
to predict accelerated rates:
af 
acceleration factor 
selogaf 
alternative specification for the se of the log(acceleration factor); see also 
addModel
predict rates based on adding rates from different models
hazFun
provides a rate function without uncertainty:
f 
rate function, possibly with 
splineFun
predicts rates using spline interpolation:
time 
exact times 
rate 
rates as per 
scale 
rate multiplier (e.g. 
The predictions are calculated using an ordinary differential equation solver. The algorithm uses a single run of the solver to calculate the state occupancy probabilities, length of stay, costs and their partial derivatives with respect to the model parameters. The predictions can also be combined to calculate differences, ratios and standardised.
The current implementation supports a list of models for each transition.
The current implementation also only allows for a vector of initial values rather than a matrix. The predictions will need to be rerun for different vectors of initial values.
For as.data.frame.markov_msm_ratio
, the data are provided in
log form, hence the default transformations and bounds are as per
as.data.frame.markov_msm_diff
, with untransformed data on the
real line.
TODO: allow for one model to predict for the different transitions.
markov_msm
returns an object of class
"markov_msm"
.
The function summary
is used to
obtain and print a summary and analysis of variance table of the
results. The generic accessor functions coef
and vcov
extract
various useful features of the value returned by markov_msm
.
An object of class "markov_msm"
is a list containing at least the
following components:
time 
a numeric vector with the times for the predictions 
P 
an 
L 
an 
Pu 
an 
Lu 
an 
newdata 
a 
vcov 
the variancecovariance matrix for the models of the transition intensities 
trans 
copy of the 
call 
the call to the function 
For debugging:
res 
data returned from the ordinary differential equation solver. This may include more information on the predictions 
Mark Clements
pmatrix.fs
, probtrans
## Not run: library(readstata13) library(mstate) library(ggplot2) library(survival) ## Two states: Initial > Final ## Note: this shows how to use markov_msm to estimate survival and risk probabilities based on ## smooth hazard models. two_states < function(model, ...) { transmat = matrix(c(NA,1,NA,NA),2,2,byrow=TRUE) rownames(transmat) < colnames(transmat) < c("Initial","Final") rstpm2::markov_msm(list(model), ..., trans = transmat) } ## Note: the first argument is the hazard model. The other arguments are arguments to the ## markov_msm function, except for the transition matrix, which is defined by the new function. death = gsm(Surv(time,status)~factor(rx), data=survival::colon, subset=(etype==2), df=3) cr = two_states(death, newdata=data.frame(rx="Obs"), t = seq(0,2500, length=301)) plot(cr,ggplot=TRUE) ## Competing risks ## Note: this shows how to adapt the markov_msm model for competing risks. competing_risks < function(listOfModels, ...) { nRisks = length(listOfModels) transmat = matrix(NA,nRisks+1,nRisks+1) transmat[1,1+(1:nRisks)] = 1:nRisks rownames(transmat) < colnames(transmat) < c("Initial",names(listOfModels)) rstpm2::markov_msm(listOfModels, ..., trans = transmat) } ## Note: The first argument for competing_risks is a list of models. Names from that list are ## used for labelling the states. The other arguments are as per the markov_msm function, ## except for the transition matrix, which is defined by the competing_risks function. recurrence = gsm(Surv(time,status)~factor(rx), data=survival::colon, subset=(etype==1), df=3) death = gsm(Surv(time,status)~factor(rx), data=survival::colon, subset=(etype==2), df=3) cr = competing_risks(list(Recurrence=recurrence,Death=death), newdata=data.frame(rx=levels(survival::colon$rx)), t = seq(0,2500, length=301)) ## Plot the probabilities for each state for three different treatment arms plot(cr, ggplot=TRUE) + facet_grid(~ rx) ## And: differences in probabilities cr_diff = diff(subset(cr,rx=="Lev+5FU"),subset(cr,rx=="Obs")) plot(cr_diff, ggplot=TRUE, stacked=FALSE) ## Extended example: Crowther and Lambert (2017) ## library(rstpm2); library(readstata13); library(ggplot2) mex.1 < read.dta13("http://fmwww.bc.edu/repec/bocode/m/multistate_example.dta") transmat < rbind("Postsurgery"=c(NA,1,2), "Relapsed"=c(NA,NA,3), "Died"=c(NA,NA,NA)) colnames(transmat) < rownames(transmat) mex.2 < transform(mex.1,osi=(osi=="deceased")+0) levels(mex.2$size)[2] < ">2050 mm" # fix typo mex < mstate::msprep(time=c(NA,"rf","os"),status=c(NA,"rfi","osi"), data=mex.2,trans=transmat,id="pid", keep=c("age","size","nodes","pr_1","hormon")) mex < transform(mex, size2=(unclass(size)==2)+0, # avoids issues with TRUE/FALSE size3=(unclass(size)==3)+0, hormon=(hormon=="yes")+0, Tstart=Tstart/12, Tstop=Tstop/12) ## c.ar < stpm2(Surv(Tstart,Tstop,status) ~ age + size2 + size3 + nodes + pr_1 + hormon, data = mex, subset=trans==1, df=3, tvc=list(size2=1,size3=1,pr_1=1)) c.ad < stpm2(Surv(Tstart, Tstop, status) ~ age + size + nodes + pr_1 + hormon, data = mex, subset=trans==2, df=1) c.rd < stpm2( Surv(Tstart,Tstop,status) ~ age + size + nodes + pr_1 + hormon, data=mex, subset=trans==3, df=3, tvc=list(pr_1=1)) ## nd < expand.grid(nodes=seq(0,20,10), size=levels(mex$size)) nd < transform(nd, age=54, pr_1=3, hormon=0, size2=(unclass(size)==2)+0, size3=(unclass(size)==3)+0) ## Predictions system.time(pred1 < rstpm2::markov_msm(list(c.ar,c.ad,c.rd), t = seq(0,15,length=301), newdata=nd, trans = transmat)) # ~2 seconds pred1 < transform(pred1, Nodes=paste("Nodes =",nodes), Size=paste("Size",size)) ## Figure 3 plot(pred1, ggplot=TRUE) + facet_grid(Nodes ~ Size) + xlab("Years since surgery") plot(pred1, ggplot=TRUE, flipped=TRUE) + facet_grid(Nodes ~ Size) + xlab("Years since surgery") plot(pred1, strata=~nodes+size, xlab="Years since surgery", lattice=TRUE) ## Figure 4 plot(subset(pred1, nodes==0 & size=="<=20 mm"), stacked=FALSE, ggplot=TRUE) + facet_grid(. ~ state) + xlab("Years since surgery") ## Figure 5 a < diff(subset(pred1,nodes==0 & size=="<=20 mm"), subset(pred1,nodes==0 & size==">2050 mm")) a < transform(a, label = "Prob(Size<=20 mm)Prob(20mm<Size<50mm)") b < ratio_markov_msm(subset(pred1,nodes==0 & size=="<=20 mm"), subset(pred1,nodes==0 & size==">2050 mm")) b < transform(b,label="Prob(Size<=20 mm)Prob(20mm<Size<50mm)") ## c < diff(subset(pred1,nodes==0 & size=="<=20 mm"), subset(pred1,nodes==0 & size==">50 mm")) c < transform(c, label = "Prob(Size<=20 mm)Prob(Size>=50mm)") d < ratio_markov_msm(subset(pred1,nodes==0 & size=="<=20 mm"), subset(pred1,nodes==0 & size==">50 mm")) d < transform(d,label= "Prob(Size<=20 mm)Prob(Size>=50mm)") ## e < diff(subset(pred1,nodes==0 & size==">2050 mm"), subset(pred1,nodes==0 & size==">50 mm")) e < transform(e,label="Prob(20mm<Size<50 mm)Prob(Size>=50mm)") f < ratio_markov_msm(subset(pred1,nodes==0 & size==">2050 mm"), subset(pred1,nodes==0 & size==">50 mm")) f < transform(f, label = "Prob(20mm<Size<50 mm)Prob(Size>=50mm)") ## combine diffs < rbind(a,c,e) ratios < rbind(b,d,f) ## Figure 5 plot(diffs, stacked=FALSE, ggplot2=TRUE) + xlab("Years since surgery") + ylim(c(0.4, 0.4)) + facet_grid(label ~ state) ## plot(ratios, stacked=FALSE, ggplot2=TRUE) + xlab("Years since surgery") + ylim(c(0, 3)) + facet_grid(label ~ state) ## Figure 6 plot(subset(pred1, nodes==0 & size=="<=20 mm"), stacked=FALSE, which="L", ggplot2=TRUE) + facet_grid(. ~ state) + xlab("Years since surgery") ## Figure 7 plot(diffs, stacked=FALSE, which="L", ggplot2=TRUE) + xlab("Years since surgery") + ylim(c(4, 4)) + facet_grid(label ~ state) plot(ratios, stacked=FALSE, which="L", ggplot2=TRUE) + xlab("Years since surgery") + ylim(c(0.1, 10)) + coord_trans(y="log10") + facet_grid(label ~ state) ## End(Not run)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.