fmsm {flexmsm}

Flexible transition intensity based models for univariate multistate processes

Description

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.

Usage

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)

Arguments

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.

Value

The function returns an object of class fmsm as described in fmsmObject.

Examples

## 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)