How to plot odds ratios instead of negative probability change after margins, dydx()

[Janine Stubbs]

Hello!

Stata's margin commands make it very easy to plot probabilities after, among other things, logistic regression.

Take for example estimating the log-odds of diabetes in females vs males and how the log-odds might differ in black and non-black ethnicities.

This code generates a graph that could naively be taken to plot a negative probability (impossible), when in fact it plots a change in probability.

webuse nhanes2f, clear
keep if !missing(diabetes, black, sex, age)
su age // mean = 47.6
logistic diabetes i.sex i.black age , nolog
// Calculate the probability of diabetes for male and females
// of African origin at mean age
margins i.sex, at(black = 1 age = 47.6 )
// Interpretation 1: after adjusting for age,
// for people of African origin at mean age,
// the Pr(diabetes) for females (.0628306) is higher than for males (.0543217)
di .0628306 - .0543217 // The difference in probability is .0085089

// estimate the uncertainty around this difference
margins , dydx(sex ) at(age = 47.6 black = 1)
// Interpretation 2: The 95% CI for the difference in probability is -.0017589 to .0187767
// so the difference in probability of diabetes for females may be a decrease of .0017589
// compared to males, to an increase of .0187767 compared to males
// Pr(diabetes for females and males of average age for different race):
margins i.black, dydx(sex ) at(age = 47.6 )
marginsplot, xscale(range(-0.5(1)1.5))

Because of this potential for confusion, I don't particularly like the output.

To be honest, I think I'd rather just plot odds ratios for sex by race. Is there a similar tool as margins to do this?

If not, I'd do linear combinations to get the ORs and 95%CIs of interest, manually type the results into Excel (yuck) and plot it up there, which is not at all satisfying.

Comments and personal opinions on the above welcome and thanks in advance for any tips and tricks!

Janine

[Joseph Coveney]

. . . I'd rather just plot odds ratios for sex by race. Is there a similar tool as margins to do this?

From the help file for margins:

After estimation by logistic, you might specify expression(exp(predict(xb))) to use relative odds rather than probabilities as the response.

So maybe something like the following?

logistic diabetes i.sex##i.black c.age, nolog
margins black, dydx(sex) at((mean) age) expression(exp(predict(xb)))

[Janine Stubbs]

Thanks Joseph,

That's it exactly.

Janine

[Janine Stubbs]

Oops, I think I spoke too soon; that margins code doesn't tally up with the manual calculation afterall. See as another example:

use https://stats.idre.ucla.edu/stat/stata/dae/binary.dta, clear
logit admit c.gre##c.gpa

/* ------------------------------------------------------------------------------
admit | Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
gre | .0174643 .0096062 1.82 0.069 -.0013635 .0362922
gpa | 3.367918 1.722688 1.96 0.051 -.0084885 6.744325
|
c.gre#c.gpa | -.0043266 .0027869 -1.55 0.121 -.0097888 .0011356
|
_cons | -13.8187 5.874827 -2.35 0.019 -25.33315 -2.304252
------------------------------------------------------------------------------ */

********* coefficient for gre at the mean of gpa? (log-linear scale)

// admit = beta_0 + 0.017*gre + 3.37*gpa - 0.004*gre*gpa

su gpa // mean = 3.3899
di .0174643 + (-.0043266 * 3.3899) // 0.00279756 ... on log-linear scale

lincom _b[gre] // 0.174
lincom _b[gre#gpa] // -0.004
lincom _b[gre] + (_b[gre#gpa]*3.3899) // 0.0028 ... on log-linear scale

su gpa
margins, dydx(gre) at(gpa=(`r(mean)')) expression(predict(xb)) // .0027976 ...coefficient for gre

********* OR for gre at the mean of gpa?

di exp(0.0028 ) // 1.003 ... exponentiated: OR for gre (relative odds)

lincom _b[gre] + (_b[gre#gpa]*3.3899), eform // 1.003

su gpa
margins, dydx(gre) at(gpa=( `r(mean)')) expression(exp(predict(xb))) // .0013802

[Andrew Musau]

su gpa
margins, dydx(gre) at(gpa=( `r(mean)')) expression(exp(predict(xb))) // .0013802

Have a look at https://www.statalist.org/forums/for...ng-marginsplot.

For plotting logistic coefficients directly, see coefplot from SSC.

[Janine Stubbs]

Thanks Musau,

True, coefplot plots OR from logistic (or coefficients from logit ) .

However, I'd like either coefplot or margins, dydx() to calculate and plot linear combinations of the regression outputs.

In #4 I show how to manually calculate a simple linear combination with lincom, (which, if you have many combinations, can be graphed using coefplot after assembling a matrix).

Happily, margins, dydx()does correctly calculate the desired coefficient, but I'd like it output the OR; the exponentiated coefficient.

Unfortunately, this is not achieved using (from #4)
margins, dydx(gre) at(gpa=( `r(mean)')) expression(exp(predict(xb))) // .0013802

From your helpful link in #5 to a previous post on this exact problem, I see that it is an outstanding problem since at least 2014.

So there is one main question:

1. how does one convince margins to exponentiate this:
margins, dydx(gre) at(gpa=(`r(mean)')) expression(predict(xb)) // .0027976 ...coefficient for gre
That is; output di exp(.0027976) // 1.0028015

And a niggling question (that probably involves a difference of probabilities - I haven't checked):

2. what is being calculated by margins, dydx(gre) at(gpa=( `r(mean)')) expression(exp(predict(xb))) // .0013802
The hoped-for answer is 1.0028015

[Andrew Musau]

1. how does one convince margins to exponentiate this:
margins, dydx(gre) at(gpa=(`r(mean)')) expression(predict(xb)) // .0027976 ...coefficient for gre
That is; output di exp(.0027976) // 1.0028015

You don't need to if your objective is to plot. coefplot will do this transformation for you.

use https://stats.idre.ucla.edu/stat/stata/dae/binary.dta, clear
logit admit c.gre##c.gpa in 1/200
margins, dydx(gre) at((mean) gpa) expression(predict(xb)) post
est sto m1

logit admit c.gre##c.gpa in 201/400
margins, dydx(gre) at((mean) gpa) expression(predict(xb)) post
est sto m2

logit admit c.gre##c.gpa 
margins, dydx(gre) at((mean) gpa) expression(predict(xb)) post
est sto m3

coefplot m1 m3 m2, eform mlabel(@b) mlabpos(12) xlab("") ylab("") ytitle(Marginal Odds Ratios (GRE Scores)) leg(order(1 "Group A" 5 "Group B"  3 "All groups"))

2. what is being calculated by margins, dydx(gre) at(gpa=( `r(mean)')) expression(exp(predict(xb))) // .0013802
The hoped-for answer is 1.0028015

Isn't this discussed in the linked thread?

[Janine Stubbs]

coefplot m1 m3 m2, eform mlabel(@b) mlabpos(12) xlab("") ylab("") ytitle(Marginal Odds Ratios (GRE Scores)) leg(order(1 "Group A" 5 "Group B" 3 "All groups"))

This is a very nice way to produce a graphical output of the exponentiated coefficients with confidence intervals (CIs). Thankyou very much!!
However, to calculate all the required information, the post-estimation routine still necessitates "manual" calculation using lincom.
To be clear, it involves:

logit
margins, dydx() ... to get the desired coefficients on the log-linear scale in tabular form (which are not needed for reporting)
coefplot (), eform ... to display the ORs and CIs
lincom ... to calculate and tabulate the ORs and CIs

By contrast, after reg, everything one wishes to report is available in tabular form after margins;

reg
margins, dydx() ... to get the desired coefficients and CIs in tabular form
marginsplot

So the question remains: can margins be used to create a tabular output of exponentiated coefficients and their CIs?
As in #4; this does not do it:
margins, dydx(gre) at(gpa=( `r(mean)')) expression(exp(predict(xb))) // .0013802 .... this is not the desired output, it would be 1.0028

[Andrew Musau]

So the question remains: can margins be used to create a tabular output of exponentiated coefficients and their CIs?

Without a doubt. The linked thread in #5 shows you how to do this in the presence of categorical RHS variables. One can work out the modification for continuous RHS variables.

However, to calculate all the required information, the post-estimation routine still necessitates "manual" calculation using lincom.

Not necessarily. Again, my reading from #1 is that you need the exponentiated estimates for purposes of graphing, and this is demonstrated in #7. If you also need to output these, use a command like estout from SSC, which can perform this transformation. Otherwise, be specific as to what your objective is.

use https://stats.idre.ucla.edu/stat/stata/dae/binary.dta, clear
logit admit c.gre##c.gpa in 1/200
margins, dydx(gre) at((mean) gpa) expression(predict(xb)) post
est sto m1

logit admit c.gre##c.gpa in 201/400
margins, dydx(gre) at((mean) gpa) expression(predict(xb)) post
est sto m2

logit admit c.gre##c.gpa
margins, dydx(gre) at((mean) gpa) expression(predict(xb)) post
est sto m3

*ssc install estout, replace
esttab m*, eform ci

Res.:

. esttab m*, eform ci

------------------------------------------------------------------------------------------
                                (1)                       (2)                       (3)  
                                                                                          
------------------------------------------------------------------------------------------
gre                           1.005**                   1.001                     1.003**
                      [1.002,1.008]             [0.998,1.004]             [1.001,1.005]  
------------------------------------------------------------------------------------------
N                               200                       200                       400  
------------------------------------------------------------------------------------------
Exponentiated coefficients; 95% confidence intervals in brackets
* p<0.05, ** p<0.01, *** p<0.001