Dear all,
A while ago, I asked the Stata list subscribers about how to
save the results from -mfx-. Kit kindly suggested me the solution.
I can get the program to run and give me the results when I use
-mfx- *once* in the program. The following program is the simplified
version of my real one. I run this with the data set containing
variables y, x1,..,x6.
This program works fine. However, I would be happy to know the shorter
way to write the code for -r( )- in line 22-26 because in my real program
I have at least 20 variables.
--------------------------------
1.capture program drop simlognorm
2.program define simlognorm, rclass
3. version 9.0
4. gen a = x1sd*invnorm(uniform())+ x1mean
5. local depvar x2 x3 x4 x5 x6
6. logit y a `depvar', robust
7. mfx, predict(p) at(mean x5=0 x6=0)
8. ret scalar pry = e(Xmfx_y)
9. ret scalar pr2 = e(r2_p)
10. ret scalar chi2 = e(chi2)
11. mat dydxmfx = e(Xmfx_dydx)
12. forval i = 1/6 {
13. ret scalar dydxmfx`i' = dydxmfx[1,`i']
14. }
15. mat sedyx = e(Xmfx_se_dydx)
16. forval i = 1/6 {
17. ret scalar sedyx`i' = sedyx[1,`i']
18. }
19. eret clear
20. drop a
21.end
22.simulate pry = r(pry) pr2 = r(pr2) chi2 = r(chi2) ///
23.dydxmfx1 = r(dydxmfx1) dydxmfx2 = r(dydxmfx2) dydxmfx3 = r(dydxmfx3) ///
24.dydxmfx4 = r(dydxmfx4) dydxmfx5 = r(dydxmfx5) dydxmfx6 = r(dydxmfx6) ///
25.sedyx1 = r(sedyx1) sedyx2 = r(sedyx2) sedyx3 = r(sedyx3) ///
26.sedyx4 = r(sedyx4) sedyx5 = r(sedyx5) sedyx6 = r(sedyx6), reps(1000) ///
27.saving("C:\data\simlognorm.dta", replace): simlognorm
--------------------------------
Here is the serious problem. I have to use -mfx- several times in the program
because I want to find the effect of variable -a- at its different values than
just its mean. The code is similar as I just add in one more section between the line 18 and the line 19 in the above program.
--------------------------------
mfx, predict(p) at(mean a=9 x5=0 x6=0)
ret scalar pry9 = e(Xmfx_y)
mat dydxmfx = e(Xmfx_dydx)
ret scalar dydxmfxp9 = dydxmfx[1,1]
mat sedyx = e(Xmfx_se_dydx)
ret scalar sedyxp9 = sedyx[1,1]
--------------------------------
And of course, I add in between line 26-27 with
--------------------------------
dydxmfxp9 = r(dydxmfxp9) sedyxp9 = r(sedyxp9) pry9 = r(pry9) ///
--------------------------------
After I ran this modified program, Stata gave me something like this
command: simlognorm
pry: r(pry)
pr2: r(pr2)
chi2: r(chi2)
df: r(df)
dydxmfx1: r(dydxmfx1)
dydxmfx2: r(dydxmfx2)
dydxmfx3: r(dydxmfx3)
dydxmfx4: r(dydxmfx4)
dydxmfx5: r(dydxmfx5)
dydxmfx6: r(dydxmfx6)
sedyx1: r(sedyx1)
sedyx2: r(sedyx2)
sedyx3: r(sedyx3)
sedyx4: r(sedyx4)
sedyx5: r(sedyx5)
sedyx6: r(sedyx6)
dydxmfxp9: r(dydxmfxp9)
sedyxp9: r(sedyxp9)
pry9: r(pry9)
Simulations (1000)
----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5
.................................................. 50
.................................................. 100
.................................................. 150
.................................................. 200
.................................................. 250
...this error is most likely due to clear being used within: sim
> lognorm
--Break--
r(1);
end of do-file
--Break--
r(1);
++++++++++++++++++++++++++++++++++
My questions are:
(1) What are the problems in my program? I am sorry that I did not
have the log file with the trace this time. But if it is necessary
I will provided this later.
(2) Before I run the complete simulation, I usually test it with -5-
repetitions so that I can see the errors faster. However, why would in
this case the error start to show after 250 reps already?
(3) Any suggestion about shortenning the codes line 23-26 would be greatly
appreciated.
Would anybody please help?
Thank you very much,
Anupit
*
* 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/
