
Last chance! 50% off unlimited learning
Sale ends in
A numerically efficient algorithm to calculate predictions from a continuous time, nonhomogeneous Markov multi-state 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=1e-10, atol=1e-10, slow=FALSE,
min.tm=1e-8,
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 markov_msm
vcov(object, ...)
# S3 method for 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 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 markov_msm_ratio
as.data.frame(x, row.names=NULL, optional=FALSE, ...)
standardise(x, ...)
# S3 method for markov_msm
standardise(x,
weights = rep(1,nrow(x$newdata)),
normalise = TRUE, ...)
# S3 method for 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 markov_msm
subset(x, subset, ...)
# S3 method for markov_msm
diff(x, y, ...)
ratio_markov_msm(x, y, ...)
# S3 method for markov_msm
rbind(..., deparse.level=1)
# S3 method for 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,...)
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:
a numeric vector with the times for the predictions
an array
for the predicted state occupancy
probabilities. The array has three dimensions: time, state, and observations.
an array
for the predicted sojourn times (or
length of stay). The array has three dimensions: time, state, and observations.
an array
for the partial derivatives of the
predicted state occupancy probabilities with respect to the model coefficients. The array has
four dimensions: time, state, coefficients, and observations.
an array
for the partial derivatives of the predicted sojourn times (or
length of stay) with respect to the model coefficients. The array has
four dimensions: time, state, coefficients, and observations.
a data.frame
with the covariates used for
the predictions
the variance-covariance matrix for the models of the transition intensities
copy of the trans
input argument
the call to the function
For debugging:
data returned from the ordinary differential equation solver. This may include more information on the predictions
For markov_msm
:
list of functions or parametric or penalised survival models. Currently
the models include
combinations of stpm2
, pstpm2
, glm
,
gam
,
survPen
or an object of class "zeroModel"
from
zeroModel
based on one of the other classes. The order in
the list matches the indexing in the trans
argument. The
functions can optionally use a t
argument for time and/or a
newdata
argument. Uncertainty in the models are incorporated into
the gradients, while uncertainty in the functions are currently not modelled.
Transition matrix describing the states and transitions
in the multi-state model. If S is the number of states in the
multi-state model, trans
should be an S x S matrix,
with (i,j)-element a positive integer if a transition from i to j
is possible in the multi-state model, NA
otherwise. In particular,
all diagonal elements should be NA
. The
integers indicating the possible transitions in the multi-state
model should be sequentially numbered, 1,...,K, with K the number
of transitions. See msprep
numerical vector for the times to evaluation the predictions. Includes the start time
data.frame
of the covariates to use in the predictions
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. "[<-"(rep(0,nrow(trans)),1,1)
).
specifies the name of the time variable. This should be set for
regression models that do not specify this (e.g. glm
) or
where the time variable is ambiguous
If there is a singularity in the observed hazard,
for example a Weibull distribution with shape < 1
has infinite
hazard at t=0
, then as a workaround, the hazard is assumed to
be a large finite number, sing.inf
, at this time. The
results should not be sensitive to the exact value assumed, but
users should make sure by adjusting this parameter in these cases.
For markov_msm
, the method used by the ordinary differential equation solver. Defaults to
Adams method ("adams"
) for non-stiff differential equations.
For splineFun
, the method jused for spline interpolation; see splinefun
.
relative error tolerance, either a
scalar or an array as long as the number of states. Passed to lsode
absolute error tolerance, either a scalar or an array as
long as the number of states. Passed to lsode
logical to show whether to use the slow R
-only
implementation. Useful for debugging. Currently needed for costs.
Minimum time used for evaluations. Avoids log(0) for some models.
a function of the form function(t)
that returns a utility for
each state at time t
for the length of stay values
a function of the form function(t)
that returns the standard
deviation for the utility for
each state at time t
for the length of stay values
logical for whether to use costs. Default: FALSE
a function of the form function(t)
that returns the cost for each transition
a function of the form function(t)
that returns the standard deviation
for the cost for each transition
a function of the form function(t)
that returns the cost per unit
time for each state
a function of the form function(t)
that returns the standard deviation
for the cost per unit time for each state
numerical value for the proportional reduction (per unit time) in the length of stay and costs
divide newdata
into blocks. Uses less memory but is slower. Reduce this number if the function call runs out of memory.
logical for whether to use spline interpolation for the transition hazards rather than the model predictions directly (default=TRUE).
logical flag for whether to keep the full output from the ordinary differential equation in the res
component (default=FALSE
).
other arguments. For markov_msm
, these are passed to the ode
solver from the
deSolve
package. For plot.markov_msm
, these arguments are passed to plot.default
For as.data.frame.markov_msm
:
add in row names to the output data-frame
(not currently used)
logical for whether to include confidence intervals. Default: TRUE
type of transformation for the confidence interval
calculation for the state occupancy probabilities. Default: log-log transformation. This is changed for
diff
and ratio_markov_msm
objects
type of transformation for the confidence interval
calculation for the length of stay calculation. Default: log transformation. This is changed for
diff
and ratio_markov_msm
objects
type of transformation for the confidence interval
calculation for the length of stay calculation. Default: log transformation. This is changed for
diff
and ratio_markov_msm
objects
valid values for the state occupancy probabilities. Default: (0,1). This is changed for
diff
and ratio_markov_msm
objects
valid values for the state occupancy probabilities. Default: (0,Inf). This is changed for
diff
and ratio_markov_msm
objects
valid values for the state occupancy probabilities. Default: (0,Inf). This is changed for
diff
and ratio_markov_msm
objects
Not currently documented
Not currently documented
For standardise.markov_msm
:
numerical vector to use in standardising the state occupancy probabilities, length of stay and costs. Default: 1 for each observation.
logical for whether to normalise the weights to 1. Default: TRUE
For plot.markov_msm
:
(currently ignored)
logical for whether to stack the plots. Default: TRUE
x-axis label
x-axis label
colours (ignored if ggplot2=TRUE
)
border colours for the polygon
(ignored if ggplot=TRUE
)
use ggplot2
alpha value for confidence bands (ggplot)
use lattice
formula for the stratification factors for the plot
For subset.markov_msm
:
expression that is evaluated on the newdata
component of the object to filter (or restrict) for the covariates used
for predictions
For transform.markov_msm
:
an object of class "markov_msm"
For rbind.markov_msm
:
not currently used
For collapse.states
:
either an index of the states to collapse or a character vector of the state names to collapse
separator to use for the collapsed state names
For zeroModel
to predict zero rates:
survival regression object to be wrapped
For hrModel
to predict rates times a hazard ratio:
hazard ratio
alternative specification for the se of the log(hazard ratio); see also ci
argument
For aftModel
to predict accelerated rates:
acceleration factor
alternative specification for the se of the log(acceleration factor); see also ci
argument
addModel
predict rates based on adding rates from different models
hazFun
provides a rate function without uncertainty:
rate function, possibly with tmvar
and/or newdata
as arguments
splineFun
predicts rates using spline interpolation:
exact times
rates as per time
rate multiplier (e.g. scale=365.25
for converting from daily rates to yearly rates)
Mark Clements
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 re-run 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.
pmatrix.fs
, probtrans
if (FALSE) {
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("Post-surgery"=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] <- ">20-50 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==">20-50 mm"))
a <- transform(a, label = "Prob(Size<=20 mm)-Prob(20mm20-50 mm"))
b <- transform(b,label="Prob(Size<=20 mm)-Prob(20mm50 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==">20-50 mm"),
subset(pred1,nodes==0 & size==">50 mm"))
e <- transform(e,label="Prob(20mm=50mm)")
f <- ratio_markov_msm(subset(pred1,nodes==0 & size==">20-50 mm"),
subset(pred1,nodes==0 & size==">50 mm"))
f <- transform(f, label = "Prob(20mm=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)
}
Run the code above in your browser using DataLab