Why Numerical Reliability Matters
Panel data record a process at visits rather than continuously. A model therefore relates the states seen at two visits to the transition probabilities over the elapsed time between them. Those probabilities may be straightforward to define but difficult to calculate repeatedly inside an optimizer, especially when there are several states, unequal visit times, or rates near zero. A numerically unstable calculation can make a sound model appear not to converge; conversely, a convergence message can accompany a poor or boundary solution.
The goal here is to assess the numerical behavior of a panel-data transition-model fit, not to change the scientific model until an optimizer succeeds. Numerical checks should be made alongside the model assumptions: the chosen state space and transition structure must be appropriate, the intensity matrix must obey its constraints, and the observation-interval likelihood must represent the recorded panel data. A stable computation cannot repair an incorrect likelihood or an unidentifiable parameterization.
Model Quantities and Numerical Constraints
For a continuous-time Markov model with \(K\) states and time-constant transition intensities over an interval, let \(Q\) be the \(K\times K\) intensity matrix. Its off-diagonal entry \(q_{rs}\) is the instantaneous transition rate from state \(r\) to state \(s\), so \(q_{rs}\geq 0\) for \(r\ne s\); its diagonal entry is \(q_{rr}=-\sum_{s\ne r}q_{rs}\). Thus every row sums to zero. Over an interval of length \(\Delta\), the transition-probability matrix is the matrix exponential \(P(\Delta)=\exp(Q\Delta)\). It is a matrix function, not an entry-by-entry exponential.
These constraints matter computationally as well as statistically. Directly optimizing every entry of \(Q\) can produce negative rates or rows that do not sum to zero. A common parameterization sets each off-diagonal rate to \(q_{rs}=\exp(\eta_{rs})\) and constructs the diagonal from the row-sum condition. This guarantees positive rates at finite \(\eta_{rs}\), but it cannot attain a rate of exactly zero: a boundary rate is approached only as \(\eta_{rs}\) tends to negative infinity. An alternative is constrained optimization on \(q_{rs}\geq0\), which can represent exact zeros but requires appropriate boundary handling.
For a fixed model, calculate \(\exp(Q\Delta)\) with a numerical matrix-exponential routine, typically one based on scaling and squaring with a Padé approximation. Do not replace it with entrywise exponentiation or assume an eigenvector decomposition will be well-conditioned: matrices with repeated or nearly repeated eigenvalues can make that route inaccurate. For large matrices or long intervals, uniformization is another option. Choose \(\nu\geq\max_r(-q_{rr})\), set \(R=I+Q/\nu\), and express the exponential as a Poisson-weighted sum of powers of the stochastic matrix \(R\). Truncation can be controlled by the omitted Poisson tail.
Worked Example: Checking a Two-State Matrix Exponential
Suppose a hypothetical process moves between state 0 and state 1 with rates \(a=0.4\) per year from 0 to 1 and \(b=0.2\) per year from 1 to 0. For an interval of \(\Delta=3\) years, form \(Q\), calculate \(P(3)\), and check that the result is a probability matrix.
The rate matrix is \(Q=\begin{pmatrix}-0.4&0.4\\0.2&-0.2\end{pmatrix}\), with rows summing to zero. For this two-state model, \(P_{00}(\Delta)=\frac{b}{a+b}+\frac{a}{a+b}e^{-(a+b)\Delta}\), \(P_{01}(\Delta)=\frac{a}{a+b}(1-e^{-(a+b)\Delta})\), \(P_{10}(\Delta)=\frac{b}{a+b}(1-e^{-(a+b)\Delta})\), and \(P_{11}(\Delta)=\frac{a}{a+b}+\frac{b}{a+b}e^{-(a+b)\Delta}\).
Here \((a+b)\Delta=0.6(3)=1.8\), and \(e^{-1.8}\approx0.1653\). Substitution gives \(P_{00}(3)=\frac{1}{3}+\frac{2}{3}(0.1653)\approx0.4435\), \(P_{01}(3)=\frac{2}{3}(1-0.1653)\approx0.5565\), \(P_{10}(3)=\frac{1}{3}(1-0.1653)\approx0.2782\), and \(P_{11}(3)=\frac{2}{3}+\frac{1}{3}(0.1653)\approx0.7218\). Thus \(P(3)\approx\begin{pmatrix}0.4435&0.5565\\0.2782&0.7218\end{pmatrix}\). Both rows sum to one to the displayed precision, and all entries are nonnegative.
The matrix exponential is essential here. Exponentiating each entry of \(Q\Delta\) would not produce the transition probabilities and would not, in general, give rows summing to one. This small model also gives a useful check against a software calculation: independently evaluate the closed-form entries, then compare them with the numerical matrix exponential.
Optimization, Convergence, and Multiple Starts
An optimizer searches for parameter values that maximize the panel log-likelihood. A convergence flag usually means that a numerical stopping rule was met; it does not prove that the reported point is a global maximum, that the likelihood was evaluated accurately, or that the data identify every rate. Inspect the fitted rates, the log-likelihood, the gradient or score, and whether any constraints are active. If the model has several free rates, repeat the fit from dispersed starting values. Agreement across runs is evidence of stability, not a mathematical proof of a global maximum.
A simple example shows how multiple starts can be assessed against an exact calculation. Suppose a hypothetical one-way process has states 0 and 1, with a single unknown rate \(q\) per year from 0 to 1 and no return transition. Twenty independent subjects are observed in state 0 at baseline and again after two years; 14 are in state 1 at follow-up and 6 remain in state 0. Assume the one-way, constant-rate model and the stated observation intervals are appropriate. Then \(P_{00}(2)=e^{-2q}\), \(P_{01}(2)=1-e^{-2q}\), and the log-likelihood is
Differentiating gives \(\ell'(q)=-12+\frac{28}{e^{2q}-1}\). Setting this to zero yields \(e^{2q}-1=28/12=7/3\), so \(e^{2q}=10/3\) and \(\hat q=\frac{1}{2}\log(10/3)\approx0.6020\) per year. The second derivative is \(\ell''(q)=-\frac{56e^{2q}}{(e^{2q}-1)^2}<0\), so this stationary point is a maximum. Substitution gives \(\ell(\hat q)=6\log(0.3)+14\log(0.7)\approx-12.22\). This exact solution supplies a benchmark for the optimizer.
Worked Example: Comparing Fits From Multiple Starting Values
Use the hypothetical data and model just described. Fit the model from starting rates \(q=0.05\), \(q=0.5\), and \(q=2.0\) per year, using a stable likelihood calculation and a positive-rate parameterization.
Step 1: Check the calculation. For each proposed \(q\), evaluate \(P_{00}(2)=e^{-2q}\) and \(P_{01}(2)=1-e^{-2q}\), then form \(\ell(q)=6\log P_{00}(2)+14\log P_{01}(2)\). The transition probabilities must lie between zero and one and sum to one.
Step 2: Compare the results. In this one-parameter example the likelihood is strictly concave for \(q>0\), as shown by \(\ell''(q)<0\). A correctly implemented optimizer started at each of the three rates should therefore reach the same maximum, approximately \(q=0.6020\) per year, with log-likelihood approximately \(-12.22\). Small differences in the final digits can arise from stopping tolerances.
Step 3: Interpret disagreement. If one run reports a materially different rate or likelihood, do not select the run with the most attractive estimate. Check the likelihood implementation, the matrix-exponential routine, parameter scaling, and convergence diagnostics; then rerun with tighter tolerances and additional starts. The analytic maximum makes clear that, in this example, materially different converged solutions indicate a computational or implementation problem rather than multiple genuine maxima.
Boundary Estimates and Weak Information
Rates can be estimated at or near zero, and other constraints in more elaborate models can also be active at a boundary. A near-zero rate may reflect genuinely little evidence of that transition, weak identification from the visit schedule, or a numerical failure. These explanations are not interchangeable. In particular, a small estimated rate is not by itself proof that a transition is absent. Examine the likelihood near the boundary, the observed transition patterns, the sensitivity to starting values, and whether the model structure permits the rate to be distinguished from other rates.
Worked Example: Diagnosing a Zero-Rate Boundary
Consider a hypothetical one-way model with rate \(q\geq0\) per year. Twenty subjects are observed in state 0 at baseline and two years later, and all 20 are still in state 0. Assume the observation intervals are correct and the one-way constant-rate model is the intended likelihood.
The likelihood is \(L(q)=[P_{00}(2)]^{20}=[e^{-2q}]^{20}\), so \(\ell(q)=-40q\). For every \(q>0\), \(\ell(q)<0\), while \(\ell(0)=0\). Therefore the maximum over \(q\geq0\) is at the boundary \(\hat q=0\) per year. This is consistent with the data under this model, but it does not establish that the true rate is exactly zero; the panel observations provide limited information about rare transitions.
If the rate is parameterized as \(q=e^\eta\), no finite \(\eta\) gives \(q=0\). An optimizer may keep decreasing \(\eta\), report an extremely small positive \(q\), or stop because its step size is small. That behavior can be a boundary signal rather than evidence of an interior optimum. A constrained fit that allows \(q=0\) can represent the maximum exactly. Compare the fitted likelihood with its boundary value, inspect the score or constraint-aware optimality conditions, and report the boundary estimate rather than implying that a tiny positive number has meaningful precision.
Diagnostics, Alternatives, and Limits
Before trusting a fit, check that \(Q\) has nonnegative off-diagonal entries and rows summing to zero, and that every calculated \(P(\Delta_i)\) has entries within numerical tolerance of the interval from zero to one and rows summing to one. Check several interval lengths, including the largest in the data, and compare selected results with an independent calculation where feasible. A tiny negative entry can arise from rounding, but it should prompt a check of the computation; silently clipping entries to zero can alter the likelihood and conceal a poorly conditioned calculation. Likewise, renormalizing rows may conceal a bug or an invalid rate matrix. If small corrections are used, quantify and document them rather than treating them as harmless by default.
Check the optimizer separately from the matrix exponential. Refit from multiple dispersed starting values; compare final log-likelihoods, estimates, gradient norms, iteration limits, and active constraints. Rescale parameters or time units if magnitudes differ greatly, and use tighter tolerances to test whether the reported estimates change materially. Underflow in very small transition probabilities can make a likelihood contribution appear to be zero; compute log probabilities carefully and investigate whether the model assigns effectively impossible probability to observed transitions. A flat likelihood, a singular or ill-conditioned information matrix, or very different parameter estimates with nearly identical likelihoods can signal weak identification rather than failed optimization.
The main numerical alternatives serve different purposes. A general-purpose matrix-exponential routine is convenient for moderate state spaces and many interval lengths. Uniformization offers a probability-based sum with a truncation error that can be controlled, though it may be costly when the Poisson mean is large. Closed forms are excellent checks for simple models but do not extend conveniently to general state spaces. Constrained rate optimization handles exact zero rates directly; log-rate optimization enforces positivity smoothly but can behave poorly near zero. Multiple starts help detect dependence on initialization, but they do not establish uniqueness or overcome nonidentifiability. No computational approach removes the need to check the model and observation assumptions.
Check Your Understanding
Answer each question, then compare your reasoning with the solutions.
- Why is exponentiating every entry of \(Q\Delta\) not a valid way to obtain \(P(\Delta)\)?
- In the one-way example with 6 subjects remaining in state 0 and 14 transitioning over two years, what is the maximum-likelihood rate, in units of inverse years?
- If all subjects remain in state 0 in the one-way model, why might an optimizer using \(q=e^\eta\) fail to report an exact zero rate?
- Name two checks that help distinguish a numerical matrix-exponential problem from an optimization problem.
Solutions. (1) The transition matrix is the matrix function \(\exp(Q\Delta)\), whose entries account for paths through the states. Entrywise exponentiation does not generally produce a matrix with nonnegative entries and rows summing to one. (2) The score equation gives \(\hat q=\frac{1}{2}\log(10/3)\approx0.6020\) per year. (3) For every finite \(\eta\), \(e^\eta>0\); the boundary value is approached only as \(\eta\) tends to negative infinity. (4) Check that \(P(\Delta)\) has nonnegative entries and row sums of one, and compare it with a closed-form or independent matrix-exponential calculation. Separately, inspect optimizer gradients, likelihood values across starts, and constraint activity.