'Manual' dy/dx using predict?

[Lawrence McKay]

Dear all,

I have a question which could not be resolved through the manuals. I am using a program called runmlwin to run a multilevel regression for an ordered categorical outcome. The key quantity of interest is the marginal effect of national economic growth (q_q_growth) on perceptions of the national economy (econGenRetroW5x_rev) at different levels of local economic prosperity (z_wages_year_no_mis_mc); i.e., comparing low-wage to high-wage areas. The runmlwin software does not post results that are manipulable via margins and thus, I am using 'predict' to do the job manually. The regression results being used are in the attached .docx file.

I first thought that I could return results for a 0.1 unit change in q_q_growth by fixing parameters at, for example, 0.2 and 0.1, and taking the difference. So my predictions look like:

      predictnl p12_wages_mar = (invlogit([FP2]cons_2 + [FP5]q_q_growth_1234*.2 + [FP5]z_wages_year_no_mis_mc_1234*z_wages_year_no_mis_mc + [FP5]growth_wages_prod_1234*z_wages_year_no_mis_mc*.2 + [FP5]dummy_wave4_1234)) - ///
(invlogit([FP2]cons_2 + [FP5]q_q_growth_1234*.1 + [FP5]z_wages_year_no_mis_mc_1234*z_wages_year_no_mis_mc + [FP5]growth_wages_prod_1234*z_wages_year_no_mis_mc*.1 + [FP5]dummy_wave4_1234)) if _est_egr_nd_pql1==1, ci(p12_wages_mar_lb p12_wages_mar_ub)

But I found the results were sensitive to where I placed the tenth-unit gap, for example, compare the results of a prediction that fixes q_q_growth at 1 and 0.9 using the following code:

         predictnl p12_wages_mar = (invlogit([FP2]cons_2 + [FP5]q_q_growth_1234*1 + [FP5]z_wages_year_no_mis_mc_1234*z_wages_year_no_mis_mc + [FP5]growth_wages_prod_1234*z_wages_year_no_mis_mc*1 + [FP5]dummy_wave4_1234)) - ///
(invlogit([FP2]cons_2 + [FP5]q_q_growth_1234*.9 + [FP5]z_wages_year_no_mis_mc_1234*z_wages_year_no_mis_mc + [FP5]growth_wages_prod_1234*z_wages_year_no_mis_mc*.9 + [FP5]dummy_wave4_1234)) if _est_egr_nd_pql1==1, ci(p12_wages_mar_lb p12_wages_mar_ub)

(Note that by using cons_2, I'm actually getting STATA to return the probability that R chooses either category 1 (Economy got a lot better) or category 2 (got a bit better) - though I don't think that's too important here.)



I'm not sure why the variation in results would occur, but whatever the case, it isn't ideal for me.

If I had got my predictions from regress, mixed etc, I would simply be able to do something like this, which would return the marginal effect of a one-unit change in q_q_growth on the probability of choosing each response category at a given level of local wages.

ologit econGenRetroW5x c.q_q_growth##c.z_wages_year_no_mis_mc
margins, dydx(q_q_growth) at(z_wages_year_no_mis_mc=(-2 2))

Is there something I can do using predict that would similarly estimate a marginal effect in this way, but was valid for the entire range of q_q_growth? After all, I understand from the manual that everything margins does is fundamentally based on STATA using 'predict' behind the scenes.

Thanks for your time.

Here is a data sample, in case it helps understand the structure, and the initial regression input, but this may be obscure for those not versed in runmlwin.

I hope I've given enough information to understand the issue, but please let me know if I can clarify anything further.

* Example generated by -dataex-. For more info, type help dataex
clear
input float n long id int panoW double q_q_growth float(z_wages_year_no_mis_mc growth_wages_prod) byte(econGenRetroW5x_rev dummy_wave4)
1522050 297 1  0 .19215575          0 . 0
1522052 297 1 .3 -.8778121 -.26334363 4 0
1522054 297 1 .3 -.8778121 -.26334363 . 0
1522059 297 1 .1 .19215575 .019215574 2 0
1522060 297 1  0 .19215575          0 2 0
end
label values econGenRetroW5x_rev reveconGenRetroW5x
label def reveconGenRetroW5x 2 "Got a little better", modify
label def reveconGenRetroW5x 4 "Got a little worse", modify

runmlwin econGenRetroW5x_rev cons (q_q_growth z_wages_year_no_mis_mc growth_wages_prod dummy_wave2 dummy_wave3 dummy_wave4 dummy_wave6 dummy_wave7 dummy_wave8 dummy_wave10 dummy_wave11 dummy_wave13 dummy_wave14 dummy_wave15 dummy_wave16 dummy_wave17, contrast(1/4)) if estsample_wages_wd_noeu==1, ///
level3(panoW: (cons q_q_growth, contrast(1/4))) ///
level2(id: (cons, contrast(1/4))) ///
level1(n) ///
discrete(distribution(multinomial) link(ologit) denominator(cons) basecategory(5) pql1) ///
mlwinpath(C:\Program Files\MLwiN v3.05\mlwin.exe)
eststo egr_nd_pql1

[Daniel Schaefer]

You are using an ordered logistic regression to predict a probability, right? On the probability scale, your line of best fit becomes a sigmoid curve. It also seems like you have some random effects in your multilevel model, because you expect the effect size of q_q_growth will depend on levels of gross pay. So, effectively, you are trying to output the slope of the line for q_q_growth at different levels of gross pay, and you are doing this by taking a difference in the predicted outcome. The problem is that the effect size of q_q_growth is not constant even given a particular level of gross pay because there are other nonlinearities in the size of the effect for q_q_growth, particularly as a result of the sigmoid curve.

Ideally what you'd do here is take the partial derivative of your regression function with respect to q_q_growth, then plug in meaningful values of your other variables to calculate predicted marginal effects. Fortunately, -dydx- is its own command. I believe you can calculate the predicted probabilities of interest, then rather than taking the difference, plug the predicted probabilities into dydx as y and q_q_growth in as x to get your marginal effects.

help dydx

*Not holding q_q_growth constant in the predicted probabilities:
predictnl p12_wages_mar = ...
dydx p12_wages_mar q_q_growth, generate(marginal_effect)

[Lawrence McKay]

Dear Daniel,

Thank you for your quick response. Your explanation of the underlying issue makes sense.

I tried implementing this solution using the following code. As you suggest, I don't hold q_q_growth constant.

          predictnl p12_wages_mar = (invlogit([FP2]cons_2 + [FP5]q_q_growth_1234*q_q_growth + [FP5]z_wages_year_no_mis_mc_1234*z_wages_year_no_mis_mc + [FP5]growth_wages_prod_1234*growth_wages_prod + [FP5]dummy_wave4_1234))
         dydx p12_wages_mar q_q_growth, generate(meffect_wages_nocon)

meffect_wages_nocon is 8 observations (1 per value of q_q_growth), And if I plot this against wages as above, using this code:

 graph twoway ///
(qfit meffect_wages_nocon z_wages_year_no_mis_mc, sort), ///
   ytitle("Marginal effect of growth???", margin(medsmall)) xtitle("Median annual gross pay in constituency, standardised", margin(medsmall)) ///
   legend(off)

I get the following graph. I think it's right that the meffect of 1-unit change in growth (ie 0% vs 1%) is something like .2 on average, but this is not showing the expected variation with wages and, for some reason, only picking up the low pay constituencies.



Do you have any further suggestions? Furthermore, is there any way to obtain the standard errors on the marginal effect? (as I don't believe the dydx function is capable of providing them).

(Also, to clarify, there are random intercepts for id and panoW (parliamentary constituency), and a random slope for q_q_growth at the constituency level.)

Best
Lawrence

[Daniel Schaefer]

Hi Lawrence,

It looks like simulating the dydx option is a little more complicated than I implied in #2. To figure this out, I created a toy example that I think serves as a simplified version of what you want to do:

clear
set obs 1000
set seed 587932
generate x1 = rnormal()
generate x2 = rnormal()
generate linear_model = (3 * x1) + (5 * x2) + (2 * x1 * x2) + rnormal()
generate probability = invlogit(linear_model)
generate yhat = probability > runiform()
sum yhat probability x1 x2
logit yhat c.x1##c.x2

So yhat is a binary outcome and there is an interaction between x1 and x2. Ultimately you want to be able to do something like this:

margins, at(x1=(-3(1)3)) dydx(x2)

I've come up with an algorithm that does something almost identical to margins with the dydx option below.

quietly capture frame drop margins_frame
* Create a new frame to store the results
frame create margins_frame at dy_dx
* Foreach value of margins "at":
forv value = -3(1)3{
    * predict values for y given a known value of x1.
    predictnl double y_pred = invlogit(_b[_cons] + _b[x1]*`value' + _b[x2]*x2 + _b[c.x1#c.x2]*`value'*x2)
    * Calculate the derrivative at each x, y point.
    dydx y_pred x2, gen(marginal_effect) double
    * calculate the mean marginal effect across all x, y points.
    quietly sum marginal_effect, meanonly
    * add the results to the margins_frame
    frame post margins_frame (`value') (r(mean))
    * prepare for next loop iteration
    drop  marginal_effect y_pred
}

I say "almost" identical because if you compare the results of this algorithm to the output of the margins command, you will see close, but slightly different values. I attribute these differences to rounding error after floating point arithmetic.

margins, at(x1=(-3(1)3)) dydx(x2)
frame change margins_frame
list, clean noobs
frame change default

------------------------------------------------------------------------------
             |            Delta-method
             |      dy/dx   std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
x2           |
         _at |
          1  |  -.0006289   .0009242    -0.68   0.496    -.0024403    .0011824
          2  |   .0052766   .0046484     1.14   0.256     -.003834    .0143872
          3  |    .223464   .0168995    13.22   0.000     .1903416    .2565865
          4  |   .3751761   .0045549    82.37   0.000     .3662486    .3841035
          5  |   .3503687   .0050618    69.22   0.000     .3404477    .3602897
          6  |   .3133302   .0146362    21.41   0.000     .2846438    .3420166
          7  |   .2721098   .0197629    13.77   0.000     .2333753    .3108443
------------------------------------------------------------------------------

. list, clean noobs

    at       dy_dx  
    -3   -.0006289  
    -2    .0052756  
    -1    .2234609  
     0    .3751761  
     1    .3503687  
     2    .3133302  
     3    .2721098

Notice as well, you can generate linear predictions of the models by removing the invlogit() function from my algorithm and adding the predict(xb) option to the margins command. This will again result in very close predictions with very small differences due to rounding error. So if you want to implement the margins functionality with the dydx option for your model, this code should help you get started.

As for acquiring the standard errors, this should be possible in principle, but that is a whole other rabbit hole to go down. If I were you, I would start by extracting the variance and covariance matrixes after estimation, do some research to find the correct formula for the standard error of a partial effect under your model, then calculate the standard error with Stata matrixes from there. Keep in mind, you are asking to reproduce some sophisticated features of Stata at a low level, which is no easy task.

Edit: I inadvertently pasted in some incorrect results when I first posted. I've fixed the results in edits.

[Lawrence McKay]

Hi Daniel,

Thank you very much for that answer. I managed to apply the solution to my code and that, indeed, gets me the expected increase in the marginal effect of X on Y as Z increases.

It seems the standard errors are indeed a stumbling block. (And I appreciate you may want to bow out if it is beyond what you're comfortable with.) However, I will post just in case.

I was trying to replicate the approach from Jeff Pitblado here, but what I am looking for is three steps abstracted from the author's example:
1) Estimate the SE for a continuous by continuous interaction.
2) Estimate the SE for different values of the marginal effect
3) Apply this to an ologit (not a logit) model.

[I can successfully follow the steps set out to correctly estimate the AME and its SE of a continuous predictor with no interaction - i.e. the estimates match those produced by margins.]

Again, you've been a great help already so there is no obligation.

Best,
Lawrence

[Lawrence McKay]

Indeed, if any other forum users feel they may have a solution, please chime in

[Lawrence McKay]

I've actually had some success now, I'll post a quick follow-up next week for those who may be interested.