# Predictions for continuous time, nonhomogeneous Markov multi-state models using parametric and penalised survival models. 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. ```r 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. ## Returns `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 ```r ## Not run: if (requireNamespace("deSolve")) { 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) } ## End(Not run) ```