problems interpreting melogit postestimation commands

[Jaye Derrick]

I am trying to interpret an interaction effect that was modeled using melogit in Stata 14.2. After the initial model, I ran a margins command and marginsplot command to graph the interaction. Then I ran a margins command to try to interpret the "simple slopes" (or re-centered predictors). It seems that the follow-up commands are very sensitive to how they are specified, which makes me think that I must not be interpreting something correctly. See the different coding used:

melogit symp_bin i.gender c.time##c.cT1anx || sub: time
    estat ic
    melogit, or
    margins, at(time=(0 1 2) cT1anx=(-0.5 0.5))
        marginsplot, noci xlabel(0 1 2) ylabel(0(.1).5) name(symptoms,        ///
        replace) title("") note("Symptoms as a function of Initial Anxiety X Time")
        
/*Estimate the significance of the time-symptom association at different levels of anxiety*/
margins, at(cT1anx=(-0.5 0.5)) dydx(time)
margins, at(cT1anx=(-0.5 0.5)) dydx(time) predict(xb)
margins, at(cT1anx=(-0.5 0.5)) dydx(time) exp(predict(xb))

Regarding the graph: I have attached the graph that Stata generates. Stata appears to be graphing the marginal mean (according to the y-axis title). Are these the probabilities? The actual values *look* like probabilities (i.e., they range from .06 to .35 or so). I have tried manually calculating the logged odds by taking the values from the original equation, substituting in low and high values for anxiety (-0.5 and 0.5) and values for time (0, 1, 2). Then I tried exponentiating them to get the odds. Then I've used exp(xb)/(1+exp(xb)) to get the probabilities. None of these values match the graph that I am getting from Stata (logged odds range from -4.61 to -1.56, odds range from .010 to .233, and probabilities range from .010 to .189). I also get very different graphs from Stata if I specify predict(xb) at the end of the margins command immediately prior to the marginsplot command (see second graph).

Regarding the simple slopes: the default is the marginal predicted mean. It appears that I am obtaining a probability? If I request predict(xb), I get very different values with very different significance levels. The exp(predict(xb)) option returns the exact same values as the predict(xb) option. If I manually exponentiate those values or use the formula exp(xb)/(1+exp(xb)) to manually turn those values into probabilities, they don't match the values from the default output. I suppose I could just report the probabilities that Stata is giving me, but I am accustomed to reporting the ORs (i.e., the exponentiated logged odds) for my simple slopes. I don't understand how to interpret probabilities in this way, particularly when they are negative. And of course, I'd still rather understand what is going on.

[Clyde Schechter]

Stata appears to be graphing the marginal mean (according to the y-axis title). Are these the probabilities?

Yes. That is exactly what they are.

I have tried manually calculating the logged odds by taking the values from the original equation, substituting in low and high values for anxiety (-0.5 and 0.5) and values for time (0, 1, 2). Then I tried exponentiating them to get the odds. Then I've used exp(xb)/(1+exp(xb)) to get the probabilities. None of these values match the graph that I am getting from Stata

Without seeing the actual code you used to do these manual calculations it is impossible to say what accounts for the difference. From a generic point of view bear in mind that you calculation has to a) include the estimated values of the random intercepts and random slopes that -melogit- estimated, and, b) for each of the data points on the graph you have to do this calculation having first modified the values of time and T1anx to the values corresponding to those points in every observation of the data set, and then average over the data set. That is a lot of hand calculation to do. It is complicated enough that few of us would be likely to get it right on any of the first few dozen tries.

I also get very different graphs from Stata if I specify predict(xb) at the end of the margins command immediately prior to the marginsplot command (see second graph).

Well, of course you do. Your plotting xb instead of the predicted probability. -xb- ignores the random effects, and also has not undergone the inverse logit transformation to get to the probability metric. There is no reason to expect those to be the same.

Regarding the simple slopes: the default is the marginal predicted mean. It appears that I am obtaining a probability? If I request predict(xb), I get very different values with very different significance levels.

When no -predict()- or -exp()- option is specified and you are estimating a slope with the dydx() option, then you are getting estimates of the derivative of the probability of symp_bin with respect to time, at the specified values of T1anx, usually more simply referred to as the marginal effect of time on the probability, conditional on the values of T1anx. Again, if you are using the -predict(xb)- option, then you are telling Stata you want the marginal effect of time on the linear predictor, not on the predicted probability, at those times. Again, there is no reason to expect them to be the same.

The exp(predict(xb)) option returns the exact same values as the predict(xb) option.

As an option in -margins-, exp(whatever) means that you want to get the marginal means of the evaluation of the expression specified by whatever. This -exp()- has nothing to do with the exponential function -exp()- used in other contexts. This notation is truly confusing and unfortunate. For this reason it is usually wise to write out the option name as -expression(predict(xb))- rather than using the allowed, but confusing, abbreviation -exp()- Perhaps StataCorp will change the name of this option in the future, or will disallow the exp abbreviation of it. If you wanted to get the marginal means of exp(xb) [meaning e^xb], the code for that would be -expression(exp(xb))-, which, if you would like to doubly confuse yourself, you could abbreviated to -exp(exp(xb))-!

I don't understand how to interpret probabilities in this way, particularly when they are negative. And of course, I'd still rather understand what is going on.

If you are getting negative values with -predict(xb)- there is nothing surprising about that--that is expected when the corresponding probability is < 0.5. You should not be getting negative predicted probabilities out of -margins- when you just use the default prediction. If you are, please show the exact code and the exact Stata output so we can try to see what is going wrong.

If I request predict(xb), I get very different values with very different significance levels.

Actually, you should just ignore the significance levels that you get from the output of -margins-, except possibly when getting a marginal effect with the dydx() option. If you are running, for example, -margins-, the p-values there are tests of the null hypothesis that the predicted probability is zero. That null hypothesis is virtually never of any interest, so the p-value tells you nothing useful. If you run -margins, predict(xb)-, you are testing the hypothesis that xb is zero, or equivalently, the hypothesis that the predicted probability is 0.5. Such a hypothesis test, again, is almost never of any use. The hypothesis tests that are usually of interest in these models are about the coefficients, marginal effects, or odds ratios (null value 1, not zero!). Those come from either the regression coefficient table itself, or the output of -test- or -testnl- commands, or -margins, dydx()-.

[Jaye Derrick]

Clyde, thank you so much for such a detailed response. A lot of this information is very helpful, but I can also see a few points where I was unclear. First, thank you for the confirmation that the default output, for the margins (with marginsplot) and the margins, dydx() commands are probabilities.

Without seeing the actual code you used to do these manual calculations it is impossible to say what accounts for the difference. From a generic point of view bear in mind that you calculation has to a) include the estimated values of the random intercepts and random slopes that -melogit- estimated, and, b) for each of the data points on the graph you have to do this calculation having first modified the values of time and T1anx to the values corresponding to those points in every observation of the data set, and then average over the data set. That is a lot of hand calculation to do. It is complicated enough that few of us would be likely to get it right on any of the first few dozen tries.

I found part of the problem. I was substituting zero for gender because I had not specified "atmeans" - but if I substitute the mean for gender, I do get the same values for predict(xb). Here is what I did. First I took the output from the command entered in my first post.

Mixed-effects logistic regression               Number of obs     =      1,680
Group variable:             sub                 Number of groups  =        560

                                                Obs per group:
                                                              min =          3
                                                              avg =        3.0
                                                              max =          3

Integration method: mvaghermite                 Integration pts.  =          7

                                                Wald chi2(4)      =      61.29
Log likelihood = -705.54111                     Prob > chi2       =     0.0000
---------------------------------------------------------------------------------
       symp_bin |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
----------------+----------------------------------------------------------------
                |
       1.gender |   .4473889   .3072977     1.46   0.145    -.1549035    1.049681
           time |   .1205314   .2545121     0.47   0.636    -.3783031    .6193659
         cT1anx |   3.150358   .4443132     7.09   0.000      2.27952    4.021195
                |
c.time#c.cT1anx |  -.4446609   .3487235    -1.28   0.202    -1.128147    .2388247
                |
          _cons |   -3.25392   .3566085    -9.12   0.000     -3.95286    -2.55498
----------------+----------------------------------------------------------------
sub             |
       var(time)|   2.454211   1.182356                      .9546226     6.30946
      var(_cons)|   4.620073   1.313481                      2.646387    8.065741
---------------------------------------------------------------------------------
LR test vs. logistic model: chi2(2) = 162.18              Prob > chi2 = 0.0000

I created the equation: -3.25392 + .4473889gender + .1205314time + 3.150358cT1anx - .4446609timeXcT1anx
Then I substituted in values of -0.5 and 0.5 for cT1anx and values of 0/1/2 for time to get six points - the log odds. These *should* match the values for the margins command if I include predict(xb):

margins, at(time=(0 1 2) cT1anx=(-0.5 0.5)) predict(xb)

And indeed, if I also substitute the mean for gender (0.5), I do obtain the same values. And so now my log-odds graph will match Stata's log-odds graph. However, I then calculated the probabilities of these values using the equation exp(xb)/(1+exp(xb)). My probabilities do not match Stata's default output. Stata gave me:

Expression   : Marginal predicted mean, predict()

1._at        : time            =           0
               cT1anx          =         -.5

2._at        : time            =           0
               cT1anx          =          .5

3._at        : time            =           1
               cT1anx          =         -.5

4._at        : time            =           1
               cT1anx          =          .5

5._at        : time            =           2
               cT1anx          =         -.5

6._at        : time            =           2
               cT1anx          =          .5

------------------------------------------------------------------------------
             |            Delta-method
             |     Margin   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         _at |
          1  |   .0504349   .0106424     4.74   0.000     .0295762    .0712935
          2  |   .2990639   .0270461    11.06   0.000     .2460546    .3520733
          3  |    .092045   .0127467     7.22   0.000     .0670619     .117028
          4  |   .3128177   .0235245    13.30   0.000     .2667104    .3589249
          5  |   .1751497   .0228361     7.67   0.000     .1303917    .2199076
          6  |   .3460419   .0281498    12.29   0.000     .2908692    .4012146
------------------------------------------------------------------------------

But I obtained:
1 | 0.010
2 | 0.189
3 | 0.014
4 | 0.174
5 | 0.019
6 | 0.160

So my probabilities graph still does not match Stata's.

Well, of course you do. Your plotting xb instead of the predicted probability. -xb- ignores the random effects, and also has not undergone the inverse logit transformation to get to the probability metric. There is no reason to expect those to be the same.

Sorry, yes, you are absolutely right that I would get different graphs when plotting the predicted probability and -xb-. I meant to state that I got a different logged odds graph when I graphed it manually than when Stata graphed it (now fixed!), and that I get a different probabilities graph when I graph it manually than when Stata graphs it (still a problem).

When no -predict()- or -exp()- option is specified and you are estimating a slope with the dydx() option, then you are getting estimates of the derivative of the probability of symp_bin with respect to time, at the specified values of T1anx, usually more simply referred to as the marginal effect of time on the probability, conditional on the values of T1anx. Again, if you are using the -predict(xb)- option, then you are telling Stata you want the marginal effect of time on the linear predictor, not on the predicted probability, at those times. Again, there is no reason to expect them to be the same.

Again, I was unclear. You are absolutely right that I should get different values when I ask for the marginal effect of time vs. the predicted probability. What I meant to ask was - if they are calculated for the same values of cT1anx, shouldn't I get the same significance levels?

Here is the predicted probabilities:

. margins, at(cT1anx=(-0.5 0.5)) dydx(time)

Average marginal effects                        Number of obs     =      1,680
Model VCE    : OIM

Expression   : Marginal predicted mean, predict()
dy/dx w.r.t. : time

1._at        : cT1anx          =         -.5

2._at        : cT1anx          =          .5

------------------------------------------------------------------------------
             |            Delta-method
             |      dy/dx   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
time         |
         _at |
          1  |   .0574286   .0119669     4.80   0.000     .0339738    .0808833
          2  |   .0166723   .0175446     0.95   0.342    -.0177144    .0510591
------------------------------------------------------------------------------

And the marginal effects:

. margins, at(cT1anx=(-0.5 0.5)) dydx(time) predict(xb)

Average marginal effects                        Number of obs     =      1,680
Model VCE    : OIM

Expression   : Linear prediction, fixed portion only, predict(xb)
dy/dx w.r.t. : time

1._at        : cT1anx          =         -.5

2._at        : cT1anx          =          .5

------------------------------------------------------------------------------
             |            Delta-method
             |      dy/dx   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
time         |
         _at |
          1  |   .3428618   .3798612     0.90   0.367    -.4016524    1.087376
          2  |  -.1017991   .2146214    -0.47   0.635    -.5224492    .3188511
------------------------------------------------------------------------------

Shouldn't I get the same significance levels? Also, if I try to calculate the probabilities by taking the marginal effect of time and using the equation exp(xb)/(1+exp(xb)), I do not get the same probabilities that Stata gave me:

. di exp(.3428618)/(1+exp(.3428618))
.58488552

. di exp(-.1017991)/(1+exp(-.1017991))
.47457218

So this is one of those places where I fear I am misinterpreting the output Stata is giving me. For the probabilities, I would say (overly simplistically) that at low levels of cT1anx, the probability of experiencing symptoms goes up over time, prob = .057, p < .001. At high levels of cT1anx, the probability of experiencing symptoms is not associated with time, prob = .017, p = .342. Is this correct? Why are the significant levels for the log odds different, and why can't I calculated the probabilities from the log odds?

As an option in -margins-, exp(whatever) means that you want to get the marginal means of the evaluation of the expression specified by whatever. This -exp()- has nothing to do with the exponential function -exp()- used in other contexts. This notation is truly confusing and unfortunate. For this reason it is usually wise to write out the option name as -expression(predict(xb))- rather than using the allowed, but confusing, abbreviation -exp()- Perhaps StataCorp will change the name of this option in the future, or will disallow the exp abbreviation of it. If you wanted to get the marginal means of exp(xb) [meaning e^xb], the code for that would be -expression(exp(xb))-, which, if you would like to doubly confuse yourself, you could abbreviated to -exp(exp(xb))-!

Thank you! Good to know! I will check this out.

If you are getting negative values with -predict(xb)- there is nothing surprising about that--that is expected when the corresponding probability is < 0.5. You should not be getting negative predicted probabilities out of -margins- when you just use the default prediction. If you are, please show the exact code and the exact Stata output so we can try to see what is going wrong.

I didn't in this particular instance, but we did just recently when my student was trying to put together a poster for a conference. I'll go back and dig that up.

Actually, you should just ignore the significance levels that you get from the output of -margins-, except possibly when getting a marginal effect with the dydx() option. If you are running, for example, -margins-, the p-values there are tests of the null hypothesis that the predicted probability is zero. That null hypothesis is virtually never of any interest, so the p-value tells you nothing useful. If you run -margins, predict(xb)-, you are testing the hypothesis that xb is zero, or equivalently, the hypothesis that the predicted probability is 0.5. Such a hypothesis test, again, is almost never of any use. The hypothesis tests that are usually of interest in these models are about the coefficients, marginal effects, or odds ratios (null value 1, not zero!). Those come from either the regression coefficient table itself, or the output of -test- or -testnl- commands, or -margins, dydx()-.

Yes, sorry, again that was me being unclear. I have only been looking at the significant levels when running -margins, dydx()-.

[Clyde Schechter]

Well, you still haven't actually shown the calculations for your hand-worked version of the mean predicted probabilities, and it's hard to know what you did in what order without seeing that. The likely problem is that logit() is a non-linear function, so the mean of invlogit(xb) is not the same thing as invlogit (mean of xb). In partcular, when -margins- calculates the marginal mean probabilities, it calculates invlogit(xb) and then averages those values over the entire data set. I'm guessing you just took invlogit() of your marginal mean xb values--which is doing it in the wrong order. [invlogit(x) is just a shorter way of writing exp(x)/(1+exp(x)), if you're not familiar with that notation. Also, Stata has a built-in invlogit() function, so you can simplify your expressions in these calculations by using it instead of the more cumbersome exponential over 1 plus exponential.]

Shouldn't I get the same significance levels? Also, if I try to calculate the probabilities by taking the marginal effect of time and using the equation exp(xb)/(1+exp(xb)), I do not get the same probabilities that Stata gave me:

No. Here you showed your calculations, and the problem is precisely what I pointed out in the paragraph just above.. You are just calculating invlogit(marginal mean xb). You have to calculate invlogit(xb) for every observation in the data set and then take the mean of that.

What I meant to ask was - if they are calculated for the same values of cT1anx, shouldn't I get the same significance levels?

No. Again, the non-linearity of invlogit() is biting you here. These two hypothesis tests are only distantly related. In particular, the standard error gets inflated or shrunk by the non-linearity of invlogit() depending on the base values you're starting from, so even though the marginal effect on xb can be exactly zero only when the marginal effect on predicted probability is zero and vice versa, the resulting z-statistics can be wildly different.

I think I've covered the remaining questions. If I missed one, do post back.

[Jaye Derrick]

Thank you again, Clyde, for your help on this issue. I think my difficulty has been, as you say, that I have been taking the invlogit on my marginal mean xb values and not on the individual values in the dataset.

Well, you still haven't actually shown the calculations for your hand-worked version of the mean predicted probabilities, and it's hard to know what you did in what order without seeing that.

I apologize, I was doing the calculations in Excel and I don't see an easy way to post those here. But my calculations were the same as the ones I posted using the -di- command from Stata, so the Excel calculations would also be in the wrong order.

The likely problem is that logit() is a non-linear function, so the mean of invlogit(xb) is not the same thing as invlogit (mean of xb). In partcular, when -margins- calculates the marginal mean probabilities, it calculates invlogit(xb) and then averages those values over the entire data set. I'm guessing you just took invlogit() of your marginal mean xb values--which is doing it in the wrong order.

My understanding was that you could use the following equations to move back and forth between the log-odds, the odds, and the probability:

logit(p) = log(p/(1-p))= β₀ + β₁*x1 + ... + β_k*xk
p= exp(β₀ + β₁*x1 + ... + β_k*xk)/(1+exp(β₀ + β₁*x1 + ... + β_k*xk))

and this is how logistic regression (but not melogit) is interpreted at, e.g., the following site: http://www.ats.ucla.edu/stat/mult_pk...odds_ratio.htm

Does this only hold for the original model, and not for the -margins- and -margins, dydx- commands? I understand that the -margins, dydx- values in particular do not describe the actual line, but rather the tangent line (i.e., the derivative of the "real" line). Is this why you can do the transformations back and forth with, say, the model intercept, but not for the values that are returned with the -margins- and the -margins, dydx- commands?

If you are getting negative values with -predict(xb)- there is nothing surprising about that--that is expected when the corresponding probability is < 0.5. You should not be getting negative predicted probabilities out of -margins- when you just use the default prediction. If you are, please show the exact code and the exact Stata output so we can try to see what is going wrong.

I have gone back to find the example of the time that I obtained a negative probability. This was actually obtained using the postestimation commands for -logit- and not -melogit-, in case that makes a difference to this problem. We ran:

logit arelapse    aageC asexC awhiteC asatC                                     ///
                c.aemoC##c.asurrogateC                                         ///
               , or vce (robust)
                    margins, at(asurrogateC=(-1.84 1.84) aemoC=(-1.20 1.20))       
                    marginsplot, noci                                    
                    margins, at(asurrogateC=(-1.84 1.84)) dydx(aemoC)            
                    margins, at(aemoC=(-1.20 1.20)) dydx(asurrogateC)

And obtained:

Logistic regression                             Number of obs     =         58
                                                Wald chi2(7)      =       9.09
                                                Prob > chi2       =     0.2462
Log pseudolikelihood = -27.722642               Pseudo R2         =     0.1885

---------------------------------------------------------------------------------------
                      |               Robust
             arelapse | Odds Ratio   Std. Err.      z    P>|z|     [95% Conf. Interval]
----------------------+----------------------------------------------------------------
                aageC |    1.08186   .0586351     1.45   0.147     .9728304    1.203108
                asexC |   .6979263   .5377129    -0.47   0.641     .1541724    3.159457
              awhiteC |   .2297783   .2465551    -1.37   0.171     .0280517    1.882173
                asatC |   .4656281   .3868329    -0.92   0.358     .0913863    2.372452
                aemoC |   .7965674   .3722589    -0.49   0.626     .3187364    1.990735
          asurrogateC |   .6380485   .1517882    -1.89   0.059     .4002742    1.017068
                      |
c.aemoC#c.asurrogateC |   1.909272   .5788254     2.13   0.033     1.053932    3.458778
                      |
                _cons |   4.924505   2.415756     3.25   0.001     1.882777    12.88031
---------------------------------------------------------------------------------------

Expression   : Pr(arelapse), predict()

1._at        : aemoC           =        -1.2
               asurrogateC     =       -1.84

2._at        : aemoC           =        -1.2
               asurrogateC     =        1.84

3._at        : aemoC           =         1.2
               asurrogateC     =       -1.84

4._at        : aemoC           =         1.2
               asurrogateC     =        1.84

------------------------------------------------------------------------------
             |            Delta-method
             |     Margin   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         _at |
          1  |   .9725043   .0377969    25.73   0.000     .8984238    1.046585
          2  |   .4205053   .1756692     2.39   0.017     .0761999    .7648107
          3  |     .62953   .1027347     6.13   0.000     .4281736    .8308863
          4  |   .8204894   .0966209     8.49   0.000     .6311159    1.009863
------------------------------------------------------------------------------

Expression   : Pr(arelapse), predict()
dy/dx w.r.t. : aemoC

1._at        : asurrogateC     =       -1.84

2._at        : asurrogateC     =        1.84

------------------------------------------------------------------------------
             |            Delta-method
             |      dy/dx   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
aemoC        |
         _at |
          1  |  -.1522301    .052205    -2.92   0.004    -.2545499   -.0499102
          2  |   .1776963   .1176537     1.51   0.131    -.0529006    .4082933
------------------------------------------------------------------------------

Expression   : Pr(arelapse), predict()
dy/dx w.r.t. : asurrogateC

1._at        : aemoC           =        -1.2

2._at        : aemoC           =         1.2

------------------------------------------------------------------------------
             |            Delta-method
             |      dy/dx   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
asurrogateC  |
         _at |
          1  |  -.1429973   .0377642    -3.79   0.000    -.2170138   -.0689809
          2  |   .0505572   .0333102     1.52   0.129    -.0147295     .115844
------------------------------------------------------------------------------

The output for the original -logit- command and the -margins- and -marginsplot- commands used to graph the interaction all make sense (we were supplied with ORs in the original output and probabilities for the graph). But then the -margins, dydx- commands return what appear to be negative probabilities. Is this just that, for those with low asurrogateC, the probability of arelapse goes down by about 15% for each one unit increase in aemoC? And for those with low aemoC, the probability of arelapse goes down by about 14% for each one unit increase in asurrogateC Because of course, you wouldn't obtain a negative OR, even if it were going down - it would just be between 0 and 1. So I had assumed that you wouldn't see a negative probability either. But I don't know how you would know that the direction of the line is negative if you only had positive probabilities.

[Clyde Schechter]

My understanding was that you could use the following equations to move back and forth between the log-odds, the odds, and the probability:

logit(p) = log(p/(1-p))= β₀ + β₁*x1 + ... + β_k*xk
p= exp(β₀ + β₁*x1 + ... + β_k*xk)/(1+exp(β₀ + β₁*x1 + ... + β_k*xk))

and this is how logistic regression (but not melogit) is interpreted at, e.g., the following site: http://www.ats.ucla.edu/stat/mult_pk...odds_ratio.htm

The free use of the logit transformation to go between odds ratios and probabilities applies to a single number. But when you are calculating average probabilities and average values of xb, then you cannot directly transform the average of one to the average of the other, because logit is non-linear. You have to laborious calculate the xb's or p's at the individual level, transform at the individual level, and the calculate the average of those. So for coefficients and odds ratios it's simple. For average marginal results, it's complicated, and this is true both for marginal values of p and xb, and also for marginal effects: sample averages are complicated and not subject to simple non-linear transformations.

Does this only hold for the original model, and not for the -margins- and -margins, dydx- commands? I understand that the -margins, dydx- values in particular do not describe the actual line, but rather the tangent line (i.e., the derivative of the "real" line). Is this why you can do the transformations back and forth with, say, the model intercept, but not for the values that are returned with the -margins- and the -margins, dydx- commands?

You are correct that the dydx calculations calculate the slope of the tangent line. But again, the main difference here is that the output of -margins- is always a sample average of individual statistics, whereas the output of -melogit- or similar regression commands is a single "holistic" estimate applicable to the sample as a whole. Of course, the situation with marginal effects is further complicated by the fact that the derivative of invlogit(xb) is not the same as invlogit(derivative of xb). So non-linearity intervenes again, this time in a slightly different way.

As for the negative values you show in your -margins- outputs, none of those are marginal probabilities. They are marginal effects, calculated using the -dydx()- option. There is no reason to expect those to be positive, nor less than 1. They can range over the entire real number line.

[Jaye Derrick]

Thank you so much, Clyde! This has been incredibly helpful very enlightening.