Numerical Integration of Cause-Specific Exit Probabilities

Why Numerical Integration Is Needed

When cause-specific decrement intensities vary with time, a cause’s probability of exit over an interval depends on two changing quantities: the intensity of that cause and the probability that a member of the group is still exposed to it. In general, the resulting cumulative incidence does not have a convenient closed form. Numerical integration provides a systematic approximation.

The key is to integrate the survival-weighted intensity, not the cause-specific intensity alone. An intensity can be high late in an interval, but relatively few members may remain in the group by then. The method below combines time-varying cause-specific intensities with the probability of remaining in the group.

Model and assumptions: At time 0, a group is in the active state. Members may leave through one of a finite set of mutually exclusive competing causes, and only the first exit is considered. Cause \(j\) has a deterministic, nonnegative intensity \(\mu_j(t)\) at time \(t\), conditional on still being active immediately before \(t\). There are no entries into the group during the interval.

Write the total intensity as \(\mu_{\text{tot}}(t)=\sum_j\mu_j(t)\), and let \(S(t)\) be the probability of remaining active through time \(t\), with \(S(0)=1\). The cumulative probability of exiting from cause \(j\) by time \(T\), also called the cumulative incidence for cause \(j\), is

F_j(T)=\int_0^T S(t)\mu_j(t)\,dt.

This formula assumes the intensities describe the first exit from the active state. It does not calculate the probability of experiencing cause \(j\) after first leaving through another cause.

Why the Integral Uses Survival Weighting

Over a short interval from \(t\) to \(t+\Delta t\), a member who is still active at \(t\) has probability approximately \(\mu_j(t)\Delta t\) of leaving from cause \(j\), conditional on being active. The probability of being both active at \(t\) and exiting from cause \(j\) in that short interval is therefore approximately \(S(t)\mu_j(t)\Delta t\). Dividing by \(\Delta t\) and taking the limit gives the cause-specific exit-time density.

f_j(t)=S(t)\mu_j(t).

Integrating this density from 0 to \(T\) gives the probability of a cause-\(j\) exit by time \(T\). The remaining probability follows from the survival differential equation:

\frac{dS(t)}{dt}=-S(t)\mu_{\text{tot}}(t), \qquad S(t)=\exp\left(-\int_0^t\mu_{\text{tot}}(u)\,du\right).

Thus, even when each \(\mu_j(t)\) is specified explicitly, the integrand for \(F_j(T)\) includes an exponential of the accumulated total intensity. Numerical quadrature approximates the integral of that full product.

Numerical Integration Methods

Choose a grid \(0=t_0<t_1<\cdots<t_n=T\). At each grid point, evaluate \(g_j(t_k)=S(t_k)\mu_j(t_k)\). Then apply a quadrature rule to the values of \(g_j\). For an interval of width \(h\), the trapezoidal rule uses the average of the integrand at its two endpoints; the midpoint rule uses the integrand at the interval midpoint. For equally spaced points, Simpson’s rule uses three consecutive points and is most directly applicable when the number of subintervals is even.

\begin{aligned} \text{Trapezoidal rule:}\quad &\int_a^b g(t)\,dt\approx \frac{b-a}{2}\bigl(g(a)+g(b)\bigr),\\ \text{Midpoint rule:}\quad &\int_a^b g(t)\,dt\approx (b-a)g\left(\frac{a+b}{2}\right),\\ \text{Simpson's rule:}\quad &\int_a^b g(t)\,dt\approx \frac{b-a}{6} \left[g(a)+4g\left(\frac{a+b}{2}\right)+g(b)\right]. \end{aligned}

For a time-varying model, calculate \(S(t)\) at the points required by the chosen rule using the accumulated total intensity up to each point. If the cumulative total intensity has no closed form, approximate that integral as well, refining the grid if needed. Applying a rule directly to \(\mu_j(t)\), without the survival factor, generally overstates the probability of exit from cause \(j\).

Worked Example: Trapezoidal Approximation Across Two Intervals

A hypothetical decrement model has two causes, A and B. On \(0\le t\le1\), their intensities vary linearly from \(\mu_A(0)=0.02\), \(\mu_A(1)=0.04\), \(\mu_B(0)=0.01\), and \(\mu_B(1)=0.02\). On \(1\le t\le2\), they vary linearly from \(\mu_A(1)=0.04\), \(\mu_A(2)=0.06\), \(\mu_B(1)=0.02\), and \(\mu_B(2)=0.04\). Time is measured in years and the rates are per year. Approximate the probability of cause A exit by time 2 using the composite trapezoidal rule with unit-width intervals.

Method. Integrate \(g_A(t)=S(t)\mu_A(t)\) by the trapezoidal rule. Because each total intensity is linear within its interval, its accumulated intensity on that interval is the interval width times the average of its endpoint rates.

The total intensity is \(0.03\) at time 0, \(0.06\) at time 1, and \(0.10\) at time 2. The accumulated total intensity to time 1 is \(H(1)=1(0.03+0.06)/2=0.045\), so \(S(1)=e^{-0.045}=0.955997\). The additional accumulated intensity from 1 to 2 is \(1(0.06+0.10)/2=0.08\). Thus \(H(2)=0.045+0.08=0.125\), and \(S(2)=e^{-0.125}=0.882497\). Also, \(S(0)=1\).

The cause A weighted intensities at the grid points are \(g_A(0)=1(0.02)=0.020000\), \(g_A(1)=0.955997(0.04)=0.038240\), and \(g_A(2)=0.882497(0.06)=0.052950\).

\begin{aligned} F_A(2) &\approx \frac{1}{2}(0.020000+0.038240) +\frac{1}{2}(0.038240+0.052950)\\ &=0.029120+0.045595\\ &=0.074715. \end{aligned}

Rounded to four decimal places, the approximate probability is 0.0747. The calculation integrates the survival-weighted cause A intensity at both intervals. A common mistake is to trapezoid-integrate \(\mu_A(t)\) alone; that ignores the decreasing number remaining exposed to cause A.

Worked Example: Refining a Midpoint Approximation

For a hypothetical model on \(0\le t\le1\), let \(\mu_A(t)=0.04+0.02t\) and \(\mu_B(t)=0.02+0.01t\), with time in years and intensities per year. Approximate the cause A cumulative incidence first with one midpoint panel and then with two equal midpoint panels.

Method. The total intensity is \(0.06+0.03t\), so the cumulative total intensity to time \(t\) is \(H(t)=0.06t+0.015t^2\), and \(S(t)=e^{-H(t)}\). For one panel, the midpoint is \(t=0.5\). Here \(H(0.5)=0.030+0.00375=0.03375\), \(S(0.5)=e^{-0.03375}=0.966813\), and \(\mu_A(0.5)=0.050\). The one-panel estimate is \(1(0.966813)(0.050)=0.048341\).

With two panels, the midpoints are \(0.25\) and \(0.75\), each with width \(0.5\). At \(t=0.25\), \(H(0.25)=0.015+0.0009375=0.0159375\), \(S(0.25)=0.984189\), and \(\mu_A(0.25)=0.045\), giving \(g_A(0.25)=0.984189(0.045)=0.044289\).

At \(t=0.75\), \(H(0.75)=0.045+0.0084375=0.0534375\), \(S(0.75)=0.947965\), and \(\mu_A(0.75)=0.055\), giving \(g_A(0.75)=0.947965(0.055)=0.052138\). Therefore,

F_A(1)\approx 0.5(0.044289+0.052138)=0.048213.

The two-panel estimate is 0.0482 to four decimal places. The one-panel estimate is 0.0483 to four decimal places. The difference illustrates why refinement is useful; it does not establish that either estimate is exact. The common mistake avoided here is using the cause A rate at the midpoint but treating the entire group as still active, rather than multiplying by \(S(t)\).

Worked Example: Simpson’s Rule and a Probability Check

For a third hypothetical model on \(0\le t\le2\), let \(\mu_A(t)=0.03+0.01t\) and \(\mu_B(t)=0.02+0.02t\), with time in years and intensities per year. Use Simpson’s rule with points 0, 1, and 2 to approximate both cause-specific cumulative incidences, then check their sum against the probability of any exit.

Method. The total intensity is \(0.05+0.03t\). Its accumulation is \(H(t)=0.05t+0.015t^2\). Hence \(S(0)=1\), \(S(1)=e^{-0.065}=0.937067\), and \(S(2)=e^{-0.16}=0.852144\).

For cause A, the weighted intensities at 0, 1, and 2 are respectively \(1(0.03)=0.030000\), \(0.937067(0.04)=0.037483\), and \(0.852144(0.05)=0.042607\). Since the full interval has width 2, Simpson’s rule gives

F_A(2)\approx \frac{2}{6} \left(0.030000+4(0.037483)+0.042607\right) =0.074179.

For cause B, the corresponding weighted intensities are \(1(0.02)=0.020000\), \(0.937067(0.04)=0.037483\), and \(0.852144(0.06)=0.051129\). Thus,

F_B(2)\approx \frac{2}{6} \left(0.020000+4(0.037483)+0.051129\right) =0.073686.

The sum of the cause-specific estimates is \(0.074179+0.073686=0.147865\). The probability of any exit by time 2 is \(1-S(2)=1-0.852144=0.147856\). These agree to four decimal places: both are 0.1479. The small difference before rounding is numerical quadrature error. Probability conservation is a useful check, but numerical estimates need not sum to the exact total exit probability unless the integrations are exact.

A common mistake is to use the intensity of each cause in the survival calculation separately. Survival depends on the total intensity from all competing causes; omitting cause B would make \(S(t)\) too large and distort both cumulative incidences.

Accuracy, Checks, and Exam Traps

Quadrature is an approximation, and its accuracy depends on how the integrand changes between grid points. A coarse grid may miss curvature in \(S(t)\mu_j(t)\), even if the cause-specific intensity itself appears nearly linear. Use more subintervals when rates change rapidly, when the interval is long, or when a probability check reveals a material discrepancy. For a fixed model, comparing estimates on successively finer grids is a practical diagnostic.

The total exit probability provides another check. Under the stated model,

\sum_jF_j(T)=1-S(T).

Exact cause-specific probabilities satisfy this identity. Separately approximated integrals may not satisfy it exactly, so compare the two sides at a precision suitable for the grid used. A large discrepancy can indicate an arithmetic error, an omitted competing cause in the total intensity, or an inadequate numerical grid.

Integrating an intensity instead of a probability density: \(\int_0^T\mu_j(t)\,dt\) is not generally the cause-\(j\) exit probability. The integrand must be \(S(t)\mu_j(t)\).
Using a cause-specific rate for survival: calculate \(S(t)\) using the sum of all competing intensities.
Mixing time units: if intensities are per year, use time in years in the integration. Convert rates or interval widths consistently.
Misapplying a quadrature rule: use each interval’s width and the correct endpoint or midpoint weights. Simpson’s rule requires equally spaced points for its standard composite form.
Treating an approximation as an exact identity: the conservation check is exact for true probabilities, but numerical quadrature introduces error.

Key takeaway: For each cause, approximate the integral of \(S(t)\mu_j(t)\). Compute survival from the accumulated total intensity, use a quadrature rule appropriate to the grid, and check the resulting cause probabilities against \(1-S(T)\).

Check Your Understanding

Use the stated numerical data and rules. Round probability answers to four decimal places unless the question asks for an expected count.

A group of 1,000 members has cause C weighted intensity values \(g_C(0)=0.015\), \(g_C(1)=0.024\), and \(g_C(2)=0.030\), in units per year. Use the trapezoidal rule on the two unit-width intervals and estimate the expected number of cause C exits. (Answer: 46.5 expected exits)
For a cause D, the survival-weighted intensity values at \(t=0,0.5,1.0,1.5\) are \(0.020,0.030,0.028,0.040\), respectively. Use the composite trapezoidal rule with width \(0.5\) to approximate its cumulative incidence by time 1.5. (Answer: 0.0440)
On \(0\le t\le2\), let \(\mu_A(t)=0.02+0.01t\) and \(\mu_B(t)=0.03+0.01t\), with time in years and rates per year. Use Simpson’s rule at 0, 1, and 2 to approximate \(F_A(2)\). (Answer: 0.0559)
For cause E, the weighted intensity values are \(g_E(0)=0.0400\), \(g_E(1)=0.1000\), and \(g_E(2)=0.0944\), per year. Apply the trapezoidal rule with unit-width intervals to estimate \(F_E(2)\). (Answer: 0.1672)

What You'll Learn