So my understanding is that the p-value is the likelihood of observing an effect at least as extreme as that shown in the sample data, if the null hypothesis is true.
But how is this a useful value? Even if the p-value is low, couldn't we have a similar, or even lower likelihood of observing an extreme effect if the null hypothesis is false?
Note this answer used to be much longer, but was also very verbose. So I have deleted all of it. If you like to read more then check version 2 of this post
Why do we need p-values? We need them for the cases when the experiment is subject to variability. Our measurements are not just clear examples of some observed effects. Often, measurements are noisy and random behavior is underlying observations. In such cases, we may wish to express/quantify how unlikely it would be that a certain effect is being observed if the null (no effect) hypothesis would be true. The p-value is an indication for the type I error probability.
In most settings, the estimated effect corresponds to the maximum likelihood estimate (or some method that is similar). So in this setting, there is no situation that "couldn't we have a similar or even lower likelihood..." because the effect is the effect/theory for which the likelihood is maximized.
The reason that we use the p-value is to express how precise the measurement/experiment is, by expressing how likely a measurement of this effect size could have occurred by pure chance. If the observed effect (or anything of similar size) could have likely occurred even when there is no effect present, then the experiment is not very accurate.
See for some related topics:
You're right - the p-value is a heuristic which needs to be interpreted in context. It should be interpreted as meaning how surprising the data would be if $H_0$ is true. So if p is small, alternative hypotheses look more attractive.
According to Bayes' Theorem:
$$ \mathbb{P}(H_0|D) = \frac{\mathbb{P}(D|H_0)\mathbb{P}(H_0)}{\mathbb{P}(D|H_0)\mathbb{P}(H_0)+\mathbb{P}(D|H_1)\mathbb{P}(H_1)} $$
So if you actually want to calculate the probability that the null hypothesis is true, your calculation has to include how likely the observed data $D$ is given $H_0$ and how likely $D$ is given $H_1$, and also a prior probability $\mathbb{P}(H_0)$ for how likely the null hypothesis is to be true.
If a study uses a p-value rather than Bayesian analysis to support rejection of a null hypothesis, that's probably fine so long as a) the data is sufficiently more likely under $H_1$ than under $H_0$ and b) $H_1$ is sufficiently plausible (so that $\mathbb{P}(H_1)$ isn't too small).
First let's take an example for reference. Assume that you are testing whether mean of a population is $0$ and sample mean is $\bar x$. Further, to keep things simple assume that the true distribution of population is $N(\mu,\sigma^2)$, with $\sigma$ known.
$$H_0: \mu = 0$$ $$H_a: \mu \ne 0$$
Now, you say:
Even if the p-value is low, couldn't we have a similar, or even lower likelihood of observing an extreme effect if the null hypothesis is false?
Hypothesis testing to check how confidently you can say that $\mu$ is equal to a certain value. That said, it is certainly possible that we may get low p-value (lower than $\alpha$) for a variety of possible values of $\mu$. In fact the set of all such values is the complement set of the confidence interval for $\mu$.
In the example above, we get confidence interval (at $5\%$ level of significance) by following equation:
\begin{align} Pr \Big(-1.96 <\frac{\bar x - \mu}{\sigma^2}<1.96\Big) &= 0.95 \\ Pr \Big(\bar x - 1.96\sigma^2 < \mu<\bar x +1.96\sigma^2\Big) &= 0.95 \end{align}
Now let, you conduct the following test: $$H_0: \mu = \theta$$ $$H_a: \mu \ne \theta$$
where, $\theta \in \mathbb R\backslash C$ and $C:= [\underline{c}, \bar c]$; $\underline{c}=\bar x - 1.96\sigma^2, \bar c =\bar x + 1.96\sigma^2$
As per definition, p-value for the above test is: \begin{align} p &= Pr \bigg(Z > \bigg|\frac{\bar x - \theta}{\sigma^2}\bigg| \bigg) + Pr \bigg(Z <- \bigg|\frac{\bar x - \theta}{\sigma^2}\bigg| \bigg) \tag{where $Z \sim N(0,1)$}\\ &= 2\bigg(1-\Phi\bigg(\bigg|\frac{\bar x - \theta}{\sigma^2}\bigg|\bigg)\bigg) \end{align}
Given the $\Phi(x)$ is increasing function and that either $\theta > \bar c$ or $\theta <\underline c$; $p<0.05$ for $\forall \theta \in \mathbb R\backslash C$.
Tests are normally defined within a model that is sufficiently flexible to accommodate all possible outcomes of the test statistic or the statistic on which the test is based. For example a model may be normal distributions with unknown mean and variance, and the null hypothesis may be "mean is zero". The test to be used here is a t-test, based on the mean standardised by its estimated standard error, and whatever value you observe will be compatible with a certain distribution in the model. So you may observe mean 15, variance 1, which is incompatible with the $H_0$: mean equals zero (i.e., producing a small p-value, unless the data set is very very small), but of course it is compatible with a normal distribution with mean 15 - and actually whatever mean you observe will be compatible with a normal distribution with that same mean value (or a mean value close to the observed one, as you could find out from computing a confidence interval), so there is always a distribution in the alternative that fits it. Meaning that if your alternative is flexible enough, a low p-value will not only indicate evidence against the $H_0$, but will also mean that the data are compatible with some other distribution in the alternative.
What can happen is that the data are not compatible with any normal distribution because their distributional shape is different. This however cannot be found out by the t-test, which by definition compares one normal distribution with other normal distributions. You'd need another test to find that out, with a nonparametric alternative including all kinds of distributional shapes - and again there will be some distributions (actually more than one, even in fact infinitely many, as the set of all distributions is very, very rich) that fits the data (of course not meaning that it is true - reality may just not be truly probabilistic).
PS: I should probably add that the involved notion of "compatibility" is relative to the test statistic, i.e., if your data look like being generated by a Gamma distribution with mean 15 and variance 1, from the point of view of the t-test (looking at means and variances) this is compatible with a normal(15,1), whereas if you look at things from the point of view of, say, a Shapiro-Wilk test for normality, this will turn out to be incompatible with normality if your sample size is large enough.
