### This script calculates the probability of a chimpanzee pulling the left lever
### Using rstanarm and leave one out cv (on actor level)

# Import libraries and data
library(rethinking)
library(rstanarm)
library(loo)
options(mc.cores = parallel::detectCores())
data(chimpanzees)
d <- chimpanzees

# How many actors?
nractors <- length(unique(d$actor))

### Define null model
m0 <- stan_glmer(pulled_left~+(1|actor),
                 data=d, family=binomial)

# Define model with one predictor
m1 <- update(m0,.~.+condition)

### Compare using "classic loo"
loo_m0 <- loo(m0)
loo_m1 <- loo(m1)
compare(loo_m0,loo_m1)

### Compare using grouped k-fold cv
folds <- kfold_split_stratified(K = nractors, x = d$actor)
kf_m0 <- kfold(m0, K = nractors, folds = folds)
kf_m1 <- kfold(m1, K = nractors, folds = folds)
compare(kf_m0,kf_m1)

### Compare respecting grouping structure
m0par<- matrix(NA,nractors,1)
m1par<- matrix(NA,nractors,1)
for(i in 1:nractors){
  m0par[i] <- sum(kf_m0$pointwise[(folds==i),])
  m1par[i] <- sum(kf_m1$pointwise[(folds==i),])
}

round(sum(m1par-m0par),2) # Diff
round(sd(m1par-m0par),2) # se