Is the t-test miscalibrated in R?

Question

When the null hypothesis is true, the p-value of a test should have the standard uniform distribution. Here is what I get with t.test(...) in R using two Gaussian samples of size 5.

set.seed(123)
p.val <- replicate(n=100000, t.test(rnorm(n=5), rnorm(n=5))$p.value)
hist(p.val, breaks=50)

You can see that there is a deficit of low p-values. Below is what I get with somewhat bigger samples of size 10.

set.seed(123)
p.val <- replicate(n=100000, t.test(rnorm(n=10), rnorm(n=10))$p.value)
hist(p.val, breaks=50)

The deficit of low p-values is gone. So what happens in the first example? Is there something wrong with t.test(...) in R for small-ish sample sizes?

---

> sessionInfo()
R version 4.2.1 (2022-06-23)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Catalina 10.15.7

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib

locale:
[1] en_CA.UTF-8/en_CA.UTF-8/en_CA.UTF-8/C/en_CA.UTF-8/en_CA.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

loaded via a namespace (and not attached):
[1] compiler_4.2.1

Answer 1

t.test performs Welch's t-test if the argument var.equal is not explicitly set to TRUE. The distribution of the test statistic (under the null hypothesis) in Welch's t-test is only approximated by a t-distribution and this approximation gets better with increasing sample sizes. Therefore, the result of your simulation is not particularly surprising.

---

Addendum
The test statistics in Welch's t-test and Student's t-test coincide if the two sample sizes $n_1$ and $n_2$ (of group $1$ and $2$, respectively) are equal. Hence, the discrepancy in the (simulated) p-value distributions of the two tests (note that the one of Student's t-test is uniform on $\left[0,1\right]$) under the null hypothesis is due to the discrepancy between the estimated degrees of freedom $\nu$ in Welch's t-test and the degrees of freedom $\tilde\nu=2\left(n-1\right)$ in Student's t-test, where $n=n_1=n_2$.

It is easy to see that, if $n_1=n_2$, the estimated degrees of freedom are given by $$ \nu = \frac{\left(n-1\right)\left(s_1^2 + s_2^2\right)^2}{s_1^4+s_2^4} = \frac{\left(n-1\right)\left(s_1^4 + s_2^4 + 2s_1^2s_2^2 \right)}{s_1^4+s_2^4}, $$ where $s_1$ and $s_2$ are the Bessel-corrected sample standard deviations.

By the AM–GM inequality and the non-negativity of sample standard deviations, $2s_1^2s_2^2 \leq s_1^4 + s_2^4$ (with equality only if $s_1 = s_2$) and $2s_1^2s_2^2 \geq 0$, therefore $n-1 \leq \nu \leq 2\left(n-1\right)=\tilde\nu$. This shows that the estimated degrees of freedom can only underestimate (or$-$but almost never$-$coincide with) the true degrees of freedom in the given situation, which leads to the conservative p-values seen in your simulation.
This behavior is nicely illustrated in Thomas Lumley's answer.

Since $s_1$ will tend to be closer to $s_2$ with increasing $n$, we can also see that $\nu$ will tend to be closer to $\tilde\nu$ as $n$ increases. Additionally, for a fixed difference $\nu - \tilde\nu$ in degrees of freedom of two t-distributions, their PDFs become increasingly similar with increasing $\nu$, and $\nu$ increases with $n$ in our case. This explains the improvement of the approximation and hence the p-value distribution with increasing group/total sample size.

Answer 2

Following up on @statmerkur's correct answer: first, here's what you get with var.equal=TRUE

which is well calibrated. Second, here's the distribution of estimated degrees of freedom for the Welch t-test

As you can see, the estimated degrees of freedom are typically near 8, but occasionally quite a bit smaller. When the df is smaller, the p-value will be conservative.

Finally, here's the p-value distribution separately for estimated df>7 and $\leq 7$