Groups keyboard shortcuts have been updated
Dismiss
See shortcuts

Setting prior when using prefixed labels in mediation models

39 views
Skip to first unread message

Yiping Zhang

unread,
May 14, 2024, 12:25:19 PM5/14/24
to blavaan
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

Ed Merkle

unread,
May 14, 2024, 12:29:38 PM5/14/24
to blavaan
Yiping,

For priors on individual parameters, you use the "prior()" modifier inside the model syntax. An example is here:


I like the idea of using the parameter labels to set priors outside of the model syntax. I will try to add this in the future.

Ed

Yiping Zhang

unread,
May 14, 2024, 12:40:52 PM5/14/24
to Ed Merkle, blavaan
Hi Ed, 
Thank you so much for your prompt response!
I tried to set prior inside the model syntax, but the output returns “normal(0,10)” Here is the non-clustered simulated data sample and code I use and a screen shot of the result for your reference: 
library(brms)
library(easystats)
library(bayestestR)


# Sample size
n<-300
set.seed(7)

# Simulate categorical predictor x
x <- sample(c("a", "b", "c"), n, replace = TRUE, prob = c(1/3, 1/3, 1/3))

# Create a numeric version of x for easier calculations

x_numeric <- as.numeric(factor(x, levels = c("a", "b", "c")))

# Define effects of x on z
effect_x_on_z <- c(0, 0.5, 1)  # Different effects for each category of x

# Simulate mediator z (normal distribution)

z <- rnorm(n, mean = effect_x_on_z[x_numeric], sd = 1)

# Define total effect of x on y
total_effect_x_on_y <- c(0, 1, 2)  # Different effects for each category of x


# Calculate direct and indirect effects
indirect_effect <- total_effect_x_on_y * 0.30
direct_effect <- total_effect_x_on_y - indirect_effect

# Simulate outcome variable y
# Using x_numeric directly for effect indexing, ensuring it fits the intended direct and indirect effect vectors
y <- rnorm(n, mean = direct_effect[x_numeric] + 0.5 * z, sd = 1)  # 0.5 is the effect of z on y

# Create a dataframe
data <- data.frame(x, z, y)

# Ensure 'x' is factor and has levels specified in case they aren't from data directly
data$x <- factor(data$x, levels = c("a", "b", "c"))

# Create dummy variables for all categories
data$x_a <- ifelse(data$x == "a", 1, 0)
data$x_b <- ifelse(data$x == "b", 1, 0)
data$x_c <- ifelse(data$x == "c", 1, 0)

library(blavaan)


mediation_model <- '
  # Mediation model with explicitly defined priors on paths and variances
  # Path specifications with priors on slopes
  z ~ 1 + a2*(x_b*prior("normal(0, 4)")) + a3*(x_c*prior("normal(0, 4)"))  # Normal(0, 0.5) priors with precision (1/0.25 = 4)
  y ~ 1 + b*(z*prior("normal(0, 4)")) + c2*(x_b*prior("normal(0, 4)")) + c3*(x_c*prior("normal(0, 4)"))  # Normal(0, 0.5) priors on regression coefficients

  # Residual variances with exponential priors
  z ~~ prior("gamma(3,3)[sd]")*z  # Exponential(1) prior for variance of z
  y ~~ prior("gamma(3,3)[sd]")*y  # Exponential(1) prior for variance of y


  # 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
'



#Fit the Bayesian SEM model with correctly specified custom priors
bmodel <- bsem(
  mediation_model,
  data = data,
  meanstructure = TRUE
)

summary(bmodel)
Screenshot 2024-05-14 at 12.37.57 PM.png

--
You received this message because you are subscribed to a topic in the Google Groups "blavaan" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/blavaan/GSVnykn16Ek/unsubscribe.
To unsubscribe from this group and all its topics, send an email to blavaan+u...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/blavaan/2b96a3e0-2709-4076-bcec-86bbe790613dn%40googlegroups.com.


--
Yiping Zhang
Associate Teaching Professor
University of Rhode Island
Address: Swan Hall 141, 60 Upper College Road
Tel: 401-874-2008   Email: yiping...@uri.edu

Ed Merkle

unread,
May 14, 2024, 1:01:11 PM5/14/24
to Yiping Zhang, blavaan
I am away from the computer right now, but I would guess that it is confused by the parentheses in your model syntax. Try putting prior() at the start of each term and/or outside parentheses.

Ed

Yiping Zhang

unread,
May 14, 2024, 3:36:06 PM5/14/24
to Ed Merkle, blavaan
Thank you!! I tried to put prior at the start of each term and inside of the parentheses (outside of the parentheses didn't work). The output is still the default prior (0, 10).

Mauricio Garnier-Villarreal

unread,
May 15, 2024, 5:21:44 AM5/15/24
to blavaan
Try writing the equation twice, in one you ad the label, and in the second one you add the prior.

Or adding the same predictor twice in the same regression, adding the label in first term, and the prior in the second one.

basically the idea is to separate the prior and label

Ed Merkle

unread,
May 20, 2024, 4:24:28 PM5/20/24
to Yiping Zhang, blavaan
Yiping,

Following up on Mauricio's post: below is some syntax that finds the correct priors. But also: note that the prior specification in Stan uses standard deviation instead of precision (which is different from JAGS). So if you want a precision of 4, you need to specify prior("normal(0, .5)").

Ed


mediation_model <- '
# Mediation model with explicitly defined priors on paths and variances
# Path specifications with priors on slopes
z ~ 1 + a2*x_b + prior("normal(0, 4)")*x_b + a3*x_c + prior("normal(0, 4)")*x_c # Normal(0, 0.5) priors with precision (1/0.25 = 4)
y ~ 1 + b*z + prior("normal(0, 4)")*z + c2*x_b + prior("normal(0, 4)")*x_b + c3*x_c + prior("normal(0, 4)")*x_c # Normal(0, 0.5) priors on regression coefficients

# Residual variances with exponential priors
z ~~ prior("gamma(3,3)[sd]")*z # Exponential(1) prior for variance of z
y ~~ prior("gamma(3,3)[sd]")*y # Exponential(1) prior for variance of y

# 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
'



Yiping Zhang

unread,
May 20, 2024, 9:26:21 PM5/20/24
to Ed Merkle, mauga...@gmail.com, blavaan
Thank you so much!!
Reply all
Reply to author
Forward
0 new messages