logit postestimation

[Nigel Moore]

I have a fairly simple dataset in which a binary response is matched against drug dose:

* Example generated by -dataex-. To install: ssc install dataex
clear
input float dose int preg
  0 1
  0 1
  0 1
  0 1
  0 1
  0 1
  0 1
  0 1
  0 1
  0 1
  0 1
  0 0
 85 1
 85 1
 85 1
 85 1
 85 1
 85 1
 85 1
 85 1
 85 1
 85 1
 85 1
 85 1
250 1
250 1
250 1
250 1
250 1
250 1
250 1
250 1
250 1
250 1
250 1
250 1
500 1
500 1
500 1
500 1
500 1
500 1
500 1
500 0
500 0
500 0
500 0
500 0
end

What I am hoping to achieve is a determination on whether any of the doses are significantly different to control (dose=0)

. logit preg dose 

Iteration 0:   log likelihood = -18.084968  
Iteration 1:   log likelihood =  -14.81676  
Iteration 2:   log likelihood = -14.003143  
Iteration 3:   log likelihood = -13.990435  
Iteration 4:   log likelihood = -13.990401  
Iteration 5:   log likelihood = -13.990401  

Logistic regression                             Number of obs     =         48
                                                LR chi2(1)        =       8.19
                                                Prob > chi2       =     0.0042
Log likelihood = -13.990401                     Pseudo R2         =     0.2264

------------------------------------------------------------------------------
        preg |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
        dose |  -.0071943    .002991    -2.41   0.016    -.0130565    -.001332
       _cons |   4.138469   1.261881     3.28   0.001     1.665227     6.61171
------------------------------------------------------------------------------

If I try margins I get an error:

. margins dose
factor 'dose' not found in list of covariates
r(322);

I tried margins, at(), but the output was nonsensical because the lowest dose groups cannot be significantly different:

. margins, at(dose=(85,250,500))

Adjusted predictions                            Number of obs     =         48
Model VCE    : OIM

Expression   : Pr(preg), predict()

1._at        : dose            =          85

2._at        : dose            =         250

3._at        : dose            =         500

------------------------------------------------------------------------------
             |            Delta-method
             |     Margin   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         _at |
          1  |   .9714451   .0286189    33.94   0.000     .9153532    1.027537
          2  |   .9121283   .0512337    17.80   0.000      .811712    1.012545
          3  |   .6321235    .136214     4.64   0.000      .365149     .899098
------------------------------------------------------------------------------

Am I just barking up the wrong tree, trying to get an outcome that's simply not possible?

[Chen Samulsion]

I understand that -margins- only accept factor variables or interactions in marginlist. for example: margins i.sex i.group i.sex#i.group

egen dose4=cut(dose), group(4)
logit preg i.dose4
margins dose4
margins i.dose4
margins, at(dose4=(0(1)3))

. egen dose4=cut(dose), group(4)

. logit preg i.dose4

note: 1.dose4 != 0 predicts success perfectly;
1.dose4 omitted and 12 obs not used.

note: 2.dose4 != 0 predicts success perfectly;
2.dose4 omitted and 12 obs not used.

Iteration 0:  Log likelihood = -13.496043  
Iteration 1:  Log likelihood = -11.702513  
Iteration 2:  Log likelihood =  -11.59278  
Iteration 3:  Log likelihood = -11.592351  
Iteration 4:  Log likelihood = -11.592351  

Logistic regression        Number of obs =     24
        LR chi2(1)    =   3.81
        Prob > chi2   = 0.0510
Log likelihood = -11.592351        Pseudo R2     = 0.1411

        
preg  Coefficient  Std. err.      z    P>z    [95% conf. interval]
        
dose4 
1            0  (empty)
2            0  (empty)
3    -2.061423     1.1974    -1.72    0.085    -4.408283    .2854374

_cons    2.397895   1.044466     2.30    0.022    .3507797    4.445011
        

. margins dose4

Adjusted predictions        Number of obs = 24
Model VCE: OIM

Expression: Pr(preg), predict()

        
Delta-method
Margin   std. err.      z    P>z    [95% conf. interval]
        
dose4 
0     .9166667   .0797856    11.49    0.000    .7602898    1.073044
3     .5833333   .1423188     4.10    0.000    .3043937     .862273
        

. margins i.dose4

Adjusted predictions        Number of obs = 24
Model VCE: OIM

Expression: Pr(preg), predict()

        
Delta-method
Margin   std. err.      z    P>z    [95% conf. interval]
        
dose4 
0     .9166667   .0797856    11.49    0.000    .7602898    1.073044
3     .5833333   .1423188     4.10    0.000    .3043937     .862273
        

. margins, at(dose4=(0(1)3))

Adjusted predictions        Number of obs = 24
Model VCE: OIM

Expression: Pr(preg), predict()
1._at: dose4 = 0
2._at: dose4 = 1
3._at: dose4 = 2
4._at: dose4 = 3

        
Delta-method
Margin   std. err.      z    P>z    [95% conf. interval]
        
_at 
1     .9166667   .0797856    11.49    0.000    .7602898    1.073044
2            .  (not estimable)
3            .  (not estimable)
4     .5833333   .1423188     4.10    0.000    .3043937     .862273

[Nigel Moore]

Thank you, Chen, I appreciate that! I managed to get the same result(s) without quite the same manipulation (which suits my simpler outlook!):

. logit preg i.dose 

note: 85.dose != 0 predicts success perfectly
      85.dose dropped and 12 obs not used

note: 250.dose != 0 predicts success perfectly
      250.dose dropped and 12 obs not used

Iteration 0:   log likelihood = -13.496043  
Iteration 1:   log likelihood = -11.702513  
Iteration 2:   log likelihood =  -11.59278  
Iteration 3:   log likelihood = -11.592351  
Iteration 4:   log likelihood = -11.592351  

Logistic regression                             Number of obs     =         24
                                                LR chi2(1)        =       3.81
                                                Prob > chi2       =     0.0510
Log likelihood = -11.592351                     Pseudo R2         =     0.1411

------------------------------------------------------------------------------
        preg |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
        dose |
         85  |          0  (empty)
        250  |          0  (empty)
        500  |  -2.061423     1.1974    -1.72   0.085    -4.408283    .2854374
             |
       _cons |   2.397895   1.044466     2.30   0.022     .3507797    4.445011
------------------------------------------------------------------------------

. margins,  at(dose=(85,250,500))

Adjusted predictions                            Number of obs     =         24
Model VCE    : OIM

Expression   : Pr(preg), predict()

1._at        : dose            =          85

2._at        : dose            =         250

3._at        : dose            =         500

------------------------------------------------------------------------------
             |            Delta-method
             |     Margin   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         _at |
          1  |          .  (not estimable)
          2  |          .  (not estimable)
          3  |   .5833333   .1423188     4.10   0.000     .3043937     .862273
------------------------------------------------------------------------------

If I read it correctly though, the outcome of the initial logit is not significant, either for the regression as a whole (p=0.0510) or the 500 dose group (p=0.085)
But the results of margins is significant for 500 (p=0.000).

Is that significance the result of not quitting at the right time? Or does it mean something more?

[Frode Andre]

I am on Stata 19.5. I wrote:

logit preg dose
margins, at(dose = (85 250 500))

And got the table:

. margins, at(dose = (85, 250, 500))

Adjusted predictions Number of obs = 48
Model VCE: OIM

Expression: Pr(preg), predict()
1._at: dose = 85
2._at: dose = 250
3._at: dose = 500

------------------------------------------------------------------------------
| Delta-method
| Margin std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
_at |
1 | .9714451 .0286189 33.94 0.000 .9153532 1.027537
2 | .9121283 .0512337 17.80 0.000 .811712 1.012545
3 | .6321235 .136214 4.64 0.000 .365149 .899098
------------------------------------------------------------------------------

[Joseph Coveney]

I have a fairly simple dataset in which a binary response is matched against drug dose:
. . . What I am hoping to achieve is a determination on whether any of the doses are significantly different to control (dose=0)

You can use

logit preg i.dose, asis nolog
margins r.dose, contrast(pveffects nowald)

for an each-dose-versus-zero set of Wald tests, or you can use

logit preg i.dose, asis nolog
estimates store Full

constraint define 1 _b[85.dose] = _b[0.dose]
constraint define 2 _b[250.dose] = _b[0.dose]
constraint define 3 _b[500.dose] = _b[0.dose]

forvalues constraint = 1/3 {
    quietly logit preg i.dose, asis constraints(`constraint')
    lrtest Full, df(1)
}

for the corresponding set of likelihood-ratio tests.

Both give you the same bottom line, viz., the two lower levels of dose don't do much, and it's the highest dose level that differs (lower) from the zero dose level. Complete do-file (with your dataset) and corresponding log file are attached for your perusal.

But you've got a couple of things going on in this dataset. First, your sample size is a bit small (n = 12) for this kind of response and, second, which is related to the first, you've got complete separation going on at the lower two dose levels. Both of these kind of complicate matters using asymptotic iterative maximum likelihood.

One alternative is to perform pairwise crosstabulations using something like

tabulate dose preg if inlist(dose, 0, 85), exact

etc., but that's pretty low-power.

Another is to use something like the user-written penlogit from Stata Journal (search penlogit), but that won't allow factor variables, and so you'll need to generate your own set of categories manually using, for example, tabulate dose, generate(dos) and omitting the first category. Even then, performing contrasts after fitting the model is a bit more involved.

A final alternative is to use official bayes: logit with, say, log-F(1,1) priors on the coefficients and follow up with visual inspection of the posterior probability density histograms or with the official bayes: postestimation commands to look at contrasts.

[Nigel Moore]

Thank you, Joseph.

I assume that the margins/contrast doesn't control for multiple comparisons, so this is also correct in order to do so:

. logit preg i.dose, asis nolog

Logistic regression                             Number of obs     =         48
                                                LR chi2(1)        =      12.99
                                                Prob > chi2       =     0.0003
Log likelihood = -11.592351                     Pseudo R2         =     0.3590

------------------------------------------------------------------------------
        preg |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
        dose |
         85  |   16.27105          .        .       .            .           .
        250  |   16.27105          .        .       .            .           .
        500  |  -2.061423     1.1974    -1.72   0.085    -4.408283    .2854374
             |
       _cons |   2.397895   1.044466     2.30   0.022     .3507797    4.445011
------------------------------------------------------------------------------
Note: 0 failures and 24 successes completely determined.

. margins r.dose, contrast(pveffects nowald) mcompare(bon)

Contrasts of adjusted predictions
Model VCE    : OIM

Expression   : Pr(preg), predict()

---------------------------
             |    Number of
             |  Comparisons
-------------+-------------
        dose |            3
---------------------------

-----------------------------------------------------
             |            Delta-method    Bonferroni
             |   Contrast   Std. Err.      z    P>|z|
-------------+---------------------------------------
        dose |
  (85 vs 0)  |   .0833333   .0797856     1.04   0.889
 (250 vs 0)  |   .0833333   .0797856     1.04   0.889
 (500 vs 0)  |  -.3333333   .1631575    -2.04   0.123
-----------------------------------------------------

@ Frode Andre: Yes, I also got that result (post #1). It was Chen who got me on the right track with i.dose, then Joseph with r.dose.