Intro to Bayesian Statistics

class: center, middle, inverse, title-slide

# Intro to Bayesian Statistics
## Part 2 <br/> Linear models
### Andrew Ellis
### Bayesian multilevel modelling workshop 2021
### 05-21-2021

---

layout: true
  

<div class="my-footer">
<span>
<a href="https://awellis.github.io/learnmultilevelmodels/" 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> Graduate School workshop 2021
</span>
</div>

---

## Estimating parameters of a normal distribution

We have 3 data points, `$y_1 = 85$`, `$y_2 = 100$`, und `$y_3 = 115$`. We can assume that they are normally distributed.

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

**Goal:** We want to estimate the mean and standard deviation of a normal distribution

.footnote[ `$Z = \sigma \sqrt{2\pi}$` is a normalising constant. 
]

We can compute the probability of one data point, and if the data points are `$i.i.d$`, then `$P(data|, \mu, \sigma)$` is the product of the indidividual probabilities.

But how do we find `$\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[Combinations of mu and 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("Which distribution?") +
  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="01-linear-models-slides_files/figure-html/unnamed-chunk-6-1.png" width="100%" />
]
]

---

## Example: t-Test as a linaer model

.alert[
We will use some simulated data:

we have two groups. One group were given a smart drug, the other is the control (placebo) group. IQ values have a mean of 100, and an SD of 15

We would like to know if the smart group is smarter than the control group.
]

.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[Mean and 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 Group

.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)))
```

]

.pull-right[
<img src="01-linear-models-slides_files/figure-html/smart-drug-1.png" width="100%" />
]

---

## Parameter estimation using 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 prior on standard deviation") +
  theme_bw(base_size = 14)
```
]
.panel[.panel-name[Plot]
<img src="01-linear-models-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 = "Normal prior on mean") +
  theme_bw(base_size = 14)
```
]
.panel[.panel-name[Plot]
<img src="01-linear-models-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()
```

---

## Sampling from the Posterior

.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="01-linear-models-slides_files/figure-html/unnamed-chunk-25-1.png" width="100%" />
]

.panel[.panel-name[Summary]

```r
summary(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[Extract samples]

```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("")
```
]

.panel[.panel-name[Plot]

<img src="01-linear-models-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()
```
]

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

<img src="01-linear-models-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)
```
]

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

<img src="01-linear-models-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()
```
]

.pull-right[

```r
post_pred_1
```

<img src="01-linear-models-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)
```

---

## Two Groups

---

## t-Test as  General Linear Model

- For simplicity, we assume equal variances

### Tradional notation

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

### Probabilistic notation

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

`$X_{ij}$` is an indicator variable.

---

## Ordinary Least Squares

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

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

Using R's formula notation:

.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

```r
knitr::include_graphics("images/two-group-iq-graphical-model.png")
```

---

## 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
```

---

## Define Priors on Group Means

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

We can sample from the prior.

```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) 
```

---

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

---

## Using default priors

```r
fit_eqvar <- brm(IQ ~ Group,
                 data = TwoGroupIQ,
                 file = "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="01-linear-models-slides_files/figure-html/unnamed-chunk-75-1.png" width="100%" />
]
]