Bayesian Structural Equation Modeling (BSEM) is increasingly used in health and social research, yet its operating characteristics under district-level constraints remain under-documented when three difficulties co-occur: small sample sizes, heavy-tailed errors, and strong correlations among exogenous constructs. We evaluate BSEM robustness using a Monte Carlo simulation with a fixed full-mediation SEM comprising four latent variables (two exogenous constructs, a mediator, and an outcome) and three reflective indicators per construct. A balanced 2×2×2 design varies sample size (n ∈ 22, 75), error family (Normal vs. Student-t with ν = 5, variance-matched), and exogenous correlation (ρ ∈ 0. 30, 0. 80). For each scenario, R = 50 independent datasets are generated and fitted using MCMC (Stan/NUTS via blavaan). We summarize global fit via posterior predictive p-values (PPP), sampling quality via R and effective sample size (ESS), and sampler diagnostics including divergent transitions, alongside parameter recovery via Monte Carlo bias, RMSE, and 95% credible-interval coverage for structural paths, mediated effects, and ρ. Across conditions, increasing n improves both parameter recovery and sampling behavior. The most fragile settings occur when n = 22 and ρ is high, where the parallel paths X1→M and X2→M become weakly identified and heavy-tailed errors can further degrade precision. These results provide practical guidance for applying BSEM to district-level studies with limited sample sizes.
Parung et al. (2026) studied this question.