Confidence intervals with margins command after logistic regression

[Derek Price]

I am using Stata 13.1 to fit a logistic model and I am getting confidence intervals below 0 and above 1 when I predict probabilities using the margins command.

MRE:

sysuse auto, clear
* simple logistic regression
logit foreign mpg
* get predicted probabilities
margins, at(mpg=(5(5)40)) predict(pr)
* same result with expression
margins, at(mpg=(5(5)40)) exp(invlogit(predict(xb)))

The output is:

------------------------------------------------------------------------------
             |            Delta-method
             |     Margin   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         _at |
          1  |   .0271183   .0252542     1.07   0.283    -.0223792    .0766157
          2  |   .0583461   .0389888     1.50   0.135    -.0180704    .1347627
          3  |   .1210596   .0509373     2.38   0.017     .0212244    .2208948
          4  |   .2344013   .0547344     4.28   0.000      .127124    .3416787
          5  |   .4049667   .0743318     5.45   0.000      .259279    .5506543
          6  |   .6020462   .1162814     5.18   0.000     .3741387    .8299536
          7  |   .7707955   .1266899     6.08   0.000     .5224879    1.019103
          8  |   .8820117   .1004224     8.78   0.000     .6851874    1.078836
------------------------------------------------------------------------------

I'm trying to figure out what is Stata doing so I extracted the results with:

* save table as matrix
mat pr = r(table)
mat p = pr' //transpose matrix

Then I predicted the log odds (linear prediction), saved the results and transformed to probability (inverse logit)

* predict log odds
margins, at(mpg=(5(5)40)) exp(predict(xb))

* save table as matrix
mat tab = r(table)
mat t = tab'

clear
svmat t

* transform logit to probability and display

gen prob = invlogit(t1)
gen problo = invlogit(t5)
gen probhi = invlogit(t6)

format %9.3f prob*
list prob* in 1/8, noobs clean // correct confidence intervals?!

prob problo probhi
0.027 0.004 0.154
0.058 0.015 0.199
0.121 0.051 0.260
0.234 0.144 0.358
0.405 0.271 0.555
0.602 0.369 0.797
0.771 0.452 0.932
0.882 0.530 0.980

Now if I take the results from the margins output in probability scale predict(pr) and (wrongly?) use the SE from that output to produce CIs I get the same as Stata:

clear

* convert results to data
svmat p

* generate confidence intervals (the wrong way?)
gen problo = p1 - (1.96*p2)
gen probhi = p1 + (1.96*p2)
format %9.3f p*
list p5 problo p6 probhi in 1/8, noobs clean

  p5   problo  p6   probhi
-0.022 -0.022 0.077 0.077
-0.018 -0.018 0.135 0.135
 0.021  0.021 0.221 0.221
 0.127  0.127 0.342 0.342
 0.259  0.259 0.551 0.551
 0.374  0.374 0.830 0.830
 0.522  0.522 1.019 1.019
 0.685  0.685 1.079 1.079

Sorry it took me so long but here is the question: Is there a way of getting the right CIs using the margins command directly? and does this problem happens with other GLMs?

Thanks

Derek Price, DVM
PhD Candidate
UPEI

[Jorge Eduardo Perez Perez]

What margins seems to be doing is using the delta method on the predicted probabilities itself, while it's better to just get the standard errors from the log odds, as stated in this FAQ.

Generating the confidence intervals for the index and then converting them to probabilities to get confidence intervals for the predicted probabilities is better than estimating the standard error of the predicted probabilities and then generating the confidence intervals directly from that standard error. The distribution of the predicted index is closer to normality than the predicted probability.

[Derek Price]

Hi All,

Thanks for the answers, I also posted the question at http://stats.stackexchange.com/quest...rgins-in-stata and got a nice answer from Jeff Pitblado who kindly wrote the command transform_margins that takes the linear prediction and transforms it to an arbitrary scale.

To install:

net from http://www.stata.com/users/jpitblado
net describe transform_margins
net install transform_margins
help transform_margins

[Dennis Foeste]

@Jeff Pitblado
Thanks for the transform_margins ado-file

I've got four questions related to the posts above:
I am using Stata 13.1 and ran into a problem using transform_margins to get the 'correct' CIs:
Comparing the margins output and the transform_margins output I recognized different Margins (not only CIs)
This problem occures when I specify the exposure option with nbreg or poisson (c. f. example below)

I guess this is due to what predict(xb) computes, when exposure or offset is specified.
This leads me to my questions:
to solve this problem (1) How can I fix this? (get the predictions and CI right before using transform_margins?)
for a better understanding (2) How does margings take offset/exposure option into account ?
(just found this http://www.stata.com/support/faqs/st...-after-offset/ for mfx )

Follow-up Question to the original post
(3) Do I have to be concerned about other options of margins like dydx or eyex, not to report the 'appropiate' CIs?
(4) Is there a way to take into account the asobserved-option of margins to use it with transform_margins?
[For example: If I wanted to state "margins agecat, asobserved"
instead of "margins agecat, atmeans" (as I did in the example below)]

Best regards & thanks in advance.

For illustration purposes I stick closely to an example given in the transform_margins help-file:



. webuse dollhill3
(Doll and Hill (1966))

. qui poisson deaths i.smokes i.agecat, exp(pyears)

. margins agecat, atmeans

Adjusted predictions Number of obs = 10
Model VCE : OIM

Expression : Predicted number of events, predict()
at : 0.smokes = .5 (mean)
1.smokes = .5 (mean)
1.agecat = .2 (mean)
2.agecat = .2 (mean)
3.agecat = .2 (mean)
4.agecat = .2 (mean)
5.agecat = .2 (mean)

------------------------------------------------------------------------------
| Delta-method
| Margin Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
agecat |
35-44 | 7.878883 1.374695 5.73 0.000 5.184531 10.57323
45-54 | 34.75047 3.484188 9.97 0.000 27.92159 41.57935
55-64 | 109.0374 8.385145 13.00 0.000 92.60285 125.472
65-74 | 224.6804 17.82199 12.61 0.000 189.7499 259.6109
75-84 | 318.7103 29.96264 10.64 0.000 259.9846 377.436
------------------------------------------------------------------------------

. qui poisson deaths i.smokes i.agecat, exp(pyears)

. qui margins agecat, predict(xb)

. transform_margins exp(@)
-----------------------------------------------
| b ll ul
-------------+---------------------------------
agecat |
35-44 | 4.627041 3.286905 6.513579
45-54 | 20.40795 16.767 24.83954
55-64 | 64.03455 55.07515 74.45143
65-74 | 131.9483 112.9498 154.1425
75-84 | 187.1694 155.6724 225.0392
-----------------------------------------------

[smith Jason]

sysuse auto, clear
* simple logistic regression logit foreign mpg
* get predicted probabilities margins, at(mpg=(5(5)40)) predict(pr)
* same result with expression
margins, at(mpg=(5(5)40)) exp(invlogit(predict(xb)))
I don't understand the last line Stata code, especially on the part "exp(invlogit(predict(xb)))".
Can someone help me with it? Thank you!

[Leonardo Guizzetti]

This is just asking margins to transform each result from the log-odds scale to the probability scale.

[smith Jason]

This is just asking margins to transform each result from the log-odds scale to the probability scale.

Could you explain this part? I really can not understand. My understanding is the result of exp() is odds ratio. How does it become predicted probability?

[Leonardo Guizzetti]

See -help margins- for details of the options. exp is short for expression, as in the expression margins is calculating. It is not exponential. Margins is computing the inverse-logit of the linear prediction of the logit model.

[smith Jason]

Thank you!