Adjusted predictions at representative values (APRs) for male/female

[Julian Pritsch]

I estimated an multilevel ordinal model and calculated the APRs (Adjusted predictions at representative values) for male and female separately. Surprisingly, the resulting margins are exactly the same for both sexes. This seems odd to me since the main effect of sex is highly significant and the strongest effect in my results.

To give you an idea about the (theoretical) model:

My DV is fear of crime and has for categories (0-3) and the IV of interest is standardized and therefore continuous. I calculated the margins for the last two outcomes (2 and 3) with the following command (the code was a bit longer since I saved the model and resulting graphs, but to keep it more simple I only show you the important code):

• armut=poverty on level 2
• sex=gender

forvalues i=2/3 {
            marginscontplot armut sex, at1(-2(0.25)2) ///
            margopts(predict(outcome(`i')fixedonly)
}

One of the resulting graphs (outcome 2) looks like this.




It seems somewhat unlikely to me that the margins for male and female are exactly the same?

What do you think?

[Richard Williams]

Why don't you show the code (and preferably the output too) from the estimation command.

Also, indicate where user-written commands come from. I am assuming you got marginscontplot (aka mcp) from Stata Journal. If not, please say where you did get it from, just in case there is some other version of it floating around.

I find that using the showmarginscmd option of mcp can help make it clearer what is being done.

[Julian Pritsch]

Hi Richard,

thanks for your answer.

You are right, I got the marginscontplot-comand from the Stata Journal. And I used the showmarginscmd option which resulted in the following output (for outcome(2)):

margins , at( armut=( -2 -1.75 -1.5 -1.25 -1 -.75 -.5 -.25 0 .25 .5 .75 1 1.25 1
> .5 1.75 2) sex=( 0 1)) predict(outcome(2)fixedonly)

Is that the output you meant? Or did you mean the code/output from the multilevel estimation?

[Richard Williams]

Show the code/output from the estimation command. I am not sure what you used - xtologit or something else?

[Julian Pritsch]

The command I used for the multilevel estimation is meologit and the output as follows:

meologit foc ib0.sex agez c.agez#c.agez ib5.bildungr ib1.erwerbr01 ib4.einkr ///
ib0.kv01 ib0.diebst01 ib0.einbr01 ib0.raub01 ib0.betrug01 ///
osten armut auslr_sqz urban c.urban#c.urban ///
|| plz: ,or vce(robust) intpoints(30) baselevel level(99)

Mixed-effects ologit regression                 Number of obs      =     31915
Group variable:             plz                 Number of groups   =      6330

                                                Obs per group: min =         1
                                                               avg =       5,0
                                                               max =        44

Integration method: mvaghermite                 Integration points =        30

                                                Wald chi2(22)      =   4942,97
Log pseudolikelihood = -31962,698               Prob > chi2        =    0,0000
                                     (Std. Err. adjusted for clustering on plz)
-------------------------------------------------------------------------------
              |               Robust
          foc | Odds Ratio   Std. Err.      z    P>|z|     [99% Conf. Interval]
--------------+----------------------------------------------------------------
          sex |
        mann  |          1  (base)
        frau  |   3,165227   ,0737817    49,43   0,000     2,980771    3,361097
              |
         agez |   1,193623   ,0165582    12,76   0,000     1,151725    1,237045
              |
c.agez#c.agez |   1,193624   ,0159851    13,22   0,000     1,153151    1,235517
              |
     bildungr |
noch schüler  |   1,639009   ,1435897     5,64   0,000     1,307909    2,053928
bis haupts..  |    1,75933    ,062664    15,86   0,000     1,605102    1,928378
realschul..)  |    1,44066   ,0428039    12,29   0,000     1,334518    1,555245
(fach-)hoc~e  |   1,168426   ,0422922     4,30   0,000     1,064413    1,282604
(fach-)hoc~s  |          1  (base)
              |
    erwerbr01 |
nicht-erwe~g  |    1,13348   ,0330882     4,29   0,000     1,051376    1,221996
erwerbstätig  |          1  (base)
              |
        einkr |
       <1000  |   1,685002   ,0817808    10,75   0,000     1,486984    1,909389
   1000-2000  |   1,424954    ,045934    10,99   0,000     1,311415    1,548323
   2000-3000  |   1,236879    ,037877     6,94   0,000     1,143063    1,338395
       >3000  |          1  (base)
keine angabe  |   1,441261   ,0672731     7,83   0,000     1,277989    1,625393
              |
         kv01 |
    kein mal  |          1  (base)
       1-mal  |   1,053738    ,084454     0,65   0,514     ,8571851     1,29536
              |
     diebst01 |
    kein mal  |          1  (base)
       1-mal  |   1,064506   ,0667287     1,00   0,319     ,9057828    1,251042
              |
      einbr01 |
    kein mal  |          1  (base)
       1-mal  |   1,983802   ,2082433     6,53   0,000     1,513804    2,599723
              |
       raub01 |
    kein mal  |          1  (base)
       1-mal  |   1,197597   ,1761374     1,23   0,220     ,8199392      1,7492
              |
     betrug01 |
    kein mal  |          1  (base)
       1-mal  |   1,284136   ,0626085     5,13   0,000     1,132583    1,455969
              |
        osten |   1,379838   ,0577357     7,69   0,000     1,238855    1,536866
        armut |   1,286003   ,0207945    15,56   0,000      1,23354    1,340697
    auslr_sqz |   1,284451   ,0258733    12,43   0,000     1,219505    1,352856
        urban |   1,223151   ,0304328     8,10   0,000      1,14722    1,304107
              |
      c.urban#|
      c.urban |   ,9302859   ,0047761   -14,08   0,000     ,9180645    ,9426699
--------------+----------------------------------------------------------------
        /cut1 |   1,162065   ,0319111    36,42   0,000     1,079868    1,244263
        /cut2 |   3,610915   ,0380961    94,78   0,000     3,512786    3,709044
        /cut3 |   5,327674   ,0475775   111,98   0,000     5,205122    5,450226
--------------+----------------------------------------------------------------
plz           |
    var(_cons)|   ,1277647   ,0151887                       ,094064    ,1735395
-------------------------------------------------------------------------------

[Richard Williams]

Thanks. Now I am wondering if you have left out something important from your intermediate steps. I adapted this from the manual:

webuse tvsfpors, clear
meologit thk prethk cc##tv || school: || class:
mcp prethk cc, margopts(predict(outcome(2)fixedonly)) show

I only get one graph with both groups together (instead of two separate graphs) and I don't get the confidence intervals included. So, I think that whatever command created your graphs is not the same as the command you show.

Why don't you run your meologit command, followed by

mcp armut sex, at1(-2(0.25)2) margopts(predict(outcome(2)fixedonly) show

Maybe you have done this, but I just want to make sure that an error didn't sneak in at some intermediate step that you aren't showing us.

Also, are you sure about the range you specified for armut? If, say, it really runs from -200 to 200, your graphs could be misleading.

[Julian Pritsch]

(1) Sorry for the confusion:

Like I wrote in my first post, I intentionally left out some lines of code to keep it simpler. Unfortunately, one of those lines included the ci-option of the command. Sorry for that.

Therefore, if you use the ci-option you should get two separate graphs like P. Roysten wrote in the Stata Journal (2013):13 Number 3, pp.520.

(2)
The "armut" variable is an z-standardized variable and I want to show the range of -2 standard deviations / +2 standard deviations.That's why I chose these values.

[Richard Williams]

I am not sure what to tell you. Tweaking my earlier code,

webuse tvsfpors, clear
meologit thk prethk i.cc i.tv || school: || class:
mcp prethk cc, margopts(predict(outcome(2)fixedonly)) show ci

The two graphs are certainly not identical.

Maybe you should specify a broader range, e.g. -4 to 4? Or run the margins command you gave earlier and look at the actual numbers? Maybe the differences are small and almost imperceptible within the range you used, but looking at the actual numbers across a broader range would reveal at least slight differences. And again, just double-check that armut really is coded the way you think it is coded, e.g. run summary statistics on it. More than once I have seen people (including me) who weren't using the variables they thought they were using,

[Julian Pritsch]

Thanks for your time anyway Richard! I will keep on trying to find out where the problem is and will post the solution (hopefully!) when I find it.

[Julian Pritsch]

Dear Richard,
dear all others who follow this threat,

I found a "solution" which is rather tedious. I estimated the multilevel models separately for male and female. After that I estimated the margins with aforementioned mcp command for outcome(2) and outcome(3) and combined the graphs.

One of the resulting graphs look like this (without going into detail here regarding its content):



In my opinion, one small problem with this approach is, that I compare two slightly different models (subgroups for male and female) in a non-linear frame which is problematic.

[Richard Williams]

By running separate models for each gender, you allow all parameters to differ by gender, whereas in the original model you only allowed the intercepts to differ. This may be ok to do, but it is a pretty major change from the original.

I don't think you should have to do it this way. I still wonder, with your original model, could the differences be small and almost imperceptible? Could mcp be making an error with the graphs? What happens if you run your original meologit command followed by

margins sex , at( armut=( -2 -1.75 -1.5 -1.25 -1 -.75 -.5 -.25 0 .25 .5 .75 1 1.25 1.5 1.75 2) ) predict(outcome(2) fixedonly)
marginsplot, noci
margins r.sex , at( armut=( -2 -1.75 -1.5 -1.25 -1 -.75 -.5 -.25 0 .25 .5 .75 1 1.25 1.5 1.75 2) ) predict(outcome(2) fixedonly)
marginsplot

This should tell you whether the estimated predictions differ by gender. Perhaps the difference is in the third decimal place, which is why the original graphs appear identical. If the results for gender seem very different, then maybe mcp broke down somewhere.

[Julian Pritsch]

Your commands worked much better. Although it still surprises me, that they won't work with the mcp command.

Here are the resulting graphs:





So it seems to me that there is a difference between male and female just as theory and empirical evidence suggests.

[Richard Williams]

That is not only what theory suggests, it is what your estimation results said should happen, since the effect of gender was highly significant.

This makes me wonder if there is a bug in mcp. But, I can't reproduce the problem with other data, e.g. this seems to work fine:

webuse tvsfpors, clear
meologit thk prethk i.cc i.tv || school: || class:
mcp prethk cc, margopts(predict(outcome(2)fixedonly)) show 
mcp prethk cc, margopts(predict(outcome(2)fixedonly)) show ci

Perhaps you want to uninstall and reinstall mcp; and also check to make sure that your version of Stata is totally up to date.

You may also want to look at Nick Winter's combomarginsplot command, available from SSC. It can be quite useful after multiple outcome commands like ologit.