Hypothesis Testing: Different p-values from different test statistics this is a question that has been lingering in my mind for a very long time, but came up again in practice and I thought I'd reach out and ask about it.
The situation I'm interested in is when you are performing a Null Hypothesis Test and you can generate multiple test statistics with your data.  You know, both analytically and confirmed via simulation, the distribution of these different test statistics.  When you compute p-values using your actual observed data against these distributions, you get different results.  Similar results, it must be said, but different.  Different enough that were one to hew to the traditional critical values to reject a null hypothesis, one might make different decisions for the different test statistics.
A little bit more concretely, consider the Likelihood Ratio Test in addition to whatever usual test statistic you are going to look at.  For example, a t-statistic when comparing sample means. -2 times the natural log of the Likelihood Ratio is going to follow a Chi-Square(1) distribution, which can be confirmed via simulation or analytically.  The t-statistic is obviously going to follow a t-distribution.  With your actual observed data, you compute p-values using both test statistics and get different results.  Which one should be used?
More concretely, I was doing homework question 8.6 in Casella and Berger.  The problem setup there is X = sum of n exponential(theta) variables, and Y = sum of m exponential(mu) variables.  The null hypothesis is that the rates theta and mu are the same.  Sums of exponentials are gamma random variables, and the ratio X / (X + Y) of gamma random variables is a beta(n,m) random variable (I have confirmed this analytically and via simulation).  So I can use a beta(n,m) to generate p-values based on observed data.  I can also generate the likelihood ratio, which follows a chi-square(1) distribution, which I have confirmed via simulation.  The two test statistics lead to different p-values, and they are different enough that I felt compelled to reach out and ask this community for advice.
To summarize my question, when different test statistics can be used to generate p-values, which test statistic and p-value should be chosen?

Depending on what you know and one the assumptions you are willing to make about the population from which your sample was taken, there may be many different ways to test whether
$H_0: \mu=\mu_0$ against $H_a: \mu=\mu_a.$
Because the two tests may use different methods or distributions their powers may differ. That is, one may be more likely to detect the difference $\delta = 1$ between $\mu_0=1$ and $\mu_a = 0.$
If you have a symmetrical sample you might feel
comfortable assuming it its population is normal and using a t test. But for another symmetrical sample, you may use a nonparametric one-sample
Wilcoxon test, based on ranks of the data.
 set.seed(2022)
 x = rnorm(25, 0, 3)

 t.test(x, mu = 1)

         One Sample t-test

 data:  x
 t = -2.2657, df = 24, p-value = 0.03277
alternative hypothesis: true mean is not equal to 1
95 percent confidence interval:
 -1.5620832  0.8805949
sample estimates:
  mean of x 
  -0.3407441 

wilcox.test(x, mu = 1)

        Wilcoxon signed rank test

 data:  x
 V = 93, p-value = 0.06263
 alternative hypothesis: true location is not equal to 1

The t test is significant at the 5% level because
its P-value is below $0.05,$ but the Wilcoxon test
is not significant at the 5% level. One test detects
the difference between $\mu = 1$ and $\mu = 0;$ the other does not. Nothing is "wrong" with either test. They are based on different assumptions and use test statistics with different distributions, so it is not surprising that they give different "answers."
The following simulations in R show that the t test has somewhat better power in this situation than does the Wilcoxon test.
 set.seed(630)
 pv.t = replicate(10^5, t.test(rnorm(25,0,3), mu = 1)$p.val)
 mean(pv.t <= 0.05)
 [1] 0.361          # Aprx power t test 36%

 pv.w = replicate(10^5, wilcox.test(rnorm(25,0,3), mu = 1)$p.val)
 mean(pv.w <= 0.05)
 [1] 0.34359b       # Aprx power Wilcoxon SR test 34%

