Bayesianische Statistik

class: center, middle, inverse, title-slide

# Bayesianische Statistik
## Teil 3<br/>Lineare Modelle mit brms
### Andrew Ellis
### Methodenkurs Neurowissenschaft im Computerlab
### 2021-03-16

---

layout: true
  

<div class="my-footer">
<span>
<a href="https://kogpsy.github.io/neuroscicomplab" target="_blank"><svg xmlns="http://www.w3.org/2000/svg" viewBox="0 0 576 512" class="rfa" style="height:0.75em;fill:#0F4C81;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>

---

## Mittelwert und Standardabweichung einer Normalverteilung schätzen

Wir haben 3 Datenpunkte mit den Werten `$y_1 = 85$`, `$y_2 = 100$`, und `$y_3 = 115$`. 
Wenn wir annehmen, dass diese Datenpunkte aus einer Normalverteilung kommen

$$ p(y | \mu, \sigma) = \frac{1}{Z} exp \left(- \frac{1}{2} \frac{(y-\mu)^2}{\sigma^2}\right)  $$

.footnote[ `$Z = \sigma \sqrt{2\pi}$` ist eine Normalisierungskonstante. 
]

können wir die Wahrscheinlichkeit eines Datenpunktes berechnen.

Unter der Annahme, dass die drei Datenpunkte unabhängig sind, ist die Wahrscheinlichkeit der Daten das Produkt der Einzelwahrscheinlichkeiten.

Wie aber finden wir die beiden Parameter der Normalverteilung, `$\mu$` und `$\sigma$`?

---

## Graphical model

<br>

.pull-left[
<div class="figure">
<img src="images/normal-graphical-model-2.png" alt="Graphical Model für normalverteilte Daten." width="50%" />
<p class="caption">Graphical Model für normalverteilte Daten.</p>
</div>
]

.pull-right[
$$ \mu \sim ???  $$
$$  \sigma \sim ??? $$
$$ y \sim N(\mu, \sigma) $$
]

---

.panelset[
.panel[.panel-name[Packages]

```r
library(tidyverse)
```
]
.panel[.panel-name[Kombinationen von mu und sigma]

```r
sequence_length <- 100

d1 <- crossing(y = seq(from = 50, to = 150, length.out = sequence_length),
              mu = c(87.8, 100, 112),
              sigma = c(7.35, 12.2, 18.4)) %>%
    mutate(density = dnorm(y, mean = mu, sd = sigma),
           mu = factor(mu, labels = str_c("mu==", c(87.8, 100, 112))),
           sigma = factor(sigma, 
                          labels = str_c("sigma==", c(7.35, 12.2, 18.4))))
```

]
.panel[.panel-name[Theme]

```r
theme_set(
  theme_bw() +
    theme(text = element_text(color = "white"),
          axis.text = element_text(color = "black"),
          axis.ticks = element_line(color = "white"),
          legend.background = element_blank(),
          legend.box.background = element_rect(fill = "white",
                                               color = "transparent"),
          legend.key = element_rect(fill = "white",
                                    color = "transparent"),
          legend.text = element_text(color = "black"),
          legend.title = element_text(color = "black"),
          panel.background = element_rect(fill = "white",
                                          color = "white"),
          panel.grid = element_blank()))
```
]
.panel[.panel-name[Plot Code]

```r
d1 %>% 
  ggplot(aes(x = y)) +
  geom_ribbon(aes(ymin = 0, ymax = density),
              fill = "steelblue") +
  geom_vline(xintercept = c(85, 100, 115), 
             linetype = 3, color = "white") +
  geom_point(data = tibble(y = c(85, 100, 115)),
             aes(y = 0.002),
             size = 2, color = "red") +
  scale_y_continuous(expression(italic(p)(italic(y)*"|"*mu*", "*sigma)), 
                     expand = expansion(mult = c(0, 0.05)), breaks = NULL) +
  ggtitle("Welche Normalverteilung?") +
  coord_cartesian(xlim = c(60, 140)) +
  facet_grid(sigma ~ mu, labeller = label_parsed) +
  theme_bw() +
  theme(panel.grid = element_blank())
```

]
.panel[.panel-name[Plot]
<img src="03-bayesian-stats-slides_files/figure-html/unnamed-chunk-6-1.png" width="100%" />
]
]

---

## Beispiel: t-Test als Lineares Modell

.alert[
Wir nehmen als Beispiel folgenden (fiktiven) Datensatz:

Wir haben IQ Scores einer Stichprobe, welche eine "Smart Drug" konsumiert hat. IQ Werte sind so normiert, dass wir in der Bevölkerung einen Mittelwert von 100, und eine SD von 15 haben. Nun wollen wir wissen, ob und wie sich die "Smart Drug" Gruppe von von diesem erwarteten Mittelwert unterscheidet. Gleichzeitig haben wir auch eine Gruppe, welche ein Placebo konsumiert hat. 
]

Die Daten sind aus dem Buch Doing Bayesian Data Analysis; wir können hier einfach die Daten von Hand eingeben, und zu einem Dataframe zusammenfügen.

.footnote[
Kruschke, John. Doing Bayesian Data Analysis (2nd Ed.). Boston: Academic Press, 2015.
]

---

.panelset[
.panel[.panel-name[2 Gruppen]

```r
smart = tibble(IQ = c(101,100,102,104,102,97,105,105,98,101,100,123,105,103,
                      100,95,102,106,109,102,82,102,100,102,102,101,102,102,
                      103,103,97,97,103,101,97,104,96,103,124,101,101,100,
                      101,101,104,100,101),
               Group = "SmartDrug")

placebo = tibble(IQ = c(99,101,100,101,102,100,97,101,104,101,102,102,100,105,
                        88,101,100,104,100,100,100,101,102,103,97,101,101,100,101,
                        99,101,100,100,101,100,99,101,100,102,99,100,99),
                 Group = "Placebo")
```
]

.panel[.panel-name[Dataframe]

```r
TwoGroupIQ <- bind_rows(smart, placebo) %>%
    mutate(Group = fct_relevel(as.factor(Group), "Placebo"))
```
]

.panel[.panel-name[Mittelwerte und SD]

```r
library(kableExtra)

TwoGroupIQ %>%
  group_by(Group) %>%
  summarise(mean = mean(IQ),
            sd = sd(IQ)) %>%
  mutate(across(where(is.numeric), round, 2)) %>% 
  kbl() %>%
  kable_styling(bootstrap_options = c("striped", "hover"))
```

<table class="table table-striped table-hover" style="margin-left: auto; margin-right: auto;">
 <thead>
  <tr>
   <th style="text-align:left;"> Group </th>
   <th style="text-align:right;"> mean </th>
   <th style="text-align:right;"> sd </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:left;"> Placebo </td>
   <td style="text-align:right;"> 100.36 </td>
   <td style="text-align:right;"> 2.52 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> SmartDrug </td>
   <td style="text-align:right;"> 101.91 </td>
   <td style="text-align:right;"> 6.02 </td>
  </tr>
</tbody>
</table>
]
]

---

## Smart Drug Gruppe

.pull-left[

```r
d <- TwoGroupIQ %>% 
  filter(Group == "SmartDrug") %>% 
  mutate(Group = fct_drop(Group))
```

```r
d %>% 
  ggplot(aes(x = IQ)) +
  geom_histogram(fill = "skyblue3", binwidth = 1) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.05))) + 
  theme_tidybayes()
```

]

.pull-right[
<img src="03-bayesian-stats-slides_files/figure-html/smart-drug-1.png" width="100%" />
]

---

## Parameter schätzen mit brms

### Priors

```r
library(brms)
```

```r
priors <- get_prior(IQ ~ 1,
          data = d)
```

---

```r
priors
```

```
##                 prior     class coef group resp dpar nlpar bound  source
##  student_t(3, 102, 3) Intercept                                  default
##    student_t(3, 0, 3)     sigma                                  default
```

---

## Priors

### `$\sigma$`

.panelset[
.panel[.panel-name[Code]

```r
tibble(x = seq(from = 0, to = 10, by = .025)) %>% 
  mutate(d = dt(x, df = 3)) %>% 
  ggplot(aes(x = x, ymin = 0, ymax = d)) +
  geom_ribbon(fill = "skyblue3") +
  scale_x_continuous(expand = expansion(mult = c(0, 0.05))) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.05)), breaks = NULL) +
  coord_cartesian(xlim = c(0, 8),
                  ylim = c(0, 0.35)) +
  xlab(expression(sigma)) +
  labs(subtitle = "Half-student-t Distribution: Prior für Standardabweichung.") +
  theme_bw(base_size = 14)
```
]
.panel[.panel-name[Plot]
<img src="03-bayesian-stats-slides_files/figure-html/unnamed-chunk-16-1.png" width="100%" />

]
]

---

## Priors

### `$\mu$`

.panelset[
.panel[.panel-name[Code]

```r
tibble(x = seq(from = 0, to = 200, by = .025)) %>% 
  mutate(d = dnorm(x, mean = 102, sd = 3)) %>% 
  ggplot(aes(x = x, ymin = 0, ymax = d)) +
  geom_ribbon(fill = "skyblue3") +
  scale_x_continuous(expand = expansion(mult = c(0, 0.05))) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.05)), breaks = NULL) +
  coord_cartesian(xlim = c(50, 150),
                  ylim = c(0, 0.15)) +
  xlab(expression(sigma)) +
  labs(subtitle = "Normalverteilter Prior für Mittelwert") +
  theme_bw(base_size = 14)
```
]
.panel[.panel-name[Plot]
<img src="03-bayesian-stats-slides_files/figure-html/unnamed-chunk-18-1.png" width="100%" />
]
]

---

## Prior Predictive Distribution

```r
m1_prior <- brm(IQ ~ 1,
          prior = priors,
          data = d,
          sample_prior = "only",
          file = "models/twogroupiq-prior-1")
```

---

## Prior Predictive Distribution

```r
summary(m1_prior)
```

```
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: IQ ~ 1 
##    Data: d (Number of observations: 47) 
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 4000
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept   101.95      6.20    92.34   111.71 1.00     2048     1510
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     3.75     17.71     0.10    12.94 1.00     1757     1199
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
```
 
---

## Prior Predictive Distribution

```r
plot(m1_prior)
```

---

## Prior Predictive Distribution

```r
library(tidybayes)

prior_pred_1 <- d %>%
  modelr::data_grid(Group) %>%
  add_predicted_draws(m1_prior) %>%
  ggplot(aes(y = Group, x = .prediction)) +
  stat_interval(.width = c(.50, .80, .95, .99)) +
  geom_point(aes(x = IQ), alpha = 0.4,  data = d) +
  scale_color_brewer() +
  theme_tidybayes()
```

---

## Vom Posterior Sampeln

.panelset[
.panel[.panel-name[Model]

```r
m1 <- brm(IQ ~ 1,
          prior = priors,
          data = d,
          file = "models/twogroupiq-1")
```
]
.panel[.panel-name[Plot]

```r
plot(m1)
```

<img src="03-bayesian-stats-slides_files/figure-html/unnamed-chunk-25-1.png" width="100%" />
]

.panel[.panel-name[Summary]

```r
print(m1)
```

```
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: IQ ~ 1 
##    Data: d (Number of observations: 47) 
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 4000
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept   101.92      0.83   100.31   103.53 1.00     3134     2075
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     6.04      0.64     4.96     7.47 1.00     3060     2458
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
```

]]

---

## Posterior Quantile Intervals

```r
library(patchwork)
mcmc_plot(m1, pars = "b") / mcmc_plot(m1, pars = "sigma")
```

---

## Posterior Samples

.panelset[
.panel[.panel-name[Samples extrahieren]

```r
samples <- posterior_samples(m1) %>% 
  transmute(mu = b_Intercept, sigma = sigma)
```
]

.panel[.panel-name[Credible Intervals]

```r
library(tidybayes)

samples %>% 
  select(mu) %>% 
  median_qi(.width = c(.50, .80, .95, .99))
```

```
##         mu    .lower   .upper .width .point .interval
## 1 101.9132 101.36505 102.4841   0.50 median        qi
## 2 101.9132 100.86039 102.9730   0.80 median        qi
## 3 101.9132 100.30502 103.5255   0.95 median        qi
## 4 101.9132  99.71696 104.1686   0.99 median        qi
```
]

.panel[.panel-name[Plot Code]

```r
samples %>% 
  select(mu) %>% 
  median_qi(.width = c(.50, .80, .95, .99)) %>% 
  ggplot(aes(x = mu, xmin = .lower, xmax = .upper)) +
  geom_pointinterval() +
  ylab("") +
  theme_tidybayes()
```
]

.panel[.panel-name[Plot]

<img src="03-bayesian-stats-slides_files/figure-html/unnamed-chunk-31-1.png" width="100%" />
]
]
---

## Posterior Samples

.panelset[
.panel[.panel-name[Density Plot Code: MW]

```r
samples %>% 
  select(mu) %>% 
  ggplot(aes(x = mu)) +
  stat_halfeye() +
  theme_tidybayes()
```
]

.panel[.panel-name[Density Plot: MW]

<img src="03-bayesian-stats-slides_files/figure-html/unnamed-chunk-33-1.png" width="100%" />
]

.panel[.panel-name[Density Plot Code: SD]

```r
samples %>% 
  select(sigma) %>% 
  ggplot(aes(x = sigma)) +
  stat_halfeye(point_interval = mode_hdi) +
  theme_tidybayes()
```
]

.panel[.panel-name[Density Plot: SD]

<img src="03-bayesian-stats-slides_files/figure-html/unnamed-chunk-35-1.png" width="100%" />
]]

---

## Posterior Predictive Distribution

.pull-left[

```r
post_pred_1 <- d %>%
  modelr::data_grid(Group) %>%
  add_predicted_draws(m1) %>%
  ggplot(aes(y = Group, x = .prediction)) +
  stat_interval(.width = c(.50, .80, .95, .99)) +
  geom_point(aes(x = IQ), alpha = 0.4,  data = d) +
  scale_color_brewer() +
  theme_tidybayes()
```
]

.pull-right[

```r
post_pred_1
```

<img src="03-bayesian-stats-slides_files/figure-html/unnamed-chunk-37-1.png" width="100%" />
]

---

## Prior vs Posterior Predictive Distribution

```r
cowplot::plot_grid(prior_pred_1, 
                   post_pred_1, 
                   labels = c('Prior predictive', 'Posterior predictive'), 
                   label_size = 12,
                   align = "h",
                   nrow = 2)
```

---

## Zwei Gruppen

---

## t-Test als Allgemeines Lineares Modell

- Annahme: Varianzgleichheit

### Tradionelle Notation

`$$y_{ij} = \alpha + \beta x_{ij} + \epsilon_{ij}$$`
`$$\epsilon \sim N(0, \sigma^2)$$`
--

### Probabilistische Notation

`$$y_{ij} \sim N(\mu, \sigma^2)$$`
`$$\mu_{ij} = \alpha + \beta x_{ij}$$`

`$X_{ij}$` ist eine Indikatorvariable.

---

## Ordinary Least Squares Schätzung

```r
levels(TwoGroupIQ$Group)
```

```
## [1] "Placebo"   "SmartDrug"
```

Mit R Formelnotation:

.panelset[

.panel[.panel-name[LM]

```r
fit_ols <- lm(IQ ~ Group,
              data = TwoGroupIQ)
```
]

.panel[.panel-name[Output]

```r
summary(fit_ols)
```

```
## 
## Call:
## lm(formula = IQ ~ Group, data = TwoGroupIQ)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -19.9149  -0.9149   0.0851   1.0851  22.0851 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    100.3571     0.7263 138.184   <2e-16 ***
## GroupSmartDrug   1.5578     0.9994   1.559    0.123    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.707 on 87 degrees of freedom
## Multiple R-squared:  0.02717,	Adjusted R-squared:  0.01599 
## F-statistic:  2.43 on 1 and 87 DF,  p-value: 0.1227
```
]
]

---

## Dummy Coding in R

.pull-left[

```r
contrasts(TwoGroupIQ$Group)
```

```
##           SmartDrug
## Placebo           0
## SmartDrug         1
```

```r
mm1 <- model.matrix(~ Group, data = TwoGroupIQ)
head(mm1)
```

```
##   (Intercept) GroupSmartDrug
## 1           1              1
## 2           1              1
## 3           1              1
## 4           1              1
## 5           1              1
## 6           1              1
```
]

.pull-right[

```r
as_tibble(mm1) %>% 
  group_by(GroupSmartDrug) %>% 
  slice_sample(n= 3)
```

```
## # A tibble: 6 x 2
## # Groups:   GroupSmartDrug [2]
##   `(Intercept)` GroupSmartDrug
##           <dbl>          <dbl>
## 1             1              0
## 2             1              0
## 3             1              0
## 4             1              1
## 5             1              1
## 6             1              1
```
]

---

## Dummy Coding in R

.pull-left[

```r
as_tibble(mm1) %>% 
  group_by(GroupSmartDrug) %>% 
  slice_sample(n= 3)
```

.pull-right[

```r
mm2 <- model.matrix(~ 0 + Group, data = TwoGroupIQ)
as_tibble(mm2) %>% 
  group_by(GroupSmartDrug) %>% 
  slice_sample(n= 3)
```

```
## # A tibble: 6 x 2
## # Groups:   GroupSmartDrug [2]
##   GroupPlacebo GroupSmartDrug
##          <dbl>          <dbl>
## 1            1              0
## 2            1              0
## 3            1              0
## 4            0              1
## 5            0              1
## 6            0              1
```
]

---

## Graphical Model

<div class="figure">
<img src="images/two-group-iq-graphical-model.png" alt="Graphical Model für 2 Gruppen." width="20%" />
<p class="caption">Graphical Model für 2 Gruppen.</p>
</div>

---

## Get Priors

```r
priors2 <- get_prior(IQ ~ 1 + Group,
                     data = TwoGroupIQ)
```

```r
priors2
```

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

---

## Get Priors

```r
priors3 <- get_prior(IQ ~ 0 + Group,
                     data = TwoGroupIQ)
```

```r
priors3
```

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

---

## Priors definieren und vom Prior Sampeln

```r
priors2_b <- prior(normal(0, 2), class = b)
```

```r
m2_prior <- brm(IQ ~ 1 + Group,
          prior = priors2_b,
          data = TwoGroupIQ,
          sample_prior = "only",
          file = "models/twogroupiq-prior-2")
```

---

```r
prior_pred_2 <- TwoGroupIQ %>%
  modelr::data_grid(Group) %>%
  add_predicted_draws(m2_prior) %>%
  ggplot(aes(y = Group, x = .prediction)) +
  stat_interval(.width = c(.50, .80, .95, .99)) +
  geom_point(aes(x = IQ), alpha = 0.4,  data = TwoGroupIQ) +
  scale_color_brewer() +
  theme_tidybayes()
```

---

## Priors definieren und vom Prior Sampeln

```r
priors3_b <- prior(normal(100, 10), class = b)
```

```r
m3_prior <- brm(IQ ~ 0 + Group,
          prior = priors3_b,
          data = TwoGroupIQ,
          sample_prior = "only",
          file = "models/twogroupiq-prior-3")
```

---

```r
prior_pred_3 <- TwoGroupIQ %>%
  modelr::data_grid(Group) %>%
  add_predicted_draws(m3_prior) %>%
  ggplot(aes(y = Group, x = .prediction)) +
  stat_interval(.width = c(.50, .80, .95, .99)) +
  geom_point(aes(x = IQ), alpha = 0.4,  data = TwoGroupIQ) +
  scale_color_brewer() +
  theme_tidybayes()
```

---

## Posterior Sampling

```r
m2 <- brm(IQ ~ 1 + Group,
          prior = priors2_b,
          data = TwoGroupIQ,
          file = "models/twogroupiq-2")
```

```r
m3 <- brm(IQ ~ 0 + Group,
          prior = priors3_b,
          data = TwoGroupIQ,
          file = "models/twogroupiq-3")
```

---

## Summary

```r
summary(m2)
```

```
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: IQ ~ 1 + Group 
##    Data: TwoGroupIQ (Number of observations: 89) 
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 4000
## 
## Population-Level Effects: 
##                Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept        100.51      0.69    99.18   101.90 1.00     3735     2921
## GroupSmartDrug     1.23      0.88    -0.49     2.95 1.00     3633     3160
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     4.72      0.36     4.09     5.52 1.00     3002     2088
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
```

---

## Summary

```r
summary(m3)
```

```
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: IQ ~ 0 + Group 
##    Data: TwoGroupIQ (Number of observations: 89) 
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 4000
## 
## Population-Level Effects: 
##                Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## GroupPlacebo     100.35      0.73    98.95   101.75 1.00     3914     3163
## GroupSmartDrug   101.92      0.69   100.62   103.28 1.00     4057     2954
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     4.71      0.36     4.07     5.45 1.00     3342     2697
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
```

---

```r
mcmc_plot(m2, "b_GroupSmartDrug")
```

---

```r
mcmc_plot(m3, "b")
```

---
## Get Posterior Samples

```r
samples_m3 <- posterior_samples(m3) %>% 
    transmute(Placebo = b_GroupPlacebo, 
              SmartDrug = b_GroupSmartDrug,
              sigma = sigma)
```

```r
samples_m3 <- samples_m3 %>% 
  mutate(diff = SmartDrug - Placebo,
         effect_size = diff/sigma)
```

---

## Get Posterior Samples

```r
samples_m3 %>% 
  select(diff) %>% 
  median_qi()
```

```
##      diff     .lower   .upper .width .point .interval
## 1 1.57189 -0.3887631 3.521777   0.95 median        qi
```

---

```r
samples_m3 %>% 
  select(diff) %>% 
  ggplot(aes(x = diff)) +
  stat_halfeye(point_interval = median_qi) +
  theme_tidybayes()
```

---

```r
samples_m3 %>% 
  select(effect_size) %>% 
  ggplot(aes(x = effect_size)) +
  stat_halfeye(point_interval = median_qi) +
  theme_tidybayes()
```

---

## Diesmal ohne eigene Priors

```r
fit_eqvar <- brm(IQ ~ Group,
                 data = TwoGroupIQ,
                 file = here::here("models/fit_eqvar"))
```

---

```r
fit_eqvar %>%
    gather_draws(b_GroupSmartDrug) %>%
    ggplot(aes(y = .variable, x = .value)) +
    stat_halfeye(fill = "Steelblue4") +
    geom_vline(xintercept = 0, color = "white", linetype = 1, size = 1) +
    ylab("") +
    xlab("Estimated difference") +
    theme_tidybayes()
```

---

.panelset[
.panel[.panel-name[Code]

```r
grid <- TwoGroupIQ %>%
    modelr::data_grid(Group)

fits_IQ <- grid %>%
    add_fitted_draws(fit_eqvar)

preds_IQ <- grid %>%
    add_predicted_draws(fit_eqvar)

pp_eqvar <- TwoGroupIQ %>%
    ggplot(aes(x = IQ, y = Group)) +
    stat_halfeye(aes(x = .value),
                  scale = 0.7,
                  position = position_nudge(y = 0.1),
                  data = fits_IQ,
                  .width = c(.66, .95, 0.99)) +
    stat_interval(aes(x = .prediction),
                   data = preds_IQ,
                   .width = c(.66, .95, 0.99)) +
    scale_x_continuous(limits = c(75, 125)) +
    geom_point(data = TwoGroupIQ) +
    scale_color_brewer() +
	labs(title = "Equal variance model predictions") +
  theme_tidybayes()
```

]
.panel[.panel-name[Plot]

<img src="03-bayesian-stats-slides_files/figure-html/unnamed-chunk-74-1.png" width="100%" />
]
]