*! version 4.2.0 16feb2026

Title

bayeshmc -- Bayesian estimation via Hamiltonian Monte Carlo
using CmdStan

Syntax

Setup (required before first use)

bayeshmc, setup path(cmdstan_path)

Estimation

bayeshmc [, options] : estimation_command

Postestimation

bayeshmc summary

bayeshmc ess

bayeshmc trace [, parameters(namelist) saving(filename)]

bayeshmc ac [, parameters(namelist) lags(#) saving(filename)]

bayeshmc density [, parameters(namelist) saving(filename)]

bayeshmc histogram [, parameters(namelist) saving(filename)]

bayeshmc waic

bayeshmc loo

MCMC
iterations(#) number of post-warmup iterations per chain;
default is iterations(2000)

warmup(#) number of warmup (adaptation) iterations per chain;
default is warmup(1000)

chains(#) number of Markov chains; default is chains(4)

seed(#) random-number seed; default is random

thin(#) thinning interval; default is thin(1)

Sampler tuning
adapt_delta(#) target acceptance rate; default is adapt_delta(0.8).
Increase toward 1.0 for difficult posteriors.

max_treedepth(#) maximum tree depth for NUTS; default is
max_treedepth(10)

Priors
normalprior(#) standard deviation of normal prior on regression
coefficients; default is normalprior(100)

prior_sd(#) synonym for normalprior()

lkjprior(#) LKJ concentration parameter eta for correlation
matrices; default is lkjprior(1) (uniform over correlations)

covprior(string) prior for unstructured covariance matrix;
one of lkj, iw, siw, huangwand, spherical;
default is covprior(lkj).
See Covariance priors below.

Credible intervals
clevel(#) set credible level; default is clevel(95)

hpd report highest posterior density (HPD) intervals instead of
equal-tailed intervals

Execution
parallel run chains in parallel (one process per chain)

threads(#) number of threads per chain for within-chain
parallelism; default is threads(1)

cmdstan(path) override the stored CmdStan installation path

nocache force recompilation of the Stan model even if
unchanged

dryrun write the Stan code and data files but do not sample

verbose display detailed progress messages

Reporting
noheader suppress the output header

saving(filename) save MCMC draws to a Stata dataset

Parameterization
reparam use non-centered parameterization for random effects
(may improve sampling in sparse data)

Description

bayeshmc fits Bayesian regression models using the No-U-Turn Sampler
(NUTS), a variant of Hamiltonian Monte Carlo (HMC), implemented through
CmdStan. It provides a Stata-native interface that automatically generates
Stan code, exports data in JSON format, invokes the CmdStan sampler, and
parses the posterior draws back into Stata for inference and visualization.

The command supports a wide range of single-level, multilevel (mixed-effects),
and panel (random-effects) models. The syntax mirrors Stata's native
estimation commands: users type bayeshmc : followed by a standard Stata
estimation command exactly as they would type it otherwise.

Compared with Stata's built-in bayes: prefix (which uses the
Metropolis--Hastings algorithm), bayeshmc uses gradient-based sampling
via NUTS/HMC, which typically delivers higher effective sample sizes per
iteration, better mixing for complex posteriors, and is the gold standard for
Bayesian computation as implemented in Stan.

Prerequisites. CmdStan must be installed on the user's system.
See mc-stan.org for
installation instructions. After installation, run
bayeshmc, setup path(cmdstan_path) once to register the
path.

Options

MCMC

iterations(#) specifies the number of post-warmup draws per chain.
The total number of draws equals iterations x chains.
Default is 2000.

warmup(#) specifies the number of warmup (adaptation) iterations during
which the sampler tunes its step size and mass matrix. These draws are
discarded. Default is 1000. For complex models (e.g., multilevel or
heteroskedastic), consider increasing to 2000--4000.

chains(#) specifies the number of independent Markov chains.
Running multiple chains enables R-hat convergence diagnostics. Default is 4.

seed(#) sets the random-number seed for reproducibility. If not
specified, a random seed is used.

thin(#) specifies the thinning interval; only every #-th draw is
retained. Thinning is generally unnecessary with HMC but may reduce storage
for very long runs. Default is 1 (no thinning).

Sampler tuning

adapt_delta(#) sets the target average acceptance probability during
warmup adaptation. Higher values (e.g., 0.95 or 0.99) yield smaller step
sizes, reducing divergent transitions at the cost of slower sampling.
Default is 0.8. If you see divergent transition warnings, increase this
value.

max_treedepth(#) sets the maximum tree depth for the NUTS algorithm.
The number of leapfrog steps per iteration is at most 2^max_treedepth.
Default is 10. Increase to 12--15 if you see maximum treedepth warnings.

Priors

normalprior(#) specifies the standard deviation of the normal prior
placed on all regression coefficients: beta_k ~ N(0, #^2). Default is
100, yielding a weakly informative prior. For standardized data, values of
1--10 may be more appropriate. prior_sd(#) is a synonym.

lkjprior(#) specifies the concentration parameter eta of the LKJ prior
on the correlation matrix in multilevel models with unstructured covariance.
When eta = 1 (default), the prior is uniform over all valid correlation
matrices. When eta > 1, the prior increasingly favors the identity matrix
(weak correlations). When eta < 1, the prior favors extreme correlations.
This option is relevant only for covprior(lkj) and
covprior(huangwand).

covprior(string) selects the prior distribution for the unstructured
covariance matrix in multilevel models.
See Covariance priors for details.

Credible intervals

clevel(#) specifies the credible interval level as a percentage.
Default is 95.

hpd requests highest posterior density (HPD) intervals. By default,
equal-tailed intervals are reported. HPD intervals are the shortest intervals
containing the specified probability mass and are preferred for asymmetric
posteriors.

Execution

parallel runs chains in parallel as separate operating-system processes.
This can substantially reduce wall-clock time on multicore machines. Each
chain runs as an independent CmdStan process.

threads(#) specifies the number of OpenMP threads per chain for
within-chain parallelism via Stan's reduce_sum mechanism.
Default is 1.

cmdstan(path) overrides the CmdStan path previously set by
bayeshmc, setup. Useful when multiple CmdStan versions are installed.

nocache forces recompilation of the Stan model. By default,
bayeshmc caches compiled models and reuses them if the Stan code has not
changed.

dryrun writes the Stan code (model.stan) and data
(data.json) files without invoking the sampler. Useful for inspecting
or customizing the generated Stan code.

verbose displays additional progress messages during model compilation
and sampling.

Subcommands

After estimation, the following postestimation subcommands are available.
All operate on the draws from the most recent bayeshmc estimation.

bayeshmc summary displays a table of posterior summaries (mean, SD, MCSE,
median, and credible interval) for all parameters.

bayeshmc ess displays effective sample sizes (ESS) and split-chain R-hat
for all parameters. ESS is computed using the variogram-based estimator of
Vehtari et al. (2021). R-hat uses the rank-normalized split-chain diagnostic.

bayeshmc trace [, parameters(namelist)
saving(filename)]
produces trace plots of MCMC draws by iteration. By default, the first four
parameters are plotted. Specify parameters() to choose specific
parameters using their display names (e.g., parameters(y.sens_ind y._cons)).

bayeshmc ac [, parameters(namelist)
lags(#) saving(filename)]
produces autocorrelation plots. Default is 40 lags.

bayeshmc density [, parameters(namelist)
saving(filename)]
produces kernel density plots of the posterior distributions.

bayeshmc histogram [, parameters(namelist)
saving(filename)]
produces histograms of the posterior distributions.

bayeshmc waic computes the Widely Applicable Information Criterion
(Watanabe 2010) using the pointwise log-likelihood.

bayeshmc loo computes the leave-one-out cross-validation information
criterion using Pareto-smoothed importance sampling (PSIS-LOO; Vehtari,
Gelman, and Gabry 2017).

Supported models

bayeshmc supports the following Stata estimation commands. The Stata
command is typed after the colon exactly as in frequentist estimation.

Continuous outcomes

regress linear regression

tobit tobit (censored) regression

truncreg truncated regression

intreg interval regression

hetregress heteroskedastic regression

betareg beta regression for fractional outcomes

heckman Heckman selection model

glm generalized linear model

Binary outcomes

logit / logistic logistic regression

probit probit regression

cloglog complementary log-log regression

hetprobit heteroskedastic probit

heckprobit probit with sample selection

Count outcomes

poisson Poisson regression

nbreg negative binomial regression

gnbreg generalized negative binomial regression

tpoisson truncated Poisson regression

zip zero-inflated Poisson regression

zinb zero-inflated negative binomial regression

Ordinal outcomes

ologit ordered logistic regression

oprobit ordered probit regression

hetoprobit heteroskedastic ordered probit

Multinomial outcomes

mlogit multinomial logistic regression

Survival outcomes

streg parametric survival regression (Weibull)

Panel (random-effects) models

These require xtset to be set before estimation. A single
random intercept per panel unit is estimated.

xtreg linear panel RE regression

xtlogit logistic panel RE regression

xtprobit probit panel RE regression

xtpoisson Poisson panel RE regression

xtnbreg negative binomial panel RE regression

xtologit ordered logistic panel RE regression

xtoprobit ordered probit panel RE regression

Multilevel (mixed-effects) models

These use the || groupvar: syntax for random effects and support
covariance(unstructured) for correlated random effects.

mixed linear mixed-effects regression

melogit mixed-effects logistic regression

meprobit mixed-effects probit regression

mecloglog mixed-effects complementary log-log

mepoisson mixed-effects Poisson regression

menbreg mixed-effects negative binomial regression

meologit mixed-effects ordered logistic regression

meoprobit mixed-effects ordered probit regression

metobit mixed-effects tobit regression

mestreg mixed-effects parametric survival regression

meglm mixed-effects generalized linear model

mehetregress mixed-effects heteroskedastic regression

mehetoprobit mixed-effects heteroskedastic ordered probit

Multilevel models with the binomial(varname) option are
supported for grouped binomial data (e.g., diagnostic test accuracy
meta-analysis).

Covariance priors for unstructured random effects

For multilevel models with covariance(unstructured), bayeshmc
supports five prior specifications for the random-effects covariance matrix.
These are selected via the covprior() option. All priors produce the
same output parameters (sigma_u, var_u, corr) regardless of
the internal parameterization.

covprior(lkj) -- LKJ decomposition (default)

Decomposes the covariance matrix into standard deviations and correlations:

tau_k ~ half-Cauchy(0, 2.5)

L_Omega ~ LKJ_corr_cholesky(eta)

Sigma = diag(tau) * Omega * diag(tau)

The concentration parameter eta is set via lkjprior(). When eta = 1,
the prior is uniform over all valid correlation matrices. This is the
recommended default and the most widely used prior in the Stan community
(Lewandowski, Kurowicka, and Joe 2009; Stan Development Team 2024).

covprior(iw) -- Inverse-Wishart

The classical conjugate prior:

Sigma ~ Inverse-Wishart(R + 1, I_R)

where R is the dimension (number of random effects) and I_R is the R x R
identity matrix. The degrees of freedom R + 1 yield the least informative
proper Inverse-Wishart. This prior couples the marginal variances and
correlations, which can be undesirable. It tends to shrink the correlation
toward zero and may underestimate variance components
(Alvarez, Niemi, and Simpson 2014).

covprior(siw) -- Scaled Inverse-Wishart

The scaled Inverse-Wishart (Gelman and Hill 2007; O'Malley and Zaslavsky 2008)
ameliorates the coupling problem of the standard IW by introducing marginal
scale parameters:

xi_k ~ half-Cauchy(0, 2.5)

S_raw ~ Inverse-Wishart(R + 1, I_R)

Sigma = diag(xi) * S_raw * diag(xi)

The scale parameters xi allow each variance component to be adjusted
independently before the IW is applied, giving better marginal behavior for
each variance.

covprior(huangwand) -- Huang-Wand

The Huang and Wand (2013) prior uses a scale-mixture representation to place
independent half-t priors on the standard deviations while retaining an LKJ
prior on the correlations:

a_k ~ Inverse-Gamma(1/2, 1/A^2) with A = 2.5

tau_k | a_k ~ half-Normal(0, a_k)

L_Omega ~ LKJ_corr_cholesky(eta)

Sigma = diag(tau) * Omega * diag(tau)

The marginal prior on each tau_k is half-t with heavy tails, providing
robustness for sparse or weakly identified variance components. This prior
is recommended for small-sample multilevel models and diagnostic test accuracy
meta-analysis (Huang and Wand 2013).

covprior(spherical) -- Spherical decomposition

The spherical (angular) decomposition parameterizes the Cholesky factor of
the correlation matrix through hyperspherical coordinates
(Pinheiro and Bates 1996; Rapisarda, Brigo, and Mercurio 2007):

tau_k ~ half-Cauchy(0, 2.5)

theta_m ~ Uniform(0, pi) for m = 1, ..., R(R-1)/2

L = f(theta)

Sigma = diag(tau) * L L' * diag(tau)

Each row of the Cholesky factor is computed from cosines and cumulative
products of sines of the angles. A Jacobian adjustment ensures the implied
prior on the correlation matrix is uniform. For R = 2 (the most common case
in bivariate meta-analysis), there is a single angle theta and the correlation
is simply cos(theta), which is uniform on (-1, 1).

This parameterization avoids the implicit regularization of LKJ and provides
a truly flat prior over the space of positive-definite correlation matrices.

Choosing a covariance prior. In practice, the choice of covariance
prior matters most for small samples, few groups, or weakly identified variance
components. For well-identified models with moderate-to-large numbers of
groups, the five priors typically give very similar results. When in doubt:

Default: covprior(lkj) with lkjprior(1) -- well-tested, the
Stan community standard.

Small samples: covprior(huangwand) -- robust half-t tails prevent
boundary estimates.

Prior sensitivity: Fit the model under two or three priors and compare
posteriors. Agreement indicates robustness.

Truly flat correlation prior: covprior(spherical) -- uniform on
correlation space without the implicit geometry of LKJ.

Diagnostic plots

Graphical MCMC diagnostics are essential for assessing convergence and
identifying sampling problems.

bayeshmc trace shows the time series of draws for each parameter.
A well-mixing chain appears as a "fuzzy caterpillar" with no trends or
periods of stasis.

bayeshmc ac shows the autocorrelation function. Low autocorrelation
indicates efficient sampling. HMC/NUTS typically exhibits much lower
autocorrelation than Metropolis--Hastings.

bayeshmc density shows kernel density estimates of the posterior.
Overlay with prior distributions to assess prior sensitivity.

bayeshmc histogram shows histograms of the posterior draws.

All diagnostic commands accept parameters(namelist) to select
specific parameters using their display names (as shown in the output table).
By default, the first four parameters are plotted. All accept
saving(filename) to save the graph.

Model comparison

bayeshmc waic computes the Widely Applicable Information Criterion:

WAIC = -2 * (lppd - p_waic)

where lppd is the log pointwise predictive density and p_waic is the
effective number of parameters (Watanabe 2010; Gelman, Hwang, and
Vehtari 2014).

bayeshmc loo computes the PSIS-LOO (Pareto-smoothed importance
sampling leave-one-out) cross-validation estimate:

LOO = -2 * sum_i log(p(y_i | y_{-i}))

estimated via importance sampling with Pareto-smoothed weights (Vehtari,
Gelman, and Gabry 2017). PSIS-LOO is generally preferred over WAIC for model
comparison. Pareto k diagnostics indicate reliability: k < 0.5 is good,
0.5 < k < 0.7 is acceptable, k > 0.7 suggests the approximation may be
unreliable.

Stored results

bayeshmc stores the following in e():

Scalars

e(N) number of observations

e(chains) number of chains

e(mcmc_size) total number of posterior draws

Macros

e(cmd) bayeshmc

e(family) estimation command (e.g., melogit)

e(depvar) name of dependent variable

Matrices

e(b) posterior means

e(V) diagonal matrix of posterior variances

Technical notes

How bayeshmc works

bayeshmc operates in five stages:

1. Parse. The Stata estimation command is parsed to identify the model
family, dependent variable, predictors, grouping structure, and covariance
type.

2. Generate. A Stan program is automatically written to the CmdStan
working directory. The program specifies the data block, parameters, priors,
likelihood, and generated quantities (including pointwise log-likelihoods for
WAIC/LOO).

3. Export. Data are exported from Stata to JSON format. Design matrices,
grouping indicators, and trial counts (for binomial models) are included.

4. Sample. CmdStan compiles the Stan program (if not cached) and runs
the NUTS sampler. Each chain produces a CSV file of posterior draws.

5. Parse results. The CSV output files are read back into Mata. Summary
statistics, ESS, and R-hat are computed. Results are displayed in a
Stata-style table and stored in e().

Effective sample size (ESS)

ESS is computed using the variogram-based estimator of Vehtari et al. (2021),
which improves upon the traditional autocorrelation-based estimator by using
the variogram V_t = mean((x_ - x_i)^2) and Geyer's initial monotone
sequence truncation. The estimator is:

rho_hat(t) = 1 - V_t / (2 * var_hat_plus)

ESS = m * n / (1 + 2 * sum_^ rho_hat(t))

where m is the number of chains and n is the draws per chain. For a single
chain, the chain is split in half and treated as two chains.

R-hat convergence diagnostic

R-hat compares within-chain and between-chain variance. Values near 1.0
indicate convergence. R-hat > 1.01 suggests the chains have not converged.
bayeshmc reports the maximum R-hat across all parameters in the header.

Divergent transitions

Divergent transitions indicate that the sampler encountered regions of high
curvature in the posterior that it could not accurately explore. If divergences
occur, increase adapt_delta() toward 0.95 or 0.99, or reparameterize
the model using reparam.

Model caching

Compiled Stan models are cached in the CmdStan working directory. If you
run the same model type with different data, the cached executable is reused,
saving compilation time (typically 30--60 seconds). Use nocache to
force recompilation.

Default priors

bayeshmc uses the following default priors unless overridden:

Parameter Default prior

Regression coefficients (beta) N(0, 100^2)

Intercept (alpha) N(0, 100^2)

Residual SD (sigma) half-Cauchy(0, 5)

Random-effect SDs (tau) half-Cauchy(0, 2.5)

Correlation matrix (Omega) LKJ(eta = 1)

Cutpoints (ordered models) N(0, 10^2)

Overdispersion (phi, nbreg) half-Cauchy(0, 5)

Weibull shape Gamma(1, 1)

Beta regression precision Gamma(0.01, 0.01)

Examples

All examples assume CmdStan has been set up:

. bayeshmc, setup path(C:\Users\username\.cmdstan\cmdstan-2.38.0)

Basic single-level models

Example 1: Linear regression

. sysuse auto, clear

. bayeshmc, iter(2000) chains(4) seed(12345) : regress price mpg weight

. bayeshmc summary

. bayeshmc ess

Example 2: Logistic regression with full diagnostics

. sysuse auto, clear

. bayeshmc, iter(2000) chains(4) seed(12345) : logit foreign mpg weight

. bayeshmc trace

. bayeshmc ac, lags(30)

. bayeshmc density

. bayeshmc histogram

Example 3: Probit regression with tighter prior

. sysuse auto, clear

. bayeshmc, iter(2000) chains(4) normalprior(10) seed(12345) :

probit foreign mpg weight length

Example 4: Ordered logistic regression with HPD intervals

. webuse fullauto, clear

. bayeshmc, iter(2000) chains(4) hpd clevel(90) seed(12345) :

ologit rep77 foreign price mpg

Example 5: Complementary log-log regression

. webuse lbw, clear

. bayeshmc, iter(2000) chains(4) seed(12345) : cloglog low age lwt smoke

Count data models

Example 6: Poisson regression

. webuse airline, clear

. bayeshmc, iter(2000) chains(4) seed(12345) : poisson injuries XYZowned

Example 7: Negative binomial regression

. webuse medpar, clear

. bayeshmc, iter(2000) chains(4) seed(12345) : nbreg los type2 type3 hmo

Example 8: Zero-inflated Poisson

. webuse fish, clear

. bayeshmc, iter(2000) chains(4) seed(12345) :

zip count persons livebait, inflate(child camper)

Example 9: Zero-inflated negative binomial

. webuse fish, clear

. bayeshmc, iter(2000) chains(4) seed(12345) :

zinb count persons livebait, inflate(child camper)

Example 10: Truncated Poisson

. webuse tpoisson1, clear

. bayeshmc, iter(2000) chains(4) seed(12345) : tpoisson count x1 x2, ll(0)

Special continuous models

Example 11: Beta regression for fractional outcomes

. webuse gastinger, clear

. bayeshmc, iter(2000) chains(4) seed(12345) :

betareg fev1rate i.center age

Example 12: Tobit regression

. sysuse auto, clear

. bayeshmc, iter(2000) chains(4) seed(12345) : tobit mpg weight foreign, ll(17)

Example 13: Truncated regression

. webuse truncreg, clear

. bayeshmc, iter(2000) chains(4) seed(12345) : truncreg wage age educ, ll(0)

Example 14: Interval regression

. webuse intregxmpl, clear

. bayeshmc, iter(2000) chains(4) seed(12345) :

intreg wage1 wage2 age nev_mar rural school tenure

Example 15: Heteroskedastic regression

. sysuse auto, clear

. bayeshmc, iter(2000) chains(4) seed(12345) : hetregress price mpg weight, het(foreign)

Example 16: Heteroskedastic probit

. webuse lbw, clear

. bayeshmc, iter(2000) chains(4) seed(12345) : hetprobit low age lwt smoke, het(race)

Selection models

Example 17: Heckman selection model

. webuse womenwk, clear

. bayeshmc, iter(2000) chains(4) seed(12345) :

heckman wage education age, select(married children education age)

Example 18: Heckman probit

. webuse womenwk, clear

. gen byte employed = (wage < .)

. bayeshmc, iter(2000) chains(4) seed(12345) :

heckprobit employed education age, select(married children education age)

Multinomial outcomes

Example 19: Multinomial logistic regression

. sysuse auto, clear

. bayeshmc, iter(2000) chains(4) seed(12345) :

mlogit rep78 price mpg weight foreign

Survival analysis

Example 20: Weibull survival regression

. webuse cancer, clear

. stset studytime, failure(died)

. bayeshmc, iter(2000) chains(4) seed(12345) : streg age drug, distribution(weibull)

Panel (random-effects) models

Example 21: Panel linear RE

. webuse nlswork, clear

. xtset idcode year

. bayeshmc, iter(1000) chains(2) seed(12345) : xtreg ln_wage age ttl_exp tenure

Example 22: Panel logistic RE

. webuse nlswork, clear

. xtset idcode year

. gen byte highwage = (ln_wage > 1.6) if ln_wage < .

. bayeshmc, iter(1000) chains(2) seed(12345) : xtlogit highwage age ttl_exp tenure

Example 23: Panel Poisson RE

. webuse nlswork, clear

. xtset idcode year

. gen hours_cat = round(hours/10) if hours < .

. bayeshmc, iter(1000) chains(2) seed(12345) : xtpoisson hours_cat age ttl_exp tenure

Example 24: Panel ordered logit RE

. webuse nlswork, clear

. xtset idcode year

. egen wage_cat = cut(ln_wage), group(4)

. replace wage_cat = wage_cat + 1

. bayeshmc, iter(1000) chains(2) seed(12345) : xtologit wage_cat age ttl_exp tenure

Multilevel (mixed-effects) models

Example 25: Linear mixed model with random intercept

. webuse nlswork, clear

. bayeshmc, iter(2000) chains(4) seed(12345) :

mixed ln_wage age ttl_exp tenure || idcode:

Example 26: Mixed-effects logistic regression

. webuse bangladesh, clear

. bayeshmc, iter(2000) chains(4) seed(12345) :

melogit c_use urban age child* || district:

Example 27: Mixed-effects Poisson

. webuse epilepsy, clear

. bayeshmc, iter(2000) chains(4) seed(12345) :

mepoisson seizures treat lbas lbas_trt lage visit || subject:

Example 28: Mixed-effects negative binomial

. webuse epilepsy, clear

. bayeshmc, iter(2000) chains(4) seed(12345) :

menbreg seizures treat lbas lbas_trt lage visit || subject:

Diagnostic test accuracy (DTA) meta-analysis

Example 29: Bivariate DTA meta-analysis (LKJ prior, default)

The bivariate random-effects model for DTA meta-analysis (Reitsma et al. 2005;
Chu and Cole 2006) is fitted using melogit with binomial data and
unstructured random effects for sensitivity and specificity:

. use dta_data, clear

. bayeshmc, iter(1000) warmup(2000) chains(4) seed(49830) :

melogit y sens_ind || study: sens_ind,

covariance(unstructured) binomial(n)

. bayeshmc summary

. bayeshmc ess

. bayeshmc trace, parameters(y.sens_ind y._cons sigma_u1._cons sigma_u2._cons)

. bayeshmc density, parameters(y.sens_ind y._cons)

Example 30: DTA meta-analysis with Inverse-Wishart prior

. bayeshmc, iter(1000) warmup(2000) chains(4) covprior(iw) seed(49830) :

melogit y sens_ind || study: sens_ind,

covariance(unstructured) binomial(n)

Example 31: DTA meta-analysis with Scaled Inverse-Wishart prior

. bayeshmc, iter(1000) warmup(2000) chains(4) covprior(siw) seed(49830) :

melogit y sens_ind || study: sens_ind,

covariance(unstructured) binomial(n)

Example 32: DTA meta-analysis with Huang-Wand prior

Recommended for meta-analyses with few studies (< 10):

. bayeshmc, iter(1000) warmup(2000) chains(4) covprior(huangwand) seed(49830) :

melogit y sens_ind || study: sens_ind,

covariance(unstructured) binomial(n)

Example 33: DTA meta-analysis with spherical decomposition prior

Provides a uniform prior on the (-1, 1) correlation space for R = 2:

. bayeshmc, iter(1000) warmup(2000) chains(4) covprior(spherical) seed(49830) :

melogit y sens_ind || study: sens_ind,

covariance(unstructured) binomial(n)

Model comparison

Example 34: Comparing Poisson vs. negative binomial with WAIC and LOO

. webuse medpar, clear

. * Model 1: Poisson

. bayeshmc, iter(2000) chains(4) seed(12345) : poisson los type2 type3 hmo

. bayeshmc waic

. bayeshmc loo

. * Model 2: Negative binomial

. bayeshmc, iter(2000) chains(4) seed(12345) : nbreg los type2 type3 hmo

. bayeshmc waic

. bayeshmc loo

Lower WAIC and LOO indicate better out-of-sample predictive performance.

Execution options

Example 35: Parallel chains for faster estimation

. bayeshmc, iter(4000) warmup(2000) chains(4) parallel seed(12345) :

mixed ln_wage age ttl_exp tenure || idcode:

Example 36: Fixing divergent transitions

If you see divergent transition warnings, increase adapt_delta:

. bayeshmc, iter(2000) chains(4) adapt_delta(0.95) max_treedepth(12)

seed(12345) : melogit y x || group: x, covariance(unstructured)

Example 37: Inspecting the generated Stan code

Use dryrun to examine and optionally customize the Stan program before
sampling:

. bayeshmc, dryrun verbose : melogit y x || group:, cov(unstruct)

The Stan file will be written to
cmdstan_path/bhmc3_work/model.stan and can be inspected with
type or edited manually before rerunning.

Example 38: Specific parameters in diagnostic plots

. bayeshmc, iter(2000) chains(4) seed(12345) : regress price mpg weight foreign

. bayeshmc trace, parameters(price.mpg price.weight)

. bayeshmc ac, parameters(price.mpg) lags(50)

. bayeshmc density, parameters(price.mpg price.weight price._cons)

. bayeshmc histogram, parameters(sigma._cons) saving(sigma_hist)

References

Hamiltonian Monte Carlo and NUTS

Betancourt, M. 2018.
A conceptual introduction to Hamiltonian Monte Carlo.
arXiv:1701.02434.

Duane, S., A. D. Kennedy, B. J. Pendleton, and D. Roweth. 1987.
Hybrid Monte Carlo.
Physics Letters B 195: 216--222.

Hoffman, M. D., and A. Gelman. 2014.
The No-U-Turn Sampler: Adaptively setting path lengths in Hamiltonian Monte
Carlo.
Journal of Machine Learning Research 15: 1593--1623.

Neal, R. M. 2011.
MCMC using Hamiltonian dynamics.
In Handbook of Markov Chain Monte Carlo, ed. S. Brooks, A. Gelman,
G. Jones, and X.-L. Meng, 113--162. Boca Raton, FL: Chapman and Hall/CRC.

Convergence diagnostics

Brooks, S. P., and A. Gelman. 1998.
General methods for monitoring convergence of iterative simulations.
Journal of Computational and Graphical Statistics 7: 434--455.

Gelman, A., and D. B. Rubin. 1992.
Inference from iterative simulation using multiple sequences.
Statistical Science 7: 457--472.

Geyer, C. J. 1992.
Practical Markov chain Monte Carlo.
Statistical Science 7: 473--483.

Vehtari, A., A. Gelman, D. Simpson, B. Carpenter, and P.-C. Buerkner. 2021.
Rank-normalization, folding, and localization: An improved R-hat for assessing
convergence of MCMC (with discussion).
Bayesian Analysis 16: 667--718.

Model comparison

Gelman, A., J. Hwang, and A. Vehtari. 2014.
Understanding predictive information criteria for Bayesian models.
Statistics and Computing 24: 997--1016.

Vehtari, A., A. Gelman, and J. Gabry. 2017.
Practical Bayesian model evaluation using leave-one-out cross-validation and
WAIC.
Statistics and Computing 27: 1413--1432.

Watanabe, S. 2010.
Asymptotic equivalence of Bayes cross validation and widely applicable
information criterion in singular learning theory.
Journal of Machine Learning Research 11: 3571--3594.

Covariance matrix priors

Alvarez, I., J. Niemi, and M. Simpson. 2014.
Bayesian inference for a covariance matrix.
Proceedings of the 26th Annual Conference on Applied Statistics in
Agriculture.

Gelman, A., and J. Hill. 2007.
Data Analysis Using Regression and Multilevel/Hierarchical Models.
Cambridge: Cambridge University Press.

Huang, A., and M. P. Wand. 2013.
Simple marginally noninformative prior distributions for covariance matrices.
Bayesian Analysis 8: 439--452.

Lewandowski, D., D. Kurowicka, and H. Joe. 2009.
Generating random correlation matrices based on vines and extended onion
method.
Journal of Multivariate Analysis 100: 1989--2001.

O'Malley, A. J., and A. M. Zaslavsky. 2008.
Domain-level covariance analysis for multilevel survey data with structured
nonresponse.
Journal of the American Statistical Association 103: 1405--1418.

Pinheiro, J. C., and D. M. Bates. 1996.
Unconstrained parameterizations for variance-covariance matrices.
Statistics and Computing 6: 289--296.

Rapisarda, F., D. Brigo, and F. Mercurio. 2007.
Parameterizing correlations: A geometric interpretation.
IMA Journal of Management Mathematics 18: 55--73.

Diagnostic test accuracy meta-analysis

Chu, H., and S. R. Cole. 2006.
Bivariate meta-analysis of sensitivity and specificity with sparse data: A
generalized linear mixed model approach.
Journal of Clinical Epidemiology 59: 1331--1332.

Dwamena, B. A. 2007.
MIDAS: Stata module for meta-analytical integration of diagnostic test
accuracy studies.
Statistical Software Components S456880, Boston College Department of
Economics.

Dwamena, B. A. 2017.
MIDAS: An update for comprehensive diagnostic test accuracy meta-analysis
in Stata.
Presented at the 2017 Stata Conference, Baltimore, MD.

Macaskill, P. 2004.
Empirical Bayes estimates generated in a hierarchical summary ROC analysis
agreed closely with those of a full Bayesian analysis.
Journal of Clinical Epidemiology 57: 925--932.

Reitsma, J. B., A. S. Glas, A. W. S. Rutjes, R. J. P. M. Scholten,
P. M. Bossuyt, and A. H. Zwinderman. 2005.
Bivariate analysis of sensitivity and specificity produces informative summary
measures in diagnostic reviews.
Journal of Clinical Epidemiology 58: 982--990.

Rutter, C. M., and C. A. Gatsonis. 2001.
A hierarchical regression approach to meta-analysis of diagnostic test
accuracy evaluations.
Statistics in Medicine 20: 2865--2884.

Bayesian methodology

Carpenter, B., A. Gelman, M. D. Hoffman, D. Lee, B. Goodrich, M. Betancourt,
M. Brambilla, J. Guo, P. Li, and A. Riddell. 2017.
Stan: A probabilistic programming language.
Journal of Statistical Software 76(1): 1--32.

Gelman, A., J. B. Carlin, H. S. Stern, D. B. Dunson, A. Vehtari, and
D. B. Rubin. 2013.
Bayesian Data Analysis. 3rd ed.
Boca Raton, FL: Chapman and Hall/CRC.

McElreath, R. 2020.
Statistical Rethinking: A Bayesian Course with Examples in R and Stan.
2nd ed. Boca Raton, FL: Chapman and Hall/CRC.

Stan Development Team. 2024.
Stan User's Guide. Version 2.35.
https://mc-stan.org/users/documentation/

Related software

Buerkner, P.-C. 2017.
brms: An R package for Bayesian multilevel models using Stan.
Journal of Statistical Software 80(1): 1--28.

Gabry, J., R. Cesnovar, and A. Johnson. 2024.
CmdStanR: R Interface to CmdStan.
https://mc-stan.org/cmdstanr/

Grant, R. L. 2017.
The Stata interface for Stan.
Stata Journal 17: 662--670.

Rue, H., S. Martino, and N. Chopin. 2009.
Approximate Bayesian inference for latent Gaussian models by using integrated
nested Laplace approximations.
Journal of the Royal Statistical Society, Series B 71: 319--392.

Authors

Ben A. Dwamena, MD

Clinical Associate Professor Emeritus of Radiology

Division of Nuclear Medicine and Molecular Imaging

University of Michigan, Ann Arbor

bdwamena@umich.edu

The MIDAS suite (midas, midas_mh, bayeshmc, and related
commands) is maintained by B. A. Dwamena. For bug reports and feature
requests, please contact the author.

Also see

Help: bayesmh, mixed, melogit, meprobit,
xtreg, xtlogit

Web: Stan website,
CmdStan documentation