- Dr Steven Bowe
- Jun 20, 2023
- 0 min read
top of page
Research Blog [ webpages currently under development ]
Search
Stata Version 18 has taken the plunge further into Meta-Analysis and Bayesian methods. So I decided to look at how we might be able to show results for a forest plot that includes a Bayesian Effect Size estimate.
The following example uses a dataset and meta and bayesmh codes from Stata 18. I have shown how to include a Bayesian mean estimate and 95% Credibility Interval as an additional diamond on the forest plot of a Random Effect model using Restricted Maximum Likeihood (REML).
Example 26: Normal–normal analysis
"Here we follow the approach of Carlin (1992) for the normal–normal analysis of the beta-blockers data. For our normal–normal analysis, we consider data in wide form and concentrate on modeling estimates of log odds-ratios from 22 studies." (Stata 18 )
use https://www.stata-press.com/data/r18/betablockers_wide, clear
(Beta-blockers data in wide form)
describe
/* Contains data from https://www.stata-press.com/data/r18/betablockers_wide.dta
Observations: 22 Beta-blockers data in wide form
Variables: 7 5 Feb 2022 19:02
(_dta has notes)
--------------------------------------------------------------------------------------------------
Variable Storage Display Value
name type format label Variable label
--------------------------------------------------------------------------------------------------
study byte %9.0g Study identifier
deaths0 int %9.0g Number of deaths in the control group
total0 int %9.0g Number of subjects in the control group
deaths1 int %9.0g Number of deaths in the treatment group
total1 int %9.0g Number of subjects in the treatment group
D double %10.0g Log odds-ratio (based on empirical logits)
var double %10.0g Squared standard error of log odds-ratio
--------------------------------------------------------------------------------------------------
Sorted by:
*/
local color3 mcolor("maroon%70")
local color2 mcolor("blue%80")
local color1 mcolor("green%80")
bayesmh D D[study], likelihood(normal(var)) noconstant ///
prior({D[study]}, normal({d},{sig2})) ///
prior({d}, normal(0,1000)) ///
prior({sig2}, igamma(0.001,0.001)) ///
block({sig2}, gibbs) ///
block({d}, gibbs) ///
rseed(17)
return list
bayesstats ess {d} {sig2}
graph export bayesdiagn.png, replace
bayesgraph diagnostics {d} {sig2}
matrix bayesES= e(mean)
matrix list bayesES
matrix EScri= e(cri)
matrix list EScri
global bayesES=round(bayesES[1,1], 0.0001)
global lower=round(EScri[1,1], 0.0001)
global upper=round(EScri[2,1], 0.0001)
di $bayesES
di $lower
di $upper
gen SE= sqrt(var)
meta set D SE
meta summarize , random(reml)
meta forestplot _id _plot _esci , random(reml) ciopts(color(navy)) omarkeropts(mcolor(maroon)) markeropts(mcolor(navy)) ///
nullrefline(lcolor(gs12) lwidth(thin)) esrefline(lcolor(maroon) lpattern(shortdash) lwidth(thin)) cibind(parentheses) xlabel(, labsize(small)) title("") columnopts(_esci , format(%9.3f)) ///
customoverall( $bayesES $lower $upper, label("Bayesian Mean [95% cred. interval] ") `color1' ) ///
note("[Random (Reml) ] " , size(small)) ///
name(graph1, replace)
Stata Outputs
bayesmh D D[study], likelihood(normal(var)) noconstant
prior({D[study]}, normal({d},{sig2}))
prior({d}, normal(0,1000))
prior({sig2}, igamma(0.001,0.001))
block({sig2}, gibbs)
block({d}, gibbs)
rseed(17)
Burn-in 2500 aaaaaaaaa1000aaaaaaaaa2000aaaaa done
Simulation 10000 .........1000.........2000.........3000.........4000.........5000.........6000.........7000.........8000.........9000.........10000 done
Model summary
------------------------------------------------------------------------------
Likelihood:
D ~ normal(xb_D,var)
Prior:
{D[study]} ~ normal({d},{sig2}) (1)
Hyperpriors:
{d} ~ normal(0,1000)
{sig2} ~ igamma(0.001,0.001)
------------------------------------------------------------------------------
(1) Parameters are elements of the linear form xb_D.
Bayesian normal regression MCMC iterations = 12,500
Metropolis–Hastings and Gibbs sampling Burn-in = 2,500
MCMC sample size = 10,000
Number of obs = 22
Acceptance rate = .7623
Efficiency: min = .02206
avg = .02348
Log marginal-likelihood max = .02491
---------------------------------------------------------------------------------------------------
| Equal-tailed
| Mean Std. dev. MCSE Median [95% cred. interval]
-----------+-------------------------------------------------------------------------------------
d | -.2537001 .0648291 .004107 -.2574083 -.371893 -.1213832
sig2 | .0191485 .0212749 .001433 .0115096 .0013426 .078143
----------------------------------------------------------------------------------------------------
bayesstats ess {d} {sig2}
Efficiency summaries MCMC sample size = 10,000
Efficiency: min = .02206
avg = .02348
max = .02491
----------------------------------------------
| ESS Corr. time Efficiency
---------+------------------------------------
d | 249.13 40.14 0.0249
sig2 | 220.55 45.34 0.0221
----------------------------------------------
ereturn list
matrix bayesES= e(mean)
matrix list bayesES
bayesES[1,2]
d sig2
Mean -.25370012 .01914852
matrix EScri= e(cri)
matrix list EScri
EScri[2,2]
d sig2
Lower -.37189299 .00134263
Upper -.12138319 .07814304
global bayesES=round(bayesES[1,1], 0.0001)
global lower=round(EScri[1,1], 0.0001)
global upper=round(EScri[2,1], 0.0001)
gen SE= sqrt(var)
di $bayesES
-.2537
di $lower
-.3719
di $upper
-.1214
meta set D SE
meta summarize , random(reml)
Effect-size label: Effect size
Effect size: D
Std. err.: SE
Meta-analysis summary Number of studies = 22
Random-effects model Heterogeneity:
Method: REML tau2 = 0.0142
I2 (%) = 19.50
H2 = 1.24
--------------------------------------------------------------------
Study | Effect size [95% conf. interval] % weight
------------------+-------------------------------------------------
Study 1 | 0.028 -1.638 1.695 0.50
Study 2 | -0.741 -1.688 0.206 1.48
Study 3 | -0.541 -1.647 0.566 1.10
Study 4 | -0.246 -0.517 0.025 11.01
Study 5 | 0.069 -0.481 0.620 3.94
Study 6 | -0.584 -1.909 0.740 0.78
Study 7 | -0.512 -0.784 -0.241 10.96
Study 8 | -0.079 -0.478 0.321 6.57
Study 9 | -0.424 -0.961 0.113 4.11
Study 10 | -0.335 -0.564 -0.105 13.14
Study 11 | -0.213 -0.595 0.169 7.02
Study 12 | -0.039 -0.489 0.411 5.48
Study 13 | -0.593 -1.427 0.240 1.88
Study 14 | 0.282 -0.121 0.684 6.50
Study 15 | -0.321 -0.905 0.262 3.56
Study 16 | -0.135 -0.647 0.376 4.45
Study 17 | 0.141 -0.573 0.854 2.50
Study 18 | 0.322 -0.761 1.405 1.15
Study 19 | 0.444 -0.960 1.849 0.69
Study 20 | -0.218 -0.727 0.292 4.48
Study 21 | -0.591 -1.095 -0.087 4.56
Study 22 | -0.608 -1.142 -0.074 4.15
------------------+-------------------------------------------------
theta | -0.247 -0.366 -0.129
--------------------------------------------------------------------
Test of theta = 0: z = -4.09 Prob > |z| = 0.0000
Test of homogeneity: Q = chi2(21) = 23.26 Prob > Q = 0.3302
meta forestplot _id _plot _esci, random(reml) ciopts(color(navy)) omarkeropts(mcolor(maroon)) /// markeropts(mcolor(navy)) nullrefline(lcolor(gs12) lwidth(thin)) esrefline(lcolor(maroon) /// lpattern(shortdash) lwidth(thin)) cibind(parentheses) xlabel(, labsize(small)) title("") columnopts(_esci , format(%9.3f)) ///
customoverall( $bayesES $lower $upper, label("Bayesian Mean [95% cred. interval] ") `color1' ) ///
note("[Random (Reml) ] " , size(small)) ///
name(graph1, replace)
bayesgraph diagnostics {d} {sig2}
Forest Plot with Bayesian Effect Size and 95% Credibility Interval included
bottom of page