A Simulation Study On The Robustness Of Bayesian Structural Equation Modeling Under Small Samples, Heavy Tails, and Collinearity

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 \hat{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.