Distribution of \(t\) statistics (one sample \(t\) test)

The distribution of \(t\) statistics under an alternative hypothesis with true effect size \(\delta\) and sample size \(n\) is \[ t \sim \mbox{noncentral}~t_{n-1}(\delta\sqrt{n}). \]

n = 50
delta = .3
sigma = 1

# simulation
tstats = replicate(100000, t.test(rnorm(n, sigma*delta, sigma))$statistic)

## theoretical
tt = seq( min(tstats), max(tstats), len = 250 )
yy = dt(tt, df = n - 1, ncp = delta * sqrt(n) )

## plot
## simulation
hist(tstats, freq = FALSE, col="gray", breaks = 20)
## theoretical
lines(tt, yy, lwd=2,col = "red")

Distribution of \(F\) statistics

If we want the two-sided \(p\) values, it is easier to use \(F = t^2\), which has an noncentral \(F_{1,n-1}(\delta^2n )\) distribution. We denote this density function \(p_f\).

## theoretical
FF = seq( 0, max(tstats^2), len = 250 )
yy = df(FF, df1 = 1, df2 = n - 1, ncp = delta^2 * n)

## plot
## simulation
hist(tstats^2, freq = FALSE, col="gray", breaks = 20)
## theoretical
lines(FF, yy, lwd=2,col = "red")

Calculation of \(p\) values

The two-tailed \(p\) value is computed from the \(F\) statistic by way of the CDF of the central \(F\) distribution (central, because we assume the null hypothesis). If we denote this CDF as \(P_0\) (in R, it is simply pf()), then \[ p = 1 - P_0(F) \] where \(F\) is the observed \(F = t^2\) statistic.

pvals = 1 - pf(tstats^2, df1 = 1, df2 = n - 1)

## plot
## simulation
hist(pvals, freq = FALSE, col="gray", breaks = 20 )

Distribution of \(p\) values

In order to obtain the theoretical distribution of \(p\) values, we need to perform a transformation (see, for instance, the details here: https://www.stat.wisc.edu/courses/st309-larget/jacobian.pdf)

Our inverse transformation for computing the \(F\) statistic from the \(p\) value is \[ F = P_0^{-1}(1 - p) = h(p) \] where \(P_0^{-1}\) is the inverse function of the CDF. Call this function \(h(p)\). The inverse function of the CDF is the quantile function of the \(F\) distribution (in R, it is qf()). If we denote the density function of the \(p\) value as \(f\), then \(f\) is \[ f(p) = p_f(h(p)) |h'(p)| \] which is \[ f(p) = \frac{p_f(F)}{p_0(F)} \] where \(p_0\) is the density function of the central \(F\) distribution.

## f(p)
## theoretical
p.val.dens = Vectorize(function(p, n, delta){
  F = qf(1 - p, df1 = 1, df2 = n - 1)
  df(F, df1 = 1, df2 = n - 1, ncp = delta^2 * n) / df(F, df1 = 1, df2 = n - 1)
},"p")

pp = seq(.001, .999, len=250)
yy = p.val.dens(pp, n, delta)

## plot
## simulation
hist(pvals, freq = FALSE, col="gray", breaks = 20 )
## theoretical
lines(pp, yy, lwd=2,col = "red")

If you look closely, you’ll notice that the density function of the p value is just the likelihood ratio, which is a neat fact (and makes sense, because the distribution of \(p\) values under the null is uniform, and the likelihood ratio is invariant under monotone transformations such as the transformation from \(F\) to \(p\)).