Behaviour of Welch's t-test with unequal group sizes

Question

It seems that in some cases Welch's t-test can give inflated p-values with imbalanced group sizes. For example, consider the following simulation of unequal group sizes sampled from the same Gaussian distribution.

set.seed(42)
sapply(2:10, function(n) {
  out <- replicate(100000, t.test(rnorm(100), rnorm(n))$p.value)
  mean(out < 0.01) / 0.01
})
# [1] 7.928 3.139 1.913 1.606 1.371 1.285 1.178 1.141 1.192

Is this conclusion correct? Are there general guidelines on group sizes for when a Welch's t-test is appropriate? Are there known ways to correct this apparent bias, or should we in general assume equal variance if groups are imbalanced in terms of size?

Update:

To clarify, I'm interested in establishing minimum sample sizes for a Welch t-test, below which we need to be cautious about the possibility of false positives (restricted to the situation where we are happy to assume equal variances).

I think this figure of a simple simulation of samples from the same distribution, with varying sample sizes more cleanly illustrates the issue. In all cases, we expect the mean number of tests with p-value < 0.05 to equal 0.05 if the test performs as expected.

I would conclude that we need to be cautious with a Welch t-test if N < ~5 for any group, particularly if N = 2.

Note, I've not investigated unequal variance, which as BruceET has pointed out can cause potentially high false positive rates when not accounted for.

dosim <- function(var.equal) {
  sapply(seq_len(nrow(g)), function(i) {
    pv <- replicate(
      reps,
      t.test(rnorm(g[i, "N1"]), rnorm(g[i, "N2"]), var.equal = var.equal)$p.value
    )
    mean(pv < 0.05)
  })
}

reps <- 10000
g <- expand.grid(list(N1 = 2:100, N2 = c(2, 3, 4, 5, 6, 10)))

g$Welch <- dosim(FALSE)
g$Student <- dosim(TRUE)
g <- tidyr::gather(g, "stat", "value", "Welch", "Student")
g$N2 <- factor(g$N2)

ggplot(g, aes(x = N1, y = value, col = N2)) +
  geom_point() +
  geom_line() +
  facet_grid(stat ~ .) +
  ylab("Mean(p < 0.05)")

Related:

• Always use Welch-t test (unequal variances t-test) instead of Student-t or Mann-Whitney test?
• Is there a minimum sample size required for the t-test to be valid?

Answer 1

This question is sufficiently broad that a comprehensive answer would require a simulation study (some key parts of which have undoubtedly been done), going far beyond our usual style of answers.

A Welch 2-sample t test can't be exactly as good as a pooled t test if we know the two populations have the same variance.

(1) A Welch test at "level 5%" with $n_1 = 5, n_2 = 100$ when variances are equal gave Type I error a little above 5% $(0.054 \pm 0.0014).$ I didn't explore samples sizes below 5.

set.seed(1234)
out = replicate(10^5, as.numeric(t.test(rnorm(5), rnorm(100))[2:3]))
mean(out[1,])
[1] 4.902672         # avg DF about 5
mean(out[2,] < .05)
[1] 0.05443          # avg P-val about 5%

In the histogram below of 100,000 simulated P-values, the maroon bar at left represents the true significance level (Type I error) of the test. Its width is 0.5 and its height is about 0.108, for a total area of about 0.054,

(2) Pooled test at level 5% with $n_1 = 5, n_2 = 100$ when variances are equal. Here, DF = $n_1 + n_2 - 2 = 103$ for all tests. The significance level is very near the intended 5%.

set.seed(1234)
out = replicate(10^5, as.numeric(t.test(rnorm(5), rnorm(100), var.eq=T)[2:3]))
mean(out[1,])
[1] 103             # DF exactly n1 + n2 - 2 = 103
mean(out[2,] < .05)
[1] 0.0493          # P-val exactly 5% from transparent theory

(3) However, if the smaller sample has variance 4 and the larger has variance 1, then the pooled tests with a nominal 5% significance level actually has significance level almost 30% (rich in false discoveries).

set.seed(1234)
out = replicate(10^5, as.numeric(t.test(rnorm(5,0,2), rnorm(100), var.eq=T)[2:3]))
mean(out[1,])
[1] 103
mean(out[2,] < .05)
[1] 0.2853

(4) The Welch test may not have exactly 5% Type I error. Even so, with average P-value near 5%, it is clearly preferable to the pooled test when variances are unequal, as in (3).

set.seed(1234)
out = replicate(10^5, as.numeric(t.test(rnorm(5,0,2), rnorm(100))[2:3]))
mean(out[1,])
[1] 4.209818
mean(out[2,] < .05)
[1] 0.05158

(5) First of two simulations focused on power. If variances are unequal and sizes are grossly unbalanced $(n_1 = 5, n_2 =100$), the power to detect that $\mu_1 = 4$ differs from $\mu_2 = 0$ is about 90%.

set.seed(1234)
out = replicate(10^5, as.numeric(t.test(rnorm(5,4,2), rnorm(100))[2:3]))
mean(out[1,])
[1] 4.209818
mean(out[2,] < .05)
[1] 0.90977

In the figure below, the green bar at left represents the power of the test. So now, a tall first bar is a good thing.

(6) If we balance the data (both samples of size 5), the power of the Welch test is reduced very little if any, compared to the previous simulation. Changing the means to be equal, leaving sample sizes at 5, and leaving variances unequal, a similar simulation (not shown) gave significance level 4.9%

set.seed(1234)
out = replicate(10^5, as.numeric(t.test(rnorm(5,4,2), rnorm(5))[2:3]))
mean(out[1,])
[1] 6.000759
mean(out[2,] < .05)
[1] 0.89743

Summary comments. None of these simulations changes the advice to prefer the Welch test to the pooled test, except for tiny sample sizes, In that case prior information whether variances are equal might help to decide between Welch test and a pooled test.

I have tried to probe in relevant directions with the few simulations above, but I certainly don't claim that they settle anything. If anyone has information from related simulations that would refine or extend my 'no new news' conclusions here, I'd be happy to see them.

Answer 2

The issue here does not appear to be due to imbalanced sample sizes per se, but rather, non-uniformity of the p-value seems to be occurring due to using very small sample sizes for one of the samples. Below I show high-resolution histograms showing simulated p-values for a comparison of one sample with $n_X = 100$ standard normal values and another sample with $n_Y = 2,...,10$ standard normal values. (Each comparison uses $K = 10^6$ simulations, and the histograms use bin widths $w = 0.01$ so you get high-resolution for the true distributions.)

As you can see, for very small samples there is a non-uniform distribution with a higher than expected probability of a small p-value. I am not entirely certain where this phenomena originates, but I have some strong suspicions. It is well-known that the standard sample variance estimator yields an estimated standard deviation that is biased for small samples (see e.g., here). For normal data the sample standard deviation is biased downward by a known "correction factor" that has quite a heavy effect for very small samples. If one were to underestimate the true standard deviation of the smaller sample in Welch's T-test, then the natural expectation would be that this would underestimate the likely different between sample means under the null hypothesis of no difference, which would bias the p-value downward. Since this is exactly what we see in the histograms, my suspicion is that this phenomena occurs due to the downward bias of sample standard deviations as an estimate of true standard deviations. I suspect that if you were to apply the standard small sample "corrections" to the sample standard deviations, in the test statistic for Welch's test, then this phenomena would be substantially ameliorated.

#Generate p-value simulations from Welch's test
#First sample is set to sample size of 100
#Second sample has sample sizes from 2, ..., 30
set.seed(97142903)
SIMS  <- 10^6;
PVALS <- matrix(NA, nrow = SIMS, ncol = 9);
for (n in 2:10)   {
for (i in 1:SIMS) {
    PVALS[i, n-1] <- t.test(rnorm(100), rnorm(n))$p.value; } }

#Generate data frame of p-value simulations
DATA <- data.frame(n = rep(2:10, each = SIMS),
                   p = as.vector(PVALS));

#Plot the KDEs of the p-value simulations
library(ggplot2);
THEME <- theme(plot.title    = element_text(hjust = 0.5, size = 14, face = 'bold'),
               plot.subtitle = element_text(hjust = 0.5, face = 'bold'));
FIGURE <- ggplot(aes(x = p), data = DATA) +
          geom_histogram(binwidth = 0.01, boundary = 0, fill = 'blue') +
          facet_wrap(~ n, ncol = 1) +
          THEME +
          ggtitle('Histograms of p-value simulations') +
          labs(subtitle = '(Sample of 100 against sample of stated size)') +
          xlab('p value') + ylab('Count');
FIGURE;