日本認知心理学会 研究法研究部会 第8回研究会(2026.03.14)

「ベータ分布を用いた割合・確率データのモデリング」サンプルコード

#パッケージ読み込み

library(brms)

library(tidyverse)

library(palmerpenguins)

#ベータ分布の描画

x <- seq(0.001, 0.999, length=1000)

par1 <- data.frame(mu=c(0.1, 0.25, 0.5, 0.75, 0.9),

                  phi=c(10, 10, 10, 10, 10))

df1 <- do.call(rbind, lapply(1:nrow(par1), function(i){

  mu <- par1$mu[i]

  phi <- par1$phi[i]

  alpha <- mu*phi

  beta <- (1-mu)*phi

  data.frame(x=x,

             density=dbeta(x, alpha, beta),

             mu=mu,

             phi=phi)

}))

ggplot(df1, aes(x, density, color=factor(mu))) +

  geom_line(linewidth=1) +

  labs(x=NULL, color="mu")

par2 <- data.frame(mu=c(0.5, 0.5, 0.5),

                   phi=c(5, 10, 50))

df2 <- do.call(rbind, lapply(1:nrow(par2), function(i){

  mu <- par2$mu[i]

  phi <- par2$phi[i]

  alpha <- mu*phi

  beta <- (1-mu)*phi

  data.frame(x=x,

             density=dbeta(x, alpha, beta),

             mu=mu,

             phi=phi)

}))

ggplot(df2, aes(x, density, color=factor(phi))) +

  geom_line(linewidth=1) +

  labs(x=NULL, color="phi")

#brmsによるベータ回帰

#palmerpenguinsパッケージのpenguinデータを使用

pen <- na.omit(penguins)

#くちばしの長さに対する太さの割合（bill_ratio）変数を作成

pen$bill_ratio <- pen$bill_depth_mm/pen$bill_length_mm

#くちばしの太さの割合を種ごとに比較

model1 <- brm(bill_ratio ~ species,

  data = pen,

  family = Beta(),

  chains = 4,

  iter = 1000,

  seed = 0314)

summary(model1)

plogis(-0.10) #Adelieの平均μ

plogis(-0.10-0.39) #Chinstrapの平均μ

plogis(-0.10-0.67) #Gentooの平均μ

model1_est <- conditional_effects(model1, effects="species")

model1_est

model1_est$species

#種間の差の事後分布を図示

nd <- data.frame(species = levels(pen$species))

mu_post <- posterior_epred(model1, newdata = nd)

mu_diff <- data.frame(Chinstrap_vs_Adelie = mu_post[,2] - mu_post[,1],

                      Gentoo_vs_Adelie = mu_post[,3] - mu_post[,1])

mu_diff |>

  pivot_longer(everything()) |>

  ggplot(aes(x = value)) +

  geom_density(fill="skyblue", alpha=.5) +

  facet_wrap(~name) +

  geom_vline(xintercept = 0, linetype="dashed") +

  labs(x = "Difference in μ", y = "Posterior density")

#μをオッズ比に変換して図示

mu_logit <- posterior_linpred(model1, newdata = nd)

OR <- data.frame(Chinstrap_vs_Adelie = exp(mu_logit[,2] - mu_logit[,1]),

                 Gentoo_vs_Adelie = exp(mu_logit[,3] - mu_logit[,1]))

OR |>

  pivot_longer(everything()) |>

  ggplot(aes(x = value)) +

  geom_density(fill="orange", alpha=.5) +

  facet_wrap(~name) +

  geom_vline(xintercept = 1, linetype="dashed") +

  labs(x = "odds rate", y = "Posterior density")

#使わないか...

#くちばしの太さの割合を翼の長さで予測

pen$flipper_length_z <- scale(pen$flipper_length_mm)

model2 <- brm(bill_ratio ~ flipper_length_z,

  data = pen,

  family = Beta(),

  chains = 4,

  iter = 1000,

  seed = 0314)

summary(model2)

model2_est <- conditional_effects(model2)

model2_est