R Bayes

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  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  0.05
pct_supported 
  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 
  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  0)
)

Odds Ratios (for logistic models)

# Convert log-odds to odds ratios
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  0)

# P(|effect| > some threshold)
prob_meaningful  0.1)

Compare Effect Magnitudes

# Test if within-person effect is larger than between-person
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 
  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  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  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 
  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

1. Define causal DAG with dagitty
2. Validate DAG against data with localTests()
3. Identify adjustment sets for target effects
4. Specify priors based on domain knowledge
5. Fit brms model with random effects for nested data
6. Check diagnostics (convergence, PPCs)
7. Extract posteriors for inference
8. Compute marginal effects on interpretable scale
9. 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)

Source & license

This open-source skill is cataloged on AgentStack and links to its original source — we do not rehost the code.

• Author: ab604
• Source: ab604/claude-code-r-skills
• License: MIT

Install and usage instructions live in the source repository linked above.