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, clear

sum female ses read math

    Variable |       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 programwill 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 miereturn 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 math   


Multiple-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 above  3.   mat b= e(b_mi)                        // save mi point estimates  4.   mat V = e(V_mi)                       // save mi vce  5.  . * run ologit and margins on the _mi_m==0 data.   quietly ologit ses female read math if _mi_m == 0  6.   quietly margins, at(female=(0 1) read=(30(10)70)) ///>     atmeans asbalanced predict(outcome(`i'))  7. .   myret                                 // myret is defined above  8.         *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 margins  9. . 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.