Graph with mean coefficients plut 95% credible intervals under - bayes - command

[Marcos Almeida]

Dear Forum Members,

I'm doing some Bayesian analysis under a hierarchical model. I wish to get a graph such as the one we get with - marginsplot - command under a frequentist frame.

I thought I could get it with the user-written (Ben Jann's) coefplot program since, according to the help files, it allows for using - matrix - command.

But I failed to get the appropriate matrix, since - bayes - provides 2 matrices: one for the mean coefficients and another for the credible intervals.

I believe I can use brutal force, I mean, I can copy and paste the information from both matrices in a single dataset, but I suspect there is a neat strategy, probably some "massage" concerning the use of - matrix - command...

Below, a toy example:

. bayes, normalprior(10) mcmcsize(1000) burnin(250): mixed mpg i.rep78 c.gear_ratio || foreign:
note: Gibbs sampling is used for regression coefficients and variance components
note: option adaptation(maxiter()) changed to 12
  
Burn-in 250 aa done
Simulation 1000 aaaaaaaaa1000 done

Multilevel structure
------------------------------------------------------------------------------
foreign
    {U0}: random intercepts
------------------------------------------------------------------------------

Model summary
------------------------------------------------------------------------------
Likelihood:
  mpg ~ normal(xb_mpg,{e.mpg:sigma2})

Priors:
  {mpg:i.rep78 gear_ratio _cons} ~ normal(0,100)                           (1)
                            {U0} ~ normal(0,{U0:sigma2})                   (1)
                  {e.mpg:sigma2} ~ igamma(.01,.01)

Hyperprior:
  {U0:sigma2} ~ igamma(.01,.01)
------------------------------------------------------------------------------
(1) Parameters are elements of the linear form xb_mpg.

Bayesian multilevel regression                   MCMC iterations  =      1,250
Metropolis-Hastings and Gibbs sampling           Burn-in          =        250
                                                 MCMC sample size =      1,000
Group variable: foreign                          Number of groups =          2

                                                 Obs per group:
                                                              min =         21
                                                              avg =       34.5
                                                              max =         48

                                                 Number of obs    =         69
                                                 Acceptance rate  =      .8539
                                                 Efficiency:  min =     .04763
                                                              avg =      .6294
Log marginal likelihood                                       max =          1
 
------------------------------------------------------------------------------
             |                                                Equal-tailed
             |      Mean   Std. Dev.     MCSE     Median  [95% Cred. Interval]
-------------+----------------------------------------------------------------
mpg          |
       rep78 |
          2  | -.6244307   2.847146   .100523  -.6267578  -6.383477   4.983676
          3  | -1.643648   2.636784   .083382  -1.592191   -6.71038   3.623224
          4  | -1.327764   2.806635   .097236  -1.374682  -6.673647    4.19083
          5  |  3.793827   3.150102     .1669   3.704065  -2.347833   9.942928
             |
  gear_ratio |  8.068264   1.716394   .148249   8.129553   4.532037   11.31953
       _cons | -2.993107   6.743514   .977144  -3.131476  -17.29503   10.58894
-------------+----------------------------------------------------------------
foreign      |
   U0:sigma2 |  76.87748   705.8558   22.3211    5.94588     .01881   311.9778
-------------+----------------------------------------------------------------
e.mpg        |
      sigma2 |  19.42023   4.544435   .154813   18.77323   13.22344   31.67032
------------------------------------------------------------------------------
Note: Default priors are used for model parameters.
Note: Adaptation continues during simulation.

. */ failing to get coefficients under - coefplot - command

. coefplot
(.: no coefficients found, all dropped, or none kept)
(nothing to plot)

. */ but they are there!

. matrix list e(mean)

e(mean)[1,9]
             mpg:        mpg:        mpg:        mpg:        mpg:        mpg:        mpg:         U0:
              1b.          2.          3.          4.          5.                                    
           rep78       rep78       rep78       rep78       rep78  gear_ratio       _cons      sigma2
Mean           0  -.62443075  -1.6436483  -1.3277645   3.7938272   8.0682642  -2.9931072   76.877476

           e.mpg:
                
          sigma2
Mean   19.420228

. matrix list e(cri)

e(cri)[2,9]
              mpg:        mpg:        mpg:        mpg:        mpg:        mpg:        mpg:         U0:
               1b.          2.          3.          4.          5.                                    
            rep78       rep78       rep78       rep78       rep78  gear_ratio       _cons      sigma2
Lower           0  -6.3834766  -6.7103804  -6.6736471  -2.3478329   4.5320369  -17.295032   .01880998
Upper           0   4.9836761    3.623224   4.1908303   9.9429277   11.319533   10.588937   311.97778

            e.mpg:
                  
           sigma2
Lower   13.223436
Upper   31.670324

. */ using - mixed - under the frequentist frame

. mixed mpg i.rep78 c.gear_ratio || foreign:

Performing EM optimization:

Performing gradient-based optimization:

Iteration 0:   log likelihood = -194.84877  
Iteration 1:   log likelihood = -194.83288  
Iteration 2:   log likelihood = -194.76731  
Iteration 3:   log likelihood = -194.75944  
Iteration 4:   log likelihood = -194.75939  
Iteration 5:   log likelihood = -194.75939  

Computing standard errors:

Mixed-effects ML regression                     Number of obs     =         69
Group variable: foreign                         Number of groups  =          2

                                                Obs per group:
                                                              min =         21
                                                              avg =       34.5
                                                              max =         48

                                                Wald chi2(5)      =      72.27
Log likelihood = -194.75939                     Prob > chi2       =     0.0000

------------------------------------------------------------------------------
         mpg |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
       rep78 |
          2  |  -.2455122   3.228211    -0.08   0.939    -6.572689    6.081665
          3  |  -1.313398   2.972638    -0.44   0.659    -7.139662    4.512866
          4  |  -1.352865   3.050797    -0.44   0.657    -7.332316    4.626587
          5  |   3.387386   3.164686     1.07   0.284    -2.815284    9.590057
             |
  gear_ratio |   7.449089   1.191231     6.25   0.000     5.114319     9.78386
       _cons |  -.6396043   4.500889    -0.14   0.887    -9.461184    8.181976
------------------------------------------------------------------------------

------------------------------------------------------------------------------
  Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
foreign: Identity            |
                  var(_cons) |   3.36e-21   1.01e-19      1.13e-46    99899.78
-----------------------------+------------------------------------------------
               var(Residual) |   16.56551    2.82032      11.86549    23.12725
------------------------------------------------------------------------------
LR test vs. linear model: chibar2(01) = 0.00          Prob >= chibar2 = 1.0000

. */ This is what I wish, but with the Credible Intervals, as shown after - bayes - command, above

. coefplot, drop (_cons _U0 _e.mpg)

This is the graph with - mixed -, but I would like to get a similar graph with - bayes: mixed -, which would provide the mean coefficients plus 95% credible intervals

Thanks in advance for any help!

[Andrew Musau]

You can try something like the following Marcos.

sysuse auto, clear
bayes, normalprior(10) mcmcsize(1000) burnin(250): mixed mpg i.rep78 c.gear_ratio || foreign:
mat R= e(mean)\e(cri)
mat list R
coefplot matrix(R), drop(_cons) ci((R[2] R[3]))  scheme(s1color)

Result:

[Marcos Almeida]

Andrew Musau : Lots of thanks, Andrew! This works like a charm!