Dear all,

I am working with cross sectional survey data from 22 countries with 16965 observations. My DVs are factor scores of latent constructs created with the SEM builder, each consisting of responses to 3 Likert style items. The model includes an interaction of an individual level continuous predictor (also a factor score, constructed from 7 Likert style items) with a country level ordered categorical variable that ranges from 1-3. The categorical variable groups observations according to their respective country's score on a metric of interest. The model includes 9 additional individual level covariates and 2 country level covariates. I am trying to evaluate the extent to which the relationships between my individual level predictors and outcomes vary across each of the 3 country categories. I had 4-6% missing data on my individual level outcomes and main predictors, and ran multiple imputation using chained equations. Following the advice of (Tilling, Kate, Williamson, Elizabeth J, Spratt, Michael, Sterne, Jonathan A.C, & Carpenter, James R. (2016). Appropriate inclusion of interactions was needed to avoid bias in multiple imputation. Journal of Clinical Epidemiology, 80, 107-115.), I imputed each group delineated by the categorical interaction term separately to avoid bias. The discussion of this approach is found on page 113.

So far I have been able to get Stata to produce and graph average marginal effects of my continuous predictor at each group level of my categorical predictor, which is helpful and does begin to tell the story the analysis seems to be indicating. The graph consists of three plots, one for each group. However, what I really want to produce are plots/graphs with predicted values of my DV on the y-axis, the range of values of my continuous predictor on the x-axis (e.g. at percentiles), with a separate line for each group. Essentially I want to do what the user generated command - marginscontplot (Royston, P. (2013). Marginscontplot: Plotting the Marginal Effects of Continuous Predictors. The Stata Journal, 13(3), 510-527.) - is able to do. However that command calls the margins command, and since I am running my analysis through the multiple imputation framework, I cannot use the margins command, and I am instead using the user generated command, mimrgns (ssc).

I've read the Stata manuals and help files for margins and predict several times, read the help file for mimrgns, and the above referenced article introducing marginscontplot, and I cannot put together an approach that works. My understanding is that at least part of the problem is related to Stata's storing of estimation results and the ways mimrgns and marginscontplot are designed to access those estimation results. But I'll leave that to the experts. My dataex(ssc) pull and code is below. This is my first post, and I really tried to do this correctly, however I apologize if this is too wordy, too vague, or difficult to read. Thank you in advance.

Code:
* Example generated by -dataex-. To install: ssc install dataex
clear
input float(PJ_c LEG_c MOR_c PL PC_YP I_status_1st I_status_2nd income01 Depravation family01 Slf_ctrl_c Age_c gender_M HDI_IA HOMICIDE) int country
          .          .  .22735445 3 1 0 1 0 4 1  -58.80923 -.7525997 1  .87  .6  310
 -.00998439   .2352361  .03168412 1 0 0 0 1 4 1   15.26484  .2474003 0 .685 1.2 3810
 -.23670854 -.13152407  .21952204 3 0 0 0 1 7 1   -21.7722  .2474003 1 .876 1.6 3580
 -.23670854 -1.2006817  -.6532929 2 0 0 0 1 6 0   -21.7722 -.7525997 0  .85  .6 4200
  .51026696  .57087344  .22343825 3 0 0 1 0 3 1   -5.10553 1.2474003 0 .861  .9  490
 -.00998439   .2352361 -.05408466 1 0 0 0 1 4 1  -6.957382  .2474003 0 .685 1.2 3810
          .          .  .07350265 3 1 0 1 0 6 0 -14.364792  .2474003 1  .87  .6  310
-.035236895   .2352361   .1596446 3 0 0 1 1 2 1  4.1537285  .2474003 1 .861  .9  490
   .5355195  .57087344   .1596446 3 0 0 0 1 5 1 -14.364792  .2474003 0 .876 1.6 3580
  .33404785  -.9315699  -.8331251 2 0 . . . . . -29.179605         . .  .85  .6 4200
 -.23670854 -.13152407   .1635608 3 0 0 0 0 2 0  .45002365 1.2474003 0 .861  .9  490
          .          .  -.4972283 3 0 0 1 0 2 0   -25.4759  .2474003 1  .87  .6  310
 -.23670854 -.13152407  .07105988 2 0 0 0 1 4 1   15.26484 1.2474003 1  .85  .6 4200
   .7117386   .9376336  .07105988 3 0 0 1 1 6 1  18.968544  .2474003 1 .876 1.6 3580
 -.43818015 -.22917242  .21952204 1 0 0 1 1 4 1 -3.2536774  .2474003 0 .685 1.2 3810
 -.23670854  .10068018 -.08279192 1 0 0 0 1 4 1  4.1537285  .2474003 0 .685 1.2 3810
          .          .  .22343825 3 0 0 0 1 6 1          . 1.2474003 1 .876 1.6 3580
  -.9836839 -1.0719105  .22735445 2 0 0 0 1 6 0 -10.661087  .2474003 1  .85  .6 4200
  .14098223   .9376336  .22735445 3 0 0 1 0 4 1  .45002365  .2474003 0 .861  .9  490
          .          . -.26474053 3 1 0 1 1 4 1  -32.88331  .2474003 0  .87  .6  310
end
label values PC_YP polcolyp
label def polcolyp 0 "no", modify
label def polcolyp 1 "yes", modify
label values Depravation revdeprfam
label def revdeprfam 2 "worse off", modify
label def revdeprfam 3 "somewhat worse off", modify
label def revdeprfam 4 "the same", modify
label def revdeprfam 5 "somewhat better off", modify
label def revdeprfam 6 "better off", modify
label def revdeprfam 7 "much better off", modify
label values gender_M male
label def male 0 "female", modify
label def male 1 "male", modify
label values country country
label def country 310 "Netherlands", modify
label def country 490 "Germany", modify
label def country 3580 "Finland", modify
label def country 3810 "Rep. Serbia", modify
label def country 4200 "Czech Rep.", modify

Code:
mi impute chained (regress) PJ_c LEG_c MOR_c Slf_ctrl_c (pmm, knn(5)) Age_c (logit) income01 PC_YP family01 gender_M P_discr1 (ologit) DUTY Depravation I_status P_corrupt P_speed = country HDI_c HOMICIDE_c ROL_c, add(20) rseed(8987) by(PL) [pweight = p_weight]

mi estimate : mixed LEG_c PL##c.MOR_c PC_YP I_status_1st I_status_2nd income01 Depravation family01 Slf_ctrl_c Age_c gender_M HDI_c HOMICIDE_c [pweight = p_weight] || country:

mimrgns, dydx(c.MOR_c) at(PL=(1 (1) 3)) atmeans cmdmargins
 
marginsplot