Nothing
R/prevalence.R
In msm: Multi-State Markov and Hidden Markov Models in Continuous Time
Defines functions
plot.prevalence.msm
print.msm.prevalence
prevalence.msm
expected.msm
observed.msm
get.covhist
Documented in
plot.prevalence.msm
prevalence.msm
## Enumerate the distinct subject covariate histories in the data
## Works for time homogeneous and inhomogeneous models
## Used to calculate expected prevalences for a population with those covariates
get.covhist <- function(x, subset=NULL) {
## Keep only times where the covariate changes, or first or last obs
mf <- x$data$mf
if (x$qcmodel$ncovs > 0) {
if (!is.null(subset)) {
subs <- mf$"(subject)" %in% subset
mf <- mf[subs,,drop=FALSE]
}
subj <- match(mf$"(subject)", unique(mf$"(subject)"))
n <- length(subj)
apaste <- do.call("paste", mf[,attr(mf,"covnames"),drop=FALSE])
first <- !duplicated(subj); last <- rev(!duplicated(rev(subj)))
keep <- (c(0, apaste[1:(n-1)]) != apaste) | first | last
## Keep and tabulate unique covariate series
covseries <- split(apaste[keep], subj[keep]) # as a list of char vectors
covseries <- sapply(covseries, paste, collapse=" , ") # as one char vector, one series per pt.
## also need p matrices for different times as well as different covs.
## but only interested in cov change times if there's more than one
## transition (at least one times change point)
change.times <- mf$"(time)"; change.times[first] <- change.times[last] <- 0
change.times <- split(change.times[keep & (!(first|last))], subj[keep & (!(first|last))])
change.times <- sapply(change.times, paste, collapse= " , ")
covseries.t <- paste(covseries, change.times, sep="; ")
ids <- unique(subj)[!duplicated(covseries.t)] # subj ids, one with each distinct series
ncombs <- table(covseries.t)[unique(covseries.t)]# how many per series
covmat <- cbind(subject=subj, time=mf$"(time)", mf[,attr(mf,"covnames"),drop=FALSE])
covmat <- covmat[(subj %in% ids) & keep,]
list(example=covmat, # rows of the original data sufficient to define the distinct histories
hist=covseries.t) # one per subject listing their covariate history as a string
}
else NULL
}
### Estimate observed state occupancies in the data at a series of times
### Assume previous observed state is retained until next observation time
### Assumes times are sorted within patient (they are in data in msm objects)
observed.msm <- function(x, times=NULL, interp=c("start","midpoint"), censtime=Inf, subset=NULL)
{
if (!inherits(x, "msm")) stop("expected x to be a msm model")
## For general HMMs use the Viterbi estimate of the observed state.
if (!is.null(x$pci)) {
state <- x$data$mf$"(state)"[!x$data$mf$"(pci.imp)"]
time <- x$data$mf$"(time)"[!x$data$mf$"(pci.imp)"]
subject <- x$data$mf$"(subject)"
subject <- subject[!x$data$mf$"(pci.imp)"]
} else {
if ((x$hmodel$hidden && !x$emodel$misc) ||
(!x$emodel$misc && x$cmodel$ncens>0) )
state <- viterbi.msm(x)$fitted
else if (x$emodel$misc && x$cmodel$ncens>0) {
vit <- viterbi.msm(x)$fitted
state <- x$data$mf$"(state)"
state[state %in% x$cmodel$censor] <- vit[state %in% x$cmodel$censor]
## TODO for misc models with censoring, impute only censored obs states from viterbi
} else
state <- x$data$mf$"(state)"
time <- x$data$mf$"(time)"; subject <- x$data$mf$"(subject)"
}
if (is.null(subset)) subset <- unique(subject)
time <- time[subject %in% subset]
state <- state[subject %in% subset]
subject <- subject[subject %in% subset]
if (is.null(times))
times <- seq(min(time), max(time), (max(time) - min(time))/10)
states.expand <- matrix(nrow=length(unique(subject)), ncol=length(times))
pts <- unique(subject)
absorb <- absorbing.msm(x)
interp <- match.arg(interp)
if (!is.numeric(censtime)) stop("censtime should be numeric")
if (length(censtime)==1) censtime <- rep(censtime, length(pts))
else if (length(censtime)!=length(pts)) stop("censtime of length ", length(censtime), ", should be 1 or ", length(pts))
for (i in seq_along(pts)){
state.i <- state[(subject==pts[i])]
time.i <- time[(subject==pts[i])]
j <- 1
while(j <= length(times)) {
if (times[j] < time.i[1]) {
mtime <- max(which(times-time.i[1] < 0))
states.expand[i, j:mtime] <- NA
j <- mtime + 1
next;
} else if (times[j] > time.i[length(time.i)]) {
if (state.i[length(time.i)] %in% absorb && (times[j] <= censtime[i])) {
states.expand[i, j:(length(times))] <- state.i[length(time.i)]
} else states.expand[i, j:(length(times))] <- NA
break;
} else {
prevtime.ind <- max(which(time.i <= times[j]))
prevtime <- time.i[prevtime.ind]
if (interp=="midpoint") {
nexttime.ind <- min(which(time.i >= times[j]))
nexttime <- time.i[nexttime.ind]
midpoint <- (prevtime + nexttime) / 2
states.expand[i,j] <- state.i[if (times[j] <= midpoint) prevtime.ind else nexttime.ind]
} else
states.expand[i,j] <- state.i[prevtime.ind]
}
j <- j+1
}
}
obstab <- t(apply(states.expand, 2, function(y) table(factor(y, levels=seq(length.out=x$qmodel$nstates)))))
obsperc <- 100*obstab / rep(rowSums(obstab), ncol(obstab))
dimnames(obstab) <- dimnames(obsperc) <- list(times, paste("State", 1:x$qmodel$nstates))
obstab <- cbind(obstab, Total=rowSums(obstab))
covhist <- get.covhist(x, subset)
covcat <- ## distinct covariate history group each subject falls into (ordinal)
if (is.null(covhist)) rep(1, length(unique(subject)))
else match(covhist$hist, unique(covhist$hist))
risk <- matrix(nrow=length(times), ncol=length(unique(covcat)), dimnames = list(times, unique(covhist$hist)))
for (i in seq_along(unique(covcat))) {
obst <- t(apply(states.expand[covcat==unique(covcat)[i],,drop=FALSE], 2,
function(y) table(factor(y, levels=seq(length.out=x$qmodel$nstates)))))
risk[,i] <- rowSums(obst)
}
list(obstab=obstab, obsperc=obsperc, risk=risk)
}
expected.msm <- function(x,
times=NULL,
timezero=NULL,
initstates=NULL,
covariates="population",
misccovariates="mean",
piecewise.times=NULL,
piecewise.covariates=NULL,
risk=NULL,
subset=NULL,
ci=c("none","normal","bootstrap"),
cl = 0.95,
B = 1000,
cores = NULL
)
{
if (!inherits(x, "msm")) stop("expected x to be a msm model")
time <- model.extract(x$data$mf, "time")
if (is.null(times))
times <- seq(min(time), max(time), (max(time) - min(time))/10)
if (is.null(timezero)) timezero <- min(time)
if (is.null(risk)) risk <- observed.msm(x, times=times, subset=subset)$risk
exptab <- matrix(0, nrow=length(times), ncol=x$qmodel$nstates)
start <- min(which(times - timezero >= 0))
if (x$emodel$misc)
initprobs <- x$emodel$initprobs
else {
if (is.null(initstates))
initstates <- observed.msm(x, times=timezero)$obstab[1:x$qmodel$nstates]
initprobs <- initstates / sum(initstates)
}
if (length(times) >= start) {
for (j in start:length(times)) {
if (x$qcmodel$ncovs>0 && isTRUE(covariates=="population")) {
covmat <- get.covhist(x, subset=subset)$example
for (i in 1:length(unique(covmat$subject))) {
## sum expected prevalences for each covariate history observed in the data
subji <- unique(covmat$subject)[i]
ni <- sum(covmat$subject==subji)
ctimes <- covmat$time[covmat$subject==subji][-c(1,ni)]
covs <- covmat[covmat$subject==subji, attr(x$data$mf,"covnames"),drop=FALSE][-ni,,drop=FALSE]
ccovs <- list()
for (k in 1:nrow(covs)) ccovs[[k]] <- as.list(covs[k,,drop=FALSE])
pmat <- pmatrix.piecewise.msm(x, t1=timezero, t2=times[j], times=ctimes, covariates=ccovs)
expji <- risk[j,i] * initprobs %*% pmat
if (x$emodel$misc) { # return expected prev of obs (not true) states
if (x$ecmodel$ncovs==0) emat <- ematrix.msm(x, ci="none")
else {
ecovs <- if(length(ctimes)==0) ccovs else ccovs[[length(ccovs)]]
emat <- ematrix.msm(x, covariates=ecovs, ci="none")
}
expji <- expji %*% emat
}
exptab[j,] <- exptab[j,] + expji
}
}
else {
pmat <-
if (is.null(piecewise.times))
pmatrix.msm(x, t=times[j] - timezero, t1=timezero, covariates=covariates)
else
pmatrix.piecewise.msm(x, timezero, times[j], piecewise.times, piecewise.covariates)
expj <- rowSums(risk)[j] * initprobs %*% pmat
if (x$emodel$misc) # return expected prev of obs (not true) states
expj <- expj %*% ematrix.msm(x, covariates=misccovariates, ci="none")
exptab[j,] <- expj
}
}
}
exptab <- cbind(exptab, apply(exptab, 1, sum))
dimnames(exptab) <- list(times, c(rownames(x$qmodel$qmatrix),"Total"))
expperc <- 100*exptab[,1:x$qmodel$nstates] / exptab[, x$qmodel$nstates+1]
ci <- match.arg(ci)
e.ci <- switch(ci,
bootstrap = expected.ci.msm(x, times, timezero, initstates, covariates, misccovariates,
piecewise.times, piecewise.covariates, risk, cl, B, cores),
normal = expected.normci.msm(x, times, timezero, initstates, covariates, misccovariates,
piecewise.times, piecewise.covariates, risk, cl, B),
none = NULL)
res <-
if (ci=="none") list(exptab=exptab, expperc=expperc)
else list(exptab=list(estimates=exptab, ci=e.ci[[1]]),
expperc=list(estimates=expperc, ci=e.ci[[2]]))
names(res) <- c("Expected","Expected percentages")
res
}
### Table of observed and expected prevalences (works for misclassification and non-misclassification models)
#' Tables of observed and expected prevalences
#'
#' This provides a rough indication of the goodness of fit of a multi-state
#' model, by estimating the observed numbers of individuals occupying each
#' state at a series of times, and comparing these with forecasts from the
#' fitted model.
#'
#' The fitted transition probability matrix is used to forecast expected
#' prevalences from the state occupancy at the initial time. To produce the
#' expected number in state \eqn{j} at time \eqn{t} after the start, the number
#' of individuals under observation at time \eqn{t} (including those who have
#' died, but not those lost to follow-up) is multiplied by the product of the
#' proportion of individuals in each state at the initial time and the
#' transition probability matrix in the time interval \eqn{t}. The proportion
#' of individuals in each state at the "initial" time is estimated, if
#' necessary, in the same way as the observed prevalences.
#'
#' For misclassification models (fitted using an \code{ematrix}), this aims to
#' assess the fit of the full model for the \emph{observed} states. That is,
#' the combined Markov progression model for the true states and the
#' misclassification model. Thus, expected prevalences of \emph{true} states
#' are estimated from the assumed proportion occupying each state at the
#' initial time using the fitted transition probabiliy matrix. The vector of
#' expected prevalences of true states is then multiplied by the fitted
#' misclassification probability matrix to obtain the expected prevalences of
#' \emph{observed} states.
#'
#' For general hidden Markov models, the observed state is taken to be the
#' predicted underlying state from the Viterbi algorithm
#' (\code{\link{viterbi.msm}}). The goodness of fit of these states to the
#' underlying Markov model is tested.
#'
#' In any model, if there are censored states, then these are replaced by
#' imputed values of highest probability from the Viterbi algorithm in order to
#' calculate the observed state prevalences.
#'
#' For an example of this approach, see Gentleman \emph{et al.} (1994).
#'
#' @param x A fitted multi-state model produced by \code{\link{msm}}.
#' @param times Series of times at which to compute the observed and expected
#' prevalences of states.
#' @param timezero Initial time of the Markov process. Expected values are
#' forecasted from here. Defaults to the minimum of the observation times given
#' in the data.
#' @param initstates Optional vector of the same length as the number of
#' states. Gives the numbers of individuals occupying each state at the initial
#' time, to be used for forecasting expected prevalences. The default is those
#' observed in the data. These should add up to the actual number of people in
#' the study at the start.
#' @param covariates Covariate values for which to forecast expected state
#' occupancy. With the default \code{covariates="population"}, expected
#' prevalences are produced by summing model predictions over the covariates
#' observed in the original data, for a fair comparison with the observed
#' prevalences. This may be slow, particularly with continuous covariates.
#'
#' Predictions for fixed covariates can be obtained by supplying covariate
#' values in the standard way, as in \code{\link{qmatrix.msm}}. Therefore if
#' \code{covariates="population"} is too slow, using the mean observed values
#' through \code{covariates="mean"} may give a reasonable approximation.
#'
#' This argument is ignored if \code{piecewise.times} is specified. If there
#' are a mixture of time-constant and time-dependent covariates, then the
#' values for all covariates should be supplied in \code{piecewise.covariates}.
#' @param misccovariates (Misclassification models only) Values of covariates
#' on the misclassification probability matrix for converting expected true to
#' expected misclassified states. Ignored if \code{covariates="population"},
#' otherwise defaults to the mean values of the covariates in the data set.
#' @param piecewise.times Times at which piecewise-constant intensities change.
#' See \code{\link{pmatrix.piecewise.msm}} for how to specify this. Ignored if
#' \code{covariates="population"}. This is only required for
#' time-inhomogeneous models specified using explicit time-dependent
#' covariates, and should not be used for models specified using "pci".
#' @param piecewise.covariates Covariates on which the piecewise-constant
#' intensities depend. See \code{\link{pmatrix.piecewise.msm}} for how to
#' specify this. Ignored if \code{covariates="population"}.
#' @param ci If \code{"normal"}, then calculate a confidence interval for the
#' expected prevalences by simulating \code{B} random vectors from the
#' asymptotic multivariate normal distribution implied by the maximum
#' likelihood estimates (and covariance matrix) of the log transition
#' intensities and covariate effects, then calculating the expected prevalences
#' for each replicate.
#'
#' If \code{"bootstrap"} then calculate a confidence interval by non-parametric
#' bootstrap refitting. This is 1-2 orders of magnitude slower than the
#' \code{"normal"} method, but is expected to be more accurate. See
#' \code{\link{boot.msm}} for more details of bootstrapping in \pkg{msm}.
#'
#' If \code{"none"} (the default) then no confidence interval is calculated.
#' @param cl Width of the symmetric confidence interval, relative to 1
#' @param B Number of bootstrap replicates
#' @param cores Number of cores to use for bootstrapping using parallel
#' processing. See \code{\link{boot.msm}} for more details.
#' @param interp Suppose an individual was observed in states \eqn{S_{r-1}} and
#' \eqn{S_r} at two consecutive times \eqn{t_{r-1}} and \eqn{t_r}, and we want
#' to estimate 'observed' prevalences at a time \eqn{t} between \eqn{t_{r-1}}
#' and \eqn{t_r}.
#'
#' If \code{interp="start"}, then individuals are assumed to be in state
#' \eqn{S_{r-1}} at time \eqn{t}, the same state as they were at \eqn{t_{r-1}}.
#'
#' If \code{interp="midpoint"} then if \eqn{t <= (t_{r-1} + t_r) / 2}, the
#' midpoint of \eqn{t_{r-1}} and \eqn{t_r}, the state at \eqn{t} is assumed to
#' be \eqn{S_{r-1}}, otherwise \eqn{S_{r}}. This is generally more reasonable
#' for "progressive" models.
#' @param censtime Adjustment to the observed prevalences to account for
#' limited follow-up in the data.
#'
#' If the time is greater than \code{censtime} and the patient has reached an
#' absorbing state, then that subject will be removed from the risk set. For
#' example, if patients have died but would only have been observed up to this
#' time, then this avoids overestimating the proportion of people who are dead
#' at later times.
#'
#' This can be supplied as a single value, or as a vector with one element per
#' subject (after any \code{subset} has been taken), in the same order as the
#' original data. This vector also only includes subjects with complete data,
#' thus it excludes for example subjects with only one observation (thus no
#' observed transitions), and subjects for whom every observation has missing
#' values. (Note, to help construct this, the complete data used for the model
#' fit can be accessed with \code{model.frame(x)}, where \code{x} is the fitted
#' model object)
#'
#' This is ignored if it is less than the subject's maximum observation time.
#' @param subset Subset of subjects to calculate observed prevalences for.
#' @param plot Generate a plot of observed against expected prevalences. See
#' \code{\link{plot.prevalence.msm}}
#' @param ... Further arguments to pass to \code{\link{plot.prevalence.msm}}.
#' @return A list of matrices, with components:
#'
#' \item{Observed}{Table of observed numbers of individuals in each state at
#' each time}
#'
#' \item{Observed percentages}{Corresponding percentage of the individuals at
#' risk at each time.}
#'
#' \item{Expected}{Table of corresponding expected numbers.}
#'
#' \item{Expected percentages}{Corresponding percentage of the individuals at
#' risk at each time.}
#'
#' Or if \code{ci.boot = TRUE}, the component \code{Expected} is a list with
#' components \code{estimates} and \code{ci}.\cr \code{estimates} is a matrix
#' of the expected prevalences, and \code{ci} is a list of two matrices,
#' containing the confidence limits. The component \code{Expected percentages}
#' has a similar format.
#' @author C. H. Jackson \email{chris.jackson@@mrc-bsu.cam.ac.uk}
#' @seealso \code{\link{msm}}, \code{\link{summary.msm}}
#' @references Gentleman, R.C., Lawless, J.F., Lindsey, J.C. and Yan, P.
#' Multi-state Markov models for analysing incomplete disease history data with
#' illustrations for HIV disease. \emph{Statistics in Medicine} (1994) 13(3):
#' 805--821.
#'
#' Titman, A.C., Sharples, L. D. Model diagnostics for multi-state models.
#' \emph{Statistical Methods in Medical Research} (2010) 19(6):621-651.
#' @keywords models
#' @export prevalence.msm
prevalence.msm <- function(x,
times=NULL,
timezero=NULL,
initstates=NULL,
covariates="population",
misccovariates="mean",
piecewise.times=NULL,
piecewise.covariates=NULL,
ci=c("none","normal","bootstrap"),
cl = 0.95,
B = 1000,
cores = NULL,
interp=c("start","midpoint"),
censtime=Inf,
subset=NULL,
plot = FALSE, ...
)
{
if (!inherits(x, "msm")) stop("expected x to be a msm model")
## Estimate observed state occupancies in the data at a series of times
time <- model.extract(x$data$mf, "time")
if (is.null(times))
times <- seq(min(time), max(time), (max(time) - min(time))/10)
obs <- observed.msm(x, times, interp, censtime, subset)
## Work out expected state occupancies by forecasting from transition probabilities
expec <- expected.msm(x, times, timezero, initstates, covariates, misccovariates,
piecewise.times, piecewise.covariates, obs$risk, subset, ci, cl, B, cores)
res <- list(observed=obs$obstab, expected=expec[[1]], obsperc=obs$obsperc, expperc=expec[[2]])
names(res) <- c("Observed", "Expected", "Observed percentages", "Expected percentages")
class(res) <- c("msm.prevalence", class(res))
if (plot) plot.prevalence.msm(x, mintime=min(times), maxtime=max(times), timezero=timezero, initstates=initstates,
interp=interp, censtime=censtime, covariates=covariates, misccovariates=misccovariates,
piecewise.times=piecewise.times, piecewise.covariates=piecewise.covariates, ...)
res
}
#' @export
print.msm.prevalence <- function(x,...){
print(unclass(x))
}
#' Plot of observed and expected prevalences
#'
#' Provides a rough indication of goodness of fit of a multi-state model, by
#' estimating the observed numbers of individuals occupying a state at a series
#' of times, and plotting these against forecasts from the fitted model, for
#' each state. Observed prevalences are indicated as solid lines, expected
#' prevalences as dashed lines.
#'
#' See \code{\link{prevalence.msm}} for details of the assumptions underlying
#' this method.
#'
#' Observed prevalences are plotted with a solid line, and expected prevalences
#' with a dotted line.
#'
#' @param x A fitted multi-state model produced by \code{\link{msm}}.
#' @param mintime Minimum time at which to compute the observed and expected
#' prevalences of states.
#' @param maxtime Maximum time at which to compute the observed and expected
#' prevalences of states.
#' @param timezero Initial time of the Markov process. Expected values are
#' forecasted from here. Defaults to the minimum of the observation times given
#' in the data.
#' @param initstates Optional vector of the same length as the number of
#' states. Gives the numbers of individuals occupying each state at the initial
#' time, to be used for forecasting expected prevalences. The default is those
#' observed in the data. These should add up to the actual number of people in
#' the study at the start.
#' @param interp Interpolation method for observed states, see
#' \code{\link{prevalence.msm}}.
#' @param censtime Subject-specific maximum follow-up times, see
#' \code{\link{prevalence.msm}}.
#' @param subset Vector of the subject identifiers to calculated observed
#' prevalences for.
#' @param covariates Covariate values for which to forecast expected state
#' occupancy. See \code{\link{prevalence.msm}} --- if this function runs too
#' slowly, as it may if there are continuous covariates, replace
#' \code{covariates="population"} with \code{covariates="mean"}.
#' @param misccovariates (Misclassification models only) Values of covariates
#' on the misclassification probability matrix. See
#' \code{\link{prevalence.msm}}.
#' @param piecewise.times Times at which piecewise-constant intensities change.
#' See \code{\link{prevalence.msm}}.
#' @param piecewise.covariates Covariates on which the piecewise-constant
#' intensities depend. See \code{\link{prevalence.msm}}.
#' @param xlab x axis label.
#' @param ylab y axis label.
#' @param lwd.obs Line width for observed prevalences. See \code{\link{par}}.
#' @param lwd.exp Line width for expected prevalences. See \code{\link{par}}.
#' @param lty.obs Line type for observed prevalences. See \code{\link{par}}.
#' @param lty.exp Line type for expected prevalences. See \code{\link{par}}.
#' @param col.obs Line colour for observed prevalences. See \code{\link{par}}.
#' @param col.exp Line colour for expected prevalences. See \code{\link{par}}.
#' @param legend.pos Vector of the \eqn{x} and \eqn{y} position, respectively,
#' of the legend.
#' @param ... Further arguments to be passed to the generic \code{\link{plot}}
#' function.
#' @seealso \code{\link{prevalence.msm}}
#' @references Gentleman, R.C., Lawless, J.F., Lindsey, J.C. and Yan, P.
#' Multi-state Markov models for analysing incomplete disease history data with
#' illustrations for HIV disease. \emph{Statistics in Medicine} (1994) 13(3):
#' 805--821.
#' @keywords models
#' @export plot.prevalence.msm
#' @export
plot.prevalence.msm <- function(x, mintime=NULL, maxtime=NULL, timezero=NULL, initstates=NULL,
interp=c("start","midpoint"), censtime=Inf, subset=NULL,
covariates="population", misccovariates="mean",
piecewise.times=NULL, piecewise.covariates=NULL, xlab="Times",ylab="Prevalence (%)",
lwd.obs=1, lwd.exp=1, lty.obs=1, lty.exp=2,
col.obs="blue", col.exp="red", legend.pos=NULL,...){
if (!inherits(x, "msm")) stop("expected x to be a msm model")
time <- model.extract(x$data$mf, "time")
if (is.null(mintime)) mintime <- min(time)
if (is.null(maxtime)) maxtime <- max(time)
t <- seq(mintime, maxtime, length.out=100)
obs <- observed.msm(x, t, interp, censtime, subset)
expec <- expected.msm(x, t, timezero=timezero, initstates=initstates, covariates=covariates, misccovariates=misccovariates,
piecewise.times=piecewise.times, piecewise.covariates=piecewise.covariates, risk=obs$risk, subset=subset, ci="none")[[2]]
states <- seq(length.out=x$qmodel$nstates)
S <- length(states)
ncols <- ceiling(sqrt(S))
nrows <- if (floor(sqrt(S))^2 < S && S <= floor(sqrt(S))*ceiling(sqrt(S))) floor(sqrt(S)) else ceiling(sqrt(S))
par(mfrow=c(nrows, ncols))
for (i in states) {
plot(t, obs$obsperc[,i], type="l", ylim=c(0, 100), xlab=xlab, ylab=ylab, lwd=lwd.obs, lty=lty.obs, col=col.obs,
main=rownames(x$qmodel$qmatrix)[i],...)
lines(t, expec[,i], lwd=lwd.exp, lty=lty.exp, col=col.exp)
}
if (!is.numeric(legend.pos) || length(legend.pos) != 2)
legend.pos <- c(0.4*maxtime, 40)
legend(x=legend.pos[1], y=legend.pos[2], legend=c("Observed","Expected"), lty=c(lty.obs,lty.exp), lwd=c(lwd.obs,lwd.exp), col=c(col.obs,col.exp))
invisible()
}
Try the msm package in your browser
Any scripts or data that you put into this service are public.
msm documentation built on Oct. 5, 2024, 1:07 a.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.
Defines functions plot.prevalence.msm print.msm.prevalence prevalence.msm expected.msm observed.msm get.covhist
Documented in plot.prevalence.msm prevalence.msm
## Enumerate the distinct subject covariate histories in the data
## Works for time homogeneous and inhomogeneous models
## Used to calculate expected prevalences for a population with those covariates
get.covhist <- function(x, subset=NULL) {
## Keep only times where the covariate changes, or first or last obs
mf <- x$data$mf
if (x$qcmodel$ncovs > 0) {
if (!is.null(subset)) {
subs <- mf$"(subject)" %in% subset
mf <- mf[subs,,drop=FALSE]
}
subj <- match(mf$"(subject)", unique(mf$"(subject)"))
n <- length(subj)
apaste <- do.call("paste", mf[,attr(mf,"covnames"),drop=FALSE])
first <- !duplicated(subj); last <- rev(!duplicated(rev(subj)))
keep <- (c(0, apaste[1:(n-1)]) != apaste) | first | last
## Keep and tabulate unique covariate series
covseries <- split(apaste[keep], subj[keep]) # as a list of char vectors
covseries <- sapply(covseries, paste, collapse=" , ") # as one char vector, one series per pt.
## also need p matrices for different times as well as different covs.
## but only interested in cov change times if there's more than one
## transition (at least one times change point)
change.times <- mf$"(time)"; change.times[first] <- change.times[last] <- 0
change.times <- split(change.times[keep & (!(first|last))], subj[keep & (!(first|last))])
change.times <- sapply(change.times, paste, collapse= " , ")
covseries.t <- paste(covseries, change.times, sep="; ")
ids <- unique(subj)[!duplicated(covseries.t)] # subj ids, one with each distinct series
ncombs <- table(covseries.t)[unique(covseries.t)]# how many per series
covmat <- cbind(subject=subj, time=mf$"(time)", mf[,attr(mf,"covnames"),drop=FALSE])
covmat <- covmat[(subj %in% ids) & keep,]
list(example=covmat, # rows of the original data sufficient to define the distinct histories
hist=covseries.t) # one per subject listing their covariate history as a string
}
else NULL
}
### Estimate observed state occupancies in the data at a series of times
### Assume previous observed state is retained until next observation time
### Assumes times are sorted within patient (they are in data in msm objects)
observed.msm <- function(x, times=NULL, interp=c("start","midpoint"), censtime=Inf, subset=NULL)
{
if (!inherits(x, "msm")) stop("expected x to be a msm model")
## For general HMMs use the Viterbi estimate of the observed state.
if (!is.null(x$pci)) {
state <- x$data$mf$"(state)"[!x$data$mf$"(pci.imp)"]
time <- x$data$mf$"(time)"[!x$data$mf$"(pci.imp)"]
subject <- x$data$mf$"(subject)"
subject <- subject[!x$data$mf$"(pci.imp)"]
} else {
if ((x$hmodel$hidden && !x$emodel$misc) ||
(!x$emodel$misc && x$cmodel$ncens>0) )
state <- viterbi.msm(x)$fitted
else if (x$emodel$misc && x$cmodel$ncens>0) {
vit <- viterbi.msm(x)$fitted
state <- x$data$mf$"(state)"
state[state %in% x$cmodel$censor] <- vit[state %in% x$cmodel$censor]
## TODO for misc models with censoring, impute only censored obs states from viterbi
} else
state <- x$data$mf$"(state)"
time <- x$data$mf$"(time)"; subject <- x$data$mf$"(subject)"
}
if (is.null(subset)) subset <- unique(subject)
time <- time[subject %in% subset]
state <- state[subject %in% subset]
subject <- subject[subject %in% subset]
if (is.null(times))
times <- seq(min(time), max(time), (max(time) - min(time))/10)
states.expand <- matrix(nrow=length(unique(subject)), ncol=length(times))
pts <- unique(subject)
absorb <- absorbing.msm(x)
interp <- match.arg(interp)
if (!is.numeric(censtime)) stop("censtime should be numeric")
if (length(censtime)==1) censtime <- rep(censtime, length(pts))
else if (length(censtime)!=length(pts)) stop("censtime of length ", length(censtime), ", should be 1 or ", length(pts))
for (i in seq_along(pts)){
state.i <- state[(subject==pts[i])]
time.i <- time[(subject==pts[i])]
j <- 1
while(j <= length(times)) {
if (times[j] < time.i[1]) {
mtime <- max(which(times-time.i[1] < 0))
states.expand[i, j:mtime] <- NA
j <- mtime + 1
next;
} else if (times[j] > time.i[length(time.i)]) {
if (state.i[length(time.i)] %in% absorb && (times[j] <= censtime[i])) {
states.expand[i, j:(length(times))] <- state.i[length(time.i)]
} else states.expand[i, j:(length(times))] <- NA
break;
} else {
prevtime.ind <- max(which(time.i <= times[j]))
prevtime <- time.i[prevtime.ind]
if (interp=="midpoint") {
nexttime.ind <- min(which(time.i >= times[j]))
nexttime <- time.i[nexttime.ind]
midpoint <- (prevtime + nexttime) / 2
states.expand[i,j] <- state.i[if (times[j] <= midpoint) prevtime.ind else nexttime.ind]
} else
states.expand[i,j] <- state.i[prevtime.ind]
}
j <- j+1
}
}
obstab <- t(apply(states.expand, 2, function(y) table(factor(y, levels=seq(length.out=x$qmodel$nstates)))))
obsperc <- 100*obstab / rep(rowSums(obstab), ncol(obstab))
dimnames(obstab) <- dimnames(obsperc) <- list(times, paste("State", 1:x$qmodel$nstates))
obstab <- cbind(obstab, Total=rowSums(obstab))
covhist <- get.covhist(x, subset)
covcat <- ## distinct covariate history group each subject falls into (ordinal)
if (is.null(covhist)) rep(1, length(unique(subject)))
else match(covhist$hist, unique(covhist$hist))
risk <- matrix(nrow=length(times), ncol=length(unique(covcat)), dimnames = list(times, unique(covhist$hist)))
for (i in seq_along(unique(covcat))) {
obst <- t(apply(states.expand[covcat==unique(covcat)[i],,drop=FALSE], 2,
function(y) table(factor(y, levels=seq(length.out=x$qmodel$nstates)))))
risk[,i] <- rowSums(obst)
}
list(obstab=obstab, obsperc=obsperc, risk=risk)
}
expected.msm <- function(x,
times=NULL,
timezero=NULL,
initstates=NULL,
covariates="population",
misccovariates="mean",
piecewise.times=NULL,
piecewise.covariates=NULL,
risk=NULL,
subset=NULL,
ci=c("none","normal","bootstrap"),
cl = 0.95,
B = 1000,
cores = NULL
)
{
if (!inherits(x, "msm")) stop("expected x to be a msm model")
time <- model.extract(x$data$mf, "time")
if (is.null(times))
times <- seq(min(time), max(time), (max(time) - min(time))/10)
if (is.null(timezero)) timezero <- min(time)
if (is.null(risk)) risk <- observed.msm(x, times=times, subset=subset)$risk
exptab <- matrix(0, nrow=length(times), ncol=x$qmodel$nstates)
start <- min(which(times - timezero >= 0))
if (x$emodel$misc)
initprobs <- x$emodel$initprobs
else {
if (is.null(initstates))
initstates <- observed.msm(x, times=timezero)$obstab[1:x$qmodel$nstates]
initprobs <- initstates / sum(initstates)
}
if (length(times) >= start) {
for (j in start:length(times)) {
if (x$qcmodel$ncovs>0 && isTRUE(covariates=="population")) {
covmat <- get.covhist(x, subset=subset)$example
for (i in 1:length(unique(covmat$subject))) {
## sum expected prevalences for each covariate history observed in the data
subji <- unique(covmat$subject)[i]
ni <- sum(covmat$subject==subji)
ctimes <- covmat$time[covmat$subject==subji][-c(1,ni)]
covs <- covmat[covmat$subject==subji, attr(x$data$mf,"covnames"),drop=FALSE][-ni,,drop=FALSE]
ccovs <- list()
for (k in 1:nrow(covs)) ccovs[[k]] <- as.list(covs[k,,drop=FALSE])
pmat <- pmatrix.piecewise.msm(x, t1=timezero, t2=times[j], times=ctimes, covariates=ccovs)
expji <- risk[j,i] * initprobs %*% pmat
if (x$emodel$misc) { # return expected prev of obs (not true) states
if (x$ecmodel$ncovs==0) emat <- ematrix.msm(x, ci="none")
else {
ecovs <- if(length(ctimes)==0) ccovs else ccovs[[length(ccovs)]]
emat <- ematrix.msm(x, covariates=ecovs, ci="none")
}
expji <- expji %*% emat
}
exptab[j,] <- exptab[j,] + expji
}
}
else {
pmat <-
if (is.null(piecewise.times))
pmatrix.msm(x, t=times[j] - timezero, t1=timezero, covariates=covariates)
else
pmatrix.piecewise.msm(x, timezero, times[j], piecewise.times, piecewise.covariates)
expj <- rowSums(risk)[j] * initprobs %*% pmat
if (x$emodel$misc) # return expected prev of obs (not true) states
expj <- expj %*% ematrix.msm(x, covariates=misccovariates, ci="none")
exptab[j,] <- expj
}
}
}
exptab <- cbind(exptab, apply(exptab, 1, sum))
dimnames(exptab) <- list(times, c(rownames(x$qmodel$qmatrix),"Total"))
expperc <- 100*exptab[,1:x$qmodel$nstates] / exptab[, x$qmodel$nstates+1]
ci <- match.arg(ci)
e.ci <- switch(ci,
bootstrap = expected.ci.msm(x, times, timezero, initstates, covariates, misccovariates,
piecewise.times, piecewise.covariates, risk, cl, B, cores),
normal = expected.normci.msm(x, times, timezero, initstates, covariates, misccovariates,
piecewise.times, piecewise.covariates, risk, cl, B),
none = NULL)
res <-
if (ci=="none") list(exptab=exptab, expperc=expperc)
else list(exptab=list(estimates=exptab, ci=e.ci[[1]]),
expperc=list(estimates=expperc, ci=e.ci[[2]]))
names(res) <- c("Expected","Expected percentages")
res
}
### Table of observed and expected prevalences (works for misclassification and non-misclassification models)
#' Tables of observed and expected prevalences
#'
#' This provides a rough indication of the goodness of fit of a multi-state
#' model, by estimating the observed numbers of individuals occupying each
#' state at a series of times, and comparing these with forecasts from the
#' fitted model.
#'
#' The fitted transition probability matrix is used to forecast expected
#' prevalences from the state occupancy at the initial time. To produce the
#' expected number in state \eqn{j} at time \eqn{t} after the start, the number
#' of individuals under observation at time \eqn{t} (including those who have
#' died, but not those lost to follow-up) is multiplied by the product of the
#' proportion of individuals in each state at the initial time and the
#' transition probability matrix in the time interval \eqn{t}. The proportion
#' of individuals in each state at the "initial" time is estimated, if
#' necessary, in the same way as the observed prevalences.
#'
#' For misclassification models (fitted using an \code{ematrix}), this aims to
#' assess the fit of the full model for the \emph{observed} states. That is,
#' the combined Markov progression model for the true states and the
#' misclassification model. Thus, expected prevalences of \emph{true} states
#' are estimated from the assumed proportion occupying each state at the
#' initial time using the fitted transition probabiliy matrix. The vector of
#' expected prevalences of true states is then multiplied by the fitted
#' misclassification probability matrix to obtain the expected prevalences of
#' \emph{observed} states.
#'
#' For general hidden Markov models, the observed state is taken to be the
#' predicted underlying state from the Viterbi algorithm
#' (\code{\link{viterbi.msm}}). The goodness of fit of these states to the
#' underlying Markov model is tested.
#'
#' In any model, if there are censored states, then these are replaced by
#' imputed values of highest probability from the Viterbi algorithm in order to
#' calculate the observed state prevalences.
#'
#' For an example of this approach, see Gentleman \emph{et al.} (1994).
#'
#' @param x A fitted multi-state model produced by \code{\link{msm}}.
#' @param times Series of times at which to compute the observed and expected
#' prevalences of states.
#' @param timezero Initial time of the Markov process. Expected values are
#' forecasted from here. Defaults to the minimum of the observation times given
#' in the data.
#' @param initstates Optional vector of the same length as the number of
#' states. Gives the numbers of individuals occupying each state at the initial
#' time, to be used for forecasting expected prevalences. The default is those
#' observed in the data. These should add up to the actual number of people in
#' the study at the start.
#' @param covariates Covariate values for which to forecast expected state
#' occupancy. With the default \code{covariates="population"}, expected
#' prevalences are produced by summing model predictions over the covariates
#' observed in the original data, for a fair comparison with the observed
#' prevalences. This may be slow, particularly with continuous covariates.
#'
#' Predictions for fixed covariates can be obtained by supplying covariate
#' values in the standard way, as in \code{\link{qmatrix.msm}}. Therefore if
#' \code{covariates="population"} is too slow, using the mean observed values
#' through \code{covariates="mean"} may give a reasonable approximation.
#'
#' This argument is ignored if \code{piecewise.times} is specified. If there
#' are a mixture of time-constant and time-dependent covariates, then the
#' values for all covariates should be supplied in \code{piecewise.covariates}.
#' @param misccovariates (Misclassification models only) Values of covariates
#' on the misclassification probability matrix for converting expected true to
#' expected misclassified states. Ignored if \code{covariates="population"},
#' otherwise defaults to the mean values of the covariates in the data set.
#' @param piecewise.times Times at which piecewise-constant intensities change.
#' See \code{\link{pmatrix.piecewise.msm}} for how to specify this. Ignored if
#' \code{covariates="population"}. This is only required for
#' time-inhomogeneous models specified using explicit time-dependent
#' covariates, and should not be used for models specified using "pci".
#' @param piecewise.covariates Covariates on which the piecewise-constant
#' intensities depend. See \code{\link{pmatrix.piecewise.msm}} for how to
#' specify this. Ignored if \code{covariates="population"}.
#' @param ci If \code{"normal"}, then calculate a confidence interval for the
#' expected prevalences by simulating \code{B} random vectors from the
#' asymptotic multivariate normal distribution implied by the maximum
#' likelihood estimates (and covariance matrix) of the log transition
#' intensities and covariate effects, then calculating the expected prevalences
#' for each replicate.
#'
#' If \code{"bootstrap"} then calculate a confidence interval by non-parametric
#' bootstrap refitting. This is 1-2 orders of magnitude slower than the
#' \code{"normal"} method, but is expected to be more accurate. See
#' \code{\link{boot.msm}} for more details of bootstrapping in \pkg{msm}.
#'
#' If \code{"none"} (the default) then no confidence interval is calculated.
#' @param cl Width of the symmetric confidence interval, relative to 1
#' @param B Number of bootstrap replicates
#' @param cores Number of cores to use for bootstrapping using parallel
#' processing. See \code{\link{boot.msm}} for more details.
#' @param interp Suppose an individual was observed in states \eqn{S_{r-1}} and
#' \eqn{S_r} at two consecutive times \eqn{t_{r-1}} and \eqn{t_r}, and we want
#' to estimate 'observed' prevalences at a time \eqn{t} between \eqn{t_{r-1}}
#' and \eqn{t_r}.
#'
#' If \code{interp="start"}, then individuals are assumed to be in state
#' \eqn{S_{r-1}} at time \eqn{t}, the same state as they were at \eqn{t_{r-1}}.
#'
#' If \code{interp="midpoint"} then if \eqn{t <= (t_{r-1} + t_r) / 2}, the
#' midpoint of \eqn{t_{r-1}} and \eqn{t_r}, the state at \eqn{t} is assumed to
#' be \eqn{S_{r-1}}, otherwise \eqn{S_{r}}. This is generally more reasonable
#' for "progressive" models.
#' @param censtime Adjustment to the observed prevalences to account for
#' limited follow-up in the data.
#'
#' If the time is greater than \code{censtime} and the patient has reached an
#' absorbing state, then that subject will be removed from the risk set. For
#' example, if patients have died but would only have been observed up to this
#' time, then this avoids overestimating the proportion of people who are dead
#' at later times.
#'
#' This can be supplied as a single value, or as a vector with one element per
#' subject (after any \code{subset} has been taken), in the same order as the
#' original data. This vector also only includes subjects with complete data,
#' thus it excludes for example subjects with only one observation (thus no
#' observed transitions), and subjects for whom every observation has missing
#' values. (Note, to help construct this, the complete data used for the model
#' fit can be accessed with \code{model.frame(x)}, where \code{x} is the fitted
#' model object)
#'
#' This is ignored if it is less than the subject's maximum observation time.
#' @param subset Subset of subjects to calculate observed prevalences for.
#' @param plot Generate a plot of observed against expected prevalences. See
#' \code{\link{plot.prevalence.msm}}
#' @param ... Further arguments to pass to \code{\link{plot.prevalence.msm}}.
#' @return A list of matrices, with components:
#'
#' \item{Observed}{Table of observed numbers of individuals in each state at
#' each time}
#'
#' \item{Observed percentages}{Corresponding percentage of the individuals at
#' risk at each time.}
#'
#' \item{Expected}{Table of corresponding expected numbers.}
#'
#' \item{Expected percentages}{Corresponding percentage of the individuals at
#' risk at each time.}
#'
#' Or if \code{ci.boot = TRUE}, the component \code{Expected} is a list with
#' components \code{estimates} and \code{ci}.\cr \code{estimates} is a matrix
#' of the expected prevalences, and \code{ci} is a list of two matrices,
#' containing the confidence limits. The component \code{Expected percentages}
#' has a similar format.
#' @author C. H. Jackson \email{chris.jackson@@mrc-bsu.cam.ac.uk}
#' @seealso \code{\link{msm}}, \code{\link{summary.msm}}
#' @references Gentleman, R.C., Lawless, J.F., Lindsey, J.C. and Yan, P.
#' Multi-state Markov models for analysing incomplete disease history data with
#' illustrations for HIV disease. \emph{Statistics in Medicine} (1994) 13(3):
#' 805--821.
#'
#' Titman, A.C., Sharples, L. D. Model diagnostics for multi-state models.
#' \emph{Statistical Methods in Medical Research} (2010) 19(6):621-651.
#' @keywords models
#' @export prevalence.msm
prevalence.msm <- function(x,
times=NULL,
timezero=NULL,
initstates=NULL,
covariates="population",
misccovariates="mean",
piecewise.times=NULL,
piecewise.covariates=NULL,
ci=c("none","normal","bootstrap"),
cl = 0.95,
B = 1000,
cores = NULL,
interp=c("start","midpoint"),
censtime=Inf,
subset=NULL,
plot = FALSE, ...
)
{
if (!inherits(x, "msm")) stop("expected x to be a msm model")
## Estimate observed state occupancies in the data at a series of times
time <- model.extract(x$data$mf, "time")
if (is.null(times))
times <- seq(min(time), max(time), (max(time) - min(time))/10)
obs <- observed.msm(x, times, interp, censtime, subset)
## Work out expected state occupancies by forecasting from transition probabilities
expec <- expected.msm(x, times, timezero, initstates, covariates, misccovariates,
piecewise.times, piecewise.covariates, obs$risk, subset, ci, cl, B, cores)
res <- list(observed=obs$obstab, expected=expec[[1]], obsperc=obs$obsperc, expperc=expec[[2]])
names(res) <- c("Observed", "Expected", "Observed percentages", "Expected percentages")
class(res) <- c("msm.prevalence", class(res))
if (plot) plot.prevalence.msm(x, mintime=min(times), maxtime=max(times), timezero=timezero, initstates=initstates,
interp=interp, censtime=censtime, covariates=covariates, misccovariates=misccovariates,
piecewise.times=piecewise.times, piecewise.covariates=piecewise.covariates, ...)
res
}
#' @export
print.msm.prevalence <- function(x,...){
print(unclass(x))
}
#' Plot of observed and expected prevalences
#'
#' Provides a rough indication of goodness of fit of a multi-state model, by
#' estimating the observed numbers of individuals occupying a state at a series
#' of times, and plotting these against forecasts from the fitted model, for
#' each state. Observed prevalences are indicated as solid lines, expected
#' prevalences as dashed lines.
#'
#' See \code{\link{prevalence.msm}} for details of the assumptions underlying
#' this method.
#'
#' Observed prevalences are plotted with a solid line, and expected prevalences
#' with a dotted line.
#'
#' @param x A fitted multi-state model produced by \code{\link{msm}}.
#' @param mintime Minimum time at which to compute the observed and expected
#' prevalences of states.
#' @param maxtime Maximum time at which to compute the observed and expected
#' prevalences of states.
#' @param timezero Initial time of the Markov process. Expected values are
#' forecasted from here. Defaults to the minimum of the observation times given
#' in the data.
#' @param initstates Optional vector of the same length as the number of
#' states. Gives the numbers of individuals occupying each state at the initial
#' time, to be used for forecasting expected prevalences. The default is those
#' observed in the data. These should add up to the actual number of people in
#' the study at the start.
#' @param interp Interpolation method for observed states, see
#' \code{\link{prevalence.msm}}.
#' @param censtime Subject-specific maximum follow-up times, see
#' \code{\link{prevalence.msm}}.
#' @param subset Vector of the subject identifiers to calculated observed
#' prevalences for.
#' @param covariates Covariate values for which to forecast expected state
#' occupancy. See \code{\link{prevalence.msm}} --- if this function runs too
#' slowly, as it may if there are continuous covariates, replace
#' \code{covariates="population"} with \code{covariates="mean"}.
#' @param misccovariates (Misclassification models only) Values of covariates
#' on the misclassification probability matrix. See
#' \code{\link{prevalence.msm}}.
#' @param piecewise.times Times at which piecewise-constant intensities change.
#' See \code{\link{prevalence.msm}}.
#' @param piecewise.covariates Covariates on which the piecewise-constant
#' intensities depend. See \code{\link{prevalence.msm}}.
#' @param xlab x axis label.
#' @param ylab y axis label.
#' @param lwd.obs Line width for observed prevalences. See \code{\link{par}}.
#' @param lwd.exp Line width for expected prevalences. See \code{\link{par}}.
#' @param lty.obs Line type for observed prevalences. See \code{\link{par}}.
#' @param lty.exp Line type for expected prevalences. See \code{\link{par}}.
#' @param col.obs Line colour for observed prevalences. See \code{\link{par}}.
#' @param col.exp Line colour for expected prevalences. See \code{\link{par}}.
#' @param legend.pos Vector of the \eqn{x} and \eqn{y} position, respectively,
#' of the legend.
#' @param ... Further arguments to be passed to the generic \code{\link{plot}}
#' function.
#' @seealso \code{\link{prevalence.msm}}
#' @references Gentleman, R.C., Lawless, J.F., Lindsey, J.C. and Yan, P.
#' Multi-state Markov models for analysing incomplete disease history data with
#' illustrations for HIV disease. \emph{Statistics in Medicine} (1994) 13(3):
#' 805--821.
#' @keywords models
#' @export plot.prevalence.msm
#' @export
plot.prevalence.msm <- function(x, mintime=NULL, maxtime=NULL, timezero=NULL, initstates=NULL,
interp=c("start","midpoint"), censtime=Inf, subset=NULL,
covariates="population", misccovariates="mean",
piecewise.times=NULL, piecewise.covariates=NULL, xlab="Times",ylab="Prevalence (%)",
lwd.obs=1, lwd.exp=1, lty.obs=1, lty.exp=2,
col.obs="blue", col.exp="red", legend.pos=NULL,...){
if (!inherits(x, "msm")) stop("expected x to be a msm model")
time <- model.extract(x$data$mf, "time")
if (is.null(mintime)) mintime <- min(time)
if (is.null(maxtime)) maxtime <- max(time)
t <- seq(mintime, maxtime, length.out=100)
obs <- observed.msm(x, t, interp, censtime, subset)
expec <- expected.msm(x, t, timezero=timezero, initstates=initstates, covariates=covariates, misccovariates=misccovariates,
piecewise.times=piecewise.times, piecewise.covariates=piecewise.covariates, risk=obs$risk, subset=subset, ci="none")[[2]]
states <- seq(length.out=x$qmodel$nstates)
S <- length(states)
ncols <- ceiling(sqrt(S))
nrows <- if (floor(sqrt(S))^2 < S && S <= floor(sqrt(S))*ceiling(sqrt(S))) floor(sqrt(S)) else ceiling(sqrt(S))
par(mfrow=c(nrows, ncols))
for (i in states) {
plot(t, obs$obsperc[,i], type="l", ylim=c(0, 100), xlab=xlab, ylab=ylab, lwd=lwd.obs, lty=lty.obs, col=col.obs,
main=rownames(x$qmodel$qmatrix)[i],...)
lines(t, expec[,i], lwd=lwd.exp, lty=lty.exp, col=col.exp)
}
if (!is.numeric(legend.pos) || length(legend.pos) != 2)
legend.pos <- c(0.4*maxtime, 40)
legend(x=legend.pos[1], y=legend.pos[2], legend=c("Observed","Expected"), lty=c(lty.obs,lty.exp), lwd=c(lwd.obs,lwd.exp), col=c(col.obs,col.exp))
invisible()
}
Try the msm package in your browser
Any scripts or data that you put into this service are public.
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.