Stata FAQ
How can I use the margins command to understand multiple interactions in regression and anova? (Stata 11)

The margins command, new in Stata 11, can be a very useful tool in understanding and interpreting interactions. We will illustrate the command in two examples using the hsbdemo dataset.

use http://www.ats.ucla.edu/stat/data/hsbdemo, clear

Example 1

We will begin with a model that has a categorical by categorical interaction (female by prog) along with a categorical by continuous interaction (honors by read). To keep things somewhat simple, the two interactions have no terms in common. We will begin by running the following regression model.

regress write female##prog honors##c.read

      Source |       SS       df       MS              Number of obs =     200
-------------+------------------------------           F(  8,   191) =   44.21
       Model |  11609.4878     8  1451.18597           Prob > F      =  0.0000
    Residual |  6269.38724   191   32.824017           R-squared     =  0.6493
-------------+------------------------------           Adj R-squared =  0.6347
       Total |   17878.875   199   89.843593           Root MSE      =  5.7292

------------------------------------------------------------------------------
       write |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
    1.female |   4.777377   1.751007     2.73   0.007     1.323584    8.231171
             |
        prog |
          2  |   2.284045   1.520393     1.50   0.135    -.7148712    5.282962
          3  |  -4.190142   1.769368    -2.37   0.019    -7.680154   -.7001312
             |
 female#prog |
        1 2  |  -2.493487   2.071853    -1.20   0.230    -6.580138    1.593164
        1 3  |   2.339986   2.387059     0.98   0.328    -2.368397    7.048369
             |
    1.honors |   28.17233   6.543047     4.31   0.000     15.26641    41.07824
        read |    .369414   .0553672     6.67   0.000     .2602043    .4786236
             |
      honors#|
      c.read |
          1  |  -.3200391   .1112185    -2.88   0.004    -.5394134   -.1006648
             |
       _cons |   28.80342   3.131345     9.20   0.000     22.62696    34.97988
------------------------------------------------------------------------------

As you can see the honor#c.read interaction is significant along with all the other one degree of freedom tests. There are two multi-degree of freedom tests that we need to follow up on using the testparm command.

testparm female#prog

 ( 1)  1.female#2.prog = 0
 ( 2)  1.female#3.prog = 0

       F(  2,   191) =    3.06
            Prob > F =    0.0493

testparm i.prog

 ( 1)  2.prog = 0
 ( 2)  3.prog = 0

       F(  2,   191) =    8.81
            Prob > F =    0.0002

The female#prog interaction is significant along with the two degree of freedom test for prog. Some people might call this the main effect for prog but that is not correct. Since we are using indicator (dummy) coding, the test for prog is really testing the effect of prog when female equals zero, that is, among males. If we multiply the F-ratio for prog by the numerator degrees of freedom, we get a value scaled like a chi-square. Thus, 2*r(F) = 17.616468 which is a value we will see again in a little while.

We can run the same model using the anova command. The anova will appear to be somewhat different because the model is parameterized differently but it is the exact same model.

anova write female##prog honors##c.read

                           Number of obs =     200     R-squared     =  0.6493
                           Root MSE      = 5.72922     Adj R-squared =  0.6347

                  Source |  Partial SS    df       MS           F     Prob > F
             ------------+----------------------------------------------------
                   Model |  11609.4878     8  1451.18597      44.21     0.0000
                         |
                  female |  922.674277     1  922.674277      28.11     0.0000
                    prog |  465.899583     2  232.949791       7.10     0.0011
             female#prog |  200.694275     2  100.347137       3.06     0.0493
                  honors |  608.523046     1  608.523046      18.54     0.0000
                    read |  424.568626     1  424.568626      12.93     0.0004
             honors#read |  271.796348     1  271.796348       8.28     0.0045
                         |
                Residual |  6269.38724   191   32.824017   
             ------------+----------------------------------------------------
                   Total |   17878.875   199   89.843593 

Note that the F-ratio for female#prog is the same as that from the testparm command and that the F-ratio for honors#read is the same as the t-value squared from the regression output ((-.3200391/.1112185)^2 = (-2.8775707)^2 = 8.2804133).

Next, we will use estimates store to save this model before using margins with the post option.

estimates store m1

We are finally ready to use the margins command to look at the female#prog interaction.

 margins female#prog, asbalanced post 

Predictive margins                                Number of obs   =        200

Expression   : Linear prediction, predict()
at           : female           (asbalanced)
               prog             (asbalanced)
               honors           (asbalanced)

------------------------------------------------------------------------------
             |            Delta-method
             |     Margin   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
 female#prog |
        0 1  |   53.82626    1.37582    39.12   0.000      51.1297    56.52281
        0 2  |    56.1103   .9608622    58.40   0.000     54.22705    57.99356
        0 3  |   49.63611   1.375616    36.08   0.000     46.93996    52.33227
        1 1  |   58.60363   1.256171    46.65   0.000     56.14158    61.06568
        1 2  |   58.39419   .9060001    64.45   0.000     56.61847    60.16992
        1 3  |   56.75348   1.213331    46.77   0.000     54.37539    59.13156
------------------------------------------------------------------------------

In case you have difficulty determining what each of the lines in the output above refers to, you can retype the margins command with the coeflegend option for more information.

 margins, coeflegend

Predictive margins                                Number of obs   =        200

Expression   : Linear prediction, predict()
at           : female           (asbalanced)
               prog             (asbalanced)
               honors           (asbalanced)

------------------------------------------------------------------------------
             |     Margin  Legend
-------------+----------------------------------------------------------------
 female#prog |
        0 1  |   53.82626  _b[0bn.female#1bn.prog]
        0 2  |    56.1103  _b[0bn.female#2.prog]
        0 3  |   49.63611  _b[0bn.female#3.prog]
        1 1  |   58.60363  _b[1.female#1bn.prog]
        1 2  |   58.39419  _b[1.female#2.prog]
        1 3  |   56.75348  _b[1.female#3.prog]
------------------------------------------------------------------------------

This margins syntax with the asbalanced option yields the "least-squares cell means" (SAS terminology), also known as the "estimated marginal cell means" (SPSS terminology), but more generally known as the adjusted cell means. And, because we used the post option, we can use the test command to compare differences in adjusted cell means.

/* test differences in prog at female==0 */

test (0.female#1.prog=0.female#2.prog)(0.female#1.prog=0.female#3.prog)

 ( 1)  0bn.female#1bn.prog - 0bn.female#2.prog = 0
 ( 2)  0bn.female#1bn.prog - 0bn.female#3.prog = 0

           chi2(  2) =   17.62
         Prob > chi2 =    0.0001

/* test differences in prog at female==1 */

test (1.female#1.prog=1.female#2.prog)(1.female#1.prog=1.female#3.prog)

 ( 1)  1.female#1bn.prog - 1.female#2.prog = 0
 ( 2)  1.female#1bn.prog - 1.female#3.prog = 0

           chi2(  2) =    1.75
         Prob > chi2 =    0.4170

The critical value of F for the per family error rate for these tests of simple main effects at alpha equals .05 is 3.71 which is equivalent to a chi-square value of 7.42. Using 7.42 as the critical value indicates that the test of differences in prog at female = 0 (males) was significant and has the same chi-square value that we calculated above in the 2*r(F). The test of prog at female equal one (females) was not significant. We should follow up on the significant test with pairwise comparisons at female equals zero.

/* pairwise comparisons for prog at female==0 */

test (0.female#1.prog=0.female#2.prog)  /* prog1 vs prog2 */

 ( 1)  0bn.female#1bn.prog - 0bn.female#2.prog = 0

           chi2(  1) =    2.26
         Prob > chi2 =    0.1330

test (0.female#1.prog=0.female#3.prog)  /* prog1 vs prog3 */

 ( 1)  0bn.female#1bn.prog - 0bn.female#3.prog = 0

           chi2(  1) =    5.61
         Prob > chi2 =    0.0179

test (0.female#2.prog=0.female#3.prog)  /* prog2 vs prog3 */

 ( 1)  0bn.female#2.prog - 0bn.female#3.prog = 0

           chi2(  1) =   17.60
         Prob > chi2 =    0.0000

These tests do not include any adjustments for multiple comparisons but we can use a Bonferroni adjustment by dividing our alpha level by the number of pairwise tests (.05/3 = .0167). With this (admittedly conservative) adjustment, only prog2 vs prog3 @ female==0 was statistically significant.

Next, we can turn our attention to the significant categorical by continuous interaction, honors by read. If you look back at the regression output you will see that the coefficient for read was .369414 with a standard error of .0553672. This value, .369414, is the slope of write on read when honors equals zero. We can easily obtain the slope when honors equals one by adding this coefficient to the coefficient for the interaction term (.369414 -.3200391 = .0493749). We can check this computation using the margins command after we use estimates restore to bring back our ANOVA/regression model.

estimates restore m1

/* slope of read for each level of honors */

margins, dydx(read) at(honors=(0 1)) vsquish

Average marginal effects                          Number of obs   =        200

Expression   : Linear prediction, predict()
dy/dx w.r.t. : read
1._at        : honors          =           0
2._at        : honors          =           1

------------------------------------------------------------------------------
             |            Delta-method
             |      dy/dx   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
read         |
         _at |
          1  |    .369414   .0553672     6.67   0.000     .2608963    .4779317
          2  |   .0493749    .099493     0.50   0.620    -.1456278    .2443776
------------------------------------------------------------------------------

These results are indeed the same as our computation of the slopes above. Of course now we also have standard errors and confidence intervals for both slopes.

Next, we will compute the predictive margins for every 10th value from 20 through 70 of read for each level of honors. The predictive margins for this model are the linear predictions of write for the six values of read for each level of honors.

 margins honors, at(read=(20(10)70)) vsquish

Predictive margins                                Number of obs   =        200

Expression   : Linear prediction, predict()
1._at        : read            =          20
2._at        : read            =          30
3._at        : read            =          40
4._at        : read            =          50
5._at        : read            =          60
6._at        : read            =          70

------------------------------------------------------------------------------
             |            Delta-method
             |     Margin   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
  _at#honors |
        1 0  |   38.53975   1.705084    22.60   0.000     35.19785    41.88165
        1 1  |   60.31129   4.090197    14.75   0.000     52.29465    68.32793
        2 0  |   42.23389   1.183996    35.67   0.000      39.9133    44.55448
        2 1  |   60.80504   3.121742    19.48   0.000     54.68654    66.92354
        3 0  |   45.92803   .7137843    64.34   0.000     44.52904    47.32702
        3 1  |   61.29879   2.177292    28.15   0.000     57.03138     65.5662
        4 0  |   49.62217   .4798269   103.42   0.000     48.68172    50.56261
        4 1  |   61.79254   1.309847    47.18   0.000     59.22529    64.35979
        5 0  |   53.31631   .7510556    70.99   0.000     51.84427    54.78835
        5 1  |   62.28629   .8188838    76.06   0.000     60.68131    63.89127
        6 0  |   57.01045   1.229244    46.38   0.000     54.60117    59.41972
        6 1  |   62.78004    1.26697    49.55   0.000     60.29682    65.26325
------------------------------------------------------------------------------

Since this is a linear model, each of the six predictive margins for honors = 0 will fall on a straight line as will the six values for honors = 1. If we want to graph these values as two lines we will need the values of the predictive margins, the values of read for which the values were computed and the value of honors to which they apply. The values for the predictive margins and for read are found in two different matrices saved by the margins command. The predictive margins are found in the matrix r(b) while the values of read are found in the matrix r(at) along with some other columns that we will discard.

Please note that if we has used the post option the two matrices would have been e(b) and e(at).

With a little bit of matrix work, we have the predictive margins and the read values in the Stata matrix b. Note the use of the Kronecker product to get two of each of the read values. The forvalues loop adds the alternating values of honors to matrix b. We finish up by saving the matrix to data with the svmat command followed by our graph twoway command.

/* graph regression lines for each level of honors */

matrix b=J(12,1,.)
matrix t=r(b)'
matrix b=t,b
matrix at=r(at)
matrix b=at[1...,8]#(1\1),b       /* here there be kronecker product */

forvalues i=1/12 {
  2.   mat b[`i',3]=~mod(`i',2)   /* 0/1 values for honors           */
  3. }     
  
matrix list b, nonames

b[12,3]
b[12,3]
         20  38.539748          0
         20  60.311292          1
         30  42.233888          0
         30  60.805041          1
         40  45.928028          0
         40   61.29879          1
         50  49.622168          0
         50  61.792539          1
         60  53.316308          0
         60  62.286288          1
         70  57.010448          0
         70  62.780037          1

svmat b, names(b)

graph twoway (connect b2 b1 if b3==0)(connect b2 b1 if b3==1), ///
   legend(order(1 "non-honors" 2 "honors")) ytitle(predicted write) xtitle(read) ///
   ylabel(20(10)70) title(Regression line by honors) 
   

After looking at the graph you might be interested in testing whether the predictive margins for honors = 0 are are different from the values for honors = 1 for each of the six values of read. If we had used the post option we could have followed up using test as a post-estimation command. However, it is easier to rerun the margins command to compute the marginal effect of honors using the dydx option. Since honors is a categorical variable margins will automatically compute the discrete change for us.

margins, dydx(honors) at(read=(20(10)70)) vsquish

Average marginal effects                          Number of obs   =        200

Expression   : Linear prediction, predict()
dy/dx w.r.t. : 1.honors
1._at        : read            =          20
2._at        : read            =          30
3._at        : read            =          40
4._at        : read            =          50
5._at        : read            =          60
6._at        : read            =          70

------------------------------------------------------------------------------
             |            Delta-method
             |      dy/dx   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
1.honors     |
         _at |
          1  |   21.77154   4.364615     4.99   0.000     13.21706    30.32603
          2  |   18.57115   3.298474     5.63   0.000     12.10626    25.03604
          3  |   15.37076   2.276819     6.75   0.000     10.90828    19.83325
          4  |   12.17037   1.400641     8.69   0.000     9.425166    14.91558
          5  |    8.96998   1.101633     8.14   0.000      6.81082    11.12914
          6  |   5.769589   1.714441     3.37   0.001     2.409347    9.129831
------------------------------------------------------------------------------
Note: dy/dx for factor levels is the discrete change from the base level.

All six of these tests would be significant using a Bonferroni adjusted critical value of .05/6= .0083.

Example 2

Our next example will be a little more complex in that it has two categorical by categorical interactions (female by prog and female by honors) with one term in common between them. In addition, there is a continuous covariate, math. This time we will start with the ANOVA model and follow it with the regression model.

anova write female##prog female##honors c.math

                           Number of obs =     200     R-squared     =  0.6304
                           Root MSE      = 5.88221     Adj R-squared =  0.6149

                  Source |  Partial SS    df       MS           F     Prob > F
           --------------+----------------------------------------------------
                   Model |  11270.2046     8  1408.77558      40.72     0.0000
                         |
                  female |  364.204296     1  364.204296      10.53     0.0014
                    prog |  407.481332     2  203.740666       5.89     0.0033
             female#prog |  228.269898     2  114.134949       3.30     0.0390
                  honors |  2398.83909     1  2398.83909      69.33     0.0000
           female#honors |  73.9785286     1  73.9785286       2.14     0.1453
                    math |  1072.64751     1  1072.64751      31.00     0.0000
                         |
                Residual |  6608.67039   191  34.6003686   
           --------------+----------------------------------------------------
                   Total |   17878.875   199   89.843593  

regress write female##prog female##honors math

      Source |       SS       df       MS              Number of obs =     200
-------------+------------------------------           F(  8,   191) =   40.72
       Model |  11270.2046     8  1408.77558           Prob > F      =  0.0000
    Residual |  6608.67039   191  34.6003686           R-squared     =  0.6304
-------------+------------------------------           Adj R-squared =  0.6149
       Total |   17878.875   199   89.843593           Root MSE      =  5.8822

------------------------------------------------------------------------------
       write |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
    1.female |   3.565842   1.786126     2.00   0.047     .0427751    7.088908
             |
        prog |
          2  |   .6948486   1.613388     0.43   0.667    -2.487499    3.877196
          3  |  -5.665851   1.785104    -3.17   0.002    -9.186901     -2.1448
             |
 female#prog |
        1 2  |   -.324647   2.152854    -0.15   0.880     -4.57107    3.921776
        1 3  |   4.833603   2.425983     1.99   0.048     .0484437    9.618763
             |
    1.honors |   11.23618   1.710336     6.57   0.000     7.862603    14.60975
             |
      female#|
      honors |
        1 1  |  -3.007322   2.056683    -1.46   0.145    -7.064052    1.049408
             |
        math |   .3262727   .0585993     5.57   0.000     .2106878    .4418577
       _cons |   31.69696   3.168456    10.00   0.000      25.4473    37.94662
------------------------------------------------------------------------------

testparm female#prog

 ( 1)  1.female#2.prog = 0
 ( 2)  1.female#3.prog = 0

       F(  2,   191) =    3.30
            Prob > F =    0.0390

testparm i.prog

 ( 1)  2.prog = 0
 ( 2)  3.prog = 0

       F(  2,   191) =    8.35
            Prob > F =    0.0003

estimates store m1

As you can see the regression and ANOVA models yield the same results for the interactions and one degree of freedom tests. The two degree of freedom test for prog is different from the anova results because regress uses indicator (dummy) coding. The testparm results for prog is actually the simple effect of prog when female is at its reference level of zero.

We will again use the margins command with the asbalanced and post options to obtain the adjusted cell means.

margins female#prog, asbalanced post

Predictive margins                                Number of obs   =        200

Expression   : Linear prediction, predict()
at           : female           (asbalanced)
               prog             (asbalanced)
               honors           (asbalanced)

------------------------------------------------------------------------------
             |            Delta-method
             |     Margin   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
 female#prog |
        0 1  |   54.49168   1.445088    37.71   0.000     51.65936      57.324
        0 2  |   55.18653   .9742899    56.64   0.000     53.27695     57.0961
        0 3  |   48.82583    1.43898    33.93   0.000     46.00548    51.64618
        1 1  |   56.55386   1.254649    45.08   0.000     54.09479    59.01292
        1 2  |   56.92406   .8203316    69.39   0.000     55.31624    58.53188
        1 3  |   55.72161   1.214488    45.88   0.000     53.34126    58.10196
------------------------------------------------------------------------------

Now we can use test commands to test the simple main effects for prog at each level of female.

/* test of prog at female==0 */
test (0.female#1.prog=0.female#2.prog)(0.female#1.prog=0.female#3.prog)

 ( 1)  0bn.female#1bn.prog - 0bn.female#2.prog = 0
 ( 2)  0bn.female#1bn.prog - 0bn.female#3.prog = 0

           chi2(  2) =   16.70
         Prob > chi2 =    0.0002

/* test of prog at female==1 */
test (1.female#1.prog=1.female#2.prog)(1.female#1.prog=1.female#3.prog)

 ( 1)  1.female#1bn.prog - 1.female#2.prog = 0
 ( 2)  1.female#1bn.prog - 1.female#3.prog = 0

           chi2(  2) =    0.66
         Prob > chi2 =    0.7174

The critical value for these tests of simple main effects is 3.76 for a per-family error rate of .05. Thus, only the test for prog at female==0 is statistically significant. We will follow up this significant test of simple main effects with pairwise comparisons among the levels of prog.

/* pairwise comparisons for prog at female==0 */
test (0.female#1.prog=0.female#2.prog)  /* prog1 vs prog2 */

 ( 1)  0bn.female#1bn.prog - 0bn.female#2.prog = 0

           chi2(  1) =    0.19
         Prob > chi2 =    0.6667

test (0.female#1.prog=0.female#3.prog)  /* prog1 vs prog3 */

 ( 1)  0bn.female#1bn.prog - 0bn.female#3.prog = 0

           chi2(  1) =   10.07
         Prob > chi2 =    0.0015

test (0.female#2.prog=0.female#3.prog)  /* prog2 vs prog3 */

 ( 1)  0bn.female#2.prog - 0bn.female#3.prog = 0

           chi2(  1) =   15.25
         Prob > chi2 =    0.0001

The Bonferroni adjusted p-values for prog1 versus prog3 and prog2 versus prog3 are .0045 and .0003, respectively. The other pairwise comparison was not significant without any adjustment.

Next, we need to look at the second interaction in the model. To do this we will use the estimates restore command. Once the estimates are restored, we will follow the same series of steps that we used for the first interaction.

estimates restore m1

margins female#honors, asbalanced post

Predictive margins                                Number of obs   =        200

Expression   : Linear prediction, predict()
at           : female           (asbalanced)
               prog             (asbalanced)
               honors           (asbalanced)

------------------------------------------------------------------------------
             |            Delta-method
             |     Margin   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
      female#|
      honors |
        0 0  |   47.21659   .7187761    65.69   0.000     45.80781    48.62536
        0 1  |   58.45276   1.573305    37.15   0.000     55.36914    61.53639
        1 0  |   52.28542   .7503273    69.68   0.000      50.8148    53.75603
        1 1  |   60.51427   1.137905    53.18   0.000     58.28402    62.74452
------------------------------------------------------------------------------

/* test of honors at female = 0 */

test 0.female#0.honors=0.female#1.honors

 ( 1)  0bn.female#0bn.honors - 0bn.female#1.honors = 0

           chi2(  1) =   43.16
         Prob > chi2 =    0.0000

/* test of honors at female = 1 */

test 1.female#0.honors=1.female#1.honors

 ( 1)  1.female#0bn.honors - 1.female#1.honors = 0

           chi2(  1) =   35.23
         Prob > chi2 =    0.0000

This time the critical value for the per family error rate is 5.10, so both tests are statistically significant.

Instead of running margins followed by test, we could have arrived at the same results by running margins with honors included in the dydx option. For categorical variables the dydx option calculates discrete change. The output for this approach is in terms of z-scores. By squaring the z-scores we can compare the results to the test command above.

estimates restore m1

margins female, dydx(honors) asbalanced post

Average marginal effects                          Number of obs   =        200

Expression   : Linear prediction, predict()
dy/dx w.r.t. : 1.honors
at           : female           (asbalanced)
               prog             (asbalanced)
               honors           (asbalanced)

------------------------------------------------------------------------------
             |            Delta-method
             |      dy/dx   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
1.honors     |
      female |
          0  |   11.23618   1.710336     6.57   0.000     7.883979    14.58837
          1  |   8.228854   1.386442     5.94   0.000     5.511478    10.94623
------------------------------------------------------------------------------
Note: dy/dx for factor levels is the discrete change from the base level.

display ([1.honors]_b[0.female]/[1.honors]_se[0.female])^2

43.159281

display ([1.honors]_b[1.female]/[1.honors]_se[1.female])^2

35.226969

Thus concludes example 2.