Saving parts of a model after Bayesian analysis

[Pantelis Kazakis]

Dear all,

Assume that one runs a model using the bayes prefix. This model includes also a large number of indicators apart from some standard controls.

I am wondering if there is a way not to have the results of the indicators shown at all in the Stata output.

I will will present a minimal working example.

input id    year    y    x1    x2
1    1990    0.765290722    0.786241262    0.783572146
1    1991    0.611380163    0.739519495    0.60693537
1    1992    0.590601317    0.490306934    0.479196734
1    1993    0.708316766    0.97180118    0.40444494
1    1994    0.587700824    0.442869206    0.399139847
1    1995    0.092595406    0.475633824    0.055023827
2    1990    0.995706082    0.356060781    0.242968298
2    1991    0.688079045    0.36475865    0.829626742
2    1992    0.564536023    0.219309875    0.59165049
2    1993    0.728154736    0.370715204    0.319482427
2    1994    0.972550408    0.331836688    0.495801468
2    1995    0.907903102    0.643605887    0.04147437
3    1990    0.666757644    0.604741524    0.892619352
3    1991    0.945322031    0.076470116    0.83826142
3    1992    0.262225766    0.622579231    0.738801605
3    1993    0.14160545    0.420033382    0.883489729
3    1994    0.674649086    0.287176323    0.798404419
3    1995    0.462505639    0.782454944    0.206593929
4    1990    0.103318329    0.860360564    0.384564281
4    1991    0.892686701    0.387916868    0.911908676
4    1992    0.23782402    0.403880512    0.962117085
4    1993    0.733468764    0.376152156    0.237543589
4    1994    0.925702313    0.438672271    0.103990407
4    1995    0.309275933    0.771466795    0.119082256
end

Suppose we run the following model:

bayes, rseed(12345) : reg y x1 x2 i.id i.year

The results we get are:

Bayesian linear regression                       MCMC iterations  =     12,500
Random-walk Metropolis-Hastings sampling         Burn-in          =      2,500
                                                 MCMC sample size =     10,000
                                                 Number of obs    =         24
                                                 Acceptance rate  =       .327
                                                 Efficiency:  min =    .003156
                                                              avg =    .006709
Log marginal likelihood = -74.144851                          max =     .02003
 
------------------------------------------------------------------------------
             |                                                Equal-tailed
             |      Mean   Std. Dev.     MCSE     Median  [95% Cred. Interval]
-------------+----------------------------------------------------------------
y            |
          x1 |  .0390891   .1586424   .025329   .0357585  -.2741725   .3541185
          x2 | -.3099317   .2938798   .039644  -.3094318  -.8912901   .2766925
             |
          id |
          2  |  .2899881   .1213703   .021604   .2966582   .0641597   .5389639
          3  |   .055742   .1603441   .022674   .0559719  -.2698375   .3579688
          4  | -.0549487   .1484017   .016786  -.0567729  -.3385171   .2308046
             |
        year |
       1991  |  .2982491   .1879188   .026803   .3004807  -.1168325   .6647947
       1992  | -.1464675   .1949344   .020653  -.1421357  -.5566844    .232561
       1993  | -.0215075   .1920174   .031623  -.0132725  -.4371793   .3662122
       1994  |  .1843418   .1978521   .020934    .188863  -.2352899   .5520474
       1995  | -.2889066   .2437549    .03487  -.2840871  -.7785481   .1575815
             |
       _cons |  .6705034   .2512025   .040967   .6678414   .1940561   1.183284
-------------+----------------------------------------------------------------
      sigma2 |  .0759295   .0308258   .002178   .0698249   .0363304   .1570932
------------------------------------------------------------------------------

Is it possible to show (or save) only the results for x1 and x2? It appears that estimates table does not work with Bayes.

Thank you.

[Phil Bromiley]

This is probably more trouble than it is worth, but you could run the regression using quietly, and then manipulate the parameter matrices. But this seems like a lot of trouble to avoid having excessl lines in the output window. If you want to save specific variables, then you can access the parameters and save them. for example local bx1=_b[x1] will normally put the parameter on x1 into the local bx1. Look at the User Manual on accessing results.

[Weiwen Ng]

Actually, I think there's a way to do this. Normally, the post-estimation command bayesstats summary will display the entire summary table. You can use that command to display only selected variables:

use http://www.stata-press.com/data/r15/womenwage
bayes, saving(bayesreg1, replace): regress wage age
bayesstats summary {wage:age} {wage:_cons}
Posterior summary statistics                      MCMC sample size =    10,000
 
------------------------------------------------------------------------------
             |                                                Equal-tailed
        wage |      Mean   Std. Dev.     MCSE     Median  [95% Cred. Interval]
-------------+----------------------------------------------------------------
         age |  .3983066   .0610067   .001723   .4002309   .2797012   .5185586
       _cons |  6.063636   1.801913   .051992    6.05618   2.496192   9.587002
------------------------------------------------------------------------------

This will be cumbersome if you have many coefficients you're trying to save. Also, I don't see a way to automatically export these results; you'll have to copy and paste them into Excel, and I know this process is error-prone.

As you noted, the Bayes prefix appears isn't compatible with estimates table. It doesn't post the matrices e(b) and e(V) - these are the means of the betas, and the variance-covariance matrix of the betas. Most estimation commands take those matrices, calculate the standard errors, test statistics, and 95% CIs, and post them to a matrix called r(table). It appears you need e(b) and e(V) to be compatible with estimates store, which would then enable you to easily export the model output using user-written commands like estout (available on SSC). It also doesn't appear there's a matrix with the mean, SD, MCSE, and the lower and upper credible interval values that you could use putexcel on, or pass to estout.

One partial workaround is that you could save the results of your Monte Carlo simulation, as I did above, to a .dta file. You can then access that file, and calcualte the mean and SD from there.

use bayesreg1
sum eq1_p1 eq1_p2 [fw=_frequency]
    Variable |        Obs        Mean    Std. Dev.       Min        Max
-------------+---------------------------------------------------------
      eq1_p1 |     10,000    .3983066    .0610067   .1818681    .622838
      eq1_p2 |     10,000    6.063636    1.801913  -1.025328   13.28191

You'll see the mean and SD of the parameter estimates correspond to the table above. I don't know how to compute the Monte Carlo SE, but I bet this is doable. The credible interval tails are more problematic, as they're the 2.5th and 97.5th quantiles of the two variables. You could do

 _pctile eq1_p1 [fw=_frequency], percentiles(2.5 97.5)
return list

scalars:
                 r(r1) =  .2797012453776874
                 r(r2) =  .5185585837551255

Again, that's the 95% credible interval for age.

However, this is much more trouble than it's probably worth. It would be better if the bayes prefix commands posted to r(table). This would enable compatibility with a bunch of commands that export regression results - if nothing else, you can save r(table) as a matrix and then use putexcel to write it. You're not the first person to raise this issue, and as Bayesian estimation gets more popular, you probably won't be the last, unless all the Bayesians jump ship to R.

[Pantelis Kazakis]

Dear all,

Thank you for the information you have provided me so far. I hope in the future things will become more automated.

[Ben A. Dwamena]

Pantelis,
See documentation for bayesstats summary. Run your model with the modelsummary option to identify how bayes labels the parameters you want to save. Then after estimation, use bayesstats summary to summarize only the parameters you want. It will be returned/saved as r(summary). HTH,
Ben