Hello,
I am trying to fit a BSEM on ordinal data. How can I set priors given the labels for the mediation analysis to fully use the strength of the Bayesian approach?
Here is the code I used to simulate data and fit the mediation model:
library(blavaan)
# Sample size
n <- 300
set.seed(7)
# Simulate cohort membership (1 to 4)
cohort <- sample(1:4, n, replace = TRUE)
# Simulate categorical predictor x
x <- sample(c("a", "b", "c"), n, replace = TRUE, prob = c(1/3, 1/3, 1/3))
# Convert x to numeric for effect calculation
x_numeric <- as.numeric(factor(x, levels = c("a", "b", "c")))
# Simulate mediator z based on effects of x
effect_x_on_z <- c(0, 0.5, 1)
z <- rnorm(n, mean = effect_x_on_z[x_numeric], sd = 1)
# Calculate direct and indirect effects of x on y
total_effect_x_on_y <- c(0, 1, 2)
indirect_effect <- total_effect_x_on_y * 0.30
direct_effect <- total_effect_x_on_y - indirect_effect
# Simulate random effects for cohorts
sigma_between <- sqrt(0.25) # sqrt of the calculated variance between groups
cohort_effect <- rnorm(4, mean = 0, sd = sigma_between) # 4 cohorts
cohort_effect <- cohort_effect[cohort] # map cohort effect to individuals
# Simulate outcome variable y incorporating cohort effects
y <- rnorm(n, mean = direct_effect[x_numeric] + 0.5 * z + cohort_effect, sd = 1)
# Create dummy variables for x
x_a <- ifelse(x == "a", 1, 0)
x_b <- ifelse(x == "b", 1, 0)
x_c <- ifelse(x == "c", 1, 0)
# Create a dataframe including all variables
data <- data.frame(cohort = factor(cohort), x_a, x_b, x_c, z, y)
# Specify two-level SEM model with random intercepts for z and y mediated by dummy variables of x
mediation_model <- '
level: within
z ~ a2*x_b + a3*x_c
y ~ c2*x_b + c3*x_c + b*z
level: between
z ~ 1
y ~ 1
# Indirect effects (a*b)
a2b := a2*b
a3b := a3*b
# Total effects
total2 := c2 + a2b
total3 := c3 + a3b
# Total proportion mediated
total_indirect := a2b + a3b
total_total := total2 + total3
total_proportion_mediated := total_indirect / total_total
'
bmodel <- bsem(mediation_model, data = data,meanstructure = TRUE,cluster = "cohort")
summary(bmodel)
Thank you in advance!!!
Yiping