My question is simple: when calculating p-values, why do we consider the distribution of the sample mean (i.e., why is the sample mean our 'test statistic')? Why not use the median, mode, lower quartile, or any other function of the data? For that matter, why not do the natural thing and just compute the joint density $f(x1,...,x_n)$?
So, you have entered into two deep questions. The first is “why the mean?” The second is implicit in your joint distribution question, “why use the sampling distribution?”
I will answer the second before the first because the answer of the first depends somewhat upon the second.
There are two broad families of statistical methodologies, the Bayesian methods, and the null hypothesis methods. Null hypothesis methods are sampling based methods. Bayesian methods are generative methods.
Because of the nature of your question, I am going to overgeneralize a little bit to make the answer shorter than a book.
There are several reasons to use one class of problem over the other, and there are circumstances where you would only use one specific class of solution. They are not interchangeable.
If one computed the joint density, $f(x_1\dots{x_n})$, one would fail because a probability can only exist with respect to a model. Implicitly, you have a model in mind with the joint density, but it does matter to call it out. In a Bayesian framework, that joint density would be the density averaged over all possible models and all possible parameters.
There are a couple of specific reasons to use Bayesian methods. The first is if you have prior knowledge. For example, let us assume that you know $\mu=1,000,000$ is almost impossible, but $\mu=5$ is quite probable. That is the definition of prior information. You know something. It might be vague. It might be unbelievably vague, but it is also not true ignorance either.
If you have prior knowledge, then non-Bayesian methods are inadmissible solutions. A solution is considered admissible if the estimator cannot be stochastically dominated by another estimator. Now admissibility is only one possible criterion. It implies that for a given data set, the Bayesian method will be intrinsically the least risky to use and the most accurate when averaged over the sample space.
A second reason to use a Bayesian method is that they generate coherent strategies. A strategy is coherent if a bookie that is taking all finite bets at the bookie’s stated prices, cannot be gamed into a sure loss by a crafty actor or set of actors. Null hypothesis methods are not coherent. If you use a null hypothesis method in a game, they can create sure losses for a bookie.
Consider the case of an inverse roulette wheel, where the ball is dropped, and the wheel spun. Instead of betting what number it will land on, the croupier sees where the ball lands, but the players cannot see it. Instead, the croupier tosses two fair coins. If a coin comes up heads, then the croupier reveals the value one place to the left of the true value; otherwise, the croupier reveals one value to the right.
The sample space is $\{\{L,L\},\{R,R\},\{R,L\},\{L,R\}\}$. To make this concrete, consider a wheel where the number came up 17. They will either be given the set, $\{16,16\}$,$\{16,18\}$, or $\{18,18\}$. Obviously, if you receive $\{16,18\}$ you should bet 17, but what should you do otherwise?
If, for example, you should receive $\{16,16\}$, Bayesian decision theory says you should toss a fair coin and either bid on 15 or 17. You will win 75% of the time. Frequentist decision theory says that the minimum variance unbiased estimator is 16. You can check that for yourself. You should bid 16. Of course, sixteen is an impossible value, known with probability one. You will win 50% of the time.
If the bookie placed odds using Frequentist decision theory but you played as a Bayesian, you would rake in the cash. The Frequentist odds are not coherent.
Although not obvious yet, this is the start of the answer to the first question on the mean. And, in fact, you do want to be careful. If you have a Bayesian model, you do not want to test it with a null hypothesis method. Conversely, if you have a Frequentist model, you do not want to test it with Bayesian methods. They would produce different predictions.
So far, we have two reasons to use the joint distribution via Bayes theorem. The first is accuracy. The second is coherence. So why use a null hypothesis method?
The first reason to choose a null hypothesis method is that you have no prior information. The nice thing about null hypothesis methods is that they work equally well regardless of the true value of the parameter. That is not true for a Bayesian method. They do not work equally well for all possible values of a parameter.
In the example above, where you felt you could all but rule out $\mu=1,000,000$ would be problematic if $\mu=1,000,000$ as you would have underweighted the possibility of the case. Bayesian methods weight the sample based on prior knowledge.
As well, they minimize the maximum possible loss you could face from drawing an unfortunate sample. The minimum variance unbiased estimator and the maximum likelihood estimator generate sampling distributions that are worst-case distributions.
Bayesian methods, generally speaking, minimize the average loss you would expect to receive from drawing an unfortunate sample. In the absence of prior information, it is challenging to minimize an average loss when you have nothing to compare to the sample.
The second reason is that you truly have a sharp null hypothesis. Bayesian methods have no clean or nice way to handle a hypothesis of the form $\mu=k$. This is because samples are not random in Bayesian thinking. You saw the sample. It is a fact. It cannot be random as there is no uncertainty in it. Parameters are uncertain in Bayesian thought. You cannot see them; you cannot know what they are. The $\Pr(\mu=k)$ is intrinsically zero. It is one point over an infinite number of points. The measure of a single point from an infinite number is zero. On the other hand, Bayesian methods handle solutions over regions fine.
The third reason is that you need an unbiased estimator for some reason. Assuring an estimator is unbiased is expensive. There can be a huge loss of information to assure it is unbiased. You may need a sample size many times larger than a biased Bayesian estimator to get the same result. All Bayesian estimators are biased estimators. Bayesian methods optimize the trade-off between bias and accuracy, though as the sample size becomes arbitrarily large, the Bayesian method will converge to the true parameters, so it becomes a just barely biased estimator.
These differences answer the question of using the joint distribution of the data to determine the probability versus using the sampling distribution of the model to determine a probability.
If you decide you want to use a Bayesian method, then you solve $\Pr(\mu\le{k}|{X})$, for example, instead of $\Pr({X}|\mu\le{k})$. The converse is true for the null hypothesis method.
Now as to using the mean for a sampling distribution. There are several reasons for this. The first is that in most standard cases, the mean is the most efficient estimator of the parameter. Secondly, all unbiased estimation has quadratic loss buried inside of it somewhere. Even if you are making an estimate based on the median, that is the mean rank. While the median minimizes absolute linear loss, it minimizes quadratic loss in rank space.
That takes us back to stochastic dominance. If you have two unbiased estimators, say the sample mean, and the sample median, but the sampling distribution of the sample mean stochastically dominates the sampling distribution of the sample median, then why would you use the median?
You see the median begin to be used, generally, when robustness is a concern rather than power.
If you want to know the probability your model is true, given the data that you saw, then you use a Bayesian method. If you want to know the probability of seeing data as extreme or more extreme than likely under your hypothesized value, then you use a null hypothesis method. If you use a Bayesian method, then it depends on the joint density of the actual sample. If you use a null hypothesis method, then you will use the sampling distribution of the model to condition the data. The mean will tend to be used due to issues of efficiency and power and not any other statistic.
EDIT
If your hypothesis concern is model selection, then you should use a Bayesian method.
Bayesian methods treat a model as a parameter. It uses the data to determine which model is most probable. You can have any number of models you like and any number of parameters.
Model selection by information criterion maps to stylized Bayesian posterior densities, but do not make their assumptions explicit. Bayesian methods allow you to let the data select the model.
EDIT
For clarity, let us imagine that we have two possible models that data are drawn from. Model one is $$\frac{1}{\pi}\frac{1}{1+(x-\mu)^2}.$$ Model two is $$\lambda\exp{[-\lambda{x}]}.$$
$$\Pr(M_1;\mu|X)=\frac{\prod_{i=1}^n\frac{1}{\pi}\frac{1}{1+(x-\mu)^2}\Pr(\mu|M_1)\Pr(M_1)}{\int_{\mu\in\Theta}\prod_{i=1}^n\frac{1}{\pi}\frac{1}{1+(x-\mu)^2}\Pr(\mu|M_1)\Pr(M_1)\mathrm{d}\mu+\int_{\lambda\in\Theta}\prod_{i=1}^n\lambda\exp{[-\lambda{x}]}\Pr(\lambda|M_2)\Pr(M_2)\mathrm{d}\lambda}$$
$$\Pr(M_2;\lambda|X)=\frac{\prod_{i=1}^n\lambda\exp{[-\lambda{x}]}\Pr(\lambda|M_2)\Pr(M_2)}{\int_{\mu\in\Theta}\prod_{i=1}^n\frac{1}{\pi}\frac{1}{1+(x-\mu)^2}\Pr(\mu|M_1)\Pr(M_1)\mathrm{d}\mu+\int_{\lambda\in\Theta}\prod_{i=1}^n\lambda\exp{[-\lambda{x}]}\Pr(\lambda|M_2)\Pr(M_2)\mathrm{d}\lambda}$$
The first probability statement $\Pr(M_1;\mu|X)$ can be reduced to $\Pr(M_1|X)$ by marginalization so that $$\Pr(M_1|X)=\int_{\mu\in\Theta}\Pr(M_1;\mu|X)\mathrm{d}\mu.$$
The same can be done for model two so that you can find the probability of model two as $$\Pr(M_2|X)=\int_{\lambda\in\Theta}\Pr(M_2;\lambda|X)\mathrm{d}\lambda.$$
Since $\Pr(M_1|X)+\Pr(M_2|X)=1$, it follows that these probability statements are the probability that these are the true models in nature given the data you have seen.
In the first model $\mu$ which is the median and the mode as this distribution has no mean or variance vanishes from the probability statement in the last line. A point value is never calculated because it is marginalized out. Likewise, $\lambda$ which is $1/mean$ is also never calculated as a point statistic.
These extra parameters are called nuisance parameters and their integral has to be solved, but only to make the integral vanish.
When solving Bayesian model selection, the parameters are nuisances and you need them as placeholders for calculation purposes, but you have to make them vanish. If you didn't, then you could have points where model one was more likely than model two, but other points where that relationship is reversed.
Imagine that $\Pr(M_1|X)=99.99\%$, it would follow that $\Pr(M_2|X)=.01\%$ and you could ignore model 2. I once tested seventy-eight models of bankruptcy. One model had a posterior probability of about 57%. the second model had a posterior probability of about 43%. The remaming 76 models had a combined probability of $1/10000$ of one percent. I could exclude the 76 models, but I couldn't exclude two models as both were approximately equally likely to be the data generating function.
Although I bypassed the problem using model averaging, that is I averaged the two regressions together by weighting one result by .57 and the other by .43, what it really meant was that I was close, but not quite there and that some similar convex combination of the two sets of variables were what mattered.
