mcpField
R
r-bayes

r-bayes

github.com/ab604データ・データベース21058 回インストール6 週間前に更新
リポジトリ

r-bayes について

ベイズモデリングや因果推論に取り組む R ユーザー向けのスキルです。DAG の設計・検証からマルチレベルモデルの推定、診断、効果の解釈までを扱い、縦断データと事後分布の可視化も取り上げます。

主な機能

  1. 01因果 DAG の作成と検証
  2. 02事前分布の設定と brms モデルの推定
  3. 03マルチレベル・縦断モデルの指定
  4. 04事後分布と限界効果の抽出
  5. 05MCMC 収束と予測適合の確認

活用シーン

  1. 01因果分析の調整集合を特定
  2. 02ランダム効果を含む回帰モデルの推定
  3. 03個人内・個人間の変化を分析
  4. 04モデルを診断し効果を解釈

インストール

コマンドをコピーしてターミナルで実行します。

$npx skills add https://github.com/ab604/claude-code-r-skills --skill r-bayes

インストール後、Claude Code で/r-bayesを入力して読み込まれたか確認します。

詳細

コアパッケージ

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

有向非巡回グラフ（DAG）

因果推論の前に、dagitty と ggdag を使って DAG を作成し、検証します。

DAG の構造を定義する
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
}
')
調整集合を特定する
# For direct effect
adjustmentSets(dag, exposure = "treatment", outcome = "outcome", effect = "direct")

# For total effect
adjustmentSets(dag, exposure = "treatment", outcome = "outcome", effect = "total")
DAG をデータと照合して検証する
# 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))
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")

brms によるベイズ回帰

標準設定
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 <- 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())
よく使われる分布族
# 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")

マルチレベルモデル

ランダム切片
# Random intercept per participant
outcome ~ predictors + (1 | participant_id)
ランダム傾き
# Random intercept and slope for time
outcome ~ time + predictors + (1 + time | participant_id)
交差ランダム効果
# Participants nested in groups, items crossed
response ~ predictors + (1 | participant_id) + (1 | item_id)

個人内中心化

縦断データでは、個人間効果と個人内効果を分けます。

# 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()
)
時間的先行性を考慮したラグ付き予測変数
# 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),
  ...
)

結果の抽出と解釈

事後サンプルを抽出する
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)
)
オッズ比（ロジスティックモデルの場合）
# Convert log-odds to odds ratios
effects_df <- effects_df |>
  mutate(
    OR = exp(estimate),
    OR_lower = exp(lower_95),
    OR_upper = exp(upper_95)
  )
方向に関する事後確率
# 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)
効果の大きさを比較する
# 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))

marginaleffects による限界効果

平均限界効果（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 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)
限界効果のプロット
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()
モデル間で傾きを比較する
# 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")
)

モデル診断

MCMC の収束を確認する
# 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
事後予測チェック
pp_check(model)
pp_check(model, type = "stat", stat = "mean")
pp_check(model, type = "stat_2d", stat = c("mean", "sd"))
事前分布と事後分布の比較
# Requires sample_prior = "yes" in brm()
prior_summary(model)

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

事後分布の操作に使う tidybayes

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

ワークフローの概要

  1. dagitty で因果 DAG を定義する
  2. localTests() でDAG をデータと照合して検証する
  3. 対象の効果に対する調整集合を特定する
  4. ドメイン知識に基づいて事前分布を指定する
  5. ネストされたデータではランダム効果を含めてbrms モデルを適合する
  6. 診断を確認する（収束、PPC）
  7. 推論のために事後分布を抽出する
  8. 解釈しやすい尺度で限界効果を計算する
  9. 不確実性を含めて効果を可視化する

避けるべきアンチパターン

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

コメント

Markdown 対応