* template do-file for two-way ANOVA
* assume a dataset with outcome y and factors factA and factB
version 15 /* also works with versions 13-14 */

* if the factors are coded by character variables, it is necessary to recode
* them as numerical variables; the new variables are called FactA and FactB
encode factA, gen(FactA)
encode factB, gen(FactB)

* standard two-way ANOVA table
anova y FactA FactB FactA#FactB
* with the ## notation, both the both main effects and the interaction are included
anova y FactA##FactB

* residuals analysis
predict pred, xb
predict stdres, rstandard
scatter stdres pred
qnorm stdres
* one example of normality test (recall that P-values are only approximate)
swilk stdres
* alternatively, the rvfplot command directly gives a plot of raw residuals versus predicted values
rvfplot

* interaction plot (factor A on x-axis), using add-on command: overlay
sort FactB FactA
overlay pred FactA, by(FactB) c(l)

* post-ANOVA inference using pwcompare command:
* means, SE and CI for any of the factors or interaction
pwcompare FactA FactB FactA#FactB, cimargins
* unadjusted pairwise comparisons for any of the factors or interaction
pwcompare FactA FactB FactA#FactB, pveffects
* bonferroni adjustment for multiple comparisons
pwcompare FactA FactB FactA#FactB, pveffects mcompare(bon)

* addition: one-way ANOVA with combined factor (for situations where the interaction is significant)
egen FactAB = group(FactA FactB)
* alternatively one may use a numeric generation of the combined factor, such as
*   FactAB=FactA+0.1*FactB
oneway y FactAB, bonferroni tabulate