Bayesianische Statistik

# Bayesianische Statistik
## Teil 5b <br/> Bayes Factors
### Andrew Ellis
### Methodenkurs Neurowissenschaft im Computerlab
### 2021-04-27

---

layout: true
  

<div class="my-footer">
<span>
<a href="https://kogpsy.github.io/neuroscicomplab" target="_blank"><svg aria-hidden="true" role="img" viewBox="0 0 576 512" style="height:1em;width:1.12em;vertical-align:-0.125em;margin-left:auto;margin-right:auto;font-size:inherit;fill:#0F4C81;overflow:visible;position:relative;"><path d="M280.37 148.26L96 300.11V464a16 16 0 0 0 16 16l112.06-.29a16 16 0 0 0 15.92-16V368a16 16 0 0 1 16-16h64a16 16 0 0 1 16 16v95.64a16 16 0 0 0 16 16.05L464 480a16 16 0 0 0 16-16V300L295.67 148.26a12.19 12.19 0 0 0-15.3 0zM571.6 251.47L488 182.56V44.05a12 12 0 0 0-12-12h-56a12 12 0 0 0-12 12v72.61L318.47 43a48 48 0 0 0-61 0L4.34 251.47a12 12 0 0 0-1.6 16.9l25.5 31A12 12 0 0 0 45.15 301l235.22-193.74a12.19 12.19 0 0 1 15.3 0L530.9 301a12 12 0 0 0 16.9-1.6l25.5-31a12 12 0 0 0-1.7-16.93z"/></svg></a> Methodenkurs Neurowissenschaft im Computerlab
</span>
</div>

---

## Schätzung von Bayes Factors

1) Savage-Dickey Density Ratio

2) `bridge sampling` Package

3) `BayesFactor` Package

4) JASP: [jasp-stats.org/download](https://jasp-stats.org/download/)

---

Wir nehme als Beispiel einen simulierten Datensatz vom Package [bcogsci](https://github.com/bnicenboim/bcogsci). Falls Sie das Package nicht installieren können, können Sie einfach das File `df_contrasts1.rda` von hier [downloaden](https://github.com/bnicenboim/bcogsci/tree/master/data).

Wir haben eine abhängige Variable `DV` (Reaktionzeit) in 2 Gruppen, F1 mit `$\mu_1=0.8$` Sekunden. und F2 mit `$\mu_2=0.4$` Sekunden. Die Daten sind von 10 simulierten Vpn.

---

```r
library(bcogsci)
data("df_contrasts1") 
df_contrasts1
```

```
## # A tibble: 10 x 3
##    F        DV    id
##    <fct> <dbl> <int>
##  1 F1    0.636     1
##  2 F1    0.841     2
##  3 F1    0.555     3
##  4 F1    1.03      4
##  5 F1    0.938     5
##  6 F2    0.123     6
##  7 F2    0.304     7
##  8 F2    0.659     8
##  9 F2    0.469     9
## 10 F2    0.444    10
```
]

```r
df_contrasts1 <- load("data/df_contrasts1.rda")
df_contrasts1
```
]

```r
df_contrasts1_sum <- df_contrasts1 %>% 
  group_by(F) %>% 
  summarise(mean = mean(DV),
            sd = sd(DV),
            se = sd/sqrt(n()))
df_contrasts1_sum
```

```
## # A tibble: 2 x 4
##   F      mean    sd     se
##   <fct> <dbl> <dbl>  <dbl>
## 1 F1      0.8 0.200 0.0894
## 2 F2      0.4 0.200 0.0894
```

]

```r
df_contrasts1_sum %>% 
  ggplot(aes(F, mean)) +
  geom_errorbar(aes(ymin = mean-se, ymax = mean+se),
                width = 0.2, size=1, color="blue") +
  geom_point(size = 4) +
  geom_line(aes(group = 1), linetype = 3)
```

]]

---

## Parameterschätzung

```r
library(brms)

get_prior(DV ~ 1 + F,
          data = df_contrasts1)
```

```
##                   prior     class coef group resp dpar nlpar bound       source
##                  (flat)         b                                       default
##                  (flat)         b  FF2                             (vectorized)
##  student_t(3, 0.6, 2.5) Intercept                                       default
##    student_t(3, 0, 2.5)     sigma                                       default
```

---

## Parameterschätzung

```r
priors <- prior(normal(0, 1), class = b)
```

```r
fit1 <- brm(DV ~ 1 + F,
            data = df_contrasts1,
            family = gaussian(),
            prior = priors,
            cores = parallel::detectCores(),
            file = "models/BF_fit_F1", 
            file_refit = "on_change") 
```

---

```r
round(fixef(fit1),3)
```

```
##           Estimate Est.Error   Q2.5  Q97.5
## Intercept    0.794     0.112  0.551  1.023
## FF2         -0.391     0.158 -0.699 -0.074
```

```r
mcmc_plot(fit1)
```

---

.pull-left-narrow[
  .huge-blue-number[1]
]
.pull-right-wide[
  .larger[
  Savage-Dickey Density Ratio
  ]
]

---

## Savage-Dickey Density Ratio

```r
fit2 <- brm(DV ~ 1 + F,
            data = df_contrasts1,
            family = gaussian(),
            prior = priors,
*           sample_prior = TRUE,
            cores = parallel::detectCores(),
            file = "models/BF_fit_F2", 
            file_refit = "on_change") 
```

---

## Savage-Dickey Density Ratio

```r
fit2 %>% 
  mcmc_plot(c("b_FF2", "prior_b"))
```

---

## Savage-Dickey Density Ratio

```r
h1 <- fit2 %>% 
  hypothesis("FF2 = 0")
h1
```

```
## Hypothesis Tests for class b:
##   Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob Star
## 1  (FF2) = 0    -0.39      0.15    -0.69     -0.1       0.26      0.21    *
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
```

---

## Savage-Dickey Density Ratio

```r
BF01 <- h1$hypothesis$Evid.Ratio
BF01
```

```
## [1] 0.2647813
```

```r
BF10 <- 1/BF01
BF10
```

```
## [1] 3.776701
```

---

---

## Bridge sampling

```r
fit3 <- brm(DV ~ 1 + F,
            data = df_contrasts1,
            family = gaussian(),
            prior = priors,
*           save_pars = save_pars(all = TRUE),
*           iter = 1e4,
            cores = parallel::detectCores(),
            file = "models/BF_fit_F3", 
            file_refit = "on_change") 
```

---

## Bridge sampling

Nullmodell (Vergleichsmodell) mit default Priors für Intercept und sigma.

```r
fit4 <- brm(DV ~ 1,
            data = df_contrasts1,
            family = gaussian(),
*           save_pars = save_pars(all = TRUE),
*           iter = 1e4,
            cores = parallel::detectCores(),
            file = "models/BF_fit_F4", 
            file_refit = "on_change") 
```

---

## Bridge sampling

```r
BF <- bayes_factor(fit3, fit4)
```

```r
BF
```

```
## Estimated Bayes factor in favor of fit3 over fit4: 3.30777
```

Bayes factor für Nullmodell: `$BF_{01}$`

```r
# Estimated Bayes factor in favor of fit4 over fit3
1/BF$bf
```

```
## [1] 0.3023181
```

---

---

## BayesFactor Package

```r
library(BayesFactor)
bf = ttestBF(formula = DV ~ F, data = df_contrasts1)

bf
```

```
## Bayes factor analysis
## --------------
## [1] Alt., r=0.707 : 4.29645 ±0%
## 
## Against denominator:
##   Null, mu1-mu2 = 0 
## ---
## Bayes factor type: BFindepSample, JZS
```

---

---

## JASP

- Daten als CSV exportieren:

```r
library(readr)

df_contrasts1 <- read_csv("https://raw.githubusercontent.com/kogpsy/neuroscicomplab/main/data/df_contrasts1.csv")

df_contrasts1 %>% 
  write_csv(file = "data/df_contrasts1.csv")
```

- Öffnen Sie den Datensatz in JASP

---

- Sie sollten zuerst den BF vom `BayesFactor` Package replizieren.
- Stellen Sie Prior und Posterior grafisch dar.
- F¨ühren Sie eine **Robustness Check** durch.
- Führen Sie eine **Sequential Analysis** durch.
]

---

## Bayesian Anova

.your-turn[
Öffnen Sie in JASP den `Bugs` Datensatz (ANOVA), und versuchen Sie Bayes Factors zu schätzen.

Vpn mussten Ratings angeben, wie sehr sie [Arthropoden](https://de.wikipedia.org/wiki/Gliederf%C3%BC%C3%9Fer) töten wollten.

- Versuchen Sie, eine repeated-measures Bayesian ANOVA durchzuführen.
- Einen Beispielbericht finden Sie hier: [https://osf.io/wae57/](https://osf.io/wae57/).
]