Help the Stat Consulting Group by giving a gift

How can I use the margins command with a 3-way anova interaction? (Stata 11)

The data will be analyzed as a 2x2x3 factorial anova. Variablesuse http://www.ats.ucla.edu/stat/stata/faq/threeway, clear

Indeed, it is the case that the 3-way interaction is statistically significant. Here is the plan for interpreting this interaction. First, we will run theanova y a##b##c /* equivalent to: anova y a b c a#b a#c b#c a#b#c */Number of obs = 24 R-squared = 0.9689 Root MSE = 1.1547 Adj R-squared = 0.9403 Source | Partial SS df MS F Prob > F -----------+---------------------------------------------------- Model | 497.833333 11 45.2575758 33.94 0.0000 | a | 150 1 150 112.50 0.0000 b | .666666667 1 .666666667 0.50 0.4930 c | 127.583333 2 63.7916667 47.84 0.0000 a#b | 160.166667 1 160.166667 120.13 0.0000 a#c | 18.25 2 9.125 6.84 0.0104 b#c | 22.5833333 2 11.2916667 8.47 0.0051 a#b#c | 18.5833333 2 9.29166667 6.97 0.0098 | Residual | 16 12 1.33333333 -----------+---------------------------------------------------- Total | 513.833333 23 22.3405797

The significant 3-way interaction suggests that there are one or more two-way interactions that are significant. We will select some two-way interactions to test for each of the levels of the third variable. If any of the two-way interactions are significant we will test the main effect of one variable at each level of the other variable. These are known as tests of simple main effects. Finally, if necessary, we follow these tests of simple main effects up with pair-wise comparisons among the appropriate means.

Here is the **margins** command.

Because we have some specific prior knowledge concerning the variables in this model we will begin by plotting the cell means for themargins a#b#c, post asbalancedAdjusted predictions Number of obs = 24 Expression : Linear prediction, predict() ------------------------------------------------------------------------------ | Delta-method | Margin Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- a#b#c | 1 1 1 | 11 .8164966 13.47 0.000 9.399696 12.6003 1 1 2 | 15 .8164966 18.37 0.000 13.3997 16.6003 1 1 3 | 19 .8164966 23.27 0.000 17.3997 20.6003 1 2 1 | 10.5 .8164966 12.86 0.000 8.899696 12.1003 1 2 2 | 10.5 .8164966 12.86 0.000 8.899696 12.1003 1 2 3 | 9.5 .8164966 11.64 0.000 7.899696 11.1003 2 1 1 | 10.5 .8164966 12.86 0.000 8.899696 12.1003 2 1 2 | 15.5 .8164966 18.98 0.000 13.8997 17.1003 2 1 3 | 18.5 .8164966 22.66 0.000 16.8997 20.1003 2 2 1 | 16.5 .8164966 20.21 0.000 14.8997 18.1003 2 2 2 | 20.5 .8164966 25.11 0.000 18.8997 22.1003 2 2 3 | 24 .8164966 29.39 0.000 22.3997 25.6003 ------------------------------------------------------------------------------

The graphs of the cell means suggest that the/* you can look at any of these matrices by using matrix list matname */ matrix means=e(b)' /* note: # is the kronecker product operator */ matrix gpha=(1\2)#(1\1\1\1\1\1) /* levels of a */ matrix gphb=(1\2\1\2)#(1\1\1) /* levels of b */ matrix gphc=(1\1\1\1)#(1\2\3) /* levels of c */ /* put them all together */ matrix gph= gpha,gphb,gphc,means matrix list gphgph[12,4] c1: c1: c1: c1 c1 c1 y1 r1:r1 1 1 1 11 r1:r2 1 1 2 15 r1:r3 1 1 3 19 r1:r4 1 2 1 10.5 r1:r5 1 2 2 10.5 r1:r6 1 2 3 9.5 r2:r7 2 1 1 10.5 r2:r8 2 1 2 15.5 r2:r9 2 1 3 18.5 r2:r10 2 2 1 16.5 r2:r11 2 2 2 20.5 r2:r12 2 2 3 24svmat gph twoway (connect gph4 gph3 if gph1==1 & gph2==1)(connect gph4 gph3 if gph1==1 & gph2==2), /// legend(order(1 "b=1" 2 "b=2")) ytitle(y) xtitle(c) xlabel(1 2 3) /// title("b#c interaction at a=1") name(bc_at_a1, replace) scheme(lean1)twoway (connect gph4 gph3 if gph1==2 & gph2==1)(connect gph4 gph3 if gph1==2 & gph2==2), /// legend(order(1 "b=1" 2 "b=2")) ytitle(y) xtitle(c) xlabel(1 2 3) /// title("b#c interaction at a=2") name(bc_at_a2, replace) scheme(lean1)

The second degree of freedom refers to this next table.c1 c2 c3 b1 + - b2 - +

These two tables translate into thec1 c2 c3 b1 + - b2 - +

To simplify adjustments for multiple testing we will use a Bonferroni adjustment across all
tests of simple main effects and pair-wise comparisons. Having already run through this
analysis before
writing this web page, I know that there will be seven such tests. Therefore I have included
a **display** command giving the adjusted p-values for each test.

Our guess from looking at the graph of the cell means was correct. The test of/* tests of b#c for each level of a */ /* test of b#c at a=1 */ test (1.a#1.b#1.c - 1.a#1.b#2.c - 1.a#2.b#1.c + 1.a#2.b#2.c = 0) /// (1.a#1.b#1.c - 1.a#1.b#3.c - 1.a#2.b#1.c + 1.a#2.b#3.c = 0)( 1) 1bn.a#1bn.b#1bn.c - 1bn.a#1bn.b#2.c - 1bn.a#2.b#1bn.c + 1bn.a#2.b#2.c = 0 ( 2) 1bn.a#1bn.b#1bn.c - 1bn.a#1bn.b#3.c - 1bn.a#2.b#1bn.c + 1bn.a#2.b#3.c = 0 chi2( 2) = 30.50 Prob > chi2 = 0.0000display "Bonferroni adjusted p-value: " %7.6f r(p)*7Bonferroni adjusted p-value: 0.000002/* test of b#c at a=2 */ test (2.a#1.b#1.c - 2.a#1.b#2.c - 2.a#2.b#1.c + 2.a#2.b#2.c = 0) /// (2.a#1.b#1.c - 2.a#1.b#3.c - 2.a#2.b#1.c + 2.a#2.b#3.c = 0)( 1) 2.a#1bn.b#1bn.c - 2.a#1bn.b#2.c - 2.a#2.b#1bn.c + 2.a#2.b#2.c = 0 ( 2) 2.a#1bn.b#1bn.c - 2.a#1bn.b#3.c - 2.a#2.b#1bn.c + 2.a#2.b#3.c = 0 chi2( 2) = 0.38 Prob > chi2 = 0.8290display "Bonferroni adjusted p-value: " %7.6f r(p)*7Bonferroni adjusted p-value: 5.803204

The test for differences among levels of/* tests of c for each level of b at a=1 */ /* test of c at b=1 & a=1 */ test (1.a#1.b#1.c - 1.a#1.b#2.c = 0)(1.a#1.b#1.c - 1.a#1.b#3.c = 0)( 1) 1bn.a#1bn.b#1bn.c - 1bn.a#1bn.b#2.c = 0 ( 2) 1bn.a#1bn.b#1bn.c - 1bn.a#1bn.b#3.c = 0 chi2( 2) = 48.00 Prob > chi2 = 0.0000display "Bonferroni adjusted p-value: " %7.6f r(p)*7Bonferroni adjusted p-value: 0.000000/* test of c at b=2 & a=2 */ test (1.a#2.b#1.c - 1.a#2.b#2.c = 0)(1.a#2.b#1.c - 1.a#2.b#3.c = 0)( 1) 1bn.a#2.b#1bn.c - 1bn.a#2.b#2.c = 0 ( 2) 1bn.a#2.b#1bn.c - 1bn.a#2.b#3.c = 0 chi2( 2) = 1.00 Prob > chi2 = 0.6065display "Bonferroni adjusted p-value: " %7.6f r(p)*7Bonferroni adjusted p-value: 4.245715

It appears that each of the pair-wise comparisons among levels of/* all pairwise comparisons between levels of c at b=1 a=1 */ test 1.a#1.b#1.c = 1.a#1.b#2.c /* c1 vs c2 */( 1) 1bn.a#1bn.b#1bn.c - 1bn.a#1bn.b#2.c = 0 chi2( 1) = 12.00 Prob > chi2 = 0.0005display "Bonferroni adjusted p-value: " %7.6f r(p)*7Bonferroni adjusted p-value: 0.003724test 1.a#1.b#1.c = 1.a#1.b#3.c /* c1 vs c3 */( 1) 1bn.a#1bn.b#1bn.c - 1bn.a#1bn.b#3.c = 0 chi2( 1) = 48.00 Prob > chi2 = 0.0000display "Bonferroni adjusted p-value: " %7.6f r(p)*7Bonferroni adjusted p-value: 0.000000test 1.a#1.b#2.c = 1.a#1.b#3.c /* c2 vs c3 */( 1) 1bn.a#1bn.b#2.c - 1bn.a#1bn.b#3.c = 0 chi2( 1) = 12.00 Prob > chi2 = 0.0005display "Bonferroni adjusted p-value: "%7.6f r(p)*7Bonferroni adjusted p-value: 0.003724

The content of this web site should not be construed as an endorsement of any particular web site, book, or software product by the University of California.