How do I reference an estimate after lincom?

[Jeff Wooldridge]

Hi Everyone: I know there must be a simple solution to this but I haven't found it. I thought that using coeflegend would tell me how I reference an estimate after lincom. But that reference gives an error when I want to display it. It also gives an error when I use it in nlcom. I'd like to insert the result from lincome into nlcom. That also fails.

. qui poisson thefts i.sameblock#c.m8 i.sameblock#c.m9 i.sameblock#c.m10 ///
>         i.sameblock#c.m11 i.sameblock#c.m12 ///
>         i.sameblock#c.m8#c.bank_dm i.sameblock#c.m9#c.bank_dm i.sameblock#c.m10#c.bank_dm ///
>         i.sameblock#c.m11#c.bank_dm i.sameblock#c.m12#c.bank_dm ///
>         i.sameblock c.m8 c.m9 c.m10 c.m11 c.m12 ///
>         i.sameblock#c.bank_dm1 c.m8#c.bank_dm c.m9#c.bank_dm c.m10#c.bank_dm ///
>         c.m11#c.bank_dm c.m12#c.bank_dm, vce(cluster block)

. lincom (1.sameblock#c.m8 + 1.sameblock#c.m9 + 1.sameblock#c.m10 ///
>         + 1.sameblock#c.m11 + 1.sameblock#c.m12)/5, coeflegend

 ( 1)  .2*[thefts]1.sameblock#c.m8 + .2*[thefts]1.sameblock#c.m9 + .2*[thefts]1.sameblock#c.m10 + .2*[thefts]1.sameblock#c.m11 + .2*[thefts]1.sameblock#c.m12 = 0

------------------------------------------------------------------------------
      thefts | Coefficient  Legend
-------------+----------------------------------------------------------------
         (1) |   -1.28483  _b[(1)]
------------------------------------------------------------------------------

. di _b[(1)]      
[(1)] not found
r(111);

[Justin Niakamal]

Try

di r(estimate)

[Jeff Wooldridge]

Thanks Justin. I did that and it properly displays. I should have mentioned that. But I cannot feed r(estimate) into nlcom, for example.

[Jeff Wooldridge]

And using nlcom directly on the average produces a puzzling outcome, too:

. nlcom exp((1.sameblock#c.m8 + 1.sameblock#c.m9 + 1.sameblock#c.m10 ///
>         + 1.sameblock#c.m11 + 1.sameblock#c.m12)/5) - 1

       _nl_1: exp((1.sameblock#c.m8 + 1.sameblock#c.m9 + 1.sameblock#c.m10         + 1.sameblock#c.m11 + 1.sameblock#c.m12)/5) - 1

------------------------------------------------------------------------------
      thefts | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
       _nl_1 |          0  (omitted)
------------------------------------------------------------------------------

[Joro Kolev]

Professor Wooldridge, you did not include in #3 an example of what goes wrong.

Try to do it in two steps, save the -lincom- estimate first into a scalar or a macro, and then plug it in the -nlcom-.

What I think is that both lincom and nlcom are r-class, so trying to plug r-class returned value does not work because the second r-class overwrites it on the fly.

Here is an example session, with some meaningless calculations:

. sysuse auto
(1978 Automobile Data)

. reg price mpg head, noheader
------------------------------------------------------------------------------
       price |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         mpg |  -259.1057   58.42485    -4.43   0.000    -375.6015   -142.6098
    headroom |  -334.0215   399.5499    -0.84   0.406    -1130.701    462.6585
       _cons |   12683.31   2074.497     6.11   0.000     8546.885    16819.74
------------------------------------------------------------------------------

. lincom mpg + headroom

 ( 1)  mpg + headroom = 0

------------------------------------------------------------------------------
       price |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         (1) |  -593.1271   427.0515    -1.39   0.169    -1444.644    258.3893
------------------------------------------------------------------------------

. nlcom r(estimate)^2 + (_b[mpg])^3

expression (r(estimate)^2 + (_b[mpg])^3) evaluates to missing
r(498);

so this does not work.

Now in two steps:

.  reg price mpg head

      Source |       SS           df       MS      Number of obs   =        74
-------------+----------------------------------   F(2, 71)        =     10.44
       Model |   144280501         2  72140250.4   Prob > F        =    0.0001
    Residual |   490784895        71  6912463.32   R-squared       =    0.2272
-------------+----------------------------------   Adj R-squared   =    0.2054
       Total |   635065396        73  8699525.97   Root MSE        =    2629.2

------------------------------------------------------------------------------
       price |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         mpg |  -259.1057   58.42485    -4.43   0.000    -375.6015   -142.6098
    headroom |  -334.0215   399.5499    -0.84   0.406    -1130.701    462.6585
       _cons |   12683.31   2074.497     6.11   0.000     8546.885    16819.74
------------------------------------------------------------------------------

. lincom mpg + headroom

 ( 1)  mpg + headroom = 0

------------------------------------------------------------------------------
       price |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         (1) |  -593.1271   427.0515    -1.39   0.169    -1444.644    258.3893
------------------------------------------------------------------------------

. sca Myestimate = r(estimate)

. nlcom Myestimate^2 + (_b[mpg])^3

       _nl_1:  Myestimate^2 + (_b[mpg])^3

------------------------------------------------------------------------------
       price |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
       _nl_1 |  -1.70e+07   1.18e+07    -1.45   0.148    -4.01e+07     6019812
------------------------------------------------------------------------------

. dis -593.1271^2 + -259.1057^3
-17747059

.

Using a local would do the trick as well:

. lincom mpg + headroom

 ( 1)  mpg + headroom = 0

------------------------------------------------------------------------------
       price |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         (1) |  -593.1271   427.0515    -1.39   0.169    -1444.644    258.3893
------------------------------------------------------------------------------

. local Myestimate = r(estimate)

. nlcom `Myestimate'^2 + (_b[mpg])^3

       _nl_1:  -593.1271329026094^2 + (_b[mpg])^3

------------------------------------------------------------------------------
       price |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
       _nl_1 |  -1.77e+07   1.18e+07    -1.51   0.132    -4.08e+07     5316212
------------------------------------------------------------------------------

. dis -593.1271^2 + -259.1057^3
-17747059

[Leonardo Guizzetti]

-lincom- does not post results like other estimation commands, so the usual _b[] and _se[] mechanisms won't work. Justin's advice is one way to grab those results. Also note that -lincom- allows for exponentiation after computing the liinear combination (e.g., -eform-, -irr-) so you can get those coefficients. However, if you were to include additive constants in your -lincom- expression, then you can no longer request exponentiated coefficients, which may or may not matter for your case.

Then, you can use -nlcom- as a replacement for -lincom-. The advantages include estimating multiple transformations simultaneously and posting the estimated covariance matrix among them. The point estimates between -lincom- and -nlcom- will be the same, but you may find differences between mathematically equivalent expressions in terms of the standard errors and associated test statistic and p-value (as warned in the documentation).

Another thing to note is that the accepted syntax between -lincom- and -nlcom- are different. You will need to use the _b[] format to reference coefficients in -nlcom-, so your syntax in #4 should be:

nlcom exp((_b[1.sameblock#c.m8] + _b[1.sameblock#c.m9] + _b[1.sameblock#c.m10]  + _b[1.sameblock#c.m11] + _b[1.sameblock#c.m12])/5) - 1

[Leonardo Guizzetti]

I like Joro's example, but numeric literal values (from a scalar or local macro) are treated as having no error. If you cared about carrying through that error you could proceed with -nlcom- in steps, as in this meaningless example.

nlcom ((_b[mpg] + _b[headroom])^2) (_b[mpg]^3)
nlcom (a: (_b[mpg] + _b[headroom])^2) (b: _b[mpg]^3), post
nlcom _b[a] + _b[b]

[Jeff Wooldridge]

Thanks Leonardo. I did try nlcom in one step, and it gives an estimate of zero -- clearly wrong -- and no standard error. Joro's solution gives an estimate but no standard error, which doesn't help. I must be using the factor notation incorrectly or something, but the estimates and standard errors from the Poisson regression look fine.

. poisson thefts i.sameblock#c.m8 i.sameblock#c.m9 i.sameblock#c.m10 ///
>         i.sameblock#c.m11 i.sameblock#c.m12 ///
>         i.sameblock#c.m8#c.bank_dm i.sameblock#c.m9#c.bank_dm i.sameblock#c.m10#c.bank_dm ///
>         i.sameblock#c.m11#c.bank_dm i.sameblock#c.m12#c.bank_dm ///
>         i.sameblock c.m8 c.m9 c.m10 c.m11 c.m12 ///
>         i.sameblock#c.bank_dm1 c.m8#c.bank_dm c.m9#c.bank_dm c.m10#c.bank_dm ///
>         c.m11#c.bank_dm c.m12#c.bank_dm, vce(cluster block)
note: you are responsible for interpretation of noncount dep. variable.

Iteration 0:   log pseudolikelihood = -2193.2568  
Iteration 1:   log pseudolikelihood =  -2193.031  
Iteration 2:   log pseudolikelihood = -2192.9905  
Iteration 3:   log pseudolikelihood = -2192.9812  
Iteration 4:   log pseudolikelihood = -2192.9793  
Iteration 5:   log pseudolikelihood = -2192.9789  
Iteration 6:   log pseudolikelihood = -2192.9788  
Iteration 7:   log pseudolikelihood = -2192.9788  

Poisson regression                                      Number of obs =  7,008
                                                        Wald chi2(20) =      .
                                                        Prob > chi2   =      .
Log pseudolikelihood = -2192.9788                       Pseudo R2     = 0.0044

                                              (Std. err. adjusted for 876 clusters in block)
--------------------------------------------------------------------------------------------
                           |               Robust
                    thefts | Coefficient  std. err.      z    P>|z|     [95% conf. interval]
---------------------------+----------------------------------------------------------------
            sameblock#c.m8 |
                        1  |  -.9509519   .5886853    -1.62   0.106    -2.104754      .20285
                           |
            sameblock#c.m9 |
                        1  |  -1.991195   .5675957    -3.51   0.000    -3.103662   -.8787279
                           |
           sameblock#c.m10 |
                        1  |  -.6905526   .5825024    -1.19   0.236    -1.832236    .4511311
                           |
           sameblock#c.m11 |
                        1  |  -1.370658   .4301334    -3.19   0.001    -2.213704   -.5276118
                           |
           sameblock#c.m12 |
                        1  |  -1.420794   .4297411    -3.31   0.001    -2.263071   -.5785166
                           |
 sameblock#c.m8#c.bank_dm1 |
                        1  |   1.416137   .7036727     2.01   0.044      .036964     2.79531
                           |
 sameblock#c.m9#c.bank_dm1 |
                        1  |   2.208649   .6871923     3.21   0.001      .861777    3.555521
                           |
sameblock#c.m10#c.bank_dm1 |
                        1  |   .5120065   .6841694     0.75   0.454     -.828941    1.852954
                           |
sameblock#c.m11#c.bank_dm1 |
                        1  |   1.571855   .6242858     2.52   0.012     .3482774    2.795433
                           |
sameblock#c.m12#c.bank_dm1 |
                        1  |   1.926151   .5602086     3.44   0.001     .8281624     3.02414
                           |
               1.sameblock |  -.7829916    .269717    -2.90   0.004    -1.311627    -.254356
                        m8 |   .1517601   .0924087     1.64   0.101    -.0293577    .3328778
                        m9 |    .037247   .0984474     0.38   0.705    -.1557063    .2302003
                       m10 |   .1228981   .0907106     1.35   0.175    -.0548914    .3006875
                       m11 |   .0557554   .0943717     0.59   0.555    -.1292098    .2407206
                       m12 |   .1058912   .0925673     1.14   0.253    -.0755373    .2873197
                           |
      sameblock#c.bank_dm1 |
                        0  |    .226014   .1824086     1.24   0.215    -.1315003    .5835283
                        1  |  -12.81342   .7611178   -16.83   0.000    -14.30518   -11.32165
                           |
           c.m8#c.bank_dm1 |  -.5486128   .3112156    -1.76   0.078    -1.158584    .0613585
                           |
           c.m9#c.bank_dm1 |  -.0876349   .3225214    -0.27   0.786    -.7197651    .5444954
                           |
          c.m10#c.bank_dm1 |   .1041835     .27923     0.37   0.709    -.4430972    .6514643
                           |
          c.m11#c.bank_dm1 |   -.144526   .4268554    -0.34   0.735    -.9811471    .6920952
                           |
          c.m12#c.bank_dm1 |  -.4988221   .3260773    -1.53   0.126    -1.137922    .1402776
                           |
                     _cons |  -2.354516   .0532507   -44.22   0.000    -2.458885   -2.250147
--------------------------------------------------------------------------------------------

.         
. nlcom exp((1.sameblock#c.m8 + 1.sameblock#c.m9 + 1.sameblock#c.m10 ///
>         + 1.sameblock#c.m11 + 1.sameblock#c.m12)/5) - 1

       _nl_1: exp((1.sameblock#c.m8 + 1.sameblock#c.m9 + 1.sameblock#c.m10         + 1.sameblock#c.m11 + 1.sameblock#c.m12)/5) - 1

------------------------------------------------------------------------------
      thefts | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
       _nl_1 |          0  (omitted)
------------------------------------------------------------------------------

.         
. lincom (1.sameblock#c.m8 + 1.sameblock#c.m9 + 1.sameblock#c.m10 ///
>         + 1.sameblock#c.m11 + 1.sameblock#c.m12)/5, coeflegend

 ( 1)  .2*[thefts]1.sameblock#c.m8 + .2*[thefts]1.sameblock#c.m9 + .2*[thefts]1.sameblock#c.m10 + .2*[thefts]1.sameblock#c.m11 + .2*[thefts]1.sameblock#c.m12 = 0

------------------------------------------------------------------------------
      thefts | Coefficient  Legend
-------------+----------------------------------------------------------------
         (1) |   -1.28483  _b[(1)]
------------------------------------------------------------------------------

. scalar avgatts = r(estimate)

. nlcom exp(avgatts) - 1

       _nl_1: exp(avgatts) - 1

------------------------------------------------------------------------------
      thefts | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
       _nl_1 |  -.7233024          .        .       .            .           .
------------------------------------------------------------------------------

[Leonardo Guizzetti]

I don't think there's anything wrong with your -poisson- model, Jeff. I can't see why you get any -nlcom- estimates at all as written, because if I try to reference coefficients without wrapping them in the -_b[]- syntax, I get the following error in Stata 16 and 17.

webuse dollhill3
poisson deaths i.smokes
nlcom (1.smokes)

Results in the error

expression (1.smokes) contains reference to X rather than _b[X]
r(198)

Can you try using the -nlcom- syntax I suggested in #6 and see if you still get non-sense results?

[Andrew Musau]

Jeff Wooldridge, I usually recommend using margins with -expression()- and -post- options in these cases. Example:

sysuse auto, clear
reg price mpg head
lincom mpg+head
margins, expression(_b[mpg]+_b[headroom]) post
nlcom _b[_cons] - 1

Res.:

. lincom mpg+head

 ( 1)  mpg + headroom = 0

------------------------------------------------------------------------------
       price |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         (1) |  -593.1271   427.0515    -1.39   0.169    -1444.644    258.3893
------------------------------------------------------------------------------

.
. margins, expression(_b[mpg]+_b[headroom]) post
Warning: expression() does not contain predict() or xb().
Warning: prediction constant over observations.

Predictive margins                              Number of obs     =         74
Model VCE    : OLS

Expression   : _b[mpg]+_b[headroom]

------------------------------------------------------------------------------
             |            Delta-method
             |     Margin   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
       _cons |  -593.1271   427.0515    -1.39   0.165    -1430.133    243.8784
------------------------------------------------------------------------------

.
. nlcom _b[_cons] - 1

       _nl_1:  _b[_cons] - 1

------------------------------------------------------------------------------
             |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
       _nl_1 |  -594.1271   427.0515    -1.39   0.164    -1431.133    242.8784
------------------------------------------------------------------------------

.

You may also be interested in the community-contributed command xlincom from SSC that can post results to e(b) and e(V).

[Jeff Wooldridge]

Many thanks again Leonardo. I got it to work by actually doing what you actually suggested rather than my own version of it. When I reference coefficients by _b[1.sameblock#c.m8] rather than just the variable name it works. I've probably run into this problem previously and forgot the solution. Age ....

Andrew: Thanks for your tip. It provides another way for me to check my results and I believe helps me solve a related problem.

[Leonardo Guizzetti]

Jeff, I'm happy it worked out. I too forget about these differences in -lincom- and -nlcom-.

Andrew, thanks for the additional ideas. The margins trick is nice to know.

[Jeff Wooldridge]

I'm still intrigued that I didn't get an error message with my incorrect usage of the nlcom command. My suggestion to Stata would be to allow the same syntax as lincom, although it was my shortcoming not to remember the proper syntax for nlcom.

[Joro Kolev]

This is a peculiar thread, Professor Wooldridge, because instead of asking about the problem you had "How to make my -nlcom- work/Why is not my -lncom- not working," you asked some other unrelated question " I'd like to insert the result from lincome into nlcom."

However it all turned out for the awesome, because we learnt lots of neat tricks. I did not know that one can use -margins- the way how Andrew showed, and although I am very familiar with -nlcom- and in principle I know of all of the features Leonardo showed, it has never come to my minds to avoid the "two stage problem" and carry forward the variance matrix of the estimate the way he did in #7.

As for why your Stata is not flagging the syntax error in -nlcom- for you, is hard to know without a reproducible example and without knowing which version of Stata you use.

Leonardo reported that he is unable to reproduce the problem on Stata 16 and 17, and I am not able to reproduce the problem on Stata 15:

.  qui reg price mpg head i1.foreign i1.foreign#c.mpg

. nlcom  exp(i1.foreign + i1.foreign#c.mpg)

expression (exp(i1.foreign + i1.foreign#c.mpg)) contains reference to X rather than _b[X]
r(198);

I'm still intrigued that I didn't get an error message with my incorrect usage of the nlcom command. My suggestion to Stata would be to allow the same syntax as lincom, although it was my shortcoming not to remember the proper syntax for nlcom.

[Jeff Wooldridge]

I should've mentioned I'm using Stata 17. My initial thought was that, because I wanted to obtain the standard error for a linear combination of coefficients followed by the exponential of those same coefficients, I should be able to use nlcom after lincom. I don't think that's such as stretch as it makes perfect sense if one remembers how the delta method works. Then, of course, I used the incorrect syntax for nlcom, but Stata didn't tell me that. Instead, it produced the wrong number and no standard error.

All I know is what I typed and what Stata gave. That it gives a different message for a different data set doesn't solve the puzzle for me.

Here's from my same example, with a simpler expression:

. nlcom exp(1.sameblock)

       _nl_1: exp(1.sameblock)

------------------------------------------------------------------------------
      thefts | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
       _nl_1 |          1          .        .       .            .           .
------------------------------------------------------------------------------

. nlcom exp(_b[1.sameblock])

       _nl_1: exp(_b[1.sameblock])

------------------------------------------------------------------------------
      thefts | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
       _nl_1 |   .4570367   .1232706     3.71   0.000     .2154307    .6986427
------------------------------------------------------------------------------