Plot observed and expected prevalences for a multi-state model

Plot observed and expected state prevalences from a fitted multi-state model. The function can also compute a rough diagnostic for where the data depart from the estimated Markov model.

Usage

prevplot(
  x,
  prev.obj,
  exacttimes = TRUE,
  M = FALSE,
  ci = FALSE,
  print_plot = TRUE,
  verbosity = getOption("msmtools.verbosity", "quiet")
)

Arguments

x: A fitted msm model object.
prev.obj: A list computed by msm::prevalence.msm(). It may include confidence intervals; prevplot() adapts automatically.
exacttimes: If TRUE (default), transition times are known and exact. This should match the value used when fitting the model with msm.
M: If TRUE, a rough indicator of deviance from the model is computed (see Details). Default is FALSE.
ci: If TRUE, confidence intervals are plotted when available. Default is FALSE.
print_plot: If TRUE (default), the plot is printed before being returned. If FALSE, the plot is returned without printing.
verbosity: Controls informational output. Use "quiet" to suppress status messages and "summary" or "progress" for high-level messages.

Value

When M = FALSE, a gg/ggplot object with observed and expected prevalences is returned. When M = TRUE, a patchwork object is returned with the prevalence plot and the deviance M plot.

The returned object also carries a $prevalence field with the long-format data.table used to build the plot. It always includes time, state, obs, and hat; it also includes lwr and upr when ci = TRUE, and M when M = TRUE. Access it directly:

p <- prevplot(model, prev_obj)
p$prevalence

print_plot only controls whether the plot is printed as a side effect. Returned objects are unchanged: use print_plot = FALSE to create the plot silently.

Details

When M = TRUE, a rough indicator of the deviance from the Markov model is computed according to Titman and Sharples (2008). A comparison at a given time t_i of a subject k in the state s between observed counts O_is and expected counts E_is is built as M_is = (O_is - E_is)^2 / E_is.

The deviance M plot is returned together with the standard prevalence plot in the second row. This layout is fixed.

When M = TRUE, the combined layout is built with patchwork, which is an optional dependency of msmtools. Install it with install.packages("patchwork") if it is not already available; prevplot() raises an informative error otherwise. The default M = FALSE path has no such requirement.

References

Titman, A. and Sharples, L.D. (2010). Model diagnostics for multi-state models, Statistical Methods in Medical Research, 19, 621-651.

Titman, A. and Sharples, L.D. (2008). A general goodness-of-fit test for Markov and hidden Markov models, Statistics in Medicine, 27, 2177-2195.

Gentleman RC, Lawless JF, Lindsey JC, Yan P. (1994). Multi-state Markov models for analysing incomplete disease data with illustrations for HIV disease. Statistics in Medicine, 13:805-821.

Jackson, C.H. (2011). Multi-State Models for Panel Data: The msm Package for R. Journal of Statistical Software, 38(8), 1-29. https://www.jstatsoft.org/v38/i08/.

Author

Francesco Grossetti francesco.grossetti@unibocconi.it.

Examples

if (FALSE) { # interactive()
data(hosp)

# augmenting the data
hosp_augmented = augment(data = hosp, data_key = subj, n_events = adm_number,
                          pattern = label_3, t_start = dateIN, t_end = dateOUT,
                          t_cens = dateCENS)

# let's define the initial transition matrix for our model
Qmat = matrix(data = 0, nrow = 3, ncol = 3, byrow = TRUE)
Qmat[1, 1:3] = 1
Qmat[2, 1:3] = 1
colnames(Qmat) = c('IN', 'OUT', 'DEAD')
rownames(Qmat) = c('IN', 'OUT', 'DEAD')

# fitting the model using
# gender and age as covariates
library(msm)
msm_model = msm(status_num ~ augmented_int, subject = subj,
                 data = hosp_augmented, covariates = ~ gender + age,
                 exacttimes = TRUE, gen.inits = TRUE, qmatrix = Qmat,
                 method = 'BFGS', control = list(fnscale = 6e+05, trace = 0,
                 REPORT = 1, maxit = 10000))

# defining the times at which compute the prevalences
t_min = min(hosp_augmented$augmented_int)
t_max = max(hosp_augmented$augmented_int)
steps = 100L

# computing prevalences
prev = prevalence.msm(msm_model, covariates = 'mean', ci = 'normal',
                       times = seq(t_min, t_max, steps))

# and plotting them using prevplot()
gof = prevplot(x = msm_model, prev.obj = prev, ci = TRUE, M = TRUE)
}