R BEST包如何按组设置贝叶斯先验

Question

我正在使用 R BEST 包来测试两组之间的均值差异，并且我想按组设置先验信念。

在下面的 R 代码中，我可以按组设置先验均值（请参阅priors2），但不能按组设置先验标准差（请参阅priors3）。
难道我做错了什么？

library(BEST)

y1 <- c(5.77, 5.33, 4.59, 4.33, 3.66, 4.48)    
y2 <- c(3.88, 3.55, 3.29, 2.59, 2.33, 3.59)

priors1 <- list(muM = 6, muSD = 2)
BESTout1 <- BESTmcmc(y1, y2, priors=priors1, parallel=FALSE)

priors2 <- list(muM = c(6, 4), muSD = 2)
BESTout2 <- BESTmcmc(y1, y2, priors=priors2, parallel=FALSE)

priors3 <- list(muM = c(6, 4), muSD = c(2, 2))
BESTout3 <- BESTmcmc(y1, y2, priors=priors3, parallel=FALSE)

Answer 1

R BEST 包是两个样本 t 检验的贝叶斯版本(John K. Kruschke) 。

从您提供的示例代码来看，您似乎正在尝试对要比较的两组使用不同的先验（在考虑数据之前合并有关参数的现有知识或信念的方法）。

然而，贝叶斯 t 检验在 BEST 包中实现的方式是，它使用单个共享先验来计算组的标准差。这是因为它假设所比较的两个组可能具有相似的方差。

在您发布的代码中，muM参数指的是组均值的先验平均值，muSD参数指的是组均值的先验标准差。虽然您可以为每个组设置不同的先验均值（如您在 中所做的那样priors2），muSD但参数在组之间共享。

不幸的是，BEST 包不支持为每个组设置不同的先验标准差。
详情请见mikemeredith/BEST第6期

如果您有兴趣为每个组的标准差设置不同的先验，则需要使用不同的包或手动执行分析。

一种替代方案是使用该paul-buerkner/brms包（本文介绍），它是 R 中用于贝叶斯分析的多功能包。它使用Stan（贝叶斯统计平台）进行计算，并支持各种模型和先验。但是，请注意，brms直接使用 或 Stan 可能比使用 BEST 包更复杂，因为它们要求您更详细地指定模型和先验。

贝叶斯 t 检验brms可能如下所示：

library(brms)

# create a data frame
dat <- data.frame(
  y = c(y1, y2),
  group = factor(rep(c("y1", "y2"), each = 6))
)

# fit the model
fit <- brm(
  y ~ group,
  prior = c(
    set_prior("normal(6, 2)", class = "Intercept", coef = "groupy1"),
    set_prior("normal(4, 2)", class = "Intercept", coef = "groupy2")
  ),
  data = dat,
  family = gaussian()
)

此代码创建一个贝叶斯线性回归模型，其中两组的均值被建模为截距，并为每组设置不同的先验。

这只是一个示例，您需要根据您的具体情况选择适当的先验。

您可以在 中设置先验brms。例如：

set_prior("normal(0,5)", class = "b", coef = "x1") 
set_prior("student_t(10, 0, 1)", class = "b", coef = "x2")

在这里，我们为系数 设定平均值为 0、标准差为 5 的正态先验x1，为 设定自由度为 10 的学生 t 先验x2。该代码需要与函数调用结合使用才能brm运行模型。

需要注意的一件关键事情是，coefin 中的参数set_prior对应于数据中的变量名称，而不是组名称。如果您想要建模特定于组的效应，则可能需要考虑在模型中使用组级效应或随机效应。

要将其与模型代码集成，您可以执行以下操作：

library(brms)
prior1 <- set_prior("normal(0,5)", class = "b", coef = "x1") 
prior2 <- set_prior("student_t(10, 0, 1)", class = "b", coef = "x2")
priors <- c(prior1, prior2)
mod_eqvar <- brm(
  IQ ~ Group,
  data = d,
  cores = 4,  # Use 4 cores for parallel processing
  prior = priors
)
summary(mod_eqvar)

（请将x1和替换x2为数据集中您希望设置先验的变量的实际名称。）

x1此代码设置和的指定先验，x2并在模型中使用它们。然后，它使用指定的先验
执行贝叶斯 t 检验，比较数据帧中的IQ不同s。 然后使用 来查看模型的结果。Groupd
summary(mod_eqvar)

为了比较组，上面的初始代码将变为：

# Load the library
library(brms)

# Set a seed for reproducibility
set.seed(42)

# Generate some data
group1 <- rnorm(50, mean = 5, sd = 1)
group2 <- rnorm(50, mean = 6, sd = 1)

# Combine into a data frame
d <- data.frame(
  IQ = c(group1, group2),
  Group = factor(rep(c("Group1", "Group2"), each = 50))
)

# Set priors for the model. This prior represents a prior belief about the difference in means between Group2 and Group1
priors <- set_prior("normal(0, 5)", class = "b", coef = "Group2")

# Fit the model
mod_eqvar <- brm(
  IQ ~ Group,
  data = d,
  prior = priors,
  cores = 4,  # Use 4 cores for parallel processing
  file = "iqgroup"  # Save results into a file
)

# View the results
summary(mod_eqvar)

在此代码中，模型正在比较dataframe 中不同 s 之间的IQ ~ Group变量。在模型公式中，是结果变量，是预测变量。预测变量是分类变量，它代表数据中的不同组。IQGroupdIQGroupGroup

（另请参阅“如何使用 R 中的稳健贝叶斯估计来比较两个组”）

在模型公式中IQ ~ Group，IQ是结果变量，Group是预测变量。预测变量Group 是分类变量，它代表数据中的不同组。

在 R 回归模型中使用因子变量时，因子的第一个水平被视为参考组，其他水平的估计值表示结果变量（在本例中为IQ）相对于参考组的差异。模型中“ ”的估计值表示与参考组（组 1）Group2之间的均值差异。Group2

因此，如果您在模型中为“Group2”的系数设置先验，如以下代码行所示：

# Load the library
library(brms)

# Set a seed for reproducibility
set.seed(42)

# Generate some data
group1 <- rnorm(50, mean = 5, sd = 1)
group2 <- rnorm(50, mean = 6, sd = 1)

# Combine into a data frame
d <- data.frame(
  IQ = c(group1, group2),
  Group = factor(rep(c("Group1", "Group2"), each = 50))
)

# Set priors for the model. This prior represents a prior belief about the difference in means between Group2 and Group1
priors <- set_prior("normal(0, 5)", class = "b", coef = "Group2")

# Fit the model
mod_eqvar <- brm(
  IQ ~ Group,
  data = d,
  prior = priors,
  cores = 4,  # Use 4 cores for parallel processing
  file = "iqgroup"  # Save results into a file
)

# View the results
summary(mod_eqvar)

您正在设定关于 和参考组之间的均值差异的先验信念Group2，而不是关于 的总体均值Group2。然后，当您拟合模型时，会根据数据更新此先验分布，从而生成两组之间均值差异的后验分布。

该summary(mod_eqvar)函数将为您提供模型的摘要，包括组 2 和组 1 之间均值的估计差异。此差异由输出中的“Group2”估计值表示。在贝叶斯背景下，我们通常会查看感兴趣参数的后验分布（在本例中为各组之间的均值差异）和 95% 可信区间。如果均值差异的 95% 可信区间不包括 0，则表明组间均值存在统计显着性差异。

该set_prior函数用于设置模型中各组之间均值差异的先验分布。coefin 中的参数对应set_prior于数据中的变量名称，而不是组名称。因此，当我们使用coef = "Group2"in时set_prior，我们实际上是为 Group2 和参考组 (Group1) 之间的均值差异设置先验，而不是为 Group2 的总体均值设置先验。

（另请参阅“ brms 模型的先前定义”）

在贝叶斯环境中，我们通常查看感兴趣参数的后验分布Group（在本例中为on的影响IQ）和 95%可信区间 ，而不是 p 值。如果可信区间不包含零，则表明该效应在 0.05 水平上具有统计显着性。

关于错误：

Error: The following priors do not correspond to any model parameter: 
  Intercept_groupy1 ~ normal(6, 2) 
  Intercept_groupy2 ~ normal(4, 2) 
Function 'get_prior' might be helpful to you. 

该错误消息表明指定的先验与任何模型参数都不匹配。这可能是因为组名称与数据中的名称不匹配。

在 brms 模型语法中，截距项表示参考组的总体平均值，这些coef项表示与其他组的该参考组平均值的差异。该coef术语应与因子变量中组级别的确切名称相匹配，在您的情况下似乎是Group1和Group2。

让我们更正函数coef中的术语set_prior以匹配数据中的组名称。您可以尝试在函数中将 and 替换Intercept_groupy1为Intercept_groupy2and Intercept_Group1（Intercept_Group2或因子变量级别的确切名称）set_prior，如下所示：

# Set priors for the model
priors <- c(
  set_prior("normal(6, 2)", class = "Intercept", coef = "Group1"),
  set_prior("normal(4, 2)", class = "Intercept", coef = "Group2")
)

# Fit the model
mod_eqvar <- brm(
  IQ ~ Group,
  data = d,
  prior = priors,
  cores = 4,  # Use 4 cores for parallel processing
  file = "iqgroup"  # Save results into a file
)

注意：这可能仍然不起作用，因为 brms 可能不允许同一模型中因子变量的不同级别具有不同的先验。在这种情况下，您可能需要为每个组拟合单独的模型，然后手动比较结果。