Use of post hoc tests following model selection based on AIC I am working on a dataset of vegetation characteristics recorded over three seasons in two different areas used by buffalo.
I have run generalised linear mixed models to determine whether vegetation biomass changes seasonally in different areas used by individual buffalo. I have run every possible combination of model and used AIC to identify the one with the best fit, which includes a fixed effect of season.
Within the best-fitting model, I would like to know which seasons differ from each other, and for this I used a post-hoc Tukey test from the multcomp package to compare pairs of seasons, and identified season pairs that differed significantly based on p-values.
I have been told that using the AIC model selection and the post-hoc tests based on p-values is mixing two types of analysis methods and is therefore not viable. Is this correct? 
I have searched the literature and cannot find any statistical papers that refer to this specific question, but I have found some research papers that have used this method. Can anyone direct me to relevant publications that can help me with this?
If indeed I cannot mix the two methods, does anyone have a suggestion for another way to look at differences between factor levels of a fixed effect?
Edit
Thank you for the response, but I just want to make sure that the question I am asking is clear. I got this response from a different forum, which may help to clarify the question that I am asking:
"This might depend on how you ran the model comparisons. If you had Season as e.g. a dummy-coded variable or something like that, and put all levels of it in as one step of the model (or e.g. just did it the R way, putting 'Season' in as its own predictor), then indeed model comparison doesn't tell you which levels of Season are different from each other, and therefore you still need some kind of post-hoc tests (or at least just inspecting the model coefficients, but indeed these should likely be corrected for multiple comparisons) to see which levels are different from which." 
I am trying to make sure that I am interpreting answers correctly, but for that to happen I have to be sure that my question is clear as well.
Many thanks for all the responses
Emily
 A: If you choose your model based on the actual observations (using aic, stepwise-regression etc.), your inferential statistics will be invalid for all the reasons gung has described here. Don't do it.

If indeed I cannot mix the two methods, does anyone have a suggestion
  for another way to look at differences between factor levels of a
  fixed effect?

Specify your model in advance, before observing any data.
Of course, you can conduct a pilot study and try different models on pilot data. This is completely valid. However, you cannot do valid inferential statistics on this dataset. Instead, specify your model based on your pilot data, and then collect independent new data on which you can perform inferential statistics with the prespecified model.
(Theoretically, you could possibly also account for your model selection in running inferential statistics on the same dataset you also selected your model for. However, I don't know how I would go about this, since this would in turn depend heavily on the specific dataset in question.)
