Logistic odds ratio in loops

[Anouk Klootwijk]

Dear statalist,

I'm using a logistic regression in a loop, which I later on use in a coefplot
First I used margins with the command below and this resulted in the first figure:

webuse lbw
global outcomes_binary "low smoke"
global outcomes_continuous "bwt"
global outcomes "$outcomes_binary $outcomes_continuous"
global covariates "age lwt ptl ht ui"

foreach y of global outcomes_binary {
logit `y' $covariates
margins, dydx(*) post
mat b_`y' = r(table)
mat b_`y' = b_`y'[1..6,1...]'
mat b_`y' = b_`y'[1...,1], b_`y'[1...,5], b_`y'[1...,6], b_`y'[1...,4]
}

foreach y of global outcomes_continuous {
reg `y' $covariates
margins, dydx(*) post
mat b_`y' = r(table)
mat b_`y' = b_`y'[1..6,1...]'
mat b_`y' = b_`y'[1...,1], b_`y'[1...,5], b_`y'[1...,6], b_`y'[1...,4]
}

local i = 1
foreach var of global covariates {
mat `var' = J(1,4,0)
foreach y of global outcomes {
mat `var' = `var' \ b_`y'[`i',1...]
}
mat `var' = `var'[2...,1...]
mat rownames `var' = $outcomes_binary $outcomes_continuous $outcomes_ordinal0 $outcomes_ordinal5 $outcomes_ordinal4
mat colnames `var' = `var' ci_l ci_h
local i = `i' + 1
}

coefplot matrix(age[,1]), ci((age[,2] age[,3])) ///
|| matrix(lwt[,1]), ci((lwt[,2] lwt[,3])) ///
|| matrix(ptl[,1]), ci((ptl[,2] ptl[,3])) ///
||, scheme(s1color) xline(0) byopts(row(1)) xlab(-200(200)200 ,format(%9.0f)) mlabel mlabposition(2) format(%9.3f)


I would like to change the marginal effects into odds ratios, but when I try this command:

foreach y of global outcomes_binary {
logistic `y' $covariates, or
mat b_`y' = r(table)
mat b_`y' = b_`y'[1..6,1...]'
mat b_`y' = b_`y'[1...,1], b_`y'[1...,5], b_`y'[1...,6], b_`y'[1...,4]
}

Something goes wrong and the figure is incomplete (smoke & bwt are missing).



Can anyone help me to fix the command for storing odds ratios in loops?
Thanks in advance!

Regards, Anouk

[William Lisowski]

Your problem occurs because in your original code a margins command follows the logit command, but a margins command does not follow the logistic command in your revised code.

[Anouk Klootwijk]

Ok thanks for your reply, do you have any suggestions how I could fix it?

[William Lisowski]

Compare your original code that used logit to your revised code that used logistic and see that the margins command is missing from your revised code. Type or paste the missing margins command into your revised code in the same place it appears in your original code.

[Anouk Klootwijk]

Thanks for your reply. I'm trying to change the margins command into the logistic command, because I'm trying to display odds ratios in the figure instead of margins. How can I make the logistic command for odds ratios work in the loop?

[William Lisowski]

The documentation provided by the help coefplot command tells us

   eform[(coeflist)] causes point estimates and confidence intervals to be exponentiated. This is
        useful if you want to plot hazard ratios (HR), incidence-rate ratios (IRR), odds ratios
        (OR), or relative-risk ratios (RRR). If eform is specified without arguments, then all
        coefficients of the model are exponentiated. To exponentiate only selected coefficients,
        specify coeflist as above for keep().

This suggests that what you want is apparently accomplished with your original code, including both the logit and margins commands, and then adding the appropriate eform options to your coefplot command. This is a start at doing so for your example.

coefplot matrix(age[,1]), ci((age[,2] age[,3])) eform($outcomes_binary) ///
|| matrix(lwt[,1]), ci((lwt[,2] lwt[,3])) eform($outcomes_binary) ///
|| matrix(ptl[,1]), ci((ptl[,2] ptl[,3])) eform($outcomes_binary) ////
|| , ...

[Anouk Klootwijk]

Thanks for your reply, unfortunately the figure does not display odds ratios when I try your command.

Is there another command to save odds ratios of logistic regression in matrix?

For example:
foreach y of global outcomes_binary {
logit `y' $covariates, or
mat b_`y' = r(table)
mat b_`y' = b_`y'[1..6,1...]'
mat b_`y' = b_`y'[1...,1], b_`y'[1...,5], b_`y'[1...,6], b_`y'[1...,4]
}

When I tried this and afterwards your command with eform in the coefplot, I had the same figure that I send earlier (figure 2) and other variables (smoke & bwt) were still missing...

Thanks in advance!

[Andrew Musau]

There is part of your code that you are not showing us because the code in #1 does not produce the same coefficients for "birthweight<2500g" and "smoked during pregnancy". In any case, running the code below following William's suggestion does result in odds ratios displayed for selected groups of coefficients. coefplot is from the Stata Journal, as you are asked to explain in FAQ Advice #12.

webuse lbw, clear
global outcomes_binary "low smoke"
global outcomes_continuous "bwt"
global outcomes "$outcomes_binary $outcomes_continuous"
global covariates "age lwt ptl ht ui"

foreach y of global outcomes_binary {
logit `y' $covariates, or
mat b_`y' = r(table)
mat b_`y' = b_`y'[1..6,1...]'
mat b_`y' = b_`y'[1...,1], b_`y'[1...,5], b_`y'[1...,6], b_`y'[1...,4]
}


coefplot matrix(age[,1]), ci((age[,2] age[,3]))  ///
|| matrix(lwt[,1]), ci((lwt[,2] lwt[,3]))  ///
|| matrix(ptl[,1]), ci((ptl[,2] ptl[,3])) ///
||, scheme(s1color) xline(0) byopts(row(1)) xlab(-200(200)200 ,format(%9.0f)) mlabel mlabposition(2) format(%9.3f) saving(gr1, replace)

coefplot matrix(age[,1]), ci((age[,2] age[,3])) eform($outcomes_binary) ///
|| matrix(lwt[,1]), ci((lwt[,2] lwt[,3])) eform($outcomes_binary) ///
|| matrix(ptl[,1]), ci((ptl[,2] ptl[,3])) eform($outcomes_binary) ///
||, scheme(s1color) xline(0) byopts(row(1)) xlab(-200(200)200 ,format(%9.0f)) mlabel mlabposition(2) format(%9.3f) saving(gr2, replace)

graph combine gr1.gph gr2.gph, col(1) scheme(s1color)

Res.:

[Anouk Klootwijk]

Thanks for your reply. I'm using Stata 15.1. When I use your code the odds ratios are still missing in the coefplot... Resulting in this figure:



What command can I use to visualize the odds ratios in the coefplot?

Thanks in advance!

[Andrew Musau]

I tried the code in Stata 15.1 and it works fine. Maybe you have an older version of coefplot, so try the following code first. If that does not fix it, update your Stata installation. Looking closer, the coefficients in my graph are off - but the fix is to extract the correct coefficients from the matrix, which you should be able to do.

ssc install coefplot, replace
webuse lbw
global outcomes_binary "low smoke"
global outcomes_continuous "bwt"
global outcomes "$outcomes_binary $outcomes_continuous"
global covariates "age lwt ptl ht ui"

foreach y of global outcomes_binary {
logit `y' $covariates
margins, dydx(*) post
mat b_`y' = r(table)
mat b_`y' = b_`y'[1..6,1...]'
mat b_`y' = b_`y'[1...,1], b_`y'[1...,5], b_`y'[1...,6], b_`y'[1...,4]
}

foreach y of global outcomes_continuous {
reg `y' $covariates
margins, dydx(*) post
mat b_`y' = r(table)
mat b_`y' = b_`y'[1..6,1...]'
mat b_`y' = b_`y'[1...,1], b_`y'[1...,5], b_`y'[1...,6], b_`y'[1...,4]
}

local i = 1
foreach var of global covariates {
mat `var' = J(1,4,0)
foreach y of global outcomes {
mat `var' = `var' \ b_`y'[`i',1...]
}
mat `var' = `var'[2...,1...]
mat rownames `var' = $outcomes_binary $outcomes_continuous $outcomes_ordinal0 $outcomes_ordinal5 $outcomes_ordinal4
mat colnames `var' = `var' ci_l ci_h
local i = `i' + 1
}

coefplot matrix(age[,1]), ci((age[,2] age[,3])) ///
|| matrix(lwt[,1]), ci((lwt[,2] lwt[,3])) ///
|| matrix(ptl[,1]), ci((ptl[,2] ptl[,3])) ///
||, scheme(s1color) xline(0) byopts(row(1)) xlab(-200(200)200 ///
,format(%9.0f)) mlabel mlabposition(2) format(%9.3f) saving(graph1, replace)

coefplot matrix(age[,1]), ci((age[,2] age[,3])) eform($outcomes_binary) ///
|| matrix(lwt[,1]), ci((lwt[,2] lwt[,3])) eform($outcomes_binary) ///
|| matrix(ptl[,1]), ci((ptl[,2] ptl[,3])) eform($outcomes_binary) ///
||, scheme(s1color) xline(0) byopts(row(1)) xlab(-200(200)200 , ///
format(%9.0f)) mlabel mlabposition(2) format(%9.3f) saving(graph2, replace)

graph combine graph1.gph graph2.gph, col(1) scheme(s1color)

[Anouk Klootwijk]

Thanks for your reply!
I was now able to make the same figure!

However, the labels of the odds ratios in the coefplot do not correspond with the actual odds ratios.
For example the odds ratios retrieved by this command show different numbers than the labels that are displayed in the figure:
foreach y of global outcomes_binary {
logit `y' $covariates, or
mat b_`y' = r(table)
mat b_`y' = b_`y'[1..6,1...]'
mat b_`y' = b_`y'[1...,1], b_`y'[1...,5], b_`y'[1...,6], b_`y'[1...,4]
}

How can I change the command to display the true odds ratios ?

Thanks in advance!

[Andrew Musau]

How can I change the command to display the true odds ratios ?

if you want odds ratios, plot the original coefficients. What you get are the exponentiated marginal effects which make no sense. As I said in #8, you do not present the exact code that generates the graph in #1. The code that you present is not it. Exit the loop and consider the following:

. webuse lbw
(Hosmer & Lemeshow data)

. logit low  age lwt ptl ht ui

Iteration 0:   log likelihood =   -117.336  
Iteration 1:   log likelihood = -105.93804  
Iteration 2:   log likelihood = -105.67556  
Iteration 3:   log likelihood = -105.67471  
Iteration 4:   log likelihood = -105.67471  

Logistic regression                             Number of obs     =        189
                                                LR chi2(5)        =      23.32
                                                Prob > chi2       =     0.0003
Log likelihood = -105.67471                     Pseudo R2         =     0.0994

------------------------------------------------------------------------------
         low |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         age |  -.0437459   .0343065    -1.28   0.202    -.1109853    .0234935
         lwt |   -.014534   .0066853    -2.17   0.030    -.0276368   -.0014311
         ptl |   .6842166   .3442857     1.99   0.047      .009429    1.359004
          ht |   1.892711   .7001172     2.70   0.007     .5205062    3.264915
          ui |   .7407733   .4479599     1.65   0.098    -.1372119    1.618759
       _cons |   1.655227   1.068401     1.55   0.121    -.4388003    3.749255
------------------------------------------------------------------------------

.
. margins, dydx(*) post

Average marginal effects                        Number of obs     =        189
Model VCE    : OIM

Expression   : Pr(low), predict()
dy/dx w.r.t. : age lwt ptl ht ui

------------------------------------------------------------------------------
             |            Delta-method
             |      dy/dx   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         age |  -.0082361   .0063741    -1.29   0.196    -.0207292    .0042569
         lwt |  -.0027363   .0012107    -2.26   0.024    -.0051094   -.0003633
         ptl |   .1288187   .0623697     2.07   0.039     .0065764    .2510611
          ht |   .3563441    .122653     2.91   0.004     .1159487    .5967396
          ui |   .1394668   .0820712     1.70   0.089    -.0213899    .3003234
------------------------------------------------------------------------------

.
. mat b_low = r(table)

.
. mat b_low = b_low[1..6,1...]'

.
. mat b_low = b_low[1...,1], b_low[1...,5], b_low[1...,6], b_low[1...,4]

.
. mat l b_low

b_low[5,4]
              b          ll          ul      pvalue
age  -.00823613  -.02072915    .0042569   .19631486
lwt  -.00273634  -.00510935  -.00036332   .02381863
ptl   .12881873   .00657639   .25106106   .03888459
 ht   .35634413   .11594868   .59673959   .00366898
 ui   .13946676  -.02138986   .30032337   .08925544

The highlighted value -0.008 is the marginal effect on age in the graph in #8. It seems to me that all you do to get the graph in #1 is to multiply the marginal effects on the binary variables by 100. If you did not write the code, I would advise you to go through it line by line (exiting the loop as I did above) to understand what it is doing. For your specific question, the following highlighted modifications will get you the graph in #1. You should extract the coefficients directly from the regression to plot the odds ratios.

webuse lbw
global outcomes_binary "low smoke"
global outcomes_continuous "bwt"
global outcomes "$outcomes_binary $outcomes_continuous"
global covariates "age lwt ptl ht ui"

foreach y of global outcomes_binary {
logit `y' $covariates
margins, dydx(*) post
mat b_`y' = r(table)
mat b_`y' = b_`y'[1..6,1...]'
mat b_`y' = b_`y'[1...,1], b_`y'[1...,5], b_`y'[1...,6], b_`y'[1...,4]
}

foreach y of global outcomes_continuous {
reg `y' $covariates
margins, dydx(*) post
mat b_`y' = r(table)
mat b_`y' = b_`y'[1..6,1...]'
mat b_`y' = b_`y'[1...,1], b_`y'[1...,5], b_`y'[1...,6], b_`y'[1...,4]
}
mat l b_bwt

local i = 1
foreach var of global covariates {
mat `var' = J(1,4,0)
foreach y of global outcomes{    
local mat= cond("`y'"!="bwt", "b_`y'[`i',1...]*100",  "b_`y'[`i',1...]")    
mat `var'=  `var'\ `mat'
}
mat `var' = `var'[2...,1...]
mat rownames `var' = $outcomes_binary $outcomes_continuous $outcomes_ordinal0 $outcomes_ordinal5 $outcomes_ordinal4
mat colnames `var' = `var' ci_l ci_h
local i = `i' + 1
}

coefplot matrix(age[,1]), ci((age[,2] age[,3])) ///
|| matrix(lwt[,1]), ci((lwt[,2] lwt[,3])) ///
|| matrix(ptl[,1]), ci((ptl[,2] ptl[,3])) ///
||, scheme(s1color) xline(0) byopts(row(1)) xlab(-200(200)200 ///
,format(%9.0f)) mlabel mlabposition(2) format(%9.3f) saving(graph1, replace)

[Anouk Klootwijk]

Thanks for your reply. The first graph is indeed showing the true margins multiplied by 100.

My question is about visualizing the odds ratios in a coefplot.

When I compute the odds ratios of low it results in this table:

logistic low age lwt ptl ht ui

Logistic regression Number of obs = 189
LR chi2(5) = 23.32
Prob > chi2 = 0.0003
Log likelihood = -105.67471 Pseudo R2 = 0.0994

------------------------------------------------------------------------------
low | Odds Ratio Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
age | .9571971 .032838 -1.28 0.202 .8949519 1.023772
lwt | .9855711 .0065888 -2.17 0.030 .9727416 .9985699
ptl | 1.982218 .6824495 1.99 0.047 1.009474 3.892316
ht | 6.637336 4.646913 2.70 0.007 1.682879 26.17789
ui | 2.097557 .9396213 1.65 0.098 .8717855 5.046821
_cons | 5.23427 5.592299 1.55 0.121 .6448095 42.48941
------------------------------------------------------------------------------

When I use your command, the odds ratio of low (birthweight <2500g) for age is displayed as 0.439 instead of the 0.957 in the table above.



What command can I use to put the true odds ratios in a coefplot?

Thanks in advance!

[Andrew Musau]

As I said in #12, it's a matter of plotting the coefficients. If I want to plot the odds ratios from your regression example, I would simply do

webuse lbw
logistic low age lwt ptl ht ui
coefplot, eform scheme(s1color) mlabel mlabposition(2) format(%9.3f)

Now, the only additional thing that the code in #1 does it to extract specific coefficients, put them in matrices and plot them. As you already do this for marginal effects, adapting it to the raw coefficients is straightforward. I do not have the time to illustrate step by step how you should adapt the code. If you are unable to figure it out, someone else in the forum may help. If not, ask a colleague who has experience with Stata.

[Andrew Musau]

Coming back to this, note that by adding the option -predict(xb)- to the margins command, you get back the coefficients. Therefore, you can use this to plot the coefficients via the margins results. That said, you often learn Stata programming by doing, so I would still encourage you to attempt to write code that extracts the coefficients directly.

webuse lbw
global outcomes_binary "low smoke"
global outcomes_continuous "bwt"
global outcomes "$outcomes_binary $outcomes_continuous"
global covariates "age lwt ptl ht ui"

foreach y of global outcomes_binary {
logit `y' $covariates
margins, dydx(*) post predict(xb)
mat b_`y' = r(table)
mat b_`y' = b_`y'[1..6,1...]'
mat b_`y' = b_`y'[1...,1], b_`y'[1...,5], b_`y'[1...,6], b_`y'[1...,4]
}

foreach y of global outcomes_continuous {
reg `y' $covariates
margins, dydx(*) post 
mat b_`y' = r(table)
mat b_`y' = b_`y'[1..6,1...]'
mat b_`y' = b_`y'[1...,1], b_`y'[1...,5], b_`y'[1...,6], b_`y'[1...,4]
}
mat l b_bwt

local i = 1
foreach var of global covariates {
mat `var' = J(1,4,0)
foreach y of global outcomes{       
mat `var'=  `var'\ b_`y'[`i',1...]
}
mat `var' = `var'[2...,1...]
mat rownames `var' = $outcomes_binary $outcomes_continuous $outcomes_ordinal0 $outcomes_ordinal5 $outcomes_ordinal4
mat colnames `var' = `var' ci_l ci_h
local i = `i' + 1
}


coefplot matrix(age[,1]), ci((age[,2] age[,3])) eform($outcomes_binary) ///
|| matrix(lwt[,1]), ci((lwt[,2] lwt[,3])) eform($outcomes_binary) ///
|| matrix(ptl[,1]), ci((ptl[,2] ptl[,3])) eform($outcomes_binary) ///
||, scheme(s1color) xline(0) byopts(row(1)) xlab(-200(200)200 , ///
format(%9.0f)) mlabel mlabposition(2) format(%9.3f)