Core Packages

library(brms)
library(cmdstanr)
library(dagitty)
library(ggdag)
library(marginaleffects)
library(tidybayes)
library(bayesplot)

Directed Acyclic Graphs (DAGs)

Prior to causal inference, create and validate DAGs with dagitty and ggdag.

Define DAG Structure

dag <- dagitty('
dag {
  # Node positions for visualization
  exposure [pos="0,1"]
  mediator [pos="1,1"]
  outcome [pos="2,1"]
  confounder [pos="1,0"]

  # Edges (arrows)
  confounder -> exposure
  confounder -> outcome
  exposure -> mediator
  mediator -> outcome
  exposure -> outcome
}
')

Identify Adjustment Sets

# For direct effect
adjustmentSets(dag, exposure = "treatment", outcome = "outcome", effect = "direct")

# For total effect
adjustmentSets(dag, exposure = "treatment", outcome = "outcome", effect = "total")

Validate DAG Against Data

# Get implied conditional independencies
implied_cis <- impliedConditionalIndependencies(dag)

# Test against data
ci_results <- localTests(dag, data = analysis_data, type = "cis")

# Assess validation
ci_df <- as.data.frame(ci_results)
ci_df$independent <- ci_df$p.value > 0.05
pct_supported <- 100 * mean(ci_df$independent, na.rm = TRUE)

cat(sprintf("DAG support: %.1f%% of implied CIs hold\n", pct_supported))

Visualize DAG

dag_tidy <- tidy_dagitty(dag)

ggplot(dag_tidy, aes(x = x, y = y, xend = xend, yend = yend)) +
  geom_dag_edges(edge_colour = "grey50") +
  geom_dag_point(size = 20) +
  geom_dag_text(size = 3.5, color = "black") +
  theme_dag() +
  labs(title = "Causal DAG")

Bayesian Regression with brms

Standard Configuration

options(mc.cores = 4)

# Standard brms model call
model <- brm(
  formula = outcome ~ predictor1 + predictor2 + (1 | group_id),
  data = model_data,
  family = bernoulli(link = "logit"),  # For binary outcomes
  prior = priors,
  sample_prior = "yes",  # For prior-posterior comparison
  chains = 4,
  cores = 4,
  iter = 4000,
  warmup = 1000,
  control = list(
    adapt_delta = 0.95,
    max_treedepth = 15
  ),
  seed = 123,  # Set seed for reproducibility
  backend = "cmdstanr",
  file = "models/model_name",         # Cache compiled model
  file_refit = "on_change"            # Only refit if formula/data change
)

Priors

Store priors separately and define explicitly:

priors <- c(
  prior(normal(0, 2), class = "Intercept"),
  prior(normal(0, 1), class = "b"),                    # Fixed effects
  prior(exponential(1), class = "sd"),                 # Random effect SD
  prior(lkj(2), class = "cor")                         # Correlation priors
)

# Get default priors for a formula
get_prior(outcome ~ predictor + (1 | id), data = data, family = bernoulli())

Common Families

# Binary outcome
family = bernoulli(link = "logit")

# Count data
family = poisson(link = "log")
family = negbinomial(link = "log")

# Continuous
family = gaussian()
family = student()  # Robust to outliers

# Ordinal
family = cumulative(link = "logit")

Multilevel Models

Random Intercepts

# Random intercept per participant
outcome ~ predictors + (1 | participant_id)

Random Slopes

# Random intercept and slope for time
outcome ~ time + predictors + (1 + time | participant_id)

Crossed Random Effects

# Participants nested in groups, items crossed
response ~ predictors + (1 | participant_id) + (1 | item_id)

Within-Person Centering

For longitudinal data, separate between-person and within-person effects:

# Create person-centered variables
model_data <- data |>
  group_by(participant_id) |>
  mutate(
    # Between-person means (stable trait)
    predictor_mean = mean(predictor, na.rm = TRUE),

    # Within-person deviations (dynamic change)
    predictor_dev = predictor - predictor_mean,

    # Volatility (person-level SD)
    predictor_sd = sd(predictor, na.rm = TRUE)
  ) |>
  ungroup() |>
  # Standardize
  mutate(
    predictor_mean_z = scale(predictor_mean)[, 1],
    predictor_dev_z = scale(predictor_dev)[, 1]
  )

# Model with both components
model <- brm(
  outcome ~ predictor_mean_z + predictor_dev_z + (1 | participant_id),
  data = model_data,
  family = bernoulli()
)

Lagged Predictors for Temporal Precedence

# Create lagged predictors within person
model_data <- data |>
  group_by(participant_id) |>
  arrange(time) |>
  mutate(
    # Lagged values (from previous timepoint)
    predictor_lag = lag(predictor, order_by = time),
    predictor_dev_lag = lag(predictor_dev, order_by = time)
  ) |>
  ungroup()

# Test if t-1 predicts outcome at t (establishes temporal precedence)
model_lagged <- brm(
  outcome ~ predictor_dev_lag_z + predictor_mean_z + (1 | participant_id),
  ...
)

Extracting and Interpreting Results

Extract Posterior Samples

posterior <- as_draws_df(model)

# Access specific parameter
samples <- posterior$b_predictor_z

# Summary statistics
tibble(
  estimate = median(samples),
  lower_95 = quantile(samples, 0.025),
  upper_95 = quantile(samples, 0.975),
  lower_80 = quantile(samples, 0.10),
  upper_80 = quantile(samples, 0.90),
  prob_negative = mean(samples < 0),
  prob_positive = mean(samples > 0)
)

Odds Ratios (for logistic models)

# Convert log-odds to odds ratios
effects_df <- effects_df |>
  mutate(
    OR = exp(estimate),
    OR_lower = exp(lower_95),
    OR_upper = exp(upper_95)
  )

Posterior Probability of Direction

# P(effect is protective)
prob_protective <- mean(posterior$b_predictor < 0)

# P(effect is harmful)
prob_harmful <- mean(posterior$b_predictor > 0)

# P(|effect| > some threshold)
prob_meaningful <- mean(abs(posterior$b_predictor) > 0.1)

Compare Effect Magnitudes

# Test if within-person effect is larger than between-person
diff <- abs(posterior$b_predictor_dev_z) - abs(posterior$b_predictor_mean_z)
prob_within_larger <- mean(diff > 0)

cat(sprintf("P(|within| > |between|) = %.1f%%\n", 100 * prob_within_larger))

Marginal Effects with marginaleffects

Average Marginal Effects (AME)

# Change in P(outcome) per 1 unit change in predictor
ame <- avg_slopes(
  model,
  variables = c("predictor1_z", "predictor2_z"),
  type = "response"  # Probability scale
)

print(ame)

Predictions at Specific Values

# Predictions at low (-1 SD), mean (0), and high (+1 SD)
predictions <- predictions(
  model,
  newdata = datagrid(
    model = model,
    predictor_z = c(-1, 0, 1)
  ),
  type = "response",
  re_formula = NA  # Population-level (ignore random effects)
)

as.data.frame(predictions) |>
  select(predictor_z, estimate, conf.low, conf.high)

Marginal Effect Plots

plot_predictions(
  model,
  by = "predictor_z",
  type = "response",
  re_formula = NA
) +
  labs(
    title = "Effect of Predictor on Outcome",
    x = "Predictor (standardized)",
    y = "P(Outcome)"
  ) +
  scale_y_continuous(labels = scales::percent) +
  theme_minimal()

Comparing Slopes Across Models

# Extract AME from multiple models
ame_model1 <- avg_slopes(model1, variables = "predictor_z", type = "response")
ame_model2 <- avg_slopes(model2, variables = "predictor_z", type = "response")

comparison <- bind_rows(
  as.data.frame(ame_model1) |> mutate(model = "Full"),
  as.data.frame(ame_model2) |> mutate(model = "Simple")
)

Model Diagnostics

Check MCMC Convergence

# Trace plots
mcmc_trace(model, pars = c("b_Intercept", "b_predictor_z"))

# R-hat (should be < 1.01)
summary(model)$fixed$Rhat

# Effective sample size (should be > 400)
summary(model)$fixed$Bulk_ESS
summary(model)$fixed$Tail_ESS

Posterior Predictive Checks

pp_check(model)
pp_check(model, type = "stat", stat = "mean")
pp_check(model, type = "stat_2d", stat = c("mean", "sd"))

Prior-Posterior Comparison

# Requires sample_prior = "yes" in brm()
prior_summary(model)

# Plot prior vs posterior
mcmc_areas(model, pars = "b_predictor_z", prob = 0.95)

tidybayes for Posterior Manipulation

# Extract draws in tidy format
draws <- model |>
  spread_draws(b_predictor1_z, b_predictor2_z) |>
  mutate(
    OR_predictor1 = exp(b_predictor1_z),
    OR_predictor2 = exp(b_predictor2_z)
  )

# Summarize
draws |>
  median_qi(OR_predictor1, OR_predictor2, .width = c(0.80, 0.95))

# Visualize
draws |>
  ggplot(aes(x = OR_predictor1)) +
  stat_halfeye() +
  geom_vline(xintercept = 1, linetype = "dashed") +
  labs(x = "Odds Ratio", y = NULL)

Workflow Summary

•Define causal DAG with dagitty
•Validate DAG against data with localTests()
•Identify adjustment sets for target effects
•Specify priors based on domain knowledge
•Fit brms model with random effects for nested data
•Check diagnostics (convergence, PPCs)
•Extract posteriors for inference
•Compute marginal effects on interpretable scale
•Visualize effects with uncertainty

Anti-Patterns to Avoid

# WRONG: Using contemporaneous predictors when temporal order matters
outcome_t ~ predictor_t  # Shows co-occurrence, not temporal precedence

# CORRECT: Use lagged predictors to establish temporal precedence
outcome_t ~ predictor_t_minus_1

# WRONG: Ignoring clustering
brm(outcome ~ predictor, data = longitudinal_data)

# CORRECT: Account for repeated measures
brm(outcome ~ predictor + (1 | participant_id), data = longitudinal_data)

# WRONG: Interpreting within-person effects from between-person variation
# Using person aggregates when you have time-varying data

# CORRECT: Person-mean centering to separate effects
outcome ~ predictor_mean_z + predictor_dev_z + (1 | id)