I have noticed a minor issue with the mode and median options - both reporting options for the MCMC mode. When I use either option the column header for the random effects is appropriate - "mode" or "median", respectively. But in both cases the column header for the fixed effects remains "mean" although the estimates listed are not, in fact, the means. (I checked the e(mode) and e(quantiles) matrices to confirm that the fixed effects listed were modes and medians, respectively.) Therefore, I think that there is a simple labeling error - the column header for the fixed effects estimates should read "mode" or "median" instead of "mean".
Not sure that this matters too much, however. My understanding is that, in most cases, the mode, median, and mean will all be nearly identical for the fixed effects. (Not so for the random effects.)
But it did make we wonder about the underlying theory. Browne and Draper suggest that the following combinations of diffuse priors and posterior point summaries for variance components: medians for gamma priors and modes for uniform priors. (Browne WJ, Draper D. A comparison of Bayesian and likelihood-based methods for fitting multilevel models. Bayesian Analysis 2006;1(3):473-514.)
Never mind the fact that the different types of posterior point summaries are likely to differ only slightly for the fixed effects, should the same following pairings be applied to them when reporting results?
Thanks. BPL
reporting fixed effects modes and medians
-
GeorgeLeckie
- Site Admin
- Posts: 432
- Joined: Fri Apr 01, 2011 2:14 pm
Re: reporting fixed effects modes and medians
Hi BPL,
Thank you for your intersting post.
We are struggling to replicate the means, models, median labelling display error that you mention. Please could you check that you are using the latest version of runmlwin
If you still get the labelling display error after this it would be great if you could post your syntax and output to see if that makes it more obvious where we have gone wrong
Best wishes
George
Thank you for your intersting post.
We are struggling to replicate the means, models, median labelling display error that you mention. Please could you check that you are using the latest version of runmlwin
Code: Select all
ssc install runmlwin, replaceBest wishes
George
-
bpeterlucas
- Posts: 8
- Joined: Sun Jan 29, 2012 5:55 pm
Re: reporting fixed effects modes and medians
George - I confirmed that I have the most recent version of runmlwin.
Below I have pasted my syntax for the ordinal model (outcome takes values 0, 1, 2) with a non-hierarchical structure (aid and resident are crossed) using uniform random effects priors.
I follow this syntax with 2 commands: "runmlwin, median" and "runmlwin, mode". You'll see that the fixed effects change for the latter 2 runmlwin commands, but the column header for the fixed effects does not change - it remains "mean".
This is the source of my question, which I suspect exposes my lack of understanding of the MCMC method. Thank you for bearing with me. If I wish to use the mode posteriors (following the pairings suggested in the Browne and Draper article that I mentioned earlier), should I report the mode fixed effects - reported as "mean" after the mode option is used? Or should i use the mean fixed effects - reported as "mean" after the mean (default) option is used?
As I mentioned, it doesn't make much substantive difference. For example, for the fixed effect "h4_01" the result after the default "mean" option is -0.9942244 whereas it is -0.997683 after the "mode" option.
Thanks. BPL
Below I have pasted my syntax for the ordinal model (outcome takes values 0, 1, 2) with a non-hierarchical structure (aid and resident are crossed) using uniform random effects priors.
I follow this syntax with 2 commands: "runmlwin, median" and "runmlwin, mode". You'll see that the fixed effects change for the latter 2 runmlwin commands, but the column header for the fixed effects does not change - it remains "mean".
This is the source of my question, which I suspect exposes my lack of understanding of the MCMC method. Thank you for bearing with me. If I wish to use the mode posteriors (following the pairings suggested in the Browne and Draper article that I mentioned earlier), should I report the mode fixed effects - reported as "mean" after the mode option is used? Or should i use the mean fixed effects - reported as "mean" after the mean (default) option is used?
As I mentioned, it doesn't make much substantive difference. For example, for the fixed effect "h4_01" the result after the default "mean" option is -0.9942244 whereas it is -0.997683 after the "mode" option.
Thanks. BPL
Code: Select all
runmlwin `var' cons (h4 hab clr2 clr3 clr4 rotc rotcsqrd, contrast(1/2)), ///
level3(aid: (cons, contrast(1/2))) ///
level2(resident: (cons, contrast(1/2))) ///
level1(eval_id) discrete(distribution(multinomial) link(ologit) denominator(cons) basecategory(2)) ///
mcmc(cc b($burnlen) c($chainlen) t($thinamt) rppriors(uniform) log) /// ///
initsb(ints) /// ///
nopause
MLwiN 2.24 multilevel model Number of obs = 688
Ordered multinomial logit response model
Estimation algorithm: MCMC
-----------------------------------------------------------
| No. of Observations per Group
Level Variable | Groups Minimum Average Maximum
----------------+------------------------------------------
aid | 60 4 11.5 30
resident | 147 1 4.7 11
-----------------------------------------------------------
----------------------------------
Contrast | Log-odds
-------------+--------------------
1 | 0 vs. 1 2
2 | 0 1 vs. 2
----------------------------------
Burnin = 500
Chain = 1000
Thinning = 1
Run time (seconds) = 17.4
Deviance (dbar) = 960.15
Deviance (thetabar) = 793.71
Effective no. of pars (pd) = 166.43
Bayesian DIC = 1126.58
------------------------------------------------------------------------------
| Mean Std. Dev. z ESS [95% Cred. Interval]
-------------+----------------------------------------------------------------
Contrast 1 |
cons_0 | .9681302 .2656094 3.64 17 .5296707 1.532364
-------------+----------------------------------------------------------------
Contrast 2 |
cons_1 | 4.296173 .4677286 9.19 13 3.497586 5.223987
-------------+----------------------------------------------------------------
h4_01 | -.9942244 .2950549 -3.37 39 -1.640761 -.4731149
hab_01 | .124938 .2858505 0.44 138 -.4156916 .692005
clr2_01 | -.1414844 .3084645 -0.46 51 -.7451616 .4699075
clr3_01 | -.4078354 .349414 -1.17 31 -1.104277 .2575084
clr4_01 | -.5087692 .3200843 -1.59 42 -1.148904 .0646879
rotc_01 | -.00222 .0298807 -0.07 121 -.0588807 .065421
rotcsqrd_01 | .0033604 .0080624 0.42 76 -.0119433 .0196369
------------------------------------------------------------------------------
------------------------------------------------------------------------------
Random-effects Parameters | Mean Std. Dev. ESS [95% Cred. Int]
-----------------------------+------------------------------------------------
Level 3: aid |
var(cons_01) | 1.090046 .5682853 26 .4442516 2.708316
-----------------------------+------------------------------------------------
Level 2: resident |
var(cons_01) | 1.834842 .9518038 10 .6043799 4.189524
------------------------------------------------------------------------------
. runmlwin, median
MLwiN 2.24 multilevel model Number of obs = 688
Ordered multinomial logit response model
Estimation algorithm: MCMC
-----------------------------------------------------------
| No. of Observations per Group
Level Variable | Groups Minimum Average Maximum
----------------+------------------------------------------
aid | 60 4 11.5 30
resident | 147 1 4.7 11
-----------------------------------------------------------
----------------------------------
Contrast | Log-odds
-------------+--------------------
1 | 0 vs. 1 2
2 | 0 1 vs. 2
----------------------------------
Burnin = 500
Chain = 1000
Thinning = 1
Run time (seconds) = 17.4
Deviance (dbar) = 960.15
Deviance (thetabar) = 793.71
Effective no. of pars (pd) = 166.43
Bayesian DIC = 1126.58
------------------------------------------------------------------------------
| Mean Std. Dev. z ESS [95% Cred. Interval]
-------------+----------------------------------------------------------------
Contrast 1 |
cons_0 | .9291338 .2656094 3.64 17 .5296707 1.532364
-------------+----------------------------------------------------------------
Contrast 2 |
cons_1 | 4.252225 .4677286 9.19 13 3.497586 5.223987
-------------+----------------------------------------------------------------
h4_01 | -.9879735 .2950549 -3.37 39 -1.640761 -.4731149
hab_01 | .1252627 .2858505 0.44 138 -.4156916 .692005
clr2_01 | -.1548411 .3084645 -0.46 51 -.7451616 .4699075
clr3_01 | -.4055779 .349414 -1.17 31 -1.104277 .2575084
clr4_01 | -.4852518 .3200843 -1.59 42 -1.148904 .0646879
rotc_01 | -.0017183 .0298807 -0.07 121 -.0588807 .065421
rotcsqrd_01 | .00338 .0080624 0.42 76 -.0119433 .0196369
------------------------------------------------------------------------------
------------------------------------------------------------------------------
Random-effects Parameters | Median Std. Dev. ESS [95% Cred. Int]
-----------------------------+------------------------------------------------
Level 3: aid |
var(cons_01) | .9719989 .5682853 26 .4442516 2.708316
-----------------------------+------------------------------------------------
Level 2: resident |
var(cons_01) | 1.651547 .9518038 10 .6043799 4.189524
------------------------------------------------------------------------------
. runmlwin, mode
MLwiN 2.24 multilevel model Number of obs = 688
Ordered multinomial logit response model
Estimation algorithm: MCMC
-----------------------------------------------------------
| No. of Observations per Group
Level Variable | Groups Minimum Average Maximum
----------------+------------------------------------------
aid | 60 4 11.5 30
resident | 147 1 4.7 11
-----------------------------------------------------------
----------------------------------
Contrast | Log-odds
-------------+--------------------
1 | 0 vs. 1 2
2 | 0 1 vs. 2
----------------------------------
Burnin = 500
Chain = 1000
Thinning = 1
Run time (seconds) = 17.4
Deviance (dbar) = 960.15
Deviance (thetabar) = 793.71
Effective no. of pars (pd) = 166.43
Bayesian DIC = 1126.58
------------------------------------------------------------------------------
| Mean Std. Dev. z ESS [95% Cred. Interval]
-------------+----------------------------------------------------------------
Contrast 1 |
cons_0 | .8442673 .2656094 3.64 17 .5296707 1.532364
-------------+----------------------------------------------------------------
Contrast 2 |
cons_1 | 4.158042 .4677286 9.19 13 3.497586 5.223987
-------------+----------------------------------------------------------------
h4_01 | -.997683 .2950549 -3.37 39 -1.640761 -.4731149
hab_01 | .1558631 .2858505 0.44 138 -.4156916 .692005
clr2_01 | -.2170185 .3084645 -0.46 51 -.7451616 .4699075
clr3_01 | -.4072731 .349414 -1.17 31 -1.104277 .2575084
clr4_01 | -.4100945 .3200843 -1.59 42 -1.148904 .0646879
rotc_01 | -.0013399 .0298807 -0.07 121 -.0588807 .065421
rotcsqrd_01 | .0018674 .0080624 0.42 76 -.0119433 .0196369
------------------------------------------------------------------------------
------------------------------------------------------------------------------
Random-effects Parameters | Mode Std. Dev. ESS [95% Cred. Int]
-----------------------------+------------------------------------------------
Level 3: aid |
var(cons_01) | .8741872 .5682853 26 .4442516 2.708316
-----------------------------+------------------------------------------------
Level 2: resident |
var(cons_01) | 1.413705 .9518038 10 .6043799 4.189524
-------------------------------------------------------------------------------
GeorgeLeckie
- Site Admin
- Posts: 432
- Joined: Fri Apr 01, 2011 2:14 pm
Re: reporting fixed effects modes and medians
Hi BPL,
Thank you for pasting the output. We have now found the labelling bug and fixed it for the next release.
In terms of your query, I am not sure what would be recommended, though as you point out it will make little difference for the fixed part parameters.
You could try emailing the authors of the article and if you get a reply, by all means post it here.
George
Thank you for pasting the output. We have now found the labelling bug and fixed it for the next release.
In terms of your query, I am not sure what would be recommended, though as you point out it will make little difference for the fixed part parameters.
You could try emailing the authors of the article and if you get a reply, by all means post it here.
George