assignment-4.Rmd

---
title: "Assignment 4"
description: | 
  Effect of tDCS stimulation on episodic memory
date: 05-28-2021
author:
  - first_name: "Andrew"
    last_name: "Ellis"
    url: https://github.com/awellis
    affiliation: Kognitive Psychologie, Wahrnehmung und Methodenlehre, Universität Bern 
    affiliation_url: https://www.kog.psy.unibe.ch
    orcid_id: 0000-0002-2788-936X

citation_url: https://awellis.github.io/learnmultilevelmodels/asssignment-4.html
bibliography: ./bibliography.bib
output: 
    distill::distill_article:
      toc: true
      toc_float: true
      toc_depth: 2
      code_folding: true
      self_contained: false
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```


```{r}
library(tidyverse)
library(brms)

theme_set(theme_grey(base_size = 14) +
            theme(panel.grid = element_blank()))
```




# Effect of tDCS on memory


You are going to analyze data from a study in which 34 elderly subjects were given anodal (activating) tDCS over their left tempero-parietal junction (TPJ). Subjects were instructed to learn association between images and pseudo-words (data were inspired by @antonenkoTDCSinducedEpisodicMemory2019).

Episodic memory is measured by percentage of correctly recalled word-image pairings. We also have response times for correct decisions.

Each subject was tested 5 times in the TPJ stimulation condition, and a further 5 times in a sham stimulation condition.

The variables in the dataset are:

```
subject: subject ID
stimulation: TPJ or sham (control)
block: block
age: age
correct: accuracy per block
rt: mean RTs for correct responses
```

You are mainly interested in whether recall ability is better during the TPJ stimulation condition than during sham.


```{r}
d <- read_csv("https://raw.githubusercontent.com/kogpsy/neuroscicomplab/main/data/tdcs-tpj.csv")
```


```{r}
glimpse(d)
```

```{r}
d <- d %>% 
  mutate(across(c(subject, stimulation, block), ~as_factor(.))) %>% 
  drop_na()
```



:::exercise

1) specify a linear model

2) specify a null model

3) compare models using approximate leave-one-out cross-validation

4) estimate a Bayes factor for the effect of TPJ stimulation
    + using the Savage-Dickey density ratio test
    + using the bridge sampling method
    
5) control for subjects' age
:::


## Linear model

```{r}
fit1 <- brm(correct ~ stimulation + (stimulation| subject),
            prior = prior(normal(0, 1), class = b),
            data = d,
            file = "models/ass4-1")
```


```{r}
summary(fit1)
```


```{r}
mcmc_plot(fit1, "b_", type = "areas")
```

## Null model

```{r}
fit2 <- brm(correct ~ 1 + (stimulation| subject),
            data = d,
            file = "models/ass4-2")
```




## Model comparison using loo

```{r}
loo_fit_1 <- loo(fit1)
loo_fit_1
```

```{r}
loo_fit_2 <- loo(fit2)
loo_fit_2
```

```{r}
loo_compare(loo_fit_1, loo_fit_2)
```

## Model comparison via Bayes factor

```{r}
fit1_bf <- brm(correct ~ stimulation + (stimulation| subject),
            prior = prior(normal(0, 1), class = b),
            data = d,
            iter = 1e4,
            save_pars = save_pars(all = TRUE),
            file = "models/ass4-bf_1")
```


```{r}
fit2_bf <- brm(correct ~ 1 + (stimulation| subject),
            prior = prior(normal(0, 1), class = b),
            data = d,
            iter = 1e4,
            save_pars = save_pars(all = TRUE),
            file = "models/ass4-bf_2")
```



```{r}
loglik_1 <- bridge_sampler(fit1_bf, silent = TRUE)
loglik_2 <- bridge_sampler(fit2_bf, silent = TRUE)
```

 
 
```{r}
BF_10 <- bayes_factor(loglik_1, loglik_2)
BF_10
```
 

## Control for age

```{r}
fit3 <- brm(correct ~ stimulation + age + (stimulation | subject),
            prior = prior(normal(0, 1), class = b),
            data = d,
            file = "models/ass4-3")
```


```{r}
fit3
```


```{r}
loo_fit_3 <- loo(fit3)
```


```{r}
loo_compare(loo_fit_1, loo_fit_3)
```



```{r}
fit3_bf <- brm(correct ~ stimulation + age + (stimulation| subject),
            prior = prior(normal(0, 1), class = b),
            data = d,
            iter = 1e4,
            save_pars = save_pars(all = TRUE),
            file = "models/ass4-bf_3")
```

```{r}
loglik_3 <- bridge_sampler(fit3_bf, silent = TRUE)
```



```{r}
BF_01 <- bayes_factor(loglik_3, loglik_1)
BF_01
```