Posterior predictive distribution after bayesmh

[Josh Doyle]

Hi, I am trying to calculate the posterior predictive distribution after fitting a model using bayesmh in Stata. Has anyone done this in Stata, yet? I have not been able to find any step-by-step examples of code online, and it seems like this should be fairly simple given Stata's margins commands. Thanks!

[Joseph Coveney]

No one's replied so far, and I'm not aware of any step-by-step examples of Stata code online, but here's my try below.

It looks long-winded, but basically it just saves the draws from the posterior distribution created by bayesmh, randomly selects (with replacement) elements of Theta from the posterior distribution, generates a predicted observation from the selected Theta using the data model (likelihood, which is conveniently shown in the header of the output of bayesmh), and repeats until you're satisfied that you've sampled the posterior predictive distribution enough.

For illustration, I used the first example from the help file for bayesmh, and chose 5000 samples as adequate to flesh out the major features of the posterior predictive distribution. For examples of graphs that you can do, I plotted the posterior predictive distribution for each observation against the corresponding observed value of the response variable in a histogram (in the manner of a visual posterior predictive check), and then taken the expectation of the predictive distribution for each observation and plotted it against the observed value of the response variable in a scatterplot. Other types of plots are possible.

version 14.2

clear *
set more off
set seed 1377321

*
* Save data for computing linear predictions from randomly drawn theta
*
webuse oxygen

tempfile data
quietly save `data'

*
* Fit model and save draws from the posterior distribution
*
tempfile posterior
bayesmh change age group, ///
    likelihood(normal({var})) ///
    prior({change:}, flat) prior({var}, jeffreys) ///
    saving(`posterior')

/* Randomly draw Theta from the posterior distribution and
   generate posterior predictive distribution from the data model */

*
* Set up posterior distribution so that it can be sampled randomly
*
use `posterior'
quietly expand _frequency
keep eq*
generate double randu = runiform()
isid randu
sort randu

// Create index on posterior distribution to allow sampling
generate long row = _n
summarize row, meanonly
local Rows = r(max)

*
* Take random samples of Theta (independently elementwise) from posterior distribution
*
// Set up file to store random samples of Theta
tempname file_handle
tempfile Theta
postfile `file_handle' double (eq1_p1 eq1_p2 eq1_p3 eq0_p1) using `Theta'

// To allow independent draws of each element of Theta
tempname eq1_p1 eq1_p2 eq1_p3 eq0_p1

forvalues rep = 1/5000 {
    foreach theta in eq1_p1 eq1_p2 eq1_p3 eq0_p1 {
        local selected = runiformint(1, `Rows')
        scalar define ``theta'' = `theta'[`selected']
    }
    post `file_handle' (`eq1_p1') (`eq1_p2') (`eq1_p3') (`eq0_p1')
}
postclose `file_handle'

*
* Generate the posterior predictive distribution (y_predicted) from Theta samples and predictors in data
*
quietly use `Theta', clear
cross using `data'
// Compute linear prediction and generate y_predicted ("change_hat") from fitted model
generate double change_hat = ///
    rnormal(eq1_p1 * age + eq1_p2 * group + eq1_p3, sqrt(eq0_p1))

// keep y_observed and y_predicted
keep change*

*
* Compare probability distribution of y_predicted against each y_observed
*
quietly levelsof change, local(changes)
pause on
foreach change in `changes' {
    quietly histogram change_hat if float(change) == float(`change'), ///
        color(none) lcolor(gs8) ///
        xline(`change', lcolor(black) lpattern(dash)) ///
        xtitle(Posterior Predictive Distribution versus Observed Response) ///
        ylabel( , angle(horizontal) nogrid) legend(off)
    pause
}

*
* Compare expectation of y_predicted to y_observed
*
set type double
collapse (mean) change_hat = change_hat, by(change)

graph twoway ///
    scatter change_hat change, mcolor(black) msize(small) || ///
    line change change, lcolor(black) lpattern(dash) ///
        ytitle("Posterior Prediction (Expectation)") ///
        ylabel( , angle(horizontal) nogrid) legend(off)

exit

[Josh Doyle]

Thanks, Joseph. This worked well.