Am trying to save some mcmc statistics in a table output from stata…and it is not quite working.
Am using estout after using mcmcsum.
I am running this code and I am getting the right ess(effective sample size) which I believe results from the mcmcsum but do not get any values for the quantiles of the chain. The other statistics (b,se,ci_l,ci_u) are fine too.
quietly runmlwin V1 cons, level2 (clid:cons) level1 (serialno:) discrete(distribution(binomial) link(logit) denominator(cons) pql2) nopause
quietly runmlwin V1 cons, level2 (clid:cons) level1 (serialno:) discrete(distribution(binomial) link(logit) denominator(cons)) mcmc(burnin(500) chain(10000)) initsprevious nopause
estimates store model
mcmcsum
estout model, cells(b se ci_l ci_u ess p50 p97_5 p2_5) stats(N), using "G:\rv1g09May9\paper 2\output5.xls"
Can anyone spot what I am doing wrong? Or suggest an alternative?
Thanks
saving statistics from mcmcsum
-
GeorgeLeckie
- Site Admin
- Posts: 432
- Joined: Fri Apr 01, 2011 2:14 pm
Re: saving statistics from mcmcsum
Hi,
This is an interesting question.
The reason why this does not work is becayse mcmcsum is not an estimation command (it is a general command) and so it doesn't post the estimation results.
estout only accesses estimation results and so the stats which are displayed come from your final runmlwin call and not the subsequent mcmcsum call.
We will have a think about how to improve mcmcsum to do what you are after
Best wishes
George
This is an interesting question.
The reason why this does not work is becayse mcmcsum is not an estimation command (it is a general command) and so it doesn't post the estimation results.
estout only accesses estimation results and so the stats which are displayed come from your final runmlwin call and not the subsequent mcmcsum call.
We will have a think about how to improve mcmcsum to do what you are after
Best wishes
George
-
GeorgeLeckie
- Site Admin
- Posts: 432
- Joined: Fri Apr 01, 2011 2:14 pm
Re: saving statistics from mcmcsum
Whoops,
I fired off the last email too quickly!
You can do exactly what you want
The quantiles of the mcmc chains are stored in e(quantiles). You can see this by typing the following after fitting your model
So to get what you want, type something like the following
Note, that you could also substitue lb and ub for p2_5 and p97_5 to get the lower and upper bounds for the 95% credible intervals.
Note that whether you run the mcmcsum command or not makes no difference. All of these stats are stored when you run the runmlwin command,
Best wishes
George
I fired off the last email too quickly!
You can do exactly what you want
The quantiles of the mcmc chains are stored in e(quantiles). You can see this by typing the following after fitting your model
Code: Select all
matrix list e(quantiles)Code: Select all
estout model, cells(b se ci_l ci_u ess quantiles[2] quantiles[5] quantiles[8]) stats(N), using "G:\rv1g09May9\paper 2\output5.xls"Note that whether you run the mcmcsum command or not makes no difference. All of these stats are stored when you run the runmlwin command,
Best wishes
George
Re: saving statistics from mcmcsum
thanks a lot for your immediate reply...managed to make that work now. thanks
Re: saving statistics from mcmcsum
Adding to the above post...is it also possible to save the DIC statistic and the variance covariance matrix of the parameter estimates using estout?
Rebecca
Rebecca
-
GeorgeLeckie
- Site Admin
- Posts: 432
- Joined: Fri Apr 01, 2011 2:14 pm
Re: saving statistics from mcmcsum
Hi Rebecca,
Anything which appears in ereturn list should be accessible by estout
So the following commands will present the parameter estimats, standard errors, sampling variances and DIC statistic
And here is the associated log
I hope this helps
George
Anything which appears in ereturn list should be accessible by estout
So the following commands will present the parameter estimats, standard errors, sampling variances and DIC statistic
Code: Select all
use http://www.bristol.ac.uk/cmm/media/runmlwin/tutorial, clear
runmlwin normexam cons standlrt, level2 (school: cons standlrt) level1 (student: cons) nopause
runmlwin normexam cons standlrt, level2 (school: cons standlrt) level1 (student: cons) mcmc(on) initsprevious nopause
ereturn list
estout, cells("b se var") stats(dic)Code: Select all
. use http://www.bristol.ac.uk/cmm/media/runmlwin/tutorial, clear
. runmlwin normexam cons standlrt, level2 (school: cons standlrt) level1 (student: cons) nopause
MLwiN 2.25 multilevel model Number of obs = 4059
Normal response model
Estimation algorithm: IGLS
-----------------------------------------------------------
| No. of Observations per Group
Level Variable | Groups Minimum Average Maximum
----------------+------------------------------------------
school | 65 2 62.4 198
-----------------------------------------------------------
Run time (seconds) = 1.62
Number of iterations = 4
Log likelihood = -4658.4351
Deviance = 9316.8701
------------------------------------------------------------------------------
normexam | Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
cons | -.0115052 .039783 -0.29 0.772 -.0894784 .066468
standlrt | .5567305 .019937 27.92 0.000 .5176548 .5958063
------------------------------------------------------------------------------
------------------------------------------------------------------------------
Random-effects Parameters | Estimate Std. Err. [95% Conf. Interval]
-----------------------------+------------------------------------------------
Level 2: school |
var(cons) | .0904446 .017924 .0553143 .1255749
cov(cons,standlrt) | .0180414 .0067229 .0048649 .031218
var(standlrt) | .0145361 .0044139 .0058851 .0231872
-----------------------------+------------------------------------------------
Level 1: student |
var(cons) | .5536575 .0124818 .5291936 .5781214
------------------------------------------------------------------------------
. runmlwin normexam cons standlrt, level2 (school: cons standlrt) level1 (student: cons) mcmc(on) initsprevious nopause
MLwiN 2.25 multilevel model Number of obs = 4059
Normal response model
Estimation algorithm: MCMC
-----------------------------------------------------------
| No. of Observations per Group
Level Variable | Groups Minimum Average Maximum
----------------+------------------------------------------
school | 65 2 62.4 198
-----------------------------------------------------------
Burnin = 500
Chain = 5000
Thinning = 1
Run time (seconds) = 9.78
Deviance (dbar) = 9122.99
Deviance (thetabar) = 9031.32
Effective no. of pars (pd) = 91.67
Bayesian DIC = 9214.65
------------------------------------------------------------------------------
normexam | Mean Std. Dev. ESS P [95% Cred. Interval]
-------------+----------------------------------------------------------------
cons | -.0056898 .0392222 281 0.432 -.0804512 .0736911
standlrt | .5582691 .0204706 806 0.000 .5176237 .5985798
------------------------------------------------------------------------------
------------------------------------------------------------------------------
Random-effects Parameters | Mean Std. Dev. ESS [95% Cred. Int]
-----------------------------+------------------------------------------------
Level 2: school |
var(cons) | .0959974 .0198513 2937 .0639578 .1399761
cov(cons,standlrt) | .0194119 .0072281 1811 .0068086 .0347633
var(standlrt) | .0154721 .0044998 1170 .008096 .025609
-----------------------------+------------------------------------------------
Level 1: student |
var(cons) | .554275 .0125677 5250 .5298432 .5796061
------------------------------------------------------------------------------
. ereturn list
scalars:
e(numlevels) = 2
e(k_f) = 2
e(k_r) = 4
e(k) = 6
e(size) = 10000
e(maxlevels) = 3
e(columns) = 1500
e(variables) = 5
e(tempmat) = 0
e(extrabinomial) = 0
e(iterations) = .
e(converged) = 0
e(burnin) = 500
e(chain) = 5000
e(thinning) = 1
e(dic) = 9214.65234375
e(pd) = 91.66587829589844
e(dbar) = 9122.986328125
e(dthetabar) = 9031.3212890625
e(mcmcdiagnostics) = 1
e(time) = 9.781000000000001
e(N) = 4059
e(level) = 95
macros:
e(chains) : "runmlwin_chains"
e(weighttype) : "standardised standardised"
e(weightvar) : ". ."
e(ivars) : "school"
e(level1var) : "student"
e(title) : "Normal response model"
e(method) : "MCMC"
e(properties) : "b V"
e(link) : "identity"
e(distribution) : "normal"
e(depvars) : "normexam"
e(version) : "2.25"
e(cmdline) : "runmlwin normexam cons standlrt, level2 (school: cons standlrt) level1 (student: cons) mcmc(on) initsprevious nopause"
e(cmd) : "runmlwin"
matrices:
e(b) : 1 x 6
e(V) : 6 x 6
e(g_max) : 1 x 1
e(g_avg) : 1 x 1
e(g_min) : 1 x 1
e(N_g) : 1 x 1
e(KD2) : 1000 x 6
e(KD1) : 1000 x 6
e(MCSE) : 1000 x 7
e(quantiles) : 9 x 7
e(PACF) : 10 x 7
e(ACF) : 100 x 7
e(pvalmode) : 1 x 6
e(pvalmedian) : 1 x 6
e(pvalmean) : 1 x 6
e(bd) : 1 x 6
e(rl2) : 1 x 6
e(rl1) : 1 x 6
e(ub) : 1 x 6
e(lb) : 1 x 6
e(ess) : 1 x 6
e(mode) : 1 x 6
e(sd) : 1 x 6
e(meanmcse) : 1 x 6
e(RP1) : 1 x 1
e(RP2) : 2 x 2
functions:
e(sample)
. estout, cells("b se var") stats(dic)
---------------------------------------------------
.
b se var
---------------------------------------------------
FP1
cons -.0056898 .0392222 .0015384
standlrt .5582691 .0204706 .000419
---------------------------------------------------
RP2
var(cons) .0959974 .0198513 .0003941
cov(cons\s~) .0194119 .0072281 .0000522
var(standl~) .0154721 .0044998 .0000202
---------------------------------------------------
RP1
var(cons) .554275 .0125677 .0001579
---------------------------------------------------
dic 9214.652
---------------------------------------------------
George