Data_simulation.R

## ----setup--------------------------------------------------
library(tidyverse)
library(kogpsy)



## ----DDM from Andrew----------------------------------------
out <- drift_diffusion(bias = 0.3, driftrate = 0.8)

ggplot2::ggplot(out, aes(time, dv)) +
  geom_line()


## ----my adaptations-----------------------------------------
my_drift_diffusion <- function(bias = 0.5,
                               driftrate = 1,
                               decision_boundary = 2,
                               ndt = 0.5,
                               diffvar = 0.1,
                               dt = 0.001) {
  assertthat::assert_that(diffvar > 0)

  bias <- as.integer(2 * decision_boundary * bias - decision_boundary)

  dv <- rnorm(1, mean = bias, sd = sqrt(dt))
  j <- 2

  while (abs(tail(dv,1)) < decision_boundary) {

      if (j <= ndt/dt) {
        dv[j] <- rnorm(1, mean = bias, sd = sqrt(dt))
      }
      else {
        error <- rnorm(1, 0, sqrt(diffvar * dt))
        dv[j] <- dv[j - 1] + driftrate * dt + error
        if (abs(dv[j]) > decision_boundary) {
          dv[j] <- ifelse(dv[j] > 0, min(dv[j], decision_boundary),
                          max(dv[j], -decision_boundary))
          (break)()
        }
      }
    j <- j+1
  }
  out <- dplyr::tibble(time = round(seq_along(dv) * dt, 2),
                       dv = dv, steps = seq_along(dv))
  return(out)
}


## ----data from one participant and one trial----------------
my_intercept = 0
my_b_congruence = 0.01
my_congruence = 1
my_b_rec = 0
my_rec = 0.3 # this will also change between participants according to my own data
my_b_interaction = -0.5 # this will be the thing that we change between participants. Also this should vary a bit within participants, so generate trial-noise

my_dr <- my_intercept + my_b_congruence*my_congruence + my_b_rec* my_rec + my_b_interaction * my_congruence * my_rec


out <- my_drift_diffusion(bias = 0.3, driftrate = my_dr, decision_boundary = 1)

ggplot2::ggplot(out, aes(time, dv)) +
  geom_line()

my_rt <- tail(out$time,1)



## ----runtrials of one Vp------------------------------------
n_trials <- 20



#this varies a bit within participants, but the means stay the same, because we assume that the effects are the same for all
my_intercept = rnorm(n = 1, mean = 0, sd = 0.1)
my_b_congruence = rnorm(n = 1, mean = 0.01, sd = 0.01)
my_congruence = 1
my_b_rec = rnorm(n = 1, mean = 0, sd = 0.01)

true_dr <-  my_intercept + my_b_congruence*my_congruence + my_b_rec* my_rec + my_b_interaction * my_congruence * my_rec
error <- abs(true_dr*0.1)

# these are the differences between participants
my_rec = rnorm(n = 1, 0.3, sd = 0.02) # this will  change between participants according to my own data
my_b_interaction = -0.2 # this will be the thing that we change between participants. Also this should vary a bit within participants, so generate trial-noise


# set up the 20 trials
rts <- rep(NA, n_trials)
for (i in 1:n_trials) {
  my_dr <- rnorm(mean = true_dr, sd = error, n = 1)
  tmpout <- my_drift_diffusion(bias = 0.3, driftrate = my_dr, decision_boundary = 1)
  rts[i] <- tail(tmpout$time, 1)
}

hist(rts)



## ----trials for some participants---------------------------
n_vp <- 100
n_trials <- 20

# prepare the output-file
result_file <- tibble(Vp = rep(1:n_vp, each = n_trials),
                      trialNum = rep(1:n_trials, n_vp),
                      congruence = rep(NA, n_vp*n_trials),
                      rec = rep(NA, n_vp*n_trials),
                      true_dr = rep(NA, n_vp*n_trials),
                      true_interaction_beta = rep(NA, n_vp*n_trials),
                      rt = rep(NA, n_vp*n_trials),
                      answer = rep(NA, n_vp*n_trials)
                      )

# vp-loop
line_counter <- 1 # to index the line in results_file

for (corrVp in 1:n_vp) {
  # set the participant-specific parameters: rec and b3
  # recurrence: sample a rec value from a distirbution that represents the rec in old experiments

   #rbeta(1000,0.6, 2) #this reflects best the distribution of the true data in the cognition papaer

  curr_true_rec <- rbeta(1, 0.6, 2)

  # interaction: sample based on djikstra: half of participants show no influence, twice as many show facilitation as show inhibition. a normal distribution with a slightly positive mean will do.
  curr_b_interaction = rnorm(1, mean = 0.2, sd = 0.3)

  #hist(rnorm(1000, mean = 0.2, sd = 0.3))

  for (currTrial in 1:n_trials) {

    #include the trial-by-trial variance
    #this varies a bit within participants, but the means stay the same, because we assume that the effects are the same for all
    my_intercept = rnorm(n = 1, mean = 0, sd = 0.1)

   my_b_congruence = rnorm(n = 1, mean = 0.1, sd = 0.01) # sligtly positive because there is an overall trend

    my_b_rec = rnorm(n = 1, mean = 0, sd = 0.01) # we don't expect a general effect

    if (currTrial <= 10) { #because in half of the trials there will be a congruent image
      my_congruence = 1 }
    else {
      my_congruence = 0}

    # set the true drift rate
    true_dr <-  my_intercept +
                my_b_congruence * my_congruence +
                my_b_rec* curr_true_rec +
                curr_b_interaction * my_congruence * curr_true_rec

    #error <- abs(true_dr*0.01)
    my_dr <- rnorm(mean = true_dr, sd = 0.01, n = 1)

    currOut <- my_drift_diffusion(bias = 0.05, driftrate = my_dr, decision_boundary = 1) # rtdists rdiffusion

    currAnswer <- ifelse(tail(currOut$dv,1) > 0, 1, -1)
    currRT <- tail(currOut$time, 1)

    result_file$congruence[line_counter] <- my_congruence
    result_file$rec[line_counter] <- curr_true_rec
    result_file$true_dr[line_counter] <- true_dr
    result_file$rt[line_counter] <-  abs(currRT)
    result_file$answer[line_counter] <- currAnswer
    result_file$true_interaction_beta[line_counter] <- curr_b_interaction

  line_counter <- line_counter + 1
  } # trial loop
} #participant loop

result_file <- result_file %>%
  mutate(congruence = as.factor(congruence))
# look at the data
result_file %>% ggplot(aes(x = rt, y = ..density.., color = congruence)) + geom_density()

#fit a brms model to recover the parameters
library(brms)

set_prior(prior = "uniform(min = -2, max", group= "Vp", coef = "congruence1:rec")

myP <- prior(student_t(3, 0.2, 8), class = sd, coef = congruence1:rec, group = Vp)

firstfit <- brm(rt|dec(answer) ~ 1 + congruence + rec + congruence*rec +
                  (1 + congruence + rec + congruence*rec | Vp),
                data = result_file,
                family = wiener(),
                prior = myP,
                iter = 2000,
                chains = 3,
                cores = 3)

prior_summary(firstfit)

conditional_effects(firstfit)
# now, recover the parameters for every individual and compare them to the real parameters (do the posterior distributions include the true b3 interaction? )

a <- ranef(firstfit)
# a$Vp[,1,"congruence:rec"] # estimates
# a$Vp[,3,"congruence:rec"] # Q2.5
# a$Vp[,4,"congruence:rec"] # Q97.5
test_result_file <- result_file %>%
  group_by(Vp) %>%
  summarize(true_b3 = unique(true_interaction_beta)) %>%
  mutate(estim_b3 = a$Vp[,1,"congruence1:rec"], # estimates,
         estim_b3Q2.5 = a$Vp[,3,"congruence1:rec"], # estimates
         estim_b3Q97.5 = a$Vp[,4,"congruence1:rec"])


test_result_file %>%
  ggplot(aes(x = Vp)) +
  geom_point(aes(y = true_b3)) +
  geom_point(aes(y = estim_b3), color = "red") +
  geom_errorbar(aes(ymin = estim_b3Q2.5,
                    y = estim_b3,
                    ymax = estim_b3Q97.5), alpha = 0.2)


test_result_file %>%
  gather(key = origin,
         value = b3size,
         2:3) %>%
  ggplot(aes(x = b3size, y= ..density.., color = origin)) + geom_density()



## ----loop over nVpn-----------------------------------------