Numerical Integration

Why numerical integration is used

In health economics, it converts survival curves into expected life-years, accumulates time-varying costs and utilities, and averages model outputs over probability distributions. This page explains the main approximation rules, their error behaviour, and the validation needed when an integral becomes a decision-model output.

Numerical integration is an approximation method, not a substitute for defining the correct quantity. The integration limits, integrand, time scale, discounting and interpolation rule must all match the clinical and economic question.

A definite integral represents accumulated quantity

The definite integral of a function (f(x)) from (a) to (b) represents signed area or accumulated contribution over that interval. Numerical methods replace the continuous area with a weighted sum of function values. Smaller subintervals or higher-order rules can improve accuracy, but only when the function is evaluated and represented appropriately.

The target quantity is:

$$I=\int_a^b f(x)\,dx$$

A general quadrature approximation is:

$$I\approx\sum_{i=0}^{n}w_i f(x_i)$$

where (x_i) are evaluation points and (w_i) are method-specific weights. The rule determines where the function is evaluated and how much each value contributes.

The rectangle rules use one height per interval

Rectangle rules approximate each subinterval by a rectangle. A left rule uses the function value at the start, a right rule uses the value at the end and a midpoint rule uses the value halfway through. These methods are simple and transparent, but left and right rules can be biased in opposite directions for monotonic functions.

For equal step size (h=(b-a)/n), the left rule is:

$$I_{L}=h\sum_{i=0}^{n-1}f(x_i)$$

The right rule is:

$$I_{R}=h\sum_{i=1}^{n}f(x_i)$$

The midpoint rule is:

$$I_{M}=h\sum_{i=0}^{n-1}f\left(\frac{x_i+x_{i+1}}{2}\right)$$

For a decreasing survival function, the left rule generally overestimates area and the right rule generally underestimates it. The midpoint rule is usually more accurate for a smooth function because it uses information from the centre of each interval.

The trapezoidal rule connects adjacent observations

The trapezoidal rule assumes that the function changes linearly between adjacent points. It is widely used for survival curves and other time-indexed data because it can accommodate unequal intervals and has a clear geometric interpretation. Its validity depends on whether straight-line interpolation is reasonable between observations.

For nodes (x_0,\ldots,x_n), the composite trapezoidal rule is:

$$I_T=\sum_{i=0}^{n-1} \frac{x_{i+1}-x_i}{2} \left[f(x_i)+f(x_{i+1})\right]$$

With equal spacing (h), this becomes:

$$I_T=h\left[ \frac{f(x_0)}{2}+ \sum_{i=1}^{n-1}f(x_i)+ \frac{f(x_n)}{2} \right]$$

The trapezoidal rule is exact for a linear function. For a smooth curved function, its error decreases as the intervals become narrower, but abrupt changes or long gaps may require more detailed modelling rather than mechanical refinement.

A worked survival-area example

Expected survival over a fixed horizon equals the area under the survival function. Suppose estimated survival is 1.00 at baseline, 0.80 at one year and 0.50 at two years. The trapezoidal rule provides an auditable approximation of restricted mean survival time over two years.

The first-year area is:

$$\frac{1}{2}(1.00+0.80)\times1=0.90$$

The second-year area is:

$$\frac{1}{2}(0.80+0.50)\times1=0.65$$

The total is:

$$RMST(2)\approx0.90+0.65=1.55\;\text{life-years}$$

This value is the area under the piecewise-linear survival curve from zero to two years. It does not include survival beyond two years and should not be interpreted as unrestricted life expectancy.

Simpson's rule uses local quadratic interpolation

Simpson's rule approximates a smooth function with quadratic segments rather than straight lines. It can achieve greater accuracy than the trapezoidal rule for the same equally spaced nodes. The standard composite form requires an even number of subintervals.

For equal spacing (h) and even (n):

$$I_S=\frac{h}{3}\left[ f(x_0)+f(x_n) +4\sum_{\substack{i=1\\i\text{ odd}}}^{n-1}f(x_i) +2\sum_{\substack{i=2\\i\text{ even}}}^{n-2}f(x_i) \right]$$

Simpson's rule is exact for polynomials up to degree three under its standard conditions. It is not automatically preferable for step functions, discontinuities, heavily censored tails or irregularly spaced empirical data.

Gaussian quadrature chooses efficient nodes and weights

Gaussian quadrature selects evaluation points and weights to integrate polynomials of a high degree exactly with relatively few function evaluations. It is useful when the integrand is smooth and can be evaluated at the required internal nodes. The familiar Gauss-Legendre rule applies to an unweighted integral over a finite interval after transformation to ([-1,1]).

An (n)-point Gauss-Legendre rule has the form:

$$\int_{-1}^{1}f(x)\,dx \approx \sum_{i=1}^{n}w_i f(x_i)$$

where the nodes (x_i) are roots of the (n)th Legendre polynomial. For integration over ([a,b]), transform:

$$x=\frac{b-a}{2}t+\frac{a+b}{2}$$

so that:

$$\int_a^b f(x)\,dx =\frac{b-a}{2} \int_{-1}^{1} f\left(\frac{b-a}{2}t+\frac{a+b}{2}\right)dt$$

Other Gaussian rules use weights tailored to infinite intervals or particular probability distributions. The selected rule must match the domain and weighting function rather than being chosen only for computational speed.

Adaptive quadrature concentrates effort where it is needed

Adaptive methods estimate local error and subdivide intervals where the function is difficult to approximate. Smooth sections receive fewer evaluations, while sections with rapid curvature or other challenging behaviour receive more. This can provide accurate results without imposing a uniformly tiny step size over the full domain.

A typical adaptive algorithm follows this logic:

1. Approximate the integral over an interval using a chosen quadrature rule.
2. Divide the interval and calculate a refined approximation.
3. Compare the coarse and refined values to estimate local error.
4. Accept the interval or subdivide it further according to absolute and relative tolerances.
5. Sum accepted interval contributions and report the estimated numerical error.

Adaptive refinement can still fail when the integrand has an undisclosed singularity, discontinuity, severe oscillation or very narrow peak. The analyst should inspect the function and not rely solely on a success flag returned by software.

Numerical integration converts survival into life-years

For a non-negative survival time (T), expected survival over horizon (\tau) can be expressed as the area under the survival curve. This links numerical integration directly to life-year estimation. The horizon and extrapolated tail assumptions are often more consequential than the quadrature rule itself.

Restricted mean survival time is:

$$RMST(\tau)=\int_0^{\tau}S(t)\,dt$$

where (S(t)=P(T>t)). Incremental restricted mean survival is:

$$\Delta RMST(\tau)= \int_0^{\tau} \left[S_1(t)-S_0(t)\right]dt$$

When a parametric survival model is used beyond observed follow-up, numerical integration includes both empirical and extrapolated contributions. Analysts should report how much of the estimated area arises beyond the observed data and test alternative survival functions.

Time-varying utility produces quality-adjusted survival

Quality-adjusted life-years can be represented as the integral of survival probability multiplied by expected health-related utility over time. This continuous-time form makes the timing of morbidity and death explicit. Utility and survival processes must be combined consistently, without applying survival twice or ignoring dependency between health state and mortality.

For expected utility (u(t)):

$$QALY(\tau)=\int_0^{\tau}S(t)u(t)\,dt$$

If utility is defined conditionally on being alive, multiplication by (S(t)) converts it to an unconditional expected contribution. If utility estimates already include zero values for death in the population mean, multiplying by survival again would double count mortality.

Discounting can be included inside the integral

Continuous-time models can accumulate discounted costs or health outcomes by multiplying the time-varying quantity by a discount function. The discount convention must match the jurisdictional guidance and the model's time unit. Mixing continuous discounting with an effective annual discrete rate without conversion changes the result.

With continuous discount rate (r), the discount factor at time (t) is:

$$d(t)=e^{-rt}$$

The present value of a time-varying cost rate (c(t)) is:

$$PV_C(\tau)=\int_0^{\tau}c(t)e^{-rt}\,dt$$

If an annual effective discrete rate (r_d) is converted to a continuous rate, then:

$$r_c=\ln(1+r_d)$$

The discount function should be evaluated at the time the cost or outcome occurs, not automatically at an interval endpoint.

Expected values are integrals over probability distributions

An expected model output can be expressed as an integral of the output over the uncertain parameter distribution. Numerical integration is useful when the parameter dimension is small and the model is smooth. In higher dimensions, simulation or specialised sparse-grid methods may be more practical because ordinary grids grow rapidly with the number of parameters.

For model output (g(\theta)) and parameter density (p(\theta)):

$$E[g(\theta)]= \int g(\theta)p(\theta)\,d\theta$$

For one parameter, deterministic quadrature may estimate this expected value with fewer model evaluations than Monte Carlo simulation. For many correlated parameters, the multidimensional integral becomes harder and the integration method must preserve the joint distribution.

Step size controls discretisation error

Discretisation error arises because a continuous area is represented with a finite number of evaluations. A convergence check repeats the calculation with progressively finer resolution and examines whether the result stabilises. The reported tolerance should be meaningful relative to the decision, not merely small in absolute numerical terms.

For a smooth function and equal spacing (h), the composite trapezoidal rule has global error of order:

$$E_T=O(h^2)$$

while composite Simpson's rule has global error of order:

$$E_S=O(h^4)$$

These orders describe asymptotic behaviour under regularity conditions. They do not guarantee accuracy when the function is discontinuous, poorly evaluated, extrapolated incorrectly or subject to large data error.

A practical refinement check may compare estimates at step sizes (h) and (h/2):

$$\Delta_h=|I_{h/2}-I_h|$$

The result should be refined until (\Delta_h) is below a prespecified tolerance and the remaining numerical uncertainty is negligible for the decision.

Interpolation determines the area between observed points

When integration begins with discrete observations, the chosen interpolation rule defines the unobserved path between them. Trapezoidal integration implies linear interpolation, while a step rule carries one observation forward or backward. Different choices can materially affect costs, survival or utility when measurement intervals are long.

Common choices include:

• Linear interpolation assumes a constant slope between adjacent observations.
• Left-continuous steps apply the earlier value across the interval.
• Right-continuous steps apply the later value across the interval.
• Spline interpolation creates a smooth curve but can overshoot natural bounds.
• Model-based interpolation uses a fitted statistical or mechanistic function between observations.

The interpolation method should reflect how the underlying process evolves. A smooth curve is not necessarily more accurate than a clinically justified step change.

Infinite and improper integrals require tail control

Some expected-value and survival calculations extend to infinity or contain difficult boundary behaviour. Numerical software transforms or truncates these domains, but the analyst must verify that the omitted tail contributes negligibly or is modelled correctly. A finite integration result does not prove that the underlying improper integral exists.

For survival:

$$E[T]=\int_0^{\infty}S(t)\,dt$$

In practice, the integral may be truncated at (\tau) and separated into observed and tail contributions:

$$E[T]= \int_0^{\tau}S(t)\,dt+ \int_{\tau}^{\infty}S(t)\,dt$$

The tail should be checked against analytic results when available, alternative extrapolation models and increasing truncation points. Arbitrary truncation can understate life-years, costs or event probabilities even when the quadrature error within ([0,\tau]) is tiny.

Numerical error and model uncertainty are different

Numerical error is the discrepancy created by the integration algorithm, finite step size and floating-point arithmetic. Model uncertainty concerns whether the integrand, parameters, extrapolation and structure correctly represent the real process. Reducing numerical error cannot correct a biased survival model or an inappropriate utility function.

The main uncertainty layers are:

Reporting only the quadrature tolerance can create false precision when the dominant uncertainty comes from the model or data.

Common implementation errors

Numerical integration can produce a stable and precise-looking number even when the wrong function or domain has been integrated. Validation must therefore connect the numerical calculation to the estimand and model. The following errors are especially relevant in health-economic work.

• Using the wrong time unit can apply annual costs or discount rates to monthly or daily integration steps without conversion.
• Omitting interval widths treats a sum of function values as an area.
• Applying a rule designed for equal spacing to irregular data assigns incorrect weights.
• Double counting interval endpoints can overstate cumulative outcomes when separate segments are added.
• Applying survival twice can understate QALYs when utilities already include death.
• Using matrix powers element by element does not represent repeated transitions.
• Treating a small quadrature error as total uncertainty ignores parameter, structural and extrapolation uncertainty.
• Truncating an infinite horizon too early can omit a material tail contribution.
• Allowing interpolation to violate natural bounds can create survival above one, negative hazards or utilities outside the intended range.
• Comparing methods on different nodes or horizons confounds the numerical rule with the data supplied to it.

A practical validation sequence

Validation should start with an analytical or transparent benchmark and then test convergence under increasingly accurate numerical settings. The contribution from each interval should remain traceable so that unexpected area can be located. The following sequence provides an auditable minimum.

1. Define the integrand, limits, units and interpretation before choosing an algorithm.
2. Calculate a function with a known exact integral to confirm the implementation and weighting rule.
3. Verify node ordering and interval widths, especially when observations are irregularly spaced.
4. State the interpolation assumption used between observed or modelled points.
5. Repeat the calculation with a finer grid or stricter tolerance and compare the results.
6. Split the domain at discontinuities, knots or known rapid changes rather than forcing one smooth rule across them.
7. Reconcile interval contributions to the reported total using unrounded values.
8. Check natural bounds and monotonicity for survival, probabilities, utilities and cumulative quantities.
9. Test longer horizons and alternative tail models when the integral includes extrapolation or an infinite domain.
10. Propagate parameter and structural uncertainty separately from numerical integration error.
11. Confirm that remaining numerical error cannot change the decision at the required reporting precision.

Numerical integration is most valuable when it turns a clearly defined continuous quantity into a reproducible model output. Accuracy depends as much on the integrand, interpolation and horizon as on the quadrature formula itself.