Simulated_contest

library(tidyverse)
library(brms)
library(patchwork)

## Define amount of data
items <- 30
chains <- 5
raters <- 30

## Define parameter values
GenerationB <- .4
sigma_raters <- .4
sigma_items <- .4
sigma_chains <- .2


# Create dataset with stimuli, with generation and chain ids
d_items <- tibble(
  item = 1:items,
  item_intercept = rnorm(items, 0, sigma_items),
  generation = rep(-1:1, 10),
  chain = rep(seq(chains), each = items/chains),
  chain_slope = NA
)

chain_slope = rnorm(chains, 0, sigma_chains)

for (i in seq(items)) {
  d_items$chain_slope[i] <- chain_slope[d_items$chain[i]]
}

# Create dataset with raters
d_raters <- tibble(
  ID = 1:raters,
  rater_slope = rnorm(raters, 0, sigma_raters)
)


## Now create all possible combinations between the stimuli
combinations <- t(combn(1:items, 2))

df_recovery <- NULL
for (rater in seq(raters)) {
  print(rater)
  temp <- tibble(
    stimulusR = combinations[,1],
    stimulusL = combinations[,2],
    chainR = NA,
    chainL = NA,
    generationR = NA,
    generationL = NA,
    intR = NA,
    intL = NA
  )
  
  for (i in seq(nrow(combinations))) {
    
    temp$chainR[i] <- d_items$chain[temp$stimulusR[i]]
    temp$chainL[i] <- d_items$chain[temp$stimulusL[i]]
    temp$generationR[i] <- d_items$generation[temp$stimulusR[i]]
    temp$generationL[i] <- d_items$generation[temp$stimulusL[i]]
        
    temp$intR[i] <- d_items$item_intercept[temp$stimulusR[i]] + 
                         (GenerationB + d_items$chain_slope[temp$stimulusR[i]] + 
                            d_raters$rater_slope[rater]) * d_items$generation[temp$stimulusR[i]]
    
    temp$intL[i] <- d_items$item_intercept[temp$stimulusL[i]] + 
                         (GenerationB + d_items$chain_slope[temp$stimulusL[i]] + 
                            d_raters$rater_slope[rater]) * d_items$generation[temp$stimulusL[i]]
    
    temp$diff[i] = temp$intR[i] - temp$intL[i]
    temp$thetaR[i] = inv_logit_scaled(temp$diff[i])
    temp$choiceR[i] = rbinom(1, 1, temp$thetaR[i])
    temp$rater[i] = rater
  }
  
  if (exists("df_recovery")) {df_recovery <- rbind(df_recovery, temp)} else {df_recovery <- temp}
  
}

p1 <- df_recovery %>% 
  group_by(generationL, generationR, stimulusR) %>% 
  summarize(x = mean(choiceR)) %>%
  mutate(generationR = as.numeric(generationR), generationL = as.numeric(generationL)) %>%
  ggplot(aes(generationR, x)) + 
  geom_point() + geom_smooth(aes(generationR, x)) + 
  facet_wrap(.~ generationL) + theme_bw() + theme(legend.position = "none")

p2 <- df_recovery %>% 
  group_by(generationL, generationR, stimulusR, rater) %>% 
  summarize(x = mean(choiceR)) %>%
  mutate(generationR = as.numeric(generationR), generationL = as.numeric(generationL)) %>%
  ggplot(aes(generationR, x, group = rater, color = rater)) + 
  geom_point() + geom_smooth(aes(generationR, x), method = lm, se = F) + 
  facet_wrap(.~ generationL) + theme_bw() + theme(legend.position = "none")

p3 <- df_recovery %>% 
  group_by(generationL, generationR, stimulusR, chainR) %>% 
  summarize(x = mean(choiceR)) %>%
  mutate(generationR = as.numeric(generationR), generationL = as.numeric(generationL)) %>%
  ggplot(aes(generationR, x, group = chainR, color = chainR)) + 
  geom_point() + geom_smooth(aes(generationR, x), method = lm, se = F) + 
  facet_wrap(.~ generationL) + theme_bw() + theme(legend.position = "none")


p1/p2/p3

write_csv(df_recovery, here("data", "exp2_simulated_data.csv"))


d1_sim <- list(
  N_choices = nrow(df_recovery),
  N_items = items,
  N_chains = chains,
  N_raters = raters,
  y = as.numeric(df_recovery$choiceR),
  item1 = as.numeric(df_recovery$stimulusR),
  item2 = as.numeric(df_recovery$stimulusL),
  rater = as.numeric(df_recovery$rater),
  generation = d_items$generation, 
  chain = d_items$chain  
)

file <- file.path(here("scripts", "Exp2_ContestModel.stan"))
mod <- cmdstan_model(file, cpp_options = list(stan_threads = TRUE), pedantic = TRUE)

# The following command calls Stan with specific options.
samples <- mod$sample(
  data = d1_sim,
  seed = 123,
  chains = 2,
  parallel_chains = 2,
  threads_per_chain = 2,
  iter_warmup = 2000,
  iter_sampling = 2000,
  refresh = 50,
  max_treedepth = 20,
  adapt_delta = 0.99,
)

samples$save_object(file = here("models", "Exp2_simulated.RDS"))