How to use coefplot margins and marginsplot in Stata.

What you see up there is (obviously) a combination of two graphs made into one, which is a thing you can do in Stata. See below.

The best way to learn this part is to just do it. I’m going to share the code that I used and you should try it out with your models and your variables. Make some interaction plots, which can be a helpful visual. Don’t limit yourself to the two that I did here. These are just examples, right?

Here is the code I used to make Figure 1. (Remember to run initial-erika-setup.do file first.) And if you are unsure of something with a built-in Stata function such as marginsplot, just type help marginsplot in the Stata command line to access the robust help documentation. My examples are using margins on a mlogit, with standardized continuous variables, and some styling is set with grstyle. This first chunk of code is long because I’m setting up your starter grstyle code, which is pretty cool and I think you will like it once you learn it. I suggest copying this into a do-file and then editing it with your models, your variable names, etc. See how I labeled my mlogit under the MLOGIT header in my code? I did that so you can hopefully readily see where to substitute your own things.

/*********************************************************
                    ACKNOWLEDGEMENTS
**********************************************************/
/*My dataset in this demo is Gallup World Poll data, which I cannot share. The post-estimation Spost commands used are courtesy of J. 
Scott Long. Many other commands are courtesy of 
Ben Jann, including coefplot, center, and grstyle.

References:
     Ben Jann, 2013. "COEFPLOT: Stata module to plot regression 
coefficients  and other results," Statistical Software Components 
S457686, Boston College Department of Economics, revised 17 Dec 2020. http://repec.sowi.unibe.ch/stata/coefplot/index.html
     J. Scott Long & Jeremy Freese. 2014. "Regression Models for 
Categorical Dependent Variables Using Stata", 3rd Edition. College 
Station, TX: Stata Press. https://jslsoc.sitehost.iu.edu/spost13.htm
*/

/*********************************************************
                        MLOGIT
**********************************************************/

/*
multinomial logit model (MNLM): 
Y = LifeEval (Life Evaluation Index, evaluative well-being)
C = income_2 (household income in international dollars)
F = Corruption (Corruption Index - perceived govt and business corruption)
Xs = gender, urbanicity, age, education, Region
I am estimating the MNLM of Y on C, F, and X. 
*/

//multinominal logit model
svy: mlogit LifeEval c.income_2 i.Corruption i.gender i.urbanicity c.age ///
i.education i.Region, base(2) noci //base(2)=suffering, middle category ref.
eststo mlogitLE
mlogtest
mlogtest, combine //Wald test to make sure this model can distinguish outcomes
/* btw mlogtest is an spost command, and it's testing all of your IVs at once
for you. Check out the output on this...you can add in a log ratio test if you
want by adding a lr */

//multinomial logit model with interaction term 
svy: mlogit LifeEval c.income_2 i.Corruption i.gender i.urbanicity c.age ///
c.income_2#i.Corruption i.education i.Region, base(2)
eststo mlogitLE_with_interact

/**********************************************************
   grstyle - for styling graphs - written by Ben Jann
**********************************************************/

*Graph Settings - set up how I like them. Click thru the docs to set your way.
//Documentation is here: http://repec.sowi.unibe.ch/stata/grstyle/
//YES, you can completely skip grstyle if you want to just use Stata defaults
grstyle clear //resets any previous grstyle in the file
grstyle init //initializes grstyle to get ready to run
grstyle set horizontal //sets y axis tick labels horizontal/readable
grstyle set ci //makes shading of CIs transparent
grstyle set grid //turns on all grid lines, horiz, vert, min and max
grstyle set legend 6 //clock position; add ", inside" to move it inside graph
grstyle color background white
grstyle yesno draw_major_hgrid yes
grstyle yesno draw_major_ygrid yes
grstyle linepattern major_grid dot
grstyle linewidth plineplot medthick
grstyle anglestyle vertical_tick horizontal
grstyle symbolsize p large
grstyle gsize axis_title_gap tiny //adds space between ticks and axis titles
grstyle color major_grid black
grstyle linewidth major_grid thin
grstyle set size small: subheading axis_title
grstyle color p1markline navy%0
grstyle color p1markfill navy%100
grstyle color p2markline maroon%0
grstyle color p2markfill maroon%100
grstyle color p3markline forest_green%0
grstyle color p3markfill forest_green%100
grstyle color p4markline dkorange%0
grstyle color p4markfill dkorange%100
grstyle color p5markline teal%0
grstyle color p5markfill teal%100
grstyle color ciarea navy%100
grstyle color ci2_area maroon%100
grstyle color ci3_area forest_green%100


/****************************************************************
 Interaction Plots - predicted probabilities - using MARGINSPLOT
*****************************************************************/

/*
PREDICTED PROBABILITIES:
Predicted probabilities refer to the predicted probability of one given outcome category, compared to all other outcome categories. Does not matter what base or reference is used. For mlogit, I will use margins to make the predictions and marginsplot to graphically display the predicted probabilities.
*/

//marginal effects of income
estimates restore mlogitLE_with_interact
mchange Corruption income_2 , amount(sd)
margins Corruption, dydx(income_2) vce(unconditional) 
/*btw that vce(unconditional) is there to take out the robust vce implied from 
using pweight in the model and this marginsplot is for the Average Marginal Effects 
(AME) of income w/95%CIs */

summarize income_2 , detail //continuous variable;to see IQR, percentiles, etc.
estimates restore mlogitLE_with_interact
margins, at(income=(280(6200)31280)) // I chose middle 98% range of HH income 
//starts at 280, goes to 31280, in increments of 6200; that's how that works

marginsplot , ///
allsimplelabels noci ///
title("MNLM: Predicted probabilities of well-being, by income") ///
subtitle("(thriving is the best category; suffering is the worst category)") ///
xtitle("1st thru 99th Percentile Annual HH Income (Int'l Dollars)") ///
ytitle("Probability of Evaluative Well-Being") ///
xlabel(280(6200)31280) ///
caption("Figure 1. Predicted Probabilities of each Life Evaluation Index Category by annual household income - full model", size(*0.5)) ///
plot1opts(lpattern(solid) ms(oh) color(purple) ) ///
plot2opts(lpattern(solid) ms(dh) color(blue) ) ///
plot3opts(lpattern(solid) ms(th) color(green) ) ///
legend(order(1 "Thriving" 2 "Struggling" 3 "Suffering") ) ///
name(MNLMincome)
graph export `pgm'-MNLM-marginsplot-1.emf, replace

And I think you should make a version of Figure 2, because then you’ll be able to try out the way I explain how to combine graphs. This is the code to produce Figure 2 as above.

/* plotting corruption by region */      
estimates restore mlogitLE_with_interact
quietly margins Corruption#Region, asbalanced
marginsplot, x(Region) noci ///
title("MNLM: Predictive Margins of Perceived Corruption by Region in India") ///
subtitle("I could put a subtitle here...") ///
legend(order(1 "No Corruption & Thriving" ///
2 "No Corruption & Struggling" ///
3 "No Corruption & Suffering" ///
4 "Some Corruption & Thriving" ///
5 "Some Corruption & Struggling" ///
6 "Some Corruption & Suffering" ///
7 "Max Corruption & Thriving" ///
8 "Max Corruption & Struggling" ///
9 "Max Corruption & Suffering")) ///
xtitle("Regions in India") ///
ytitle("Probability of Evaluative Well-Being") ///
caption("Figure 2. Predicted Probabilities of each Life Evaluation Index Category by perceived corruption level and National Region - full model", size(*0.5)) ///
name(MNLMcorruption)
graph export `pgm'-MNLM-marginsplot-2.emf, replace

And if you made figure 1 and figure 2, in some form, the following code will let you combine them to have one graph output file, which is useful. You should also see how to modify this for other use cases.

/*To combine marginsplots together into one graph export, I do it this way. There are other ways - I just like this way so I share it, if you want to use: 
*/

graph combine MNLMincome MNLMcorruption, xsize(6.5) ysize(2.7) ///
iscale(.8) name(MarginsPlotsCombo) ///
title("These Things are Usually Called Interaction Plots") ///
subtitle("Made with marginsplots and love, in Stata...") ///
caption("Note: You can put a caption or figure notation here too...it will go all the way across")
graph close MNLMincome MNLMcorruption
graph export `pgm'-MNLM-marginsplot-combo.emf, replace

okay, okay. Let’s use coefplot now

Okay. So what you see in Figure 3 above is the most basic use for coefplot and that is to plot coefficients. I mean, you know what nonlinear regression output looks like, yeah?

I should make a nice regression table and put it here. Because we absolutely do need to make those too, for the sake of precision, and who doesn’t like seeing little * significance stars ** on our work, right? But if you just want a sense of which coefficients are having the most influence in your model and how? You could use coefplot. Seriously. It’s also really useful for comparing several models, because you can put plots side by side, or on top of one another for that matter, to readily see the similarities and differences between models. Super helpful, and you can make some great looking visuals for your posters and publications.

Here is the code I used to make Figure 3 as above.

/* In statistics, we can standardize things on different scales so that 
we can make comparisons about relative magnitude. 
We're going to do that now using Ben Jann's center command to do it.
It standardizes continuous variables to a mean of zero and SD of 1. Like this:
*/

center income_2 age, standardize prefix(std_)
//then check them (mean is not zero only due to rounding errors)
summarize income_2 std_income_2 age std_age

/*
BTW to get the exact names of coefficients stata is using, so that you can 
rename them like I am doing here, use this:
mlogit, coeflegend
*/

//super basic coefplot example to visualize coefficients of a single model
//running mlogit with standardized variables 
mlogit LifeEval c.std_income_2 i.Corruption i.gender ///
i.urbanicity c.std_age c.std_income_2#i.Corruption i.education ///
i.Region if YEAR_WAVE==2019
eststo std_mlogitfull2019

//coefplot as seen in Figure 3
coefplot (std_mlogitfull2019), ///
drop(_cons) xlabel(, grid) baselevels xline(0) headings(, labsize(vsmall)) ///
ciopts(recast(rcap) color(pink)) citop ///
msymbol(D) mcolor(pink) msize(medium) ///
title("Beta Coefficients in 2019 Model") ///
subtitle("Multinominal Logistic Regression Model of Evaluative Well-Being on Income," ///
"Corruption, and Region, with continuous variables standardized", size(*0.8)) ///
coeflabels(std_income_2="Income (std)" ///
std_age="age in years" ///
O.Corruption="No Corruption" ///
5O.Corruption="Some Corruption" ///
10O.Corruption="Max Corruption" ///
1.gender = "Male" ///
2.gender = "Female" ///
0.Corruption#c.std_income_2 = "No Corruption x Income (std)" ///
50.Corruption#c.std_income_2 = "Some Corruption x Income" ///
100.Corruption#c.std_income_2 = "Max Corruption x Income (std)" ///
1.education = "Up to 8 years edu." ///
2.education = "9-15 years of education" ///
3.education = "College grad or HS+4 edu" ///
1.Region = "Central India" ///
2.Region = "East India" ///
3.Region = "West India" ///
4.Region = "North India" ///
5.Region = "South India" ///
1.urbanicity = "Urban" ///
2.urbanicity = "Rural") ///
caption("Note: I have shown the base levels here for clarity but they are at zero by definition.", size(*0.5)) ///
note("Figure 3. Coefplot demo to look at some coefficients on a single model...", span) ///
name(coefplot1)
graph export `pgm'-MNLM-coefplot-1.emf, replace

What if you do not want all of those coefficients in the plot? You can use drop or keep. Keep commands tend to freak me out (because if you omit something you want to keep, it’s gone, whereas if you forget to drop something you just have an extra thing to drop). I just find it emotionally less riskly to use drop commands, so that’s what you’ll see in the code below, to make Figure 7, which is the same as Figure 3, plus some drops. See below.

//coefplot for figure 7, which is Figure 3 with some coefficients removed
//Note structure of the drop command, and the difference between the above two visuals

estimates restore std_mlogitfull2019
coefplot std_mlogitfull2019, ///
drop(_cons || *.Region || *.urbanicity || 100.Corruption) ///
xlabel(, grid) baselevels xline(0) headings(, labsize(vsmall)) ///
//snip...the rest is the same as for Figure 3