Search
   >> Home >> Resources & support >> FAQs >> Performing multiple comparison tests following two-way ANOVA

Stata 7: Does Stata support any multiple comparison tests following two-way ANOVA?

Title   Stata 7: Performing multiple comparison tests following two-way ANOVA
Author Kenneth Higbee, StataCorp
Date August 2001

Question:

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?

Answer:

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:

  1. Run your anova command.
  2. Use test to perform one of your single degree-of-freedom tests. The reported p-value is not corrected for multiple comparisons.
  3. 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.)
  4. 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.

Example

Here is a made-up dataset:

 . list
 
              a          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*b
        
                    Number 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

 . list
 
              a          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, noconst
 
                    Number 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 list
 
 scalars:
               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.

The Stata Blog: Not Elsewhere Classified Find us on Facebook Follow us on Twitter LinkedIn Google+ Watch us on YouTube