// #######################################################################

// bayes basic example - drug study
set seed 7431
use "impact.dta", clear
replace time=610 if time==.
stset time, failure(relapse)
streg i.treat i.site, distribution(weibull) nohr
bayes: streg i.treat i.site, distribution(weibull) nohr

bayesgraph diagnostics _all
// mixing looks ok but autocorrelated - we would like more posterior draws with the mcmcsize option
// we'll also get dots to show progress with the dots option
bayes, mcmcsize(50000) dots(1000): streg i.treat i.site, distribution(weibull)
/* ok, just a little slow. There are options for thinning but we will side with 
   Christian Robert (Chapter 1, Handbook of MCMC) and keep them all.
*/
bayesstats ess // it's slow to mix but not terrible; we might even want more than 50,000 draws

// we can find the probability for the HR for long treatment < 0.85 (i.e. {_t:1.treat} < -0.1625)
bayestest interval {_t:1.treat}, upper(-0.1625)

// let's save the draws to a new file
bayes, mcmcsize(50000) dots(1000) saving(impact_bayes, replace): streg i.treat i.site, distribution(weibull)
use impact_bayes.dta, clear
summarize [fw=_frequency]
// to understand the varnames:
display e(parnames)
display e(varnames)

// or expand it
expand _frequency
sort _index
generate newindex = _n

// joint distribution of the two log-hazard ratios
scatter eq1_p2 eq1_p4, msize(vsmall) mlalign(outside) mlwidth(none) mcolor(%6) ///
	xtitle(site) ytitle(duration) xline(0) yline(0)

gen sitecat=floor(eq1_p4*20)/20
gen duracat=floor(eq1_p2*20)/20
collapse (count) n=newindex, by(sitecat duracat)
twoway contour n duracat sitecat, levels(10) scolor(blue) ecolor(yellow) ///
	xtitle(site) ytitle(duration) ztitle("Posterior draws")





// #######################################################################

// auto example - specifying priors and why they matter

set seed 7431
sysuse auto, clear
regress price mpg
bayes: regress price mpg

use `e(filename)', clear
qui summ eq1_p1
global badmpg=r(mean)
qui summ eq1_p2
global badcons=r(mean)
sysuse auto, clear

// change priors & initials
bayes, prior({price: mpg}, normal(0, 4000000)) ///
	   prior({price: _cons}, uniform(-10000,25000)) ///
	   prior({sigma2}, uniform(0.1, 5000)) ///
	   initial({price: _cons} 12000 {price: mpg} -300 {sigma2} 500)  ///
	   : regress price mpg
use `e(filename)', clear
qui summ eq1_p1
global goodmpg=r(mean)
qui summ eq1_p2
global goodcons=r(mean)
sysuse auto, clear


// the default priors are highly informative!
gen badpred = $badcons + $badmpg * mpg
gen goodpred = $goodcons + $goodmpg * mpg

sort mpg
twoway (scatter price mpg, ytitle("Price (USD)")) ///
	   (lfit price mpg, lcolor(navy)) ///
	   (line badpred mpg, lcolor(cranberry) lpattern(dash)) ///
	   (line goodpred mpg, lcolor(cranberry) lwidth(thick) lpattern(dash)) ///
	   , legend(order(1 "Data" 2 "OLS" 3 "Bayes - default priors" 4 " Bayes - WIP priors"))