markov_msm: Predictions for continuous time, nonhomogeneous Markov...
In mclements/rstpm2: Smooth Survival Models, Including Generalized Survival Models
markov_msm R Documentation
Predictions for continuous time, nonhomogeneous Markov multi-state
models using parametric and penalised survival models.
Description
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.
Usage
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 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,...)
Arguments
For markov_msm:
x
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.
trans
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
t
numerical vector for the times to evaluation the predictions. Includes
the start time
newdata
data.frame of the covariates to use in the predictions
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. "[<-"(rep(0,nrow(trans)),1,1)).
tmvar
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
sing.inf
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.
method
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.
rtol
relative error tolerance, either a
scalar or an array as long as the number of states. Passed to lsode
atol
absolute error tolerance, either a scalar or an array as
long as the number of states. Passed to lsode
slow
logical to show whether to use the slow R-only
implementation. Useful for debugging. Currently needed for costs.
min.tm
Minimum time used for evaluations. Avoids log(0) for some models.
utility
a function of the form function(t) that returns a utility for
each state at time t for the length of stay values
utility.sd
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
use.costs
logical for whether to use costs. Default: FALSE
transition.costs
a function of the form function(t) that returns the cost for each transition
transition.costs.sd
a function of the form function(t) that returns the standard deviation
for the cost for each transition
state.costs
a function of the form function(t) that returns the cost per unit
time for each state
state.costs.sd
a function of the form function(t) that returns the standard deviation
for the cost per unit time for each state
discount.rate
numerical value for the proportional reduction (per unit time) in the length of stay
and costs
block.size
divide newdata into blocks. Uses less memory but is slower. Reduce this number if the function call runs out of memory.
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 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:
row.names
add in row names to the output data-frame
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: log-log transformation. This is changed for
diff and ratio_markov_msm objects
L.conf.type
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
C.conf.type
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
P.range
valid values for the state occupancy probabilities. Default: (0,1). This is changed for
diff and ratio_markov_msm objects
L.range
valid values for the state occupancy probabilities. Default: (0,Inf). This is changed for
diff and ratio_markov_msm objects
C.range
valid values for the state occupancy probabilities. Default: (0,Inf). This is changed for
diff and ratio_markov_msm objects
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
x-axis label
ylab
x-axis label
col
colours (ignored if ggplot2=TRUE)
border
border colours for the polygon (ignored if ggplot=TRUE)
ggplot2
use ggplot2
alpha
alpha value for confidence bands (ggplot)
lattice
use lattice
strata
formula for the stratification factors for the plot
For subset.markov_msm:
subset
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:
_data
an object of class "markov_msm"
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 ci argument
For aftModel to predict accelerated rates:
af
acceleration factor
selogaf
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:
f
rate function, possibly with tmvar and/or newdata as arguments
splineFun predicts rates using spline interpolation:
time
exact times
rate
rates as per time
scale
rate multiplier (e.g. scale=365.25 for converting from daily rates to yearly rates)
Details
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.
Value
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 array for the predicted state occupancy
probabilities. The array has three dimensions: time, state, and observations.
L
an array for the predicted sojourn times (or
length of stay). The array has three dimensions: time, state, and observations.
Pu
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.
Lu
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.
newdata
a data.frame with the covariates used for
the predictions
vcov
the variance-covariance matrix for the models of the
transition intensities
trans
copy of the trans input argument
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
Author(s)
Mark Clements
See Also
pmatrix.fs, probtrans
Examples
## 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("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(20mm<Size<50mm)")
b <- ratio_markov_msm(subset(pred1,nodes==0 & size=="<=20 mm"),
subset(pred1,nodes==0 & size==">20-50 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==">20-50 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==">20-50 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)
mclements/rstpm2 documentation built on Feb. 5, 2023, 8 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
| markov_msm | R Documentation |
Predictions for continuous time, nonhomogeneous Markov multi-state models using parametric and penalised survival models.
Description
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.
Usage
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 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,...)
Arguments
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 multi-state model. If S is the number of states in the
multi-state 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 data-frame |
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: log-log 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 |
x-axis label |
ylab |
x-axis 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. |
Details
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.
Value
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 variance-covariance 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 |
Author(s)
Mark Clements
See Also
pmatrix.fs, probtrans
Examples
## 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("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(20mm<Size<50mm)")
b <- ratio_markov_msm(subset(pred1,nodes==0 & size=="<=20 mm"),
subset(pred1,nodes==0 & size==">20-50 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==">20-50 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==">20-50 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)
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.