dojo_multilevel.Rmd

---
title: "Multilevel models and metaanalysis"
author: "Martin Modrák"
#date: "2016/12/12 (updated: `r Sys.Date()`)"
output:
  xaringan::moon_reader:
    lib_dir: libs
    nature:
      highlightStyle: github
      highlightLines: true
      countIncrementalSlides: false
---

# Intro

- Goes by many names: 
  - Multilevel models 
  - Random effects 
  - Varying intercept/slope 
  - Mixed models 
  - Metaanalysis
- Packages:
  - `lme4` (frequentist)
  - `rstanarm` (Bayesian, easy to install)
  - `brms` (Bayesian, very flexible)

```{r setup, include=FALSE}
options(htmltools.dir.version = FALSE)
knitr::opts_chunk$set(echo = FALSE, fig.height = 6, fig.width = 9)
library(rstanarm)
library(tidybayes)
library(tidyverse)
library(cowplot)
library(scico)
library(patchwork)
theme_set(theme_cowplot())
output_dir <- here::here("local_temp_data")
if(!dir.exists(output_dir)) {
  dir.create(output_dir)
}
options(mc.cores = parallel::detectCores(), brms.backend = "cmdstanr")
```

---
class: inverse, center, middle



# Linear models recap

---

## Model formulation

$$\mu = \beta X$$
--
$$\mu_i = \beta_0 + \beta_{age} X_{i,age} + \beta_{treatment} X_{i,treatment}$$

--

$$y_i \sim N(\mu_i, \sigma)$$

---

## Model formulation (R)


$$\mu_i = \beta_0 + \beta_{age} X_{i,age} + \beta_{treatment} X_{i,treatment}$$

--

```{r}
mock_lin_data <- data.frame(y = round(rnorm(3, 40, 10), 1), age = c(23, 28, 30), group = c("control", "control", "treatment"))
mock_lin_data
```

--

R formula 

```
y ~ 1 + age + group
```

--

```{r, echo = TRUE}
model.matrix(y ~ 1 + age + group, mock_lin_data)
```


---

## Everything is a linear model!

`t.test(groupA, groupB)` = `lm(y ~ group)` 

---
class:inverse, middle, center

# A simple example

```{r}
# Helper plotting functions
base_plot <- function(data) {
  set.seed(9882264) #Force seed to maintain identical jitter
  data %>% 
    ggplot(aes(x = group, y = log_relative_expression, color = group, fill = group)) +
    scale_y_continuous("Log(Relative expression)") +
    # ggplot(aes(x = group, y = log_relative_expression, color = group, fill = group)) + 
    # scale_y_continuous("Log(Relative expression)") +
    theme(axis.title.x = element_blank())
}

add_my_point <- function(plot) {
  plot + 
    geom_point( position = position_jitter(width = 0.05), size = 4, shape = 21) +
    scale_color_manual(values = c("#c46e9b","#05a8aa"), guide = FALSE) + 
    scale_fill_manual(values = c("#f49ac8", "#5fc7c8"), guide = FALSE)
}

mean_se_log <- function(x) {
  exp(mean_se(log(x)))
}

add_error_bars <- function(plot, nudge = -0.3, width = 0.2) {
  position = position_nudge(x = nudge)
  plot +   
    stat_summary(fun.data = "mean_se", geom = "errorbar", width = width, color = "#707070", position = position) + 
    stat_summary(fun = "mean", geom = "point", shape = 3,  size = 3, stroke = 3, color = "#707070", position = position) 
}

```



---

## Single group of data

```{r}
batches <- c("A","B", "C")
set.seed(5432188)
base_data_single <- tibble(batch = batches, var_intercept = c(0.2, 0, -0.2), group ="") %>% 
  crossing(tibble(replicate = 1:6)) %>%
  filter( (batch != "B" | replicate < 5) & (batch != "C" | replicate < 4)) %>%
  mutate(log_relative_expression = var_intercept + rnorm(n(), 0, 0.25), 
         relative_expression = exp(log_relative_expression))


base_data_single %>% 
  base_plot() %>%
  add_my_point()


```

---

## Breaking it down

```{r}
base_data_single %>% base_plot() %>%
  add_my_point() +
  facet_wrap(~batch)
```

---

## A set of bad options

A) Ignore batches

$$\mu_i = \alpha$$

```{r, echo = TRUE}
ignore <- t.test(log_relative_expression ~ 1, 
                data = base_data_single)
ignore$conf.int
```

--

```{r, fig.height = 3.5, fig.width = 6}
plot_single_all_together <-  base_data_single %>% 
  base_plot() %>%
  add_error_bars() %>%
  add_my_point()

plot_single_all_together
```


---

## A set of bad options

B) Treat batches as separate

$$\mu_i = \alpha + \beta_1 \mathrm{isB}_i + \beta_2 \mathrm{isC}_i$$

```{r, echo = TRUE}
separate <- lm(log_relative_expression ~ batch, 
               data = base_data_single) 
confint(separate)
```

--

```{r, fig.height = 3.5, fig.width = 6}
base_data_single %>% base_plot() %>%
  add_error_bars() %>% add_my_point() +
  facet_wrap(~batch)
```


---

## A set of bad options

C) Take averages of the batches

```{r}
batch_means <- base_data_single  %>%
  group_by(batch, group) %>%
  summarise(log_relative_expression = mean(log_relative_expression), .groups = "drop") 

```

```{r, echo = TRUE}
take_means <- t.test(log_relative_expression ~ 1, 
                     data = batch_means)
take_means$conf.int
```

--

```{r, fig.height = 3.5, fig.width = 6}

plot_batch_means <-  batch_means %>%
  base_plot() %>%
  add_error_bars() %>%
  add_my_point()

plot_batch_means

```

---

## Varying intercept

The batches are not identical, but also not independent

$$\mu_i = \alpha + \beta_{batch(i)} \\
\beta_{batch(1,2,3)} \sim N(0, \sigma_{batch})$$

--

```{r varying_intercept_rstanarm, cache=TRUE, echo = TRUE}
varying_intercept <- rstanarm::stan_lmer(
  log_relative_expression ~ 1 + (1 | batch), 
  data = base_data_single)
summ_varying_intercept <- 
  summary(varying_intercept, probs = c(0.025,0.975))
summ_varying_intercept[1:6, c("2.5%", "97.5%")]

```

???

More parameters, less flexible

---

## A middle ground

```{r}
predict_data <- tibble(batch = batches) %>% crossing(group = unique(base_data_single$group))
s <- posterior_linpred(varying_intercept, newdata = predict_data)
varying_intercept_to_plot <- predict_data %>% mutate(log_mean = colMeans(s), mean = exp(log_mean), sd = sqrt(diag(cov(s))), low = log_mean - sd, high = log_mean + sd)


add_my_crossbar <- function(plot, crossbar_data) { 
  plot + geom_crossbar(aes(y = log_mean, ymin = low, ymax = high), data = crossbar_data, color = "#5d2e8c", fill = "white", width = 0.2, size = 2, position = position_nudge(x = 0.3)) 
}

plot_multilevel <- base_data_single %>% 
  base_plot() %>%
  add_error_bars() %>%
  add_my_crossbar(varying_intercept_to_plot)  %>% add_my_point() + facet_wrap(~batch)
  

plot_multilevel
```

---

## Partial pooling

```{r}
hide_axis_theme <- theme(axis.text.y = element_blank(), axis.title.y = element_blank(), axis.line.y = element_blank(),
                         axis.ticks.y = element_blank())
(((plot_batch_means + (plot_single_all_together+ hide_axis_theme)) + (plot_multilevel + hide_axis_theme)) & expand_limits( y = c(-0.55, 0.5))) + plot_layout(widths = c(0.25,0.25,0.5))
```



---

## Different prediction tasks

```{r}
predict_data_new_level <- tibble(group = unique(base_data_single$group), batch = "NEW_BATCH")

s_n <- posterior_linpred(varying_intercept, newdata = predict_data_new_level)
varying_intercept_to_plot_n <- predict_data_new_level %>% mutate(log_mean = colMeans(s_n), mean = exp(log_mean), sd = sqrt(diag(cov(s_n))), low = log_mean - sd, high = log_mean + sd)


plot_multilevel_intercept <- base_data_single %>% base_plot() %>%
  add_error_bars() %>%
  add_my_crossbar(varying_intercept_to_plot_n)  %>%
  add_my_point()

plot_multilevel_intercept + plot_multilevel & expand_limits( y = c(-0.55, 0.5))

```



---
class: inverse, middle, center

# More complex example

---

## All data together

```
log_relative_expression ~ group
```


$$\mu_i = \alpha + \beta_1 \mathrm{isTreatment}_i$$


```{r}
set.seed(6546542)
batches <- c("A","B", "C")
base_data <- tibble(batch = batches, var_intercept = c(0.2, 0, -0.2), var_effect = c(-0.3,0.1, 0.2)) %>% 
  crossing(tibble(group = c("control", "treatment"), b = c(0, 1))) %>%
  crossing(tibble(replicate = 1:10)) %>%
  filter( (batch != "B" | replicate < 7) & (batch != "C" | replicate < 4)) %>%
  #filter(batch == "A" | replicate < 3) %>%
  mutate(log_relative_expression = var_intercept + b* (0.2 + var_effect)  + rnorm(n(), 0, 0.25), 
         relative_expression = exp(log_relative_expression))


plot_together_points <- base_data %>% 
  base_plot() %>%
  add_error_bars() %>%
  add_my_point()

plot_together_points

```


---

## Break it down

```
log_relative_expression ~ batch * group
```

$$\mu_i = \alpha + \beta_1 \mathrm{isTreatment}_i + \beta_2 \mathrm{isB}_i + \beta_3 \mathrm{isC}_i + \\
\beta_4 \mathrm{isTreatment}_i \mathrm{isB}_i + \beta_5  \mathrm{isTreatment}_i \mathrm{isC}_i$$



```{r, fig.width = 6, fig.height = 5}
plot_separate <- base_data %>% base_plot() %>%
  add_error_bars() %>% add_my_point() +
  facet_wrap(~batch)

plot_separate
```

???

I can also make the effect be the same in all groups

---

## Varying intercept + varying effect

```
log_relative_expression ~ group + (1 + group || batch)
```


$$\mu_i = \alpha + \beta_1 \mathrm{isTreatment}_i + \gamma_{batch(i)} + \delta_{batch(i)} \mathrm{isTreatment}_i \\
\gamma_{batch(1,2,3)} \sim N(0, \sigma_{\mathrm{batch\_intercept}}) \\ 
\delta_{batch(1,2,3)} \sim N(0, \sigma_{\mathrm{batch\_effect}})$$

```{r}
# \begin{pmatrix}\gamma_{batch(1,2,3)} \\ \delta_{batch(1,2,3)} \end{pmatrix} 
# \sim MVN(0, \Sigma)
```


```{r varying_effect_brms, cache=TRUE, echo = TRUE, results = "hide", message = FALSE, warning = FALSE}
varying_effect <- brms::brm(
  log_relative_expression ~ group + (1 + group || batch), 
  data = base_data,
  prior = c(brms::set_prior("normal(0,1)", class = "b"),
            brms::set_prior("normal(0,1)", class = "sd")),
  control = list(adapt_delta = 0.99))
```

---

## Varying intercept + varying effect

```{r}
summ_varying_effect <- summary(varying_effect, probs = c(0.025,0.975))
summ_varying_effect
```
---


## Varying intercept + varying effect

```{r}
predict_data <- tibble(batch = batches) %>% crossing(group = unique(base_data$group))
s <- posterior_linpred(varying_effect, newdata = predict_data)
varying_effect_to_plot <- predict_data %>% mutate(log_mean = colMeans(s), mean = exp(log_mean), sd = sqrt(diag(cov(s))), low = log_mean - sd, high = log_mean + sd)


plot_multilevel_effect <- base_data %>% 
  base_plot() %>%
  add_error_bars() %>%
  add_my_crossbar(varying_effect_to_plot)  %>% add_my_point() + facet_wrap(~batch)
  

plot_multilevel_effect
```



```{r}
# # Partial pooling
# hide_axis_theme <- theme(axis.text.y = element_blank(), axis.title.y = element_blank(), axis.line.y = element_blank(),
#                          axis.ticks.y = element_blank())
# (((plot_all_together + hide_axis_theme)) + (plot_multilevel_effect + hide_axis_theme)) & expand_limits( y = c(-0.55, 0.5))) + plot_layout(widths = c(0.25,0.25,0.5))
```



---
class: inverse, center, middle

# Meta-analysis

---

## Now it's simple

That's exactly what meta-analysis is! 

```{r}
metaanalysis_data <- data.frame(estimate = c(-0.7, -1.5, -0.4, 0.1), std_error = c(1, 0.8, 0.4, 3)/3, study = LETTERS[1:4] ) %>% mutate(ci_low = estimate - 1.96 * std_error, ci_high = estimate + 1.96 * std_error)

plot_metaanalysis <- metaanalysis_data %>%
  ggplot(aes(x = study, y = estimate, ymin = ci_low, ymax = ci_high)) +
  geom_hline(yintercept = 0, color = "blue") + 
  geom_linerange() + 
  geom_point() + coord_flip() 

plot_metaanalysis

```


---

## Runing a meta-analysis 

```{r brm_meta, echo=TRUE, results="hide", cache=TRUE}
metaanalysis <- brms::brm(
  estimate | se(std_error) ~ 1 + (1 | study), 
  data = metaanalysis_data, 
  control = list(adapt_delta = 0.99))
```

---

```{r}
summ_meta <- summary(metaanalysis)
metaanalysis_data %>%
  select(-std_error) %>%
  rbind(data.frame(estimate = summ_meta$fixed[,"Estimate"], ci_low = summ_meta$fixed[,"l-95% CI"], ci_high = summ_meta$fixed[,"u-95% CI"], study = "META")) %>%
  ggplot(aes(x = study, y = estimate, ymin = ci_low, ymax = ci_high)) +
  geom_hline(yintercept = 0, color = "blue") + 
  geom_linerange() + 
  geom_point() + coord_flip() 

```


---


# Further reading

- How to develop a model without overfitting?
  - Bayesian Workflow https://arxiv.org/abs/2011.01808

## Thanks!