2
$\begingroup$

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

$\endgroup$
  • 4
    $\begingroup$ There's clearly a software component in this question, but I think we can keep this open since it goes deeper into the question of how to compute standard errors and confidence intervals for bounded range statistics (predicted probabilities). $\endgroup$ – StasK Jul 31 '15 at 14:35
5
$\begingroup$

margins computes standard errors from nonlinear predictions using the delta-method, and as donlelek points out, it also uses a normal approximation for computing confidence intervals.

Consider a slight variation on donlelek's example.

sysuse auto, clear
logistic foreign mpg
margins, predict(pr) nopvalues

The result from margins is

Predictive margins                              Number of obs     =         74
Model VCE    : OIM

Expression   : Pr(foreign), predict()

--------------------------------------------------------------
             |            Delta-method
             |     Margin   Std. Err.     [95% Conf. Interval]
-------------+------------------------------------------------
       _cons |   .2972973   .0487662      .2017172    .3928773
--------------------------------------------------------------

Let's call this marginal prediction $m$.

Here we verify by hand the confidence limits produced by margins.

. display .2972973 - invnormal(.975)*.0487662
.2017173

. display .2972973 + invnormal(.975)*.0487662
.3928773

$m$ is the mean of the predicted probabilities. The predicted probability for observation $i$ is

$$ p_i = \frac{1}{1 + \exp(-z_i)} $$

where $z_i$ is the corresponding linear prediction

$$ z_i = b_0 + b_1*x_i $$

$x_i$ is the corresponding value of mpg, and $b$ is the vector of regression coefficients. Thus $m$ is

$$ m = \frac{1}{N}\sum_{i=1}^N p_i $$

$m$ is not a simple transformation of the marginal linear prediction. This is the general situation that margins was developed to handle. As such, the only method for computing a confidence interval for $m$ is the normal approximation using standard errors computed via the delta-method.

margins does not provide an option that will transform the confidence limits of a marginal linear prediction.

However, donlelek's question inspired me to write transform_margins, which will transform the point estimates and confidence limits from margins results with a user specified expression. Here is donlelek's example using transform_margins:

sysuse auto
logistic foreign mpg
margins, at(mpg=(5(5)40)) predict(xb)
transform_margins invlogit(@)

Here is the output from transform_margins

----------------------------------------------
             |         b         ll         ul
-------------+--------------------------------
         _at |
          1  |  .0271183   .0042517    .153952
          2  |  .0583461   .0151856   .1993467
          3  |  .1210596   .0511398   .2603462
          4  |  .2344013    .144129   .3575909
          5  |  .4049667   .2710298   .5547246
          6  |  .6020462   .3688264   .7966118
          7  |  .7707955   .4519781     .93203
          8  |  .8820117   .5300384   .9802168
----------------------------------------------

Type the following commands in Stata to install transform_margins.

net from http://www.stata.com/users/jpitblado
net describe transform_margins
net install transform_margins
help transform_margins
$\endgroup$
  • $\begingroup$ Thanks for your answer @Jeff Pitblado, just to clarify, there is no way to do this directly with margins, and the correct way of doing it would be to get the linear predictions and then transform to the desired scale with your transform_margins command. $\endgroup$ – donlelek Aug 1 '15 at 23:45
  • 1
    $\begingroup$ margins does not have an option that provides the kind of confidence interval calculations you want. The normal approximation, using delta method standard errors, is the only option. In your example, I can see why you do not like this, so I provided transform_margins as a convenience. $\endgroup$ – Jeff Pitblado Aug 3 '15 at 14:44
  • $\begingroup$ Thanks again @Jeff Pitblado, as an alternative; would it be ok to use the confidence intervals provided by margins and just: replace ll = 0 if ll < 0 and replace ll = 1 if ll > 1, and the same for ul? $\endgroup$ – donlelek Aug 3 '15 at 14:51
  • $\begingroup$ I do not see why not, provided you document that the confidence intervals are based on a normal approximation and that you are cropping at the end-points. $\endgroup$ – Jeff Pitblado Aug 3 '15 at 22:23
1
$\begingroup$

Confidence intervals coming from predicted values in logistic regression models should be nonsymmetric and bounded between 0 and 1. The key is that the inverse link (the expit function) should be applied at the very end of calculation. Try the lincom command instead.

Margins uses marginal standardization which employs the delta method approximation to standard errors for predicted values on the probability scale.

$\endgroup$
  • $\begingroup$ Thanks for your answer @AdamO but my question is if there is a way of doing it right with margins. Also, this MRE is not my actual model and using lincom to do what I want would be tedious (I would if it's the only option) but it seems to me like a serious bug if Stata get the CIs wrong for all non-linear predictions. $\endgroup$ – donlelek Jul 31 '15 at 12:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.