update data

study1/analysis/1_cleaning.qmd

@@ -0,0 +1,193 @@
+---
+title: "Doggo/Nogo (Study 1) - Data Cleaning"
+editor: source
+editor_options: 
+  chunk_output_type: console
+format:
+  html:
+    code-fold: true
+    self-contained: true
+---
+
+## Data Preparation
+
+```{r}
+#| message: false
+#| warning: false
+
+library(tidyverse)
+library(easystats)
+library(patchwork)
+library(ggside)
+```
+
+
+```{r}
+#| code-fold: false
+
+dfsub <- read.csv("../data/rawdata_participants.csv") 
+df <- read.csv("../data/rawdata_game.csv") 
+```
+
+The initial sample consisted of `r report::report_participants(df, age="Age", gender="Gender", education="Education")`.
+
+
+```{r}
+# Feedback computation
+newdf <- data.frame()
+for(i in unique(df$Participant)) {
+  for(s in unique(df$Session)) {
+    dat <- df[df$Participant == i & df$Session == s, ]
+    dat$Feedback <- c("Positive", ifelse(dat$RT <= dat$Threshold, "Positive", "Negative")[2:nrow(dat)])
+  
+    n_positive <- c(0)
+    n_negative <- c(0)
+    counter_positive <- 0
+    counter_negative <- 0
+    for(j in 2:nrow(dat)) {
+      prev_feedback <- dat$Feedback[j-1]
+      current_feedback <- dat$Feedback[j]
+  
+      if(prev_feedback == "Positive") {
+        counter_positive <- counter_positive + 1
+      } else {
+        counter_positive <- 0
+      }
+      n_positive <- c(n_positive, counter_positive)
+  
+      if(prev_feedback == "Negative") {
+        counter_negative <- counter_negative + 1
+      } else {
+        counter_negative <- 0
+      }
+      n_negative <- c(n_negative, counter_negative)
+    }
+    dat$Feedback_N_Positive <- n_positive
+    dat$Feedback_N_Negative <- n_negative
+    newdf <- rbind(newdf, dat)
+  }
+}
+df <- newdf
+```
+
+### Recruitment History
+
+```{r}
+#| message: false
+
+# Consecutive count of participants per day (as area)
+dfsub |>
+  mutate(Date = as.Date(Date, format = "%d/%m/%Y")) |> 
+  group_by(Date) |> 
+  summarize(N = n()) |> 
+  ungroup() |>
+  mutate(N = cumsum(N)) |>
+  ggplot(aes(x = Date, y = N)) +
+  geom_area() +
+  scale_y_continuous(expand = c(0, 0)) +
+  labs(
+    title = "Recruitment History",
+    x = "Date",
+    y = "Total Number of Participants"
+  ) +
+  see::theme_modern()
+```
+
+
+## Manipulation Check
+
+### ISI
+
+```{r}
+#| code-fold: false
+
+df |> 
+  ggplot(aes(x=ISI, y=RT)) +
+  geom_point2(alpha=1/10, size=3) +
+  geom_smooth(method = "gam", formula = y ~ s(x, bs = "cs"), se = FALSE, linewidth=2) +
+  scale_y_log10() +
+  theme_minimal() +
+  coord_cartesian(ylim = c(0.1, 1))
+```
+
+### Threshold
+
+```{r}
+#| code-fold: false
+
+cumulative_median <- function(x) {
+  sapply(seq_along(x), function(i) median(x[1:i]))
+}
+
+df |> 
+  arrange(Participant, Session, Trial) |> 
+  group_by(Participant, Session) |> 
+  mutate(CumMedian = cumulative_median(RT)) |> 
+  ungroup() |> 
+  ggplot(aes(x=Threshold, y=CumMedian)) +
+  geom_abline(intercept = 0, slope = 1) +
+  geom_point(aes(color=Participant)) +
+  facet_grid(~Session) +
+  scale_x_log10() +
+  scale_y_log10() 
+```
+
+
+### Feedback
+
+
+
+
+
+```{r}
+#| code-fold: false
+
+df |>
+  filter(RT < 1) |> 
+  mutate(Feedback_N_Positive = ifelse(Feedback_N_Positive > 5, 6, Feedback_N_Positive),
+         Feedback_N_Negative = ifelse(Feedback_N_Negative > 3, 4, Feedback_N_Negative)) |>
+  ggplot(aes(x=Feedback_N_Positive, y=RT)) +
+  geom_smooth(method="lm", formula="y~poly(x, 2)") +
+  ggdist::stat_halfeye(aes(group=interaction(as.factor(Feedback_N_Positive), as.factor(Feedback_N_Negative)),
+                           fill=as.factor(Feedback_N_Negative)), alpha=1/3) +
+  scale_y_log10() +
+  coord_cartesian(ylim = c(0.15, 0.85))
+
+# t.test(RT ~ FeedbackN1, data = filter(!is.na(FeedbackN1)))
+```
+
+
+## Exclusion
+
+### Trials
+
+```{r}
+#| code-fold: false
+
+d <- bayestestR::estimate_density(df$RT, method="KernSmooth")
+dsub <- bayestestR::estimate_density(df, select="RT", at=c("Participant", "Session"), method="KernSmooth")
+
+d |>
+  ggplot(aes(x=x, y=y)) +
+  geom_vline(xintercept = c(0.15, 0.85), linetype = "dashed") +
+  geom_area(fill="grey") +
+  geom_line(data=dsub, aes(color=Participant, linetype=Session), linewidth=1) +
+  coord_cartesian(xlim = c(0, 1)) +
+  theme_minimal()
+
+outliers <- (df$RT > 0.85) | (df$RT < 0.15)
+df <- df[!outliers, ]
+```
+
+Outlier trials were removed from the dataset (`r sum(outliers)`).
+
+
+
+## Save 
+
+```{r}
+#| code-fold: false
+
+write.csv(df, "../data/data_game.csv", row.names = FALSE)
+write.csv(dfsub, "../data/data_participants.csv", row.names = FALSE)
+```

study1/analysis/2_optimization.qmd

@@ -0,0 +1,72 @@
+---
+title: "Doggo/Nogo (Study 1) - Task Optimization"
+editor: source
+editor_options: 
+  chunk_output_type: console
+format:
+  html:
+    code-fold: true
+    self-contained: true
+---
+
+## Data Preparation
+
+```{r}
+#| message: false
+#| warning: false
+
+library(tidyverse)
+library(easystats)
+library(patchwork)
+library(ggside)
+```
+
+
+```{r}
+#| code-fold: false
+
+dfsub <- read.csv("../data/data_participants.csv") 
+df <- read.csv("../data/data_game.csv") 
+```
+
+## ISI
+
+```{r}
+#| code-fold: false
+
+m <- glmmTMB::glmmTMB(RT ~ poly(ISI, 2) + (poly(ISI, 2)|Participant), 
+                      data = df)
+
+
+modelbased::estimate_relation(m, length=20) |> 
+  ggplot(aes(x = ISI, y = Predicted)) +
+  geom_ribbon(aes(ymin = CI_low, ymax = CI_high), alpha = 0.2) +
+  geom_line(data=modelbased::estimate_relation(m, length=20, include_random=TRUE), aes(color=Participant)) +
+  geom_line(linewidth=1) +
+  geom_hline(yintercept = 0.325, linetype = "dashed") +
+  theme_minimal() +
+  coord_cartesian(xlim = c(0, 4)) 
+```
+
+## Effect on Total Duration
+
+```{r}
+#| code-fold: false
+
+dat <- rbind(
+  data.frame(ISI = seq(1, 4, length.out = 500), Type = "1-4"),
+  data.frame(ISI = seq(1.5, 3.5, length.out = 500), Type = "1.5-3.5"),
+  data.frame(ISI = seq(2, 3.25, length.out = 500), Type = "2-3.25")
+)
+
+
+dat$RT <- insight::get_predicted(m, dat) 
+dat$Duration <- dat$RT + dat$ISI
+
+dat |> 
+  group_by(Type) |>
+  summarise(Duration = sum(Duration)) |> 
+  ggplot(aes(x=Type, y=Duration)) +
+  geom_line(aes(group=1)) +
+  theme_minimal()
+```

study1/analysis/best_model.jl

@@ -0,0 +1,86 @@
+using CSV
+using DataFrames
+using Turing
+using SequentialSamplingModels
+using StatsModels
+using StatsPlots
+using GLMakie
+using RCall
+
+
+# Predictions ===============================================================================
+
+cd(@__DIR__)  # pwd()
+include("fun_datagrid.jl")
+include("fun_data_poly.jl")
+
+df = CSV.read("../data/data_game.csv", DataFrame)
+
+
+# https://cosmicmar.com/MuseInference.jl/latest/
+# LBA model
+@model function model_lba(data; min_rt=minimum(data.rt), isi=nothing, participant=nothing)
+
+  # Transform ISI into polynomials
+  isi = data_poly(isi, 2; orthogonal=true)
+  ppt = unique(participant)
+
+  # Priors for coefficients
+  drift_intercept ~ filldist(truncated(Normal(3, 5), 0.0, Inf), 1)
+  drift_isi1 ~ filldist(Normal(0, 1), 1)
+  drift_isi2 ~ filldist(Normal(0, 1), 1)
+
+  # Prior for random intercepts (requires thoughtful specification)
+  # Participant-level intercepts' SD
+  drift_intercept_ppt_sd ~ truncated(Normal(0, 0.1), 0.0, Inf)
+  # Participant-level intercepts
+  drift_intercept_ppt ~ filldist(
+    Normal(0, drift_intercept_ppt_sd),
+    length(ppt)
+  )
+
+  σ ~ filldist(truncated(Normal(0, 1), 0.0, Inf), 1)
+  A ~ truncated(Normal(0.4, 0.4), 0.0, Inf)
+  k ~ truncated(Normal(0.2, 0.2), 0.0, Inf)
+  τ ~ truncated(Normal(0.2, 0.05), 0.0, min_rt)
+
+  for i in 1:length(data)
+    drift = drift_intercept .+ drift_intercept_ppt[findfirst(s -> s == participant[i], ppt)]
+    drift .+= drift_isi1 * isi[i, 1]
+    drift .+= drift_isi2 * isi[i, 2]
+    data[i] ~ LBA(drift, A, k, τ, σ)
+  end
+end
+
+# Fit 
+dat = [(choice=1, rt=df.RT[i]) for i in 1:nrow(df)]
+chain_lba = sample(model_lba(dat, min_rt=minimum(df.RT), isi=df.ISI, participant=df.Participant), NUTS(0.65, max_depth=8), 100)
+# StatsPlots.plot(chain_lba; size=(600, 2000))
+# summarystats(chain_wald)
+
+
+
+# Predictions
+grid = datagrid(df.ISI)
+pred = predict(model_lba([(missing) for i in 1:length(grid)]; min_rt=minimum(df.RT), isi=grid), chain_lba)
+pred = Array(pred)[:, 2:2:end]
+# Remove extreme
+pred[pred.>1] .= NaN
+pred = DataFrame(hcat(grid, transpose(pred)), vcat(:ISI, [Symbol("iter_$i") for i in 1:500]))
+
+
+@rput pred
+
+R"""
+library(tidyverse)
+library(bayestestR)
+library(ggdist)
+
+pred <- reshape_iterations(pred)
+# head(pred)
+pred |>
+  mutate(ISI = ISI) |> 
+  ggplot(aes(x = ISI, y = iter_value)) +
+  stat_halfeye() +
+  coord_cartesian(ylim = c(0.2, 0.6)) 
+"""