Bayesian Structural Equation Modeling (BSEM) allows for the formal integration of prior knowledge into model estimation (Lee, 1981).
Rather than relying solely on the likelihood of the current data, Bayesian updating creates a posterior distribution: \(p(\theta \mid y) \propto p(y \mid \theta)\, p(\theta)\)
Why this matters:
- Stabilizes estimation in complex models or small samples.
- Prevents overfitting through regularization.
- Accumulates scientific knowledge across studies.
Unidimensional (5-indicator) Model
We will use a 5-item unidimensional scale (SIJS, Sinval & Marôco, 2020) to demonstrate prior specification.
By default, blavaan uses diffuse, weakly informative priors (e.g., \(\mathcal{N}(0, 10)\) for factor loadings \(\left[\lambda\right]\))) (Merkle & Rosseel, 2016). This allows the data to dominate the posterior.
01-default-priors.R
# Estimate model using standard defaultsset.seed(1987)fit_def <-bcfa(model, data = ds, std.ov = T, std.lv=T, n.chains =3, bcontrol =list(cores =3))# Extract posterior samples samp_def <-blavInspect(fit_def, "mcmc")
Default Priors
\(Std. Est.\)
\(HPD 95\%\)
\(\hat{R}\)
\(N_{eff}\)
Hyperparameters
Autocorrelation
Trace
Density
\(\lambda_{11}\)
0.59
(0.36; 0.76)
1
3368.01
\(\mathcal{N}(0, 10)\)
\(\lambda_{21}\)
0.95
(0.84; 1.00)
1
1920.45
\(\mathcal{N}(0, 10)\)
\(\lambda_{31}\)
0.57
(0.31; 0.75)
1
2694.65
\(\mathcal{N}(0, 10)\)
\(\lambda_{41}\)
0.85
(0.71; 0.94)
1
2126.40
\(\mathcal{N}(0, 10)\)
\(\lambda_{51}\)
0.58
(0.34; 0.74)
1
2785.95
\(\mathcal{N}(0, 10)\)
Regularizing Priors
We can apply small-variance priors (shrinkage) to parameters to anchor estimates toward zero. In blavaan, we can use the dpriors() function to specify regularizing priors for specific parameters. For example, we can apply a \(\mathcal{N}\left(0, .1\right)\) to the \(\lambda\).
02-regularization.R
# Syntax with regularizing priors on specific loadingsmodel_reg <-'js =~ sijs1 + sijs2 + sijs3 + sijs4 + sijs5'reg_priors <-dpriors(lambda ="normal(0, .1)")set.seed(1987)fit_reg <-bcfa(model_reg, data = ds, dp = reg_priors, std.ov = T, std.lv=T, n.chains =3, bcontrol =list(cores =3))
Regularizing Priors
\(Std. Est.\)
\(HPD 95\%\)
\(\hat{R}\)
\(N_{eff}\)
Hyperparameters
Autocorrelation
Trace
Density
\(\lambda_{11}\)
0.13
( 0.01; 0.36)
1
1787.44
\(\mathcal{N}(0, .1)\)
\(\lambda_{21}\)
0.17
(-0.19; 0.65)
1
1044.77
\(\mathcal{N}(0, .1)\)
\(\lambda_{31}\)
0.10
(-0.18; 0.36)
1
1454.53
\(\mathcal{N}(0, .1)\)
\(\lambda_{41}\)
0.15
(-0.18; 0.59)
1
1166.04
\(\mathcal{N}(0, .1)\)
\(\lambda_{51}\)
0.10
(-0.17; 0.37)
1
1462.21
\(\mathcal{N}(0, .1)\)
Informative Priors
Instead of shrinkage, we can extract empirical effect sizes from published studies, systematic reviews, or Meta-Analytic Structural Equation Modeling (MASEM) studies (Depaoli, 2021; vandeSchoot2021?). This allows us to directly build upon cumulative scientific knowledge.
For this example, suppose a prior MASEM study examining Job Satisfaction indicators reported that the standardized factor loading for sijs2 is 0.70 with a standard error of 0.08(Depaoli, 2021). In blavaan using the Stan backend, we translate this directly into a normal prior where the second parameter is treated as the standard deviation.
03-informative-priors.R
# Parameters extracted from prior literature/MASEM:mu <-0.70sigma <-0.08# Injecting the literature-based informative prior directly into the model syntaxmodel_obj <-sprintf('js =~ sijs1 + prior("normal(%s, %s)")*sijs2 + sijs3 + sijs4 + sijs5', mu, sigma)set.seed(1987)fit_inf <-bcfa(model_obj, data = ds, std.ov =TRUE, std.lv =TRUE, n.chains =3, bcontrol =list(cores =3))
Informative Priors
Informative Priors
\(Std. Est.\)
\(HPD 95\%\)
\(\hat{R}\)
\(N_{eff}\)
Hyperparameters
Autocorrelation
Trace
Density
\(\lambda_{11}\)
0.55
(0.32; 0.72)
1
3056.07
\(\mathcal{N}(0, 10)\)
\(\lambda_{21}\)
0.89
(0.77; 1.00)
1
2972.27
\(\mathcal{N}(0.7, 0.08)\)
\(\lambda_{31}\)
0.53
(0.29; 0.71)
1
3705.56
\(\mathcal{N}(0, 10)\)
\(\lambda_{41}\)
0.83
(0.69; 0.94)
1
2004.17
\(\mathcal{N}(0, 10)\)
\(\lambda_{51}\)
0.53
(0.31; 0.71)
1
3540.42
\(\mathcal{N}(0, 10)\)
Subjective Priors
We can translate expert beliefs into parametric distributions using the SHELF protocol (O’Hagan, 2019). For standardized factor loadings, we ask experts to provide their \(P_{25}\), \(P_{50}\) (median), and \(P_{75}\) quantiles. This captures both the expected magnitude of the loading and the uncertainty around it.
03-shelf-elicitation.R
library(SHELF)# Columns represent the 25%, 50%, and 75% quantiles of the expert's belief.# Example: For sijs1, the expert's median belief is 0.70, and they are 50% sure the true standardized loading falls between 0.50 and 0.85.elicited_vals <-matrix(c(0.50, 0.70, 0.85, # sijs1: Standard prior0.40, 0.65, 0.85, # sijs2: Greater uncertainty (wider interval)0.60, 0.75, 0.90, # sijs3: Expected to be a strong indicator0.45, 0.60, 0.80, # sijs4: Expected to be a weaker indicator0.55, 0.70, 0.85# sijs5: Standard prior ), byrow =TRUE, ncol =3)# Define the probabilities corresponding to the elicited quantilesprobs <-c(0.25, 0.50, 0.75)# Extract Normal hyperparameters into a structured data framepriors <-data.frame(mu =numeric(5), sigma =numeric(5))for(i in1:5) {# lower=0 and upper=1 truncate the distribution for standardized limits fit <-fitdist(vals = elicited_vals[i, ], probs = probs, lower =0, upper =1) priors$mu[i] <-round(fit$Normal[1], 2) priors$sigma[i] <-round(fit$Normal[2], 2)}
Subjective Priors
We inject the elicited hyperparameters directly into specific parameter paths within the model syntax.
04-shelf-model.R
# Construct syntax dynamically using the data frame columnssyntax_items <-sprintf('prior("normal(%.2f, %.2f)")*sijs%d', priors$mu, priors$sigma, 1:5)model_shelf <-paste0("js =~ ", paste(syntax_items, collapse =" + "))cat(model_shelf)set.seed(1987)fit_shelf <-bcfa(model_shelf, data = ds, std.ov =TRUE, std.lv =TRUE, n.chains =3, bcontrol =list(cores =3))
\(Std. Est.\)
\(HPD 95\%\)
\(\hat{R}\)
\(N_{eff}\)
Hyperparameters
Autocorrelation
Trace
Density
\(\lambda_{11}\)
0.58
(0.39; 0.73)
1
3249.16
\(\mathcal{N}(0.69, 0.26)\)
\(\lambda_{21}\)
0.94
(0.84; 1.00)
1
1848.68
\(\mathcal{N}(0.64, 0.33)\)
\(\lambda_{31}\)
0.59
(0.41; 0.73)
1
2820.68
\(\mathcal{N}(0.75, 0.22)\)
\(\lambda_{41}\)
0.83
(0.71; 0.92)
1
2332.28
\(\mathcal{N}(0.61, 0.26)\)
\(\lambda_{51}\)
0.59
(0.41; 0.73)
1
3555.12
\(\mathcal{N}(0.70, 0.22)\)
Sensitivity Analysis
A single prior specification is rarely sufficient. Sensitivity analysis is required to test how varying the prior distributions impacts the posterior inferences (Kaplan, 2023).
Key steps:
- Test a range of plausible prior specifications (diffuse, regularizing, expert-elicited [subjective]).
- Compare the posterior means, medians, and credible intervals (e.g., \(HPD\)) across models.
- If conclusions change drastically depending on the prior, the data carries less weight than the prior, which must be transparently reported.
Sinval, J., & Marôco, J. (2020). Short Index of Job Satisfaction: Validity evidence from Portugal and Brazil. PLoS ONE, 15(4), 1–21. https://doi.org/10.1371/journal.pone.0231474
