For a general Markov chain multi-state model with interval censored
transitions calculate the NPMLE of the transition intensities. The estimates
are returned as an msfit object.
npmsm(
gd,
tmat,
method = c("multinomial", "poisson"),
inits = c("equalprob", "homogeneous", "unif", "beta"),
beta_params,
support_manual,
exact,
maxit = 100,
tol = 1e-04,
conv_crit = c("haz", "prob", "lik"),
verbose = FALSE,
manual = FALSE,
newmet = FALSE,
include_inf = FALSE,
checkMLE = TRUE,
checkMLE_tol = 1e-04,
prob_tol = tol/10,
remove_redundant = TRUE,
remove_bins = FALSE,
estimateSupport = FALSE,
init_int = c(0, 0),
...
)A list with the following entries:
A: A list of class msfit containing
the cumulative intensities for each transition and the transition matrix used;
Ainit: Initial intensities, in an object of class msfit;
gd: Data used for the estimation procedure;
ll: Log-likelihood value of the procedure at the last iteration;
delta: Change in log-likelihood value at the last iteration;
it: Number of iterations of the procedure;
taus: Unique time points of the data, the cumulative intensity only makes jumps at these time points.;
tmat: The transition matrix used, see transMat;
tmat2: A summary of the transitions in the model, see to.trans2;
ll_history: The log-likelihood value at every iteration of the procedure;
KKT_violated: How often were KKT conditions violated during maximisation of the likelihood? In other words, how often did we hit the optimization boundary during the procedure?;
min_time: The smallest time of an observation in the used data. Note that the smallest time in the data is used as zero reference;
reduced_gradient: The reduced gradient at the last iteration. Rows indicate the transitions and columns the unique observation times;
isMLE: Has the procedure converged to the NPMLE? Checked
using checkMLE_tol;
langrangemultiplier: The lagrange multipliers at the last iteration;
aghmat: A matrix representation of the transition intensities in A.
Rows represent transitions and columns unique observation times;
Ygk: The summed at-risk indicator for all subjects in the data at the last iteration. Rows represent transitions and columns unique observation times;
Dmk: The summed probability of making a transition for all subjects at the last iteration. Rows represent transitions and columns unique observation times;
method: Method used for the optimization procedure;
maxit: Maximum number of allowed iterations;
tol: Tolerance of the convergence procedure;
conv_crit: Convergence criterion of the procedure;
checkMLE: Was the reduced gradient checked at the last iteration to determine convergence?;
checkMLE_tol: The tolerance of the checkMLE procedure;
prob_tol: Tolerance for probabilities to be set to zero;
remove_redundant: Were redundant observations removed before performing the procedure?;
A data.frame with the following named columns
id:Subject idenitifier;
state:State at which the subject is observed at time;
time:Time at which the subject is observed;
The true transition time between states is then interval censored between the times.
A transition matrix as created by transMat
Which method should be used for the EM algorithm. Choices are
c("multinomial", "poisson"), with multinomial the default. Multinomial
will use the EM algorithm described in Gomon and Putter (2024) and Poisson
will use the EM algorithm described in Gu et al. (2023).
Which distribution should be used to generate the initial estimates
of the intensities in the EM algorithm. One of c("equalprob", "unif", "beta"),
with "equalprob" assigning 1/K to each intensity, with K the number of distinct
observation times (length(unique(gd[, "time"]))). For "unif", each
intensity is sampled from the Unif[0,1]
distribution and for "beta" each intensity is sampled from the Beta(a, b) distribution.
If "beta" is chosen, the argument beta_params must be specified as a
vector of length 2 containing the parameters of the beta distribution.
Default = "equalprob".
A vector of length 2 specifying the beta distribution parameters
for initial distribution generation. First entry will be used as shape1
and second entry as shape2. See help(rbeta). Only used if inits = "beta".
Used for specifying a manual support region for the transitions.
A list of length the number of transitions in tmat,
each list element containing a data frame with 2 named columns L and R indicating the
left and right values of the support intervals. When specified, all intensities
outside of these intervals will be set to zero for the corresponding transitions.
Intensities set to zero cannot be changed by the EM algorithm. Will use inits = "equalprob".
Numeric vector indicating to which states transitions are observed at exact times.
Must coincide with the column number in tmat.
Maximum number of iterations. Default = 100.
Tolerance of the convergence procedure. A change in the value of
conv_crit in an iteration of less than tol will make the procedure stop.
Convergence criterion. Stops procedure when the difference
in the chosen quantity between two consecutive iterations is smaller
than the tolerance level tol. One of the following:
Stop when change in maximum estimated intensities (hazards) < tol.
Stop when change in estimated probabilities < tol.
Stop when change in observed-data likelihood < tol.
Default is "haz". The options "haz" and "lik" can be compared across different
methods, but "prob" is dependent on the chosen method. Most
conservative (requiring most iterations) is "prob", followed by "haz" and finally "lik".
Should iteration messages be printed? Default is FALSE
Manually specify starting transition intensities? If TRUE,
the transition intensity for each bin for each transition must be entered manually.
DO NOT USE for large data sets, and in general it is not adviced to use this.
Should contributions after last observation time also be used in the likelihood? Default is FALSE.
Should an additional bin from the largest observed time to infinity be included in the algorithm? Default is FALSE.
Should a check be performed whether the estimate has converged
towards a local maximum? This is done by comparing
the reduced gradient to the value of checkMLE_tol. Default is TRUE.
Tolerance for checking whether the estimate has converged to local maximum.
Whenever an estimated transition intensity is smaller than prob_tol, it is assumed
to be zero and its reduced gradient is not considered for determining whether
the maximum has been reached. Default = 1e-4.
If an estimated probability is smaller than prob_tol,
it will be set to zero during estimation. Default value is tol/10.
Should redundant observations be removed before running the algorithm? An observation is redundant when the same state has been observed more than 3 times consecutively, or if it is a repeat observation of an absorbing state. Default is TRUE.
Should a bin be removed during the algorithm if all
estimated intensities are zero for a single bin? Can improve
computation speed for large data sets. Note that zero means the estimated intensities
are smaller than prob_tol. Default is FALSE.
Should the support of the transitions be estimated using
the result of Hudgens (2005)? Currently produces incorrect support sets -
DO NOT USE. Default = FALSE
A vector of length 2, with the first entry indicating what percentage of mass should be distributed over (second entry) what percentage of all first bins. Default is c(0, 0), in which case the argument is ignored. This argument has no practical uses and only exists for demonstration purposes in the related article.
Further arguments to estimate_support_msm
Denote the unique observation times in the data as \(0 = \tau_0, \tau_1, \ldots, \tau_K\) Let \(g, h \in H\) denote the possible states in the model and \(X(t)\) the state of the process at time t.
Then this function can be used to estimate the transition intensities \(\alpha_{gh}^k = \alpha_{gh}(\tau_k)\).
Having obtained these estimated, it is possible to recover the transition probabilities
\(\mathbf{P}(X(t) = h | X(s) = g)\) for \(t > s\) using
the transprob functions.
D. Gomon and H. Putter, Non-parametric estimation of transition intensities in interval censored Markov multi-state models without loops, arXiv, 2024, tools:::Rd_expr_doi("10.48550/arXiv.2409.07176")
Y. Gu, D. Zeng, G. Heiss, and D. Y. Lin, Maximum likelihood estimation for semiparametric regression models with interval-censored multistate data, Biometrika, Nov. 2023, tools:::Rd_expr_doi("10.1093/biomet/asad073")
Michael G. Hudgens, On Nonparametric Maximum Likelihood Estimation with Interval Censoring and Left Truncation, Journal of the Royal Statistical Society Series B: Statistical Methodology, Volume 67, Issue 4, September 2005, Pages 573-587, tools:::Rd_expr_doi("10.1111/j.1467-9868.2005.00516.x")
transprob for calculating transition probabilities,
plot.npmsm for plotting the cumulative intensities,
print.npmsm for printing some output summaries,
visualise_msm for visualising data,
msfit for details on the output object.
#Create transition matrix using mstate functionality: illness-death
if(require(mstate)){
tmat <- mstate::trans.illdeath()
}
#Write a function for evaluation times: observe at 0 and uniform inter-observation times.
eval_times <- function(n_obs, stop_time){
cumsum( c( 0, runif( n_obs-1, 0, 2*(stop_time-4)/(n_obs-1) ) ) )
}
#Use built_in function to simulate illness-death data
#from Weibull distributions for each transition
sim_dat <- sim_id_weib(n = 50, n_obs = 6, stop_time = 15, eval_times = eval_times,
start_state = "stable", shape = c(0.5, 0.5, 2), scale = c(5, 10, 10/gamma(1.5)))
tmat <- mstate::trans.illdeath()
#Fit the model using method = "multinomial"
mod_fit <- npmsm(gd = sim_dat, tmat = tmat, tol = 1e-2)
#Plot the cumulative intensities for each transition
plot(mod_fit)
Run the code above in your browser using DataLab