Add hddm scripts

Signed-off-by: Liang Zhang <[email protected]>

DESCRIPTION

@@ -13,6 +13,7 @@ Depends:
 Imports:
     Amelia,
     arrow,
+    brms,
     clustermq,
     corrplot,
     cpmr,

_scripts_hddm/compile_models.R

@@ -0,0 +1,66 @@
+library(tidyverse)
+library(brms)
+# retest models
+## congruency effect/switch cost
+brm(
+  bf(
+    rt | dec(acc) ~ 0 + ocassion:type + (0 + ocassion:type | p | user_id),
+    bs ~ 0 + ocassion + (0 + ocassion | p | user_id),
+    ndt ~ 0 + ocassion + (0 + ocassion | p | user_id),
+    bias = 0.5
+  ),
+  withr::with_seed(
+    1,
+    expand_grid(
+      user_id = 0,
+      ocassion = c("test", "retest"),
+      type = c("inc", "con") # could be changed to switch related levels
+    ) |>
+      mutate(
+        rt = rexp(n()),
+        acc = sample(c(0, 1), n(), replace = TRUE)
+      )
+  ),
+  family = wiener(
+    link_bs = "identity",
+    link_ndt = "identity"
+  ),
+  prior = c(
+    set_prior("cauchy(0, 5)", class = "b"),
+    set_prior("gamma(3, .5)", class = "b", dpar = "bs", lb = 0),
+    set_prior("gamma(1, .5)", class = "b", dpar = "ndt", lb = 0)
+  ),
+  file = "data/model_diff_retest",
+  chains = 0
+)
+## simple intercept models
+brm(
+  bf(
+    rt | dec(acc) ~ 0 + ocassion + (0 + ocassion | p | user_id),
+    bs ~ 0 + ocassion + (0 + ocassion | p | user_id),
+    ndt ~ 0 + ocassion + (0 + ocassion | p | user_id),
+    bias = 0.5
+  ),
+  withr::with_seed(
+    1,
+    expand_grid(
+      user_id = 0,
+      ocassion = c("test", "retest")
+    ) |>
+      mutate(
+        rt = rexp(n()),
+        acc = sample(c(0, 1), n(), replace = TRUE)
+      )
+  ),
+  family = wiener(
+    link_bs = "identity",
+    link_ndt = "identity"
+  ),
+  prior = c(
+    set_prior("cauchy(0, 5)", class = "b"),
+    set_prior("gamma(3, .5)", class = "b", dpar = "bs", lb = 0),
+    set_prior("gamma(1, .5)", class = "b", dpar = "ndt", lb = 0)
+  ),
+  file = "data/model_simple_retest",
+  chains = 0
+)

_scripts_hddm/hddm.qsub

@@ -0,0 +1,10 @@
+#!/usr/bin/env bash
+#$ -cwd
+#$ -S /bin/bash
+#$ -N hddm
+#$ -o logs/qsub
+#$ -q withpe.q
+#$ -pe ompi 4
+#$ -j y
+
+Rscript $SCRIPT

_scripts_hddm/retest_alert.R

@@ -0,0 +1,62 @@
+suppressPackageStartupMessages({
+  library(tidyverse)
+  library(brms)
+})
+requireNamespace("bit64", quietly = TRUE)
+projects <- targets::tar_config_yaml()
+game_id <- "380173315257221"
+file_save <- str_glue("data-raw/hddm-retest/game-{game_id}_type-cong_wiener")
+cat("Processing: ", data.iquizoo::match_info(game_id, "game_name"), "\n\n")
+data <- targets::tar_read_raw(
+  paste0("data_valid_", game_id),
+  store = projects$prepare_source_data_retest$store
+) |>
+  filter(n() == 2, .by = user_id) |>
+  mutate(
+    ocassion = if_else(row_number(game_time) == 1, "test", "retest"),
+    .by = user_id
+  ) |>
+  unnest(raw_parsed) |>
+  mutate(
+    RT = RT / 1000,
+    # unify levels
+    StimType = case_match(
+      StimType,
+      "congruent" ~ "con",
+      "incongruent" ~ "inc"
+    )
+  ) |>
+  filter(RT > 0.25) |>
+  select(user_id, ocassion, type = StimType, acc = ACC, rt = RT)
+model <- readRDS("data/model_diff_retest.rds")
+
+if (interactive()) {
+  View(standata(update(model, newdata = data, chains = 0)))
+}
+
+initfun <- function() {
+  K <- 4 # 2 by 2
+  K_bs <- 2
+  K_ndt <- 2
+  M_1 <- K + K_bs + K_ndt
+  N_1 <- 189 # number of subjects
+  list(
+    b = as.array(rnorm(K)),
+    b_bs = as.array(runif(K_bs, 1, 2)),
+    b_ndt = as.array(runif(K_ndt, 0.05, 0.1)),
+    sd_1 = as.array(runif(M_1, 0.5, 1)),
+    z_1 = matrix(rnorm(M_1 * N_1, 0, 0.01), M_1, N_1),
+    L_1 = diag(M_1)
+  )
+}
+fit_wiener <- update(
+  model,
+  newdata = data,
+  init = initfun,
+  file = file_save,
+  iter = 1000,
+  warmup = 500,
+  chains = 4,
+  cores = 4,
+  control = list(max_treedepth = 15, adapt_delta = 0.9)
+)

_scripts_hddm/retest_antisac.R

@@ -0,0 +1,54 @@
+suppressPackageStartupMessages({
+  library(tidyverse)
+  library(brms)
+})
+requireNamespace("bit64", quietly = TRUE)
+projects <- targets::tar_config_yaml()
+game_id <- "238239294447813"
+file_save <- str_glue("data-raw/hddm-retest/game-{game_id}_type-simple_wiener")
+cat("Processing: ", data.iquizoo::match_info(game_id, "game_name"), "\n\n")
+data <- targets::tar_read_raw(
+  paste0("data_valid_", game_id),
+  store = projects$prepare_source_data_retest$store
+) |>
+  filter(n() == 2, .by = user_id) |>
+  mutate(
+    ocassion = if_else(row_number(game_time) == 1, "test", "retest"),
+    .by = user_id
+  ) |>
+  unnest(raw_parsed) |>
+  mutate(RT = RT / 1000) |>
+  filter(RT > 0.25 & RT < 20) |>
+  select(user_id, ocassion, acc = ACC, rt = RT)
+model <- readRDS("data/model_simple_retest.rds")
+
+if (interactive()) {
+  View(standata(update(model, newdata = data, chains = 0)))
+}
+
+initfun <- function() {
+  K <- 2 # 2 by 2
+  K_bs <- 2
+  K_ndt <- 2
+  M_1 <- K + K_bs + K_ndt
+  N_1 <- 194 # number of subjects
+  list(
+    b = as.array(rnorm(K)),
+    b_bs = as.array(runif(K_bs, 1, 2)),
+    b_ndt = as.array(runif(K_ndt, 0.05, 0.1)),
+    sd_1 = as.array(runif(M_1, 0.5, 1)),
+    z_1 = matrix(rnorm(M_1 * N_1, 0, 0.01), M_1, N_1),
+    L_1 = diag(M_1)
+  )
+}
+fit_wiener <- update(
+  model,
+  newdata = data,
+  init = initfun,
+  file = file_save,
+  iter = 1000,
+  warmup = 500,
+  chains = 4,
+  cores = 4,
+  control = list(max_treedepth = 15, adapt_delta = 0.9)
+)

_scripts_hddm/retest_birds_cong.R

@@ -0,0 +1,62 @@
+suppressPackageStartupMessages({
+  library(tidyverse)
+  library(brms)
+})
+requireNamespace("bit64", quietly = TRUE)
+projects <- targets::tar_config_yaml()
+game_id <- "224377200398533"
+file_save <- str_glue("data-raw/hddm-retest/game-{game_id}_type-cong_wiener")
+cat("Processing: ", data.iquizoo::match_info(game_id, "game_name"), "\n\n")
+data <- targets::tar_read_raw(
+  paste0("data_valid_", game_id),
+  store = projects$prepare_source_data_retest$store
+) |>
+  filter(n() == 2, .by = user_id) |>
+  mutate(
+    ocassion = if_else(row_number(game_time) == 1, "test", "retest"),
+    .by = user_id
+  ) |>
+  unnest(raw_parsed) |>
+  mutate(
+    RT = RT / 1000,
+    # unify levels
+    StimType = case_match(
+      StimType,
+      "Congruent" ~ "con",
+      "Incongruent" ~ "inc"
+    )
+  ) |>
+  filter(RT > 0.15) |>
+  select(user_id, ocassion, type = StimType, acc = ACC, rt = RT)
+model <- readRDS("data/model_diff_retest.rds")
+
+if (interactive()) {
+  View(standata(update(model, newdata = data, chains = 0)))
+}
+
+initfun <- function() {
+  K <- 4 # 2 by 2
+  K_bs <- 2
+  K_ndt <- 2
+  M_1 <- K + K_bs + K_ndt
+  N_1 <- 383 # number of subjects
+  list(
+    b = as.array(rnorm(K)),
+    b_bs = as.array(runif(K_bs, 1, 2)),
+    b_ndt = as.array(runif(K_ndt, 0.05, 0.1)),
+    sd_1 = as.array(runif(M_1, 0.5, 1)),
+    z_1 = matrix(rnorm(M_1 * N_1, 0, 0.01), M_1, N_1),
+    L_1 = diag(M_1)
+  )
+}
+fit_wiener <- update(
+  model,
+  newdata = data,
+  init = initfun,
+  file = file_save,
+  iter = 1000,
+  warmup = 500,
+  chains = 4,
+  cores = 4,
+  control = list(max_treedepth = 15, adapt_delta = 0.9)
+)

_scripts_hddm/retest_birds_switch.R

@@ -0,0 +1,62 @@
+suppressPackageStartupMessages({
+  library(tidyverse)
+  library(brms)
+})
+requireNamespace("bit64", quietly = TRUE)
+projects <- targets::tar_config_yaml()
+game_id <- "224377200398533"
+file_save <- str_glue("data-raw/hddm-retest/game-{game_id}_type-switch_wiener")
+cat("Processing: ", data.iquizoo::match_info(game_id, "game_name"), "\n\n")
+data <- targets::tar_read_raw(
+  paste0("data_valid_", game_id),
+  store = projects$prepare_source_data_retest$store
+) |>
+  filter(n() == 2, .by = user_id) |>
+  mutate(
+    ocassion = if_else(row_number(game_time) == 1, "test", "retest"),
+    .by = user_id
+  ) |>
+  unnest(raw_parsed) |>
+  mutate(
+    RT = RT / 1000,
+    # unify levels
+    TaskType = case_match(
+      TaskType,
+      "Switch" ~ "switch",
+      "Repeat" ~ "repeat"
+    )
+  ) |>
+  filter(RT > 0.15, !is.na(TaskType)) |>
+  select(user_id, ocassion, type = TaskType, acc = ACC, rt = RT)
+model <- readRDS("data/model_diff_retest.rds")
+
+if (interactive()) {
+  View(standata(update(model, newdata = data, chains = 0)))
+}
+
+initfun <- function() {
+  K <- 4 # 2 by 2
+  K_bs <- 2
+  K_ndt <- 2
+  M_1 <- K + K_bs + K_ndt
+  N_1 <- 383 # number of subjects
+  list(
+    b = as.array(rnorm(K)),
+    b_bs = as.array(runif(K_bs, 1, 2)),
+    b_ndt = as.array(runif(K_ndt, 0.05, 0.1)),
+    sd_1 = as.array(runif(M_1, 0.5, 1)),
+    z_1 = matrix(rnorm(M_1 * N_1, 0, 0.01), M_1, N_1),
+    L_1 = diag(M_1)
+  )
+}
+fit_wiener <- update(
+  model,
+  newdata = data,
+  init = initfun,
+  file = file_save,
+  iter = 1000,
+  warmup = 500,
+  chains = 4,
+  cores = 4,
+  control = list(max_treedepth = 15, adapt_delta = 0.9)
+)

_scripts_hddm/retest_cardsort.R

@@ -0,0 +1,62 @@
+suppressPackageStartupMessages({
+  library(tidyverse)
+  library(brms)
+})
+requireNamespace("bit64", quietly = TRUE)
+projects <- targets::tar_config_yaml()
+game_id <- "224377872388293"
+file_save <- str_glue("data-raw/hddm-retest/game-{game_id}_type-switch_wiener")
+cat("Processing: ", data.iquizoo::match_info(game_id, "game_name"), "\n\n")
+data <- targets::tar_read_raw(
+  paste0("data_valid_", game_id),
+  store = projects$prepare_source_data_retest$store
+) |>
+  filter(n() == 2, .by = user_id) |>
+  mutate(
+    ocassion = if_else(row_number(game_time) == 1, "test", "retest"),
+    .by = user_id
+  ) |>
+  unnest(raw_parsed) |>
+  mutate(
+    RT = RT / 1000,
+    # unify levels
+    Type = case_match(
+      Type,
+      "Switch" ~ "switch",
+      "Repeat" ~ "repeat"
+    )
+  ) |>
+  filter(RT > 0.15, !is.na(Type)) |>
+  select(user_id, ocassion, type = Type, acc = ACC, rt = RT)
+model <- readRDS("data/model_diff_retest.rds")
+
+if (interactive()) {
+  View(standata(update(model, newdata = data, chains = 0)))
+}
+
+initfun <- function() {
+  K <- 4 # 2 by 2
+  K_bs <- 2
+  K_ndt <- 2
+  M_1 <- K + K_bs + K_ndt
+  N_1 <- 183 # number of subjects
+  list(
+    b = as.array(rnorm(K)),
+    b_bs = as.array(runif(K_bs, 1, 2)),
+    b_ndt = as.array(runif(K_ndt, 0.05, 0.1)),
+    sd_1 = as.array(runif(M_1, 0.5, 1)),
+    z_1 = matrix(rnorm(M_1 * N_1, 0, 0.01), M_1, N_1),
+    L_1 = diag(M_1)
+  )
+}
+fit_wiener <- update(
+  model,
+  newdata = data,
+  init = initfun,
+  file = file_save,
+  iter = 1000,
+  warmup = 500,
+  chains = 4,
+  cores = 4,
+  control = list(max_treedepth = 15, adapt_delta = 0.9)
+)