贝叶斯分析Article
---
title: "贝叶斯分析Article"
author: "xuefliang"
date: "7/11/2019"
output:
html_document:
theme: readable
highlight: haddock
df_print: paged
code_folding: show
toc: true
toc_float: TRUE
number_sections: false
fig_width: 10.08
fig_height: 6
tags: [r, bayesian, posterior, test]
editor_options:
chunk_output_type: console
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE,warning = FALSE,message = FALSE,comment=">")
options(knitr.kable.NA = '',digits=2)
set.seed(333)
library(tidyverse)
library(bayestestR)
library(rstanarm)
library(insight)
```
## 可信区间（CI）
### 什么是可信区间？
可信区间是贝叶斯统计中的一个重要概念。 其核心目的是描述和总结与您的参数相关的不确定性 。 在这方面，它可能看起来与频率主义的置信区间很相似。然而，虽然他们的目标相似，但他们的统计定义和意义却截然不同。实际上，虽然后者是通过一个充满罕见测试假设和近似的复杂算法获得的，但可信区间的计算相当简单。
由于贝叶斯推断返回了可能的效应值（后验）的分布，可信区间只是包含特定百分比的可能值的范围。 例如，95％可信区间仅仅是后验分布的中心部分，其包含95％的值。
请注意，与频率主义区间相比，这大大提高了贝叶斯区间的可解释性。事实上，贝叶斯框架允许说“给定观察到的数据，效果有95％的概率落在这个范围内” ，而频率较低的直接替代方案（95％ 置信区间 ）将“ 有95％的可能性当从这种数据计算置信区间时，效果落在这个范围内 “。
### 为什么默认为89％？
当然，当选择CI级别进行默认报告时，人们开始使用95％，这是在频率论世界中使用的任意约定。 然而，一些作者认为95％可能不是贝叶斯后验分布最合适的，如果没有足够的后验样本可能缺乏稳定性（Kruschke 2014） 。该命题是使用90％而不是95％。 然而，最近，McElreath（2014年，2018年）建议，如果我们首先使用任意阈值，为什么不使用89％，因为这个值具有作为素数的附加参数。
因此，默认情况下，CI以89％的间隔（ci=0.89）计算，被认为比例如95％的间隔更稳定（Kruschke 2014）。如果应计算95％的间隔，则建议有效样本量至少为10.000（Kruschke，2014，p.183ff）。 此外，89是最高素数，不超过已经不稳定的95％阈值。什么与任何事情有关？ 没什么 ，但它提醒我们任何这些约定的完全任意性（McElreath 2014，@ mcelreath2018statistical） 。
### 不同类型的CI
读者可能会注意到bayestestR提供了两种计算可信区间的方法，即最高密度区间（HDI）（hdi()）和等尾区间（ETI） （eti()）。也可以通过ci()函数的method参数更改这些方法。有什么不同？ 让我们来看看：
```{r}
# Generate a normal distribution
posterior <- distribution_normal(1000)
# Compute HDI and ETI
ci_hdi <- ci(posterior, method = "HDI")
ci_eti <- ci(posterior, method = "ETI")
# Plot the distribution and add the limits of the two CIs
posterior %>%
estimate_density(extend=TRUE) %>%
ggplot(aes(x=x, y=y)) +
geom_area(fill="orange") +
theme_classic() +
# HDI in blue
geom_vline(xintercept=ci_hdi$CI_low, color="royalblue", size=3) +
geom_vline(xintercept=ci_hdi$CI_high, color="royalblue", size=3) +
# Quantile in red
geom_vline(xintercept=ci_eti$CI_low, color="red", size=1) +
geom_vline(xintercept=ci_eti$CI_high, color="red", size=1)
```
这些几乎完全一样......
但其他类型的分布也是如此吗？
```{r}
# Generate a beta distribution
posterior <- distribution_beta(1000, 6, 2)
# Compute HDI and Quantile CI
ci_hdi <- ci(posterior, method = "HDI")
ci_eti <- ci(posterior, method = "ETI")
# Plot the distribution and add the limits of the two CIs
posterior %>%
estimate_density(extend=TRUE) %>%
ggplot(aes(x=x, y=y)) +
geom_area(fill="orange") +
theme_classic() +
# HDI in blue
geom_vline(xintercept=ci_hdi$CI_low, color="royalblue", size=3) +
geom_vline(xintercept=ci_hdi$CI_high, color="royalblue", size=3) +
# Quantile in red
geom_vline(xintercept=ci_eti$CI_low, color="red", size=1) +
geom_vline(xintercept=ci_eti$CI_high, color="red", size=1)
```
这个差异很大。
与HDI相反，对于该HDI ，区间内的所有点具有比区间外的点更高的概率密度， ETI是等尾的 。 这意味着90％的间隔在其限制的任一侧具有5％的分布。 它表示第5百分位数和95h百分位数。 在对称分布中，计算可靠区间的两种方法（ETI和HDI）返回类似的结果。
对于偏斜的分布，情况并非如此。实际上，ETI中的参数值可能比ETI之外的参数值具有更低的可信度（更不可能）。 作为分布中可靠值的摘要，此属性似乎不合需要。
另一方面，当将变换应用于分布时，ETI范围确实会发生变化（例如，对于概率的对数几率缩放）：变换分布的下限和上限将对应于变换的下限和上限。相反，将变换应用于分布将改变最终的HDI。
## 方向概率（pd）
### 什么是pd？
方向概率（pd）是效应存在的指数，范围从50％到100％，表示效果在特定方向上的确定性（即 正或负）。除了解释，理解和计算的简单性之外，该索引还提供了其他有趣的属性：
它独立于模型 ：它完全基于后验分布，不需要来自数据或模型的任何附加信息。它对响应变量和预测变量的规模都很稳健。它与频率论p值密切相关，因此可用于绘制相似之处，并为不熟悉贝叶斯统计的读者提供参考。然而，该指数与评估影响的程度和重要性无关，通过其他指数（如ROPE百分比）可以更好地实现。事实上，重要性和存在的指数是完全独立的。您可以使用99.99％的pd效果，其中整个后验分布集中在[0.0001,0.002]范围内。在这种情况下，效果是肯定的，具有高确定性 ，但也不显著 （即非常小）。
有效指数存在，例如pd ，在探索性研究或临床研究中特别有用，其重点是确保感兴趣的效果不是相反的方向（对于临床研究，如治疗无害）。然而，一旦效果的方向得到确认，焦点就应转向其重要性，包括对其大小，相关性和重要性的精确估计。
### 与p值的关系
在大多数情况下，似乎pd通过以下公式与频率为单侧的p值直接对应：$\ [p_ {one sided} = 1-p_d / 100 ]$同样， 双边p值 （最常报道的一个）通过以下公式是等价的： $\ [p_ {双面} = 2 *（1-p_d / 100）]$因此，双边p值分别为.1 ， .05 , . 01和.001大致对应于95％ ， 97.5％ ， 99.5％和99.95％的pd 。 但如果它像p值一样，那一定是坏的，因为p值很差[ 参见重复性危机]。
事实上，重复性危机这一方面可能被误解了。实际上，并不是p值是一个本质上错误或错误的指数。 相反，它的误用 ，误解和误解会助长局势的衰退。例如，pd与p值高度相关的事实表明后者更多地是效应存在的指数而非重要性（即，“感兴趣的价值”）。贝叶斯版本pd具有直观的含义，并且明显表明所有阈值都是任意的。然而，pd的数学和解释透明度及其作为效应存在指数的重新概念化，为贝叶斯结果的表征提供了有价值的见解。此外，它与频率主义p值的接近使其成为一个完美的衡量标准，可以简化心理学研究向贝叶斯框架的采用。
### 计算方法
计算pd的最简单和直接的方法是1）查看中位数的符号，2）选择相同符号的后部分，3）计算该部分代表的百分比。 这种“简单”的方法是最困难的方法，但其精确度与后验抽取的数量直接相关。第二种方法依赖于密度估计。它首先估算密度函数（有许多方法可用），然后计算0的另一侧密度曲线的曲线下面积（AUC）。基于密度的方法可以假设更精确，但很大程度上取决于用于估算密度函数的方法。
### 方法比较
让我们比较4种可用方法，直接方法和3种基于密度的方法，它们的密度估计算法不同（见estimate_density ）。
### 关联
让我们首先测试通过不同方法获得的结果的接近度和相似性。
```{r}
library(bayestestR)
library(logspline)
library(KernSmooth)
# Compute the correlations
data <- data.frame()
for(the_mean in runif(25, 0, 4)){
for(the_sd in runif(25, 0.5, 4)){
x <- rnorm(100, the_mean, abs(the_sd))
data <- rbind(data,
data.frame("direct" = pd(x),
"kernel" = pd(x, method="kernel"),
"logspline" = pd(x, method="logspline"),
"KernSmooth" = pd(x, method="KernSmooth")
))
}
}
data <- as.data.frame(sapply(data, as.numeric))
# Visualize the correlations
library(GGally)
GGally::ggpairs(data) +
theme_classic()
```
所有给出的方法都是高度相关的，并给出非常相似的结果。这意味着方法选择并不激烈，不能用于过多地调整结果。
### 准确性
为了测试每种方法的准确性，我们将从非常密集的分布（具有大量观察值）计算直接pd开始。 这将是我们的基线，或“真实”pd。然后，我们将迭代地从该父分布中绘制较小的样本，并且我们将使用不同的方法计算pd 。 这个估计值与参考值越接近越好。
```{r}
data <- data.frame()
for(i in 1:25){
the_mean <- runif(1, 0, 4)
the_sd <- abs(runif(1, 0.5, 4))
parent_distribution <- rnorm(100000, the_mean, the_sd)
true_pd <- pd(parent_distribution)
for(j in 1:25){
sample_size <- round(runif(1, 25, 5000))
subsample <- sample(parent_distribution, sample_size)
data <- rbind(data,
data.frame("sample_size" = sample_size,
"true" = true_pd,
"direct" = pd(subsample) - true_pd,
"kernel" = pd(subsample, method="kernel")- true_pd,
"logspline" = pd(subsample, method="logspline") - true_pd,
"KernSmooth" = pd(subsample, method="KernSmooth") - true_pd
))
}
}
data <- as.data.frame(sapply(data, as.numeric))
data %>%
tidyr::gather(Method, Distance, -sample_size, -true) %>%
ggplot(aes(x=sample_size, y = Distance, color = Method, fill= Method)) +
geom_point(alpha=0.3, stroke=0, shape=16) +
geom_smooth(alpha=0.2) +
geom_hline(yintercept=0) +
theme_classic() +
xlab("\nDistribution Size")
```
基于“核心”的密度方法似乎一直低估了pd 。有趣的是，即使在少量后抽取的情况下，“直接”方法看起来也更可靠。
## 实用对等区（ROPE）
### 什么是ROPE？
与频率论方法不同，贝叶斯推断不是基于统计显著性，其中效果针对“零”进行测试。实际上，贝叶斯框架提供了参数的概率视图，允许评估与它们相关的不确定性。因此，我们可以得出结论，在超出特定范围之外的可能被视为“实际上没有效果”（即，可忽略不计的幅度）的概率就足够了，而不是断定效果存在。 该范围称为实际等效区域（ROPE） 。
实际上，从统计上来说，后验分布与0不同的概率没有多大意义（它与单个点的无限概率不同）。 因此，强调ROPE的想法是让用户定义空值周围的区域，其中包含等于实际目的的空值的值（Kruschke 2010,2014; Kruschke，Aguinis和Joo 2012） 。
### 等效性测试
ROPE是对应于“无效”假设的区域，用于等效性测试，以测试参数是否显著（在重要程度上足以被关注）。 该测试通常基于“HDI + ROPE决策规则” 来检查是否应该针对明确制定的“零假设”（ 即 ，ROPE）接受或拒绝参数值。换句话说，它检查作为空区域（ROPE）的可信区间（CI）的百分比。 如果该百分比足够低，则拒绝原假设。如果该百分比足够高，则接受零假设。ROPE中的
### 可信区间与ROPE中的完全后验区间使用ROPE和HDI作为可信区间
Kruschke（2018）建议使用落在ROPE内的95％HDI的百分比作为决策规则。然而，由于89％的HDI被认为是更好的选择,因此bayestestR默认提供了落在ROPE内的89％HDI的百分比。然而，初步数据倾向于表明使用完全后验分布的百分比而不是CI可能更敏感（特别是描绘非常显著的影响）。因此，我们建议用户考虑使用完整的ROPE百分比（通过设置ci=1），这将返回ROPE中整个后验分布的部分。
### ROPE接受或拒绝的百分比是多少？
如果HDI完全在ROPE之外，则该参数的“零假设”被“拒绝”。如果ROPE完全涵盖HDI，即参数的所有最可靠值都在实际等价区域内，则接受零假设。 否则，尚未决定是接受还是拒绝零假设。
如果使用完整的ROPE（即100％的HDI），则如果ROPE中后验的百分比小于2.5％或大于97.5％，则拒绝或接受零假设。 理想的结果是ROPE内部的比例较低（越接近于零越好）。
### 如何定义ROPE范围？
Kruschke（2018）提出，默认情况下，ROPE可以设置为标准化参数的-0.1到0.1的范围（根据Cohen1988可忽略的影响大小）。
* 对于线性模型（lm） ，这可以推广为：$\ [ - 0.1 * SD_ {y}，0.1 * SD_ {y}] $ 。
* 对于逻辑模型 ，以对数优势比表示的参数可通过以下公式转换为标准差： $\ [\ sqrt {3} / \ \pi ] $，从而产生-0.055到-0.055的范围。对于具有二进制结果的其他模型，强烈建议手动指定rope参数。 目前，对于逻辑模型应用相同的默认值。
* 对于t检验 ，使用响应的标准偏差，类似于线性模型（见上文）。
* 对于相关性 ，使用-0.05, 0.05 即 Cohen（1988）经验法则所建议的可忽略相关值的一半。
* 对于所有其他模型， -0.1, 0.1用于确定ROPE限制，但强烈建议手动指定。
### 对参数刻度的灵敏度
当使用基于ROPE的指数时，考虑预测变量的单位（即尺度）是重要的，因为ROPE的正确解释表示实际等效于零的区域取决于预测变量的尺度。实际上，与其他指数（例如pd）不同，ROPE中的百分比取决于其参数的单位。换句话说，由于ROPE表示响应比例的固定部分，因此其与系数的接近程度取决于系数本身的比例。
例如，如果我们考虑一个简单的回归growth ~ time，模拟Wookies婴儿的发展，可忽略的变化（ROPE）小于54厘米。如果我们的time变量用天表示 ，我们会发现系数（代表白天的增长）大约是10厘米（系数后验的中位数是10）。我们认为这可以忽略不计。但是，如果我们决定以年为单位表示time变量，系数将通过此变换进行缩放（因为它现在将代表按年增长）。系数现在约为3550厘米（10*355），我们现在认为这是重要的 。
```{r}
library(see)
data <- iris # Use the iris data
model <- stan_glm(Sepal.Length ~ Sepal.Width, data=data) # Fit model
# Compute indices
pd <- p_direction(model)
percentage_in_rope <- rope(model, ci=1)
# Visualise the pd
plot(pd)
pd
# Visualise the percentage in ROPE
plot(percentage_in_rope)
percentage_in_rope
```
我们可以看出， Sepal.Length和Sepal.Width之间线性关系的pd和ROPE百分比分别约为92.95%和15.95% ，对应于不确定和不显著的影响。 如果我们扩展预测器会发生什么？
```{r}
data$Sepal.Width_scaled <- data$Sepal.Width / 100 # Divide predictor by 100
model <- stan_glm(Sepal.Length ~ Sepal.Width_scaled, data=data)
# Compute indices
pd <- p_direction(model)
percentage_in_rope <- rope(model, ci=1)
# Visualise the pd
plot(pd)
pd
# Visualise the percentage in ROPE
plot(percentage_in_rope)
percentage_in_rope
```
正如您所看到的，通过简单地将预测器除以100，我们彻底改变了与ROPE中的百分比相关的结论（其变得非常接近0）：现在可以将效果解释为重要。因此，我们建议在选择ROPE范围时密切关注预测变量的单位（例如，哪个系数对应于一个小影响？），以及报告或读取ROPE结果时。
### 多重共线性：非独立协变量
当参数显示强相关性时，即，当协变量不是独立的时，联合参数分布可以朝向或远离ROPE移位。 共线性使基于单变量边缘的ROPE和假设检验无效，因为概率是以独立性为条件的。最有问题的是仅与ROPE区域部分重叠的参数。在共线性的情况下，这些参数的（联合）分布可以获得增加或减少的ROPE，这意味着基于ROPE的推断是不合适的（Kruschke 2014）。equivalence_test()和rope()函数对参数之间的成对相关性执行简单检查，但由于两个以上变量之间可能存在共线性，因此检查此假设检验假设的第一步是查看不同的对图。 更复杂的检查是投影预测变量选择（Piironen和Vehtari 2017） 。
## 贝叶斯因子
贝叶斯应用统计框架的采用，特别是在社会或心理科学领域，似乎正朝着两个不同的方向发展。 标记他们分离的关键主题之一是他们对贝叶斯因子的看法。简而言之，一些作者（例如，由Wagenmakers领导的“阿姆斯特丹学校”）提倡使用并强调其作为统计指标的质量，而其他作者则指出其局限性，而更倾向于对后验分布的精确描述（使用CI）， ROPEs等）。
bayestestR在本次辩论中没有参与，而是提供工具来帮助您进行任何想要实现的分析。相反，它强烈支持一个明智的选择的概念： 发现方法，尝试它们，理解它们，了解它们，并自己决定 。
话虽如此，这里是贝叶斯因子的介绍:
### 贝叶斯因子
贝叶斯因子（BF）是一个“模型”相对于另一个“模型”的相对证据的指数，其在贝叶斯推断中用作经典（频率论）假设检验指数（例如p-值）的替代。在贝叶斯框架中，贝叶斯因子也可以被认为是根据观察到的数据更新关于两个模型的相对几率的一些先验信念的量。
根据贝叶斯定理：
$$[P（M | D）= \frac {P（D | M）P（M）} {P（D）} ]$$
然后通过比较两个模型，我们得到：
$$ [\frac {P（M_1 | D）} {P（M_2 | D）} = \frac {P（D | M_1）} {P（D | M_2）} \times \frac {P（M_1）} { P（M_2）} ]$$
其中中间项是贝叶斯因子：
$$[BF_ {12} = \frac {P（D | M_1）} {P（D | M_2）}]$$ 这些通常按比例计算两个竞争假设/模型的边际可能性，但正如我们从上面的等式中可以看到的，它们也可以通过将后验概率除以先验概率来计算。重要的是，贝叶斯因子涵盖了广泛的指数和应用 ，并有不同的风格。
### 比较模型
一般来说，贝叶斯因子用于比较模型并回答问题：在哪种模型下观察到的数据更可能？换句话说，哪种模型更有可能产生观察到的数据？这通常通过计算两个模型的边际可能性来完成。在这种情况下，贝叶斯因子是一个比较模型相对于另一个模型的相对证据的度量。
### 贝叶斯模型（ brms和rstanarm ）
注意：为了计算模型的贝叶斯因子，必须在拟合时添加非默认参数：
* brmsfit模型必须符合save_all_pars = TRUE
* stanreg模型必须已经安装了定义的diagnostic_file 。
让我们首先用5个贝叶斯回归与brms来预测Sepal.Length ：
```{r}
library(brms)
m0 <- brm(Sepal.Length ~ 1, data = iris, save_all_pars = TRUE)
m1 <- brm(Sepal.Length ~ Petal.Length, data = iris, save_all_pars = TRUE)
m2 <- brm(Sepal.Length ~ Species, data = iris, save_all_pars = TRUE)
m3 <- brm(Sepal.Length ~ Species + Petal.Length, data = iris, save_all_pars = TRUE)
m4 <- brm(Sepal.Length ~ Species * Petal.Length, data = iris, save_all_pars = TRUE)
```
我们现在可以将这些模型与bayesfactor_models()函数进行比较，使用denominator参数指定将与所有模型进行比较的模型（在本例中为仅截距模型）：
```{r}
comparison <- bayesfactor_models(m1, m2, m3, m4, denominator = m0)
comparison
```
我们可以看到完整模型是最好的模型 - 与$（BF _ {\text {m0}} = 9 \times 10 ^ {55} ）$与null（仅截距）相比。我们还可以将参考模型更改为主效应模型：
```{r}
update(comparison, reference = 3)
```
我们可以看到，尽管完整模型是最好的，但几乎没有证据表明它比主效应模型更可取。我们还可以将参考模型更改为Species模型：
```{r}
update(comparison, reference = 2)
```
请注意，在贝叶斯框架中，比较模型不需要是嵌套模型，就像我们将Petal.Length模型与Species Petal.Length模型进行比较时所发生的那样（在频率论者框架中无法做到的事情，比较模型必须互相嵌套）。
注意：为了正确和精确地估计贝叶斯因子，你总是需要4P: Proper Priors (1, 2, 3), and a Plentiful Posterior (4)。
### 频率论模型的BIC近似
也可以为频率模型计算贝叶斯因子。这是通过比较BIC测量来完成的，允许对非嵌套频率模型进行贝叶斯比较（Wagenmakers 2007） 。 让我们试试一些线性混合模型 ：
```{r}
library(lme4)
m0 <- lmer(Sepal.Length ~ (1 | Species), data = iris)
m1 <- lmer(Sepal.Length ~ Petal.Length + (1 | Species), data = iris)
m2 <- lmer(Sepal.Length ~ Petal.Length + (Petal.Length | Species), data = iris)
m3 <- lmer(Sepal.Length ~ Petal.Length + Petal.Width + (Petal.Length | Species), data = iris)
m4 <- lmer(Sepal.Length ~ Petal.Length * Petal.Width + (Petal.Length | Species), data = iris)
bayesfactor_models(m1, m2, m3, m4, denominator = m0)
```
### 通过贝叶斯模型平均包含贝叶斯因子
包含贝叶斯因子回答了这个问题：在具有特定效果的模型下观察到的数据是否比在没有特定效果的模型下更可能？ 换句话说，平均而言-具有效果的模型$X$更有可能产生观察到的数据而不是没有效果的模型$X$？让我们使用上面的brms示例：
```{r}
bayesfactor_inclusion(comparison)
```
如果我们检查交互项的包含贝叶斯因子，我们可以看到，在所有5个模型中，具有交互项（Species:Petal.Length）的模型比没有交互项的模型的可能性高5倍。
我们也可以仅比较匹配的模型-即，没有效果的模型$A$将仅与具有效果$A$的模型进行比较，但不与具有更高级别交互的模型进行比较（请参阅解释为什么您可能需要在这里这样做）。
```{r}
bayesfactor_inclusion (comparison, match_models = TRUE )
```
#### 与JASP比较
bayesfactor_inclusion()旨在提供跨模型平均值的贝叶斯因子，类似于JASP的效果选项。 让我们比较两者：
比较所有模型
```{r}
library(BayesFactor)
ToothGrowth$dose <- as.factor(ToothGrowth$dose)
BF_ToothGrowth <- anovaBF(len ~ dose*supp, ToothGrowth)
bayesfactor_inclusion(BF_ToothGrowth)
```
比较匹配的模型
```{r}
bayesfactor_inclusion(BF_ToothGrowth, match_models = TRUE)
```
使用Nuisance Effects
我们将在JASP中为null模型添加剂量，并在R执行相同的操作：
```{r}
BF_ToothGrowth_against_dose <- BF_ToothGrowth[3:4]/BF_ToothGrowth[2] # OR:
# update(bayesfactor_models(BF_ToothGrowth), subset = c(4,5), reference = 3)
BF_ToothGrowth_against_dose
bayesfactor_inclusion(BF_ToothGrowth_against_dose)
```
### Savage-Dickey密度比率贝叶斯因子
Savage-Dickey密度比可用于回答以下问题：鉴于观察到的数据，null更多还是更不可能？
这是通过比较先前和后验分布之间的空值密度来完成的，并且是比较提供的（点）空模型的贝叶斯因子的近似值：
$H_0$ 与$H_1$的贝叶斯因子可以通过分析积分模型参数$\theta$获得 。 然而，贝叶斯因子同样可以通过仅考虑$H_1$，并且将后部的高度除以$\theta$的高度除以先前的高度$\theta$ ，兴趣。“ （Wagenmakers et al.2010）
让我们使用学生睡眠数据，并尝试回答这个问题： 根据观察到的数据，药物（变量group ）或多或少可能对额外睡眠时数（可变extra）没有影响？
blopplot表明第二组的额外睡眠时间更长 。 多少钱？ 让我们拟合一个简单的贝叶斯线性模型 。
```{r}
library(rstanarm)
model <- stan_glm(extra ~ group, data = sleep)
```
我们可以在此模型上使用as.data.frame()来提取与group2的效果相关的后验分布，以及来自insight包的 get_priors() ，以查看使用的先前分布：
```{r}
posterior <- as.data.frame(model)$group2
insight::get_priors(model)
```
对于group2参数，使用的先验是平均值（位置）0和SD（比例）5.044799。我们可以模拟这个先前的分布如下：
```{r}
prior <- distribution_normal(length(posterior), mean = 0, sd = 5.044799)
```
我们现在可以绘制先前和后验分布：观察分布，我们可以看到后验中心位于1.56。 但这并不意味着0的影响必然不太可能。 为了测试它，我们将使用bayesfactor_savagedickey() 。
```{r}
test_group2 <- bayesfactor_savagedickey(posterior = posterior, prior = prior)
test_group2
```
该BF表示在给定数据的情况下，0（点空效应模型）的效果不太可能为0.88 。 或者，当解释为$H_0$与$H_1$的贝叶斯因子时，我们可以说观察到的数据在空值下比在具有效果的模型下更可能是1 / 0.88=1.14倍！因此，尽管分布中心偏离了零点，但它仍然在零点周围非常密集。请注意，可以在此处找到贝叶斯因子的解释指南 。
#### 定向测试
如果我们先前有关于效果方向的假设，我们也可以进行方向测试（“单侧”或“单尾”测试）。这是通过对先前和后验分布设置订单限制来完成的（Morey和Wagenmakers 2014; Morey 2015） 。
```{r}
test_group2_right <- bayesfactor_savagedickey(posterior = posterior, prior = prior, direction = ">")
test_group2_right
```
#### 测试所有模型参数
或者，我们也可以将模型直接传递给bayesfactor_savagedickey()以同时测试所有模型的参数：
```{r}
bayesfactor_savagedickey(model)
```
#### 测试对比（与emmeans ）
我们也可以使用bayesfactor_savagedickey()和emmeans ，让我们测试贝叶斯对比 。
```{r}
library(emmeans)
group_diff <- pairs(emmeans(model, ~ group))
group_diff
bayesfactor_savagedickey(group_diff, prior = model)
```
## 深度1：点估计的比较
### 贝叶斯框架中的效应点估计
### 介绍
贝叶斯框架和频率框架之间的主要区别之一是前者返回每个效应的概率分布（即模型的感兴趣参数，例如回归斜率）而不是单个值。然而，对于报告或用于进一步分析的需求和需求，仍然需要能够最好地表征潜在后验分布的单个值（点估计）。在文献中使用三个主要指数用于效果估计：平均值，中值或MAP（最大后验）估计（大致对应于分布的模式“峰值”）。不幸的是，没有就使用哪一个达成共识，因为从未进行过系统的比较。在目前的工作中，我们将比较这三个点之间的效果估计值，以及从可比较的频率论模型中提取的广为人知的beta 。 通过这种比较，我们期望在两个框架之间架起桥梁和关系，帮助和缓解公众的过渡。
### 实验1：与误差（噪声）和样本大小的关系
#### 方法
模拟旨在调整以下特征：
模型的类型 ：线性或逻辑。
“真实”效果 （绘制数据的原始回归系数）：可以是1或0（无效果）。
样本量 ：从20到100步长为10。
误差 ：应用于预测器的高斯噪声，SD均匀分布在0.33和6.66之间（具有1000个不同的值）。
我们为这些特征的每个组合生成了一个数据集，从而产生了总共2*2*9*1000=36000贝叶斯和频率模型。这里可以使用用于生成的代码（请注意，通常需要几天/几周才能完成）。
```{r}
library(see)
df <- read.csv("https://raw.github.com/easystats/circus/master/data/bayesSim_study1.csv")
```
#### 结果
对噪音的敏感度
```{r}
df %>%
select(error, true_effect, outcome_type, beta, Median, Mean, MAP) %>%
gather(estimate, value, -error, -true_effect, -outcome_type) %>%
mutate(temp = as.factor(cut(error, 10, labels = FALSE))) %>%
group_by(temp) %>%
mutate(error_group = round(mean(error), 1)) %>%
ungroup() %>%
filter(value < 6) %>%
ggplot(aes(x = error_group, y = value, fill = estimate, group = interaction(estimate, error_group))) +
# geom_hline(yintercept = 0) +
# geom_point(alpha=0.05, size=2, stroke = 0, shape=16) +
# geom_smooth(method="loess") +
geom_boxplot(outlier.shape=NA) +
theme_modern() +
scale_fill_manual(values = c("beta" = "#607D8B", "MAP" = "#795548", "Mean" = "#FF9800", "Median" = "#FFEB3B"),
name = "Index") +
ylab("Point-estimate") +
xlab("Noise") +
facet_wrap(~ outcome_type * true_effect, scales="free")
```
对样本量的敏感度
```{r}
df %>%
select(sample_size, true_effect, outcome_type, beta, Median, Mean, MAP) %>%
gather(estimate, value, -sample_size, -true_effect, -outcome_type) %>%
mutate(temp = as.factor(cut(sample_size, 10, labels = FALSE))) %>%
group_by(temp) %>%
mutate(size_group = round(mean(sample_size))) %>%
ungroup() %>%
filter(value < 6) %>%
ggplot(aes(x = size_group, y = value, fill = estimate, group = interaction(estimate, size_group))) +
# geom_hline(yintercept = 0) +
# geom_point(alpha=0.05, size=2, stroke = 0, shape=16) +
# geom_smooth(method="loess") +
geom_boxplot(outlier.shape=NA) +
theme_modern() +
scale_fill_manual(values = c("beta" = "#607D8B", "MAP" = "#795548", "Mean" = "#FF9800", "Median" = "#FFEB3B"),
name = "Index") +
ylab("Point-estimate") +
xlab("Sample size") +
facet_wrap(~ outcome_type * true_effect, scales="free")
```
统计建模
我们拟合（频率论）多元线性回归来统计测试预测效果的存在与否，以及它们与噪声和样本大小的相互作用。
```{r}
df %>%
select(sample_size, error, true_effect, outcome_type, beta, Median, Mean, MAP) %>%
gather(estimate, value, -sample_size, -error, -true_effect, -outcome_type) %>%
glm(true_effect ~ outcome_type / estimate / value, data=., family="binomial") %>%
broom::tidy() %>%
select(term, estimate, p=p.value) %>%
filter(stringr::str_detect(term, 'outcome_type'),
stringr::str_detect(term, ':value')) %>%
arrange(desc(estimate)) %>%
knitr::kable(digits=2)
```
这表明，与频率论者的beta相比，为了描述效果的存在与否;
* 对于线性模型， 均值是更好的预测因子，紧随其后的是中位数 ， MAP和频率主义beta 。
* 对于逻辑模型， MAP是更好的预测因子，其次是中位数 ， 均值和频率测试 。
总的来说，中位数似乎是一个恰当的选择，在不同类型的模型中保持高性能。
### 实验2：与采样特征的关系
#### 方法
模拟旨在调整以下特征：
* 模型类型 ：线性或逻辑。
* “真实”效果 （绘制数据的原始回归系数）：可以是1或0（无效果）。
* 绘制 ：从5到1000（步长为5）（1000次迭代）。
* 预热 ：预热 迭代的比率。 从1/10到9/10，步长为0.1（9次迭代）。
我们为这些特征的每个组合生成了3个数据集，从而得到总共2*2*8*40*9*3=34560贝叶斯和频率模型。这里可以使用用于生成的代码（请注意，通常需要几天/几周才能完成）。
```{r}
df <- read.csv("https://raw.github.com/easystats/circus/master/data/bayesSim_study2.csv")
```
#### 结果
对迭代次数的敏感性
```{r}
df %>%
select(iterations, true_effect, outcome_type, beta, Median, Mean, MAP) %>%
gather(estimate, value, -iterations, -true_effect, -outcome_type) %>%
mutate(temp = as.factor(cut(iterations, 5, labels = FALSE))) %>%
group_by(temp) %>%
mutate(iterations_group = round(mean(iterations), 1)) %>%
ungroup() %>%
filter(value < 6) %>%
ggplot(aes(x = iterations_group, y = value, fill = estimate, group = interaction(estimate, iterations_group))) +
geom_boxplot(outlier.shape=NA) +
theme_classic() +
scale_fill_manual(values = c("beta" = "#607D8B", "MAP" = "#795548", "Mean" = "#FF9800", "Median" = "#FFEB3B"),
name = "Index") +
ylab("Point-estimate of the true value 0\n") +
xlab("\nNumber of Iterations") +
facet_wrap(~ outcome_type * true_effect, scales="free")
```
对预热比率的敏感度
```{r}
df %>%
mutate(warmup = warmup / iterations) %>%
select(warmup, true_effect, outcome_type, beta, Median, Mean, MAP) %>%
gather(estimate, value, -warmup, -true_effect, -outcome_type) %>%
mutate(temp = as.factor(cut(warmup, 3, labels = FALSE))) %>%
group_by(temp) %>%
mutate(warmup_group = round(mean(warmup), 1)) %>%
ungroup() %>%
filter(value < 6) %>%
ggplot(aes(x = warmup_group, y = value, fill = estimate, group = interaction(estimate, warmup_group))) +
geom_boxplot(outlier.shape=NA) +
theme_classic() +
scale_fill_manual(values = c("beta" = "#607D8B", "MAP" = "#795548", "Mean" = "#FF9800", "Median" = "#FFEB3B"),
name = "Index") +
ylab("Point-estimate of the true value 0\n") +
xlab("\nNumber of Iterations") +
facet_wrap(~ outcome_type * true_effect, scales="free")
```
### 实验3：与Priors规范的关系
结论可以在指南部分找到。
## 深度2：效应存在指数和显著指数的比较
### 贝叶斯框架中效应存在的指标及其意义
现在普遍认为贝叶斯统计框架是进行心理科学研究的正确方法。然而，它的灵活性是它的力量和弱点，因为没有就应该计算或报告哪些指数达成一致。此外，缺乏一个共识的效应存在指数，例如频率论的p值，可能导致许多非熟悉的读者在贝叶斯统计中感知到的不必要的模糊。因此，本研究描述并比较了几种效应存在指数，提供了这些指数“行为”的直观视觉表示，与传统指标（如样本量和频率显著性）相关。结果有助于培养对研究人员报告的价值观的直观理解，并允许为贝叶斯统计描述提出建议，这对科学报告的标准化至关重要。
#### 介绍
贝叶斯框架在心理学家和神经科学家中迅速普及。更喜欢这种方法的原因是可靠性，噪声数据的准确度更高，小样本的估计更好，I类错误的倾向性更小，将先验知识引入分析的可能性，以及关键的结果直观性和直接的解释。频率论的方法与对零假设检验的关注有关，并且p值的误用已被证明对心理科学的再现性危机起着关键作用。人们普遍认为，贝叶斯方法的推广是一种克服这一问题的方法 。
虽然贝叶斯框架推动的概率推理遍及大多数数据科学方面，但它已经为统计建模奠定了基础。心理学大量依赖的这个方面可以大致分为两个软边类别;预测和结构建模。虽然统计模型可以（通常）同时满足这两个目的，但预测建模致力于构建并找到准确预测给定结果的最佳模型。 它以适合度量，预测准确性和模型比较等概念为中心。在这个维度的极端情况下，机器和深度学习模型被用于强大的预测能力，通常以牺牲人类可读性为代价（这些模型通常被称为“黑盒子”，强调评估其内部功能的难度;另一方面，心理学家经常使用更简单的模型（例如与一般线性框架相关）来探索他们的数据。在此框架内，目标从构建最佳模型转变为理解模型内的参数。实际上，方法论管道通常从涉及模型比较的预测建模开始“什么是世界上最好的模型（即观察到的变量）”）后看似过渡到结构建模：“给定这个世界模型，如何效果（即模型参数）正在影响结果“。 对于最后一部分，他们通常依赖于效果指数“存在”。
事实上，虽然贝叶斯框架的优势之一是其概率参数估计，允许量化与每个估计相关的固有不确定性，但心理学家也对参数存在感兴趣：一个决策标准，允许他们得出结论，如果效果是“不同的从0“（统计上对应于”不可忽略的“或”不是相反的方向“）。换句话说，在了解效果的重要性，相关性或强度之前，要知道它是否与给定方向的结果有关。这种需要导致了频率主义p值的广泛采用，用作效应存在的指标，并且它的接受伴随着为其分类创建任意聚类（.05，.01和.001）。不幸的是，这些启发式方法严重僵化，成为达到目标和门槛而不是理解数据的工具 。
因此，贝叶斯框架在不需要这种零假设检验指标的情况下回答心理问题的能力通常被宣传为“没有p值的新世界”的承诺，而不是旧的和有缺陷的频率问题。尽管如此，似乎“效应存在”指数和标准对于人类获得对其数据的相互作用和结构的直观理解是有用的。因此，贝叶斯用户友好型实施的发展伴随着贝叶斯因子（BF）的推广，这一指数反映了模型对另一个模型的预测性能 。它提供了超过p值的许多优点，具有直接的解释（“数据是在替代下比零假设更容易发生3次（ BF = 3）”并且允许对替代方案做出陈述，而不仅仅是零假设。此外，最近的数学发展允许其计算复杂模型。 尽管与p值相比，BF符合预期的稳定，有效，直观和更好的指数，但其用于模型选择仍然是一个争论的问题。实际上，由于用于其计算的预测是从模型参数的先前分布生成的，因此它高度依赖于先验规范。重要的是，对于本文的目的，其用于估计较大模型中参数的效果存在的用途仍然是有限的，其计算在技术上是困难的并且其解释不像诸如t检验或相关性的“简单”测试那样简单。然而，考虑到对这些信息的需求，研究人员基于后验分布的特征开发了其他指数，其表示给定观察数据的不同参数值的概率分布。例如，通过呈现中心点（平均值，中位数，......）和离散度（标准偏差，中位数绝对偏差，......）的点估计值，可以总结这种不确定性，通常伴随一个百分比（89％，90％或95）最高密度区间的百分比（HDI;称为可信区间 - CI）。尽管贝叶斯框架提供了计算许多效应存在指数的可能性，但尚未对使用的存在指标达成共识，因为从未进行过任何比较。对于有兴趣采用贝叶斯框架的科学家来说，这可能是一个反驳。 此外，这个灰色区域会增加不熟悉贝叶斯框架的读者或评论者遵循假设和结论的难度，这可能反过来对整个研究产生不必要的怀疑。虽然我们认为这些指数及其解释指南（以经验法则的形式）在实践中是有用的，但我们也坚信这些指数应该伴随着与样本大小和效应大小相关的“行为”的知识。 。 这些知识对于人们隐含和直观地评估他们报告的数学价值的含义和含义非常重要。反过来，这可以防止可能的启发式和从这些指数得到的类别的结晶。
因此，基于多元线性回归（最广泛使用的模型之一）的模拟，本工作旨在比较仅从后验分布得出的几个效应存在指数，提供这些指数的“行为”的视觉表示。与样本大小，噪声，先验的关系以及频率论p值（超出其众多缺陷的指数，众所周知，可用作贝叶斯新手的参考），并为贝叶斯统计报告提出建议。
对于所有模拟模型，我们计算了以下指数：
* p值 ：基于频率回归，该指数表示给定统计模型的概率，当零假设为真时，统计汇总（例如两个比较组之间的样本均值差异）将大于或等于实际观察结果（Wasserstein，Lazar等人，2016年） 。
* “ 方向概率”（pd）：对应于效果严重为正或负的概率（以百分比表示）（与中位数符号一致）。
* p -MAP ：基于MAP的p值与参数对零假设的几率相关 。
* p -ROPE ：基于ROPE的p值，表示在ROPE定义的可忽略值空间中不包含（正值）或完全包含（负值）的HDI的最大百分比。
* ROPE（89％） ：该指数指的是ROPE中89％HDI的百分比，该区间被认为比基于频率的95％更稳定 。
* ROPE（95％） ：该指数是指ROPE中95％HDI的百分比。该间隔对应于Kruschke（2014）最初建议的间隔 。
* ROPE（完整） ：该指数是指ROPE内整个后验分布的百分比。
* 贝叶斯因子 ：TODO。
实用等效区域（ROPE）定义为-0.1至0.1。
### 研究1：与噪音和样本量的关系
#### 方法
模拟旨在调整以下特征：
* 模型类型 ：线性或逻辑。
* “真实”效果 （绘制数据的原始回归系数）：可以是1或0（无效果）。
* 样本量 ：从20到100步长为10。
* 误差 ：应用于预测器的高斯噪声，SD均匀分布在0.33和6.66之间（具有1000个不同的值）。
我们为这些特征的每个组合生成了一个数据集，从而产生了总共2*2*9*1000=36000贝叶斯和频率模型。 这里可以使用用于生成的代码（请注意，通常需要几天/几周才能完成）。
```{r}
df <- read.csv("https://raw.github.com/easystats/circus/master/data/bayesSim_study1.csv")
```
#### 结果
效果检测
对噪音的敏感度
```{r}
df %>%
select(outcome_type, true_effect, error, sample_size, p_value, p_direction, p_MAP, p_ROPE, ROPE_89, ROPE_95, ROPE_full, BF_log) %>%
gather(index, value, -error, -sample_size, -true_effect, -outcome_type) %>%
mutate(true_effect = as.factor(true_effect),
index = factor(index, levels=c("p_value", "p_direction", "p_MAP", "p_ROPE", "ROPE_89", "ROPE_95", "ROPE_full", "BF_log"))) %>%
mutate(temp = as.factor(cut(error, 10, labels = FALSE))) %>%
group_by(temp) %>%
mutate(error_group = round(mean(error), 1)) %>%
ungroup() %>%
ggplot(aes(x = error_group, y = value, fill = index, colour=true_effect, group = interaction(index, true_effect, error_group))) +
# geom_jitter(shape=16, alpha=0.02) +
geom_boxplot(outlier.shape = NA) +
facet_wrap(~index * outcome_type, scales = "free", ncol=2) +
theme_modern() +
scale_color_manual(values = c(`0` = "#f44336", `1` = "#8BC34A", name="Effect")) +
scale_fill_manual(values = c("p_value"="#607D8B", "p_MAP" = "#4CAF50", "p_direction" = "#2196F3",
"ROPE_89" = "#FFC107", "ROPE_95" = "#FF9800", "ROPE_full" = "#FF5722",
"p_ROPE"="#E91E63", "BF_log"="#9C27B0"), guide=FALSE) +
ylab("Index Value") +
xlab("Noise")
```
对样本量的敏感度
```{r}
df %>%
select(outcome_type, true_effect, error, sample_size, p_value, p_direction, p_MAP, p_ROPE, ROPE_89, ROPE_95, ROPE_full, BF_log) %>%
gather(index, value, -error, -sample_size, -true_effect, -outcome_type) %>%
mutate(true_effect = as.factor(true_effect),
index = factor(index, levels=c("p_value", "p_direction", "p_MAP", "p_ROPE", "ROPE_89", "ROPE_95", "ROPE_full", "BF_log"))) %>%
mutate(temp = as.factor(cut(sample_size, 10, labels = FALSE))) %>%
group_by(temp) %>%
mutate(size_group = round(mean(sample_size))) %>%
ungroup() %>%
ggplot(aes(x = size_group, y = value, fill = index, colour=true_effect, group = interaction(index, true_effect, size_group))) +
# geom_jitter(shape=16, alpha=0.02) +
geom_boxplot(outlier.shape = NA) +
facet_wrap(~index * outcome_type, scales = "free", ncol=2) +
theme_modern() +
scale_color_manual(values = c(`0` = "#f44336", `1` = "#8BC34A", name="Effect")) +
scale_fill_manual(values = c("p_value"="#607D8B", "p_MAP" = "#4CAF50", "p_direction" = "#2196F3",
"ROPE_89" = "#FFC107", "ROPE_95" = "#FF9800", "ROPE_full" = "#FF5722",
"p_ROPE"="#E91E63", "BF_log"="#9C27B0"), guide=FALSE) +
ylab("Index Value") +
xlab("Sample Size")
```
统计建模
我们拟合（频率）多元线性回归来预测不同指数的影响的存在与否，以及它们与噪声和样本大小的相互作用。
与频率p值的关系
```{r}
df %>%
select(outcome_type, true_effect, error, sample_size, p_value, p_direction, p_MAP, p_ROPE, ROPE_89, ROPE_95, ROPE_full, BF_log) %>%
gather(index, value, -error, -sample_size, -true_effect, -outcome_type, -p_value) %>%
mutate(true_effect = as.factor(true_effect),
index = factor(index, levels=c("p_direction", "p_MAP", "p_ROPE", "ROPE_89", "ROPE_95", "ROPE_full", "BF_log"))) %>%
mutate(temp = as.factor(cut(sample_size, 3, labels = FALSE))) %>%
group_by(temp) %>%
mutate(size_group = as.character(round(mean(sample_size)))) %>%
ungroup() %>%
ggplot(aes(x = p_value, y = value, color = true_effect, shape=size_group)) +
geom_point(alpha=0.025, stroke = 0, shape=16) +
facet_wrap(~index * outcome_type, scales = "free", ncol=2) +
theme_modern() +
scale_color_manual(values = c(`0` = "#f44336", `1` = "#8BC34A"), name="Effect") +
guides(colour = guide_legend(override.aes = list(alpha = 1)),
shape = guide_legend(override.aes = list(alpha = 1), title="Sample Size"))
```
与频率的基于p的任意聚类的关系
```{r}
df$sig_1 <- factor(ifelse(df$p_value >= .1, "n.s.", "-"), levels=c("n.s.", "-"))
df$sig_05 <- factor(ifelse(df$p_value >= .05, "n.s.", "*"), levels=c("n.s.", "*"))
df$sig_01 <- factor(ifelse(df$p_value >= .01, "n.s.", "**"), levels=c("n.s.", "**"))
df$sig_001 <- factor(ifelse(df$p_value >= .001, "n.s.", "***"), levels=c("n.s.", "***"))
get_data <- function(predictor, outcome, lbound=0, ubound=0.3){
fit <- suppressWarnings(glm(paste(outcome, "~ outcome_type * ", predictor), data=df, family = "binomial"))
# data <- data.frame(x=rep(1:100, 2))
data <- data.frame(outcome_type=rep(c("linear", "binary"), each=100))
data[predictor] <- rep(seq(lbound, ubound, length.out = 100), 2)
data$index <- predictor
predict_fit <- predict(fit, newdata=data, type="response", se.fit = TRUE)
data[outcome] <- predict_fit$fit
data$CI_lower <- predict_fit$fit - (qnorm(0.99) * predict_fit$se.fit)
data$CI_upper <- predict_fit$fit + (qnorm(0.99) * predict_fit$se.fit)
data <-
select(
data,
"value" := !!predictor,
!!outcome,
.data$outcome_type,
.data$index,
.data$CI_lower,
.data$CI_upper
)
return(data)
}
rbind(
rbind(
get_data(predictor="p_direction", outcome="sig_001", lbound=99.5, ubound=100),
get_data(predictor="p_MAP", outcome="sig_001", lbound=0, ubound=0.01),
get_data(predictor="p_ROPE", outcome="sig_001", lbound=97, ubound=100),
get_data(predictor="ROPE_89", outcome="sig_001", lbound=0, ubound=0.5),
get_data(predictor="ROPE_95", outcome="sig_001", lbound=0, ubound=0.5),
get_data(predictor="ROPE_full", outcome="sig_001", lbound=0, ubound=0.5),
get_data(predictor="BF_log", outcome="sig_001", lbound=0, ubound=10)
) %>%
rename("sig"=sig_001) %>%
mutate(threshold="p < .001"),
rbind(
get_data(predictor="p_direction", outcome="sig_01", lbound=98, ubound=100),
get_data(predictor="p_MAP", outcome="sig_01", lbound=0, ubound=0.1),
get_data(predictor="p_ROPE", outcome="sig_01", lbound=85, ubound=100),
get_data(predictor="ROPE_89", outcome="sig_01", lbound=0, ubound=2),
get_data(predictor="ROPE_95", outcome="sig_01", lbound=0, ubound=2),
get_data(predictor="ROPE_full", outcome="sig_01", lbound=0, ubound=2),
get_data(predictor="BF_log", outcome="sig_01", lbound=0, ubound=5)
) %>%
rename("sig"=sig_01) %>%
mutate(threshold="p < .01"),
rbind(
get_data(predictor="p_direction", outcome="sig_05", lbound=95, ubound=100),
get_data(predictor="p_MAP", outcome="sig_05", lbound=0, ubound=0.3),
get_data(predictor="p_ROPE", outcome="sig_05", lbound=50, ubound=100),
get_data(predictor="ROPE_89", outcome="sig_05", lbound=0, ubound=10),
get_data(predictor="ROPE_95", outcome="sig_05", lbound=0, ubound=10),
get_data(predictor="ROPE_full", outcome="sig_05", lbound=0, ubound=10),
get_data(predictor="BF_log", outcome="sig_05", lbound=0, ubound=2)
) %>%
rename("sig"=sig_05) %>%
mutate(threshold="p < .05"),
rbind(
get_data(predictor="p_direction", outcome="sig_1", lbound=90, ubound=100),
get_data(predictor="p_MAP", outcome="sig_1", lbound=0, ubound=0.5),
get_data(predictor="p_ROPE", outcome="sig_1", lbound=25, ubound=100),
get_data(predictor="ROPE_89", outcome="sig_1", lbound=0, ubound=20),
get_data(predictor="ROPE_95", outcome="sig_1", lbound=0, ubound=20),
get_data(predictor="ROPE_full", outcome="sig_1", lbound=0, ubound=20),
get_data(predictor="BF_log", outcome="sig_1", lbound=0, ubound=1)
) %>%
rename("sig"=sig_1) %>%
mutate(threshold="p < .1")
) %>%
mutate(index = as.factor(index)) %>%
ggplot(aes(x=value, y=sig)) +
geom_ribbon(aes(ymin=CI_lower, ymax=CI_upper), alpha=0.1) +
geom_line(aes(color=index), size=1) +
facet_wrap(~ index * threshold * outcome_type, scales = "free", ncol=8) +
theme_modern() +
scale_color_manual(values = c("p_value"="#607D8B", "p_MAP" = "#4CAF50", "p_direction" = "#2196F3",
"ROPE_89" = "#FFC107", "ROPE_95" = "#FF9800", "ROPE_full" = "#FF5722",
"p_ROPE"="#E91E63", "BF_log"="#9C27B0"), guide=FALSE) +
ylab("Probability of being significant") +
xlab("Index Value")
```
与等效性测试的关系
```{r}
df$equivalence_95 <- factor(ifelse(df$ROPE_95 == 0, "significant", "n.s."), levels=c("n.s.", "significant"))
df$equivalence_89 <- factor(ifelse(df$ROPE_89 == 0, "significant", "n.s."), levels=c("n.s.", "significant"))
rbind(
rbind(
get_data(predictor="p_direction", outcome="equivalence_95", lbound=97.5, ubound=100),
get_data(predictor="p_MAP", outcome="equivalence_95", lbound=0, ubound=0.2),
get_data(predictor="p_ROPE", outcome="equivalence_95", lbound=92.5, ubound=97.5),
get_data(predictor="ROPE_89", outcome="equivalence_95", lbound=0, ubound=3),
get_data(predictor="ROPE_full", outcome="equivalence_95", lbound=0, ubound=3),
get_data(predictor="BF_log", outcome="equivalence_95", lbound=0.5, ubound=2.5)
) %>%
rename("equivalence"=equivalence_95) %>%
mutate(level="95 HDI"),
rbind(
get_data(predictor="p_direction", outcome="equivalence_89", lbound=99.5, ubound=100),
get_data(predictor="p_MAP", outcome="equivalence_89", lbound=0, ubound=0.02),
get_data(predictor="p_ROPE", outcome="equivalence_89", lbound=0, ubound=100),
get_data(predictor="ROPE_95", outcome="equivalence_89", lbound=0, ubound=0.5),
get_data(predictor="ROPE_full", outcome="equivalence_89", lbound=0, ubound=0.5),
get_data(predictor="BF_log", outcome="equivalence_89", lbound=0, ubound=3)
) %>%
rename("equivalence"=equivalence_89) %>%
mutate(level="90 HDI")
) %>%
ggplot(aes(x=value, y=equivalence)) +
geom_ribbon(aes(ymin=CI_lower, ymax=CI_upper), alpha=0.1) +
geom_line(aes(color=index), size=1) +
facet_wrap(~ index * level * outcome_type, scales = "free", nrow=7) +
theme_modern() +
scale_color_manual(values = c("p_value"="#607D8B", "p_MAP" = "#4CAF50", "p_direction" = "#2196F3",
"ROPE_89" = "#FFC107", "ROPE_95" = "#FF9800", "ROPE_full" = "#FF5722",
"p_ROPE"="#E91E63", "BF_log"="#9C27B0"), guide=FALSE) +
ylab("Probability of rejecting H0 with the equivalence test") +
xlab("Index Value")
```
ROPE（完整）与p（方向）之间的关系
```{r}
df %>%
mutate(true_effect = as.factor(true_effect)) %>%
ggplot(aes(x=p_direction, y=ROPE_full, color=true_effect)) +
geom_point(alpha=0.025, stroke = 0, shape=16) +
facet_wrap(~ outcome_type, scales = "free", ncol=2) +
theme_modern() +
scale_color_manual(values = c(`0` = "#f44336", `1` = "#8BC34A"), name="Effect") +
ylab("ROPE (full)") +
xlab("Probability of Direction (pd)")
```
### 研究2：与采样特征的关系
### 研究3：与Priors规范的关系
讨论
基于多元线性回归的模拟，目前的工作旨在比较几个仅来自后验分布的效应存在指数，提供这些指数的“行为”与样本大小，噪声，先验和频率论的关系的视觉表示。 p值。
虽然这种与频率指数的比较可能看似违反直觉或错误（因为贝叶斯思想与频率论框架本质上不同），但我们认为这种比较因教学原因而引人注目。频率的p值对许多人“说话”，因此可以被视为参考和促进向贝叶斯框架转变的一种方式。然而，这并不排除在必要时寻求效应存在的一般范式的变化，并且贝叶斯指数与频率论p基本不同，而不仅仅是近似或等价。 至关重要的是，我们非常同意效果的存在与意义之间的区别和可能的分离。尽管如此，我们认为评估影响是否“有意义”在很大程度上依赖于文献，先验，新颖性，背景或领域，并且不能仅仅基于统计指标进行评估（即使某些指数，如ROPE相关的，尝试以有意义的方式弥合存在）。因此，研究人员应该依靠统计数据来评估效应的存在（以及大小和方向估计），并且系统地，但在上下文中，从更大的角度讨论其意义和重要性。
结论可以在指南部分找到。
## 报告指南
### 如何描述和报告模型的参数
贝叶斯分析返回每个参数（或效果）的后验分布。为了最低限度地描述这些分布，我们建议报告中心点的点估计以及表征估计不确定性（分散）的信息。另外，还可以报告效果存在和/或重要性的指数。基于先前对点估计和效应存在指数的比较 ，我们可以得出以下建议。
#### 确定性
我们建议将中位数报告为中心性指数，因为与平均值或MAP估计值相比，它更加稳健。然而，在严重偏斜的后部分布的情况下，MAP估计可能是一个很好的选择。
#### 不确定性
89％可信区间（CI）似乎是表征与估计相关的不确定性的合理范围，比更高的阈值（例如90％和95％）更稳定。 我们还建议基于HDI而不是分位数计算CI，其倾向于可能的值而不是中心值。注意，基于分位数（等尾区间）的CI在转换的情况下可能更合适（例如，当将对数转换为概率时）。 否则，最初未覆盖null的间隔可能会在转换后覆盖它。
#### 存在
贝叶斯框架允许整齐地描绘和量化假设检验的不同方面，例如效应存在和重要性。描述效应存在的最直接的指标是方向概率（pd），表示与效果的最可能方向（正面或负面）相关联的确定性。 该索引易于理解，易于解释，易于计算，对模型特征具有鲁棒性，并且与数据规模无关。
此外，它与频率论p值密切相关，因此可用于绘制相似之处，并为不熟悉贝叶斯统计的读者提供参考。 分别为.1， .05，.01和.001的双侧p值大致对应于95％，97.5％，99.5％和99.95％的pd 。 因此，为方便起见，我们建议将以下参考值作为解释助手：
* pd <= 95％ ~ p > .1：不确定
* pd > 95％ ~ p <.1：有可能存在
* pd > 97％ ：可能存在
* pd > 99％ ：极可能存在
* pd > 99.9％：肯定存在
#### 意义
ROPE中的百分比是显著性指数（在其主要含义中），告知我们参数是否与结果中的不可忽略的变化（就数量而言）相关 - 或不相关。 我们建议在ROPE中报告完全后验分布的百分比（完整的 ROPE）而不是给定比例的CI，这似乎更敏感（特别是描绘高度显著的影响）。我们建议使用百分比作为连续的显著性指数，而不是将其用作二进制，全有或全无的决策标准，例如原始等效性测试所建议的。 但是，根据模拟数据 ，我们建议将以下参考值作为解释助手：
* ROPE > 99％ ：可以忽略不计（我们可以接受零假设）
* ROPE > 97.5％ ：可能微不足道
* ROPE<= 97.5％ ＆ > = 2.5％ ：未定意义
* ROPE <2.5％ ：可能很重要
* ROPE <1％ ：显著（我们可以拒绝零假设）
请注意，需要格外小心，因为它的解释很大程度上取决于其他参数，如样本大小和ROPE范围（见此处 ） 。
#### 模板
基于这些建议，基于其后验分布的参数的最小报告的模板句子可以是：“X的影响具有pd的概率为负,中位数= ，89％CI[HDI低，HDI高]并且可以被认为是显著的（ROPE% in ROPE）。”
### 如何比较不同的模型
尽管它也可用于评估效应的存在和重要性，贝叶斯因子（BF）是一种通用的指数，可用于直接比较不同的模型（或数据生成过程）。贝叶斯因子是一个比率，通知我们观察到的数据在两个比较模型下的可能性有多大（或更少）-通常是具有效果的模型与没有效果的模型。根据空模型的规范（无论是点估计（例如，0）还是间隔），贝叶斯因子可以在效应存在和重要性的上下文中使用。
一般来说，贝叶斯因子大于1，证明有利于其中一个模型，而贝叶斯因子小于1，证明有利于另一个模型。 有几条经验法则可以帮助解释（见这里），>3是一个常用的门槛，可以对非轶事证据进行分类。
#### 模板
报告贝叶斯因子（BF）时，可以使用以下句子：“有中等证据支持x（BF = BF） 不存在影响。”
title: "贝叶斯分析Article"
author: "xuefliang"
date: "7/11/2019"
output:
html_document:
theme: readable
highlight: haddock
df_print: paged
code_folding: show
toc: true
toc_float: TRUE
number_sections: false
fig_width: 10.08
fig_height: 6
tags: [r, bayesian, posterior, test]
editor_options:
chunk_output_type: console
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE,warning = FALSE,message = FALSE,comment=">")
options(knitr.kable.NA = '',digits=2)
set.seed(333)
library(tidyverse)
library(bayestestR)
library(rstanarm)
library(insight)
```
## 可信区间（CI）
### 什么是可信区间？
可信区间是贝叶斯统计中的一个重要概念。 其核心目的是描述和总结与您的参数相关的不确定性 。 在这方面，它可能看起来与频率主义的置信区间很相似。然而，虽然他们的目标相似，但他们的统计定义和意义却截然不同。实际上，虽然后者是通过一个充满罕见测试假设和近似的复杂算法获得的，但可信区间的计算相当简单。
由于贝叶斯推断返回了可能的效应值（后验）的分布，可信区间只是包含特定百分比的可能值的范围。 例如，95％可信区间仅仅是后验分布的中心部分，其包含95％的值。
请注意，与频率主义区间相比，这大大提高了贝叶斯区间的可解释性。事实上，贝叶斯框架允许说“给定观察到的数据，效果有95％的概率落在这个范围内” ，而频率较低的直接替代方案（95％ 置信区间 ）将“ 有95％的可能性当从这种数据计算置信区间时，效果落在这个范围内 “。
### 为什么默认为89％？
当然，当选择CI级别进行默认报告时，人们开始使用95％，这是在频率论世界中使用的任意约定。 然而，一些作者认为95％可能不是贝叶斯后验分布最合适的，如果没有足够的后验样本可能缺乏稳定性（Kruschke 2014） 。该命题是使用90％而不是95％。 然而，最近，McElreath（2014年，2018年）建议，如果我们首先使用任意阈值，为什么不使用89％，因为这个值具有作为素数的附加参数。
因此，默认情况下，CI以89％的间隔（ci=0.89）计算，被认为比例如95％的间隔更稳定（Kruschke 2014）。如果应计算95％的间隔，则建议有效样本量至少为10.000（Kruschke，2014，p.183ff）。 此外，89是最高素数，不超过已经不稳定的95％阈值。什么与任何事情有关？ 没什么 ，但它提醒我们任何这些约定的完全任意性（McElreath 2014，@ mcelreath2018statistical） 。
### 不同类型的CI
读者可能会注意到bayestestR提供了两种计算可信区间的方法，即最高密度区间（HDI）（hdi()）和等尾区间（ETI） （eti()）。也可以通过ci()函数的method参数更改这些方法。有什么不同？ 让我们来看看：
```{r}
# Generate a normal distribution
posterior <- distribution_normal(1000)
# Compute HDI and ETI
ci_hdi <- ci(posterior, method = "HDI")
ci_eti <- ci(posterior, method = "ETI")
# Plot the distribution and add the limits of the two CIs
posterior %>%
estimate_density(extend=TRUE) %>%
ggplot(aes(x=x, y=y)) +
geom_area(fill="orange") +
theme_classic() +
# HDI in blue
geom_vline(xintercept=ci_hdi$CI_low, color="royalblue", size=3) +
geom_vline(xintercept=ci_hdi$CI_high, color="royalblue", size=3) +
# Quantile in red
geom_vline(xintercept=ci_eti$CI_low, color="red", size=1) +
geom_vline(xintercept=ci_eti$CI_high, color="red", size=1)
```
这些几乎完全一样......
但其他类型的分布也是如此吗？
```{r}
# Generate a beta distribution
posterior <- distribution_beta(1000, 6, 2)
# Compute HDI and Quantile CI
ci_hdi <- ci(posterior, method = "HDI")
ci_eti <- ci(posterior, method = "ETI")
# Plot the distribution and add the limits of the two CIs
posterior %>%
estimate_density(extend=TRUE) %>%
ggplot(aes(x=x, y=y)) +
geom_area(fill="orange") +
theme_classic() +
# HDI in blue
geom_vline(xintercept=ci_hdi$CI_low, color="royalblue", size=3) +
geom_vline(xintercept=ci_hdi$CI_high, color="royalblue", size=3) +
# Quantile in red
geom_vline(xintercept=ci_eti$CI_low, color="red", size=1) +
geom_vline(xintercept=ci_eti$CI_high, color="red", size=1)
```
这个差异很大。
与HDI相反，对于该HDI ，区间内的所有点具有比区间外的点更高的概率密度， ETI是等尾的 。 这意味着90％的间隔在其限制的任一侧具有5％的分布。 它表示第5百分位数和95h百分位数。 在对称分布中，计算可靠区间的两种方法（ETI和HDI）返回类似的结果。
对于偏斜的分布，情况并非如此。实际上，ETI中的参数值可能比ETI之外的参数值具有更低的可信度（更不可能）。 作为分布中可靠值的摘要，此属性似乎不合需要。
另一方面，当将变换应用于分布时，ETI范围确实会发生变化（例如，对于概率的对数几率缩放）：变换分布的下限和上限将对应于变换的下限和上限。相反，将变换应用于分布将改变最终的HDI。
## 方向概率（pd）
### 什么是pd？
方向概率（pd）是效应存在的指数，范围从50％到100％，表示效果在特定方向上的确定性（即 正或负）。除了解释，理解和计算的简单性之外，该索引还提供了其他有趣的属性：
它独立于模型 ：它完全基于后验分布，不需要来自数据或模型的任何附加信息。它对响应变量和预测变量的规模都很稳健。它与频率论p值密切相关，因此可用于绘制相似之处，并为不熟悉贝叶斯统计的读者提供参考。然而，该指数与评估影响的程度和重要性无关，通过其他指数（如ROPE百分比）可以更好地实现。事实上，重要性和存在的指数是完全独立的。您可以使用99.99％的pd效果，其中整个后验分布集中在[0.0001,0.002]范围内。在这种情况下，效果是肯定的，具有高确定性 ，但也不显著 （即非常小）。
有效指数存在，例如pd ，在探索性研究或临床研究中特别有用，其重点是确保感兴趣的效果不是相反的方向（对于临床研究，如治疗无害）。然而，一旦效果的方向得到确认，焦点就应转向其重要性，包括对其大小，相关性和重要性的精确估计。
### 与p值的关系
在大多数情况下，似乎pd通过以下公式与频率为单侧的p值直接对应：$\ [p_ {one sided} = 1-p_d / 100 ]$同样， 双边p值 （最常报道的一个）通过以下公式是等价的： $\ [p_ {双面} = 2 *（1-p_d / 100）]$因此，双边p值分别为.1 ， .05 , . 01和.001大致对应于95％ ， 97.5％ ， 99.5％和99.95％的pd 。 但如果它像p值一样，那一定是坏的，因为p值很差[ 参见重复性危机]。
事实上，重复性危机这一方面可能被误解了。实际上，并不是p值是一个本质上错误或错误的指数。 相反，它的误用 ，误解和误解会助长局势的衰退。例如，pd与p值高度相关的事实表明后者更多地是效应存在的指数而非重要性（即，“感兴趣的价值”）。贝叶斯版本pd具有直观的含义，并且明显表明所有阈值都是任意的。然而，pd的数学和解释透明度及其作为效应存在指数的重新概念化，为贝叶斯结果的表征提供了有价值的见解。此外，它与频率主义p值的接近使其成为一个完美的衡量标准，可以简化心理学研究向贝叶斯框架的采用。
### 计算方法
计算pd的最简单和直接的方法是1）查看中位数的符号，2）选择相同符号的后部分，3）计算该部分代表的百分比。 这种“简单”的方法是最困难的方法，但其精确度与后验抽取的数量直接相关。第二种方法依赖于密度估计。它首先估算密度函数（有许多方法可用），然后计算0的另一侧密度曲线的曲线下面积（AUC）。基于密度的方法可以假设更精确，但很大程度上取决于用于估算密度函数的方法。
### 方法比较
让我们比较4种可用方法，直接方法和3种基于密度的方法，它们的密度估计算法不同（见estimate_density ）。
### 关联
让我们首先测试通过不同方法获得的结果的接近度和相似性。
```{r}
library(bayestestR)
library(logspline)
library(KernSmooth)
# Compute the correlations
data <- data.frame()
for(the_mean in runif(25, 0, 4)){
for(the_sd in runif(25, 0.5, 4)){
x <- rnorm(100, the_mean, abs(the_sd))
data <- rbind(data,
data.frame("direct" = pd(x),
"kernel" = pd(x, method="kernel"),
"logspline" = pd(x, method="logspline"),
"KernSmooth" = pd(x, method="KernSmooth")
))
}
}
data <- as.data.frame(sapply(data, as.numeric))
# Visualize the correlations
library(GGally)
GGally::ggpairs(data) +
theme_classic()
```
所有给出的方法都是高度相关的，并给出非常相似的结果。这意味着方法选择并不激烈，不能用于过多地调整结果。
### 准确性
为了测试每种方法的准确性，我们将从非常密集的分布（具有大量观察值）计算直接pd开始。 这将是我们的基线，或“真实”pd。然后，我们将迭代地从该父分布中绘制较小的样本，并且我们将使用不同的方法计算pd 。 这个估计值与参考值越接近越好。
```{r}
data <- data.frame()
for(i in 1:25){
the_mean <- runif(1, 0, 4)
the_sd <- abs(runif(1, 0.5, 4))
parent_distribution <- rnorm(100000, the_mean, the_sd)
true_pd <- pd(parent_distribution)
for(j in 1:25){
sample_size <- round(runif(1, 25, 5000))
subsample <- sample(parent_distribution, sample_size)
data <- rbind(data,
data.frame("sample_size" = sample_size,
"true" = true_pd,
"direct" = pd(subsample) - true_pd,
"kernel" = pd(subsample, method="kernel")- true_pd,
"logspline" = pd(subsample, method="logspline") - true_pd,
"KernSmooth" = pd(subsample, method="KernSmooth") - true_pd
))
}
}
data <- as.data.frame(sapply(data, as.numeric))
data %>%
tidyr::gather(Method, Distance, -sample_size, -true) %>%
ggplot(aes(x=sample_size, y = Distance, color = Method, fill= Method)) +
geom_point(alpha=0.3, stroke=0, shape=16) +
geom_smooth(alpha=0.2) +
geom_hline(yintercept=0) +
theme_classic() +
xlab("\nDistribution Size")
```
基于“核心”的密度方法似乎一直低估了pd 。有趣的是，即使在少量后抽取的情况下，“直接”方法看起来也更可靠。
## 实用对等区（ROPE）
### 什么是ROPE？
与频率论方法不同，贝叶斯推断不是基于统计显著性，其中效果针对“零”进行测试。实际上，贝叶斯框架提供了参数的概率视图，允许评估与它们相关的不确定性。因此，我们可以得出结论，在超出特定范围之外的可能被视为“实际上没有效果”（即，可忽略不计的幅度）的概率就足够了，而不是断定效果存在。 该范围称为实际等效区域（ROPE） 。
实际上，从统计上来说，后验分布与0不同的概率没有多大意义（它与单个点的无限概率不同）。 因此，强调ROPE的想法是让用户定义空值周围的区域，其中包含等于实际目的的空值的值（Kruschke 2010,2014; Kruschke，Aguinis和Joo 2012） 。
### 等效性测试
ROPE是对应于“无效”假设的区域，用于等效性测试，以测试参数是否显著（在重要程度上足以被关注）。 该测试通常基于“HDI + ROPE决策规则” 来检查是否应该针对明确制定的“零假设”（ 即 ，ROPE）接受或拒绝参数值。换句话说，它检查作为空区域（ROPE）的可信区间（CI）的百分比。 如果该百分比足够低，则拒绝原假设。如果该百分比足够高，则接受零假设。ROPE中的
### 可信区间与ROPE中的完全后验区间使用ROPE和HDI作为可信区间
Kruschke（2018）建议使用落在ROPE内的95％HDI的百分比作为决策规则。然而，由于89％的HDI被认为是更好的选择,因此bayestestR默认提供了落在ROPE内的89％HDI的百分比。然而，初步数据倾向于表明使用完全后验分布的百分比而不是CI可能更敏感（特别是描绘非常显著的影响）。因此，我们建议用户考虑使用完整的ROPE百分比（通过设置ci=1），这将返回ROPE中整个后验分布的部分。
### ROPE接受或拒绝的百分比是多少？
如果HDI完全在ROPE之外，则该参数的“零假设”被“拒绝”。如果ROPE完全涵盖HDI，即参数的所有最可靠值都在实际等价区域内，则接受零假设。 否则，尚未决定是接受还是拒绝零假设。
如果使用完整的ROPE（即100％的HDI），则如果ROPE中后验的百分比小于2.5％或大于97.5％，则拒绝或接受零假设。 理想的结果是ROPE内部的比例较低（越接近于零越好）。
### 如何定义ROPE范围？
Kruschke（2018）提出，默认情况下，ROPE可以设置为标准化参数的-0.1到0.1的范围（根据Cohen1988可忽略的影响大小）。
* 对于线性模型（lm） ，这可以推广为：$\ [ - 0.1 * SD_ {y}，0.1 * SD_ {y}] $ 。
* 对于逻辑模型 ，以对数优势比表示的参数可通过以下公式转换为标准差： $\ [\ sqrt {3} / \ \pi ] $，从而产生-0.055到-0.055的范围。对于具有二进制结果的其他模型，强烈建议手动指定rope参数。 目前，对于逻辑模型应用相同的默认值。
* 对于t检验 ，使用响应的标准偏差，类似于线性模型（见上文）。
* 对于相关性 ，使用-0.05, 0.05 即 Cohen（1988）经验法则所建议的可忽略相关值的一半。
* 对于所有其他模型， -0.1, 0.1用于确定ROPE限制，但强烈建议手动指定。
### 对参数刻度的灵敏度
当使用基于ROPE的指数时，考虑预测变量的单位（即尺度）是重要的，因为ROPE的正确解释表示实际等效于零的区域取决于预测变量的尺度。实际上，与其他指数（例如pd）不同，ROPE中的百分比取决于其参数的单位。换句话说，由于ROPE表示响应比例的固定部分，因此其与系数的接近程度取决于系数本身的比例。
例如，如果我们考虑一个简单的回归growth ~ time，模拟Wookies婴儿的发展，可忽略的变化（ROPE）小于54厘米。如果我们的time变量用天表示 ，我们会发现系数（代表白天的增长）大约是10厘米（系数后验的中位数是10）。我们认为这可以忽略不计。但是，如果我们决定以年为单位表示time变量，系数将通过此变换进行缩放（因为它现在将代表按年增长）。系数现在约为3550厘米（10*355），我们现在认为这是重要的 。
```{r}
library(see)
data <- iris # Use the iris data
model <- stan_glm(Sepal.Length ~ Sepal.Width, data=data) # Fit model
# Compute indices
pd <- p_direction(model)
percentage_in_rope <- rope(model, ci=1)
# Visualise the pd
plot(pd)
pd
# Visualise the percentage in ROPE
plot(percentage_in_rope)
percentage_in_rope
```
我们可以看出， Sepal.Length和Sepal.Width之间线性关系的pd和ROPE百分比分别约为92.95%和15.95% ，对应于不确定和不显著的影响。 如果我们扩展预测器会发生什么？
```{r}
data$Sepal.Width_scaled <- data$Sepal.Width / 100 # Divide predictor by 100
model <- stan_glm(Sepal.Length ~ Sepal.Width_scaled, data=data)
# Compute indices
pd <- p_direction(model)
percentage_in_rope <- rope(model, ci=1)
# Visualise the pd
plot(pd)
pd
# Visualise the percentage in ROPE
plot(percentage_in_rope)
percentage_in_rope
```
正如您所看到的，通过简单地将预测器除以100，我们彻底改变了与ROPE中的百分比相关的结论（其变得非常接近0）：现在可以将效果解释为重要。因此，我们建议在选择ROPE范围时密切关注预测变量的单位（例如，哪个系数对应于一个小影响？），以及报告或读取ROPE结果时。
### 多重共线性：非独立协变量
当参数显示强相关性时，即，当协变量不是独立的时，联合参数分布可以朝向或远离ROPE移位。 共线性使基于单变量边缘的ROPE和假设检验无效，因为概率是以独立性为条件的。最有问题的是仅与ROPE区域部分重叠的参数。在共线性的情况下，这些参数的（联合）分布可以获得增加或减少的ROPE，这意味着基于ROPE的推断是不合适的（Kruschke 2014）。equivalence_test()和rope()函数对参数之间的成对相关性执行简单检查，但由于两个以上变量之间可能存在共线性，因此检查此假设检验假设的第一步是查看不同的对图。 更复杂的检查是投影预测变量选择（Piironen和Vehtari 2017） 。
## 贝叶斯因子
贝叶斯应用统计框架的采用，特别是在社会或心理科学领域，似乎正朝着两个不同的方向发展。 标记他们分离的关键主题之一是他们对贝叶斯因子的看法。简而言之，一些作者（例如，由Wagenmakers领导的“阿姆斯特丹学校”）提倡使用并强调其作为统计指标的质量，而其他作者则指出其局限性，而更倾向于对后验分布的精确描述（使用CI）， ROPEs等）。
bayestestR在本次辩论中没有参与，而是提供工具来帮助您进行任何想要实现的分析。相反，它强烈支持一个明智的选择的概念： 发现方法，尝试它们，理解它们，了解它们，并自己决定 。
话虽如此，这里是贝叶斯因子的介绍:
### 贝叶斯因子
贝叶斯因子（BF）是一个“模型”相对于另一个“模型”的相对证据的指数，其在贝叶斯推断中用作经典（频率论）假设检验指数（例如p-值）的替代。在贝叶斯框架中，贝叶斯因子也可以被认为是根据观察到的数据更新关于两个模型的相对几率的一些先验信念的量。
根据贝叶斯定理：
$$[P（M | D）= \frac {P（D | M）P（M）} {P（D）} ]$$
然后通过比较两个模型，我们得到：
$$ [\frac {P（M_1 | D）} {P（M_2 | D）} = \frac {P（D | M_1）} {P（D | M_2）} \times \frac {P（M_1）} { P（M_2）} ]$$
其中中间项是贝叶斯因子：
$$[BF_ {12} = \frac {P（D | M_1）} {P（D | M_2）}]$$ 这些通常按比例计算两个竞争假设/模型的边际可能性，但正如我们从上面的等式中可以看到的，它们也可以通过将后验概率除以先验概率来计算。重要的是，贝叶斯因子涵盖了广泛的指数和应用 ，并有不同的风格。
### 比较模型
一般来说，贝叶斯因子用于比较模型并回答问题：在哪种模型下观察到的数据更可能？换句话说，哪种模型更有可能产生观察到的数据？这通常通过计算两个模型的边际可能性来完成。在这种情况下，贝叶斯因子是一个比较模型相对于另一个模型的相对证据的度量。
### 贝叶斯模型（ brms和rstanarm ）
注意：为了计算模型的贝叶斯因子，必须在拟合时添加非默认参数：
* brmsfit模型必须符合save_all_pars = TRUE
* stanreg模型必须已经安装了定义的diagnostic_file 。
让我们首先用5个贝叶斯回归与brms来预测Sepal.Length ：
```{r}
library(brms)
m0 <- brm(Sepal.Length ~ 1, data = iris, save_all_pars = TRUE)
m1 <- brm(Sepal.Length ~ Petal.Length, data = iris, save_all_pars = TRUE)
m2 <- brm(Sepal.Length ~ Species, data = iris, save_all_pars = TRUE)
m3 <- brm(Sepal.Length ~ Species + Petal.Length, data = iris, save_all_pars = TRUE)
m4 <- brm(Sepal.Length ~ Species * Petal.Length, data = iris, save_all_pars = TRUE)
```
我们现在可以将这些模型与bayesfactor_models()函数进行比较，使用denominator参数指定将与所有模型进行比较的模型（在本例中为仅截距模型）：
```{r}
comparison <- bayesfactor_models(m1, m2, m3, m4, denominator = m0)
comparison
```
我们可以看到完整模型是最好的模型 - 与$（BF _ {\text {m0}} = 9 \times 10 ^ {55} ）$与null（仅截距）相比。我们还可以将参考模型更改为主效应模型：
```{r}
update(comparison, reference = 3)
```
我们可以看到，尽管完整模型是最好的，但几乎没有证据表明它比主效应模型更可取。我们还可以将参考模型更改为Species模型：
```{r}
update(comparison, reference = 2)
```
请注意，在贝叶斯框架中，比较模型不需要是嵌套模型，就像我们将Petal.Length模型与Species Petal.Length模型进行比较时所发生的那样（在频率论者框架中无法做到的事情，比较模型必须互相嵌套）。
注意：为了正确和精确地估计贝叶斯因子，你总是需要4P: Proper Priors (1, 2, 3), and a Plentiful Posterior (4)。
### 频率论模型的BIC近似
也可以为频率模型计算贝叶斯因子。这是通过比较BIC测量来完成的，允许对非嵌套频率模型进行贝叶斯比较（Wagenmakers 2007） 。 让我们试试一些线性混合模型 ：
```{r}
library(lme4)
m0 <- lmer(Sepal.Length ~ (1 | Species), data = iris)
m1 <- lmer(Sepal.Length ~ Petal.Length + (1 | Species), data = iris)
m2 <- lmer(Sepal.Length ~ Petal.Length + (Petal.Length | Species), data = iris)
m3 <- lmer(Sepal.Length ~ Petal.Length + Petal.Width + (Petal.Length | Species), data = iris)
m4 <- lmer(Sepal.Length ~ Petal.Length * Petal.Width + (Petal.Length | Species), data = iris)
bayesfactor_models(m1, m2, m3, m4, denominator = m0)
```
### 通过贝叶斯模型平均包含贝叶斯因子
包含贝叶斯因子回答了这个问题：在具有特定效果的模型下观察到的数据是否比在没有特定效果的模型下更可能？ 换句话说，平均而言-具有效果的模型$X$更有可能产生观察到的数据而不是没有效果的模型$X$？让我们使用上面的brms示例：
```{r}
bayesfactor_inclusion(comparison)
```
如果我们检查交互项的包含贝叶斯因子，我们可以看到，在所有5个模型中，具有交互项（Species:Petal.Length）的模型比没有交互项的模型的可能性高5倍。
我们也可以仅比较匹配的模型-即，没有效果的模型$A$将仅与具有效果$A$的模型进行比较，但不与具有更高级别交互的模型进行比较（请参阅解释为什么您可能需要在这里这样做）。
```{r}
bayesfactor_inclusion (comparison, match_models = TRUE )
```
#### 与JASP比较
bayesfactor_inclusion()旨在提供跨模型平均值的贝叶斯因子，类似于JASP的效果选项。 让我们比较两者：
比较所有模型
```{r}
library(BayesFactor)
ToothGrowth$dose <- as.factor(ToothGrowth$dose)
BF_ToothGrowth <- anovaBF(len ~ dose*supp, ToothGrowth)
bayesfactor_inclusion(BF_ToothGrowth)
```
比较匹配的模型
```{r}
bayesfactor_inclusion(BF_ToothGrowth, match_models = TRUE)
```
使用Nuisance Effects
我们将在JASP中为null模型添加剂量，并在R执行相同的操作：
```{r}
BF_ToothGrowth_against_dose <- BF_ToothGrowth[3:4]/BF_ToothGrowth[2] # OR:
# update(bayesfactor_models(BF_ToothGrowth), subset = c(4,5), reference = 3)
BF_ToothGrowth_against_dose
bayesfactor_inclusion(BF_ToothGrowth_against_dose)
```
### Savage-Dickey密度比率贝叶斯因子
Savage-Dickey密度比可用于回答以下问题：鉴于观察到的数据，null更多还是更不可能？
这是通过比较先前和后验分布之间的空值密度来完成的，并且是比较提供的（点）空模型的贝叶斯因子的近似值：
$H_0$ 与$H_1$的贝叶斯因子可以通过分析积分模型参数$\theta$获得 。 然而，贝叶斯因子同样可以通过仅考虑$H_1$，并且将后部的高度除以$\theta$的高度除以先前的高度$\theta$ ，兴趣。“ （Wagenmakers et al.2010）
让我们使用学生睡眠数据，并尝试回答这个问题： 根据观察到的数据，药物（变量group ）或多或少可能对额外睡眠时数（可变extra）没有影响？
blopplot表明第二组的额外睡眠时间更长 。 多少钱？ 让我们拟合一个简单的贝叶斯线性模型 。
```{r}
library(rstanarm)
model <- stan_glm(extra ~ group, data = sleep)
```
我们可以在此模型上使用as.data.frame()来提取与group2的效果相关的后验分布，以及来自insight包的 get_priors() ，以查看使用的先前分布：
```{r}
posterior <- as.data.frame(model)$group2
insight::get_priors(model)
```
对于group2参数，使用的先验是平均值（位置）0和SD（比例）5.044799。我们可以模拟这个先前的分布如下：
```{r}
prior <- distribution_normal(length(posterior), mean = 0, sd = 5.044799)
```
我们现在可以绘制先前和后验分布：观察分布，我们可以看到后验中心位于1.56。 但这并不意味着0的影响必然不太可能。 为了测试它，我们将使用bayesfactor_savagedickey() 。
```{r}
test_group2 <- bayesfactor_savagedickey(posterior = posterior, prior = prior)
test_group2
```
该BF表示在给定数据的情况下，0（点空效应模型）的效果不太可能为0.88 。 或者，当解释为$H_0$与$H_1$的贝叶斯因子时，我们可以说观察到的数据在空值下比在具有效果的模型下更可能是1 / 0.88=1.14倍！因此，尽管分布中心偏离了零点，但它仍然在零点周围非常密集。请注意，可以在此处找到贝叶斯因子的解释指南 。
#### 定向测试
如果我们先前有关于效果方向的假设，我们也可以进行方向测试（“单侧”或“单尾”测试）。这是通过对先前和后验分布设置订单限制来完成的（Morey和Wagenmakers 2014; Morey 2015） 。
```{r}
test_group2_right <- bayesfactor_savagedickey(posterior = posterior, prior = prior, direction = ">")
test_group2_right
```
#### 测试所有模型参数
或者，我们也可以将模型直接传递给bayesfactor_savagedickey()以同时测试所有模型的参数：
```{r}
bayesfactor_savagedickey(model)
```
#### 测试对比（与emmeans ）
我们也可以使用bayesfactor_savagedickey()和emmeans ，让我们测试贝叶斯对比 。
```{r}
library(emmeans)
group_diff <- pairs(emmeans(model, ~ group))
group_diff
bayesfactor_savagedickey(group_diff, prior = model)
```
## 深度1：点估计的比较
### 贝叶斯框架中的效应点估计
### 介绍
贝叶斯框架和频率框架之间的主要区别之一是前者返回每个效应的概率分布（即模型的感兴趣参数，例如回归斜率）而不是单个值。然而，对于报告或用于进一步分析的需求和需求，仍然需要能够最好地表征潜在后验分布的单个值（点估计）。在文献中使用三个主要指数用于效果估计：平均值，中值或MAP（最大后验）估计（大致对应于分布的模式“峰值”）。不幸的是，没有就使用哪一个达成共识，因为从未进行过系统的比较。在目前的工作中，我们将比较这三个点之间的效果估计值，以及从可比较的频率论模型中提取的广为人知的beta 。 通过这种比较，我们期望在两个框架之间架起桥梁和关系，帮助和缓解公众的过渡。
### 实验1：与误差（噪声）和样本大小的关系
#### 方法
模拟旨在调整以下特征：
模型的类型 ：线性或逻辑。
“真实”效果 （绘制数据的原始回归系数）：可以是1或0（无效果）。
样本量 ：从20到100步长为10。
误差 ：应用于预测器的高斯噪声，SD均匀分布在0.33和6.66之间（具有1000个不同的值）。
我们为这些特征的每个组合生成了一个数据集，从而产生了总共2*2*9*1000=36000贝叶斯和频率模型。这里可以使用用于生成的代码（请注意，通常需要几天/几周才能完成）。
```{r}
library(see)
df <- read.csv("https://raw.github.com/easystats/circus/master/data/bayesSim_study1.csv")
```
#### 结果
对噪音的敏感度
```{r}
df %>%
select(error, true_effect, outcome_type, beta, Median, Mean, MAP) %>%
gather(estimate, value, -error, -true_effect, -outcome_type) %>%
mutate(temp = as.factor(cut(error, 10, labels = FALSE))) %>%
group_by(temp) %>%
mutate(error_group = round(mean(error), 1)) %>%
ungroup() %>%
filter(value < 6) %>%
ggplot(aes(x = error_group, y = value, fill = estimate, group = interaction(estimate, error_group))) +
# geom_hline(yintercept = 0) +
# geom_point(alpha=0.05, size=2, stroke = 0, shape=16) +
# geom_smooth(method="loess") +
geom_boxplot(outlier.shape=NA) +
theme_modern() +
scale_fill_manual(values = c("beta" = "#607D8B", "MAP" = "#795548", "Mean" = "#FF9800", "Median" = "#FFEB3B"),
name = "Index") +
ylab("Point-estimate") +
xlab("Noise") +
facet_wrap(~ outcome_type * true_effect, scales="free")
```
对样本量的敏感度
```{r}
df %>%
select(sample_size, true_effect, outcome_type, beta, Median, Mean, MAP) %>%
gather(estimate, value, -sample_size, -true_effect, -outcome_type) %>%
mutate(temp = as.factor(cut(sample_size, 10, labels = FALSE))) %>%
group_by(temp) %>%
mutate(size_group = round(mean(sample_size))) %>%
ungroup() %>%
filter(value < 6) %>%
ggplot(aes(x = size_group, y = value, fill = estimate, group = interaction(estimate, size_group))) +
# geom_hline(yintercept = 0) +
# geom_point(alpha=0.05, size=2, stroke = 0, shape=16) +
# geom_smooth(method="loess") +
geom_boxplot(outlier.shape=NA) +
theme_modern() +
scale_fill_manual(values = c("beta" = "#607D8B", "MAP" = "#795548", "Mean" = "#FF9800", "Median" = "#FFEB3B"),
name = "Index") +
ylab("Point-estimate") +
xlab("Sample size") +
facet_wrap(~ outcome_type * true_effect, scales="free")
```
统计建模
我们拟合（频率论）多元线性回归来统计测试预测效果的存在与否，以及它们与噪声和样本大小的相互作用。
```{r}
df %>%
select(sample_size, error, true_effect, outcome_type, beta, Median, Mean, MAP) %>%
gather(estimate, value, -sample_size, -error, -true_effect, -outcome_type) %>%
glm(true_effect ~ outcome_type / estimate / value, data=., family="binomial") %>%
broom::tidy() %>%
select(term, estimate, p=p.value) %>%
filter(stringr::str_detect(term, 'outcome_type'),
stringr::str_detect(term, ':value')) %>%
arrange(desc(estimate)) %>%
knitr::kable(digits=2)
```
这表明，与频率论者的beta相比，为了描述效果的存在与否;
* 对于线性模型， 均值是更好的预测因子，紧随其后的是中位数 ， MAP和频率主义beta 。
* 对于逻辑模型， MAP是更好的预测因子，其次是中位数 ， 均值和频率测试 。
总的来说，中位数似乎是一个恰当的选择，在不同类型的模型中保持高性能。
### 实验2：与采样特征的关系
#### 方法
模拟旨在调整以下特征：
* 模型类型 ：线性或逻辑。
* “真实”效果 （绘制数据的原始回归系数）：可以是1或0（无效果）。
* 绘制 ：从5到1000（步长为5）（1000次迭代）。
* 预热 ：预热 迭代的比率。 从1/10到9/10，步长为0.1（9次迭代）。
我们为这些特征的每个组合生成了3个数据集，从而得到总共2*2*8*40*9*3=34560贝叶斯和频率模型。这里可以使用用于生成的代码（请注意，通常需要几天/几周才能完成）。
```{r}
df <- read.csv("https://raw.github.com/easystats/circus/master/data/bayesSim_study2.csv")
```
#### 结果
对迭代次数的敏感性
```{r}
df %>%
select(iterations, true_effect, outcome_type, beta, Median, Mean, MAP) %>%
gather(estimate, value, -iterations, -true_effect, -outcome_type) %>%
mutate(temp = as.factor(cut(iterations, 5, labels = FALSE))) %>%
group_by(temp) %>%
mutate(iterations_group = round(mean(iterations), 1)) %>%
ungroup() %>%
filter(value < 6) %>%
ggplot(aes(x = iterations_group, y = value, fill = estimate, group = interaction(estimate, iterations_group))) +
geom_boxplot(outlier.shape=NA) +
theme_classic() +
scale_fill_manual(values = c("beta" = "#607D8B", "MAP" = "#795548", "Mean" = "#FF9800", "Median" = "#FFEB3B"),
name = "Index") +
ylab("Point-estimate of the true value 0\n") +
xlab("\nNumber of Iterations") +
facet_wrap(~ outcome_type * true_effect, scales="free")
```
对预热比率的敏感度
```{r}
df %>%
mutate(warmup = warmup / iterations) %>%
select(warmup, true_effect, outcome_type, beta, Median, Mean, MAP) %>%
gather(estimate, value, -warmup, -true_effect, -outcome_type) %>%
mutate(temp = as.factor(cut(warmup, 3, labels = FALSE))) %>%
group_by(temp) %>%
mutate(warmup_group = round(mean(warmup), 1)) %>%
ungroup() %>%
filter(value < 6) %>%
ggplot(aes(x = warmup_group, y = value, fill = estimate, group = interaction(estimate, warmup_group))) +
geom_boxplot(outlier.shape=NA) +
theme_classic() +
scale_fill_manual(values = c("beta" = "#607D8B", "MAP" = "#795548", "Mean" = "#FF9800", "Median" = "#FFEB3B"),
name = "Index") +
ylab("Point-estimate of the true value 0\n") +
xlab("\nNumber of Iterations") +
facet_wrap(~ outcome_type * true_effect, scales="free")
```
### 实验3：与Priors规范的关系
结论可以在指南部分找到。
## 深度2：效应存在指数和显著指数的比较
### 贝叶斯框架中效应存在的指标及其意义
现在普遍认为贝叶斯统计框架是进行心理科学研究的正确方法。然而，它的灵活性是它的力量和弱点，因为没有就应该计算或报告哪些指数达成一致。此外，缺乏一个共识的效应存在指数，例如频率论的p值，可能导致许多非熟悉的读者在贝叶斯统计中感知到的不必要的模糊。因此，本研究描述并比较了几种效应存在指数，提供了这些指数“行为”的直观视觉表示，与传统指标（如样本量和频率显著性）相关。结果有助于培养对研究人员报告的价值观的直观理解，并允许为贝叶斯统计描述提出建议，这对科学报告的标准化至关重要。
#### 介绍
贝叶斯框架在心理学家和神经科学家中迅速普及。更喜欢这种方法的原因是可靠性，噪声数据的准确度更高，小样本的估计更好，I类错误的倾向性更小，将先验知识引入分析的可能性，以及关键的结果直观性和直接的解释。频率论的方法与对零假设检验的关注有关，并且p值的误用已被证明对心理科学的再现性危机起着关键作用。人们普遍认为，贝叶斯方法的推广是一种克服这一问题的方法 。
虽然贝叶斯框架推动的概率推理遍及大多数数据科学方面，但它已经为统计建模奠定了基础。心理学大量依赖的这个方面可以大致分为两个软边类别;预测和结构建模。虽然统计模型可以（通常）同时满足这两个目的，但预测建模致力于构建并找到准确预测给定结果的最佳模型。 它以适合度量，预测准确性和模型比较等概念为中心。在这个维度的极端情况下，机器和深度学习模型被用于强大的预测能力，通常以牺牲人类可读性为代价（这些模型通常被称为“黑盒子”，强调评估其内部功能的难度;另一方面，心理学家经常使用更简单的模型（例如与一般线性框架相关）来探索他们的数据。在此框架内，目标从构建最佳模型转变为理解模型内的参数。实际上，方法论管道通常从涉及模型比较的预测建模开始“什么是世界上最好的模型（即观察到的变量）”）后看似过渡到结构建模：“给定这个世界模型，如何效果（即模型参数）正在影响结果“。 对于最后一部分，他们通常依赖于效果指数“存在”。
事实上，虽然贝叶斯框架的优势之一是其概率参数估计，允许量化与每个估计相关的固有不确定性，但心理学家也对参数存在感兴趣：一个决策标准，允许他们得出结论，如果效果是“不同的从0“（统计上对应于”不可忽略的“或”不是相反的方向“）。换句话说，在了解效果的重要性，相关性或强度之前，要知道它是否与给定方向的结果有关。这种需要导致了频率主义p值的广泛采用，用作效应存在的指标，并且它的接受伴随着为其分类创建任意聚类（.05，.01和.001）。不幸的是，这些启发式方法严重僵化，成为达到目标和门槛而不是理解数据的工具 。
因此，贝叶斯框架在不需要这种零假设检验指标的情况下回答心理问题的能力通常被宣传为“没有p值的新世界”的承诺，而不是旧的和有缺陷的频率问题。尽管如此，似乎“效应存在”指数和标准对于人类获得对其数据的相互作用和结构的直观理解是有用的。因此，贝叶斯用户友好型实施的发展伴随着贝叶斯因子（BF）的推广，这一指数反映了模型对另一个模型的预测性能 。它提供了超过p值的许多优点，具有直接的解释（“数据是在替代下比零假设更容易发生3次（ BF = 3）”并且允许对替代方案做出陈述，而不仅仅是零假设。此外，最近的数学发展允许其计算复杂模型。 尽管与p值相比，BF符合预期的稳定，有效，直观和更好的指数，但其用于模型选择仍然是一个争论的问题。实际上，由于用于其计算的预测是从模型参数的先前分布生成的，因此它高度依赖于先验规范。重要的是，对于本文的目的，其用于估计较大模型中参数的效果存在的用途仍然是有限的，其计算在技术上是困难的并且其解释不像诸如t检验或相关性的“简单”测试那样简单。然而，考虑到对这些信息的需求，研究人员基于后验分布的特征开发了其他指数，其表示给定观察数据的不同参数值的概率分布。例如，通过呈现中心点（平均值，中位数，......）和离散度（标准偏差，中位数绝对偏差，......）的点估计值，可以总结这种不确定性，通常伴随一个百分比（89％，90％或95）最高密度区间的百分比（HDI;称为可信区间 - CI）。尽管贝叶斯框架提供了计算许多效应存在指数的可能性，但尚未对使用的存在指标达成共识，因为从未进行过任何比较。对于有兴趣采用贝叶斯框架的科学家来说，这可能是一个反驳。 此外，这个灰色区域会增加不熟悉贝叶斯框架的读者或评论者遵循假设和结论的难度，这可能反过来对整个研究产生不必要的怀疑。虽然我们认为这些指数及其解释指南（以经验法则的形式）在实践中是有用的，但我们也坚信这些指数应该伴随着与样本大小和效应大小相关的“行为”的知识。 。 这些知识对于人们隐含和直观地评估他们报告的数学价值的含义和含义非常重要。反过来，这可以防止可能的启发式和从这些指数得到的类别的结晶。
因此，基于多元线性回归（最广泛使用的模型之一）的模拟，本工作旨在比较仅从后验分布得出的几个效应存在指数，提供这些指数的“行为”的视觉表示。与样本大小，噪声，先验的关系以及频率论p值（超出其众多缺陷的指数，众所周知，可用作贝叶斯新手的参考），并为贝叶斯统计报告提出建议。
对于所有模拟模型，我们计算了以下指数：
* p值 ：基于频率回归，该指数表示给定统计模型的概率，当零假设为真时，统计汇总（例如两个比较组之间的样本均值差异）将大于或等于实际观察结果（Wasserstein，Lazar等人，2016年） 。
* “ 方向概率”（pd）：对应于效果严重为正或负的概率（以百分比表示）（与中位数符号一致）。
* p -MAP ：基于MAP的p值与参数对零假设的几率相关 。
* p -ROPE ：基于ROPE的p值，表示在ROPE定义的可忽略值空间中不包含（正值）或完全包含（负值）的HDI的最大百分比。
* ROPE（89％） ：该指数指的是ROPE中89％HDI的百分比，该区间被认为比基于频率的95％更稳定 。
* ROPE（95％） ：该指数是指ROPE中95％HDI的百分比。该间隔对应于Kruschke（2014）最初建议的间隔 。
* ROPE（完整） ：该指数是指ROPE内整个后验分布的百分比。
* 贝叶斯因子 ：TODO。
实用等效区域（ROPE）定义为-0.1至0.1。
### 研究1：与噪音和样本量的关系
#### 方法
模拟旨在调整以下特征：
* 模型类型 ：线性或逻辑。
* “真实”效果 （绘制数据的原始回归系数）：可以是1或0（无效果）。
* 样本量 ：从20到100步长为10。
* 误差 ：应用于预测器的高斯噪声，SD均匀分布在0.33和6.66之间（具有1000个不同的值）。
我们为这些特征的每个组合生成了一个数据集，从而产生了总共2*2*9*1000=36000贝叶斯和频率模型。 这里可以使用用于生成的代码（请注意，通常需要几天/几周才能完成）。
```{r}
df <- read.csv("https://raw.github.com/easystats/circus/master/data/bayesSim_study1.csv")
```
#### 结果
效果检测
对噪音的敏感度
```{r}
df %>%
select(outcome_type, true_effect, error, sample_size, p_value, p_direction, p_MAP, p_ROPE, ROPE_89, ROPE_95, ROPE_full, BF_log) %>%
gather(index, value, -error, -sample_size, -true_effect, -outcome_type) %>%
mutate(true_effect = as.factor(true_effect),
index = factor(index, levels=c("p_value", "p_direction", "p_MAP", "p_ROPE", "ROPE_89", "ROPE_95", "ROPE_full", "BF_log"))) %>%
mutate(temp = as.factor(cut(error, 10, labels = FALSE))) %>%
group_by(temp) %>%
mutate(error_group = round(mean(error), 1)) %>%
ungroup() %>%
ggplot(aes(x = error_group, y = value, fill = index, colour=true_effect, group = interaction(index, true_effect, error_group))) +
# geom_jitter(shape=16, alpha=0.02) +
geom_boxplot(outlier.shape = NA) +
facet_wrap(~index * outcome_type, scales = "free", ncol=2) +
theme_modern() +
scale_color_manual(values = c(`0` = "#f44336", `1` = "#8BC34A", name="Effect")) +
scale_fill_manual(values = c("p_value"="#607D8B", "p_MAP" = "#4CAF50", "p_direction" = "#2196F3",
"ROPE_89" = "#FFC107", "ROPE_95" = "#FF9800", "ROPE_full" = "#FF5722",
"p_ROPE"="#E91E63", "BF_log"="#9C27B0"), guide=FALSE) +
ylab("Index Value") +
xlab("Noise")
```
对样本量的敏感度
```{r}
df %>%
select(outcome_type, true_effect, error, sample_size, p_value, p_direction, p_MAP, p_ROPE, ROPE_89, ROPE_95, ROPE_full, BF_log) %>%
gather(index, value, -error, -sample_size, -true_effect, -outcome_type) %>%
mutate(true_effect = as.factor(true_effect),
index = factor(index, levels=c("p_value", "p_direction", "p_MAP", "p_ROPE", "ROPE_89", "ROPE_95", "ROPE_full", "BF_log"))) %>%
mutate(temp = as.factor(cut(sample_size, 10, labels = FALSE))) %>%
group_by(temp) %>%
mutate(size_group = round(mean(sample_size))) %>%
ungroup() %>%
ggplot(aes(x = size_group, y = value, fill = index, colour=true_effect, group = interaction(index, true_effect, size_group))) +
# geom_jitter(shape=16, alpha=0.02) +
geom_boxplot(outlier.shape = NA) +
facet_wrap(~index * outcome_type, scales = "free", ncol=2) +
theme_modern() +
scale_color_manual(values = c(`0` = "#f44336", `1` = "#8BC34A", name="Effect")) +
scale_fill_manual(values = c("p_value"="#607D8B", "p_MAP" = "#4CAF50", "p_direction" = "#2196F3",
"ROPE_89" = "#FFC107", "ROPE_95" = "#FF9800", "ROPE_full" = "#FF5722",
"p_ROPE"="#E91E63", "BF_log"="#9C27B0"), guide=FALSE) +
ylab("Index Value") +
xlab("Sample Size")
```
统计建模
我们拟合（频率）多元线性回归来预测不同指数的影响的存在与否，以及它们与噪声和样本大小的相互作用。
与频率p值的关系
```{r}
df %>%
select(outcome_type, true_effect, error, sample_size, p_value, p_direction, p_MAP, p_ROPE, ROPE_89, ROPE_95, ROPE_full, BF_log) %>%
gather(index, value, -error, -sample_size, -true_effect, -outcome_type, -p_value) %>%
mutate(true_effect = as.factor(true_effect),
index = factor(index, levels=c("p_direction", "p_MAP", "p_ROPE", "ROPE_89", "ROPE_95", "ROPE_full", "BF_log"))) %>%
mutate(temp = as.factor(cut(sample_size, 3, labels = FALSE))) %>%
group_by(temp) %>%
mutate(size_group = as.character(round(mean(sample_size)))) %>%
ungroup() %>%
ggplot(aes(x = p_value, y = value, color = true_effect, shape=size_group)) +
geom_point(alpha=0.025, stroke = 0, shape=16) +
facet_wrap(~index * outcome_type, scales = "free", ncol=2) +
theme_modern() +
scale_color_manual(values = c(`0` = "#f44336", `1` = "#8BC34A"), name="Effect") +
guides(colour = guide_legend(override.aes = list(alpha = 1)),
shape = guide_legend(override.aes = list(alpha = 1), title="Sample Size"))
```
与频率的基于p的任意聚类的关系
```{r}
df$sig_1 <- factor(ifelse(df$p_value >= .1, "n.s.", "-"), levels=c("n.s.", "-"))
df$sig_05 <- factor(ifelse(df$p_value >= .05, "n.s.", "*"), levels=c("n.s.", "*"))
df$sig_01 <- factor(ifelse(df$p_value >= .01, "n.s.", "**"), levels=c("n.s.", "**"))
df$sig_001 <- factor(ifelse(df$p_value >= .001, "n.s.", "***"), levels=c("n.s.", "***"))
get_data <- function(predictor, outcome, lbound=0, ubound=0.3){
fit <- suppressWarnings(glm(paste(outcome, "~ outcome_type * ", predictor), data=df, family = "binomial"))
# data <- data.frame(x=rep(1:100, 2))
data <- data.frame(outcome_type=rep(c("linear", "binary"), each=100))
data[predictor] <- rep(seq(lbound, ubound, length.out = 100), 2)
data$index <- predictor
predict_fit <- predict(fit, newdata=data, type="response", se.fit = TRUE)
data[outcome] <- predict_fit$fit
data$CI_lower <- predict_fit$fit - (qnorm(0.99) * predict_fit$se.fit)
data$CI_upper <- predict_fit$fit + (qnorm(0.99) * predict_fit$se.fit)
data <-
select(
data,
"value" := !!predictor,
!!outcome,
.data$outcome_type,
.data$index,
.data$CI_lower,
.data$CI_upper
)
return(data)
}
rbind(
rbind(
get_data(predictor="p_direction", outcome="sig_001", lbound=99.5, ubound=100),
get_data(predictor="p_MAP", outcome="sig_001", lbound=0, ubound=0.01),
get_data(predictor="p_ROPE", outcome="sig_001", lbound=97, ubound=100),
get_data(predictor="ROPE_89", outcome="sig_001", lbound=0, ubound=0.5),
get_data(predictor="ROPE_95", outcome="sig_001", lbound=0, ubound=0.5),
get_data(predictor="ROPE_full", outcome="sig_001", lbound=0, ubound=0.5),
get_data(predictor="BF_log", outcome="sig_001", lbound=0, ubound=10)
) %>%
rename("sig"=sig_001) %>%
mutate(threshold="p < .001"),
rbind(
get_data(predictor="p_direction", outcome="sig_01", lbound=98, ubound=100),
get_data(predictor="p_MAP", outcome="sig_01", lbound=0, ubound=0.1),
get_data(predictor="p_ROPE", outcome="sig_01", lbound=85, ubound=100),
get_data(predictor="ROPE_89", outcome="sig_01", lbound=0, ubound=2),
get_data(predictor="ROPE_95", outcome="sig_01", lbound=0, ubound=2),
get_data(predictor="ROPE_full", outcome="sig_01", lbound=0, ubound=2),
get_data(predictor="BF_log", outcome="sig_01", lbound=0, ubound=5)
) %>%
rename("sig"=sig_01) %>%
mutate(threshold="p < .01"),
rbind(
get_data(predictor="p_direction", outcome="sig_05", lbound=95, ubound=100),
get_data(predictor="p_MAP", outcome="sig_05", lbound=0, ubound=0.3),
get_data(predictor="p_ROPE", outcome="sig_05", lbound=50, ubound=100),
get_data(predictor="ROPE_89", outcome="sig_05", lbound=0, ubound=10),
get_data(predictor="ROPE_95", outcome="sig_05", lbound=0, ubound=10),
get_data(predictor="ROPE_full", outcome="sig_05", lbound=0, ubound=10),
get_data(predictor="BF_log", outcome="sig_05", lbound=0, ubound=2)
) %>%
rename("sig"=sig_05) %>%
mutate(threshold="p < .05"),
rbind(
get_data(predictor="p_direction", outcome="sig_1", lbound=90, ubound=100),
get_data(predictor="p_MAP", outcome="sig_1", lbound=0, ubound=0.5),
get_data(predictor="p_ROPE", outcome="sig_1", lbound=25, ubound=100),
get_data(predictor="ROPE_89", outcome="sig_1", lbound=0, ubound=20),
get_data(predictor="ROPE_95", outcome="sig_1", lbound=0, ubound=20),
get_data(predictor="ROPE_full", outcome="sig_1", lbound=0, ubound=20),
get_data(predictor="BF_log", outcome="sig_1", lbound=0, ubound=1)
) %>%
rename("sig"=sig_1) %>%
mutate(threshold="p < .1")
) %>%
mutate(index = as.factor(index)) %>%
ggplot(aes(x=value, y=sig)) +
geom_ribbon(aes(ymin=CI_lower, ymax=CI_upper), alpha=0.1) +
geom_line(aes(color=index), size=1) +
facet_wrap(~ index * threshold * outcome_type, scales = "free", ncol=8) +
theme_modern() +
scale_color_manual(values = c("p_value"="#607D8B", "p_MAP" = "#4CAF50", "p_direction" = "#2196F3",
"ROPE_89" = "#FFC107", "ROPE_95" = "#FF9800", "ROPE_full" = "#FF5722",
"p_ROPE"="#E91E63", "BF_log"="#9C27B0"), guide=FALSE) +
ylab("Probability of being significant") +
xlab("Index Value")
```
与等效性测试的关系
```{r}
df$equivalence_95 <- factor(ifelse(df$ROPE_95 == 0, "significant", "n.s."), levels=c("n.s.", "significant"))
df$equivalence_89 <- factor(ifelse(df$ROPE_89 == 0, "significant", "n.s."), levels=c("n.s.", "significant"))
rbind(
rbind(
get_data(predictor="p_direction", outcome="equivalence_95", lbound=97.5, ubound=100),
get_data(predictor="p_MAP", outcome="equivalence_95", lbound=0, ubound=0.2),
get_data(predictor="p_ROPE", outcome="equivalence_95", lbound=92.5, ubound=97.5),
get_data(predictor="ROPE_89", outcome="equivalence_95", lbound=0, ubound=3),
get_data(predictor="ROPE_full", outcome="equivalence_95", lbound=0, ubound=3),
get_data(predictor="BF_log", outcome="equivalence_95", lbound=0.5, ubound=2.5)
) %>%
rename("equivalence"=equivalence_95) %>%
mutate(level="95 HDI"),
rbind(
get_data(predictor="p_direction", outcome="equivalence_89", lbound=99.5, ubound=100),
get_data(predictor="p_MAP", outcome="equivalence_89", lbound=0, ubound=0.02),
get_data(predictor="p_ROPE", outcome="equivalence_89", lbound=0, ubound=100),
get_data(predictor="ROPE_95", outcome="equivalence_89", lbound=0, ubound=0.5),
get_data(predictor="ROPE_full", outcome="equivalence_89", lbound=0, ubound=0.5),
get_data(predictor="BF_log", outcome="equivalence_89", lbound=0, ubound=3)
) %>%
rename("equivalence"=equivalence_89) %>%
mutate(level="90 HDI")
) %>%
ggplot(aes(x=value, y=equivalence)) +
geom_ribbon(aes(ymin=CI_lower, ymax=CI_upper), alpha=0.1) +
geom_line(aes(color=index), size=1) +
facet_wrap(~ index * level * outcome_type, scales = "free", nrow=7) +
theme_modern() +
scale_color_manual(values = c("p_value"="#607D8B", "p_MAP" = "#4CAF50", "p_direction" = "#2196F3",
"ROPE_89" = "#FFC107", "ROPE_95" = "#FF9800", "ROPE_full" = "#FF5722",
"p_ROPE"="#E91E63", "BF_log"="#9C27B0"), guide=FALSE) +
ylab("Probability of rejecting H0 with the equivalence test") +
xlab("Index Value")
```
ROPE（完整）与p（方向）之间的关系
```{r}
df %>%
mutate(true_effect = as.factor(true_effect)) %>%
ggplot(aes(x=p_direction, y=ROPE_full, color=true_effect)) +
geom_point(alpha=0.025, stroke = 0, shape=16) +
facet_wrap(~ outcome_type, scales = "free", ncol=2) +
theme_modern() +
scale_color_manual(values = c(`0` = "#f44336", `1` = "#8BC34A"), name="Effect") +
ylab("ROPE (full)") +
xlab("Probability of Direction (pd)")
```
### 研究2：与采样特征的关系
### 研究3：与Priors规范的关系
讨论
基于多元线性回归的模拟，目前的工作旨在比较几个仅来自后验分布的效应存在指数，提供这些指数的“行为”与样本大小，噪声，先验和频率论的关系的视觉表示。 p值。
虽然这种与频率指数的比较可能看似违反直觉或错误（因为贝叶斯思想与频率论框架本质上不同），但我们认为这种比较因教学原因而引人注目。频率的p值对许多人“说话”，因此可以被视为参考和促进向贝叶斯框架转变的一种方式。然而，这并不排除在必要时寻求效应存在的一般范式的变化，并且贝叶斯指数与频率论p基本不同，而不仅仅是近似或等价。 至关重要的是，我们非常同意效果的存在与意义之间的区别和可能的分离。尽管如此，我们认为评估影响是否“有意义”在很大程度上依赖于文献，先验，新颖性，背景或领域，并且不能仅仅基于统计指标进行评估（即使某些指数，如ROPE相关的，尝试以有意义的方式弥合存在）。因此，研究人员应该依靠统计数据来评估效应的存在（以及大小和方向估计），并且系统地，但在上下文中，从更大的角度讨论其意义和重要性。
结论可以在指南部分找到。
## 报告指南
### 如何描述和报告模型的参数
贝叶斯分析返回每个参数（或效果）的后验分布。为了最低限度地描述这些分布，我们建议报告中心点的点估计以及表征估计不确定性（分散）的信息。另外，还可以报告效果存在和/或重要性的指数。基于先前对点估计和效应存在指数的比较 ，我们可以得出以下建议。
#### 确定性
我们建议将中位数报告为中心性指数，因为与平均值或MAP估计值相比，它更加稳健。然而，在严重偏斜的后部分布的情况下，MAP估计可能是一个很好的选择。
#### 不确定性
89％可信区间（CI）似乎是表征与估计相关的不确定性的合理范围，比更高的阈值（例如90％和95％）更稳定。 我们还建议基于HDI而不是分位数计算CI，其倾向于可能的值而不是中心值。注意，基于分位数（等尾区间）的CI在转换的情况下可能更合适（例如，当将对数转换为概率时）。 否则，最初未覆盖null的间隔可能会在转换后覆盖它。
#### 存在
贝叶斯框架允许整齐地描绘和量化假设检验的不同方面，例如效应存在和重要性。描述效应存在的最直接的指标是方向概率（pd），表示与效果的最可能方向（正面或负面）相关联的确定性。 该索引易于理解，易于解释，易于计算，对模型特征具有鲁棒性，并且与数据规模无关。
此外，它与频率论p值密切相关，因此可用于绘制相似之处，并为不熟悉贝叶斯统计的读者提供参考。 分别为.1， .05，.01和.001的双侧p值大致对应于95％，97.5％，99.5％和99.95％的pd 。 因此，为方便起见，我们建议将以下参考值作为解释助手：
* pd <= 95％ ~ p > .1：不确定
* pd > 95％ ~ p <.1：有可能存在
* pd > 97％ ：可能存在
* pd > 99％ ：极可能存在
* pd > 99.9％：肯定存在
#### 意义
ROPE中的百分比是显著性指数（在其主要含义中），告知我们参数是否与结果中的不可忽略的变化（就数量而言）相关 - 或不相关。 我们建议在ROPE中报告完全后验分布的百分比（完整的 ROPE）而不是给定比例的CI，这似乎更敏感（特别是描绘高度显著的影响）。我们建议使用百分比作为连续的显著性指数，而不是将其用作二进制，全有或全无的决策标准，例如原始等效性测试所建议的。 但是，根据模拟数据 ，我们建议将以下参考值作为解释助手：
* ROPE > 99％ ：可以忽略不计（我们可以接受零假设）
* ROPE > 97.5％ ：可能微不足道
* ROPE<= 97.5％ ＆ > = 2.5％ ：未定意义
* ROPE <2.5％ ：可能很重要
* ROPE <1％ ：显著（我们可以拒绝零假设）
请注意，需要格外小心，因为它的解释很大程度上取决于其他参数，如样本大小和ROPE范围（见此处 ） 。
#### 模板
基于这些建议，基于其后验分布的参数的最小报告的模板句子可以是：“X的影响具有pd的概率为负,中位数= ，89％CI[HDI低，HDI高]并且可以被认为是显著的（ROPE% in ROPE）。”
### 如何比较不同的模型
尽管它也可用于评估效应的存在和重要性，贝叶斯因子（BF）是一种通用的指数，可用于直接比较不同的模型（或数据生成过程）。贝叶斯因子是一个比率，通知我们观察到的数据在两个比较模型下的可能性有多大（或更少）-通常是具有效果的模型与没有效果的模型。根据空模型的规范（无论是点估计（例如，0）还是间隔），贝叶斯因子可以在效应存在和重要性的上下文中使用。
一般来说，贝叶斯因子大于1，证明有利于其中一个模型，而贝叶斯因子小于1，证明有利于另一个模型。 有几条经验法则可以帮助解释（见这里），>3是一个常用的门槛，可以对非轶事证据进行分类。
#### 模板
报告贝叶斯因子（BF）时，可以使用以下句子：“有中等证据支持x（BF = BF） 不存在影响。”
评论
发表评论