Worked Vector Delta-Method Calculations

From a Vector Limit to a New Estimator

Many estimators are not averages themselves but smooth functions of averages. A ratio, a covariance, or a vector of probability contrasts may all be obtained by applying a function to a vector of sample statistics. The vector delta method turns a central limit theorem for that underlying vector into an asymptotic distribution for the transformed estimator.

The calculation is systematic: identify the vector being estimated, write the transformation, evaluate its Jacobian at the population value, and multiply the covariance matrix on both sides by the Jacobian. The examples below emphasize the details that are easiest to get wrong: derivative signs, off-diagonal covariance terms, and the distinction between a limiting covariance and a finite-sample approximation.

Here the \(i,j\) entry of the Jacobian is the derivative of output coordinate \(i\) with respect to input coordinate \(j\). The matrix multiplication retains dependence among the input estimates. In particular, it is generally incorrect to calculate each transformed variance using only the corresponding input variances.

Worked Example: A Ratio of Two Sample Means

Let \((X_i,Y_i)\) be independent and identically distributed pairs with the following joint distribution, where each entry is a probability:

Consider estimating the ratio \(\mu_X/\mu_Y\) by \(\bar X/\bar Y\), defining the ratio when \(\bar Y>0\). The marginal means are \(\mu_X=0.50\) and \(\mu_Y=0.40\). Also, \(E[XY]=0.30\), so \(\operatorname{Var}(X)=0.50(0.50)=0.25\), \(\operatorname{Var}(Y)=0.40(0.60)=0.24\), and \(\operatorname{Cov}(X,Y)=0.30-(0.50)(0.40)=0.10\).

The input vector is \((\bar X,\bar Y)^{\mathsf T}\), whose limiting covariance for root-\(n\) scaling is \(\Sigma=\begin{pmatrix}0.25&0.10\\0.10&0.24\end{pmatrix}\). For \(g(a,b)=a/b\), the population ratio is \(g(0.50,0.40)=1.25\), and the Jacobian is the row vector

Thus the limiting variance is

Consequently, \(\sqrt n(\bar X/\bar Y-1.25)\xrightarrow{d}N(0,75/32)\). Since \(Y_i\) is Bernoulli with positive mean 0.40, the law of large numbers gives \(\bar Y\xrightarrow{p}0.40\); in particular, \(P(\bar Y>0)\to1\). The ratio is therefore defined with probability tending to one. For large \(n\), its approximate variance is \((75/32)/n\), not \(75/32\).

Worked Example: A Smooth Functional for Covariance

The covariance of \(X\) and \(Y\) can be written as \(E[XY]-E[X]E[Y]\), so it is a smooth function of three expectations. Use the same paired Bernoulli distribution as above and let \(W=XY\). The vector of empirical moments is \((\bar X,\bar Y,\bar W)\), with population value \((0.50,0.40,0.30)\). Define \(h(a,b,c)=c-ab\). Its gradient, viewed as a one-row Jacobian, is

To obtain the input covariance matrix, note that \(W\) is Bernoulli with mean 0.30. Further, \(\operatorname{Cov}(X,W)=E[X^2Y]-E[X]E[W]=0.30-0.15=0.15\), because \(X^2=X\), and \(\operatorname{Cov}(Y,W)=E[XY^2]-E[Y]E[W]=0.30-0.12=0.18\), because \(Y^2=Y\). Together with the previous calculations, this gives

The limiting variance of the estimated covariance is the scalar \(J_h\Sigma J_h^{\mathsf T}\). Writing out all diagonal and cross terms shows exactly how each covariance contributes:

This result can also be checked by deriving the first-order error directly. Write \(\mu_X=E[X]\), \(\mu_Y=E[Y]\), and \(\mu_W=E[W]\). The empirical covariance functional satisfies the exact identity

The three empirical means have root-\(n\) fluctuations by the multivariate central limit theorem. Each centered mean is therefore \(O_p(n^{-1/2})\), so their product in the last line is \(O_p(n^{-1})\), which becomes \(o_p(1)\) after multiplication by \(\sqrt n\). The leading term is the sample average of the centered variable \(\psi=W-\mu_W-\mu_Y(X-\mu_X)-\mu_X(Y-\mu_Y)\). The central limit theorem gives limiting variance \(\operatorname{Var}(\psi)=J_h\Sigma J_h^{\mathsf T}=0.05\). Thus the empirical covariance functional, centered at the population covariance \(0.10\), has root-\(n\) limiting variance \(0.05\). This direct expansion also explains why the product of the two empirical-mean errors does not contribute at first order.

The result applies to the plug-in covariance with divisor \(n\). The familiar sample covariance with divisor \(n-1\) differs by the factor \(n/(n-1)\); since this factor tends to one, that version has the same root-\(n\) limit.

Worked Example: Transforming a Multinomial Probability Vector

Suppose \(n\) independent categorical observations have three category probabilities \(p=(p_1,p_2,p_3)^{\mathsf T}=(0.20,0.30,0.50)^{\mathsf T}\), and let \(\hat p\) be the vector of observed proportions. Consider the two ratios \(g(p)=(p_1/p_3,\ p_2/p_3)^{\mathsf T}\). The population transformed value is \((0.40,0.60)^{\mathsf T}\).

For probabilities in the interior of the three-category simplex, the two-by-three Jacobian is

The multinomial covariance of the root-\(n\) limit for \(\hat p\) is \(\Sigma=\operatorname{diag}(p)-pp^{\mathsf T}\), which here equals

For clarity, calculate the entries of \(J_g\Sigma J_g^{\mathsf T}\). Its first diagonal entry is \(4(0.16)+(-0.8)^2(0.25)+2(2)(-0.8)(-0.10)=0.64+0.16+0.32=1.12\). The second is \(4(0.21)+(-1.2)^2(0.25)+2(2)(-1.2)(-0.15)=0.84+0.36+0.72=1.92\). The off-diagonal entry is \(2(2)(-0.06)+2(-1.2)(-0.10)+(-0.8)(2)(-0.15)+(-0.8)(-1.2)(0.25)\) \(=-0.24+0.24+0.24+0.24=0.48\). Therefore

The probability vector is constrained to sum to one, so its multinomial covariance matrix is singular. This causes no difficulty: the delta-method calculation uses the covariance on the actual probability-vector fluctuations, and the displayed transformed covariance is well-defined. The denominator \(p_3\) is positive, and \(\hat p_3\xrightarrow{p}0.50\), so the ratios are defined with probability tending to one.

Common Mistakes and Practical Checks

A delta-method calculation is short only after its inputs have been specified carefully. Before accepting a result, check the following:

• Dropping covariance terms. For a scalar transformation of a vector, expand \(J\Sigma J^{\mathsf T}\) or multiply the matrices explicitly. Off-diagonal terms can increase or decrease the variance.
• Changing a derivative sign. For a ratio \(a/b\), the derivative with respect to \(b\) is \(-a/b^2\), not \(a/b^2\). That negative sign affects the covariance contribution.
• Using only marginal variances. Numerator and denominator, or the empirical moments defining a covariance, are generally dependent. Their covariance is part of the calculation.
• Confusing covariance scales. The matrix in the root-\(n\) limiting distribution is the large-sample covariance for \(\sqrt n\) times the estimation error. The estimator itself has approximate covariance equal to that matrix divided by \(n\).
• Ignoring the domain of the transformation. A ratio requires a nonzero denominator at the population value and a denominator that is nonzero with probability tending to one. Smoothness at the population value is essential.

A useful final check is dimension: if the input has \(k\) coordinates and the output has \(r\), then \(J\) is \(r\)-by-\(k\), \(\Sigma\) is \(k\)-by-\(k\), and \(J\Sigma J^{\mathsf T}\) must be \(r\)-by-\(r\). For a scalar output it must reduce to a single variance.