Why doesn't the TOST equivalence testing procedure use non-central $t$ distribution to determine the $p$ value?

I'm looking at Lakens' (2017) primer and he tests the hypothesis that there is a difference between two groups $X_1$and $X_2$of magnitude $\Delta=E[X_1]-E[X_2] \neq0$ by subtracting $\Delta$ from the difference between the sample means $d=M_1-M_2$ and then he computes the $t$ value and $p$ value based on this difference score with a procedure analogous to Welch t-test (Eqs 3 and 4 on p. 357).

This looks wrong to me. IMO to determine the p value, one should be using non-central $t$ distribution with non-centrality parameter based on $\Delta$ and t value based on $d$. This is because under the hypothesis, $d$ and hence the test statistic $t=d/h$ (where we would substitute $h=\sqrt{s_1^2/n_1+s_2^2/n_2}$ following Welch) does not have a symmetric distribution because $E[X_1]-E[X_2]=\Delta\neq 0$. Subtracting $\Delta$ from $d$ does merely shift the distribution, but it can't make it symmetric, since $\Delta$ is not a random variable.

Using the zero centered T distribution here would test the hypothesis assuming that X1 - X2 -3 is symmetric about zero.

(X1-X2-3 corresponds to $X_1-X_2-\Delta$). I don't see how it is possible to obtain such symmetric test statistic distribution unless $E[X_1]=E[X_2]$ which under the hypothesis is not the case ($E[X_1]-E[X_2]=\Delta\neq 0$).

So, which is the correct procedure to conduct equivalence test?

Literature

Lakens, D. (2017). Equivalence tests: a practical primer for t tests, correlations, and meta-analyses. Social Psychological and Personality Science, 8(4), 355-362.

You are correct, in the sense that the uniformly most powerful (UMP) test for equivalence for the one-sample, paired, and two-sample $$T$$-test does involve a non-central $$t$$-distribution.

The two one-sided tests (TOST) procedure with both rejection regions with size $$\alpha$$, or equivalently the $$1 - 2 \alpha$$ confidence-interval inclusion procedure, has size $$\alpha$$ when the null hypothesis is true. However, it is not as powerful as the UMP test for small sample sizes. One can show that, as the sample size grows, the TOST procedure is asymptotically equivalent to the UMP test. But a great deal of power can be lost for small sample sizes.

This code simulates from a standard Normal distribution and performs the test for equivalence in the one-sample setting using both the TOST test (via Daniel Lakens TOSTER package) and the UMP test, where the P-value for the UMP test is given by $$P(|T| < |t_{\text{obs}}|; \texttt{df} = n - 1, \texttt{ncp} = \sqrt{n}\cdot\delta/\sigma)$$ where $$T$$ is a non-central $$t$$-random variable and $$t_{\text{obs}} = \frac{\bar{x}}{s/\sqrt{n}}$$.

n <- 10
alpha <- 0.05

theta <- 0.5 # delta / sigma

ncp <- theta*sqrt(n) # delta/sigma*sqrt(n)

# Simulate 10000 samples, and compute power.
S <- 10000

# theta.true <- theta # Null is true
theta.true <- 0 # Alternative is true

Pvals.tost <- rep(NA, length(S))
Pvals.ump  <- rep(NA, length(S))

for (s in 1:S){
Xs <- rnorm(n, mean = theta.true)

xbar <- mean(Xs)

tobs <- xbar/(sd(Xs)/sqrt(n))

Pvals.ump[s] <- pt(abs(tobs), df = n-1, ncp = ncp) - pt(-abs(tobs), df = n-1, ncp = ncp)

out <- TOSTER::TOSTone(m = xbar,
mu = 0,
sd = 1,
n = n,
low_eqbound_d = -theta,
high_eqbound_d = theta,
alpha = 0.05,
plot = FALSE,
verbose = FALSE)

Pvals.tost[s] <- max(out$$TOST_p1, out$$TOST_p2)
}

mean(Pvals.tost <= 0.05)
mean(Pvals.ump <= 0.05)


With $$n = 10$$, the TOST procedure has 0 power to detect the discrepancy from the null, while the UMP test has power ~ 0.17 to detect it.

Reference: Stefan Wellek, Testing Statistical Hypotheses of Equivalence and Noninferiority, pages 64 and 92.

What you are proposing does not appear to be any different from what actually occurs in this test. Remember that in a classical hypothesis test, the p-value is calculated using the null distribution of the test statistic. This calculation assumes that the null hypothesis is true. You suggest that the calculation of the p-value should use the non-central T-distribution with non-centrality parameter $$\Delta$$. However, the null hypothesis for this test is that $$\Delta = 0$$, so under that condition the non-central T-distribution simplifies down to the standard (centralised) T-distribution. (Note that the noncentral T-distribution with a non-zero noncentrality parameter would arise when you compute the power function for the test.)

I want to expand upon / comment on the other answers.

First, note that if you are using the noncentral t-distribution, you need to specify your equivalence bounds as a standardized mean difference $$\delta$$, since the noncentrality parameter will be $$\delta\sqrt{n}$$. If you are using TOST (i.e. central t-distribution), you need to specify it on the raw scale. In either case, if the equivalence bound is not on the right scale, you'll have to estimate it using the sample standard deviation, which will affect the error rates of your test.

Second, I do not follow David how using the noncentral t-distribution yields more power than TOST. It seems to me that (for a given sample size), this depends on the equivalence bounds (which the shape of the noncentral t depends on).

Here is the result of a simulation:

• x-axis shows the population effect size, y-axis is power
• solid curves are for TOST, dashed curves for noncentral t
• the colors correspond to different equivalence bounds of 0.5 (navy), 1 (purple), 1.5 (red), 2 (orange)
• vertical dotted lines show the corresponding equivalence bounds
• horizontal dotted line is at $$\alpha=0.05$$
• note that you can see David's result in the plot: the dashed navy line crosses the point (0,0.17) and the solid navy line (0,0)

In this simulation, equivalence bounds are specified as a standardized mean difference. For the TOST, the raw bounds are estimated using the sample SD. Hence, the max false positive error rate is not actually equal to $$\alpha$$ (solid curves do not necessarily go through where the dotted lines cross). However, if I didn't make an error somewhere, power for the noncnetral t is not necessarily higher than for TOST.

Code:

# sample size
n <- 10
# point null on raw scale
mu0 <- 0
# population sd
sigma <- 1.5
# equivalence bounds (magnitudes)
delta_eq <- c(0.5, 1, 1.5, 2)
# population ES = (mu-mu0)/sigma
delta <- seq(-max(delta_eq)*1.1, max(delta_eq)*1.1, length.out=200)

colors <- c("navy", "purple", "red", "orange")

par(mfrow=c(1,1))
nreps <- 10000

power.ump  <- matrix(nrow = length(delta), ncol = length(delta_eq))
power.tost <- matrix(nrow = length(delta), ncol = length(delta_eq))

for (i in seq_along(delta)) {

# population mean
mu <- mu0 + delta[i]*sigma

# each col hold p-values for a given delta_eq
p.tost <- matrix(nrow = nreps, ncol = length(delta_eq))
p.ump  <- matrix(nrow = nreps, ncol = length(delta_eq))

for (j in 1:nreps) {
# sample
x <- rnorm(n, mu, sigma)
m <- mean(x)
s <- sd(x)

# UMP

# non-centrality parameters (note: do not depend on s)
ncp.lower <- -delta_eq * sqrt(n)
ncp.upper <- +delta_eq * sqrt(n)

# observed t-statistic (two-sided test)
t <- abs((m-mu0) * sqrt(n) / s)

# p-value UMP
p.ump.lower <- pt(-t, n-1, ncp.lower, lower.tail=FALSE) - pt(t, n-1, ncp.lower, lower.tail=FALSE)
p.ump.upper <- pt(+t, n-1, ncp.upper, lower.tail=TRUE) - pt(-t, n-1, ncp.upper, lower.tail=TRUE)
p.ump[j, ] <- pmax(p.ump.lower, p.ump.upper)

# TOST

# lower and upper equivalence bounds TOST (raw scale, estimated using s)
mu0.tost.lower <- mu0 - delta_eq * s
mu0.tost.upper <- mu0 + delta_eq * s

# observed t-statistic TOST
t.tost.lower <- (m - mu0.tost.lower) * sqrt(n) / s
t.tost.upper <- (m - mu0.tost.upper) * sqrt(n) / s

# p-value tost
p.tost.lower <- pt(t.tost.lower, n-1, ncp=0, lower.tail=FALSE)
p.tost.upper <- pt(t.tost.upper, n-1, ncp=0, lower.tail=TRUE)
p.tost[j, ] <- pmax(p.tost.lower, p.tost.upper)

}

power.ump[i, ]  <- colMeans(p.ump < 0.05)
power.tost[i, ] <- colMeans(p.tost < 0.05)

}

plot(0, 0, type="n", xlim=range(delta), ylim=c(0,1), xlab=expression(delta), ylab="Power")
for (i in seq_along(delta_eq)) {
lines(delta, power.tost[ ,i], lty="solid", col=colors[i])
lines(delta, power.ump[ ,i], lty="dashed", col=colors[i])
abline(v = -delta_eq[i], col=colors[i], lty="dotted")
abline(v = delta_eq[i], col=colors[i], lty="dotted")
}
abline(h=0.05, lty="dotted")
`