Importance Sampling for Random-Effects Integration

Why Importance Sampling Helps

In a joint model, a subject’s longitudinal and time-to-event observations are conditionally modeled given that subject’s random effects. The subject-specific marginal likelihood contribution requires integrating over those unobserved effects. When this integral is difficult to evaluate accurately by a fixed grid, importance sampling replaces the integral with a weighted average over random draws.

The key is not merely to draw random effects. The proposal distribution must put substantial probability wherever the subject’s data make the integrand substantial. In particular, it should cover the subject-specific posterior region, which can be far from the center of the population random-effects distribution when a subject’s observations are informative.

Definition: For subject \(i\), let \(b_i\) denote the random-effects vector, let \(y_i\) denote the longitudinal data, and let \((T_i,\delta_i)\) denote the observed event-time information. At a fixed parameter value \(\theta\), write \(h_i(b_i;\theta)=f(y_i,T_i,\delta_i\mid b_i;\theta)p(b_i;\theta)\), where the conditional likelihood and random-effects density are both included. The subject’s marginal likelihood contribution is \(L_i(\theta)=\int h_i(b_i;\theta)\,db_i\).

The conditional likelihood here is the one specified by the joint model, including its longitudinal and event-time components. The likelihood contribution also depends on the model’s treatment of censoring and observation. Importance sampling approximates this integral; it does not correct a misspecified conditional model or an unsuitable censoring assumption.

The Importance-Sampling Identity

Choose a proposal density \(q_i(b_i)\) from which it is straightforward to draw. It may depend on the subject’s observed data and on the current parameter values. Wherever \(h_i(b_i;\theta)\) is nonzero, the proposal must also be positive. Then multiplying and dividing the integrand by the proposal gives

L_i(\theta) = \int \frac{h_i(b_i;\theta)}{q_i(b_i)}q_i(b_i)\,db_i = E_{q_i}\left[\frac{h_i(B_i;\theta)}{q_i(B_i)}\right].

If \(B_i^{(1)},\ldots,B_i^{(M)}\) are independent draws from \(q_i\), define the importance weights \(w_i^{(m)}=h_i(B_i^{(m)};\theta)/q_i(B_i^{(m)})\). The estimator is

\widehat L_i(\theta)=\frac{1}{M}\sum_{m=1}^{M}w_i^{(m)}.

The weights are ratios of densities or likelihood densities evaluated at the same random effect. They are not generally probabilities and need not sum to one. If \(E_{q_i}[w_i^2]\) is finite, the estimator has variance \(\operatorname{Var}_{q_i}(w_i)/M\), and the central limit theorem gives the usual large-sample approximation for its sampling error. If this second moment is infinite, the estimator can remain unbiased under the support condition, but its variability may be so large that ordinary error estimates and practical convergence are unreliable.

Key condition: The proposal must cover the support of the integrand: if \(h_i(b_i;\theta)>0\), then \(q_i(b_i)>0\). For stable estimation, it is also important that the proposal not be too small in regions where \(h_i\) is appreciable. Support coverage prevents a structural omission; good coverage prevents a few extreme weights from dominating the estimate.

Worked Examples

Worked Example: A Three-Draw Integral Estimate

Consider a one-dimensional, hypothetical integrand \(h(b)=0.08\exp\{-\tfrac12[(b-1.2)/0.4]^2\}\). Its integral is \(0.08(0.4)\sqrt{2\pi}\approx0.08021\). Use the proposal \(q(b)\) for \(B\sim N(1.2,0.25)\), whose density is \(q(b)=0.79788456\exp\{-(b-1.2)^2/0.5\}\). Suppose three independent proposal draws, rounded for display, are \(0.8,1.2,1.6\).

At \(b=1.2\), \(h(1.2)=0.08000\) and \(q(1.2)=0.79788456\). Thus \(w(1.2)=0.08000/0.79788456\approx0.100265\).

At either \(b=0.8\) or \(b=1.6\), the distance from the mode is \(0.4\). Therefore \(h(b)=0.080e^{-0.5}\approx0.04852245\), while \(q(b)=0.79788456e^{-0.32}\approx0.579383\). Each off-center weight is \(0.04852245/0.579383\approx0.083748\).

The estimate is

\widehat L = \frac{0.083748+0.100265+0.083748}{3} \approx 0.08925.

This estimate exceeds the exact integral, about \(0.08021\), but that is possible with only three draws: unbiasedness is a repeated-sampling property, not a guarantee that every estimate is close. The proposal is centered at the integrand’s mode and is wider than the integrand, so it covers the important region. More draws would be needed to assess the estimate’s precision.

Worked Example: Building a Subject-Specific Estimate

For a hypothetical subject in a joint model, suppose the random effect is scalar and the observed longitudinal and event-time data produce the following integrand values. For illustration, a non-Gaussian proposal centered near the subject-specific posterior mode is used. The \(q(b)\) column gives its density evaluated at each draw.

Draw \(b^{(m)}\)	\(h(b^{(m)};\theta)\)	\(q(b^{(m)})\)	Weight \(h/q\)
0.8	0.010	0.20	0.050
1.0	0.040	0.40	0.100
1.2	0.080	0.40	0.200
1.4	0.040	0.20	0.200

For example, the draw at \(b=1.2\) contributes weight \(0.080/0.40=0.200\). Averaging all four weights gives

\widehat L_i(\theta) = \frac{0.050+0.100+0.200+0.200}{4} = 0.1375.

This is an estimate of the marginal likelihood contribution for this subject at the stated parameter values. The table is illustrative: to judge its quality, one would repeat the computation with more draws or independent batches and examine the variability. The proposal’s center should be chosen using the subject’s data, not automatically set to the population mean of the random effects.

Choosing and Checking a Proposal

The ideal proposal for estimating \(L_i\) is the normalized integrand \(h_i(b_i;\theta)/L_i(\theta)\), which is the subject-specific posterior density of the random effects. Using it would make every importance weight equal to \(L_i(\theta)\), but its normalizing constant is the integral being estimated. In practice, a proposal can approximate this posterior using a mode and local curvature, a Gaussian approximation, or information from a previous iteration. A proposal with somewhat heavier tails than the posterior approximation can protect against under-covering uncertainty.

The following checks are useful before trusting a subject-level estimate:

Support: verify that the proposal can generate every region where the integrand is positive.
Weight concentration: inspect the largest weights and the distribution of weights. A few weights much larger than the rest signal poor coverage or a rare high-contribution region.
Replication: compare estimates from independent batches or increasing numbers of draws. Large changes indicate Monte Carlo error or instability.
Numerics: evaluate likelihoods and weights on the log scale when products of densities are very small. A log weight is \(\log h_i(b_i;\theta)-\log q_i(b_i)\).

One summary is the effective sample size, calculated from normalized weights \(\tilde w_m=w_m/\sum_r w_r\) as \(ESS=1/\sum_m\tilde w_m^2\). It ranges from 1 to \(M\): equal weights give \(M\), while a single dominant weight gives an ESS close to 1. This is a diagnostic of weight concentration, not a substitute for checking repeated estimates or tail coverage.

Worked Example: Detecting Weight Concentration

Suppose four positive weights for one subject are \(0.02, 0.03, 0.05,\) and \(0.90\). Their sum is \(1.00\), so the normalized weights are the same values. The effective sample size is

ESS = \frac{1}{0.02^2+0.03^2+0.05^2+0.90^2} = \frac{1}{0.8138} \approx1.229.

Although four draws were taken, the weights behave roughly like just over one equally weighted draw. This is a warning that the estimate may be highly variable. Increasing the draw count can help, but if the proposal repeatedly misses an important part of the posterior, a better-centered or heavier-tailed proposal is needed as well.

Justification, Alternatives, and Limitations

The identity behind the estimator also explains its unbiasedness. Provided the proposal covers the integrand’s support and the expectation exists,

E_{q_i}\left[\frac{h_i(B_i;\theta)}{q_i(B_i)}\right] = \int \frac{h_i(b_i;\theta)}{q_i(b_i)}q_i(b_i)\,db_i = \int h_i(b_i;\theta)\,db_i = L_i(\theta).

Independence of the \(M\) draws then makes the average unbiased, and, when the weight variance is finite, gives \(\operatorname{Var}(\widehat L_i)=\operatorname{Var}_{q_i}(w_i)/M\). This derivation depends on using the actual proposal density in the denominator. An incorrect density changes the expectation and generally yields the wrong integral.

Importance sampling complements rather than replaces other integration strategies. The earlier material on numerical integration and Monte Carlo EM provides the context: quadrature evaluates a deterministic weighted grid, while Monte Carlo EM uses random-effects draws to approximate E-step expectations. Importance sampling instead reweights draws from a selected proposal to estimate an integral such as a subject’s marginal likelihood. It can be effective when the posterior is concentrated in a region that a data-adapted proposal captures well. Fixed quadrature can be more stable in low dimensions when a suitable grid is feasible; Monte Carlo methods are more adaptable to complex integrands but can be noisy.

A narrow proposal can miss posterior tails; a proposal centered at the population mean can miss a subject-specific mode. Either can produce very large or unstable weights. Conversely, a proposal that is excessively diffuse may spend most draws where the integrand is negligible. Increasing \(M\) reduces Monte Carlo variance at the usual \(1/M\) rate only when the weight variance is finite and the proposal is held fixed. It cannot repair a missing-support error or guarantee practical accuracy when weight variance is extreme. For an EM E-step, importance-sampling approximations to posterior expectations require normalized weights; those self-normalized estimates are generally not exactly unbiased at finite \(M\).

Common Mistakes

Using the prior as the proposal by default. The prior covers population variation, but informative subject data can concentrate the posterior elsewhere. Check whether the proposal represents the subject-specific posterior region.
Forgetting the proposal density. Averaging \(h_i(B_i)\) alone estimates an expectation under \(q_i\), not the required integral. Use \(h_i(B_i)/q_i(B_i)\).
Assuming a large draw count guarantees accuracy. A large number of draws can still yield an unstable estimate if a few weights dominate or important regions are poorly covered.
Confusing normalized and unnormalized weights. The marginal likelihood estimate averages unnormalized \(h_i/q_i\) weights. Normalized weights are used for posterior expectations and diagnostics such as ESS.

Key takeaway: Importance sampling estimates a marginal likelihood by averaging integrand-to-proposal ratios. Its validity requires support coverage; its practical usefulness depends on a proposal that tracks the subject-specific posterior closely enough to keep the weights stable.

Check Your Understanding

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

State the support condition required for an importance-sampling estimate of a subject’s marginal likelihood.
If every importance weight is multiplied by the same positive constant, what happens to the marginal likelihood estimate?
Four normalized weights are \(0.25,0.25,0.25,0.25\). Calculate the ESS and interpret it.
Why can a proposal centered at the population random-effects mean be poor for an individual subject?
What does a single very large weight suggest, and what are two useful responses?

Solutions. (1) Wherever \(h_i(b_i;\theta)>0\), the proposal must satisfy \(q_i(b_i)>0\). (2) The estimate is multiplied by the same constant because it is the average of the weights; it no longer estimates the original integral unless the constant is one. (3) \(ESS=1/[4(0.25^2)]=4\), the maximum possible for four draws, indicating equal weight contribution. (4) The subject’s observations may shift or sharpen the random-effects posterior away from the population mean, so that proposal may undersample regions of high integrand. (5) It suggests weight concentration, possibly from inadequate proposal coverage. Compare independent or larger batches, and revise the proposal to better cover the posterior, for example by recentering it or allowing heavier tails.