Importance Sampling with Correlated and High-Dimensional Effects

Why Correlation and Dimension Matter

In a nonlinear mixed-effects model, an individual’s random effects are often correlated. For example, clearance and volume effects may move together under the population covariance model, while the individual likelihood may concentrate on a narrow, tilted region of random-effects space. An importance proposal that treats coordinates independently can miss that geometry. Even when it covers the important region, small mismatches can accumulate across many dimensions and produce highly variable weights.

The goal here is to construct and assess multivariate and block proposals for integrating over random effects. The marginal likelihood integral and the importance-sampling identity are standard earlier material; the focus now is the geometry of the proposal and how its dimension and dependence affect the sampling weights.

Definition: For individual $i$, let $\eta_i\in\mathbb{R}^d$ denote the random-effects vector, with population distribution $p(\eta_i)$, and let $L_i(\eta_i)$ be the conditional likelihood of that individual’s observations given $\eta_i$. The marginal likelihood is $$ L_i=\int L_i(\eta_i)p(\eta_i)\,d\eta_i. $$ If $\eta_i^{(m)}$ are independently drawn from a proposal density $q_i$, the importance contribution is $W_i^{(m)}=L_i(\eta_i^{(m)})p(\eta_i^{(m)})/q_i(\eta_i^{(m)})$, and the Monte Carlo estimate is $\widehat L_i=M^{-1}\sum_{m=1}^M W_i^{(m)}$. The proposal must be positive wherever the integrand $L_i(\eta_i)p(\eta_i)$ is positive.

A common population model is $\eta_i\sim N(0,\Omega)$, where $\Omega$ is a positive-definite covariance matrix. If $\Omega=CC^\mathsf{T}$ is a Cholesky factorization and $z\sim N(0,I_d)$, then $\eta_i=Cz$ has the required correlated distribution. In these standard-normal coordinates, the integral becomes $$ L_i=\int L_i(Cz)\phi_d(z)\,dz, $$ where $\phi_d$ is the $d$-dimensional standard-normal density. A proposal can be defined on either the $\eta_i$ or $z$ scale, but its density and the target integrand must be expressed on the same scale, including the change-of-variables factor when needed.

Correlated and Block Proposals

A multivariate Gaussian proposal $q_i(\eta_i)=N(m_i,C_i)$ can represent both shifts in location and dependence among effects. Its covariance $C_i$ should be positive definite. A block proposal groups selected random effects and allows dependence within each group. Blocks are especially useful when a few effects are strongly coupled but a full covariance update is unstable or expensive. A fully independent proposal is the special case with diagonal covariance.

For a product proposal $q(z)=\prod_{b=1}^B q_b(z_b)$, the proposal assumes independence between blocks, while each $q_b$ can retain dependence among the coordinates inside its block. If the target integrand also factors by block, the total importance ratio is a product of block ratios. In general NLME likelihoods, the conditional likelihood may couple effects from different blocks, so this factorization of the proposal does not imply that the target factors.

Worked Example: Mapping Correlated Effects from Standard-Normal Draws

Consider four hypothetical random effects. Suppose the first pair has correlation $0.8$, the second pair has correlation $0.5$, and the pairs are independent of one another. Let $z=(z_1,z_2,z_3,z_4)^\mathsf{T}\sim N(0,I_4)$. The Cholesky transformation is $$ \eta_1=z_1,\qquad \eta_2=0.8z_1+0.6z_2,\qquad \eta_3=z_3,\qquad \eta_4=0.5z_3+\sqrt{0.75}\,z_4. $$

The coefficient $0.6=\sqrt{1-0.8^2}$ makes $\operatorname{Var}(\eta_2)=0.8^2+0.6^2=1$; likewise, $0.75=1-0.5^2$ keeps $\operatorname{Var}(\eta_4)=1$. The covariance of the first pair is $\operatorname{Cov}(\eta_1,\eta_2)=0.8$, and that of the second pair is $\operatorname{Cov}(\eta_3,\eta_4)=0.5$.

For the hypothetical draw $z=(0.5,-1,1,0)$, substitution gives $$ \eta_1=0.5,\qquad \eta_2=0.8(0.5)+0.6(-1)=0.4-0.6=-0.2,\qquad \eta_3=1,\qquad \eta_4=0.5(1)+\sqrt{0.75}(0)=0.5. $$ The resulting draw is $(0.5,-0.2,1,0.5)$. Sampling independent standard-normal coordinates and then applying this transformation generates correlated effects with the specified covariance.

How Weight Variability Compounds with Dimension

A useful way to study proposal mismatch is the relative second moment of the importance ratio. For a normalized target density $p$ and proposal $q$, define $W(x)=p(x)/q(x)$. Since $E_q[W]=1$, the quantity $R=E_q[W^2]$ is at least one, and $\operatorname{Var}_q(W)=R-1$. A large $R$ indicates variable weights; under standard large-sample conditions, the effective sample size fraction is approximately $1/R$. This approximation is a diagnostic, not a guarantee about a particular finite simulation.

The effect of dimension can be shown exactly when the target and proposal both factor across coordinates. Write $p(x)=\prod_{j=1}^d p_j(x_j)$ and $q(x)=\prod_{j=1}^d q_j(x_j)$. The ratio is $W(x)=\prod_{j=1}^d p_j(x_j)/q_j(x_j)$. Under the proposal, the coordinates are independent, so $$ E_q[W^2] =\prod_{j=1}^d E_{q_j}\left[\left(\frac{p_j(X_j)}{q_j(X_j)}\right)^2\right]. $$ If every coordinate has the same relative second moment $r$, then $R=r^d$. Thus, even a modest mismatch per coordinate can produce a large second moment as $d$ grows.

Worked Example: A Small Gaussian Mismatch in Many Dimensions

Let the normalized target be $p=N(0,I_d)$, and use the independent proposal $q=N(0,1.5I_d)$. In one dimension, for $p=N(0,1)$ and $q=N(0,\lambda)$, the ratio’s second moment is finite when $0<\lambda<2$, and direct integration gives $$ E_q[W^2]=\int\frac{p(x)^2}{q(x)}\,dx =\frac{1}{\sqrt{\lambda(2-\lambda)}}. $$ For $\lambda=1.5$, this is $1/\sqrt{1.5(0.5)}=1/\sqrt{0.75}\approx1.1547$. The proposal is broader than the target, but its one-coordinate weight variation is already greater than one.

Independence across dimensions makes the $d$-dimensional relative second moment $R=(1/\sqrt{0.75})^d$. Using the unrounded one-coordinate factor, $$ R_{10}=(1/\sqrt{0.75})^{10}\approx4.214,\qquad 1/R_{10}\approx0.2373, $$ and $$ R_{30}=(1/\sqrt{0.75})^{30}\approx74.83,\qquad 1/R_{30}\approx0.01336. $$ Under the large-sample ESS approximation, this corresponds to an effective sample size fraction of about $23.73\%$ at $d=10$, but only $1.336\%$ at $d=30$. These figures illustrate how a fixed, moderate per-coordinate mismatch can become consequential in higher dimensions.

Correlated Gaussian Mismatch

Matching each random effect’s marginal variance is not enough if the target and proposal have different dependence. For a normalized Gaussian target $p=N(0,\Sigma)$ and proposal $q=N(0,C)$, the relative second moment is $$ E_q[W^2]=\frac{|C|}{\sqrt{|\Sigma|\,|2C-\Sigma|}}, $$ provided $2C-\Sigma$ is positive definite. This condition ensures that the Gaussian integral for the second moment is finite. It is stronger than merely requiring both $\Sigma$ and $C$ to be positive definite.

To justify the formula, square the target density and divide by the proposal density. The exponential term in the resulting integrand is $\exp\{-\tfrac12x^\mathsf{T}(2\Sigma^{-1}-C^{-1})x\}$. It is integrable precisely when $A=2\Sigma^{-1}-C^{-1}$ is positive definite. Applying the multivariate Gaussian integral yields $$ E_q[W^2]=|\Sigma|^{-1}|C|^{1/2}|A|^{-1/2}. $$ Since $A=\Sigma^{-1}(2C-\Sigma)C^{-1}$, its determinant is $|2C-\Sigma|/(|\Sigma||C|)$; substituting gives the stated expression. If $C=\Sigma$, then $W=1$ everywhere and the formula gives $E_q[W^2]=1$, as it should.

Worked Example: The Cost of Ignoring Correlation

Suppose the hypothetical target is bivariate normal with unit marginal variances and correlation $0.8$: $$ \Sigma=\begin{pmatrix}1&0.8\\0.8&1\end{pmatrix}. $$ Use the independent proposal $C=I_2$, which matches both marginal variances but ignores the correlation. Here $|\Sigma|=1-0.8^2=0.36$, $|C|=1$, and $$ 2C-\Sigma=\begin{pmatrix}1&-0.8\\-0.8&1\end{pmatrix},\qquad |2C-\Sigma|=0.36. $$ This matrix is positive definite because its eigenvalues are $0.2$ and $1.8$, both positive. The second moment is therefore finite and equals $$ E_q[W^2]=\frac{1}{\sqrt{0.36(0.36)}}=\frac{1}{0.36}\approx2.778. $$ The corresponding asymptotic ESS fraction is approximately $1/2.778=0.3600$. A proposal with $C=\Sigma$ instead has constant weight and fraction one for this Gaussian target. In an NLME integral, the posterior-like integrand is not generally Gaussian, so matching this way cannot be exact in general; the calculation isolates the cost of this particular covariance mismatch.

Assessing Blocks and Diagnosing Failure

A practical proposal can use a full covariance matrix when the number of effects and available samples support a stable estimate. If only certain effects show strong dependence, blocks can represent those relationships while avoiding the burden of estimating every cross-block covariance. A coordinate-wise proposal is simpler, but it can be inefficient when important regions are tilted or curved relative to the coordinate axes. A single multivariate Gaussian can capture linear correlation, but it cannot in general reproduce a curved, skewed, or multimodal target.

Worked Example: Combining Block Ratios and Checking ESS

Suppose a hypothetical proposal samples two blocks independently. For four draws, the computed target-to-proposal density ratios for block A are $1,2,0.5,1.5$, and those for block B are $1,0.5,2,1$. Because the joint proposal is the product of the two block proposals, multiply the ratios within each draw to obtain joint weights $1,1,1,1.5$.

The importance-sampling effective sample size is $\mathrm{ESS}=(\sum_m W^{(m)})^2/\sum_m(W^{(m)})^2$. Here the sum of weights is $4.5$, and the sum of squared weights is $1+1+1+2.25=5.25$. Thus $$ \mathrm{ESS}=\frac{4.5^2}{5.25}=\frac{20.25}{5.25}\approx3.857. $$ The estimate is based on four draws, so the ESS is below four because the weights are unequal. This small example demonstrates how block ratios combine; it does not establish that a block proposal will always improve ESS. That depends on how well the proposal covers the integrand.

Useful diagnostics include comparing ESS across individuals and across proposal choices, inspecting the largest normalized weights, and checking whether a small number of draws dominate the estimate. Repeated very large weights can signal poor overlap, an inadequately represented dependence pattern, or tails that the proposal does not cover well. A covariance estimate may also be unreliable when there are too few effective draws, even if the raw Monte Carlo sample count is large. For correlated random effects, inspect whether the proposal captures the important joint directions, not only whether its marginal variances look plausible.

A defensive mixture can help guard against a proposal that undercovers parts of the target, but it does not eliminate dimension-related weight variability. Likewise, an adaptive proposal can improve its location and covariance using weighted draws, but adaptation cannot guarantee adequate coverage or finite weight variance. These methods should be assessed empirically using independent or suitably validated simulations, not judged by the proposal’s appearance alone.

Common Pitfall and Practical Limits

A common mistake is to conclude that a proposal is adequate because it matches each random effect’s marginal variance. The bivariate example shows why this reasoning fails: the target and proposal had identical marginal variances, yet ignoring correlation produced a relative second moment of about $2.778$. The correction is to assess the joint dependence and the resulting weights, using full or block covariance structures where justified.

There are limits to what covariance matching can accomplish. A high-dimensional full covariance matrix takes many parameters to estimate, and noisy estimates can make the proposal less stable rather than more effective. Block proposals trade some flexibility for more stable, manageable estimation, but can miss important cross-block dependence. A Gaussian proposal may also fit poorly when the integrand has separated modes or non-Gaussian tails. In those cases, weight diagnostics and comparison of plausible proposals are essential; no single covariance choice guarantees reliable Monte Carlo integration.

Key takeaway: Importance proposals for multivariate random effects should represent consequential dependence as well as marginal scale. Under factorization, second-moment mismatch compounds multiplicatively across dimensions; block or correlated proposals can address dependence, but their adequacy must be checked through weight behavior and coverage.

Check Your Understanding

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

If two independent target coordinates each have importance-ratio second moment $1.2$, what is the joint relative second moment?
For a bivariate target with correlation $0.6$ and unit marginal variances, what covariance should a Gaussian proposal use to match the target dependence exactly?
Why does a positive-definite proposal covariance not by itself guarantee a finite importance-weight second moment?
Four draws have weights $1,1,2,0$. Calculate their ESS.
When might a block proposal be preferable to a full covariance proposal, and what dependence might it fail to represent?

Solutions. 1. Independence makes the joint second moment the product, $1.2\times1.2=1.44$. 2. Use $C=\Sigma=\begin{pmatrix}1&0.6\\0.6&1\end{pmatrix}$; this matches the target covariance, including its correlation. 3. The second moment requires $2C-\Sigma$ to be positive definite in the Gaussian case; positive definiteness of $C$ alone does not ensure that condition. More generally, the proposal must cover the target and its weight ratio must have a finite second moment for that diagnostic. 4. The weights sum to $4$ and their squares sum to $1+1+4+0=6$, so $\mathrm{ESS}=4^2/6=16/6\approx2.667$. 5. A block proposal can be preferable when dependence is concentrated within a few groups and a full covariance estimate would be unstable or costly. It may fail to represent important dependence between different blocks.