Page 1 of 1
Average Marginal Effects via runmlwin
Posted: Wed Feb 12, 2014 10:04 am
by leap
Dear everyone,
My question is the following:
How can I get average marginal effects and theirs standard errors from MLwiN via runmlwin?
A previous post on how to get predictions from MLwiN via runmlwin has already been answered. The post is extremely useful to obtain marginal effects at average and/or specific values of the covariates (
http://www.cmm.bristol.ac.uk/forum/view ... a2112#p452).
Using a similar method than the one explained in the post above, I can easily obtained the average marginal effects but not the standard errors. Without being able to use of the command margins, my understanding is that I would need to calculate the Jacobian matrix. However, this is not so straightforward in a software such as STATA, which is not meant for mathematics but applied statistics.
I have been struggling with this issue and still cannot find a solution, so any help would be gladly appreciated.
Thank you,
Léa Pessin
Re: Average Marginal Effects via runmlwin
Posted: Wed Feb 12, 2014 11:49 am
by GeorgeLeckie
Hi Léa,
Calculating the AME is not so hard (see below). But I don't think there is any easy way to calculate the standard error of this average.
Sorry not to be of more help
Best wishes
George
Code: Select all
use "http://www.bristol.ac.uk/cmm/media/runmlwin/tutorial.dta", clear
generate pass = (normexam>0)
melogit pass standlrt i.girl
margins, dydx(girl)
runmlwin pass cons standlrt girl, ///
level1(student:) ///
discrete(distribution(binomial) link(logit) denom(cons) pql2) ///
nopause
predictnl ame = ///
invlogit(_b[cons] + _b[standlrt]*standlrt + _b[girl]*1) ///
- invlogit(_b[cons] + _b[standlrt]*standlrt + _b[girl]*0)
tabstat ame, stat(mean)
with output...
Code: Select all
. use "http://www.bristol.ac.uk/cmm/media/runmlwin/tutorial.dta", clear
.
. generate pass = (normexam>0)
.
. melogit pass standlrt i.girl
Iteration 0: log likelihood = -2274.4478
Iteration 1: log likelihood = -2271.5217
Iteration 2: log likelihood = -2271.5149
Iteration 3: log likelihood = -2271.5149
Logistic regression Number of obs = 4059
Wald chi2(2) = 734.81
Log likelihood = -2271.5149 Prob > chi2 = 0.0000
------------------------------------------------------------------------------
pass | Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
standlrt | 1.28272 .0475777 26.96 0.000 1.18947 1.375971
1.girl | .2207443 .0737735 2.99 0.003 .0761509 .3653376
_cons | -.0851353 .0575463 -1.48 0.139 -.1979239 .0276533
------------------------------------------------------------------------------
.
. margins, dydx(girl)
Average marginal effects Number of obs = 4059
Model VCE : OIM
Expression : Predicted mean, predict()
dy/dx w.r.t. : 1.girl
------------------------------------------------------------------------------
| Delta-method
| dy/dx Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
1.girl | .042072 .0140612 2.99 0.003 .0145126 .0696314
------------------------------------------------------------------------------
Note: dy/dx for factor levels is the discrete change from the base level.
.
. runmlwin pass cons standlrt girl, ///
> level1(student:) ///
> discrete(distribution(binomial) link(logit) denom(cons) pql2) ///
> nopause
MLwiN 2.29 multilevel model Number of obs = 4059
Binomial logit response model
Estimation algorithm: IGLS, PQL2
Run time (seconds) = 1.39
Number of iterations = 5
------------------------------------------------------------------------------
pass | Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
cons | -.0851353 .0575431 -1.48 0.139 -.1979177 .0276472
standlrt | 1.28272 .0475705 26.96 0.000 1.189484 1.375956
girl | .2207442 .0737696 2.99 0.003 .0761584 .3653299
------------------------------------------------------------------------------
.
. predictnl ame = ///
> invlogit(_b[cons] + _b[standlrt]*standlrt + _b[girl]*1) ///
> - invlogit(_b[cons] + _b[standlrt]*standlrt + _b[girl]*0)
.
. tabstat ame, stat(mean)
variable | mean
-------------+----------
ame | .042072
------------------------
Best wishes
George
Re: Average Marginal Effects via runmlwin
Posted: Wed Feb 12, 2014 11:59 am
by GeorgeLeckie
Hi Léa,
Just had another thought.
If you can specify the same model using one of Stata's own multilevel modelling commands, but are avoiding Stata because it is computationally slow with your particular model and data, then you could use runmlwin to obtain the point estimates then plug them into the Stata multilevel modelling command as starting values. Then fit the model. The Stata multilevel modelling command should take far fewer iterations than normal and so converge far more quickly. So really this is just away of giving sensible starting values to Stata. You can then of course use the margins command however you want.
Example given below
George
Commands:
Code: Select all
* Load the data
use "http://www.bristol.ac.uk/cmm/media/runmlwin/tutorial.dta", clear
* Generate a binary response
generate pass = (normexam>0)
* Fit the model in runmlwin
runmlwin pass cons standlrt girl, ///
level1(student:) ///
discrete(distribution(binomial) link(logit) denom(cons) pql2) ///
nopause
* Save the vector of runmlwin point estimates
matrix b1 = e(b)
* Specify model using melogit, but do not estiamte iit
melogit pass standlrt i.girl, noestimate
* Save the vector of starting values
matrix b = e(b)
* List the runmlwin point estimates
matrix list b1
* List the melogit starting values
matrix list b
* Replace the melogit starting values with the runmlwin point estimates
matrix b[1,1] = b1[1,2]
matrix b[1,2] = 0
matrix b[1,3] = b1[1,3]
matrix b[1,4] = b1[1,1]
* List the updated melogit starting values
matrix list b
* Note model fits in one iteration, because we provided starting values which were effectively the estimates
melogit pass standlrt i.girl, from(b)
* Note that had we fitted the model with no starting values, it would have taken three iterations
melogit pass standlrt i.girl
* Compute the average marginal effect
margins, dydx(girl)
Output:
Code: Select all
. * Load the data
. use "http://www.bristol.ac.uk/cmm/media/runmlwin/tutorial.dta", clear
.
. * Generate a binary response
. generate pass = (normexam>0)
.
. * Fit the model in runmlwin
. runmlwin pass cons standlrt girl, ///
> level1(student:) ///
> discrete(distribution(binomial) link(logit) denom(cons) pql2) ///
> nopause
MLwiN 2.29 multilevel model Number of obs = 4059
Binomial logit response model
Estimation algorithm: IGLS, PQL2
Run time (seconds) = 1.48
Number of iterations = 5
------------------------------------------------------------------------------
pass | Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
cons | -.0851353 .0575431 -1.48 0.139 -.1979177 .0276472
standlrt | 1.28272 .0475705 26.96 0.000 1.189484 1.375956
girl | .2207442 .0737696 2.99 0.003 .0761584 .3653299
------------------------------------------------------------------------------
.
. * Save the vector of runmlwin point estimates
. matrix b1 = e(b)
.
. * Specify model using melogit, but do not estiamte iit
. melogit pass standlrt i.girl, noestimate
Posting starting values:
Logistic regression Number of obs = 4059
( 1) [pass]0b.girl = 0
------------------------------------------------------------------------------
pass | Coef. Legend
-------------+----------------------------------------------------------------
standlrt | 1.170814 _b[standlrt]
1.girl | .1997406 _b[1.girl]
_cons | -.0626776 _b[_cons]
------------------------------------------------------------------------------
Note: The above coefficient values are starting values and not the result of a fully fitted model.
.
. * Save the vector of starting values
. matrix b = e(b)
.
. * List the runmlwin point estimates
. matrix list b1
b1[1,4]
FP1: FP1: FP1: OD:
cons standlrt girl bcons_1
y1 -.08513525 1.2827201 .22074416 1
.
. * List the melogit starting values
. matrix list b
b[1,4]
pass: pass: pass: pass:
0b. 1.
standlrt girl girl _cons
y1 1.1708138 0 .19974063 -.06267763
.
. * Replace the melogit starting values with the runmlwin point estimates
. matrix b[1,1] = b1[1,2]
. matrix b[1,2] = 0
. matrix b[1,3] = b1[1,3]
. matrix b[1,4] = b1[1,1]
.
. * List the updated melogit starting values
. matrix list b
b[1,4]
pass: pass: pass: pass:
0b. 1.
standlrt girl girl _cons
y1 1.2827201 0 .22074416 -.08513525
.
. * Note model fits in one iteration, because we provided starting values which were effectively the estimates
. melogit pass standlrt i.girl, from(b)
Iteration 0: log likelihood = -2271.5149
Iteration 1: log likelihood = -2271.5149
Logistic regression Number of obs = 4059
Wald chi2(2) = 734.81
Log likelihood = -2271.5149 Prob > chi2 = 0.0000
------------------------------------------------------------------------------
pass | Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
standlrt | 1.28272 .0475777 26.96 0.000 1.18947 1.375971
1.girl | .2207443 .0737735 2.99 0.003 .0761509 .3653376
_cons | -.0851353 .0575463 -1.48 0.139 -.1979239 .0276533
------------------------------------------------------------------------------
.
. * Note that had we fitted the model with no starting values, it would have taken three iterations
. melogit pass standlrt i.girl
Iteration 0: log likelihood = -2274.4478
Iteration 1: log likelihood = -2271.5217
Iteration 2: log likelihood = -2271.5149
Iteration 3: log likelihood = -2271.5149
Logistic regression Number of obs = 4059
Wald chi2(2) = 734.81
Log likelihood = -2271.5149 Prob > chi2 = 0.0000
------------------------------------------------------------------------------
pass | Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
standlrt | 1.28272 .0475777 26.96 0.000 1.18947 1.375971
1.girl | .2207443 .0737735 2.99 0.003 .0761509 .3653376
_cons | -.0851353 .0575463 -1.48 0.139 -.1979239 .0276533
------------------------------------------------------------------------------
.
. * Compute the average marginal effect
. margins, dydx(girl)
Average marginal effects Number of obs = 4059
Model VCE : OIM
Expression : Predicted mean, predict()
dy/dx w.r.t. : 1.girl
------------------------------------------------------------------------------
| Delta-method
| dy/dx Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
1.girl | .042072 .0140612 2.99 0.003 .0145126 .0696314
------------------------------------------------------------------------------
Note: dy/dx for factor levels is the discrete change from the base level.
Re: Average Marginal Effects via runmlwin
Posted: Wed Feb 12, 2014 12:54 pm
by leap
Hello George,
Thank you for your prompt and detailed answer.
Unfortunately (or fortunately thanks to the existence of MLwiN), I am fully using the functionality of MLwiN (MCMC estimation, etc.) so I cannot use STATA as an alternative.
Is there any hope that runmlwin will one day be compatible with margins? It is really a shame because margins and marginsplot really offer some nice ways to illustrate regression output.
Thank you again,
Léa
Re: Average Marginal Effects via runmlwin
Posted: Wed Feb 12, 2014 1:11 pm
by GeorgeLeckie
Hi Léa,
I'm afraid we are not planning to make runmlwin compatible with margins. Not because it would not be useful, but because it would require a lot of programming.
If you are using MCMC estimation then things might actually be easier...
(1) Run the model for 1000 iterations and the save the chains as data. Or run it for 10000 iterations and set thinning to 10.
(2) Go back to the original data and expand the estimation sample by 1000 and then merge in the MCMC chains
(3) Calculate the AME in the usual way, but do it separately at each iteration of the MCMC chain
You will get 1000 values for the AME. The mean is your point estimate for the AME, the sd is your standard error.
Best wishes
George
variate normal.