The** margins** and **marginsplot ** commands**,** introduced in Stata 11 and Stata 12, respectively,
are very popular post-estimation commands. However, they can be tricky to use in
conjunction with multiple imputation.

Let’s begin by looking at the data.

use https://stats.idre.ucla.edu/stat/data/hsbmar_s10, clearsum female ses read mathVariable | Obs Mean Std. Dev. Min Max -------------+-------------------------------------------------------- female | 185 .5459459 .4992356 0 1 ses | 200 2.055 .7242914 1 3 read | 185 51.61622 10.19104 28 76 math | 190 52.17895 9.246168 33 75

As you can see from the table above, all of
the variables except for **ses** have missing values.

Running **margins** and **marginsplot** after multiple imputation involves
a multi-step process. We will demonstrate this process using an ordered logit model with
**ses** as the response variable. It can take on the values 1, 2 or 3. It’s not
a great response variable from a theoretical standpoint, but at least it is
ordinal.
The predictor variables are **female**, **read** and **math**.

## So, what’s the problem?

Why not just impute the data and run the analyses. Well, we can impute the data, but we
need a way to run both **ologit** and **margins** on each imputed dataset and
then combine the **margins** results into a single output. The problem
is that **margins (**an rclass command**)** does not work with **mi estimate
(**an eclass command). Additionally, since we are
looking for predicted probabilities, we need to compute them for each of the three
response values.

We can accomplish this by writing a wrapper program called **emargins.ado**.
It contains both the **ologit** and **margins** commands. By setting the option
**properties** to **mi**, **emargins** can be used with **mi estimate**.
We also need to declare **emargins** to be an **eclass** program.

Here is what the program looks like.

program emargins, eclass properties(mi) version 15 args outcome ologit ses female read math margins, at(female=(0 1) read=(30(10)70)) atmeans asbalanced /// post predict(outcome(`outcome')) end

The **emargins **program
will run the **ologit** and then estimate **margins** to give
the predicted probabilities for each
level of ses.

The important part to notice is that the program is marked
“eclass” and we have use the “post”
option on the **margins **statement. This is done so that the
predicted probabilities and variance-covariance matrix estimated for each
imputed dataset will be saved
correctly in the **mi** **ereturn** list where **mi estimate**
can access the estimates (not the **return **list where it would normally go).

Here is how you use **emargins** program with **mi estimate**:

mi estimate, cmdok: emargins 1

The **cmdok** is needed because Stata does not recognize **emargins** as an
**mi estimable**
program. The value one (1) after **emargins** is passed to **margins** indicating which
response value is being predicted.

Once **mi estimate** has combined the **margins** from each of the imputed datasets using Rubin’s
rules into one table, how do we get the **marginsplot** to run? If you try running **marginsplot**
after **mi estimate** you get the error message, “previous command was not margins.” This happens because **mi
estimate** does not leave the results in the right place for **marginsplot** to
find them.

Remember, **mi estimate** is an **eclass **command
saves results in the
**ereturn** list but **margins** is **rclass** and saves
its results in the **return** list. Somehow we need to move the
information from the **mi estimate ereturn** list to the **
margins return** list.

The solution to this problem is to save the combined margins predicted
probabilities **e(b_mi)** and variance-covariance matrix **e(V_mi)**
produced by **mi estimate **into matrices **b** and **V**, run
a standard **margins** on the
**_mi_m == 0** (non-imputed) data, and then repost the results from **b** and **V** back into
the** margins return** list
**r(b)** and **r(V) **where **marginsplot **can access them. We do the last part with a program called **myret.ado** which
looks like this.

program myret, rclass return add return matrix b = b return matrix V= V end

Now putting all of the pieces together into a do-file we get…

set seed 1234543 mi set mlong mi register imputed female math read science socst mi impute mvn female math read science socst = /// ses write awards, add(10) * this is to get the ologit coefficients and standard errors mi estimate: ologit ses female read math * loop once for each of the response values of ses forvalues i=1/3 { mi estimate, cmdok: emargins `i' // emargins is defined above mat b= e(b_mi) // save mi point estimates mat V = e(V_mi) // save mi vce * run ologit and margins on the _mi_m==0 data quietly ologit ses female read math if _mi_m == 0 quietly margins, at(female=(0 1) read=(30(10)70)) /// atmeans asbalanced predict(outcome(`i')) myret // myret is defined above*Technically we ran the program myret between margins and marginsplot.*E(cmd) is the eclass scalar that tells Stata what the previous command was.* So we have to set that to "margins" for marginsplot to work correctly. mata: st_global("e(cmd)", "margins") // set previous cmd to margins marginsplot, x(read) recast(line) noci name(ologit`i', replace) }

Here is what the output looks like when we run the do-file.

. mi register imputed female math read science socst(51 m=0 obs. now marked as incomplete). mi impute mvn female math read science socst = /// > ses write awards, add(10)Performing EM optimization: observed log likelihood = -1814.7997 at iteration 10 Performing MCMC data augmentation ... Multivariate imputation Imputations = 10 Multivariate normal regression added = 10 Imputed: m=1 through m=10 updated = 0 Prior: uniform Iterations = 1000 burn-in = 100 between = 100 ------------------------------------------------------------------ | Observations per m |---------------------------------------------- Variable | Complete Incomplete Imputed | Total -------------------+-----------------------------------+---------- female | 185 15 15 | 200 math | 190 10 10 | 200 read | 185 15 15 | 200 science | 193 7 7 | 200 socst | 188 12 12 | 200 ------------------------------------------------------------------ (complete + incomplete = total; imputed is the minimum across m of the number of filled-in observations.). * this is to get the ologit coefficients and standard errors. mi estimate: ologit ses female read mathMultiple-imputation estimates Imputations = 10 Ordered logistic regression Number of obs = 200 Average RVI = 0.0497 Largest FMI = 0.1225 DF adjustment: Large sample DF: min = 628.16 avg = 10,711.70 max = 25,804.25 Model F test: Equal FMI F( 3, 3581.4) = 6.37 Within VCE type: OIM Prob > F = 0.0003 ------------------------------------------------------------------------------ ses | Coef. Std. Err. t P>|t| [95% Conf. Interval] -------------+---------------------------------------------------------------- female | -.3637883 .2892283 -1.26 0.209 -.9317596 .2041831 read | .0469533 .0179698 2.61 0.009 .0117115 .082195 math | .0189451 .0197378 0.96 0.337 -.0197633 .0576536 -------------+---------------------------------------------------------------- /cut1 | 1.965144 .8645886 .2704931 3.659795 /cut2 | 4.230366 .9113859 2.443998 6.016733 ------------------------------------------------------------------------------. * loop once for each of the response values of ses. forvalues i=1/3 {2.. mi estimate, cmdok: emargins `i' // emargins is defined above3. mat b= e(b_mi) // save mi point estimates4. mat V = e(V_mi) // save mi vce5.. * run ologit and margins on the _mi_m==0 data. quietly ologit ses female read math if _mi_m == 06. quietly margins, at(female=(0 1) read=(30(10)70)) ///> atmeans asbalanced predict(outcome(`i'))7.. myret // myret is defined above8. *Technically we ran the program myret between margins and margins> plot.. *E(cmd) is the eclass scalar that tells Stata what the previous command was.. * So we have to set that to "margins" for marginsplot to work correctly... mata: st_global("e(cmd)", "margins") // set previous cmd to margins9.. marginsplot, x(read) recast(line) noci name(ologit`i', replace)10. }Multiple-imputation estimates Imputations = 10 Adjusted predictions Number of obs = 200 Average RVI = 0.0512 Largest FMI = 0.0912 DF adjustment: Large sample DF: min = 1,122.48 avg = 3,787.74 Within VCE type: Delta-method max = 7,748.29 ------------------------------------------------------------------------------ | Coef. Std. Err. t P>|t| [95% Conf. Interval] -------------+---------------------------------------------------------------- _at | 1 | .3920011 .1073537 3.65 0.000 .1815432 .6024589 2 | .2872842 .0643243 4.47 0.000 .1611907 .4133777 3 | .2013181 .0389331 5.17 0.000 .1249888 .2776474 4 | .1363001 .0339945 4.01 0.000 .0696162 .202984 5 | .0900615 .0341854 2.63 0.009 .0229872 .1571359 6 | .4811309 .1088753 4.42 0.000 .2675365 .6947253 7 | .3671169 .0681464 5.39 0.000 .2334318 .5008019 8 | .2661107 .041778 6.37 0.000 .1841926 .3480288 9 | .1848442 .039817 4.64 0.000 .106792 .2628964 10 | .1243357 .0433637 2.87 0.004 .0393239 .2093476 ------------------------------------------------------------------------------ Variables that uniquely identify margins: female read Multiple-imputation estimates Imputations = 10 Adjusted predictions Number of obs = 200 Average RVI = 0.0325 Largest FMI = 0.0653 DF adjustment: Large sample DF: min = 2,164.73 avg = 234,713.28 Within VCE type: Delta-method max = 1103688.17 ------------------------------------------------------------------------------ | Coef. Std. Err. t P>|t| [95% Conf. Interval] -------------+---------------------------------------------------------------- _at | 1 | .4689071 .0618147 7.59 0.000 .347738 .5900763 2 | .507704 .0394224 12.88 0.000 .4304361 .5849718 3 | .506791 .0382264 13.26 0.000 .4318685 .5817136 4 | .4663628 .0427754 10.90 0.000 .3825192 .5502065 5 | .3968078 .0655825 6.05 0.000 .2681965 .5254191 6 | .4175739 .0727861 5.74 0.000 .2748366 .5603111 7 | .4806358 .044185 10.88 0.000 .3940283 .5672433 8 | .5110769 .0378897 13.49 0.000 .4368143 .5853395 9 | .5008984 .0394922 12.68 0.000 .4234918 .5783051 10 | .4527477 .0577852 7.84 0.000 .3394723 .5660232 ------------------------------------------------------------------------------ Variables that uniquely identify margins: female read Multiple-imputation estimates Imputations = 10 Adjusted predictions Number of obs = 200 Average RVI = 0.0487 Largest FMI = 0.0966 DF adjustment: Large sample DF: min = 1,001.59 avg = 3,854.63 Within VCE type: Delta-method max = 9,547.55 ------------------------------------------------------------------------------ | Coef. Std. Err. t P>|t| [95% Conf. Interval] -------------+---------------------------------------------------------------- _at | 1 | .1390918 .0564432 2.46 0.014 .0284384 .2497452 2 | .2050118 .0526088 3.90 0.000 .1018873 .3081363 3 | .2918909 .0467674 6.24 0.000 .2002166 .3835652 4 | .397337 .0599498 6.63 0.000 .2797618 .5149122 5 | .5131306 .0929592 5.52 0.000 .3307635 .6954977 6 | .1012952 .0431331 2.35 0.019 .0166569 .1859335 7 | .1522473 .0417033 3.65 0.000 .0704114 .2340832 8 | .2228125 .0389796 5.72 0.000 .146346 .2992789 9 | .3142574 .0536114 5.86 0.000 .2091475 .4193672 10 | .4229165 .0904808 4.67 0.000 .245528 .6003051 ------------------------------------------------------------------------------ Variables that uniquely identify margins: female read

Please note: The values in the tables and graphs above are predicted probabilities. The column heading for the margins tables, Coef., is incorrect.

In case you lose track of which values in the **margins** output are which, you
can list the **r(at)** matrix.

matrix list r(at)r(at)[10,3] female read math 1._at 0 30 51.691358 2._at 0 40 51.691358 3._at 0 50 51.691358 4._at 0 60 51.691358 5._at 0 70 51.691358 6._at 1 30 51.691358 7._at 1 40 51.691358 8._at 1 50 51.691358 9._at 1 60 51.691358 10._at 1 70 51.691358

This may not be the most transparent process ever but, in the end, we got the plots
of the predicted probabilities. Of course, the technique shown here is not restricted
to **ologit** but generalizes to many other estimation procedures for use with
**margins** and **marginsplot**.

## Acknowledgements

Thanks to Isabel Cannette of Stata Corp for the suggestion to use **myret** to repost
the **margins** results.