Flexible transition intensity based models for univariate multistate processes
Fits a flexible multistate survival model. Any type of process is supported, including both forward and backward transitions, and must be specified by providing a list of equations, one for each transition intensity allowed. Any type of observation scheme is allowed: the process can be observed in continuous time, intermittently at fixed times, there can be an absorbing state as well as censored states. Virtually any type of covariate effects are supported and can be specified by means of splines, with the same syntax used to specify Generalised Additive Models (GAMs) in R.
fmsm(formula, data, id, state, death, pmethod = 'eigendecomp', aggregate = TRUE, params.0 = NULL, sp.0 = NULL, constraint = NULL, sp.method = 'perf', iterlimsp = 50, Q.diagnostics = TRUE, fit = TRUE, iterlim = 100, tolsp = 1e-7, tolsp.EFS = 0.1, parallel = FALSE, no_cores = 2, cens.state = NULL, living.exact = NULL, verbose = FALSE, justComp = NULL, approxHess = FALSE)
formula: Model specification for the transition intensities.
data: Dataset.
id: Name of the variable in the dataset representing the unique code associated with each patient.
state: Name of the variable in the dataset representing the state occupied by the patient at the given time.
death: TRUE if the last state is an absorbing state, FALSE otherwise.
pmethod: Which method should be used for the computation of the transition probability matrix. Available options are
'eigendecomp' (default): this method is based on the eigendecomposition of the transition intensity matrix (from Kalbfleisch & Lawless 1985);'analytic': uses analytic expressions of the transition probabilities, obtained by solving the Kolmogorov forward differential equation, only implemented for IDMs for now;'scaling&squaring': this is the scaling and squaring method implemented as proposed in Fung (2004).This is inefficient, so its use is not recommended. Can be used to investigate convergence errors.aggregate: Whether or not data should be aggregated (this slightly improves efficiency as redundancies in the data are eliminated). The default is TRUE.
params.0: Starting values for the model parameters. Defaults to NULL, i.e. they are computed internally.
sp.0: Starting values for the smoothing parameters. Defaults to NULL, i.e. they are computed internally.
constraint: A list containing the constraints to be applied to the model parameters. For example, assuming a process with three transitions, constraint = list(x1 = c(1,1,1), x2 = c(1,2,2)) means that the effect of covariate x1 is constrained to be equal across the tree transitions and that the effects of x2 on the second and third transitions are constrained to be equal, but the effect on the first transition is left unconstrained.
sp.method: Method to be used for smoothing parameter estimation. The default is magic, the automatic multiple smoothing parameter selection algorithm. Alternatively, efs can be used for the Fellner-Schall method. To suppress the smoothing parameter estimation set this to NULL.
iterlimsp: Maximum allowed iterations for smoothing parameter estimation.
Q.diagnostics: If TRUE, diagnostics information on the Q matrix are saved. The default TRUE.
fit: If FALSE, fitting is not carried. May be useful to extract model setup quantities.
iterlim: Maximum allowed iterations for trust region algorithm.
tolsp: Convergence criterion used in magic based smoothing parameter estimation.
tolsp.EFS: Convergence criterion used in efs based smoothing parameter estimation.
parallel: If TRUE parallel computing is used during estimation. This can only be used by Windows users for now.
no_cores: Number of cores used if parallel computing chosen. The default is 2. If NULL, all available cores are used.
cens.state: Code used in the dataset to indicate the censored states.
living.exact: Name of the variable in the dataset indicating whether an observation is exactly observed or not.
verbose: If TRUE, prints the convergence criterion obtained at each iteration of the full algorithm. The default is FALSE.
justComp: Can be c('lik', 'grad', 'hess') (or a subvector of this) to compute the log-likelihood, gradient and Hessian at the starting parameter value, without carrying out fitting. Defaults to NULL.
approxHess: If TRUE an approximation of the Hessian based on the outer product of the gradient vector is computed. The default is FALSE.
The function returns an object of class fmsm as described in fmsmObject.
## Not run: ################################################## # MULTISTATE SURVIVAL MODELLING with CAV DATA #### ################################################## library(flexmsm) # Import data Data <- IDM_cav # MODEL SPECIFICATION #### formula <- list(years ~ s(years, bs = 'cr', k = 10) + dage + pdiag, # 1-2 years ~ s(years, bs = 'cr', k = 10) + dage + pdiag, # 1-3 0, # 2-1 years ~ s(years, bs = 'cr', k = 10) + dage + pdiag, # 2-3 0, # 3-1 0 # 3-2 ) # Counts of pairs of consecutive states observed (C = counts, T = times) counts.CT <- state.pairs.CT(formula = formula, data = Data, time = 'years', state = 'state', id = 'PTNUM') counts.CT$counts # MODEL FITTING ### # NOTE *** # Takes about 18 minutes on a machine with Windows 10, # Intel 2.20 GHz core, 16 GB of RAM and 8 cores, using all cores. # The default is to use 2 cores, this takes about 26 minutes. # To use all available cores on your device input no_cores = NULL. # **** fmsm.out <- fmsm(formula = formula, data = Data, id = 'PTNUM', state = 'state', death = TRUE, fit = TRUE, parallel = TRUE, no_cores = 2, pmethod = 'analytic') print(fmsm.out) AIC(fmsm.out) BIC(fmsm.out) # FITTING SUMMARY #### summary(fmsm.out) conv.check(fmsm.out) #################### # VISUALISATION #### #################### # PLOT THE SMOOTHS OF TIME FOR EACH TRANSITION #### # par(mfrow = c(1,3)) plot(fmsm.out) # Consider a patient with: dage.pred <- 16 # - 16 year old donor pdiag.pred <- 0 # - IDC as principal diagnosis start.pred <- 0 # - start observation at time t = 0 stop.pred <- 15 # - t = 15 years for time horizon n.pred <- 21 # - 21 time points no.state.pred <- -13 # - (because we don't need this, so anything is fine) newdata <- data.frame(PTNUM = rep(1, n.pred), years = seq(start.pred, stop.pred, length.out = n.pred), state = rep(no.state.pred, n.pred), dage = rep(dage.pred, n.pred), pdiag = rep(pdiag.pred, n.pred)) # ESTMATED TRANSITION INTENSITIES #### # Plot of estimated transition intensities # par(mfrow = c(1,3)) Q.hat <- Q.pred(fmsm.out, newdata = newdata, get.CI = TRUE, plot.Q = TRUE, rug = TRUE, ylim = c(0, 1.5)) # Estimated transition intensity matrix at, e.g., t = 0 round(Q.hat$Q.hist[,,1], 3) # ESTMATED TRANSITION PROBABILITIES #### # Plot of estimated transition probabilities # par(mfrow = c(2,3)) P.hat <- P.pred(fmsm.out, newdata = newdata, get.CI = TRUE, plot.P = TRUE, rug = TRUE) # Estimated 15 year transition probability matrix round(P.hat$P.pred, 3) # e.g., there is a 6.2% chance of observing CAV onset 15 years after transplant ## End(Not run)
Useful links