我正在从事一个科学出版项目（在社会科学领域），我正在用研究 1、研究 2 等格式构建我的手稿。基本上，我有一个令人信服的发现，其中一个数据集本身在我的领域中可能足以发表论文。但是，我知道如果将其复制到单独的数据集中，我的同事会认为这会更有趣，所以这就是我正在做的事情。
因此，让我们将此作为一个更具体的示例（这是生成的数据，但在统计上与我的真实数据相似）。我预测一个人每周食用的甜甜圈的平均数量与体重增加有关。出于问题的目的，让我们假设这种影响独立于我们通常在更复杂的模型中控制的所有其他人口统计数据。
下面是数据生成过程：
set.seed(3)
donuts1 <- runif(100, 0, 15)
weight1 <- 5 * donuts1 + rnorm(100, 0, 50) + 185
正如预期的那样，跑步lm(weight1 ~ donuts1)揭示了甜甜圈消费对体重的显着影响。
伟大的！现在我将去另一个人群收集更多数据。虽然在这个新人群中效果是恒定的，但随机误差更大。
set.seed(8)
donuts2 <- runif(100, 0, 15)
weight2 <- 5 * donuts2 + rnorm(100, 0, 100) + 185
现在，如果我将我的研究 2 报告为lm(weight2 ~ donuts2)，我们将看到 p 值约为 0.22，因此无法复制。现在，一个理性的人可以观察研究 2 的数据并说，尽管统计上不显着，但研究 2 的数据与研究 1 中的数据是一致的。不过，观察它并说“足够接近”通常不是我们想要的当涉及到统计推断时。
在我看来，有几个选择：
- 将两个样本合二为一，并对组合数据进行分析。这在我的领域和许多其他领域会被认为是奇怪的。
- 像对待荟萃分析一样对待它们，按精度分析效应大小和加权。
- 将研究 1 的信息作为贝叶斯方法分析研究 2 的信息先验。
我对最后一个选项特别感兴趣，即使只是为了演示——我正在学习贝叶斯分析，想看看我是否理解正确。
那么问题来了：我如何将研究 1 的结果作为研究 2 分析的信息先验？
我倾向于这样做，下面用在 JAGS 中运行分析的代码进行了演示。
library(rjags)
library(runjags)
model.inform <- '
model {
for (i in 1:N) {
weight2[i] ~ dnorm(y.hat[i], tau)
y.hat[i] <- a + b * donuts2[i]
}
a ~ dnorm(183.41, 9.11)
b ~ dnorm(5.07, 1.08)
tau <- pow(sigma, -2)
sigma ~ dunif(0, 100)
}
'
我已经将aandb项的先验（假设一个线性模型：）设置为研究 1 的参数估计值。先验的标准偏差设置为 OLS 派生的标准误差。如果您运行分析（问题末尾的代码），您会发现这更好地反映了我们从两个数据集中知道的情况，因为 95% 的可信区间不包括零，中值估计接近真实效果。我会注意到，使用平坦的先验会导致较低的估计效果和更宽的置信区间。
话虽如此，我在已发表的研究中从未见过这样的事情，尽管我承认我很少在我的领域看到任何贝叶斯分析。据我了解，这类似于经验贝叶斯，但在这种情况下，我没有使用相同的数据来估计先验和后验。令我震惊的是，也许我所做的并没有表达出对先验的足够不确定性，但我不确定。
我知道许多人称这些主观先验是有原因的，但我认为这个问题本身并不是那么主观，即使可能有多种方法可以解决这个问题并模拟相同的一般情况，也无法在交叉验证上回答方法。
在 JAGS 中拟合上述模型的代码：
fit.inform <- run.jags(model.inform, c('a','b'),
data = list('donuts2' = donuts2,
'weight2' = weight2,
'N' = 100),
inits=list(list(.RNG.name="lecuyer::RngStream"),
list(.RNG.name="lecuyer::RngStream"),
list(.RNG.name="lecuyer::RngStream")),
burnin = 1000, sample = 5000, thin = 15, summarise = TRUE,
n.chains = 3)