Capturing Direct, Indirect, and Mediation Effects with Bayesian Statistics - Private-Projects237/Statistics GitHub Wiki
This wikipage will be recalculating estimates from the examples used in Intermediate Path Analysis Examples in R. The goal is to understand how to capture direct, indirect, and mediation effects using Bayesian statistics. Specifically, the brm() function from the 'brms' package. To validate that the estimate we are obtaining are actually representing what we think they are, we will be comparing the outputs specifically to the examples from the path analysis wiki.
# These are the same examples from the 'Intermediate Path Analysis Examples.R'
# Load required packages
library(brms)
library(tidyverse)
library(MASS)
library(kableExtra)
# Set seed for reproducibility
set.seed(123)
# Set population parameters
b0 = .5
b1 = -.05
b2 = -.06
b3 = -.04
sig = .12 # more appropriate for logistic model
# Step 1: Generate synthetic data
n <- 5000 # Sample size
# Create a correlation matrix (moderate correlations)
cov_matrix <- matrix(c(
1.0, 0.5, 0.4,
0.5, 1.0, 0.55,
0.4, 0.55, 1.0
), nrow = 3)
# Generate predictors from multivariate normal distribution
predictors <- mvrnorm(n, mu = c(0, 0, 0), Sigma = cov_matrix)
# Create data frame with new variable names
data <- data.frame(
x1 = predictors[, 1], # anxiety
x2 = predictors[, 2], # depression
x3 = predictors[, 3] # stress
)
# Generate a linear predictor with noise
data$y <- b0 + b1 * data$x1 + b2 * data$x2 + b3 * data$x3 + rnorm(n, mean =0, sd = sig)
# Force the outcome to be between 0 and 100
data$y <- pmax(0, pmin(1, data$y ))
# Run the Bayesian Model
mul_bay <- brm(
y ~ x1 + x2 + x3, # Same formula as lm()
data = data,
family = gaussian(), # Continuous outcome (default link = identity)
chains = 4,
cores = 4,
seed = 123
)
# Model output summary
summary(mul_bay)
# Extract the parameters
round(fixef(mul_bay),2) # fixed effects
round(VarCorr(mul_bay)$residual__$sd,2) # sigma
round(bayes_R2(mul_bay), 2) # r^2
check_collinearity(mul_bay) # Checks multicollinearity
suppressMessages(suppressWarnings(
standardize_parameters(mul_bay, method = "refit", verbose = FALSE))) # Standardized betas- removes annoying message
There is a lot of good information that is presented from the model summary. Starting at the top, we see that the model was y ~ x1 + x2 + x3, which is the multiple regression we intended to model. Then we see that there are 5000 observations within the dataset, which is correct. Moving to the 'Regression Coefficients', here we see all of the unstandardized beta coefficients and the y-intercept. There are no p-values but there are clearly 95% CI and we can see that they along with the standard errors are very very small, indicating that we can have a lot of confidence that are our estimates represent the true population estimates (which they do!). The in the 'Further Distributional Parameters', we see that the residual standard error is also present. Thus the bayesian model was able to reproduce the fixed effects along with the one random effect for the errors. In comparison to the path analysis model, we see that the estimates are identical!
| Model Summary | Path Analysis Model (Unstandardized) |
|---|---|
 |
 |
There are very cool functions that we can use to extract specific information from our Bayesian model. We can use functions like fixef() and VarCorr() to extract the unstandardized estimates. bayes_R2() will produce the coefficient of determination, which is nice (unsure if its for the outcome or for the full model), there is a test for vif through check_collinearity(), which is awesome! Lastly, we can convert our regressions into standardized beta coefficients using standardize_parameters(), and we see they also match that of the path analysis example.
| Extracting Info From Model | Path Analysis (Standardized ) |
|---|---|
 |
 |
# Load required packages
library(tidyverse)
library(brms)
library(MASS)
library(kableExtra)
library(parameters) # standardize_parameters
library(performance) # check_collinearity
# Set seed for reproducibility
set.seed(123)
# Set population parameters
b0 = .5
b1 = -.05
b2 = -.06
b3 = -.04
sig = .12 # more appropriate for logistic model
# Step 1: Generate synthetic data
n <- 5000 # Sample size
# Create a correlation matrix (moderate correlations)
cov_matrix <- matrix(c(
1.0, 0.75, 0.80,
0.75, 1.0, 0.83,
0.80, 0.83, 1.0
), nrow = 3)
# Generate predictors from multivariate normal distribution
predictors <- mvrnorm(n, mu = c(0, 0, 0), Sigma = cov_matrix)
# Create data frame with new variable names
data <- data.frame(
x1 = predictors[, 1], # anxiety
x2 = predictors[, 2], # depression
x3 = predictors[, 3] # stress
)
# Generate a linear predictor with noise
data$y <- b0 + b1 * data$x1 + b2 * data$x2 + b3 * data$x3 + rnorm(n, mean =0, sd = sig)
# Force the outcome to be between 0 and 100
data$y <- pmax(0, pmin(1, data$y ))
# Run the Bayesian Model
mul2_bay <- brm(
y ~ x1 + x2 + x3, # Same formula as lm()
data = data,
family = gaussian(), # Continuous outcome (default link = identity)
chains = 4,
cores = 4,
seed = 123
)
# Model output summary
summary(mul2_bay)
# Extract the parameters
round(fixef(mul2_bay),2) # fixed effects
round(VarCorr(mul2_bay)$residual__$sd,2) # sigma
round(bayes_R2(mul2_bay), 2) # r^2
check_collinearity(mul2_bay) # Checks multicollinearity
suppressMessages(suppressWarnings(
standardize_parameters(mul2_bay, method = "refit", verbose = FALSE))) # Standardized betas- removes annoying message
We can skip the model output this time and directly focus on the information extracted from the model. We can see that we have extracted the unstandardized beta coefficients and the residual standard error along with
| Extracting Info From Model | Path Analysis (Standardized ) |
|---|---|
 |
 |
Example 3: Introducing a Mediator (No Multicollinearity)
# Load required packages
library(tidyverse)
library(brms)
library(MASS)
library(kableExtra)
library(tidybayes) # For spread_draws(), mean_qi()
# Set seed for reproducibility
set.seed(123)
# Set population parameters
b0_Y = .6
b0_M = .4
b1 = -.05 # Keep these the same (beta estimates; they are the total effect)
b2 = -.06
b3 = -.04
a1 = -.07 # pathways a
a2 = -.06
a3 = -.03
b = 0.2 # pathway b (only 1)
e_sig = .08
y_sig = .08 # more appropriate for this type of model
c_prime1 = b1 - (a1 * b)
c_prime2 = b2 - (a2 * b)
c_prime3 = b3 - (a3 * b)
# Step 1: Generate synthetic data
n <- 5000 # Sample size
# Create a correlation matrix (moderate correlations)
cov_matrix <- matrix(c(
1.0, 0.5, 0.4,
0.5, 1.0, 0.55,
0.4, 0.55, 1.0
), nrow = 3)
# Generate predictors from multivariate normal distribution
predictors <- mvrnorm(n, mu = c(0, 0, 0), Sigma = cov_matrix)
# Create data frame with new variable names
data <- data.frame(
x1 = predictors[, 1], # anxiety
x2 = predictors[, 2], # depression
x3 = predictors[, 3] # stress
)
# Generate a Mediator variable
data$m <- b0_M + a1 * data$x1 + a2* data$x2 + a3 * data$x3 + rnorm(n, 0, sd = e_sig)
# Generate a linear predictor with noise + Mediator
data$y <- b0_Y + b1 * data$x1 + b2 * data$x2 + b3 * data$x3 + b * data$m + rnorm(n, 0, sd = y_sig)
# Force the outcome to be between 0 and 1
data$y <- pmax(0, pmin(1, data$y))
bayes_mediation <- brm(
bf(m ~ x1 + x2 + x3) + # Path a: for all predictors
bf(y ~ x1 + x2 + x3 + m) + # Paths c' for all predictors + mediator
set_rescor(FALSE), # No residual correlation; SUPER IMPORTANT!
data = data,
family = gaussian(),
chains = 4,
cores = 4,
seed = 123
)
# Model output summary
summary(bayes_mediation)
# Extract the parameters
round(fixef(bayes_mediation),2) # fixed effects
round(VarCorr(bayes_mediation)$residual__$sd,2) # sigma
round(bayes_R2(bayes_mediation), 2) # r^2
# Check multicollinearity manually
lm1 <- lm(m ~ x1 + x2 + x3, data = data)
lm2 <- lm(y ~ x1 + x2 + x3 + m, data = data)
check_collinearity(lm1)
check_collinearity(lm2)
# Check the standardized beta coefficients
# Rerun the model fully standardized
data_std <- data %>% mutate(across(c(x1, x2, x3, m, y), scale))
bayes_mediation_std <- suppressMessages(suppressWarnings(
update(bayes_mediation, newdata = data_std)))
round(fixef(bayes_mediation_std),2)
This time let's look at the model output- we see that this is a mediation analysis and a more complex one since we have three variables (x1, x2, x3) predicting a mediator 'm' and then those same predictors plus 'm' are being used to predict the outcome 'y'. Below we see the unstandardized beta coefficients, which include the relationships between x and m, x and y, and the b pathway which is m and y. Below that we see the residual error distribution for both endogenous variables 'm' and 'y'- extremely interesting is that the standard deviation for both of these variables is exactly as what was programmed! Additionally all of the unstandardized beta coefficients are identical to what was produced by the path analysis version.
| Model Output |
|---|
 |
Okay so we needed to make some changes here since the previous functions we were using no longer work with these more complex models. To check vif, it is now recommended to run multiple regressions for each endogenous variables and check the variance inflation factor through the vif() function. The mediator is getting pretty high but I think that is okay and it makes sense considering that it is a composite of the three predictors. Additionally, to get the standardized beta coefficients, it is know recommended to use the scale() function on all the variables to standardized them. Aftewards, rerun the model using the update() function, this will produced the standardized estimates that we can then extract. We see that also from the bayes_R2() function that
Extracting Info From Model + lm() models |
Path Analysis (Standardized ) |
|---|---|
 |
 |
We are not done just yet- we should also look at the indirect effects possible through the mediation in the model. This next part, of the code to produce it is very tedious in Bayesian. We basically need to use the spread_draws() function to obtain thousands of simulated estimates, those for the slope of a predictor (x1, x2, x3) and the mediator ('m'), which are our pathway a's, and the slope of 'm' with 'y', which is our pathway b. Once these values are drawn we need to multiply them together to produce the indirect effect and then do some additional algebra to calculate the total proportion of mediation. we then use the median_qi() function to produce the means of these calculations plus 95% confidence intervals so we know if they are statistically significant or not. Afterwards we display this information.
# Extract relevant posterior draws
draws <- bayes_mediation_std %>%
spread_draws(b_m_x1, b_m_x2, b_m_x3, b_y_m, b_y_x1, b_y_x2, b_y_x3)
# Compute indirect effects and mediation proportions
indirect_effects <- draws %>%
mutate(
indirect_x1 = b_m_x1 * b_y_m,
indirect_x2 = b_m_x2 * b_y_m,
indirect_x3 = b_m_x3 * b_y_m,
total_x1 = indirect_x1 + b_y_x1,
total_x2 = indirect_x2 + b_y_x2,
total_x3 = indirect_x3 + b_y_x3,
prop_med_x1 = indirect_x1 / total_x1,
prop_med_x2 = indirect_x2 / total_x2,
prop_med_x3 = indirect_x3 / total_x3
)
# Summarize the indirect effects and mediation proportions
indirect_summary <- indirect_effects %>%
median_qi(indirect_x1, indirect_x2, indirect_x3,
prop_med_x1, prop_med_x2, prop_med_x3) %>%
as.data.frame() %>%
dplyr::select(-.width, -.point, -.interval) %>%
round(3) %>% c()
# Print results as transposed table
indirect_summary_matrix <- matrix(indirect_summary, nrow = 2,byrow = TRUE)
row.names(indirect_summary_matrix) <- c("Indirect Effect", "Mediation")
colnames(indirect_summary_matrix) <- sub(".*_", "", names(indirect_summary))[1:(length(indirect_summary)/2)]
indirect_summary_matrix
As we can see from below- the indirect effects and the mediation proportion match. Also, these are the standardized indirect effects not the unstandardized ones.
| Indirect effects | Indirect Effects Path Analysis |
|---|---|
 |
 |