Graphing results from the margins command can help in the interpretation of your model. Almost all of the needed results will be found in various matrices saved by margins. If you do not use the post option the matrices of interest are r(b) which has the point estimates, r(V) the variance-covariance matrix (standard errors are squares roots of the diagonal values), and r(at) which has the values used in the at option. If you do use the post option then the relevant matrices are e(b), e(V) and e(at).
Generally, the point estimates will be plotted on the y-axis. The x-axis will either be the at values or levels of a categorical predictor. The general strategy shown here will show how to get the necessary values for graphing into a Stata matrix which we will then save to data in memory using the svmat command. Once the values are part of the data we can plot the points using standard graph commands.
Example 1
Let’s start with a model that has a categorical by categorical interaction along with two continuous covariates. We will run the model as an anova but we would get the same results if we ran it as a regression.
use https://stats.idre.ucla.edu/stat/data/hsbdemo, clear
anova write prog##female math read
Number of obs = 200 R-squared = 0.6932
Root MSE = 6.59748 Adj R-squared = 0.5155
Source | Partial SS df MS F Prob > F
------------+----------------------------------------------------
Model | 12394.507 73 169.787767 3.90 0.0000
|
prog | 22.3772184 2 11.1886092 0.26 0.7737
female | 1125.44586 1 1125.44586 25.86 0.0000
prog#female | 287.95987 2 143.979935 3.31 0.0398
math | 2449.19165 39 62.799786 1.44 0.0670
read | 2015.49976 29 69.4999916 1.60 0.0411
|
Residual | 5484.36799 126 43.5267301
------------+----------------------------------------------------
Total | 17878.875 199 89.843593
Next, we run the margins command to get the six adjust cell means from the 2 by 3 interaction. These adjusted cells means are called lsmeans in SAS or estimated marginal means in SPSS. After the margins command we will list the r(b) matrix.
margins prog#female, asbalanced
Adjusted predictions Number of obs = 200
Expression : Linear prediction, predict()
at : prog (asbalanced)
female (asbalanced)
math (asbalanced)
read (asbalanced)
------------------------------------------------------------------------------
| Delta-method
| Margin Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
prog#female |
1 0 | 50.43791 1.969439 25.61 0.000 46.57788 54.29794
1 1 | 56.9152 1.948113 29.22 0.000 53.09697 60.73343
2 0 | 52.58387 1.408514 37.33 0.000 49.82323 55.34451
2 1 | 55.05205 1.41875 38.80 0.000 52.27136 57.83275
3 0 | 47.81982 2.108078 22.68 0.000 43.68806 51.95158
3 1 | 57.29057 1.785628 32.08 0.000 53.79081 60.79034
------------------------------------------------------------------------------
matrix list r(b)
r(b)[1,6]
1.prog# 1.prog# 2.prog# 2.prog# 3.prog# 3.prog#
0.female 1.female 0.female 1.female 0.female 1.female
r1 50.437911 56.9152 52.583868 55.052054 47.81982 57.290573
Note that r(b) is a 1×6 matrix. We will want the 6×1 transpose of r(b) which we will call b.
matrix b=r(b)'
matrix list b
b[6,1]
r1
1.prog#
0.female 50.437911
1.prog#
1.female 56.9152
2.prog#
0.female 52.583868
2.prog#
1.female 55.052054
3.prog#
0.female 47.81982
3.prog#
1.female 57.290573
If you look back at the margins output you will see a column for prog that has two 1’s followed by two 2’s and then two 3’s. We will create this array using a Kronecker product of a column vector (123) and a column vector (11).
matrix prg=(123)#(11) /* # is the kronecker product operator */
matrix list prg
prg[6,1]
c1:
c1
r1:r1 1
r1:r2 1
r2:r3 2
r2:r4 2
r3:r5 3
r3:r6 3
We will use this same trick again to create an series of alternating 0’s and 1’s. Next, we will combine the three columns together into a matrix called c.
matrix fem=(111)#(01) /* # is the kronecker product operator */
matrix list fem
fem[6,1]
c1:
c1
r1:r1 0
r1:r2 1
r2:r3 0
r2:r4 1
r3:r5 0
r3:r6 1
matrix c=prg,fem,b
matrix list c
c[6,3]
c1: c1:
c1 c1 r1
r1:r1 1 0 50.437911
r1:r2 1 1 56.9152
r2:r3 2 0 52.583868
r2:r4 2 1 55.052054
r3:r5 3 0 47.81982
r3:r6 3 1 57.290573
Now we can save c to data using the svmat command.
svmat c, names(c)
clist c1-c3 in 1/7
c1 c2 c3
1. 1 0 50.43791
2. 1 1 56.9152
3. 2 0 52.58387
4. 2 1 55.05206
5. 3 0 47.81982
6. 3 1 57.29057
7. . . .
We can now plot the adjusted cell means for prog by female using graph twoway
/* plot prog by female */
graph twoway (connect c3 c1 if c2==0)(connect c3 c1 if c2==1), ///
xlabel(1 2 3) legend(order(1 "female=0" 2 "female=1")) ///
xtitle(prog) ytitle(ajdusted cell means) ///
name(prog_by_female, replace)
We can also graph the results for female by prog just by reversing c1 and c2 in the graph twoway command.
/* plot female by prog */
graph twoway (connect c3 c2 if c1==1)(connect c3 c2 if c1==2) ///
(connect c3 c2 if c1==3), xlabel(0 1) ///
legend(order(1 "prog=1" 2 "prog=2" 3 "prog=3")) ///
xtitle(female) ytitle(ajdusted cell means) ///
name(female_by_prog, replace)
Example 2
For our second example we will graph the results of a categorical by continuous interaction from a logistic regression model.
use https://stats.idre.ucla.edu/stat/data/logitcatcon, clear
logit y i.f##c.s, nolog
Logistic regression Number of obs = 200
LR chi2(3) = 71.01
Prob > chi2 = 0.0000
Log likelihood = -96.28586 Pseudo R2 = 0.2694
------------------------------------------------------------------------------
y | Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
1.f | 5.786811 2.302518 2.51 0.012 1.273959 10.29966
s | .1773383 .0364362 4.87 0.000 .1059248 .2487519
|
f#c.s |
1 | -.0895522 .0439158 -2.04 0.041 -.1756255 -.0034789
|
_cons | -9.253801 1.94189 -4.77 0.000 -13.05983 -5.447767
------------------------------------------------------------------------------
We will use the margins command to get the predicted probabilites for 11 values of s from 20 to 70 for both f equal zero and f equal one. The vsquish option just reduces the number of blank lines in the output.
margins f, at(s=(20(5)70)) vsquish
Adjusted predictions Number of obs = 200
Model VCE : OIM
Expression : Pr(y), predict()
1._at : s = 20
2._at : s = 25
3._at : s = 30
4._at : s = 35
5._at : s = 40
6._at : s = 45
7._at : s = 50
8._at : s = 55
9._at : s = 60
10._at : s = 65
11._at : s = 70
------------------------------------------------------------------------------
| Delta-method
| Margin Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
_at#f |
1 0 | .0033115 .0040428 0.82 0.413 -.0046123 .0112353
1 1 | .1529993 .098668 1.55 0.121 -.0403864 .3463851
2 0 | .0079995 .0083179 0.96 0.336 -.0083033 .0243023
2 1 | .2188574 .1104097 1.98 0.047 .0024583 .4352564
3 0 | .0191964 .0164503 1.17 0.243 -.0130456 .0514383
3 1 | .3029251 .1126363 2.69 0.007 .082162 .5236882
4 0 | .0453489 .0304422 1.49 0.136 -.0143166 .1050145
4 1 | .4026401 .1026147 3.92 0.000 .201519 .6037612
5 0 | .1033756 .0500784 2.06 0.039 .0052238 .2015275
5 1 | .5111116 .0827069 6.18 0.000 .349009 .6732142
6 0 | .2186457 .0674192 3.24 0.001 .0865065 .3507849
6 1 | .6185467 .0611384 10.12 0.000 .4987176 .7383758
7 0 | .4044675 .0706157 5.73 0.000 .2660634 .5428717
7 1 | .7155135 .0476424 15.02 0.000 .6221361 .8088909
8 0 | .622414 .0675171 9.22 0.000 .4900828 .7547452
8 1 | .795962 .0437198 18.21 0.000 .7102727 .8816513
9 0 | .8000327 .0610245 13.11 0.000 .6804269 .9196385
9 1 | .8581703 .0421992 20.34 0.000 .7754613 .9408793
10 0 | .9066322 .0443791 20.43 0.000 .8196508 .9936136
10 1 | .9037067 .0387213 23.34 0.000 .8278143 .9795991
11 0 | .9592963 .0267857 35.81 0.000 .9067973 1.011795
11 1 | .9357181 .0332641 28.13 0.000 .8705218 1.000914
------------------------------------------------------------------------------
In total, there are 22 values in the r(b) matrix. Once again we will put these values into a matrix called b.
matrix b=r(b)'
matrix list b
b[22,1]
r1
1._at#0.f .00331151
1._at#1.f .15299932
2._at#0.f .00799952
2._at#1.f .21885736
3._at#0.f .01919637
3._at#1.f .30292513
4._at#0.f .04534893
4._at#1.f .4026401
5._at#0.f .10337563
5._at#1.f .51111162
6._at#0.f .21864569
6._at#1.f .6185467
7._at#0.f .40446752
7._at#1.f .71551351
8._at#0.f .622414
8._at#1.f .79596199
9._at#0.f .8000327
9._at#1.f .85817031
10._at#0.f .9066322
10._at#1.f .9037067
11._at#0.f .95929634
11._at#1.f .93571812
Next, we need the 11 values of read found in the r(at) matrix. However, r(at) is an 11×3 matrix for which we only need the column labeled s. We also need to have two of each of the values of read which we will get by again using the Kronecker product.
matrix at=r(at)
matrix list at
at[11,3]
0b. 1.
f f s
r1 . . 20
r2 . . 25
r3 . . 30
r4 . . 35
r5 . . 40
r6 . . 45
r7 . . 50
r8 . . 55
r9 . . 60
r10 . . 65
r11 . . 70
matrix at=at[1...,"s"]#(11) /* # is the kronecker operator */
matrix list at
at[22,1]
s:
c1
r1:r1 20
r1:r2 20
r2:r3 25
r2:r4 25
r3:r5 30
r3:r6 30
r4:r7 35
r4:r8 35
r5:r9 40
r5:r10 40
r6:r11 45
r6:r12 45
r7:r13 50
r7:r14 50
r8:r15 55
r8:r16 55
r9:r17 60
r9:r18 60
r10:r19 65
r10:r20 65
r11:r21 70
r11:r22 70
We will put the at and b matrices together into matrix p and save the values to data. In this case it is easier to create the alternating series of 0’s and 1’s once the data are in memory using the line number modulo 2 with the not prefix so that the series starts with zero.
matrix p=at,b
matrix list p
p[22,2]
s:
c1 r1
r1:r1 20 .00331151
r1:r2 20 .15299932
r2:r3 25 .00799952
r2:r4 25 .21885736
r3:r5 30 .01919637
r3:r6 30 .30292513
r4:r7 35 .04534893
r4:r8 35 .4026401
r5:r9 40 .10337563
r5:r10 40 .51111162
r6:r11 45 .21864569
r6:r12 45 .6185467
r7:r13 50 .40446752
r7:r14 50 .71551351
r8:r15 55 .622414
r8:r16 55 .79596199
r9:r17 60 .8000327
r9:r18 60 .85817031
r10:r19 65 .9066322
r10:r20 65 .9037067
r11:r21 70 .95929634
r11:r22 70 .93571812
svmat p, names(p)
generate p3 = ~mod(_n,2) in 1/22
clist p1-p3 in 1/23
p1 p2 p3
1. 20 .0033115 0
2. 20 .1529993 1
3. 25 .0079995 0
4. 25 .2188574 1
5. 30 .0191964 0
6. 30 .3029251 1
7. 35 .0453489 0
8. 35 .4026401 1
9. 40 .1033756 0
10. 40 .5111116 1
11. 45 .2186457 0
12. 45 .6185467 1
13. 50 .4044675 0
14. 50 .7155135 1
15. 55 .622414 0
16. 55 .795962 1
17. 60 .8000327 0
18. 60 .8581703 1
19. 65 .9066322 0
20. 65 .9037067 1
21. 70 .9592963 0
22. 70 .9357181 1
23. . . .
Now we can go ahead and graph the probabilities for f equal zero and f equal one using the graph twoway command.
graph twoway (line p2 p1 if p3==0)(line p2 p1 if p3==1), ///
legend(order(1 "f=0" 2 "f=1")) xtitle(variable s) ///
ytitle(probability) name(probability, replace)
The graph of the probabilities above is nice as far as it goes but the presentation of the results might be clearer if we were to graph the difference in probabilities between f equal zero and f equal one. To do this we will need to rerun the margins command computing the discrete change for f at each value of read. We can get the difference using the dydx (derivative) option.
margins, dydx(f) at(s=(20(5)70)) vsquish
Conditional marginal effects Number of obs = 200
Model VCE : OIM
Expression : Pr(y), predict()
dy/dx w.r.t. : 1.f
1._at : s = 20
2._at : s = 25
3._at : s = 30
4._at : s = 35
5._at : s = 40
6._at : s = 45
7._at : s = 50
8._at : s = 55
9._at : s = 60
10._at : s = 65
11._at : s = 70
------------------------------------------------------------------------------
| Delta-method
| dy/dx Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
1.f |
_at |
1 | .1496878 .0987508 1.52 0.130 -.0438602 .3432358
2 | .2108578 .1107226 1.90 0.057 -.0061545 .4278701
3 | .2837288 .1138312 2.49 0.013 .0606236 .5068339
4 | .3572912 .107035 3.34 0.001 .1475064 .567076
5 | .407736 .0966865 4.22 0.000 .2182339 .597238
6 | .399901 .0910124 4.39 0.000 .22152 .578282
7 | .311046 .0851843 3.65 0.000 .1440878 .4780042
8 | .173548 .0804362 2.16 0.031 .0158959 .3312001
9 | .0581376 .0741941 0.78 0.433 -.0872801 .2035553
10 | -.0029255 .0588969 -0.05 0.960 -.1183612 .1125102
11 | -.0235782 .042708 -0.55 0.581 -.1072843 .0601279
------------------------------------------------------------------------------
Note: dy/dx for factor levels is the discrete change from the base level.
Once again we will need the r(b) and r(at) values. But, we need to be careful this time because of the way that margins parameterizes these estimates there are 11 zero’s in the r(b) matrix before the 11 point estimates of interest. We will need to drop all of the rows with zero values.
matrix b=r(b)'
matrix list b
b[22,1]
r1
0b.f:1._at 0
0b.f:2._at 0
0b.f:3._at 0
0b.f:4._at 0
0b.f:5._at 0
0b.f:6._at 0
0b.f:7._at 0
0b.f:8._at 0
0b.f:9._at 0
0b.f:10._at 0
0b.f:11._at 0
1.f:1._at .14968781
1.f:2._at .21085784
1.f:3._at .28372876
1.f:4._at .35729117
1.f:5._at .40773598
1.f:6._at .399901
1.f:7._at .31104599
1.f:8._at .17354799
1.f:9._at .05813761
1.f:10._at -.00292551
1.f:11._at -.02357822
matrix b=b[12...,1]
matrix list b
b[11,1]
r1
1.f:1._at .14968781
1.f:2._at .21085784
1.f:3._at .28372876
1.f:4._at .35729117
1.f:5._at .40773598
1.f:6._at .399901
1.f:7._at .31104599
1.f:8._at .17354799
1.f:9._at .05813761
1.f:10._at -.00292551
1.f:11._at -.02357822
matrix at=r(at)
matrix at=at[1...,"s"]
matrix list at
at[11,1]
s
r1 20
r2 25
r3 30
r4 35
r5 40
r6 45
r7 50
r8 55
r9 60
r10 65
r11 70
For this graph of the differences in probabilities we will also include the lines for the 95% confidence interval. We will need to compute the upper and lower confidence values manually. To do this we will need the standard errors from the margins command. We will begin by getting the variance-covariance matrix of the estimates. Once again we will have to drop all of the rows and columns that are all zeros.
matrix v=r(V)
matrix v=v[12...,12...]
matrix list v
symmetric v[11,11]
1.f: 1.f: 1.f: 1.f: 1.f: 1.f: 1.f: 1.f:
1. 2. 3. 4. 5. 6. 7. 8.
_at _at _at _at _at _at _at _at
1.f:1._at .00975172
1.f:2._at .01090899 .01225949
1.f:3._at .01106858 .01252966 .01295755
1.f:4._at .00988489 .01133267 .01196523 .0114565
1.f:5._at .00746032 .00876206 .00961621 .00982147 .00934828
1.f:6._at .00437948 .00540815 .00638981 .00727808 .00804735 .00828325
1.f:7._at .00138994 .00201386 .00284747 .0039552 .0053869 .00685686 .00725637
1.f:8._at -.00090133 -.00072541 -.00028361 .00050889 .00177443 .00364155 .00568431 .00646999
1.f:9._at -.00219008 -.00231219 -.00218045 -.00173163 -.00085595 .00078305 .0033247 .00541876
1.f:10._at -.00260473 -.00282965 -.00279956 -.00246344 -.00177185 -.00047237 .00167411 .00372089
1.f:11._at -.00249143 -.00271383 -.00268483 -.00235557 -.00172561 -.00069732 .00085579 .00235456
1.f: 1.f: 1.f:
9. 10. 11.
_at _at _at
1.f:9._at .00550476
1.f:10._at .00422715 .00346884
1.f:11._at .00284913 .00245636 .00182397
The standard errors are found by taking the square roots of the digonal values of v. This actually invloves a series of steps; getting the diagonal values (vecdiag), putting them into a diagonal matrix (diag), getting the square roots (cholesky) and finally putting the values into a column vector (vecdiag). The one step that may not have been obvious is using the Cholesky decompostion on a diagonal matrix which gives us the square roots of the diagonal values. These four matrix functions are included in a single command shown below.
matrix se=vecdiag(cholesky(diag(vecdiag(v))))'
matrix list se
se[11,1]
r1
1.f:1._at .09875081
1.f:2._at .1107226
1.f:3._at .11383123
1.f:4._at .10703502
1.f:5._at .0966865
1.f:6._at .09101237
1.f:7._at .08518432
1.f:8._at .08043624
1.f:9._at .07419408
1.f:10._at .05889686
1.f:11._at .04270799
Next, we will put the at, b and se matrices together into matrix d and save them to the data in memory. Now we can generate the values for the upper and lower 95% confidence intervals.
matrix d=at,b,se
matrix list d
d[11,3]
s r1 r1
r1 20 .14968781 .09875081
r2 25 .21085784 .1107226
r3 30 .28372876 .11383123
r4 35 .35729117 .10703502
r5 40 .40773598 .0966865
r6 45 .399901 .09101237
r7 50 .31104599 .08518432
r8 55 .17354799 .08043624
r9 60 .05813761 .07419408
r10 65 -.00292551 .05889686
r11 70 -.02357822 .04270799
svmat d, names(d)
generate ul = d2 + 1.96*d3
generate ll = d2 - 1.96*d3
clist d1-ll in 1/12
d1 d2 d3 ul ll
1. 20 .1496878 .0987508 .3432394 -.0438638
2. 25 .2108578 .1107226 .4278741 -.0061585
3. 30 .2837287 .1138312 .506838 .0606195
4. 35 .3572912 .107035 .5670798 .1475025
5. 40 .407736 .0966865 .5972415 .2182304
6. 45 .399901 .0910124 .5782852 .2215168
7. 50 .311046 .0851843 .4780073 .1440847
8. 55 .173548 .0804362 .331203 .015893
9. 60 .0581376 .0741941 .203558 -.0872828
10. 65 -.0029255 .0588969 .1125123 -.1183634
11. 70 -.0235782 .042708 .0601294 -.1072859
12. . . . . .
Here is the graph twoway command to plot the differeces in probability with lines for the upper and lower confidence intervals.
graph twoway (line d2 ul ll d1), yline(0) legend(off) ///
xtitle(variable s) ytitle(difference in probability) ///
name(difference1, replace)
As nice as the above graph is, it might look better done as a range plot with area shading between the upper and lower confidence bounds.
graph twoway (rarea ul ll d1, color(gs13) lcolor(gs13)) ///
(line d2 d1), yline(0) legend(off) ///
xtitle(variable s) ytitle(difference in probability) ///
name(difference2, replace)
If you want the lines in these graphs to be smoother, just include more values in the at option, say (20(2)70) instead of (20(5)70).
Below you will find all of the code from both examples without the commentary.
/* example 1 -- plot adjusted cell means */
use https://stats.idre.ucla.edu/stat/data/hsbdemo, clear
anova write prog##female math read
margins prog#female, asbalanced
mat b=r(b)'
mat prg=(123)#(11)
mat fem=(111)#(01)
mat c=prg,fem,b
svmat c, names(c)
/* plot prog by female */
graph twoway (connect c3 c1 if c2==0)(connect c3 c1 if c2==1), ///
xlabel(1 2 3) legend(order(1 "female=0" 2 "female=1")) ///
xtitle(prog) ytitle(ajdusted cell means) ///
name(prog_by_female, replace)
/* plot female by prog */
graph twoway (connect c3 c2 if c1==1)(connect c3 c2 if c1==2)(connect c3 c2 if c1==3), ///
xlabel(0 1) legend(order(1 "prog=1" 2 "prog=2" 3 "prog=3")) ///
xtitle(female) ytitle(ajdusted cell means) ///
name(female_by_prog, replace)
/* example 2 -- categorical by continuous */
use https://stats.idre.ucla.edu/stat/data/logitcatcon, clear
logit y i.f##c.s, nolog
/* plot probabilities by categorical variable */
margins f, at(s=(20(5)70)) vsquish
mat b=r(b)'
mat at=r(at)
mat at=at[1...,3]#(11)
mat p=at,b
svmat p, names(p)
generate p3 = ~mod(_n,2) in 1/52
graph twoway (line p2 p1 if p3==0)(line p2 p1 if p3==1), ///
legend(order(1 "f=0" 2 "f=1")) xtitle(variable s) ///
ytitle(probability) name(probability, replace)
/* difference in probabilities */
margins, dydx(f) at(s=(20(5)70)) vsquish
mat b=r(b)'
mat at=r(at)
mat at=at[1...,3]
local nrow=rowsof(at)
mat b=b[`nrow'+1...,1]
mat se=r(V)
mat se=se[`nrow'+1...,`nrow'+1...]
mat se=vecdiag(cholesky(diag(vecdiag(se))))'
mat d=at,b,se
mat lis d
svmat d, names(d)
gen ul = d2 + 1.96*d3
gen ll = d2 - 1.96*d3
graph twoway (line d2 ul ll d1), yline(0) legend(off) ///
xtitle(variable s) ytitle(difference in probability) ///
name(difference1, replace)
graph twoway (rarea ul ll d1, color(gs13) lcolor(gs13)) ///
(line d2 d1), yline(0) legend(off) ///
xtitle(variable s) ytitle(difference in probability) ///
name(difference2, replace)