Welch-James Approximation for Two Multivariate Means

Why Adjust the Wald Test?

When two groups have different covariance matrices, pooling their observations to estimate a common covariance can distort inference. The unequal-covariance Wald test avoids that pooling: it estimates the covariance of the difference in sample means separately from each group. Under the null hypothesis, its quadratic form has an asymptotic chi-square calibration. In finite samples, however, the covariance estimate is itself variable, and treating it as fixed can make the chi-square reference distribution inaccurate.

The Welch–James approximation addresses this finite-sample issue by using an approximate F reference distribution with adjusted degrees of freedom. The adjustment reflects how much each group contributes to the estimated covariance of the mean difference. It is not an exact finite-sample result for arbitrary unequal covariance matrices; its value is a more tailored calibration than the unadjusted chi-square test.

Statistic and Approximate Calibration

Suppose two independent groups have sample sizes $n_1,n_2$, sample mean vectors $\bar X_1,\bar X_2$, and unbiased sample covariance matrices $S_1,S_2$, each of dimension $p\times p$. The target is equality of the two population mean vectors. The estimated covariance of the difference in sample means is

\widehat V=\frac{S_1}{n_1}+\frac{S_2}{n_2}.

Write $d=\bar X_1-\bar X_2$. When $\widehat V$ is nonsingular, the unequal-covariance Wald quadratic form is $Q=d^\mathsf{T}\widehat V^{-1}d$. The Welch–James approximation scales $Q$ and assigns it an F distribution with numerator degrees of freedom $p$ and an estimated denominator degrees of freedom.

Definition: Let $A_i=S_i/n_i$, so that $\widehat V=A_1+A_2$, and define the symmetric matrices $H_i=\widehat V^{-1/2}A_i\widehat V^{-1/2}$, for $i=1,2$. The Welch–James estimated degrees of freedom are $$ \hat\nu=\frac{p(p+1)} {\displaystyle\sum_{i=1}^2\frac{(\operatorname{tr}H_i)^2+\operatorname{tr}(H_i^2)}{n_i-1}}. $$ The approximate F statistic is $$ F_{\mathrm{WJ}}=\frac{\hat\nu-p+1}{\hat\nu p}\,Q, $$ calibrated against $F_{p,\,\hat\nu-p+1}$. This calibration requires $\hat\nu>p-1$.

The matrices $H_i$ describe each group’s estimated contribution to the total covariance of the mean difference, on a scale where the two contributions sum to the identity matrix. Thus the degrees of freedom depend not just on sample sizes, but also on how the two groups’ covariance contributions are distributed across directions.

Why This Degrees-of-Freedom Formula Appears

The adjustment can be derived by approximating the random estimated covariance matrix with a scaled Wishart matrix. For the derivation, suppose the two groups are independent multivariate normal samples with population covariance matrices $\Sigma_i$. The sample covariance matrix satisfies $(n_i-1)S_i\sim W_p(\Sigma_i,n_i-1)$, where $W_p$ denotes the $p$-dimensional Wishart distribution. Let $V=\Sigma_1/n_1+\Sigma_2/n_2$ be the true covariance of the mean difference, and define the population contribution matrices $C_i=V^{-1/2}(\Sigma_i/n_i)V^{-1/2}$. They sum to the identity matrix.

Result: A moment-matched Wishart approximation to $V^{-1/2}\widehat V V^{-1/2}$ has degrees of freedom $$ \nu=\frac{p(p+1)} {\displaystyle\sum_{i=1}^2\frac{(\operatorname{tr}C_i)^2+\operatorname{tr}(C_i^2)}{n_i-1}}. $$ Replacing the unknown population contributions $C_i$ by the corresponding sample contributions $H_i$ gives $\hat\nu$.

Here is the covariance calculation behind the result. For a Wishart sample covariance, the covariance between entries $a,b$ and $c,d$ of the corresponding sample covariance contribution is proportional to $C_{i,ac}C_{i,bd}+C_{i,ad}C_{i,bc}$, with multiplier $1/(n_i-1)$. Summing the variances of all entries of the whitened estimated covariance gives

\sum_{a=1}^p\sum_{b=1}^p \operatorname{Var}\!\left\{\big(V^{-1/2}\widehat V V^{-1/2}\big)_{ab}\right\} = \sum_{i=1}^2\frac{(\operatorname{tr}C_i)^2+\operatorname{tr}(C_i^2)}{n_i-1}.

For a Wishart matrix divided by its degrees of freedom and having mean $I_p$, each entry covariance is $(\delta_{ac}\delta_{bd}+\delta_{ad}\delta_{bc})/\nu$, where $\delta_{ab}$ is 1 when $a=b$ and 0 otherwise. Summing the entry variances yields $p(p+1)/\nu$. Equating this total variability with the preceding expression gives the stated $\nu$. This matches aggregate covariance variability; it does not make the actual covariance estimate exactly Wishart when the two population covariance matrices differ.

Under the matching model, the classical Hotelling–Wishart relationship gives a scaled F distribution for the quadratic form. Specifically, if the whitened covariance estimate were exactly Wishart divided by $\nu$, then

Q\ \sim\ \frac{p\nu}{\nu-p+1}F_{p,\,\nu-p+1}.

Rearranging this relationship gives $F_{\mathrm{WJ}}=(\nu-p+1)Q/(\nu p)$, which motivates the stated test statistic. In practice $\nu$ is unknown, so the sample-based $\hat\nu$ is used. The argument explains the calibration; because the covariance estimate is not generally exactly Wishart with a single scale, the resulting F test remains approximate.

Proof of the degrees-of-freedom adjustment. Whiten each group’s covariance contribution by $V^{-1/2}$. Its expectation is then $C_i$, and the contributions sum to $I_p$. The Wishart covariance formula gives the total entrywise variance for group $i$ as $\{(\operatorname{tr}C_i)^2+\operatorname{tr}(C_i^2)\}/(n_i-1)$; independence means the two groups’ variances add. A mean-$I_p$ Wishart matrix divided by its degrees of freedom has total entrywise variance $p(p+1)/\nu$. Equating the two totals and solving for $\nu$ gives the formula. Substitution of the sample contributions gives $\hat\nu$, and the Hotelling–Wishart relationship gives the F scaling. $\square$

Worked Applications

Worked Example: Reduction to Welch’s Unequal-Variance Test

Take $p=1$, $n_1=10$, $n_2=20$, and suppose the estimated variance contributions to the mean difference are $A_1=0.4$ and $A_2=0.1$. Then $\widehat V=0.5$, and the one-dimensional contribution matrices are $H_1=0.4/0.5=0.8$ and $H_2=0.1/0.5=0.2$.

For $p=1$, $(\operatorname{tr}H_i)^2+\operatorname{tr}(H_i^2)=2H_i^2$. Therefore

\hat\nu =\frac{2}{2(0.8^2/9+0.2^2/19)} =\frac{1}{0.64/9+0.04/19} =\frac{4275}{313} \approx 13.66.

This is the Welch–Satterthwaite degrees-of-freedom formula for a two-sample unequal-variance comparison. For example, if the observed mean difference is $1.2$, then $Q=1.2^2/0.5=2.88$. With $p=1$, the scaling factor $(\hat\nu-p+1)/(\hat\nu p)$ equals 1, so $F_{\mathrm{WJ}}=2.88$, calibrated against $F_{1,13.66}$. Equivalently, this is the square of the unequal-variance t statistic.

Worked Example: A Two-Dimensional Heteroscedastic Comparison

Suppose $p=2$, $n_1=10$, $n_2=12$, and the estimated covariance contributions are

A_1= \begin{pmatrix}0.4&0\\0&0.2\end{pmatrix}, \qquad A_2= \begin{pmatrix}0.1&0\\0&0.3\end{pmatrix}.

Thus $\widehat V=0.5I_2$, so $H_1=\operatorname{diag}(0.8,0.4)$ and $H_2=\operatorname{diag}(0.2,0.6)$. For group 1, $\operatorname{tr}H_1=1.2$, and $\operatorname{tr}(H_1^2)=0.8$; its contribution to the denominator is $(1.2^2+0.8)/9=2.24/9$. For group 2, $\operatorname{tr}H_2=0.8$, and $\operatorname{tr}(H_2^2)=0.4$; its contribution is $(0.8^2+0.4)/11=1.04/11$. Consequently,

\hat\nu =\frac{2(2+1)}{2.24/9+1.04/11} =\frac{6}{34/99} =\frac{297}{17} \approx 17.47.

If the observed mean difference is $d=(1,0.5)^\mathsf{T}$, then

Q=d^\mathsf{T}(0.5I_2)^{-1}d =2(1^2+0.5^2)=2.5.

The adjusted statistic is therefore

F_{\mathrm{WJ}} =\frac{17.47-2+1}{17.47(2)}(2.5) \approx 1.178,

with approximate degrees of freedom $2$ and $17.47-2+1\approx16.47$. A test decision would compare this value with an upper-tail critical value from that F distribution, or report the corresponding upper-tail probability.

Worked Example: Checking the F Scaling

Suppose a two-dimensional analysis has $Q=6$ and its sample-based adjustment gives $\hat\nu=15.34$. The numerator degrees of freedom are $p=2$, and the denominator degrees of freedom are $15.34-2+1=14.34$. The scaling calculation is

F_{\mathrm{WJ}} =\frac{15.34-2+1}{15.34\cdot2}\,6 =\frac{14.34\cdot6}{30.68} =\frac{86.04}{30.68} \approx 2.804.

Thus the appropriate reference distribution is $F_{2,14.34}$, not a chi-square distribution and not an F distribution with denominator degrees of freedom equal to $\hat\nu$. The $-p+1$ adjustment is part of the conversion from the Hotelling-type quadratic form.

Assumptions and Practical Limits

The finite-sample moment calculation above assumes independent multivariate normal samples within the two groups. The approximation is most credible when the group sample sizes are adequate, the covariance estimates are stable, and the estimated degrees of freedom are comfortably above $p-1$. Both sample covariance matrices must be estimable, and $\widehat V$ must be nonsingular for $Q$ to be defined in the ordinary way.

For a fixed dimension $p$, a large-sample justification is also available under weaker distributional assumptions: observations are independent within and between groups, each group has a finite fourth moment, the group sample sizes grow with neither group becoming negligible, and the limiting covariance of the mean difference is positive definite. Under these conditions, the sample covariance estimates are consistent and the Wald quadratic form has its asymptotic chi-square calibration under the null. As sample sizes grow, the degrees-of-freedom adjustment grows as well, and the scaled F reference approaches that chi-square calibration. This does not guarantee good small-sample performance for strongly nonnormal data or unstable covariance estimates.

A practical limitation is that the approximation compresses the variability of a matrix estimate into a single degrees-of-freedom number. When the group contributions differ greatly across directions, that one-number summary may not capture every feature of the true finite-sample distribution. If $\hat\nu\le p-1$, the stated F calibration is unavailable; this is a warning that the approximation is not usable as written, rather than a reason to replace the value with an arbitrary positive number.

Common Pitfall

Using the wrong denominator degrees of freedom. The denominator is $\hat\nu-p+1$, not $\hat\nu$. The numerator remains $p$.
Pooling the covariance matrices. The test uses $S_1/n_1+S_2/n_2$, not a pooled covariance estimate that imposes equality of the group covariances.
Calling the F calibration exact. The degrees-of-freedom construction matches aggregate covariance variability to a Wishart model. With unequal covariance matrices, that match does not generally reproduce the full finite-sample distribution.
Dropping the trace terms. The formula contains both $(\operatorname{tr}H_i)^2$ and $\operatorname{tr}(H_i^2)$. Omitting either changes the adjustment and fails, in general, to recover Welch’s univariate degrees of freedom.

Key takeaway: The Welch–James test starts with the unequal-covariance Wald quadratic form and uses a moment-based degrees-of-freedom adjustment to obtain an approximate F calibration. The adjustment reflects both sample sizes and the groups’ relative covariance contributions; it improves calibration without making the test exact.

Check Your Understanding

Use the definitions and derivation above to answer these questions.

Why does the Welch–James statistic use $S_1/n_1+S_2/n_2$ rather than a pooled covariance estimate?
For a problem with $p=3$ and $\hat\nu=20$, what are the numerator and denominator degrees of freedom for the approximate F reference distribution?
Explain why matching total entrywise variance does not make the unequal-covariance estimate exactly Wishart.
Show how the Welch–James degrees-of-freedom formula reduces to the Welch–Satterthwaite formula when $p=1$.
What does it indicate if the estimated value satisfies $\hat\nu\le p-1$?