Title | Stata 7: Performing multiple comparison tests following two-way ANOVA | |

Author | Kenneth Higbee, StataCorp |

I am performing two-way ANOVA and wish to proceed with a multiple comparison test. Does Stata support any multiple comparison tests following two-way ANOVA?

Currently there is no official Stata command to perform multiple comparison tests for two-way (and higher order) ANOVA. It is, however, easy to construct the tests in Stata following a few simple steps:

- Run your
**anova**command. - Use
**test**to perform one of your single degree-of-freedom tests. The reported*p*-value is not corrected for multiple comparisons. - Apply a correction to account for the number of multiple comparisons
you are performing. See the "Methods and Formulas" section of
[R]
**oneway**for the appropriate correction. For instance, to obtain the Bonferroni adjusted*p*-value, multiply the uncorrected*p*-value by the total number of comparisons. (If the answer exceeds 1.0, the corrected*p*-value is reported as 1.0.) - Repeat steps (2) and (3) for each of your remaining tests.

You can either do these steps interactively in Stata, or you can write a do-file or ado-file to help automate the process.

Step (3), the correction for the multiple testing, is easy. The hardest
part is step (2), performing the single degree-of-freedom test. Sometimes
it can be difficult to figure out which linear combination of coefficients
you are going to test. Take a look at the section titled "Testing
coefficients" in [R] **anova** (page 59 -- Version 7 manual) and the FAQ
“How can I form various tests
comparing the different levels of a categorical variable after anova or
regress?” for guidance on obtaining single degree-of-freedom tests.

The simplest case is when each of your single degree-of-freedom tests use
the residual error from the ANOVA model as the denominator in the *F*
tests. When your single degree-of-freedom tests do not use residual error
as the denominator, step (2) requires more than one **test** command.
See the FAQ “How can you
specify a term other than residual error as the denominator in a single
degree-of-freedom *F* test after ANOVA” for guidance.

Here is a made-up dataset:

. lista b y 1. 1 1 4 2. 1 1 19 3. 1 1 17 4. 1 2 13 5. 1 2 9 6. 1 2 8 7. 2 1 34 8. 2 1 8 9. 2 1 11 10. 2 2 30 11. 2 2 38 12. 2 2 39 13. 3 1 39 14. 3 1 14 15. 3 1 10 16. 3 2 52 17. 3 2 51 18. 3 2 43 19. 4 1 36 20. 4 1 32 21. 4 1 50 22. 4 2 45 23. 4 2 35 24. 4 2 29

and here is the corresponding two-way ANOVA:

. anova y a b a*bNumber of obs = 24 R-squared = 0.7401 Root MSE = 9.52628 Adj R-squared = 0.6264 Source | Partial SS df MS F Prob > F -----------+---------------------------------------------------- Model | 4134.50 7 590.642857 6.51 0.0010 | a | 2470.16667 3 823.388889 9.07 0.0010 b | 580.166667 1 580.166667 6.39 0.0224 a*b | 1084.16667 3 361.388889 3.98 0.0270 | Residual | 1452.00 16 90.75 -----------+---------------------------------------------------- Total | 5586.50 23 242.891304

The factors a and b and their interaction are significant.

Pretend that when we set up this experiment, we had decided that we wished to examine four particular individual degree-of-freedom tests on factor a. The four tests of interest for factor a are

1 and 2 versus 3 and 4 1 versus 2 3 versus 4 1, 2, and 3 versus 4

While not related to the topic of single degree-of-freedom and multiple
comparison tests, looking at an interaction plot can be helpful in
understanding the significant interaction from this ANOVA. One way to
obtain an interaction plot is to use the
**graph** command after
the **adjust** command
with the **replace** option.
(**preserve** and
**restore** make it
easy to go to the cell mean data and then back to the original.)

. preserve . adjust , by(a b) replace------------------------------------------------------------------------- Dependent variable: y Command: anova ------------------------------------------------------------------------- ---------------------------- | b a | 1 2 ----------+----------------- 1 | 13.3333 10 2 | 17.6667 35.6667 3 | 21 48.6667 4 | 39.3333 36.3333 ---------------------------- Key: Linear Prediction

After the **adjust** command with the **replace** option, the data now
look like

. lista b y xb 1. 1 1 4 13.33333 2. 1 2 13 10 3. 2 1 34 17.66667 4. 2 2 30 35.66667 5. 3 1 39 21 6. 3 2 52 48.66667 7. 4 1 36 39.33333 8. 4 2 45 36.33333

One way to graph the interaction plot from this data is

. sort b a . graph xb a, c(L) s([b]) xlab ylab

but I prefer

. sort b a . gen b1 = xb if b==1 . gen b2 = xb if b==2 . graph b1 b2 a, c(ll) xlab ylab psize(200)

which produces the following Stata graph:

It appears that when factor a is equal to 2 or 3, the factor b effect is relatively different compared to when factor a is equal to 1 or 4.

The **restore** command brings us back to the original data.

. restore

Now, back to our original problem of producing multiple comparison tests. Remember that the four tests of interest for factor a are

1 and 2 versus 3 and 4 1 versus 2 3 versus 4 1, 2, and 3 versus 4

I find that using a cell means model is the easiest method for obtaining single degree-of-freedom tests. (See the FAQ How can I form various tests comparing the different levels of a categorical variable after ANOVA or regress? for details.)

The **group()** function of
**egen** produces a
variable indicating the cells of our design. The
**table** command then
shows the relationship between this new variable and the original a and b
variables for our example.

. egen z = group(a b) . table a b, c(mean z)---------------------- | b a | 1 2 ----------+----------- 1 | 1 2 2 | 3 4 3 | 5 6 4 | 7 8 ----------------------

Now, I run **anova** using the new variable, and I make sure to specify
the **noconstant** option.

. anova y z, noconstNumber of obs = 24 R-squared = 0.9397 Root MSE = 9.52628 Adj R-squared = 0.9095 Source | Partial SS df MS F Prob > F -----------+---------------------------------------------------- Model | 22616.00 8 2827.00 31.15 0.0000 | z | 22616.00 8 2827.00 31.15 0.0000 | Residual | 1452.00 16 90.75 -----------+---------------------------------------------------- Total | 24068.00 24 1002.83333

(Also of interest: I could now run **adjust, by(z)** and see that the
numbers agree with our original results.)

If I desire a Bonferroni correction on the four single degree-of-freedom
tests, I simply adjust the *p*-values for the tests by multiplying by
4. Here is the result for comparing levels 1 and 2 against levels 3 and 4
for factor a. (Remember how the z variable relates to the original factors
a and b? Level 1 of a corresponds to levels 1 and 2 of z; level 2 of a
corresponds to levels 3 and 4 of z; level 3 of a corresponds to levels 5 and
6 of z; and level 4 of a corresponds to levels 7 and 8 of z.)

. test _b[z[1]]+_b[z[2]]+_b[z[3]]+_b[z[4]] = _b[z[5]]+_b[z[6]]+_b[z[7]]+_b[z[8]]( 1) z[1] + z[2] + z[3] + z[4] - z[5] - z[6] - z[7] - z[8] = 0.0 F( 1, 16) = 19.48 Prob > F = 0.0004. return listscalars: r(drop) = 0 r(df_r) = 16 r(F) = 19.48393021120293 r(df) = 1 r(p) = .0004343772345068. di "corrected p-value = " min(1,r(p)*4)corrected p-value = .00173751

You can see what the **test** command returns in **r()** with the
command **return list**. We can take the uncorrected *p*-value of
0.00043 and produce the Bonferonni corrected *p*-value of .00174 by
multiplying **r(p)** (the uncorrected *p*-value) by 4 (the number of
tests to be performed).
The test for level 1 compared to level 2 of factor a is obtained next.

. test _b[z[1]]+_b[z[2]] = _b[z[3]]+_b[z[4]]( 1) z[1] + z[2] - z[3] - z[4] = 0.0 F( 1, 16) = 7.44 Prob > F = 0.0149. di "corrected p-value = " min(1,r(p)*4)corrected p-value = .05965724

Similarly, the test for level 3 compared to level 4 of factor a is obtained with

. test _b[z[5]]+_b[z[6]] = _b[z[7]]+_b[z[8]]( 1) z[5] + z[6] - z[7] - z[8] = 0.0 F( 1, 16) = 0.30 Prob > F = 0.5930. di "corrected p-value = " min(1,r(p)*4)corrected p-value = 1

Notice that in this case the corrected *p*-value is reported as 1 since
0.5930*4 is greater than 1.

Finally, the fourth test comparing levels 1, 2, and 3 against level 4 for factor a can be obtained with the following:

. test (_b[z[1]]+_b[z[2]]+_b[z[3]]+_b[z[4]]+_b[z[5]]+_b[z[6]])/3=_b[z[7]]+_b[z[8]]( 1) .3333333 z[1] + .3333333 z[2] + .3333333 z[3] + .3333333 z[4] + > .3333333 z[5] + .3333333 z[6] - z[7] - z[8] = 0.0 F( 1, 16) = 8.96 Prob > F = 0.0086. di "corrected p-value = " min(1,r(p)*4)corrected p-value = .03435797

Other multiple comparison tests could have been performed. One popular set of multiple comparison tests for a factor with four levels are the six tests of 1 versus 2, 1 versus 3, 1 versus 4, 2 versus 3, 2 versus 4, and 3 versus 4. The same approach will work for this situation and others.