MCMC posterior distribution
MCMC posterior distribution
I used the MCMC estimation method to get an estimate for a parameter of a 2 level logistic regression analysis. Based on the regression diagnostics of a parameter I would like to obtain a 99% prediction interval for this parameter instead of a 95% prediction interval. As I understand it I need the 1% and 99% quantile values of the posterior distribution. Default the program shows me the 2.5% and 97.5% quatile values. Is it possible to change the quatiles of the posterior distribution presented in the MCMC diagnostics from 2.5% and 97.5% into 1% and 99%. I cannot find in the manual how to do this. The diagnostic setting button also offers no solution here.
-
ChrisCharlton
- Posts: 1384
- Joined: Mon Oct 19, 2009 10:34 am
Re: MCMC posterior distribution
Unfortunately these values are not changeable within MLwiN, however you can calculate them yourselves from the chains. An example of the commands for this is in the quantiles.txt macro provided with MLwiN:
The chains of the parameters are stored stacked in c1090, so the first command generates a column indicating which rows correspond to each parameter. To adapt this to your model you will need to change the 3 to the total number of parameters in your model and the 10001 to the number of iterations that you ran the model for. You will also need to change c20 to an empty column in your worksheet to store this indicator.
The "split" line creates a column for each parameter chain. You will need to change c20 to the column used above and change the c21-c23 to refer to a list of empty columns, one for each parameter.
As in the above code you will then need to sort each of these parameter chain columns generated in the previous step.
To get your required quantile value you will then need to calculate where this lies given the number of iterations that the model has run for. In the above example the chain is 10001 long, so the 2.5% quantile will be in row 251 and the 97.5% quantile will be in row 9751. You will need to extract these for each parameter that you are interested in. Note that these values may be slightly different to what MLwiN reports in the MCMC diagnostics window for some chain lengths, as this method is doing no interpolation if a row doesn't correspond exactly to a quantile.
The final line above simply joins the extracted quantiles together into a single column.
Alternatively if you were to call MLwiN from Stata via runmlwin (http://www.bris.ac.uk/cmm/software/runmlwin/) you should be able to retrieve the required quantiles directly through the level() option. Unfortunately this functionality is currently broken, but will be fixed in the next release.
Code: Select all
note macro returns the 2.5% and 97.5% quantile estimates
code 3 1 10001 c20
split c1090 c20 c21-c23
sort c21 c21
sort c22 c22
sort c23 c23
pick 251 c21 b21
pick 9751 c21 b22
pick 251 c22 b23
pick 9751 c22 b24
pick 251 c23 b25
pick 9751 c23 b26
join b21-b26 c15 c15
The "split" line creates a column for each parameter chain. You will need to change c20 to the column used above and change the c21-c23 to refer to a list of empty columns, one for each parameter.
As in the above code you will then need to sort each of these parameter chain columns generated in the previous step.
To get your required quantile value you will then need to calculate where this lies given the number of iterations that the model has run for. In the above example the chain is 10001 long, so the 2.5% quantile will be in row 251 and the 97.5% quantile will be in row 9751. You will need to extract these for each parameter that you are interested in. Note that these values may be slightly different to what MLwiN reports in the MCMC diagnostics window for some chain lengths, as this method is doing no interpolation if a row doesn't correspond exactly to a quantile.
The final line above simply joins the extracted quantiles together into a single column.
Alternatively if you were to call MLwiN from Stata via runmlwin (http://www.bris.ac.uk/cmm/software/runmlwin/) you should be able to retrieve the required quantiles directly through the level() option. Unfortunately this functionality is currently broken, but will be fixed in the next release.
Re: MCMC posterior distribution
thank you very much for your clear reply
Maarten
Maarten