How to do post-estimation plot of discrete time analysis?

[Man Yang]

Hello, I am using discrete time analysis on a project to examine the first occurrence of an event. Below are my syntax

. logit EVENT PERIOD pr2 pr3 i.BirthYear ib5.ethnicrace ELL female MRLEVEL PracIndp PersSocSkll ChallBeh In
> tegration WellBeing if MR==1, iterate(50)       

Iteration 0:   log likelihood = -17774.871  
Iteration 1:   log likelihood = -16968.503  
Iteration 2:   log likelihood = -15866.331  
Iteration 3:   log likelihood = -15860.161  
Iteration 4:   log likelihood =  -15860.15  
Iteration 5:   log likelihood =  -15860.15  

Logistic regression                               Number of obs   =     252934
                                                  LR chi2(36)     =    3829.44
                                                  Prob > chi2     =     0.0000
Log likelihood =  -15860.15                       Pseudo R2       =     0.1077

-------------------------------------------------------------------------------------------
                    EVENT |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
--------------------------+----------------------------------------------------------------
                   PERIOD |  -.5272217   .0454329   -11.60   0.000    -.6162686   -.4381749
                      pr2 |   .0112451   .0062711     1.79   0.073    -.0010461    .0235363
                      pr3 |  -.0002635   .0002424    -1.09   0.277    -.0007386    .0002116
                          |
                BirthYear |
                    1982  |  -.3731549   .1285857    -2.90   0.004    -.6251782   -.1211316
                    1983  |  -1.002686   .1392994    -7.20   0.000    -1.275708   -.7296639
                    1984  |  -1.057614   .1334958    -7.92   0.000     -1.31926   -.7959666
                    1985  |  -1.427331    .138246   -10.32   0.000    -1.698288   -1.156374
                    1986  |  -1.953165   .1454343   -13.43   0.000    -2.238211   -1.668119
                    1987  |  -2.034847   .1447087   -14.06   0.000    -2.318471   -1.751223
                    1988  |  -2.358076   .1497545   -15.75   0.000    -2.651589   -2.064562
                    1989  |  -2.725663    .154599   -17.63   0.000    -3.028671   -2.422654
                    1990  |  -3.068246   .1587418   -19.33   0.000    -3.379375   -2.757118
                    1991  |  -3.478518   .1631533   -21.32   0.000    -3.798293   -3.158743
                    1992  |  -3.999524   .1671678   -23.93   0.000    -4.327167   -3.671881
                    1993  |  -4.207497   .1694396   -24.83   0.000    -4.539593   -3.875402
                    1994  |  -4.361056    .172867   -25.23   0.000    -4.699869   -4.022243
                    1995  |  -4.668352   .1796521   -25.99   0.000    -5.020464   -4.316241
                    1996  |  -4.518742   .1802471   -25.07   0.000     -4.87202   -4.165464
                    1997  |  -4.887824   .1928031   -25.35   0.000    -5.265711   -4.509937
                    1998  |  -4.858595   .2020433   -24.05   0.000    -5.254592   -4.462597
                    1999  |  -4.844775   .2265945   -21.38   0.000    -5.288892   -4.400658
                    2000  |  -4.496112   .3219503   -13.97   0.000    -5.127123   -3.865101
                          |
               ethnicrace |
        African American  |    .786026   .0460669    17.06   0.000     .6957366    .8763154
          Asian American  |  -.4359661   .0953101    -4.57   0.000    -.6227705   -.2491617
                Latino/a  |   .0629631   .0508165     1.24   0.215    -.0366354    .1625615
Native American-Hawaiian  |   .5642024    .193542     2.92   0.004     .1848671    .9435376
                   Other  |  -.0456702   .1145911    -0.40   0.690    -.2702646    .1789243
                 Unknown  |   .9808834   .1137712     8.62   0.000      .757896    1.203871
                          |
                      ELL |  -.9358915   .0602884   -15.52   0.000    -1.054055   -.8177284
                   female |  -.1778972   .0369286    -4.82   0.000    -.2502758   -.1055185
                  MRLEVEL |   -.034015   .0318231    -1.07   0.285    -.0963871    .0283571
                 PracIndp |  -.0938565   .0280988    -3.34   0.001    -.1489292   -.0387839
              PersSocSkll |    .329711   .0292657    11.27   0.000     .2723512    .3870707
                 ChallBeh |  -.3234386   .0158676   -20.38   0.000    -.3545385   -.2923388
              Integration |  -.0569525   .0229395    -2.48   0.013    -.1019131   -.0119919
                WellBeing |  -.0546078   .0097765    -5.59   0.000    -.0737694   -.0354462
                    _cons |   2.783287   .2250984    12.36   0.000     2.342102    3.224472
-------------------------------------------------------------------------------------------

in which EVENT is the outcome, PERIOD is the time variable, pr2 and pr3 are the quadratic and cubic term of PERIOD, the rest of the variables are just covariates I throw in the model. if MR=1 limits the original sample to those who have MR equals to 1.

I have two questions:
1. since all I need is to control for birth year but I don't really need the model output for each birth year. Is there a way for me to let STATA omit the output for birth year?
2. after getting the result, now I want to make plot of fitted hazard probability and I thought about two ways: one is using predict and the other is using margins. However, neither one was able to produce a meaningful plot that I am expecting. For instance, below I included the plots generated using predict and margins specifically.

predict p_33

line p_33 PERIOD, sort ///
xtitle("Age in years") ytitle("Hazard Probability") ///
xlabel(1 "15" 2 "16" 3 "17" 4 "18" 5 "19" 6 "20" 7 "21" 8 "22" 9 "23" 10 "24" 11 "25" 12 "26" 13 "27" 14 "28" 15 "29" 16 "30" 17 "31" 18 "32" 19 "33") ///
region(lstyle(none) lcolor(none) color(none)))

The above plot is really strange and I wonder whether it is because of the covariates I added to the model.

Or if using margins, I got the following plot

margins, at(PERIOD=(1(1)19)) nose
marginsplot, recast(line) xlabel(1 "15" 2 "16" 3 "17" 4 "18" 5 "19" 6 "20" 7 "21" 8 "22" 9 "23" 10 "24" 11 "25" 12 "26" 13 "27" 14 "28" 15 "29" 16 "30" 17 "31" 18 "32" 19 "33") noci

Any clues please? Thanks!

[Oded Mcdossi]

As for you first question, try:

estimates table, drop(i.BirthYear)  b(%4.2fc) star

As for the second, I don't understand why margins results are not meaningful.

[Man Yang]

Hi Oded, thanks for your answer. For the second one, the reason why I am suspicious about margins plot is because the curve is too smooth. For instance, below I ran the unadjusted model using the same limiter.

 logit EVENT PERIOD pr2 pr3 if MR==1

Iteration 0:   log likelihood = -25456.813  
Iteration 1:   log likelihood = -25215.932  
Iteration 2:   log likelihood = -25195.311  
Iteration 3:   log likelihood = -25195.245  
Iteration 4:   log likelihood = -25195.245  

Logistic regression                               Number of obs   =     462322
                                                  LR chi2(3)      =     523.14
                                                  Prob > chi2     =     0.0000
Log likelihood = -25195.245                       Pseudo R2       =     0.0103

------------------------------------------------------------------------------
       EVENT |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
      PERIOD |  -.3246224   .0308826   -10.51   0.000    -.3851512   -.2640936
         pr2 |   .0339812   .0045182     7.52   0.000     .0251256    .0428368
         pr3 |  -.0012593   .0001848    -6.81   0.000    -.0016214   -.0008971
       _cons |  -3.734809   .0543489   -68.72   0.000    -3.841331   -3.628288
------------------------------------------------------------------------------

Once I got the model result, I use -predict- (instead of margins)

predict p_1

line p_1 PERIOD, sort ///
xtitle("Age in years") ytitle("Hazard Probability") ///
xlabel(1 "15" 2 "16" 3 "17" 4 "18" 5 "19" 6 "20" 7 "21" 8 "22" 9 "23" 10 "24" 11 "25" 12 "26" 13 "27" 14 "28" 15 "29" 16 "30" 17 "31" 18 "32" 19 "33")

I got the following plot

And then I tried to use -margins-

logit EVENT PERIOD pr2 pr3 if MR==1
margins, at(PERIOD=(1(1)19)) nose
marginsplot, recast(line) xlabel(1 "15" 2 "16" 3 "17" 4 "18" 5 "19" 6 "20" 7 "21" 8 "22" 9 "23" 10 "24" 11 "25" 12 "26" 13 "27" 14 "28" 15 "29" 16 "30" 17 "31" 18 "32" 19 "33") noci title("AUT only (Full)")

You see I got a different plot? In my perspective, I think the first plot is more accurate, in other words, the plot generated using -predict- is the correct one. That's why I am thinking of using -predict- instead of -margins-. But the problem is -predict- works in the unadjusted model but not working properly in the adjusted model (shown in the very first plot of this post). Do you know what is happening here? Thanks.

[Oded Mcdossi]

You should use Stata factor variable notation ( help fvvarlist ) so margins will recognize the qubic and quadratic terms of period.

logit EVENT c.PERIOD##c.PERIOD##c.PERIOD if MR==1
margins, at(PERIOD=(1(1)19)) nose
marginsplot, recast(line) xlabel(1 "15" 2 "16" 3 "17" 4 "18" 5 "19" 6 "20" 7 "21" 8 "22" 9 "23" 10 "24" 11 "25" 12 "26" 13 "27" 14 "28" 15 "29" 16 "30" 17 "31" 18 "32" 19 "33") noci title("AUT only (Full)")

[Man Yang]

Hi Oded, yes! That works! Thanks for the response. But I still found that once adding all the other covariates (in the words, in the adjusted model) using margins after the model still produces the second plot in #1. Does it mean that in the full model, the curve is supposed to be like that? Or do you think there is a way for me to use -predict- to generate plot after getting the adjusted model results? I just want to test whether the -predict- and -margins- will produce the same plot.

[Oded Mcdossi]

I can't say much without seeing the code and the output, but your interpretation might be correct if the effects of the qubic and the quadratic terms become zero.

[Man Yang]

When you say the effects of cubic and quadratic terms become zero, do you mean their coefficients = 0 in the model output?

[Oded Mcdossi]

Yes, this is exactly what I meant. Anyway, predict and margins would produce the same results under the same conditions.

Added: Since you are using the if qualifier in the model you should also use it in the predict post estimation command. Margins, however, take the if qualifier into account when calculating the marginal effects so I would recommend using margins.