brmsで多項ロジットモデルのベイズ推定(交互作用あり)

brmsパッケージによるベイズモデリングについてはdas_kino先生のページで詳しく解説されています。

# パッケージロード -----------

library(brms)

# サンプルデータ ------------

N <- 15

dat <- data.frame(

  y1 = rbinom(N, 10, 0.3), y2 = rbinom(N, 10, 0.5), 

  y3 = rbinom(N, 10, 0.7), x1 = rnorm(N), x2 = rnorm(N)

)

dat$size <- with(dat, y1 + y2 + y3)

dat$y <- with(dat, cbind(y1, y2, y3))

# サンプリング ------------

fit <- brm(bf(y | trials(size)  ~  x1*x2), data = dat, 

           family = multinomial())

# 結果 -------------

summary(fit)

# 可視化 ---------

conditions = brms::make_conditions(fit, "x1")

brms::conditional_effects(fit, effects = "x2", conditions = conditions, categorical = TRUE)

ggplot2を用いて以下のように交互作用の可視化もできます

 便宜上, 変数名は隠しています