bayesian_p2_ANOVA.Rmd

---
title: "Bayesian effect size for ANOVA designs"
author: "Fabian Dablander, Maarten Marsman"
date: "May 22, 2016"
output: html_document
---

## 1) Example: One-way design
Download the data from [https://osf.io/zcipy/](https://osf.io/zcipy/), save it under
"example1.sav" in the current directory.

```{r, message = FALSE, warning = FALSE}
library('haven')
library('dplyr')
library('BayesFactor')
NITER = 100000

dat = read_spss('example1.sav') %>% 
  filter(ThrowOut == 2, !is.na(Accuracy)) %>% 
  select(Accuracy, IdentityCorrectSalience) %>% 
  rename(y = Accuracy, x = IdentityCorrectSalience) %>% 
  mutate(x = factor(x, levels = 1:3, labels = c('Asian', 'Control', 'Female')))
  # note that control is number 2, not 3 as in the manuscript

head(dat)
```

Convenience and squared multiple correlation function.
```{r}
# gets the posterior of beta weights for each factor level
get_posterior = function(post, levls) {
  cnames = colnames(post)
  select = sapply(levls, function(lv) grep(lv, cnames))
  post[, select]
}

# don't call this directly
.partial2 = function(betas, psig, w) {
  # easier to read than matrix operations
  beta = sum(w * betas)
  num = sum(w * (betas - beta)^2)
  num / (num + psig)
}

# calls the p2 function for each sample of the posterior
partial2 = function(betas, psig, w) {
  n1 = ncol(betas)
  joint = cbind(betas, psig)
  
  apply_partial2 = function(vec) {
    .partial2(vec[seq_len(n1)], vec[n1 + 1], w)
  }
  
  unname(apply(joint, 1, apply_partial2))
}
```


### 1a) Estimation and effect size (unrestricted)
```{r}
bf = anovaBF(y ~ x, data = data.frame(dat))
post = posterior(bf, iterations = 2*NITER)
apply(post, 2, mean)
```

```{r}
levls = levels(dat$x)

psig = get_posterior(post, 'sig2')
betas = get_posterior(post, levls)
weights = table(dat$x) / sum(table(dat$x))

p2 = partial2(betas, psig, weights)
mean(p2)
```

### 1b) Estimation and effect size (restricted)
```{r}
select = betas[, 3] < betas[, 2] & betas[, 2] < betas[, 1]
restr_post = post[select, ]
apply(restr_post, 2, mean)
```

```{r}
restr_psig = get_posterior(restr_post, 'sig2')
restr_betas = get_posterior(restr_post, levls)

restr_p2 = partial2(restr_betas, restr_psig, weights)
mean(restr_p2)
```

```{r}
hist(p2, breaks = 40, main = 'unrestricted')
hist(restr_p2, breaks = 40, main = 'restricted')
```

## 2) Example: 2 x 3 mixed design
Download the data from [https://osf.io/4qejs/](https://osf.io/4qejs/), save it under
"example2.sav" in the current directory.

```{r, warning = FALSE, message = FALSE}
library('tidyr')

dat = read_spss('example2.sav') %>%
  filter(Orig_rep == 2) %>%
  select(-Orig_rep, -mean_judge)

dat = mutate(dat, ID = 1:nrow(dat)) %>% 
  gather(., condition, y, incest, dog, flag) %>% 
  arrange(y, condition, time) %>% 
  mutate(y = as.numeric(y),
         condition = as.factor(condition),
         time = as.factor(time), ID = as.factor(ID))

head(dat)
```

Functions for squared semipartial correlations (pchange).
```{r}
# don't call this directly
.p2change = function(betas1, betas2, psig, tau, w1, w2) {
  beta1 = sum(w1 * betas1)
  beta2 = sum(w2 * betas2)
  num1 = sum(w1 * (betas1 - beta1)^2)
  num2 = sum(w2 * (betas2 - beta2)^2)
  
  denom = num1 + num2 + psig + tau
  c(num1, num2, tau) / denom
}

p2change = function(betas1, betas2, psig, tau, w1, w2) {
  n1 = ncol(betas1)
  n2 = ncol(betas2)
  joint = cbind(betas1, betas2, psig, tau)
  
  apply_p2change = function(vec) {
    .p2change(vec[seq_len(n1)], vec[seq_len(n2) + n1],
              vec[length(vec) - 1], vec[length(vec)], w1, w2)
  }
  
  unname(t(apply(joint, 1, apply_p2change)))
}
```

### Estimation and effect size
```{r}
bf = lmBF(y ~ time + condition + ID, whichRandom = 'ID', data = data.frame(dat))
post = posterior(bf, iterations = NITER)

levls2 = levels(dat$condition)
levls1 = paste0('time-', levels(dat$time))

psig = get_posterior(post, 'sig2')
betas1 = get_posterior(post, levls1)
betas2 = get_posterior(post, levls2)

# there should be a better way to get the random effects ..
prandom = post[, grepl('^ID-*', colnames(post))]
tau = apply(prandom, 1, var)
apply(cbind(betas1, betas2, psig, tau), 2, mean)
```

```{r}
w1 = with(dat, tapply(y, time, length))
w1 = w1 / sum(w1)
w2 = with(dat, tapply(y, condition, length))
w2 = w2 / sum(w2)

sp2 = p2change(betas1, betas2, psig, tau, w1, w2)
apply(sp2, 2, mean)
```

```{r}
hist(sp2[, 1], main = 'time (x1)', xlab = '',
     breaks = 40, xlim = c(0, 1))
hist(sp2[, 1] + sp2[, 2], main = 'time (x1) + moral (x2)',
     xlab = '', breaks = 40, xlim = c(0, 1))
hist(sp2[, 1] + sp2[, 2] + sp2[, 3], main = 'time (x1) + moral (x2) + ID (x3)',
     xlab = '', breaks = 40, xlim = c(0, 1))
```