Sparse-Grid Quadrature for Higher-Dimensional Effects

Why Sparse Grids Are Useful

A subject’s contribution to a joint-model likelihood is usually obtained by integrating a conditional contribution over that subject’s random effects. If the random-effects vector has several components, a tensor-product quadrature rule can become expensive quickly: with \(m\) nodes in each of \(d\) dimensions, it requires \(m^d\) function evaluations. Sparse-grid quadrature reduces this cost by combining selected lower-dimensional rules rather than using every combination of nodes.

The reduction is useful, not automatic proof of accuracy. Sparse grids work particularly well when the transformed integrand is reasonably smooth and its important structure can be represented by contributions involving relatively few coordinates. Strong interactions or a sharply concentrated integrand can be missed at a low grid level. The approximation therefore needs to be checked by increasing the grid level and examining the behavior of the resulting likelihood contribution.

From a Random-Effects Integral to a Sparse Grid

Let \(b_i\) be the \(d\)-dimensional random-effects vector for subject \(i\), with \(b_i\sim N(0,D)\), where \(D\) is positive definite. Write \(g_i(b_i;\theta)\) for the subject’s conditional joint contribution given the random effects, including the observed longitudinal and event-time data. The marginal contribution is

Here \(\theta\) collects the model parameters, and \(\phi_d(\,\cdot\,;0,D)\) is the multivariate normal density. Factor \(D=CC^\mathsf{T}\), for example by a Cholesky decomposition, and set \(b_i=Cz\), where \(z\sim N(0,I_d)\). The change of variables converts the integral to a standard-normal expectation:

This transformation is valid when \(D\) is positive definite and the integral exists. It is especially helpful because a product standard-normal measure can be approximated using one-dimensional Gaussian quadrature rules. The same transformation and covariance convention are used for tensor-product quadrature; sparse grids change how those one-dimensional rules are combined.

The Smolyak Construction

Let \(Q_\ell\) denote a one-dimensional quadrature rule at level \(\ell\), and set \(Q_0=0\). For instance, \(Q_\ell f\) is a weighted sum of evaluations of \(f\) at the nodes for that level. Define the difference rule \(\Delta_\ell=Q_\ell-Q_{\ell-1}\). With dimension \(d\) and level parameter \(q\), one common Smolyak convention is

The product \(\otimes\) means that the one-dimensional rules are applied across the indicated coordinates. The level convention matters: other references may index Smolyak rules differently, so a reported level should always be accompanied by the convention or rule definition. Higher \(q\) includes more levels and generally improves the approximation, at greater computational cost.

For \(d=2\) and \(q=2\), the allowable level pairs are \((1,1)\), \((1,2)\), and \((2,1)\). Since \(\Delta_1=Q_1\) and \(\Delta_2=Q_2-Q_1\), expansion gives

This identity is a short derivation of the combination form: the two terms with one level-2 rule each add the refinements along either coordinate, while the level-1 product is subtracted because it was counted in both. The subtraction is essential; simply adding the two refined products would double-count their shared coarse contribution.

A Numerical Likelihood-Integral Calculation

The following hypothetical calculation uses a smooth, nonnegative integrand to illustrate the approximation itself. Take independent standard-normal coordinates \(Z_1,Z_2\), and let \(h(z_1,z_2)=\exp(-0.1z_1^2-0.2z_2^2)\). This function is at most one, so it can serve as a simple illustration of a conditional contribution; general likelihood integrands need not be bounded by one.

Use the one-node rule \(Q_1f=f(0)\), and the three-node standard-normal rule \(Q_2\), with nodes \(0,\sqrt{3},-\sqrt{3}\) and respective weights \(2/3,1/6,1/6\). This rule integrates polynomials through degree five exactly. Applying the two-dimensional combination formula gives

For the first term, the second coordinate is zero, so the evaluations are \(1\) at \(z_1=0\) and \(\exp(-0.3)\) at each of \(z_1=\pm\sqrt{3}\). Thus \((Q_2\otimes Q_1)h=2/3+(1/3)\exp(-0.3)\). Similarly, \((Q_1\otimes Q_2)h=2/3+(1/3)\exp(-0.6)\), and the coarse term is one. Substitution yields

For comparison, the exact expectation factors because the coordinates are independent. For \(Z\sim N(0,1)\), direct Gaussian integration gives \(E\{\exp(-aZ^2)\}=1/\sqrt{1+2a}\) for \(a\geq0\). Therefore,

The sparse approximation is about \(0.008307\) below this exact value. It is a useful approximation in this illustration, but its error is not zero. In an actual joint model, the exact expectation is not usually available; comparing neighboring sparse-grid levels is one practical way to assess numerical stability.

When a Sparse Grid Misses an Interaction

The low-level example above puts its nonzero evaluations on the coordinate axes: one random effect varies while the other is fixed at zero. Consequently, it provides little information about how the integrand behaves when both coordinates are away from zero. This can matter in joint models, where correlated random effects and association structures can produce substantial joint variation.

Checking the Approximation and Choosing an Alternative

A practical sparse-grid calculation should be treated as a numerical approximation whose stability must be evaluated. Useful checks include:

• Increase the level. Recalculate selected subject contributions at a higher \(q\). Compare the values and, importantly, the total log likelihood; small changes in many contributions can accumulate.
• Inspect the scale and positivity. Sparse-grid combination weights can be negative, even when the one-dimensional quadrature weights are positive. The resulting sum approximates a positive integral but may be nonpositive from numerical error or an inadequate grid. A nonpositive likelihood contribution cannot be logged and calls for investigation, not silent truncation.
• Check for concentrated or irregular integrands. Plot or evaluate the transformed integrand along coordinate directions and selected combinations of coordinates. A narrow peak, steep tails, or strong cross-coordinate pattern may require a finer grid or a different integration strategy.
• Compare where feasible. In a low-dimensional check, compute a more extensive tensor-product approximation. Agreement across methods and levels is reassuring, though it does not prove the approximation error is negligible.

Tensor-product quadrature is the natural direct alternative: it is simple to define and can be highly accurate at moderate dimension, but its node count grows exponentially with dimension. Sparse grids trade this exhaustive coverage for a structured selection of combinations. Monte Carlo integration avoids a deterministic node grid and can be useful at high dimension, but it introduces sampling error. No method is uniformly best: the choice depends on dimension, smoothness, concentration, computational budget, and the accuracy needed for estimation.

Sparse-grid accuracy is not guaranteed merely by smoothness, and there is no universal error bound that can be inferred from the level alone. The most favorable situation is a smooth transformed integrand whose important behavior is captured by low-order coordinate interactions. A low-level sparse grid may perform poorly for strong interactions, narrow posterior concentration, or substantial tail contributions. If refinement does not stabilize the likelihood, using more nodes in the same sparse-grid family may not be enough; the transformation, integration method, or model computation may need reconsideration.

Common Pitfall: Treating the Sparse Sum as a Probability Average

A common mistake is to assume that a sparse-grid estimate must be a weighted average with nonnegative weights summing to one, just like a probability expectation. The combination formula includes subtraction, so some combined weights can be negative. The correction is to view the sparse sum as a quadrature approximation to an integral, not as a discrete probability distribution. Check the computed integral’s positivity and its stability under refinement, and do not interpret individual sparse-grid weights as probabilities.