The Multivariate Behrens-Fisher Problem: Comparing Mean Vectors under Unequal Covariances

In agricultural and biological research, experiments often involve measuring multiple correlated traits on each experimental unit. For instance, when comparing two crop varieties, researchers may analyze grain yield, protein content, and plant height simultaneously. When the data meet the assumptions of multivariate normality and equal covariance matrices (homoscedasticity) across both groups, Hotelling’s two-sample T² test is the standard method to compare mean vectors. However, covariance matrices are frequently unequal due to differences in genetic variability or environmental conditions. This scenario represents the Multivariate Behrens-Fisher problem. Using Hotelling’s pooled T² test under unequal covariance matrices, especially with unequal sample sizes, alters the Type I error rate and leads to incorrect statistical conclusions.

This guide explains the mathematical formulation of the Multivariate Behrens-Fisher problem, the consequences of ignoring unequal covariances, the Yao, Nel-Van der Merwe, and Johansen approximations, and a real-life agricultural example. We also provide Python and R implementations for analysis.

1. Mathematical Formulation of the Multivariate Behrens-Fisher Problem

Let X₁₁, X₁₂, …, X_1n1 be a random sample of p-dimensional vectors of size n₁ from a multivariate normal distribution with mean vector μ₁ and covariance matrix Σ₁. Let X₂₁, X₂₂, …, X_2n2 be an independent random sample of size n₂ from a multivariate normal distribution with mean vector μ₂ and covariance matrix Σ₂:

X_1i ~ N_p(μ₁, Σ₁),     i = 1, …, n₁

X_2j ~ N_p(μ₂, Σ₂),     j = 1, …, n₂

The parameters (μ₁, μ₂, Σ₁, Σ₂) are unknown. We want to test the null hypothesis that the population mean vectors are equal against the alternative hypothesis that they differ:

H₀: μ₁ – μ₂ = 0    vs.    H₁: μ₁ – μ₂ ≠ 0

When the covariance matrices are equal (Σ₁ = Σ₂ = Σ), we use Hotelling’s standard pooled T² statistic:

T² = [ (n₁ n₂) / (n₁ + n₂) ] (X̄₁ – X̄₂)’ S_p^-1 (X̄₁ – X̄₂)

Where X̄₁ and X̄₂ are the sample mean vectors, and S_p is the pooled sample covariance matrix:

S_p = [ (n₁ – 1)S₁ + (n₂ – 1)S₂ ] / (n₁ + n₂ – p – 1)

Here, S₁ and S₂ are the sample covariance matrices. Under the null hypothesis, [ (n₁ + n₂ – p – 1) / ((n₁ + n₂ – 2)p) ] T² follows an F-distribution with p and n₁ + n₂ – p – 1 degrees of freedom.

When Σ₁ ≠ Σ₂, the pooled covariance matrix S_p is a biased estimator of the covariance of the difference between the mean vectors. The components of the covariance matrices function as nuisance parameters. In classical frequentist statistics, no exact test exists where the test statistic follows a standard probability distribution independent of the nuisance parameters for finite samples. This difficulty is the Multivariate Behrens-Fisher problem.

2. Consequences of Ignoring the Multivariate Behrens-Fisher Problem

Ignoring unequal covariance matrices and using the pooled T² test leads to incorrect rates of Type I error (rejecting a true null hypothesis). The severity depends on the relation between sample sizes and covariance matrices:

• Equal Sample Sizes (n₁ = n₂): The pooled T² test is relatively robust. The actual Type I error rate remains close to the nominal level (e.g., 5%), though it slightly exceeds it as the differences between the covariance matrices increase. However, the test suffers a minor loss of power compared to tests designed for unequal covariance matrices.
• Larger Covariance in the Smaller Sample (n₁ < n₂ and Σ₁ > Σ₂ in terms of matrix ordering): The pooled covariance estimator understates the true covariance of the mean difference. The calculated T² value is artificially inflated, causing the actual Type I error rate to exceed the nominal level. A nominal 5% test can have an actual Type I error rate of 15% or higher, leading to excessive false positives.
• Larger Covariance in the Larger Sample (n₁ < n₂ and Σ₁ < Σ₂): The pooled covariance estimator overstates the true covariance. The test statistic is biased toward zero, making the test conservative. The actual Type I error rate falls below the nominal level (often under 2%), reducing the statistical power to detect real differences.

This problem is a direct multivariate extension of the univariate case. The consequences and general principles of sample size selection under heteroscedasticity are detailed in our guide on the Two-Sample Behrens-Fisher Problem and can be evaluated using a sample size calculator or detailed studies of sample size determination formulas.

3. Key Statistical Approximations

Because no exact frequentist solution exists, statisticians use approximations. The three most common frequentist methods are Yao’s test, Nel and Van der Merwe’s test, and Johansen’s F-approximation.

Yao’s Test

Yao (1965) proposed an approximation based on the statistic T_Y², defined as:

T_Y² = (X̄₁ – X̄₂)’ (S₁/n₁ + S₂/n₂)^-1 (X̄₁ – X̄₂)

Yao approximated the distribution of T_Y² under the null hypothesis by a Hotelling’s T² distribution with degrees of freedom f, estimated as:

1/f = ∑_i=1² (1 / (n_i – 1)) [ ( (X̄₁ – X̄₂)’ S̃^-1 (S_i/n_i) S̃^-1 (X̄₁ – X̄₂) ) / T_Y² ]²

Where S̃ = S₁/n₁ + S₂/n₂. The critical value is then obtained from the F-distribution: [ (f – p + 1) / (f p) ] T_Y² ~ F_{p, f–p+1}.

Nel and Van der Merwe’s Test

Nel and Van der Merwe (1986) used the same test statistic T_Y² but approximated the degrees of freedom f by equating the first two moments of the estimated covariance matrices to a Wishart distribution. The degrees of freedom f is calculated as:

f = [ tr(S̃²) + [tr(S̃)]² ] / ∑_i=1² (1 / (n_i – 1)) [ tr((S_i/n_i)²) + [tr(S_i/n_i)]² ]

Where tr(A) denotes the trace of matrix A. Under the null hypothesis, the statistic is approximated by:

[ (f – p + 1) / (f p) ] T_Y² ~ F_{p, f–p+1}

Johansen’s F-Approximation

Johansen (1980) proposed an alternative F-approximation for the test statistic. The test statistic is scaled by a constant c, such that:

T_Y² / c ~ F_{p, ν}

Where the scaling constant c and the degrees of freedom ν are calculated as:

c = p + 2 A – 6 A / (p (p + 2))

ν = p (p + 2) / (3 A)

The term A is computed from the sample covariance matrices and the sum matrix S̃:

A = (1 / 2) ∑_i=1² (1 / (n_i – 1)) [ tr( ( (S_i/n_i) S̃^-1 )² ) + [tr( (S_i/n_i) S̃^-1 )]² ]

Johansen’s test provides a robust correction when covariance matrices differ substantially, maintaining the Type I error rate close to the nominal level.

4. Real-Life Agricultural Example

To illustrate the calculations, we present a multivariate agricultural research example. An agronomist compares two wheat cultivars (Cultivar A and Cultivar B) grown under different conditions. Cultivar A is grown in a greenhouse using drip irrigation ($n_1$ = 10 plots), resulting in highly uniform conditions. Cultivar B is grown in open fields using flood irrigation ($n_2$ = 12 plots), resulting in higher variability. For each plot, two variables are measured: grain yield (tons per hectare) and grain protein content (%). The raw data are shown in Table 1. This type of analysis is detailed in the textbook Elements of Statistics for Agriculture and Forestry.

Plot No.	Cultivar A (Drip Irrigation)		Cultivar B (Flood Irrigation)
	Yield (tons/ha)	Protein (%)	Yield (tons/ha)	Protein (%)
1	5.2	12.1	4.5	13.5
2	5.5	12.3	4.8	14.2
3	5.1	11.9	4.2	12.8
4	5.4	12.2	5.1	15.1
5	5.3	12.0	4.6	13.8
6	5.2	12.1	4.3	13.1
7	5.6	12.4	4.9	14.5
8	5.3	12.2	4.5	13.4
9	5.1	11.9	4.1	12.6
10	5.3	12.1	4.8	14.0
11	–	–	4.4	13.3
12	–	–	4.6	13.7

Table 2 shows the summary statistics calculated from the raw values.

Statistic	Cultivar A (Group 1)	Cultivar B (Group 2)
Sample Size (n)	10	12
Mean Vector (X̄)	[5.3000, 12.1200]	[4.5667, 13.6667]
Variance-Covariance Matrix (S)	[ [0.0267, 0.0244], [0.0244, 0.0262] ]	[ [0.0879, 0.2097], [0.2097, 0.5097] ]

The variance of protein content for Cultivar B (0.5097) is approximately 19.4 times larger than that of Cultivar A (0.0262). The sample sizes are also unequal ($n_1$ = 10, $n_2$ = 12). This represents a severe case of heteroscedasticity requiring multivariate approximations.

Step 1: Calculate the difference vector and the sum matrix S̃

d = X̄₁ – X̄₂ = [5.3000 – 4.5667,   12.1200 – 13.6667]’ = [0.7333,   -1.5467]’

S̃ = S₁/10 + S₂/12 = [ [0.00267, 0.00244], [0.00244, 0.00262] ] + [ [0.00732, 0.01747], [0.01747, 0.04247] ] = [ [0.00999, 0.01992], [0.01992, 0.04510] ]

Step 2: Calculate the inverse S̃^-1 and the test statistic T_Y²

The determinant of S̃ is: |S̃| = (0.00999 × 0.04510) – (0.01992)² = 0.0004505 – 0.0003968 = 0.0000537.

S̃^-1 = (1 / 0.0000537) [ [0.04510, -0.01992], [-0.01992, 0.00999] ] = [ [839.85, -370.95], [-370.95, 186.03] ]

T_Y² = d‘ S̃^-1 d = [0.7333, -1.5467] [ [839.85, -370.95], [-370.95, 186.03] ] [0.7333, -1.5467]’ = 1736.80

Step 3: Calculate degrees of freedom and scaling factor (Johansen’s Method)

Using matrix algebra, we compute the term A for the sample. For $n_1$ = 10 and $n_2$ = 12:

A = 0.2825

The scaling constant c is:

c = 2 + 2(0.2825) – 6(0.2825) / 8 = 2 + 0.5650 – 0.2119 = 2.2354

The degrees of freedom ν is:

ν = 8 / (3 × 0.2825) = 8 / 0.8475 = 9.44 ≈ 9

The calculated F-statistic is:

F = T_Y² / c = 1736.80 / 2.2354 = 776.95

For p = 2 and ν = 9 degrees of freedom, the critical value of the F-distribution at significance level α = 0.05 is 4.256. Since F = 776.95 > 4.256, we reject the null hypothesis of equal mean vectors, indicating a significant difference in traits between Cultivar A and Cultivar B.

Nel and Van der Merwe’s Test Comparison:
Using Nel’s formula, the calculated degrees of freedom is f = 13.08. The corresponding F-statistic is:

F = [ (f – p + 1) / (f p) ] T_Y² = [ (13.08 – 2 + 1) / 26.16 ] × 1736.80 = 0.4618 × 1736.80 = 801.96

Comparing this to the critical F-value for 2 and 12 degrees of freedom (F_{0.05, 2, 12} = 3.885), the result is highly significant, confirming a difference between cultivar trait profiles.

5. R and Python Implementation

Here are the scripts to analyze the data from Example 1 using matrix calculations for the Nel & Van der Merwe and Johansen approximations.

R Implementation

# Raw data input
x1 <- matrix(c(
    5.2, 12.1, 5.5, 12.3, 5.1, 11.9, 5.4, 12.2, 5.3, 12.0,
    5.2, 12.1, 5.6, 12.4, 5.3, 12.2, 5.1, 11.9, 5.3, 12.1
), ncol = 2, byrow = TRUE)

x2 <- matrix(c(
    4.5, 13.5, 4.8, 14.2, 4.2, 12.8, 5.1, 15.1, 4.6, 13.8,
    4.3, 13.1, 4.9, 14.5, 4.5, 13.4, 4.1, 12.6, 4.8, 14.0,
    4.4, 13.3, 4.6, 13.7
), ncol = 2, byrow = TRUE)

n1 <- nrow(x1)
n2 <- nrow(x2)
p <- ncol(x1)

m1 <- colMeans(x1)
m2 <- colMeans(x2)

S1 <- cov(x1)
S2 <- cov(x2)

# Sum covariance matrix
S_tilde <- S1/n1 + S2/n2
S_tilde_inv <- solve(S_tilde)

diff <- m1 - m2
T2 <- as.numeric(t(diff) %*% S_tilde_inv %*% diff)

# Nel & Van der Merwe degrees of freedom
tr <- function(M) sum(diag(M))
f_nel <- (tr(S_tilde %*% S_tilde) + tr(S_tilde)^2) / (
    (1/(n1-1)) * (tr((S1/n1) %*% (S1/n1)) + tr(S1/n1)^2) +
    (1/(n2-1)) * (tr((S2/n2) %*% (S2/n2)) + tr(S2/n2)^2)
)

# Johansen parameters
A <- 0.5 * (
    (1/(n1-1)) * (tr((S1/n1 %*% S_tilde_inv) %*% (S1/n1 %*% S_tilde_inv)) + tr(S1/n1 %*% S_tilde_inv)^2) +
    (1/(n2-1)) * (tr((S2/n2 %*% S_tilde_inv) %*% (S2/n2 %*% S_tilde_inv)) + tr(S2/n2 %*% S_tilde_inv)^2)
)
c_joh <- p + 2*A - 6*A/(p*(p+2))
nu_joh <- p*(p+2)/(3*A)

cat("Hotelling T2:", T2, "\n")
cat("Nel df:", f_nel, "\n")
cat("Johansen c:", c_joh, "nu:", nu_joh, "F-stat:", T2/c_joh, "\n")

Python Implementation

import numpy as np

# Raw data input
x1 = np.array([
    [5.2, 12.1], [5.5, 12.3], [5.1, 11.9], [5.4, 12.2], [5.3, 12.0],
    [5.2, 12.1], [5.6, 12.4], [5.3, 12.2], [5.1, 11.9], [5.3, 12.1]
])

x2 = np.array([
    [4.5, 13.5], [4.8, 14.2], [4.2, 12.8], [5.1, 15.1], [4.6, 13.8],
    [4.3, 13.1], [4.9, 14.5], [4.5, 13.4], [4.1, 12.6], [4.8, 14.0],
    [4.4, 13.3], [4.6, 13.7]
])

n1 = len(x1)
n2 = len(x2)
p = x1.shape[1]

m1 = np.mean(x1, axis=0)
m2 = np.mean(x2, axis=0)

S1 = np.cov(x1, rowvar=False)
S2 = np.cov(x2, rowvar=False)

S_tilde = S1/n1 + S2/n2
S_tilde_inv = np.linalg.inv(S_tilde)

diff = m1 - m2
T2 = diff.dot(S_tilde_inv).dot(diff)

# Nel & Van der Merwe df
tr_St = np.trace(S_tilde)
tr_St2 = np.trace(S_tilde.dot(S_tilde))

denom = 0.0
for S_i, n_i in [(S1, n1), (S2, n2)]:
    S_n = S_i / n_i
    tr_Sn = np.trace(S_n)
    tr_Sn2 = np.trace(S_n.dot(S_n))
    denom += (1.0 / (n_i - 1)) * (tr_Sn2 + tr_Sn**2)

f_nel = (tr_St2 + tr_St**2) / denom

# Johansen parameters
A = 0.0
for S_i, n_i in [(S1, n1), (S2, n2)]:
    S_n = S_i / n_i
    mat = S_n.dot(S_tilde_inv)
    A += (1.0 / (n_i - 1)) * (np.trace(mat.dot(mat)) + np.trace(mat)**2)
A = A / 2.0

c_joh = p + 2*A - 6*A/(p*(p+2))
nu_joh = p*(p+2)/(3.0*A)

print(f"Hotelling T2: {T2:.4f}")
print(f"Nel df: {f_nel:.4f}")
print(f"Johansen c: {c_joh:.4f}, nu: {nu_joh:.4f}, F-stat: {T2/c_joh:.4f}")

6. Conclusion and Recommendations

The Multivariate Behrens-Fisher problem highlights the importance of checking statistical assumptions in multi-trait comparisons. Relying on standard pooled Hotelling's T² test under unequal covariance matrices alters Type I error rates, leading to false discoveries or missed findings. When covariance matrices differ, Yao's test, Nel & Van der Merwe's test, or Johansen's F-approximation should be used to protect the integrity of the analysis.