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:
- 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.
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.
|
