Hypothesis test: numeric vs. ranked

Question

I believed that the most powerful hypothesis test for judging whether a single sample comes from $N(0,1)$ or from $N(1,1)$ uses the average value as test statistics. Thus, I calculate the sample size using a Monte Carlo simulation and compared two methods: (a) using the numeric average value as test statistic, and (b) transforming the data first to ranks and then calculating the average value. What I found is that the ranked data yields a smaller sample size. How is this possible?

What I am finally interested in is to determine the sample size for a future experiment. I like to take $n$ samples from an industrial process with a known distribution, analyse the samples and show that the process is well within its specification limits. This allows me to produce the product without measuring each part.

There are many subtle difficulties involved in the hypothesis test. That's why I setup the following toy model to focus the attention on the above described "numeric vs. ranked" question.

set.seed(2020)
transformToRank = FALSE
alpha           = 0.05
betaTarget      = 0.20

T = matrix(nrow = 1e4, ncol = 2) # alloc space for test statistic
for ( nSample in 3:100 ) {
    for ( i in 1:1e4 ){
        # combine data:
        random = c(rnorm(nSample,0,1), rnorm(nSample,1,1))
        if ( transformToRank ){
            # rank data (combined):
            random = rank(random)
        } 
        # Calc test statistic:
        T[i,] = c(mean(random[1:nSample]), mean(random[(nSample+1):(2*nSample)]) );
    }
    cutoff1 = quantile(T[, 1], alpha/2)
    cutoff2 = quantile(T[, 1], 1-alpha/2)
    beta    = sum(cutoff1 < T[, 2] & T[, 2] < cutoff2) / 1e4
    if ( beta <= betaTarget ){
        break 
    } else {
        nSample = nSample + 1
    }
}
print(beta)
print(nSample)

The numerical case can be checked by using the following code:

library(BSDA)
set.seed(2020)
pValue = replicate(1e4,z.test(rnorm(8,1,1), 
                    alternative="two.sided", mu=0, sigma.x=1)$p.value); 
power = mean(pValue <= 0.05)
beta  = 1 - power
print(beta)

If I use the numeric value, the sample size $n=8$ satisfies the $\beta$-risk condition. In contrast, if I transform the random numbers to ranks, I only need $n=4$. Thus, if I am willing to use a randomly generated dataset to analyse (rank) the experimental data, the power of the test increases significantly. This logic also applies if I sample from a location scaled version of the $t$-distribution. What am I missing?

My key question is, why is the ranked version superior in power by such a huge amount? I believed that I'll get approximately the same sample size, because it is known that the rank transformation provides a bridge to the non-parametric hypothesis tests, see e.g. Ref1, or Ref2.

Answer 1

Here are simulations comparing two samples of size 15 from $\mathsf{Norm}(0,1)$ and $\mathsf{Norm}(1,1),$ respectively. My simulation shows that the pooled t test has better power than the two-sample Wilcoxon test, which is well-known, and that neither test has power $0.8.$

set.seed(2020)
pv = replicate(10^4, t.test(rnorm(15,0,1),
                            rnorm(15,1,1), var.eq=T)$p.val)
mean(pv <= 0.05)
[1] 0.7525

set.seed(911)
pv = replicate(10^4, wilcox.test(rnorm(15,0,1),
                                 rnorm(15,1,1))$p.val)
mean(pv <= 0.05)
[1] 0.7118

It seems that I have misunderstood what you are doing, that your simulation code is wrong, or both. It might be helpful to have a clearer explanation of what you are doing with ranks, and to see the inner loop of your program where you compute power.

[It makes no sense to take averages of ranks for the two samples separately: for example, if $n=15,$ then both sets of ranks would run from 1 through 15 and both sets of ranks would always sum to 120. You might want to look at what the Wilcoxon rank-sum test does with ranks of the two samples.]

Here are simulations with sample sizes $n=25$ and difference $0.5$ in population means. In neither case is power anywhere near 80%.

set.seed(1066)
pv = replicate(10^4, t.test(rnorm(25,0,1),
                            rnorm(25,.5,1), var.eq=T)$p.val)
mean(pv <= 0.05)
[1] 0.3978

set.seed(1776)
pv = replicate(10^4, wilcox.test(rnorm(25,0,1),
                                 rnorm(25,.5,1))$p.val)
mean(pv <= 0.05)
[1] 0.3867

Note: For pooled t tests, here is an online 'power and sample size' calculator, that works for reasonable parameters.