How to report two-way anova results?

After finding a significant treatment in two-way anova, how do I report where the differences are? Every text I have read has left me at: Ah, the results are significant, so let's move on. In one-way anova I could use Tukey's HSD to find the means that differed, making my report much more effective. Right now I don't know how to do that with two-factor anova.

Here is an example from a text done in R:

Is there a difference in Friday tardy rates at different plants?

> mlate
day  plant absences
1   march-4 plant1       19
2  march-11 plant1       22
3  march-18 plant1       20
4   march-4 plant2       18
5  march-11 plant2       20
6  march-18 plant2       16
7   march-4 plant3       27
8  march-11 plant3       32
9  march-18 plant3       28
10  march-4 plant4       22
11 march-11 plant4       27
12 march-18 plant4       26

> anova(lm(absences ~ plant+day, data=mlate))
Analysis of Variance Table

Response: absences
Df  Sum Sq Mean Sq F value    Pr(>F)
plant      3 216.250  72.083  41.191 0.0002134 ***
day        2  30.167  15.083   8.619 0.0172128 *
Residuals  6  10.500   1.750
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


So I can see that both the plant and the day have an effect on worker's being late, but how do I compare the means like I would with a tukey's hsd in a one-way test?

-
In your example, you don't check for an interaction. If there is a strong theoretical reason that establishes there is no interaction, or if the data are sufficiently strong to establish that there is no interaction, then the effects of the two factors are independent of each other, and you can just as well run 2 one-way ANOVAs, with Tukey tests etc. –  gung Dec 28 '11 at 7:45