贝叶斯分析Example
---
title: "Untitled"
author: "xuefliang"
date: "7/11/2019"
output:
html_document:
theme: readable
highlight: haddock
df_print: paged
code_folding: show
toc: true
number_sections: true
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)
```
## 示例1：启动贝叶斯模型
一旦安装了必要的软件包，我们就可以加载rstanarm（拟合模型），bayestestR（计算有用的指标）和访问参数。
### 简单的线性模型（又名回归）
我们将从一个简单的回归开始，测试Petal.Length（预测变量或独立变量）与Sepal.Length（响应变量或依赖变量）之间的关系，它们来自iris数据集，默认包含在R中。
#### 拟合模型
让我们从适合模型的频率版本开始，只是为了得到一些参考：
```{r}
model <- lm(Sepal.Length ~ Petal.Length, data=iris)
summary(model)
```
在该模型中，Petal.Length和Sepal.Length之间的线性关系是正的和显着的（β= 0.41，t= 21.6，p <.001）。这意味着，对于Petal.Length（预测值）中每增加1，Sepal.Length（响应）增加0.41。 这可以通过使用ggplot2包绘制x轴上的预测值和响应值y来显示：
```{r}
## The ggplot function takes the data as argument, and then the variables
## related to aesthetic features such as the x and y axes.
ggplot(iris, aes(x=Petal.Length, y=Sepal.Length)) +
geom_point() + ## This adds the points
geom_smooth(method="lm") ## This adds a regression line
```
让我们拟合该模型的贝叶斯版本：
```{r}
model <- stan_glm(Sepal.Length ~ Petal.Length, data=iris)
```
您可以看到正在运行的采样算法。
### 后验概率
完成后，让我们提取模型的参数（即系数）。
```{r}
posteriors <- insight::get_parameters(model)
head(posteriors) ## Show the first 6 rows
```
我们可以看到，参数采用具有两列的长数据的形式，对应于截距和Petal.Length的效果。 这些列包含这两个参数的后验分布，即一组不同的合理值。
### 关于 posterior draws
让我们来看看后验的长度。
```{r}
nrow(posteriors) ## Size (number of rows)
```
#### 为什么大小为4000，而不是更多或更少？
首先，这些观察值（行）通常被称为posterior draws。基本思想是贝叶斯采样算法（如蒙特卡罗马尔可夫链-MCMC）将从隐藏的真实后验分布中得出。因此，通过这些posterior draws，我们可以估计潜在的真实后验分布。你得到的posterior draws越多，你对后验分布的估计就越好。但是计算也需要更长的时间。
如果我们查看上面模型中默认使用的rstanarm “sampling”算法的文档，我们可以看到几个影响posterior draws次数的参数。默认情况下，有4个链（您可以将其视为不同的采样运行），每个链创建2000 iter（draws）。但是，这些迭代中只有一半被保留，因为一半用于预热（算法的收敛）。总数是4链*（2000次迭代-1000次预热）=4000次后抽取。我们可以改变它，例如：
```{r}
model <- stan_glm(Sepal.Length ~ Petal.Length, data=iris, chains = 2, iter = 1000, warmup = 250)
nrow(insight::get_parameters(model)) ## Size (number of rows)
```
在这种情况下，我们确实有2个链*（1000次迭代-250次预热）= 1500次posterior draws。让我们使用默认设置保留我们的第一个模型。
### 可视化后验分布
现在我们已经了解了这些值来自何处，我们将首先可视化我们感兴趣的参数的分布，Petal.Length的影响。
```{r}
ggplot(posteriors, aes(x = Petal.Length)) +
geom_density(fill = "orange")
```
该分布表示不同效果（x轴）的概率（y轴）。中心值比极值更可能。正如您所看到的，这种分布范围从大约0.35到0.50，其中大部分在0.41左右。恭喜！你刚刚描述了你的后验分布。这是贝叶斯分析的核心。 我们不需要p值，t值或自由度，一切都在在这个后验分布中。如您所见，我们的描述与频率派回归的值（其值为0.41）一致。这让人放心！事实上，在大多数情况下，贝叶斯分析并没有大幅改变结果或解释。相反，它使结果更具解释性，直观性，易于理解和描述。我们现在可以继续并精确地描述这种后验分布。
### 描述后验
不幸的是，将整个后验分布报告为图表通常是不切实际的。我们需要找到一种简明的方法来总结它。 我们建议用3个元素描述后验分布。
1、点估计，单值摘要（类似于频率论回归中的beta）。 2、一个可靠的区间，代表相关的不确定性。 3、一些重要指数，提供有关这种影响重要性的信息。
#### 点估计
什么单一值最能代表我的后验分布？
中心化指标，例如均值，中位数或模式通常用作点估计。 有什么不同？ 我们来检查一下这个意思：
```{r}
mean(posteriors$Petal.Length)
```
这接近频率论的版本。 但正如我们所知，平均值对异常值或极值值非常敏感。也许中位数可能更强劲？
```{r}
median(posteriors$Petal.Length)
```
这接近平均值。 也许我们可以采取后验分布的高峰？在贝叶斯框架中，该值称为最大后验（MAP）。
```{r}
map_estimate(posteriors$Petal.Length)
```
他们都非常接近！ 让我们在后验分布上可视化这些值：
```{r}
ggplot(posteriors, aes(x = Petal.Length)) +
geom_density(fill = "orange") +
## The mean in blue
geom_vline(xintercept=mean(posteriors$Petal.Length), color="blue", size=1) +
## The median in red
geom_vline(xintercept=median(posteriors$Petal.Length), color="red", size=1) +
## The MAP in purple
geom_vline(xintercept=map_estimate(posteriors$Petal.Length), color="purple", size=1)
```
所有这些值都给出了非常相似的结果。因此，我们将选择中位数，因为这个值从概率的角度来看具有直接意义：有50％的可能性比真实效果低或高（因为它在两个相等的部分削减分布）。
#### 不确定性
现在有了点估计，我们必须描述不确定性。 我们可以计算范围：
```{r}
range(posteriors$Petal.Length)
```
但是包含所有这些极端值是否有意义？ 可能不是。 因此，我们将计算可靠的间隔。 长话短说，它类似于频率派的置信区间，但它更容易解释，更容易计算,它更有意义。
我们将根据最高密度区间（HDI）计算此可靠区间。它将为我们提供包含89％更可能的效果值的范围。 请注意，我们将使用89％CI而不是95％CI（如在频率框架中），因为89％的级别提供更稳定的结果（Kruschke 2014）并且提醒我们这些约定的任意性（McElreath 2018）。
```{r}
hdi(posteriors$Petal.Length, ci=0.89)
```
很好，所以我们可以得出结论，该效应有89％的可能性落在[0.38,0.44]范围内。我们刚刚计算了两个最重要的信息来描述我们的效果。
#### 效果意义
然而，在许多科学领域，仅仅描述效果是不够的。 读者也想知道这种效果是否意味着什么。 那很重要么？ 它与0不同吗？ 换句话说，评估效果的重要性。 我们应该怎么做？
嗯，在这种特殊情况下，它是非常雄辩的：所有可能的效果值（即整个后验分布）是正的并且优于0.35，这已经很多了。
但是，我们仍然需要一些客观的决策标准，如果是或否，效果是显着的。 我们可以类似于频率框架，并查看可信区间是否包含0这里，情况并非如此，这意味着我们的影响很大。
但这个指数不是很细，不是吗？ 我们可以做得更好吗？ 是。
### 具有分类预测器的线性模型
让我们对鸡的重量感兴趣，这取决于2种饲料类型。 我们将首先从chickwts数据集（也可在基础R中获得）中选择我们感兴趣的2种饲料类型（因为原因）：肉粉和向日葵
#### 数据准备和模型拟合
```{r}
## We keep only rows for which feed is meatmeal or sunflower
data <- chickwts %>%
filter(feed %in% c("meatmeal", "sunflower"))
```
让我们运行另一个贝叶斯回归来预测两种饲料类型的权重。
```{r}
model <- stan_glm(weight ~ feed, data=data)
```
#### 后验描述
```{r}
posteriors <- insight::get_parameters(model)
ggplot(posteriors, aes(x=feedsunflower)) +
geom_density(fill = "red")
```
这代表了肉粉和向日葵之间差异的后验分布。似乎差异是相当积极的...吃向日葵会让你更胖（如果你是一只鸡）。 多少钱？ 让我们计算中位数和CI：
```{r}
median(posteriors$feedsunflower)
hdi(posteriors$feedsunflower)
```
它使你的脂肪达到约51（中位数）。 然而，不确定性非常高：两种饲料类型之间的差异有89％的可能性介于7.77和87.66之间。
这个效果与0有什么不同？
#### ROPE百分比
测试此分布是否与0不同是没有意义的，因为0是单个值（并且任何分布与单个值不同的概率是无限的）。然而，评估重要性的一种方法可能是定义0左右的区域，这将被视为等于零（不，-或可忽略不计，效果）。这被称为实用等价区（ROPE），是测试参数重要性的一种方法。
我们如何定义这个区域？现在我们知道我们可以将ROPE定义为[-20,20]范围。此范围内的所有效果均视为无效（可忽略不计）。我们现在可以计算89％最可能值（89％CI）的比例，这些值不为空，即超出此范围。
```{r}
rope(posteriors$feedsunflower, range = c(-20, 20), ci=0.89)
```
89％CI中的7.75％可被视为无效。 那是很多吗？ 根据我们的指导方针，是的，这太过分了。 基于ROPE的这一特定定义，我们得出结论，这种影响并不显着（可忽略的概率太高）。
虽然，老实说，我对霍华德教授有些怀疑。 我真的不相信他对ROPE的定义。 有更客观的方法来定义它吗？
是。 其中一种做法是例如使用响应变量的标准偏差（SD）的十分之一（1/10 = 0.1），这可以被认为是“可忽略的”效应大小。
```{r}
rope_value <- 0.1 * sd(data$weight)
rope_range <- c(-rope_value, rope_value)
rope_range
```
让我们重新定义我们的ROPE作为[-6.2,6.2]范围内的区域。注意，这可以通过rope_range函数直接获得:
```{r}
rope_value <- rope_range(model)
rope_range
```
让我们重新计算ROPE中的百分比：
```{r}
rope(posteriors$feedsunflower, range = rope_range, ci=0.89)
```
根据ROPE的这种合理定义，我们观察到效应的89％的后验分布与ROPE不重叠。因此，我们可以得出结论，这种影响是显着的（在重要性的意义上要注意）。
### Probability of Direction方向概率（pd）
也许我们对这种影响是否可以忽略不感兴趣。也许我们只想知道这种影响是积极的还是消极的。 在这种情况下，无论效果的“大小”如何，我们都可以简单地计算后验的正比例。
```{r}
n_positive <- posteriors %>%
filter(feedsunflower > 0) %>% ## select only positive values
nrow() ## Get length
n_positive / nrow(posteriors) * 100
```
我们可以得出结论，效果是积极的，概率为97.82％。 我们称这个指数为方向概率（pd）。 事实上，它可以通过以下方式更容易地计算：
```{r}
p_direction(posteriors$feedsunflower)
```
有趣的是，这个指数通常与频率论的p值高度相关。我们几乎可以通过简单的转换大致推断出相应的p值：
```{r}
pd <- 97.82
onesided_p <- 1 - pd / 100
twosided_p <- onesided_p * 2
twosided_p
```
如果我们在频率框架中运行我们的模型，我们应该近似观察p值为0.04的效果。 真的吗？
### 与频率论者的比较
```{r}
lm(weight ~ feed, data=data) %>%
summary()
```
频率论模型告诉我们差异是积极的和显着的（$\beta$= 52，$p$ = 0.04）。
虽然我们得出了类似的结论，但贝叶斯框架使我们能够对我们的影响及其估计的不确定性有更深刻和直观的理解。
### 一个实现所有需求的功能
然而，我同意，提取和计算所有指数有点乏味。 但是如果我告诉你我们只用一个功能就可以完成所有这些工作呢？
此函数计算所有提到的索引，并且可以直接在模型上运行：
```{r}
describe_posterior(model)
```
我们有它！ 中位数，CI，pd和ROPE百分比！
理解和描述后验分布只是贝叶斯建模的一个方面......你准备好了吗？
## 示例2：确认贝叶斯技能
现在描述和理解线性回归的后验分布对你来说没有任何秘密，让我们回过头来研究一些更简单的模型： 相关性和t-测试 。
但在我们这样做之前，让我们花点时间提醒自己，并了解所有基本的统计程序 ，如相关性， t-测试，ANOVA或Chisquare测试都是线性回归 （我们强烈推荐这个优秀的演示 ）。 但是，这些简单模型仍将是引入更复杂指数的机会，例如贝叶斯因子 。
### 相关性
#### 频率版本
两个连续变量之间的频率相关性 ，一些花的萼片的宽度和长度 。 数据在R中可用作iris数据集（与前一个教程中使用的相同）。让我们计算一个Pearson相关性测试，将结果存储在一个名为result的对象中，然后显示它：
```{r}
result <- cor.test (iris $ Sepal.Width, iris $ Sepal.Length)
result
```
正如您在输出中所看到的，我们所做的测试实际上比较了两个假设： 零假设 （无相关）与备选假设 （非零相关）。基于p值，不能拒绝零假设：两个变量之间的相关性是负的但不显着 （r = -12，p> .05）。
#### 贝叶斯相关
要计算贝叶斯相关性测试，我们需要BayesFactor包，使用correlationBF()函数计算相关correlationBF()并以类似的方式存储结果。
```{r}
library (BayesFactor)
result <- correlationBF (iris $ Sepal.Width, iris $ Sepal.Length)
#让我们运行describe_posterior()函数：
describe_posterior (result)
```
我们在这里再看到很多东西，但现在重要的指数是后验分布的中位数为-0.11 。这非常接近频率相关性。 如前所述，我们可以描述可信区间 ， pd或ROPE百分比 ，但我们将重点关注贝叶斯框架贝叶斯因子提供的另一个指数。
### 贝叶斯因子（BF）
我们说，相关性实际上比较了两个假设，一个无效（没有效果）和一个altnernative（存在效应）。 贝叶斯因子（BF）允许相同的比较并且确定观察数据更可能在两个模型中的哪个模型下 ：具有感兴趣效果的模型和没有感兴趣效果的空模型。 我们可以使用bayesfactor()来专门计算贝叶斯因子，比较这些模型
```{r}
bayesfactor (result)
```
我们得到了BF=0.51 。 这是什么意思？
贝叶斯因子是相对证据的连续测量，贝叶斯因子大于1，证明有利于其中一个模型（通常称为分子），贝叶斯因子小于1，证明有利于另一个模型（分母）。
是的，你听到了正确的事，支持null的证据！这就是贝叶斯框架有时被认为优于频率框架的原因之一。 请记住，从您的统计数据课程中， p值只能用于拒绝h0 ，但不能接受它。 使用贝叶斯因子,您可以衡量证据支持null假设的情况 。
表示替代null假设的证据的BF可以使用来反转，$BF_{01}=1/BF_{10}$以提供替代null的证据。 在针对null假设的替代的BF小于1（支持null）的情况下，这提高了人的可读性。
在我们的例子中， BF = 1/0.51 = 2 ，表明与替代假设相比，数据null假设下的概率是2倍 ，尽管有利于空，但只被认为是反对空的轶事证据。因此，我们可以得出结论，有传闻证据表明两个变量之间缺乏相关性（r中位数 = 0.11，BF = 0.51） ，这是一个更具信息性的陈述，表明我们可以用频率统计数据做些什么。
这还不是全部！
## 示例3：成为贝叶斯主人
未完待续。
title: "Untitled"
author: "xuefliang"
date: "7/11/2019"
output:
html_document:
theme: readable
highlight: haddock
df_print: paged
code_folding: show
toc: true
number_sections: true
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)
```
## 示例1：启动贝叶斯模型
一旦安装了必要的软件包，我们就可以加载rstanarm（拟合模型），bayestestR（计算有用的指标）和访问参数。
### 简单的线性模型（又名回归）
我们将从一个简单的回归开始，测试Petal.Length（预测变量或独立变量）与Sepal.Length（响应变量或依赖变量）之间的关系，它们来自iris数据集，默认包含在R中。
#### 拟合模型
让我们从适合模型的频率版本开始，只是为了得到一些参考：
```{r}
model <- lm(Sepal.Length ~ Petal.Length, data=iris)
summary(model)
```
在该模型中，Petal.Length和Sepal.Length之间的线性关系是正的和显着的（β= 0.41，t= 21.6，p <.001）。这意味着，对于Petal.Length（预测值）中每增加1，Sepal.Length（响应）增加0.41。 这可以通过使用ggplot2包绘制x轴上的预测值和响应值y来显示：
```{r}
## The ggplot function takes the data as argument, and then the variables
## related to aesthetic features such as the x and y axes.
ggplot(iris, aes(x=Petal.Length, y=Sepal.Length)) +
geom_point() + ## This adds the points
geom_smooth(method="lm") ## This adds a regression line
```
让我们拟合该模型的贝叶斯版本：
```{r}
model <- stan_glm(Sepal.Length ~ Petal.Length, data=iris)
```
您可以看到正在运行的采样算法。
### 后验概率
完成后，让我们提取模型的参数（即系数）。
```{r}
posteriors <- insight::get_parameters(model)
head(posteriors) ## Show the first 6 rows
```
我们可以看到，参数采用具有两列的长数据的形式，对应于截距和Petal.Length的效果。 这些列包含这两个参数的后验分布，即一组不同的合理值。
### 关于 posterior draws
让我们来看看后验的长度。
```{r}
nrow(posteriors) ## Size (number of rows)
```
#### 为什么大小为4000，而不是更多或更少？
首先，这些观察值（行）通常被称为posterior draws。基本思想是贝叶斯采样算法（如蒙特卡罗马尔可夫链-MCMC）将从隐藏的真实后验分布中得出。因此，通过这些posterior draws，我们可以估计潜在的真实后验分布。你得到的posterior draws越多，你对后验分布的估计就越好。但是计算也需要更长的时间。
如果我们查看上面模型中默认使用的rstanarm “sampling”算法的文档，我们可以看到几个影响posterior draws次数的参数。默认情况下，有4个链（您可以将其视为不同的采样运行），每个链创建2000 iter（draws）。但是，这些迭代中只有一半被保留，因为一半用于预热（算法的收敛）。总数是4链*（2000次迭代-1000次预热）=4000次后抽取。我们可以改变它，例如：
```{r}
model <- stan_glm(Sepal.Length ~ Petal.Length, data=iris, chains = 2, iter = 1000, warmup = 250)
nrow(insight::get_parameters(model)) ## Size (number of rows)
```
在这种情况下，我们确实有2个链*（1000次迭代-250次预热）= 1500次posterior draws。让我们使用默认设置保留我们的第一个模型。
### 可视化后验分布
现在我们已经了解了这些值来自何处，我们将首先可视化我们感兴趣的参数的分布，Petal.Length的影响。
```{r}
ggplot(posteriors, aes(x = Petal.Length)) +
geom_density(fill = "orange")
```
该分布表示不同效果（x轴）的概率（y轴）。中心值比极值更可能。正如您所看到的，这种分布范围从大约0.35到0.50，其中大部分在0.41左右。恭喜！你刚刚描述了你的后验分布。这是贝叶斯分析的核心。 我们不需要p值，t值或自由度，一切都在在这个后验分布中。如您所见，我们的描述与频率派回归的值（其值为0.41）一致。这让人放心！事实上，在大多数情况下，贝叶斯分析并没有大幅改变结果或解释。相反，它使结果更具解释性，直观性，易于理解和描述。我们现在可以继续并精确地描述这种后验分布。
### 描述后验
不幸的是，将整个后验分布报告为图表通常是不切实际的。我们需要找到一种简明的方法来总结它。 我们建议用3个元素描述后验分布。
1、点估计，单值摘要（类似于频率论回归中的beta）。 2、一个可靠的区间，代表相关的不确定性。 3、一些重要指数，提供有关这种影响重要性的信息。
#### 点估计
什么单一值最能代表我的后验分布？
中心化指标，例如均值，中位数或模式通常用作点估计。 有什么不同？ 我们来检查一下这个意思：
```{r}
mean(posteriors$Petal.Length)
```
这接近频率论的版本。 但正如我们所知，平均值对异常值或极值值非常敏感。也许中位数可能更强劲？
```{r}
median(posteriors$Petal.Length)
```
这接近平均值。 也许我们可以采取后验分布的高峰？在贝叶斯框架中，该值称为最大后验（MAP）。
```{r}
map_estimate(posteriors$Petal.Length)
```
他们都非常接近！ 让我们在后验分布上可视化这些值：
```{r}
ggplot(posteriors, aes(x = Petal.Length)) +
geom_density(fill = "orange") +
## The mean in blue
geom_vline(xintercept=mean(posteriors$Petal.Length), color="blue", size=1) +
## The median in red
geom_vline(xintercept=median(posteriors$Petal.Length), color="red", size=1) +
## The MAP in purple
geom_vline(xintercept=map_estimate(posteriors$Petal.Length), color="purple", size=1)
```
所有这些值都给出了非常相似的结果。因此，我们将选择中位数，因为这个值从概率的角度来看具有直接意义：有50％的可能性比真实效果低或高（因为它在两个相等的部分削减分布）。
#### 不确定性
现在有了点估计，我们必须描述不确定性。 我们可以计算范围：
```{r}
range(posteriors$Petal.Length)
```
但是包含所有这些极端值是否有意义？ 可能不是。 因此，我们将计算可靠的间隔。 长话短说，它类似于频率派的置信区间，但它更容易解释，更容易计算,它更有意义。
我们将根据最高密度区间（HDI）计算此可靠区间。它将为我们提供包含89％更可能的效果值的范围。 请注意，我们将使用89％CI而不是95％CI（如在频率框架中），因为89％的级别提供更稳定的结果（Kruschke 2014）并且提醒我们这些约定的任意性（McElreath 2018）。
```{r}
hdi(posteriors$Petal.Length, ci=0.89)
```
很好，所以我们可以得出结论，该效应有89％的可能性落在[0.38,0.44]范围内。我们刚刚计算了两个最重要的信息来描述我们的效果。
#### 效果意义
然而，在许多科学领域，仅仅描述效果是不够的。 读者也想知道这种效果是否意味着什么。 那很重要么？ 它与0不同吗？ 换句话说，评估效果的重要性。 我们应该怎么做？
嗯，在这种特殊情况下，它是非常雄辩的：所有可能的效果值（即整个后验分布）是正的并且优于0.35，这已经很多了。
但是，我们仍然需要一些客观的决策标准，如果是或否，效果是显着的。 我们可以类似于频率框架，并查看可信区间是否包含0这里，情况并非如此，这意味着我们的影响很大。
但这个指数不是很细，不是吗？ 我们可以做得更好吗？ 是。
### 具有分类预测器的线性模型
让我们对鸡的重量感兴趣，这取决于2种饲料类型。 我们将首先从chickwts数据集（也可在基础R中获得）中选择我们感兴趣的2种饲料类型（因为原因）：肉粉和向日葵
#### 数据准备和模型拟合
```{r}
## We keep only rows for which feed is meatmeal or sunflower
data <- chickwts %>%
filter(feed %in% c("meatmeal", "sunflower"))
```
让我们运行另一个贝叶斯回归来预测两种饲料类型的权重。
```{r}
model <- stan_glm(weight ~ feed, data=data)
```
#### 后验描述
```{r}
posteriors <- insight::get_parameters(model)
ggplot(posteriors, aes(x=feedsunflower)) +
geom_density(fill = "red")
```
这代表了肉粉和向日葵之间差异的后验分布。似乎差异是相当积极的...吃向日葵会让你更胖（如果你是一只鸡）。 多少钱？ 让我们计算中位数和CI：
```{r}
median(posteriors$feedsunflower)
hdi(posteriors$feedsunflower)
```
它使你的脂肪达到约51（中位数）。 然而，不确定性非常高：两种饲料类型之间的差异有89％的可能性介于7.77和87.66之间。
这个效果与0有什么不同？
#### ROPE百分比
测试此分布是否与0不同是没有意义的，因为0是单个值（并且任何分布与单个值不同的概率是无限的）。然而，评估重要性的一种方法可能是定义0左右的区域，这将被视为等于零（不，-或可忽略不计，效果）。这被称为实用等价区（ROPE），是测试参数重要性的一种方法。
我们如何定义这个区域？现在我们知道我们可以将ROPE定义为[-20,20]范围。此范围内的所有效果均视为无效（可忽略不计）。我们现在可以计算89％最可能值（89％CI）的比例，这些值不为空，即超出此范围。
```{r}
rope(posteriors$feedsunflower, range = c(-20, 20), ci=0.89)
```
89％CI中的7.75％可被视为无效。 那是很多吗？ 根据我们的指导方针，是的，这太过分了。 基于ROPE的这一特定定义，我们得出结论，这种影响并不显着（可忽略的概率太高）。
虽然，老实说，我对霍华德教授有些怀疑。 我真的不相信他对ROPE的定义。 有更客观的方法来定义它吗？
是。 其中一种做法是例如使用响应变量的标准偏差（SD）的十分之一（1/10 = 0.1），这可以被认为是“可忽略的”效应大小。
```{r}
rope_value <- 0.1 * sd(data$weight)
rope_range <- c(-rope_value, rope_value)
rope_range
```
让我们重新定义我们的ROPE作为[-6.2,6.2]范围内的区域。注意，这可以通过rope_range函数直接获得:
```{r}
rope_value <- rope_range(model)
rope_range
```
让我们重新计算ROPE中的百分比：
```{r}
rope(posteriors$feedsunflower, range = rope_range, ci=0.89)
```
根据ROPE的这种合理定义，我们观察到效应的89％的后验分布与ROPE不重叠。因此，我们可以得出结论，这种影响是显着的（在重要性的意义上要注意）。
### Probability of Direction方向概率（pd）
也许我们对这种影响是否可以忽略不感兴趣。也许我们只想知道这种影响是积极的还是消极的。 在这种情况下，无论效果的“大小”如何，我们都可以简单地计算后验的正比例。
```{r}
n_positive <- posteriors %>%
filter(feedsunflower > 0) %>% ## select only positive values
nrow() ## Get length
n_positive / nrow(posteriors) * 100
```
我们可以得出结论，效果是积极的，概率为97.82％。 我们称这个指数为方向概率（pd）。 事实上，它可以通过以下方式更容易地计算：
```{r}
p_direction(posteriors$feedsunflower)
```
有趣的是，这个指数通常与频率论的p值高度相关。我们几乎可以通过简单的转换大致推断出相应的p值：
```{r}
pd <- 97.82
onesided_p <- 1 - pd / 100
twosided_p <- onesided_p * 2
twosided_p
```
如果我们在频率框架中运行我们的模型，我们应该近似观察p值为0.04的效果。 真的吗？
### 与频率论者的比较
```{r}
lm(weight ~ feed, data=data) %>%
summary()
```
频率论模型告诉我们差异是积极的和显着的（$\beta$= 52，$p$ = 0.04）。
虽然我们得出了类似的结论，但贝叶斯框架使我们能够对我们的影响及其估计的不确定性有更深刻和直观的理解。
### 一个实现所有需求的功能
然而，我同意，提取和计算所有指数有点乏味。 但是如果我告诉你我们只用一个功能就可以完成所有这些工作呢？
此函数计算所有提到的索引，并且可以直接在模型上运行：
```{r}
describe_posterior(model)
```
我们有它！ 中位数，CI，pd和ROPE百分比！
理解和描述后验分布只是贝叶斯建模的一个方面......你准备好了吗？
## 示例2：确认贝叶斯技能
现在描述和理解线性回归的后验分布对你来说没有任何秘密，让我们回过头来研究一些更简单的模型： 相关性和t-测试 。
但在我们这样做之前，让我们花点时间提醒自己，并了解所有基本的统计程序 ，如相关性， t-测试，ANOVA或Chisquare测试都是线性回归 （我们强烈推荐这个优秀的演示 ）。 但是，这些简单模型仍将是引入更复杂指数的机会，例如贝叶斯因子 。
### 相关性
#### 频率版本
两个连续变量之间的频率相关性 ，一些花的萼片的宽度和长度 。 数据在R中可用作iris数据集（与前一个教程中使用的相同）。让我们计算一个Pearson相关性测试，将结果存储在一个名为result的对象中，然后显示它：
```{r}
result <- cor.test (iris $ Sepal.Width, iris $ Sepal.Length)
result
```
正如您在输出中所看到的，我们所做的测试实际上比较了两个假设： 零假设 （无相关）与备选假设 （非零相关）。基于p值，不能拒绝零假设：两个变量之间的相关性是负的但不显着 （r = -12，p> .05）。
#### 贝叶斯相关
要计算贝叶斯相关性测试，我们需要BayesFactor包，使用correlationBF()函数计算相关correlationBF()并以类似的方式存储结果。
```{r}
library (BayesFactor)
result <- correlationBF (iris $ Sepal.Width, iris $ Sepal.Length)
#让我们运行describe_posterior()函数：
describe_posterior (result)
```
我们在这里再看到很多东西，但现在重要的指数是后验分布的中位数为-0.11 。这非常接近频率相关性。 如前所述，我们可以描述可信区间 ， pd或ROPE百分比 ，但我们将重点关注贝叶斯框架贝叶斯因子提供的另一个指数。
### 贝叶斯因子（BF）
我们说，相关性实际上比较了两个假设，一个无效（没有效果）和一个altnernative（存在效应）。 贝叶斯因子（BF）允许相同的比较并且确定观察数据更可能在两个模型中的哪个模型下 ：具有感兴趣效果的模型和没有感兴趣效果的空模型。 我们可以使用bayesfactor()来专门计算贝叶斯因子，比较这些模型
```{r}
bayesfactor (result)
```
我们得到了BF=0.51 。 这是什么意思？
贝叶斯因子是相对证据的连续测量，贝叶斯因子大于1，证明有利于其中一个模型（通常称为分子），贝叶斯因子小于1，证明有利于另一个模型（分母）。
是的，你听到了正确的事，支持null的证据！这就是贝叶斯框架有时被认为优于频率框架的原因之一。 请记住，从您的统计数据课程中， p值只能用于拒绝h0 ，但不能接受它。 使用贝叶斯因子,您可以衡量证据支持null假设的情况 。
表示替代null假设的证据的BF可以使用来反转，$BF_{01}=1/BF_{10}$以提供替代null的证据。 在针对null假设的替代的BF小于1（支持null）的情况下，这提高了人的可读性。
在我们的例子中， BF = 1/0.51 = 2 ，表明与替代假设相比，数据null假设下的概率是2倍 ，尽管有利于空，但只被认为是反对空的轶事证据。因此，我们可以得出结论，有传闻证据表明两个变量之间缺乏相关性（r中位数 = 0.11，BF = 0.51） ，这是一个更具信息性的陈述，表明我们可以用频率统计数据做些什么。
这还不是全部！
## 示例3：成为贝叶斯主人
未完待续。
评论
发表评论