如何在R中对含缺失单元格的贝叶斯重复测量ANOVA进行分析？

1. 安装并加载所需包 #

首先安装分析需要的工具包：

install.packages(c("tidyverse", "brms", "emmeans"))

加载包：

library(tidyverse)
library(brms)
library(emmeans)

2. 数据准备 #

2.1 创建数据框 #

将提供的数据导入R：

data_wide <- tibble(
  id = 1:12,
  group = c("B", "A", "A", "C", "B", "C", "A", "B", "C", "B", "A", "C"),
  m_c = c(NA, 1.15, 6.92, NA, NA, NA, 5.38, NA, NA, NA, 5.88, NA),
  m_s = c(NA, NA, NA, 4.38, NA, 2.42, NA, NA, 6.2, NA, NA, 2.21),
  m_v = c(4.97, NA, NA, NA, 2.71, NA, NA, 5.1, NA, 6.73, NA, NA),
  a_c = c(2.58, NA, NA, NA, 3.98, NA, NA, 1.37, NA, 1.76, NA, NA),
  a_s = c(NA, 2.26, 1.24, NA, NA, NA, 3.06, NA, NA, NA, 3.45, NA),
  a_v = c(NA, NA, NA, 5.43, NA, 5.64, NA, NA, 6.15, NA, NA, 5.65),
  e_c = c(NA, NA, NA, 6.83, NA, 6.15, NA, NA, 5.97, NA, NA, 5.95),
  e_s = c(2.9, NA, NA, NA, 5.72, NA, NA, 6.26, NA, 6.93, NA, NA),
  e_v = c(NA, 6.18, 4.14, NA, NA, NA, 6.85, NA, NA, NA, 5.01, NA)
)

2.2 转换为长格式 #

贝叶斯重复测量模型需要长格式数据，将宽格式转换为长格式：

data_long <- data_wide %>%
  pivot_longer(
    cols = m_c:e_v,
    names_to = "condition",
    values_to = "happiness"
  ) %>%
  separate(condition, into = c("time_abbr", "flavor_abbr"), sep = "_") %>%
  mutate(
    time_of_day = recode(time_abbr, m = "Morning", a = "Afternoon", e = "Evening"),
    flavor = recode(flavor_abbr, c = "Chocolate", s = "Strawberry", v = "Vanilla")
  ) %>%
  select(id, group, time_of_day, flavor, happiness)

3. 拟合贝叶斯重复测量ANOVA模型 #

使用brms包拟合包含组、时间段、口味及其交互项的模型，同时包含参与者的随机效应以控制重复测量：

full_model <- brm(
  formula = happiness ~ group + time_of_day * flavor + (1 | id),
  data = data_long,
  family = gaussian(),
  chains = 4,
  iter = 2000,
  warmup = 1000,
  cores = 4,
  seed = 123
)

查看模型摘要：

summary(full_model)

4. 检验各效应的显著性（贝叶斯因子） #

通过比较不同模型的贝叶斯因子来检验时间段、口味及交互项的效应：

4.1 检验时间段效应 #

拟合不含时间段的模型并计算贝叶斯因子：

no_time_model <- brm(
  formula = happiness ~ group + flavor + (1 | id),
  data = data_long,
  family = gaussian(),
  chains = 4,
  iter = 2000,
  warmup = 1000,
  cores = 4,
  seed = 123
)

bf_time <- bayes_factor(full_model, no_time_model)
print(bf_time)

贝叶斯因子>3表示中等强度支持包含时间段的模型，即时间段对幸福感有显著影响；>10表示强支持。

4.2 检验口味效应 #

拟合不含口味的模型并计算贝叶斯因子：

no_flavor_model <- brm(
  formula = happiness ~ group + time_of_day + (1 | id),
  data = data_long,
  family = gaussian(),
  chains = 4,
  iter = 2000,
  warmup = 1000,
  cores = 4,
  seed = 123
)

bf_flavor <- bayes_factor(full_model, no_flavor_model)
print(bf_flavor)

4.3 检验交互效应 #

拟合不含交互项的模型并计算贝叶斯因子：

no_interaction_model <- brm(
  formula = happiness ~ group + time_of_day + flavor + (1 | id),
  data = data_long,
  family = gaussian(),
  chains = 4,
  iter = 2000,
  warmup = 1000,
  cores = 4,
  seed = 123
)

bf_interaction <- bayes_factor(full_model, no_interaction_model)
print(bf_interaction)

5. 事后检验 #

使用emmeans包进行事后 pairwise 比较，结果包含后验均值和95%可信区间：

5.1 时间段的事后检验 #

posthoc_time <- emmeans(full_model, pairwise ~ time_of_day, adjust = "holm")
print(posthoc_time, digits = 3)

5.2 口味的事后检验 #

posthoc_flavor <- emmeans(full_model, pairwise ~ flavor, adjust = "holm")
print(posthoc_flavor, digits = 3)

5.3 交互项的事后检验 #

posthoc_interaction <- emmeans(full_model, pairwise ~ time_of_day:flavor, adjust = "holm")
print(posthoc_interaction, digits = 3)

说明 #

• brms包通过MCMC自然处理缺失数据，无需额外插补，解决了JASP中遇到的缺失数据问题。
• 事后检验使用Holm方法调整多重比较，可信区间不包含0表示两组差异显著。

内容的提问来源于stack exchange，提问作者justAnewbie