Worked Analysis of Two Groups With Unequal Covariances

Why the Covariance Assumption Matters

Comparing two multivariate means requires more than subtracting sample means. The uncertainty of that difference depends on the variability and correlation of outcomes in each group. If the groups have different covariance matrices, combining their covariance estimates as though they shared one common matrix can misstate that uncertainty. This worked analysis keeps the group covariances separate and uses the Welch-James approximation for calibration.

The method is approximate in finite samples. Its test statistic uses a covariance estimate appropriate to unequal group covariances, while its reference F distribution approximates the statistic’s null distribution. Accordingly, conclusions and confidence regions should be reported as approximate rather than exact.

Model and notation: Let $X_{ij}$ be the $p$-dimensional response vector for individual $j$ in group $i$, where $i=1,2$ and $j=1,\ldots,n_i$. Assume the observations are independent, with $X_{ij}\sim N_p(\mu_i,\Sigma_i)$ within group $i$. The two groups are independent of one another, and $\Sigma_1$ and $\Sigma_2$ need not be equal. Write $\bar X_i$ for the group sample mean and $S_i$ for the usual unbiased sample covariance matrix, with divisor $n_i-1$.

The inferential target is the difference $\delta=\mu_1-\mu_2$. To test $H_0:\delta=0$, the estimated difference is $d=\bar X_1-\bar X_2$. Because the group samples are independent, its estimated covariance is $V=S_1/n_1+S_2/n_2$. The separate terms retain each group’s own estimated covariance; no pooled covariance matrix is needed.

Worked Example: Specifying the Model and Estimating Its Components

Suppose two independent groups of patients are measured on two outcomes, systolic blood pressure and a symptom score. There are $n_1=12$ patients in group 1 and $n_2=10$ in group 2. Their observed summaries, in the outcome order stated, are

\bar X_1=\begin{pmatrix}130\\8\end{pmatrix},\qquad \bar X_2=\begin{pmatrix}127\\6\end{pmatrix},\qquad S_1=\begin{pmatrix}4&0\\0&9\end{pmatrix},\qquad S_2=\begin{pmatrix}16&0\\0&4\end{pmatrix}.

The model treats each patient’s two measurements as one response vector. It assumes independent patients within and between groups and a bivariate normal distribution within each group. The unequal estimated covariance matrices are retained separately. In particular, the estimated variance of systolic blood pressure is 4 in group 1 and 16 in group 2, whereas the estimated symptom-score variance is 9 and 4, respectively.

The estimated mean difference and its covariance estimate are

d=\bar X_1-\bar X_2=\begin{pmatrix}3\\2\end{pmatrix},\qquad V=\frac{S_1}{12}+\frac{S_2}{10} =\begin{pmatrix}\frac{4}{12}+\frac{16}{10}&0\\0&\frac{9}{12}+\frac{4}{10}\end{pmatrix} =\begin{pmatrix}\frac{29}{15}&0\\0&\frac{23}{20}\end{pmatrix}.

Thus, $V$ reflects the variance contributed by each group at its own sample size. The larger variance for systolic pressure in group 2 contributes $16/10$, not $16/12$ or a pooled estimate.

The Wald Statistic and Welch-James Calibration

Provided $V$ is nonsingular, the two-group Wald statistic is $Q=d'V^{-1}d$. Under the null hypothesis its reference distribution is not generally an exact chi-square distribution when the covariance matrices are estimated separately. The Welch-James procedure approximates its null distribution by a scaled F distribution, with degrees of freedom determined by the relative covariance contributions.

Welch-James two-group calibration: For $p$-dimensional outcomes, set $A_i=S_i/n_i$, so $V=A_1+A_2$, and define $C_i=V^{-1/2}A_iV^{-1/2}$. The Welch-James approximate denominator degrees of freedom are $$ \nu=\frac{p(p+1)}{\displaystyle\sum_{i=1}^2 \frac{\operatorname{tr}(C_i^2)+\{\operatorname{tr}(C_i)\}^2}{n_i-1}}. $$ For testing $H_0:\mu_1=\mu_2$, use $$ F_{\mathrm{WJ}}=\frac{\nu-p+1}{\nu p}\,Q \quad\text{and compare it with }F_{p,\,\nu-p+1}. $$ This is an approximate calibration, not an exact finite-sample result. It requires an invertible $V$ and positive denominator degrees of freedom.

The matrices $C_i$ describe how much each group contributes to the estimated covariance of the difference, after putting that covariance on a common scale. Their traces and squared traces, together with the group covariance degrees of freedom $n_i-1$, determine the approximation. When the groups contribute very differently, the effective degrees of freedom can be substantially smaller than the total sample size.

Worked Example: Testing the Mean Difference

Use the summaries above to test $H_0:\mu_1=\mu_2$ against the two-sided multivariate alternative. Since $V$ is diagonal, its inverse is

V^{-1}=\begin{pmatrix}\frac{15}{29}&0\\0&\frac{20}{23}\end{pmatrix}, \qquad Q=d'V^{-1}d =\frac{3^2}{29/15}+\frac{2^2}{23/20} =\frac{135}{29}+\frac{80}{23} \approx 8.133.

For the degrees of freedom, $C_1$ and $C_2$ are diagonal. Their diagonal entries are the corresponding entries of $A_i$ divided by those of $V$:

C_1=\begin{pmatrix}5/29&0\\0&15/23\end{pmatrix}, \qquad C_2=\begin{pmatrix}24/29&0\\0&8/23\end{pmatrix}.

For a diagonal matrix with entries $a,b$, the trace expression is $a^2+b^2+(a+b)^2=2(a^2+b^2+ab)$. Thus the two contributions to the denominator of $\nu$ are

\frac{2\{(5/29)^2+(15/23)^2+(5/29)(15/23)\}}{11} \approx 0.10318,

\frac{2\{(24/29)^2+(8/23)^2+(24/29)(8/23)\}}{9} \approx 0.24305.

Adding them gives approximately $0.34623$. With $p=2$, the numerator of $\nu$ is $2(3)=6$, so $\nu\approx 6/0.34623=17.33$. The approximate F statistic is therefore

F_{\mathrm{WJ}}=\frac{17.33-2+1}{17.33(2)}(8.133)\approx 3.83, \qquad F_{\mathrm{WJ}}\ \dot\sim\ F_{2,\,16.33}.

For an $F$ variable with 2 numerator degrees of freedom and $v$ denominator degrees of freedom, the upper-tail probability at $f$ is $\{v/(v+2f)\}^{v/2}$. Substituting $v=16.33$ and $f=3.83$ gives an approximate p-value of

\left(\frac{16.33}{16.33+2(3.83)}\right)^{16.33/2}\approx 0.043.

At the 5% level, this approximate test rejects equality of the two mean vectors. The result concerns the joint vector of outcomes; it does not by itself identify which individual outcome accounts for the evidence.

Confidence Interpretation and Two Useful Justifications

The same calibration gives an approximate joint confidence region for $\delta$. It is a region of candidate mean differences whose Wald distances from the observed difference do not exceed the critical value. For confidence level $1-\alpha$, its boundary is specified by

(\delta-d)'V^{-1}(\delta-d) \leq \frac{\nu p}{\nu-p+1}F_{p,\,\nu-p+1;\,1-\alpha}.

For the example, $F_{2,16.33;0.95}$ is approximately 3.62. The corresponding cutoff for the quadratic distance is $17.33(2)3.62/16.33\approx7.68$. The observed difference has distance $Q\approx8.13$ from zero, just above that cutoff. Equivalently, zero is outside this approximate 95% region. This interpretation remains joint: the region describes plausible pairs of outcome differences, not two separately calibrated 95% intervals.

Result: Covariance of the sample mean difference. If observations are independent across groups and within each group, with finite covariance matrices $\Sigma_i$, then $\operatorname{Cov}(\bar X_1-\bar X_2)=\Sigma_1/n_1+\Sigma_2/n_2$. If $S_i$ is the usual unbiased sample covariance, $V=S_1/n_1+S_2/n_2$ is an unbiased estimator of this covariance.

To verify the result, independence implies that the covariance of the difference is the sum of the two sample-mean covariances: the cross-covariance terms are zero. Within group $i$, independence gives

\operatorname{Cov}(\bar X_i) =\operatorname{Cov}\left(\frac{1}{n_i}\sum_{j=1}^{n_i}X_{ij}\right) =\frac{1}{n_i^2}\sum_{j=1}^{n_i}\operatorname{Cov}(X_{ij}) =\frac{\Sigma_i}{n_i}.

Adding these expressions gives the stated covariance. Since $E[S_i]=\Sigma_i$, linearity of expectation gives $E[V]=\Sigma_1/n_1+\Sigma_2/n_2$. This justifies using the two separate covariance estimates in the Wald statistic; it does not establish that its Welch-James F calibration is exact. $\square$

Result: Invariance to a change of outcome coordinates. For any nonsingular $p\times p$ matrix $B$, replacing every response vector by $BX_{ij}$ leaves the Wald statistic $Q$ unchanged.

Under this transformation the estimated difference becomes $Bd$, and its estimated covariance becomes $BVB'$. Its inverse is $B'^{-1}V^{-1}B^{-1}$. Consequently,

(Bd)'(BVB')^{-1}(Bd) =d'B'B'^{-1}V^{-1}B^{-1}Bd =d'V^{-1}d =Q.

Thus, changing units or applying another nonsingular linear re-expression of all outcomes does not change the Wald distance. This fact applies to $Q$; it does not make marginal confidence intervals or individual outcome interpretations invariant. $\square$

Assessing Power

Power is the probability of rejecting the null under a specified alternative. For planning, one needs an assumed mean difference and covariance structure, not just the observed sample difference. A useful large-sample approximation is that the Wald statistic under a fixed alternative has approximately a noncentral chi-square distribution with $p$ degrees of freedom and noncentrality parameter $\lambda=\delta'V_0^{-1}\delta$, where $V_0$ is the planned covariance of the estimated difference. This approximation ignores the additional variability from estimating the covariance matrices, so it is not the finite-sample Welch-James power calculation.

Worked Example: A Large-Sample Power Assessment

Suppose a planned two-outcome study has covariance of the estimated difference $V_0=I_2$, and the scientifically relevant difference is $\delta=(2,2)'$. The noncentrality parameter is

\lambda=\delta'V_0^{-1}\delta =(2,2)\begin{pmatrix}1&0\\0&1\end{pmatrix}\begin{pmatrix}2\\2\end{pmatrix} =2^2+2^2=8.

Using the large-sample 5% critical value $\chi^2_{2;0.95}\approx5.991$, approximate power is $P\{\chi^2_2(8)>5.991\}$. For two degrees of freedom, the noncentral chi-square upper tail has the Poisson-mixture representation

P\{\chi^2_2(\lambda)>x\} =\sum_{k=0}^{\infty}e^{-\lambda/2}\frac{(\lambda/2)^k}{k!} P\{\chi^2_{2+2k}>x\}.

Here the Poisson mean is $\lambda/2=4$, and the central chi-square tails in the sum can be calculated as $P\{\chi^2_{2+2k}>x\}=e^{-x/2}\sum_{j=0}^{k}(x/2)^j/j!$. At $x=5.991$, the Poisson weights for $k=0,\ldots,10$ are approximately $0.0183, 0.0733, 0.1465, 0.1954, 0.1954, 0.1563, 0.1042, 0.0595, 0.0298, 0.0132, 0.0053$. The corresponding central chi-square tails are approximately $0.0500, 0.1998, 0.4241, 0.6480, 0.8157, 0.9162, 0.9663, 0.9877, 0.9957, 0.9984, 0.9992$. Multiplying each pair and summing gives approximately $0.716$; the remaining terms contribute less than $0.001$. The large-sample power estimate is therefore about $0.72$.

This calculation is a planning illustration, not a claim that a finite-sample Welch-James test has exactly 72% power. For a modest sample, power should instead be evaluated using the noncentral version of the chosen approximate F calibration, or by simulation under the specified group distributions. The assumptions about the effect and covariance should also be varied to see how sensitive the planned power is.

Common Mistake and Practical Interpretation

Wrong reasoning: “The groups have different covariance matrices, but I can pool $S_1$ and $S_2$, since pooling gives one more stable covariance estimate.” Pooling imposes a common-covariance structure. If that structure is not justified, it changes the estimated uncertainty of the difference and no longer matches the unequal-covariance analysis.

Correction: Estimate the covariance of the difference as $S_1/n_1+S_2/n_2$, and use a calibration designed for unequal group covariances, such as the Welch-James approximation. Report the approximation and its degrees of freedom. Check that the covariance estimate is nonsingular, and remember that small samples, strong departures from multivariate normality, or influential observations can make the approximation less dependable.

The practical value of the procedure is that it separates two questions that are easy to conflate: how far apart the estimated mean vectors are, and how uncertain that difference is under the two groups’ distinct variability patterns. Reporting the estimated difference alongside the test and its confidence region keeps the direction and scale of the estimated effects visible, even when the joint test gives a single decision.

Check Your Understanding

Use the model and calculations in this tutorial to check your reasoning.

Why is the estimated covariance of the mean difference $S_1/n_1+S_2/n_2$, rather than a pooled covariance matrix?
What does the Welch-James denominator degrees of freedom represent, and why is its calibration approximate?
For the worked test, how is the F statistic obtained from $Q$, $p$, and $\nu$?
What does it mean for zero to lie outside the approximate joint confidence region?
Why does a power calculation require a specified alternative and covariance structure?