Bayesian SEM Priors: Moving Beyond Defaults

1. BSEM and the Usefulness of Priors

Bayesian Structural Equation Modeling (BSEM) allows for the formal integration of prior knowledge into model estimation (Depaoli 2021).

Rather than relying solely on the likelihood of the current data, Bayesian updating creates a posterior distribution: Posterior ∝ Likelihood × Prior

Why this matters:
- Stabilizes estimation in complex models or small samples.
- Prevents overfitting through regularization.
- Accumulates scientific knowledge across studies.

2. A 5-Indicator Unidimensional Model

We will use a 5-item unidimensional scale (SIJS) to demonstrate prior specification.

00-data-setup.R

pacman::p_load(blavaan, dplyr, bayesplot, ggplot2, semPlot)

# Import the SIJS Dataset
ds <- rio::import('https://ndownloader.figshare.com/files/22299075') |> 
  select(paste0("SIJS", 1:5)) |> 
  janitor::clean_names()

# Base model syntax
model <- 'js =~ sijs1 + sijs2 + sijs3 + sijs4 + sijs5'

2. A 5-Indicator Unidimensional Model

3. Default (Flat) Priors

By default, blavaan uses diffuse, weakly informative priors (e.g., normal(0, 10) for factor loadings) (Merkle and Rosseel 2016). This allows the data to completely dominate the posterior.

01-default-priors.R

# Estimate model using standard defaults
set.seed(123)
fit_def <- bcfa(model, data = ds, n.chains = 3, bcontrol = list(cores = 3))

# Extract posterior samples 
samp_def <- blavInspect(fit_def, "mcmc")

4. Results: Default Priors

mcmc_areas(samp_def, pars = "js=~sijs2", prob = 0.95) +
  ggtitle("Posterior Density: Default Diffuse Prior", 
          subtitle = "Loading of sijs2 on sijs") +
  theme_minimal(base_size = 16) +
  theme(plot.title = element_text(face = "bold"),
        panel.grid.minor = element_blank())

5. Regularizing Priors

We can apply small-variance priors (shrinkage) to parameters to anchor estimates toward zero. This acts as a penalty against overfitting, similar to Ridge or Lasso regression.

02-regularization.R

# Syntax with regularizing priors on specific loadings
model_reg <- '
  js =~ sijs1 + 
          prior("normal(0, 0.5)")*sijs2 + 
          prior("normal(0, 0.5)")*sijs3 + 
          prior("normal(0, 0.5)")*sijs4 + 
          prior("normal(0, 0.5)")*sijs5
'

set.seed(123)
fit_reg <- bcfa(model_reg, data = ds, n.chains = 3, bcontrol = list(cores = 3))
samp_reg <- blavInspect(fit_reg, "mcmc")

6. Results: Regularizing Priors

mcmc_areas(samp_reg, pars = "js=~sijs2", prob = 0.95) +
  ggtitle("Posterior Density: Regularizing Prior", 
          subtitle = "Loading of sijs2 on sijs (shrunk toward zero)") +
  theme_minimal(base_size = 16) +
  theme(plot.title = element_text(face = "bold"),
        panel.grid.minor = element_blank())

7. Subjective / Informative Priors

Instead of arbitrary shrinkage, we can elicit expert knowledge using the SHELF protocol (O’Hagan 2019) or extract effect sizes from meta-analyses (O’Hagan et al. 2006).

03-objective-priors.R

library(SHELF)

# Expert elicitation for sijs2 loading:
# "25% sure it is < 0.40, 50% sure < 0.65, 75% sure < 0.85"
fit_dist <- fitdist(vals = c(0.40, 0.65, 0.85), 
                    probs = c(0.25, 0.50, 0.75), 
                    lower = 0, upper = 1)

mu <- round(fit_dist$Normal[1], 2)
sigma <- round(fit_dist$Normal[2], 2)

model_obj <- sprintf('
  js =~ sijs1 + prior("normal(%s, %s)")*sijs2 + 
          sijs3 + sijs4 + sijs5
', mu, sigma)

set.seed(123)
fit_obj <- bcfa(model_obj, data = ds, n.chains = 3, bcontrol = list(cores = 3))
samp_obj <- blavInspect(fit_obj, "mcmc")

8. Results: Objective Priors

mcmc_areas(samp_obj, pars = "js=~sijs2", prob = 0.95) +
  ggtitle(sprintf("Posterior Density: Objective Prior ~ N(%s, %s)", mu, sigma), 
          subtitle = "Loading of sijs2 on sijs") +
  theme_minimal(base_size = 16) +
  theme(plot.title = element_text(face = "bold"),
        panel.grid.minor = element_blank())

9. Always Conduct a Sensitivity Analysis

A single prior specification is rarely sufficient. Methodological rigor requires sensitivity analysis to test how varying the prior distributions impacts the posterior inferences (Kaplan 2023).

Key steps:
- Test a range of plausible prior specifications (diffuse, regularizing, objective).
- Compare the posterior means, medians, and highest posterior density (HPD) intervals across models.
- If conclusions change drastically depending on the prior, the data carries less weight than the prior, which must be transparently reported.

10. References

Depaoli, Sarah. 2021. Bayesian structural equation modeling. The Guilford Press. https://www.guilford.com/books/Bayesian-Structural-Equation-Modeling/Sarah-Depaoli/9781462547746/contents.
Kaplan, David. 2023. Bayesian statistics for the social sciences. 2nd ed. Guilford Press. https://www.guilford.com/books/Bayesian-Statistics-for-the-Social-Sciences/David-Kaplan/9781462553549.
Merkle, Edgar C., and Yves Rosseel. 2016. blavaan: Bayesian structural equation models via parameter expansion. http://arxiv.org/abs/1511.05604v2.
O’Hagan, Anthony. 2019. Expert knowledge elicitation: Subjective but scientific.” The American Statistician 73 (sup1): 69–81. https://doi.org/10.1080/00031305.2018.1518265.
O’Hagan, Anthony, Caitlin E. Buck, Alireza Daneshkhah, et al. 2006. Uncertain judgements: Eliciting experts’ probabilities. Wiley. https://doi.org/10.1002/0470033312.