Re: st: ratio of marginal effects when using two -margins- commands

[jpitblado@stata.com (Jeff Pitblado, StataCorp LP)]

Mirko <mirko.moro@gmail.com> asks how to use -nlcom- with results from
multiple calls to -margins-:

> I am wondering whether is possible to obtain standard errors and
> confidence intervals of ratio of marginal effects when they are
> obtained by running two (or more) times the command -margins-.
> -margins- makes life a lot easier especially after estimating models
> with interaction terms or nonlinear models:
> 
> * simple case
> use http://www.stata-press.com/data/r11/margex
> logistic outcome sex##group age
> margins sex, post
> nlcom (risk_ratio: _b[1.sex] / _b[0.sex])
> 
> 
> However, there may be some cases in which one needs to compute two
> separate -margins- commands to obtain the marginal effects of
> interest. For example, after a model of Y on X, Z and the interaction
> term XZ:
> 
> Y = a + bX + bZ + bXZ + e,
> 
> I'd like to obtain the statistical significance of the ratio of
> 
> (marginal effects of variable X conditional on variable Z)/ (marginal
> effects of variable Z conditional on X) =
> 
> ME1/ME2
> 
> ****begin example****
> * -margins*
> 
> sysuse auto, clear
> set more off
> qui regress mpg foreign i.rep78##c.weight headroom
> * ME1
> margins, dydx(rep78) atmeans
> * ME2
> margins, dydx(weight) over(rep78)
> 
> ****end example****
> 
> Is there an easy way to get standard errors and confidence intervals
> of ratios of two -margins- ME1/ME2?
> 
> Or do I need to use -nlcom- like below?
> 
> 
> ****begin example****
> * ratio of marginal effects with -nlcom-
> sysuse auto, clear
> set more off
> qui tab rep78, gen(rep)
> forval i=2/5{
> 	qui gen rep`i'Xweight = rep`i'*weight
> }
> regress mpg foreign rep2-rep5 rep2Xweight-rep5Xweight weight headroom
> qui sum weight if e(sample)
> local meanw = r(mean)
> * ratio of ME1/ME2
> nlcom (_b[rep2] + _b[rep2Xweight]*`meanw')/(_b[weight] + _b[rep2Xweight])
> nlcom (_b[rep3] + _b[rep3Xweight]*`meanw')/(_b[weight] + _b[rep3Xweight])
> nlcom (_b[rep4] + _b[rep4Xweight]*`meanw')/(_b[weight] + _b[rep4Xweight])
> nlcom (_b[rep5] + _b[rep5Xweight]*`meanw')/(_b[weight] + _b[rep5Xweight])
> 
> ****end example****

The -margins- command specifically posts the Jacobian matrix in -r(Jacobian)-
so that results from different -margins- calls on the same model fit can be
combined.  Mirko just needs to be careful about equation names in the combined
results before -ereturn post-ing them.

We'll use Mirko's second example to illustrate how to combine the results
from two separate calls to -margins-.

First let's respecify the model using factor variables notation so that we can
use -margins- to get the numerator and denominator marginal effects:

. regress mpg for rep78##c.weight headroom
(output omitted)

The equivalent -nlcom- command to Mirko's is:

***** BEGIN:
. nlcom	(R2: (_b[2.r] + _b[2.r#w]*`meanw')/(_b[w] + _b[2.r#w]))	///
>	(R3: (_b[3.r] + _b[3.r#w]*`meanw')/(_b[w] + _b[3.r#w]))	///
>	(R4: (_b[4.r] + _b[4.r#w]*`meanw')/(_b[w] + _b[4.r#w]))	///
>	(R5: (_b[5.r] + _b[5.r#w]*`meanw')/(_b[w] + _b[5.r#w]))

          R2:  (_b[2.r] + _b[2.r#w]*3032)/(_b[w] + _b[2.r#w])
          R3:  (_b[3.r] + _b[3.r#w]*3032)/(_b[w] + _b[3.r#w])
          R4:  (_b[4.r] + _b[4.r#w]*3032)/(_b[w] + _b[4.r#w])
          R5:  (_b[5.r] + _b[5.r#w]*3032)/(_b[w] + _b[5.r#w])

------------------------------------------------------------------------------
         mpg |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
          R2 |  -33.06524   311.6278    -0.11   0.916    -657.0891    590.9586
          R3 |   166.0963   502.9187     0.33   0.742     -840.981    1173.174
          R4 |   76.42387   460.7441     0.17   0.869    -846.2003     999.048
          R5 |   363.8427   148.9601     2.44   0.018     65.55527    662.1302
------------------------------------------------------------------------------
***** END:

Now that we have the -regress- model using factor variables notation, we can
get the two sets of marginal effects results.  We'll need to grab the point
estimates from -r(b)- and corresponding Jacobian matrix from -r(Jacobian)-.
We can use the Jacobian matrix to reproduce the variance estimates by the
following matrix product

	r(V) = r(Jacobian)*e(V)*r(Jacobian)'

***** BEGIN:
. * reference contrasts on the margins of rep78, i.e. effects of factor rep78
. margins, dydx(rep78)
(output omitted)
. matrix b_num = r(b)
. matrix colna b_num = num:
. matrix J_num = r(Jacobian)
. matrix rowna J_num = num:

. * marginal effects of weight for each level of rep78
. margins, dydx(weight) over(rep78) 
(output omitted)
. matrix b_den = r(b)
. matrix colna b_den = den:
. matrix J_den = r(Jacobian)
. matrix rowna J_den = den:
***** END:

Notice that we added our own equation name to each set of results we pulled
from -margins-.  We used 'num' for the numerator results, and 'den' for the
denominator.

All we need to do now is post the combined results, then use -nlcom-:

***** BEGIN:
. * combine the results
. matrix b = b_num, b_den
. matrix J = J_num \ J_den
. matrix V = J*e(V)*J'
. ereturn post b V
. * check that the combined results match those from the original
. ereturn display
(output omitted)
***** END:

***** BEGIN:
. nlcom   (R2: [num]_b[2.r]/[den]_b[2.r]) ///
>         (R3: [num]_b[3.r]/[den]_b[3.r]) ///
>         (R4: [num]_b[4.r]/[den]_b[4.r]) ///
>         (R5: [num]_b[5.r]/[den]_b[5.r])

          R2:  [num]_b[2.r]/[den]_b[2.r]
          R3:  [num]_b[3.r]/[den]_b[3.r]
          R4:  [num]_b[4.r]/[den]_b[4.r]
          R5:  [num]_b[5.r]/[den]_b[5.r]

------------------------------------------------------------------------------
             |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
          R2 |  -33.06645   311.6233    -0.11   0.915    -643.8368    577.7039
          R3 |   166.0723   502.9125     0.33   0.741     -819.618    1151.763
          R4 |   76.40598   460.7396     0.17   0.868    -826.6271     979.439
          R5 |   363.8585   148.9596     2.44   0.015     71.90304     655.814
------------------------------------------------------------------------------
***** END:

Note that -nlcom- after -regress- reports 't' statistics, p-values, and CIs;
but -margins- and our combined results report 'z' statistics.  This is because
we didn't post the -e(df_r)- from the -regress- results in our combined
results.

--Jeff
jpitblado@stata.com
*
*   For searches and help try:
*   http://www.stata.com/help.cgi?search
*   http://www.stata.com/support/statalist/faq
*   http://www.ats.ucla.edu/stat/stata/