Why is ANOVA not p-hacking? Say we have some data with many parameters. As an example let's say I'm an not-so-ethical journalist working for a food website and I'm looking to write some clickbait article "backed by science" about how some food or lifestyle is good/bad for you.
My data might contain a set of thousands of people, their socioeconomic status, what they had for breakfast, whether they're vegetarian, whether they prefer tea or coffee, etc. and also what their scores are for an IQ test.
If I split the dataset into many (100+) groups, where group one might split vegetarians/meat eaters, group two might split coffee drinkers and tea drinkers, group 3 might split female coffee drinkers and male tea drinkers etc. If then use ANOVA to compare all groups to see if there's a statistically significant difference, is this not effectively p-hacking?
I understand that ANOVA won't tell us which groups are different, but it seems like it would be a quick way to confirm that there IS a (false) positive in there. And then it's a case of finding out which groups are different.

There is no doubt that ANOVA--and any many other statistical tests--can be used for P-hacking, if used repeatedly on many random datasets, as you describe.
At the 5% level, about 1 in 20 of the tests will have a P-value below 0.05.
Under $H_0$ for an exact test based on a continuous test statistic, the distribution of the P-value is standard uniform, and thus has probability 0.05, below 0.05.
Shown below are simulations for (a) t.tests with normal data, (b) exact binomial tests with binomial counts, and (c) one-sample
Wilcoxon tests with normal data.
When $H_0$ is true: Only the t.test shows a standard uniform distribution for the P-value. The binomial test has P-values as close
to 5% as possible (without exceeding 5%), so there is no (nonrandomized) binomial test at exactly the 5% level. The Wilcoxon test is based on ranks; hence the test statistic is not continuous. But a test near the 5% level is available.

R code for figure:
set.seed(1234)
pv.t = replicate(10^5, 
         t.test(rnorm(25, 100, 15), mu = 100)$p.val)
mean(pv.t <= 0.05)
[1] 0.05017       # 5% within margin of sim error


set.seed(1235)
pv.b = replicate(10^5, 
        binom.test(rbinom(1, 10, .5), 10)$p.val)
mean(pv.b <= 0.05)
[1] 0.02173     # No test at exactly 5% available
                #  This test rejects for x=0,1,9, or 10
set.seed(1236)
pv.w = replicate(10^5, 
        wilcox.test(rnorm(25, 100, 15), mu = 100)$p.val)           
mean(pv.w <= 0.05)
[1] 0.04905       # test very nearly at 5%

par(mfrow=c(1,3))
 hist(pv.t, prob=T, col="skyblue2", main="T Test")
  curve(dunif(x), add=T, col="brown", lwd=2)
  abline(v=.05, col="red", lty="dotted", lwd=2)
 hist(pv.b, prob=T, col="skyblue2", main="Exact Binomial Test")
  curve(dunif(x), add=T, col="brown", lwd=2)
  abline(v=.05, col="red", lty="dotted", lwd=2)
 hist(pv.w, prob=T, col="skyblue2", main="Wilcoxon Test")
  curve(dunif(x), add=T, col="brown", lwd=2)
  abline(v=.05, col="red", lty="dotted", lwd=2)
par(mfrow=c(1,1))

