Predicting Transition Probabilities with probtrans

From Fitted Hazards to State Probabilities

A fitted transition-specific hazard describes instantaneous movement out of one state and into another. It is not itself the probability of occupying a state at a later time: that probability depends on the possible paths through the state graph, including remaining in a state and moving through intermediate states. The mstate function probtrans carries out this propagation from fitted transition hazards.

This tutorial assumes the transition structure has been encoded and transition-specific Cox models have been fitted. The introductory material on Aalen–Johansen state occupancy establishes the probability-matrix idea; here the focus is how to calculate and interpret model-based transition probabilities with probtrans.

Definition: Let \(A(t)\) be the cumulative transition-hazard matrix. For distinct states \(r\) and \(s\), \(dA_{rs}(t)\) is the estimated hazard increment for transition \(r\to s\) at time \(t\). The diagonal increment is \(dA_{rr}(t)=-\sum_{s\ne r}dA_{rs}(t)\), summing over transitions out of state \(r\). The transition-probability matrix \(P(s,t)\) has entry \(P_{rs}(s,t)=P(X(t)=s\mid X(s)=r)\), for a process starting in state \(r\) at time \(s\).

In the product-integral calculation, each hazard update contributes a matrix \(I+dA(u)\), where \(I\) is the identity matrix. If probabilities are represented as a row vector, then the update is \(p(u)=p(u-) [I+dA(u)]\). The order of updates follows time: the state distribution after one update is the input to the next.

P(s,t)=\prod_{s<u\leq t}[I+dA(u)]

The result depends on the fitted transition hazards and on the model assumptions behind them. In particular, the transition intensities must describe the process on the chosen time scale and, for the usual Markov interpretation, the current state and time (and any covariates held fixed for the prediction) must suffice to describe the future transition rates. Cox models also require their hazard-model assumptions to be reasonable, including proportional hazards for any fitted proportional-hazards terms. Censoring must be adequately handled by the estimation procedure, typically by being independent of the event process conditional on modeled information.

Calling probtrans and Reading Its Output

The usual workflow is to obtain estimated cumulative transition hazards with mstate::msfit and pass that fitted-hazard object to mstate::probtrans. The transition map used for the fitted hazards must match the transition structure used in the analysis. The newdata supplied to msfit must provide covariate information in the form required by the fitted model; for a transition-stratified fit, it must also identify the transitions appropriately. The exact model-fitting setup may differ, but the prediction step follows the same pattern.

For example, once coxfit, trans, and an appropriate newdata object nd have been prepared, the calls are:

msf <- mstate::msfit(coxfit, trans, newdata = nd)
pt <- mstate::probtrans(msf, predt = 0, variance = TRUE)

The argument predt sets the prediction start time. With predt = 0, the calculation starts at time zero. The returned object is a list: pt[[1]] contains probabilities for a process starting in state 1, pt[[2]] for a process starting in state 2, and so on, according to the state numbering in the transition map. Each table gives times and state probabilities, commonly named pstate1, pstate2, and so forth. When variance estimation is requested, standard-error columns are also returned. Check the column names and times in the actual object rather than assuming a particular row order.

Worked Example: Two-State Probabilities from Hazard Increments

Suppose a hypothetical fitted model has two states: state 1 is alive without the event, and state 2 is an absorbing event state. Assume the estimated transition-hazard increments from state 1 to state 2 are 0.20 at time 1 and 0.30 at time 2. There are no other transitions.

At time 1, the update matrix is \(I+dA(1)\). The outgoing increment is 0.20, so the diagonal increment for state 1 is \(-0.20\). State 2 has no outgoing transition.

I+dA(1)= \begin{pmatrix} 0.80 & 0.20\\ 0 & 1 \end{pmatrix}

Starting in state 1 gives the row vector \((1,0)\). After the first update, \((1,0)(I+dA(1))=(0.80,0.20)\). Thus the model-based probabilities at time 1 are 0.80 in state 1 and 0.20 in state 2.

At time 2, the second update matrix is \(\begin{pmatrix}0.70&0.30\\0&1\end{pmatrix}\). Multiplying the current probabilities by it gives \((0.80,0.20)\begin{pmatrix}0.70&0.30\\0&1\end{pmatrix}=(0.56,0.44)\). The probability of occupying state 2 by time 2 is therefore 0.44, not the sum 0.20 + 0.30: some individuals have already entered state 2 by time 1 and cannot transition into it a second time.

In an applied analysis, probtrans performs this propagation from the cumulative hazards in the msfit object. If the hypothetical increments above were the increments supplied to that calculation, the table for a starting state of 1 would report pstate1 = 0.80 and pstate2 = 0.20 at time 1, then pstate1 = 0.56 and pstate2 = 0.44 at time 2.

Why the Matrix Updates Need Checking

Each row of \(I+dA(u)\) sums to one by construction: the diagonal subtracts the total outgoing increment, while the off-diagonal entries add those increments back. That property alone does not guarantee that the update is a valid probability update. For example, if a state has one outgoing increment of 2, its update row is \((-1,2)\). The row sums to one, but the negative entry can produce a negative state probability.

For an update to preserve nonnegative probabilities, every outgoing increment sum for a state must be at most 1. Then the diagonal entry, one minus that sum, is nonnegative; if the off-diagonal increments are also nonnegative, all entries of the update matrix are nonnegative and its row sums are one. Under those conditions, multiplying a probability row vector by the update gives another probability row vector. If an outgoing sum exceeds 1, do not interpret the resulting product-integral values as valid probabilities merely because the rows sum to one.

Worked Example: A Three-State Illness–Death Model

Consider a hypothetical model with state 1 = well, state 2 = ill, and state 3 = dead. Allowed transitions are 1 to 2, 1 to 3, and 2 to 3. Suppose the fitted hazards yield these increments at time 1: 0.20 for 1 to 2, 0.10 for 1 to 3, and 0.25 for 2 to 3. The outgoing sums are 0.30 from state 1 and 0.25 from state 2, each no greater than 1, so the update has nonnegative entries.

I+dA(1)= \begin{pmatrix} 0.70 & 0.20 & 0.10\\ 0 & 0.75 & 0.25\\ 0 & 0 & 1 \end{pmatrix}

For a process starting in state 1, the initial row vector is \((1,0,0)\). After time 1 it is \((1,0,0)(I+dA(1))=(0.70,0.20,0.10)\). The entries sum to one: 0.70 remains well, 0.20 is ill, and 0.10 is dead.

Now suppose the increments at time 2 are 0.30 for 1 to 2, 0.10 for 1 to 3, and 0.20 for 2 to 3. The outgoing sums are 0.40 and 0.20, again satisfying the nonnegativity condition. The update matrix is \(\begin{pmatrix}0.60&0.30&0.10\\0&0.80&0.20\\0&0&1\end{pmatrix}\). The probabilities become \((0.70,0.20,0.10)\begin{pmatrix}0.60&0.30&0.10\\0&0.80&0.20\\0&0&1\end{pmatrix}=(0.42,0.37,0.21)\). At time 2, the predicted probabilities of being well, ill, and dead are 0.42, 0.37, and 0.21, respectively. They sum to one and are all nonnegative.

In R, these are the kinds of probabilities to inspect in pt[[1]], the table for a process starting in state 1. The arithmetic also shows why the probability of being ill at time 2 is not just the new 1-to-2 increment: it includes both paths that remain in state 2 and those that move from well to ill during the second update.

Extracting Probabilities for a Chosen Starting State

The list index identifies the starting state; the pstate columns identify the state occupied at each output time. For instance, pt[[2]] is the table for a process starting in state 2, while pstate3 within that table gives the probability of occupying state 3. These are conditional transition probabilities from the specified starting state, not a population-average distribution unless the starting-state distribution has been combined appropriately.

Worked Example: Probability of Death Starting from the Ill State

Use the three-state model above, but start in state 2 (ill) at time 0. There is no transition from state 2 back to state 1. The state-2-to-state-3 increments are 0.25 at time 1 and 0.20 at time 2. The probability of remaining ill after the first update is \(1-0.25=0.75\); after the second it is \(0.75(1-0.20)=0.60\). Consequently, the state-3 probability at time 2 is \(1-0.60=0.40\).

The corresponding extraction is pt[[2]], followed by reading pstate2 and pstate3 at time 2. In this hypothetical example, those values are 0.60 and 0.40. The interpretation is: among processes represented as starting in the ill state at time 0, the fitted model predicts probability 0.40 of occupying the dead state at time 2. It is not a claim that an individual has a 40% chance of dying regardless of their starting state, covariates, or time scale.

Confirm the state map.
Check that the state numbers and allowed transitions in the fitted hazards are the intended ones. A mismatch changes the meaning of every probability column.

Compute fitted hazards and probabilities.
Pass the fitted transition hazards to probtrans and choose predt to represent the prediction start time.

Select the starting state.
Use the list element corresponding to that state, then locate the requested time and pstate column.

Check the result.
At each reported time, inspect whether the probabilities are between 0 and 1 and whether they sum to approximately 1. Also inspect the hazard increments for outgoing sums above 1.

Diagnostics, Alternatives, and Limitations

Probability tables should be checked, not merely printed and accepted. At each time, the probabilities for all states should sum to approximately one, allowing for rounding. Values outside the interval from 0 to 1, or totals that depart materially from one, are warnings. Inspect the increments for each state and time: negative hazard increments, incompatible transition coding, or outgoing sums greater than 1 can invalidate the update. Also check that the prediction times are within the range where fitted hazards are supported and that the requested starting state has a meaningful interpretation.

The model-based calculations depend on the transition-specific hazard estimates. Poor Cox-model fit, sparse risk sets, violations of proportional hazards, or an inadequate Markov assumption can make the resulting probabilities unreliable even when the matrix arithmetic is valid. Standard errors quantify estimation uncertainty under the fitted procedure; they do not correct misspecification or poor support. When comparing predictions with observed state occupancy, keep in mind that observed occupancy and model-based predictions are different quantities and should be assessed in a way that respects their respective assumptions.

The main alternative is the empirical Aalen–Johansen estimator, which estimates state occupancy directly from observed transition histories without requiring transition-specific Cox regression. It is a natural choice for describing occupancy in the observed population, while probtrans combines fitted transition hazards to produce model-based probabilities. The model-based approach can use regression structure and specified covariate values, but depends on the hazard models and their assumptions. Neither approach makes unsupported extrapolation safe; predictions far beyond observed follow-up can be especially unstable.

Common mistake: “Each update matrix has rows summing to one, so it must always return valid probabilities.” This is false: a row such as \((-1,2)\) sums to one but has a negative entry. The correction is to check that every outgoing increment sum is at most 1, as well as confirming nonnegative off-diagonal increments, before relying on the update to preserve probabilities.

Check Your Understanding

Answer each question, then compare your reasoning with the short solutions.

In a three-state model, what does pt[[3]] represent, and what does pstate2 in that table represent?
A state has outgoing increments 0.35 and 0.25 at one update. What is its diagonal entry in \(I+dA\), and does this row meet the nonnegativity condition?
Why is a 1-to-2 increment at a later time not, by itself, the probability of occupying state 2 at that time?
What is one key distinction between probabilities from probtrans and empirical Aalen–Johansen occupancy estimates?

Solutions. (1) pt[[3]] gives probabilities for a process starting in state 3; pstate2 is its probability of occupying state 2 at each listed time. (2) The outgoing sum is 0.35 + 0.25 = 0.60, so the diagonal entry is 1 - 0.60 = 0.40; the row meets the condition because the outgoing sum is at most 1 and both increments are nonnegative. (3) The probability also depends on the probability distribution just before that update and on all earlier paths, including prior entry into state 2 and transitions out of it. (4) probtrans propagates model-based fitted transition hazards, whereas Aalen–Johansen estimates occupancy empirically from observed transition histories.