Manual calculation of marginal effects by plug-ing parameter estimates

[Petar Peric]

Dear all,
I would really appreciate it if anyone could help me out with a code for my calculations of marginal effects.

The data I am working on can be obtained only with another statistical program which isn't able of calculating marginal effects.
The model I am working with is Probit, which gives me the regression coefficients I would like to use in manual calculation of marginal effects in Stata. As I am not used to working with Stata this calculations seem very abstract to me. How do I get the marginal effects by manually plug-ing coefficients and means of predictors?

So far I have read page 28 from article - https://www3.nd.edu/~rwilliam/stats/Margins01.pdf , but it seems not sufficient and its for logit output.

I tried using the method described at - http://www.stata.com/support/faqs/st...-after-offset/ where it describes the calculation of marginal effects of two dichotomous variables by hand, but even that couldn't get me to replicate marginal effects from article - https://www3.nd.edu/~rwilliam/stats3/Margins02.pdf which I want to do for exercise before conducting work on my data.

I would really appreciate every or any suggestion you might have.
I thank you in advance.
Kind regards,
Petar

[Scott Merryman]

Do you means something like this:

. sysuse auto,clear
(1978 Automobile Data)

. qui probit foreign  mpg

. margins, atmeans

Adjusted predictions                            Number of obs     =         74
Model VCE    : OIM

Expression   : Pr(foreign), predict()
at           : mpg             =     21.2973 (mean)

------------------------------------------------------------------------------
             |            Delta-method
             |     Margin   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
       _cons |   .2777808   .0550856     5.04   0.000      .169815    .3857465
------------------------------------------------------------------------------

. qui sum mpg

. local m = r(mean)

. disp "myprediction = " normal(_b[_cons] +  `m'*_b[mpg])
myprediction = .27778078

Or perhaps this:

. margins, dydx(mpg) atmean

Conditional marginal effects                    Number of obs     =         74
Model VCE    : OIM

Expression   : Pr(foreign), predict()
dy/dx w.r.t. : mpg
at           : mpg             =     21.2973 (mean)

------------------------------------------------------------------------------
             |            Delta-method
             |      dy/dx   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         mpg |   .0322111   .0101675     3.17   0.002     .0122833     .052139
------------------------------------------------------------------------------

. disp "mymarginaleffect = " (_b[mpg])*normalden(_b[_cons] +  `m'*_b[mpg])
mymarginaleffect = .03221113

.

[Petar Peric]

Scott, thank You!
Those are the basics I had to know, yes.
Maybe you can help me also with resolving a problem of how to, also manually with only (probit) parameter estimates at my hand, obtain marginal effect which will show me the difference in the predicted probabilities for cases in one category of predictor relative to the reference category.
Let us say in case of predictor 'rank', how much more/less likely is for a subject with rank 1 to get into graduate school then rank 2, or 3, or 4.
I am referring to example at:
http://www.ats.ucla.edu/stat/stata/dae/binary.dta

[Scott Merryman]

Isn't that just the difference in the adjusted predicted values?

use http://www.ats.ucla.edu/stat/stata/dae/binary.dta, clear
qui probit admit gre gpa i.rank
qui sum gre
local mgre = r(mean)
qui sum gp
local mgpa = r(mean)
disp "Rank 2  vs Rank 1 = " normal(_b[_cons] +  `mgre'*_b[gre] + `mgpa'*_b[gpa]) - normal(_b[_cons] +  _b[2.rank] + `mgre'*_b[gre] + `mgpa'*_b[gpa])

[Petar Peric]

Thank You again Scott!
That is exactly what I needed to know.

[Maxence Morlet]

Dear Scott Merryman,

In #2, if the model were a logit, would the following be correct to calculate the marginal effect "manually" of the variable mpg:

logit foreign mpg
predict xb_logit, xb
g marg_eff_mpg_logit=(exp(xb_logit)/1+exp(xb_logit))^2*_b[mpg]

And then, in order to get only one number, an average partial effect, that does not vary with each unit, write

egen avg_marg_mpg_logit=mean(marg_eff_mpg_logit)

?

If mpg had been binary instead of continuous, could we just predict the outcome when mpg is one, then predict the outcome when mpg is 0, and then take the difference?

[Scott Merryman]

For logit model the conditional and average marginal effects can be calculated by hand:

. sysuse auto,clear
(1978 automobile data)

. qui logit foreign  mpg

. //Conditional marginal effect | mpg = mean
. qui sum mpg

. gen pr1 =  exp(_b[_cons]+ r(mean)*_b[mpg])/(1 + exp(_b[_cons]+ r(mean)*_b[mpg]))

. gen pr0 = 1- pr1

. gen pr1pr0 = pr1*pr0

. qui sum pr1pr0

. disp  "Conditional marginal effect = " r(mean)*_b[mpg]
Conditional marginal effect = .03175261

. margins, dydx(mpg) atmean 

Conditional marginal effects                                Number of obs = 74
Model VCE: OIM

Expression: Pr(foreign), predict()
dy/dx wrt:  mpg
At: mpg = 21.2973 (mean)

------------------------------------------------------------------------------
             |            Delta-method
             |      dy/dx   std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
         mpg |   .0317526   .0103945     3.05   0.002     .0113798    .0521254
------------------------------------------------------------------------------

. 
. //Average marginal effect
. gen pr_1 =  exp(_b[_cons]+ mpg*_b[mpg])/(1 + exp(_b[_cons]+ mpg*_b[mpg]))

. gen pr_0 = 1- pr_1

. gen pr_1pr_0 = pr_1*pr_0

. qui sum pr_1pr_0

. disp  "Average marginal effect = " r(mean)*_b[mpg]
Average marginal effect = .02811538

. margins, dydx(mpg) 

Average marginal effects                                    Number of obs = 74
Model VCE: OIM

Expression: Pr(foreign), predict()
dy/dx wrt:  mpg

------------------------------------------------------------------------------
             |            Delta-method
             |      dy/dx   std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
         mpg |   .0281154   .0071019     3.96   0.000     .0141959    .0420349
------------------------------------------------------------------------------

If mpg had been binary instead of continuous, could we just predict the outcome when mpg is one, then predict the outcome when mpg is 0, and then take the difference?

Yes, but keep in mind it is when all mpg values are one minus when all mpg values are zero:

. sysuse auto,clear
(1978 automobile data)

. replace mpg = (mpg < 20)
(74 real changes made)

. qui logit foreign  i.mpg

. margins, dydx(mpg) 

Conditional marginal effects                                Number of obs = 74
Model VCE: OIM

Expression: Pr(foreign), predict()
dy/dx wrt:  1.mpg

------------------------------------------------------------------------------
             |            Delta-method
             |      dy/dx   std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
       1.mpg |  -.2930403   .0990123    -2.96   0.003    -.4871008   -.0989797
------------------------------------------------------------------------------
Note: dy/dx for factor levels is the discrete change from the base level.

. replace mpg = 1
(39 real changes made)

. gen pr_1 =  exp(_b[_cons]+ mpg*_b[1.mpg])/(1 + exp(_b[_cons]+ mpg*_b[1.mpg]))

. replace mpg = 0
(74 real changes made)

. gen pr_0  =  exp(_b[_cons]+ mpg*_b[1.mpg])/(1 + exp(_b[_cons]+ mpg*_b[1.mpg]))

. gen delta  = pr_1 - pr_0

. qui sum delta

. 
. disp "Discrete change = " r(mean)
Discrete change = -.29304028

[Maxence Morlet]

Dear Scott Merryman,

Thank you very much for your prompt response. I used your reasoning in a Monte Carlo context, with a multivariate regression. Is the following then correct? I am interested in the marginal effect of being female on the outcome.

set seed 1
set obs 1000
    gen age = runiformint(18,65) 

    gen union_unbiased = rbinomial(1,0.3) 
    rename union_unbiased union_latent_unbiased
    g union_unbiased=union_latent_unbiased>0
    gen female_unbiased = rbinomial(1, 0.1) 
    rename female_unbiased female_latent_unbiased
    g female_unbiased=female_latent_unbiased>0

    g highschool_unbiased=rbinomial(1,0.8)>0

g y_latent_unbiased = 0.03 + 0.4*age - 0.02*age^2 + 3*highschool_unbiased -0.06*union_unbiased + 5*female_unbiased 

g y_logit_unbiased = (y_latent_unbiased+rlogistic())>0

g y_latent_unbiased1= 0.03 + 0.4*age - 0.02*age^2 + 3*highschool_unbiased -0.06*union_unbiased +5 //Supposing all individuals are female, 5*1 = 5. Correct?
g y_latent_unbiased0= 0.03 + 0.4*age - 0.02*age^2 + 3*highschool_unbiased -0.06*union_unbiased //Now no one is female
g y_logit_unbiased_fem1=exp(y_latent_unbiased1)/(1+exp(y_latent_unbiased1))
g y_logit_unbiased_fem0=exp(y_latent_unbiased0)/(1+exp(y_latent_unbiased0))

g true_marg_logit_unbiased=y_logit_unbiased_fem1 - y_logit_unbiased_fem0 //Finally, is this correct?

[FernandoRios]

Hi Maxence

Two points,
First, I should have sent you this in the earlier interactions we have:
https://blog.stata.com/2016/01/14/re...obit-or-logit/
https://blog.stata.com/2016/01/07/pr...k-your-weapon/
Second, I think getting marginal effects estimation can be done slightly different (for the sake of speed)

First. Create all your variables and unbiased latent
call it XB

You can estimate marginal effects for logit using this xb and the formulas you have above.
however, for the montecarlo analysis you can focus on creating variation through the error.

in other words, XB will be fixed across all simulations (so you don't need to reestimate the marginal effects)
the error will change every time.

Here an example:

clear
set obs 100
gen x1 = rchi2(2)-2
gen x2 = rchi2(2)-2
gen x3 = rnormal()>.75
gen y_latent = 1 + x1+x2+x3
gen mfx = logistic(1+x1+x2+1)-logistic(1+x1+x2)
gen y = .
set seed 10101
capture program drop simx
program simx, eclass
    sum mfx, meanonly
    matrix b=r(mean) 
    replace y =(y_latent+rlogistic())>0
    logit y x1 x2 i.x3
    margins, dydx(x3) post
    matrix b=b,_b[1.x3]
    probit y x1 x2 i.x3
    margins, dydx(x3) post
    matrix b=b,_b[1.x3]
    reg y x1 x2 i.x3
    matrix b=b,_b[1.x3]
    ereturn post b
end

simulate, reps(200): simx
sum


    Variable |        Obs        Mean    Std. dev.       Min        Max
-------------+---------------------------------------------------------
       _b_c1 |        200    .1116338           0   .1116338   .1116338
       _b_c2 |        200    .1181155    .0694874  -.0607635   .3355345
       _b_c3 |        200    .1179743     .069984  -.0641696   .3216946
       _b_c4 |        200    .1296684    .0772143  -.0573053   .3722865

Of course, you will need far more repetitions, but this should be a start

Fernando

[Maxence Morlet]

Our previous exchanges have really helped me, without them I wouldn't have gotten very far!

Thank you for pointing out the computational inefficiencies in my code in #8.

So just to be sure, one point you made in #9 was that the code in #8 is correct, but not efficient and may take unnecessarily long to run due to the re-estimation of marginal effects at each iteration?

Thank you very much again for all your help!

[FernandoRios]

I think its just a different way of thinking about the model.
If you randomize X's, your marginal effects will vary as well (due to sampling). However, we usually say lets consider X as fixed.
In that case, the only think changing are the errors, which is what I suggest here.
You could still do the analysis sampling everything, and your "true" marginal effect will just be different on each iteration.
In my example, the "true" marginal effect has zero variance. If you randomize everything, it will have Some variance (sampling variation)

Best wishes
Fernando