OSCR

Computational mechanisms of perception in autism revealed using games inspired by rodent operant tasks.

Code ↔ Paper

9 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 9 matches · 2 of them tie a paragraph to a whole file, not to given lines: weak matches, whose lines are not tinted
  1. [1] § MATERIALS AND METHODS › Data analysis › Nonlinear mixed-effects learning curve models ↔ scripts/run_nlme_learningcurve.R, lines 1–35 · score 0.70 · nonlinear mixed, fit separately, 11–200, 1–10, curve, model
  2. [2] § MATERIALS AND METHODS › Data analysis › Analysis of survey scores ↔ cave/pixi/pixi.min.js, the whole file · a weak match · score 0.62 · ABC, Composite, reverse, domains
  3. [3] § MATERIALS AND METHODS › Data analysis › Analysis of survey scores ↔ data_analysis_final.R, lines 233–283 · score 0.59 · sensory seeking, sensory avoiding, reverse, raw, AASP
  4. [4] § MATERIALS AND METHODS › Video game development ↔ cave/js/plugins.js, the whole file · a weak match · score 0.56 · 2.5 s, modified, touches, ms, video, clicks
  5. [5] § MATERIALS AND METHODS › Participants ↔ scripts/run_learningcurve_3group_ks.R, lines 1–47 · score 0.54 · missing VABS, 11–200, accuracy, SRS, score, game
  6. [6] § MATERIALS AND METHODS › Data analysis › SDT model ↔ data_analysis_final.R, lines 813–888 · score 0.53 · lose switch, win stay, bs, SRS, fit, Model
  7. [7] § MATERIALS AND METHODS › Data analysis › SDT model ↔ data_analysis_final.R, lines 813–888 · score 0.53 · lose switch, win stay, bs, SRS, fit, Model
  8. [8] § RESULTS › Integration noise correlates with symptom severity in social and adaptive domains ↔ data_analysis_final.R, lines 2221–2261 · score 0.52 · sensation seeking, perceptual noise, BIS, AASP, VABS, k1
  9. [9] § RESULTS › Integration noise correlates with symptom severity in social and adaptive domains ↔ data_analysis_final.R, lines 2222–2262 · score 0.52 · sensation seeking, perceptual noise, BIS, AASP, VABS, k1

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

R · 2,591 lines · 90 KB · CC0-1.0 · 3 matches

  1. # Created: 2020-07-07 15:00:00
  2. # Last modified: 2020-07-07 15:00:00
  3. setwd("/Users/suchetachakravarty/Documents/projects/geode/")
  4. rm(list=ls())
  5. # Load the required libraries
  6. library(data.table)
  7. library(ggplot2)
  8. library(dplyr)
  9. library(ggpubr)
  10. library(rstatix)
  11. library(gridExtra)
  12. library(Hmisc)
  13. library(rstan)
  14. library(zoo)
  15. library(latex2exp)
  16. library(afex)
  17. library(emmeans)
  18. library(minpack.lm)
  19. library(tidyverse)
  20. library(glmnet)
  21. library(lme4)
  22. library(broom.mixed)
  23. # load the data
  24. load("data/preprocessed_data.RData")
  25. # create a binary group variable
  26. alldata[, binary_group := fifelse(grepl("ASD", group, ignore.case = TRUE), "ASD", "TD")]
  27. # create a verbal group variable
  28. alldata[, group_verbal := fcase(
  29. is.na(vinabcstd), "flag",
  30. binary_group == "ASD" & vinabcstd >= 75, "ASD",
  31. binary_group == "TD" & vinabcstd >= 75, "TD",
  32. binary_group == "ASD" & vinabcstd < 75, "ASD (low verbal)",
  33. binary_group == "TD" & vinabcstd < 75, "flag_low_td",
  34. default = "Uncategorized"
  35. )]
  36. # create a srs normal group variable
  37. alldata[, srs_normal := fifelse(
  38. (binary_group == "ASD" & srs2total > 20) | (binary_group == "TD" & srs2total < 120),1,0)]
  39. # alldata[, table(srs_normal)]
  40. # create a above criteria group variable
  41. alldata[, above_crit := fifelse(
  42. (nTrials == 200 & perf > .616),1,0)]
  43. # data summary
  44. summary <- unique(alldata[, .(username, subID, game_version, source, binary_group, group_verbal, srs_normal, gender, race, age, ethnicity, vinabcstd, srs2total, bistotal,aasp_low_reg_raw,
  45. aasp_sen_seek_raw,aasp_sen_ses_raw,aasp_sen_avoid_raw, perf, trainperf, nTrials, above_crit)])
  46. ###### table for subject numbers ####
  47. summary[srs_normal == 1 & game_version=="ft", table(binary_group)]
  48. summary[srs_normal == 1 & game_version=="ft" & group_verbal!="flag", table(binary_group)]
  49. summary[above_crit==1 & srs_normal == 1 & game_version=="ft" & group_verbal!="flag" & !is.na(bistotal) & !is.na(aasp_low_reg_raw) & !is.na(aasp_sen_seek_raw) & !is.na(aasp_sen_ses_raw) & !is.na(aasp_sen_avoid_raw), table(binary_group)]
  50. ############################################################################################################
  51. # Data for Supplemenatry Figure 1
  52. # total number of participants in each binary group for game_version = "ft"
  53. summary[game_version == "ft", table(binary_group)] #+1 ASD subject whose game data we could not track
  54. # total number of participants in each binary group for game_version = "ft" and srs_normal = 1
  55. summary[game_version == "ft" & srs_normal == 1, table(binary_group)]
  56. # total number of participants in each binary group for game_version = "ft" and srs_normal = 1 and above_crit = 1
  57. summary[game_version == "ft" & srs_normal == 1 & above_crit == 1, table(binary_group)]
  58. # number of participants in binary_group = "ASD" for game_version = "ft" and srs_normal = 1 and above_crit = 0
  59. summary[game_version == "ft" & srs_normal == 1 & above_crit == 0, table(binary_group)]
  60. # subIDs of the ASD participants who did not meet the criteria
  61. id1 = summary[game_version == "ft" & above_crit == 0 & binary_group == "ASD", .(subID)]
  62. # subIDs ofparticipants who played game again
  63. id2 = summary[game_version != "ft", .(subID)]
  64. # remove the "_2" and "_3" at the end of the subID for id2 only where it is present
  65. id2[, subID := gsub("_2$|_3$", "", subID)]
  66. # common subIDs between id1 and id2
  67. common_ids = id1[subID %in% id2$subID, .(subID)]
  68. # subIDs of the ASD participants who played game_version = "ft-2" and met the criteria
  69. id3 = summary1[game_version == "ft-2" & srs_normal == 1 & above_crit == 1 & binary_group == "ASD", .(subID)]
  70. # clear id1, id2, common_ids
  71. rm(id1,id2,common_ids)
  72. ############################################################################################################
  73. # Supplemenatry Figure 2
  74. # only include subjects with game_version = "ft" and srs_normal = 1
  75. df = summary[game_version == "ft" & srs_normal==1]
  76. # gender distribution
  77. # If gender = "M", Male, if "F", Female, otherwise call it Other/Unknown
  78. df[,gender := fcase(gender == "M", "Male",gender == "F","Female", default = "Other/Unknown")]
  79. p1<-
  80. ggplot(df,aes(x=gender, fill=binary_group, color=binary_group)) +
  81. geom_bar(position="dodge") +
  82. theme_bw()+labs(fill="",color="",xlab="",
  83. # ylab="Frequency"
  84. )+
  85. theme(legend.position = c(.7,.8),
  86. panel.grid.major = element_blank(),
  87. panel.grid.minor = element_blank(),
  88. axis.text.x = element_text(angle = 90, vjust = 1, hjust=1),
  89. legend.background=element_blank())
  90. # age distribution
  91. p2 <-
  92. ggplot(df,aes(x=age,fill=binary_group,color=binary_group))+
  93. geom_bar(width=.5,position=position_dodge(width = .5),alpha=.7)+
  94. theme_bw()+
  95. scale_x_continuous(breaks = c(11, 12, 13,14,15,16,17),limits = c(10,18))+
  96. # ylab("Frequency")+
  97. xlab("Age (y)")+
  98. labs(fill="",color="")+
  99. theme(legend.position = "None",
  100. panel.grid.major = element_blank(),
  101. panel.grid.minor = element_blank(),
  102. legend.background=element_blank())
  103. # race distribution
  104. unique(df$race)
  105. df[,race := fcase(race=="LA6156-9","Asian",
  106. race=="LA10610-6","Black or African American",
  107. race=="LA10611-4","Native Hawaiian or Other Pacific Islander",
  108. race=="LA4457-3","White",
  109. race=="LA4489-6","Native Hawaiian or Other Pacific Islander",default="Unknown")]
  110. df$race<-factor(df$race,levels = c("White","Black or African American", "Asian",
  111. "Native Hawaiian or Other Pacific Islander","Unknown"))
  112. p3 <-
  113. ggplot(df,aes(x=race,fill=binary_group,color=binary_group))+
  114. geom_bar(width=.5,position=position_dodge(width = .5),alpha=.7)+
  115. theme_bw()+
  116. # ylab("Frequency")+
  117. xlab("")+
  118. labs(fill="",color="")+
  119. scale_x_discrete(labels = c("White","Black/African\nAmerican","Asian",
  120. "Native Hawaiian/\nPacific Islander","Unknown"))+
  121. theme(legend.position = "None",
  122. panel.grid.major = element_blank(),
  123. axis.text.x = element_text(angle = 90, vjust = 1, hjust=1),
  124. panel.grid.minor = element_blank())
  125. # ethnicity distribution
  126. unique(df$ethnicity)
  127. df[,ethnicity := fcase(ethnicity=="2135-2","Hispanic or Latino",
  128. ethnicity=="2186-5","Not Hispanic or Latino",
  129. default="Unknown")]
  130. df$ethnicity<-factor(df$ethnicity,levels = c("Hispanic or Latino","Not Hispanic or Latino","Unknown"))
  131. p4 <-
  132. ggplot(df,aes(x=ethnicity,fill=binary_group,color=binary_group))+
  133. geom_bar(width=.5,position=position_dodge(width = .5),alpha=.7)+
  134. theme_bw()+
  135. # ylab("Frequency")+
  136. xlab("")+
  137. labs(fill="",color="")+
  138. scale_x_discrete(labels = c("Hispanic/Latino","Not Hispanic/Latino","Unknown"))+
  139. theme(legend.position = "None",
  140. panel.grid.major = element_blank(),
  141. axis.text.x = element_text(angle = 90, vjust = 1, hjust=1),
  142. panel.grid.minor = element_blank())
  143. # srs2total distribution
  144. p5 <-
  145. ggplot(df,
  146. aes(x=srs2total))+
  147. geom_histogram(aes(#y = after_stat(density),
  148. fill=binary_group,
  149. color=binary_group),
  150. position="identity",
  151. alpha=.5)+
  152. theme_bw()+
  153. xlab("SRS-2")+
  154. geom_vline(xintercept=70,linetype="dashed")+
  155. # ylab("Frequency")+
  156. labs(fill="",color="")+
  157. theme(legend.position = "None",
  158. panel.grid.major = element_blank(),
  159. panel.grid.minor = element_blank())
  160. # bistotal distribution
  161. p6 <-
  162. ggplot(df,
  163. aes(x=bistotal))+
  164. geom_histogram(aes(#y = after_stat(density),
  165. fill=binary_group,
  166. color=binary_group),
  167. position="identity",
  168. alpha=.5)+
  169. theme_bw()+
  170. xlab("BIS")+
  171. # ylab("Frequency")+
  172. labs(fill="",color="")+
  173. theme(legend.position = "None",
  174. panel.grid.major = element_blank(),
  175. panel.grid.minor = element_blank())
  176. # reverse scale vinabcstd distribution
  177. df[,vinabcstd := 140 - vinabcstd]
  178. # vinabcstd distribution
  179. p7 <-
  180. ggplot(df,
  181. aes(x=vinabcstd))+
  182. geom_histogram(aes(#y = after_stat(density),
  183. fill=binary_group,
  184. color=binary_group),
  185. position="identity",
  186. alpha=.5)+
  187. theme_bw()+
  188. xlab("VABS-3")+
  189. # ylab("Frequency")+
  190. labs(fill="",color="")+
  191. theme(legend.position = "None",
  192. panel.grid.major = element_blank(),
  193. panel.grid.minor = element_blank())
  194. # aasp_low_reg_raw distribution
  195. p8 <-
  196. ggplot(df,
  197. aes(x=aasp_low_reg_raw))+
  198. geom_histogram(aes(#y = after_stat(density),
  199. fill=binary_group,
  200. color=binary_group),
  201. position="identity",
  202. alpha=.5)+
  203. theme_bw()+
  204. xlab("AASP\n(Low Registration)")+
  205. # ylab("Frequency")+
  206. labs(fill="",color="")+
  207. theme(legend.position = "None",
  208. panel.grid.major = element_blank(),
  209. panel.grid.minor = element_blank())
  210. # aasp_sen_avoid_raw distribution
  211. p9 <-
  212. ggplot(df,
  213. aes(x=aasp_sen_avoid_raw))+
  214. geom_histogram(aes(#y = after_stat(density),
  215. fill=binary_group,
  216. color=binary_group),
  217. position="identity",
  218. alpha=.5)+
  219. theme_bw()+
  220. xlab("AASP\n(Sensory Avoiding)")+
  221. # ylab("Frequency")+
  222. labs(fill="",color="")+
  223. theme(legend.position = "None",
  224. panel.grid.major = element_blank(),
  225. panel.grid.minor = element_blank())
  226. # reverse scale aasp_sen_seek_raw distribution
  227. df[,aasp_sen_seek_raw := 75 - aasp_sen_seek_raw]
  228. p10 <-
  229. ggplot(df,
  230. aes(x=aasp_sen_seek_raw))+
  231. geom_histogram(aes(#y = after_stat(density),
  232. fill=binary_group,
  233. color=binary_group),
  234. position="identity",
  235. alpha=.5)+
  236. theme_bw()+
  237. xlab("AASP\n(Sensory Seeking)")+
  238. # ylab("Frequency")+
  239. labs(fill="",color="")+
  240. theme(legend.position = "None",
  241. panel.grid.major = element_blank(),
  242. panel.grid.minor = element_blank())
  243. # aasp_sen_ses_raw distribution
  244. p11 <-
  245. ggplot(df,
  246. aes(x=aasp_sen_ses_raw))+
  247. geom_histogram(aes(#y = after_stat(density),
  248. fill=binary_group,
  249. color=binary_group),
  250. position="identity",
  251. alpha=.5)+
  252. theme_bw()+
  253. xlab("AASP\n(Sensory Sensitivity)")+
  254. # ylab("Frequency")+
  255. labs(fill="",color="")+
  256. theme(legend.position = "None",
  257. panel.grid.major = element_blank(),
  258. panel.grid.minor = element_blank())
  259. # combine all the plots into publication ready format, add labels and save the plot
  260. ggarrange(p1,p2,p3,p4,p5,p6,p7,p8,p9,p10,p11,ncol=4,nrow=3,common.legend = TRUE,align = "v",
  261. labels = c("A","B","C","D","E","F","G","H","I","J","K")) %>% ggexport(filename = "results/figures/geodems/demographics_v2.pdf")
  262. # remove plot objects
  263. rm(p1,p2,p3,p4,p5,p6,p7,p8,p9,p10,p11)
  264. ############################################################################################################
  265. # Main Figure 2
  266. # distribution of vinabcstc for group_verbal (remove flagged subjects) and for srs_normal=1 and game_version = "ft"
  267. df = summary[game_version == "ft" & srs_normal==1 & group_verbal != "flag" & group_verbal != "flag_low_td"]
  268. df$group_verbal <- factor(df$group_verbal, levels = c("TD","ASD","ASD (low verbal)"),ordered = TRUE)
  269. # plot distribution of vinabcstd for group_verbal
  270. p1 <-
  271. ggplot(df,
  272. aes(x=vinabcstd))+
  273. geom_histogram(aes(
  274. fill=group_verbal,color=group_verbal),
  275. position="identity",
  276. alpha=.5)+
  277. theme_bw()+
  278. scale_fill_manual(values = c("blue","orange","red"))+
  279. scale_color_manual(values = c("blue","orange","red"))+
  280. xlab("VABS-3")+
  281. ylab("Frequency")+
  282. labs(fill="",color="")+
  283. theme(legend.position = "None",
  284. panel.grid.major = element_blank(),
  285. panel.grid.minor = element_blank(),
  286. legend.key = element_rect(colour = "transparent"),
  287. legend.text=element_text(size=4))
  288. # scatter plot perf by vinabcstd
  289. p2 <-
  290. ggplot(df,aes(x=vinabcstd,y=perf*100))+
  291. geom_point(aes(color=group_verbal))+
  292. ylab("%Correct")+
  293. xlab("VABS-3")+
  294. scale_color_manual(values = c("blue","orange","red"))+
  295. geom_hline(yintercept = 62,linetype="dashed",size=.5)+
  296. theme_bw()+
  297. labs(color="")+
  298. theme(legend.position = "None",
  299. panel.grid.major = element_blank(),
  300. panel.grid.minor = element_blank())
  301. # plot percentage of subjects above the criteria for each group_verbal, show the percentage on the plot
  302. df = df[, .(above_crit = mean(above_crit)), by = .(group_verbal)]
  303. df$group_verbal <- factor(df$group_verbal, levels = c("TD","ASD","ASD (low verbal)"),ordered = TRUE)
  304. p3 <-
  305. ggplot(df,
  306. aes(x=group_verbal,y=above_crit))+
  307. geom_bar(stat="identity",aes(fill=group_verbal),position="dodge")+
  308. geom_text(aes(label = paste0(round(above_crit*100,1),"%")),position = position_stack(vjust = 0.5), col = "white")+
  309. theme_bw()+
  310. scale_fill_manual(values = c("blue","orange","red"))+
  311. xlab("")+
  312. ylab("%above criteria")+
  313. labs(fill="")+
  314. theme(legend.position = "None",
  315. axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1),
  316. panel.grid.major = element_blank(),
  317. panel.grid.minor = element_blank(),
  318. legend.key = element_rect(colour = "transparent"),
  319. legend.text=element_text(size=4))
  320. # score for trials 1 through 10 for each group_verbal from alldata, for game_version = "ft" and srs_normal = 1 and above_crit = 1
  321. df = alldata[game_version == "ft" & srs_normal == 1 & above_crit == 1 & group_verbal != "flag" & group_verbal != "flag_low_td"]
  322. df = df[,c("subID","group_verbal","trial","score")]
  323. df = df[trial %in% 1:10]
  324. # plot the average scores and se for each group_verbal
  325. df = df[, .(mean_score = mean(score), se = sd(score)/sqrt(.N)), by = .(group_verbal,trial)]
  326. df$group_verbal <- factor(df$group_verbal, levels = c("TD","ASD","ASD (low verbal)"),ordered = TRUE)
  327. pairwise_test <- df %>%
  328. pairwise_t_test(mean_score ~ group_verbal, p.adjust.method = "bonferroni")
  329. pairwise_test <- pairwise_test %>%
  330. add_xy_position(x = "group_verbal")
  331. fit_learning_curve <- function(data) {
  332. nls(mean_score ~ pinf - (pinf - p0) * exp(-alpha * trial),
  333. data = data,
  334. start = list(pinf = 1, p0 = 0.5, alpha = 0.1),
  335. control = nls.control(maxiter = 100))
  336. }
  337. td_fit <- df %>%
  338. filter(group_verbal == "TD") %>%
  339. fit_learning_curve()
  340. asd_fit <- df %>%
  341. filter(group_verbal == "ASD") %>%
  342. fit_learning_curve()
  343. asd_low_verbal_fit <- df %>%
  344. filter(group_verbal == "ASD (low verbal)") %>%
  345. fit_learning_curve()
  346. new_data <- data.frame(trial = seq(1, 10, 0.1))
  347. td_pred <- data.frame(
  348. trial = new_data$trial,
  349. mean_score = predict(td_fit, newdata = new_data),
  350. group_verbal = "TD"
  351. )
  352. asd_pred <- data.frame(
  353. trial = new_data$trial,
  354. mean_score = predict(asd_fit, newdata = new_data),
  355. group_verbal = "ASD"
  356. )
  357. asd_low_verbal_pred <- data.frame(
  358. trial = new_data$trial,
  359. mean_score = predict(asd_low_verbal_fit, newdata = new_data),
  360. group_verbal = "ASD (low verbal)"
  361. )
  362. pred_data <- rbind(td_pred, asd_pred)
  363. pred_data <- rbind(pred_data, asd_low_verbal_pred)
  364. p4 <-
  365. ggplot() +
  366. geom_point(data = df, aes(x = trial, y = 100*mean_score, color = group_verbal)) +
  367. geom_errorbar(data = df, aes(x = trial, ymin = 100*(mean_score - se), ymax = 100*(mean_score + se), color = group_verbal), width = 0.2) +
  368. geom_line(data = pred_data, aes(x = trial, y = 100*mean_score, color = group_verbal)) +
  369. # stat_pvalue_manual(pairwise_test, label = "p.adj.signif",
  370. # y.position = max(df$mean_score*100) + 0.1) +
  371. scale_x_continuous(breaks = 1:10) +
  372. labs(title = "Level 1",
  373. x = "Trial",
  374. y = "%Correct",
  375. color = "Group") +
  376. theme_bw() +
  377. ylim(60,105)+
  378. scale_color_manual(values = c("TD" = "blue", "ASD" = "orange", "ASD (low verbal)" = "red"))+
  379. theme(legend.position = "None",
  380. panel.grid.major = element_blank(),
  381. panel.grid.minor = element_blank(),
  382. legend.key = element_rect(colour = "transparent"),
  383. legend.text=element_text(size=4))
  384. # rolling average (trials 11 through 200) for each group_verbal
  385. df = alldata[game_version == "ft" & srs_normal == 1 & above_crit == 1 & group_verbal != "flag" & group_verbal != "flag_low_td"]
  386. df = df[,c("subID","group_verbal","trial","score")]
  387. df = df[trial %in% 11:200]
  388. window_size <- 50
  389. df_rolling <- df %>%
  390. group_by(subID, group_verbal) %>%
  391. arrange(trial) %>%
  392. mutate(rolling_avg = rollmean(score, k = window_size, fill = NA, align = "left")) %>% ungroup()
  393. df_summary <- df_rolling %>%
  394. group_by(group_verbal, trial) %>%
  395. summarise(
  396. avg_score = mean(rolling_avg, na.rm = TRUE),
  397. se = sd(rolling_avg, na.rm = TRUE) / sqrt(n()),
  398. .groups = "drop"
  399. )
  400. # fit_learning_curve <- function(data) {
  401. # nls(avg_score ~ pinf - (pinf - p0) * exp(-alpha * trial),
  402. # data = data,
  403. # start = list(pinf = 1, p0 = 0.5, alpha = 0.1),
  404. # control = nls.control(maxiter = 100))
  405. # }
  406. # td_fit <- df_summary %>%
  407. # filter(group_verbal == "TD") %>%
  408. # fit_learning_curve()
  409. # asd_fit <- df_summary %>%
  410. # filter(group_verbal == "ASD") %>%
  411. # fit_learning_curve()
  412. # asd_low_verbal_fit <- df_summary %>%
  413. # filter(group_verbal == "ASD (low verbal)") %>%
  414. # fit_learning_curve()
  415. p5 <-
  416. ggplot(df_summary, aes(x = trial, y = 100*avg_score, color = group_verbal)) +
  417. # geom_point() +
  418. geom_smooth(method = "loess",span=1,se=F,size=.5) +
  419. geom_ribbon(aes(ymin = 100*(avg_score - se), ymax = 100*(avg_score + se), fill = group_verbal), alpha = .4) +
  420. labs(title = "Rolling Average Score by Group (Trials 11-200)",
  421. x = "Trial",
  422. y = "Rolling Average Score",
  423. color = "Group",
  424. fill = "Group") +
  425. theme_bw() +
  426. ylim(75,90)+
  427. labs(title = "Level 2-10",
  428. x = "Trial",
  429. y = "%Correct",
  430. color = "Group") +
  431. scale_x_continuous(breaks = seq(10, 150, 20),limits = c(10,150)) +
  432. scale_color_manual(values = c("ASD" = "orange", "ASD (low verbal)" = "red", "TD" = "blue")) +
  433. scale_fill_manual(values = c("ASD" = "orange", "ASD (low verbal)" = "red", "TD" = "blue"))+
  434. theme(legend.position = "None",
  435. panel.grid.major = element_blank(),
  436. panel.grid.minor = element_blank(),
  437. legend.key = element_rect(colour = "transparent"),
  438. legend.text=element_text(size=4))
  439. # psychometric curve (unsigned)
  440. df = alldata[game_version == "ft" & srs_normal == 1 & above_crit == 1 & group_verbal != "flag" & group_verbal != "flag_low_td"]
  441. df = df[,c("subID","group_verbal","trial","score","dflash")]
  442. df = df[trial %in% 101:200 & dflash != 0]
  443. df$dflash <- abs(df$dflash)
  444. df_subject_avg <- df %>%
  445. group_by(subID, group_verbal, dflash) %>%
  446. summarise(avg_score = mean(score, na.rm = TRUE), .groups = "drop")
  447. df_long <- df_subject_avg %>%
  448. select(subID, group_verbal, dflash, avg_score) %>%
  449. mutate(dflash = factor(dflash,levels = c(1:10))) # Convert dflash to a factor
  450. aov_model <- aov_ez(
  451. id = "subID",
  452. dv = "avg_score",
  453. data = df_long,
  454. between = "group_verbal",
  455. within = "dflash",
  456. type = 3
  457. )
  458. print(aov_model)
  459. posthoc <- emmeans(aov_model, ~ group_verbal | dflash)
  460. group_comparisons <- pairs(posthoc)
  461. print(group_comparisons)
  462. significant_dflash <- group_comparisons %>%
  463. as.data.frame() %>%
  464. filter(p.value < 0.05) %>%
  465. pull(dflash) %>%
  466. unique()
  467. df_group_avg <- df_subject_avg %>%
  468. group_by(group_verbal, dflash) %>%
  469. summarise(
  470. mean_score = mean(avg_score, na.rm = TRUE),
  471. se = sd(avg_score, na.rm = TRUE) / sqrt(n()),
  472. .groups = "drop"
  473. )
  474. p6<-
  475. ggplot(df_group_avg, aes(x = dflash, y = mean_score, color = group_verbal)) +
  476. geom_point(size=.5)+
  477. geom_smooth(se=F,method="glm",method.args = list(family = "binomial"),size=.5) +
  478. # geom_ribbon(aes(ymin = mean_score - se, ymax = mean_score + se, fill = group_verbal), alpha = 0.2) +
  479. geom_errorbar(aes(ymin = mean_score - se, ymax = mean_score + se), alpha = 0.8, size=.5,width=.5) +
  480. labs(title = "Trial 101-200",
  481. x = "Cue difference",
  482. y = "Accuracy",
  483. color = "Group",
  484. fill = "Group") +
  485. ylim(.65,1.05)+
  486. theme_bw() +
  487. scale_x_continuous(breaks = seq(1, 10),limits = c(.5,10.5)) +
  488. scale_color_manual(values = c("ASD" = "orange", "ASD (low verbal)" = "red", "TD" = "blue")) +
  489. scale_fill_manual(values = c("ASD" = "orange", "ASD (low verbal)" = "red", "TD" = "blue"))+
  490. theme(legend.position = "None",
  491. panel.grid.major = element_blank(),
  492. panel.grid.minor = element_blank())
  493. # p <- p + geom_point(data = df_group_avg %>%
  494. # filter(dflash %in% significant_dflash) %>%
  495. # group_by(dflash) %>%
  496. # summarise(max_y = max(mean_score + se)),
  497. # aes(x = dflash, y = max_y + 0.05),
  498. # shape = 8, size = 3, color = "black")
  499. ggarrange(p1,p2,p3,p4,p5,p6,ncol=3,nrow=3,common.legend = TRUE,align = "v",
  500. labels = c("A","B","C","D","E","F")) %>% ggexport(filename = "results/figures/geodems/descriptives_v2.pdf")
  501. # remove variables
  502. rm(df,df_group_avg,df_long,df_rolling,df_subject_avg,group_comparisons,pairwise_test,posthoc,significant_dflash,td_fit,asd_fit,asd_low_verbal_fit)
  503. rm(td_pred,asd_pred,asd_low_verbal_pred,df_summary,new_data,pred_data)
  504. rm(aov_model,window_size,fit_learning_curve)
  505. # remove plot objects
  506. rm(p1,p2,p3,p4,p5,p6)
  507. ####### Supplementary Figure 3 ######################
  508. # the full psychometric curve
  509. df = alldata[game_version == "ft" & srs_normal == 1 & above_crit == 1 & group_verbal != "flag" & group_verbal != "flag_low_td"]
  510. df = df[,c("subID","group_verbal","trial","choice","dflash","score")]
  511. df = df[trial %in% 11:200]
  512. df_subject_avg <- df %>%
  513. group_by(subID, group_verbal, dflash) %>%
  514. summarise(avg_choice = mean(choice, na.rm = TRUE), .groups = "drop")
  515. df_group_avg <- df_subject_avg %>% group_by(group_verbal, dflash) %>% summarise(
  516. mean_choice = mean(avg_choice, na.rm = TRUE),
  517. se = sd(avg_choice, na.rm = TRUE) / sqrt(n()),
  518. .groups = "drop")
  519. df_group_avg$group_verbal <- factor(df_group_avg$group_verbal, levels = c("TD","ASD","ASD (low verbal)"),ordered = TRUE)
  520. p1<-
  521. ggplot(df_group_avg, aes(x = dflash, y = mean_choice, color = group_verbal)) +
  522. geom_point(size=.5)+
  523. geom_smooth(se=F,method="glm",method.args = list(family = "binomial"),size=.5) +
  524. # geom_ribbon(aes(ymin = mean_choice - se, ymax = mean_choice + se, fill = group_verbal), alpha = 0.2) +
  525. geom_errorbar(aes(ymin = mean_choice - se, ymax = mean_choice + se), alpha = 0.8, size=.5,width=.5) +
  526. labs(title = "Trial 11-200",
  527. x = "Cue difference",
  528. y = "%Went right",
  529. color = "",
  530. fill = ""
  531. ) +
  532. ylim(-.05,1.05)+
  533. theme_bw() +
  534. scale_x_continuous(breaks = seq(-10, 10,3),limits = c(-11,11)) +
  535. scale_color_manual(values = c("ASD" = "orange", "ASD (low verbal)" = "red", "TD" = "blue")) +
  536. scale_fill_manual(values = c("ASD" = "orange", "ASD (low verbal)" = "red", "TD" = "blue"))+
  537. theme(legend.position = "None",
  538. panel.grid.major = element_blank(),
  539. panel.grid.minor = element_blank())
  540. df$dflash <- abs(df$dflash)
  541. df = df[df$dflash != 0]
  542. df_subject_avg <- df %>%
  543. group_by(subID, group_verbal, dflash) %>%
  544. summarise(avg_score = mean(score, na.rm = TRUE), .groups = "drop")
  545. df_group_avg <- df_subject_avg %>% group_by(group_verbal, dflash) %>% summarise(
  546. mean_score = mean(avg_score, na.rm = TRUE),
  547. se = sd(avg_score, na.rm = TRUE) / sqrt(n()),
  548. .groups = "drop")
  549. df_group_avg$group_verbal <- factor(df_group_avg$group_verbal, levels = c("TD","ASD","ASD (low verbal)"),ordered = TRUE)
  550. p2<-
  551. ggplot(df_group_avg, aes(x = dflash, y = mean_score, color = group_verbal)) +
  552. geom_point(size=.5)+
  553. geom_smooth(se=F,method="glm",method.args = list(family = "binomial"),size=.5) +
  554. # geom_ribbon(aes(ymin = mean_choice - se, ymax = mean_choice + se, fill = group_verbal), alpha = 0.2) +
  555. geom_errorbar(aes(ymin = mean_score - se, ymax = mean_score + se), alpha = 0.8, size=.5,width=.5) +
  556. labs(title = "Trial 11-200",
  557. x = "Cue difference (absolute)",
  558. y = "Accuracy",
  559. color = "",
  560. fill = ""
  561. ) +
  562. ylim(.65,1.05)+
  563. theme_bw() +
  564. scale_x_continuous(breaks = seq(1, 10),limits = c(.5,10.5)) +
  565. scale_color_manual(values = c("ASD" = "orange", "ASD (low verbal)" = "red", "TD" = "blue")) +
  566. scale_fill_manual(values = c("ASD" = "orange", "ASD (low verbal)" = "red", "TD" = "blue"))+
  567. theme(legend.position = "None",
  568. panel.grid.major = element_blank(),
  569. panel.grid.minor = element_blank())
  570. df <- alldata[trial>10 & srs_normal == 1 & above_crit == 1 & game_version == "ft",c("username","subID","choice","dflash","winstay","loseswitch","binary_group","group_verbal")]
  571. df <- df[group_verbal != "flag" & group_verbal != "flag_low_td"]
  572. # create trial bins
  573. df <- df %>%
  574. group_by(subID) %>%
  575. mutate(trial= row_number(),
  576. trial_bin = cut(trial, breaks = 19, labels = FALSE)) %>%
  577. ungroup()
  578. # function to fit a glm for each trial bin
  579. fit_bin_glm <- function(bin_data) {
  580. glm(choice ~ dflash + winstay + loseswitch,
  581. data = bin_data,
  582. family = binomial(link = "logit"))
  583. }
  584. # fit the models
  585. model_fits <- df %>%
  586. group_by(group_verbal, trial_bin) %>%
  587. nest() %>%
  588. mutate(model = map(data, fit_bin_glm))
  589. # get parameter estimates for each bin and group
  590. param_estimates <- model_fits %>%
  591. mutate(params = map(model, broom::tidy)) %>%
  592. unnest(params) %>%
  593. select(group_verbal, trial_bin, term, estimate, std.error)
  594. param_estimates <- param_estimates %>%
  595. mutate(significant = abs(estimate / std.error) > 1.96) # p < 0.05
  596. # Define the learning curve function
  597. learning_curve <- function(n, p0, pinf, alpha) {
  598. pinf - (pinf - p0) * exp(-alpha * n)
  599. }
  600. # Filter data for dflash and split by group
  601. dflash_data <- param_estimates %>%
  602. filter(term == "dflash") %>%
  603. split(.$group_verbal)
  604. # Function to fit the learning curve for a group
  605. fit_learning_curve <- function(data) {
  606. nlsLM(estimate ~ learning_curve(trial_bin, p0, pinf, alpha),
  607. data = data,
  608. start = list(p0 = min(data$estimate),
  609. pinf = max(data$estimate),
  610. alpha = 0.1),
  611. lower = c(0.5, 0.5, 0.1),
  612. upper = c(1, 1, 1))
  613. }
  614. # Fit the model for each group
  615. fits <- map(dflash_data, safely(fit_learning_curve))
  616. # Extract parameters and create predictions
  617. results <- imap_dfr(fits, function(fit, group_verbal) {
  618. if (is.null(fit$error)) {
  619. params <- coef(fit$result)
  620. new_data <- data.frame(trial_bin = seq(1, 19, length.out = 100))
  621. new_data$predicted <- predict(fit$result, newdata = new_data)
  622. new_data$group_verbal <- group_verbal
  623. return(new_data)
  624. } else {
  625. return(NULL)
  626. }
  627. })
  628. # Print the fitted parameters
  629. walk2(fits, names(fits), function(fit, group_verbal) {
  630. if (is.null(fit$error)) {
  631. params <- coef(fit$result)
  632. cat(sprintf("%s: p0 = %.4f, pinf = %.4f, alpha = %.4f\n",
  633. group_verbal, params["p0"], params["pinf"], params["alpha"]))
  634. } else {
  635. cat(sprintf("%s: Fitting failed\n", binary_group))
  636. }
  637. })
  638. # Plot
  639. p3 <-
  640. ggplot() +
  641. geom_point(data = param_estimates %>% filter(term == "dflash"),
  642. aes(x = trial_bin*10, y = estimate, color = group_verbal),size=.5) +
  643. geom_ribbon(data = param_estimates %>% filter(term == "dflash"),
  644. aes(x = trial_bin*10, ymin = estimate - std.error, ymax = estimate + std.error, fill = group_verbal),
  645. alpha = 0.2) +
  646. geom_errorbar(data = param_estimates %>% filter(term == "dflash"),
  647. aes(x = trial_bin*10, ymin = estimate - std.error, ymax = estimate + std.error, color = group_verbal),
  648. width = .5,size=.5) +
  649. geom_line(data = results,
  650. aes(x = trial_bin*10, y = predicted, color = group_verbal)) +
  651. labs(title = "",
  652. x = "Trial",
  653. y = "Psychometric Slope",
  654. color="") +
  655. theme_bw()+
  656. scale_x_continuous(breaks = seq(11, 200,30),limits = c(5,205)) +
  657. scale_color_manual(values = c("ASD" = "orange", "ASD (low verbal)" = "red", "TD" = "blue")) +
  658. scale_fill_manual(values = c("ASD" = "orange", "ASD (low verbal)" = "red", "TD" = "blue"))+
  659. theme(legend.position = "None",
  660. panel.grid.major = element_blank(),
  661. panel.grid.minor = element_blank())
  662. p4 <-
  663. ggplot() +
  664. geom_point(data = param_estimates %>% filter(term == "(Intercept)"), aes(x = trial_bin*10, y = estimate, color = group_verbal),size=.5) +
  665. # geom_line(data = new_data, aes(x = trial_bin, y = predicted), color = "red") +
  666. geom_ribbon(data = param_estimates %>% filter(term == "(Intercept)"),
  667. aes(x = trial_bin*10, ymin = estimate - std.error, ymax = estimate + std.error, fill = group_verbal),
  668. alpha = 0.2) +
  669. geom_errorbar(data = param_estimates %>% filter(term == "(Intercept)"),
  670. aes(x = trial_bin*10, ymin = estimate - std.error, ymax = estimate + std.error, color = group_verbal),
  671. width = .5,size=.5) +
  672. geom_smooth(data = param_estimates %>% filter(term == "(Intercept)"),
  673. aes(x = trial_bin*10, y = estimate, color = group_verbal),
  674. method = "lm",
  675. # formula = y ~ s(x, bs = "cs"),
  676. se = F,size=.5) +
  677. # facet_wrap(~ term, scales = "free_y", ncol = 3, nrow = 3) +
  678. labs(title = "",
  679. x = "Trial",
  680. y = "Side-bias",
  681. color="",
  682. fill="") +
  683. ylim(-.8,.8)+
  684. theme_bw()+
  685. scale_x_continuous(breaks = seq(11, 200,30),limits = c(5,205)) +
  686. scale_color_manual(values = c("ASD" = "orange", "ASD (low verbal)" = "red", "TD" = "blue")) +
  687. scale_fill_manual(values = c("ASD" = "orange", "ASD (low verbal)" = "red", "TD" = "blue"))+
  688. theme(legend.position = "None",
  689. panel.grid.major = element_blank(),
  690. panel.grid.minor = element_blank())
  691. p5 <-
  692. ggplot() +
  693. geom_point(data = param_estimates %>% filter(term == "winstay"), aes(x = trial_bin*10, y = estimate, color = group_verbal),size=.5) +
  694. # geom_line(data = new_data, aes(x = trial_bin, y = predicted), color = "red") +
  695. geom_ribbon(data = param_estimates %>% filter(term == "winstay"),
  696. aes(x = trial_bin*10, ymin = estimate - std.error, ymax = estimate + std.error, fill = group_verbal),
  697. alpha = 0.2) +
  698. geom_errorbar(data = param_estimates %>% filter(term == "winstay"),
  699. aes(x = trial_bin*10, ymin = estimate - std.error, ymax = estimate + std.error, color = group_verbal),
  700. width = .5,size=.5) +
  701. geom_smooth(data = param_estimates %>% filter(term == "winstay"),
  702. aes(x = trial_bin*10, y = estimate, color = group_verbal),
  703. method = "lm",
  704. # formula = y ~ s(x, bs = "cs"),
  705. se = F,size=.5) +
  706. # facet_wrap(~ term, scales = "free_y", ncol = 3, nrow = 3) +
  707. labs(title = "",
  708. x = "Trial",
  709. y = "Win-stay",
  710. color="",
  711. fill="") +
  712. ylim(-.8,.8)+
  713. theme_bw()+
  714. scale_x_continuous(breaks = seq(11, 200,30),limits = c(5,205)) +
  715. scale_color_manual(values = c("ASD" = "orange", "ASD (low verbal)" = "red", "TD" = "blue")) +
  716. scale_fill_manual(values = c("ASD" = "orange", "ASD (low verbal)" = "red", "TD" = "blue"))+
  717. theme(legend.position = "None",
  718. panel.grid.major = element_blank(),
  719. panel.grid.minor = element_blank())
  720. p6 <-
  721. ggplot() +
  722. geom_point(data = param_estimates %>% filter(term == "loseswitch"), aes(x = trial_bin*10, y = estimate, color = group_verbal),size=.5) +
  723. # geom_line(data = new_data, aes(x = trial_bin, y = predicted), color = "red") +
  724. geom_ribbon(data = param_estimates %>% filter(term == "loseswitch"),
  725. aes(x = trial_bin*10, ymin = estimate - std.error, ymax = estimate + std.error, fill = group_verbal),
  726. alpha = 0.2) +
  727. geom_errorbar(data = param_estimates %>% filter(term == "loseswitch"),
  728. aes(x = trial_bin*10, ymin = estimate - std.error, ymax = estimate + std.error, color = group_verbal),
  729. width = .5,size=.5) +
  730. geom_smooth(data = param_estimates %>% filter(term == "loseswitch"),
  731. aes(x = trial_bin*10, y = estimate, color = group_verbal),
  732. method = "lm",
  733. # formula = y ~ s(x, bs = "cs"),
  734. se = F,size=.5) +
  735. # facet_wrap(~ term, scales = "free_y", ncol = 3, nrow = 3) +
  736. labs(title = "",
  737. x = "Trial",
  738. y = "Lose-switch",
  739. color="",
  740. fill="") +
  741. ylim(-.8,1.8)+
  742. theme_bw()+
  743. scale_x_continuous(breaks = seq(11, 200,30),limits = c(5,205)) +
  744. scale_color_manual(values = c("ASD" = "orange", "ASD (low verbal)" = "red", "TD" = "blue")) +
  745. scale_fill_manual(values = c("ASD" = "orange", "ASD (low verbal)" = "red", "TD" = "blue"))+
  746. theme(legend.position = "None",
  747. panel.grid.major = element_blank(),
  748. panel.grid.minor = element_blank())
  749. ggarrange(p1,p2,p3,p4,p5,p6,labels = c("A","B","C","D","E","F"),ncol=3,nrow=3,common.legend = TRUE,align = "v") %>% ggexport(filename = "results/figures/geodems/suppl_figure3.pdf")
  750. # remove variables, functions and objects
  751. rm(df,df_group_avg,df_subject_avg,dflash_data,fits,results)
  752. rm(param_estimates,fit_bin_glm,model_fits)
  753. rm(p1,p2,p3,p4,p5,p6)
  754. rm(learning_curve,fit_learning_curve)
  755. ########### Figure 3 (flash weights) ######################
  756. df = alldata[game_version == "ft" & srs_normal == 1 & above_crit == 1 & group_verbal != "flag" & group_verbal != "flag_low_td"]
  757. df = df[,c("subID","group_verbal","trial","choice","dflash","lbin1","lbin2","lbin3","lbin4","lbin5","lbin6","lbin7","lbin8","lbin9","lbin10","rbin1","rbin2","rbin3","rbin4","rbin5","rbin6","rbin7","rbin8","rbin9","rbin10")]
  758. df = df[trial %in% 101:200]
  759. # Fit logistic regression model for each group
  760. fit_model <- function(data) {
  761. glm(choice ~ lbin1 + lbin2 + lbin3 + lbin4 + lbin5 + lbin6 + lbin7 + lbin8 + lbin9 + lbin10 +
  762. rbin1 + rbin2 + rbin3 + rbin4 + rbin5 + rbin6 + rbin7 + rbin8 + rbin9 + rbin10,
  763. family = binomial,
  764. data = data)
  765. }
  766. models <- df %>%
  767. group_by(group_verbal) %>%
  768. nest() %>%
  769. mutate(model = map(data, fit_model))
  770. # Extract coefficients
  771. coef_data <- models %>%
  772. mutate(coef = map(model, tidy)) %>%
  773. select(group_verbal, coef) %>%
  774. unnest(coef) %>%
  775. filter(term != "(Intercept)") %>%
  776. mutate(
  777. side = ifelse(grepl("^l", term), "left", "right"),
  778. bin = as.numeric(gsub("^[lr]bin", "", term)),
  779. weight = estimate * ifelse(side == "left", -1, 1)
  780. ) %>%
  781. group_by(group_verbal, bin) %>%
  782. reframe(avg_weight = mean(weight),se_avg_weight = sqrt(sum(std.error^2)) / 2) # Using propagation of error method)
  783. model <- aov(avg_weight ~ group_verbal * bin, data = coef_data)
  784. summary(model)
  785. # For group comparisons
  786. emmeans(model, pairwise ~ group_verbal)
  787. # For bin comparisons
  788. emmeans(model, pairwise ~ bin)
  789. # For side comparisons
  790. # emmeans(lm_model, pairwise ~ side)
  791. # For group:bin interaction
  792. emmeans(model, pairwise ~ group_verbal | bin)
  793. # For group:side interaction
  794. # emmeans(lm_model, pairwise ~ group | side)
  795. # pairwise_test <- coef_data %>%
  796. # pairwise_t_test(avg_weight ~ group_verbal, p.adjust.method = "bonferroni")
  797. # pairwise_test <- pairwise_test %>%
  798. # add_xy_position(x = "group_verbal")
  799. # Plot model weights
  800. # p <-
  801. # ggplot(coef_data, aes(x = bin, y = estimate, color = group_verbal,shape=side)) +
  802. # geom_line() +
  803. # geom_point() +
  804. # labs(x = "Bin", y = "Weight", color = "Group") +
  805. # theme_minimal()
  806. coef_data$group_verbal <- factor(coef_data$group_verbal, levels = c("TD","ASD","ASD (low verbal)"),ordered = TRUE)
  807. p1<-
  808. ggplot(coef_data, aes(x = bin/4, y = avg_weight, color = group_verbal),size=.5) +
  809. # geom_ribbon(aes(ymin=weight-std.error,ymax=weight+std.error,fill=group_verbal),alpha=.5) +
  810. # geom_point() +
  811. geom_errorbar(aes(ymin=avg_weight-se_avg_weight,ymax=avg_weight+se_avg_weight,color=group_verbal),size=.5,width=.1) +
  812. geom_smooth(span=1,se=F,size=.5)+
  813. ylim(0.25,1)+
  814. labs(x = "Flash time (s)", y = TeX("$\\beta$"), color = "") +
  815. theme_bw()+
  816. scale_x_continuous(breaks = seq(0, 2.5, .5),limits = c(0,2.8)) +
  817. scale_color_manual(values = c("ASD" = "orange", "ASD (low verbal)" = "red", "TD" = "blue")) +
  818. scale_fill_manual(values = c("ASD" = "orange", "ASD (low verbal)" = "red", "TD" = "blue"))+
  819. theme(legend.position = "None",
  820. panel.grid.major = element_blank(),
  821. panel.grid.minor = element_blank())
  822. # Function to group bins
  823. # group_bins <- function(data) {
  824. # data %>%
  825. # mutate(
  826. # l_early = as.integer(lbin1 + lbin2 + lbin3>0),
  827. # l_middle = as.integer(lbin4 + lbin5 + lbin6 + lbin7>0),
  828. # l_late = as.integer(lbin8 + lbin9 + lbin10>0),
  829. # r_early = as.integer(rbin1 + rbin2 + rbin3>0),
  830. # r_middle = as.integer(rbin4 + rbin5 + rbin6 + rbin7>0),
  831. # r_late = as.integer(rbin8 + rbin9 + rbin10>0)
  832. # ) %>%
  833. # select(-starts_with("lbin"), -starts_with("rbin"))
  834. # }
  835. group_bins <- function(data) {
  836. data %>%
  837. mutate(
  838. l_early = lbin1 + lbin2 + lbin3,
  839. l_middle = lbin4 + lbin5 + lbin6 + lbin7,
  840. l_late = lbin8 + lbin9 + lbin10,
  841. r_early = rbin1 + rbin2 + rbin3,
  842. r_middle = rbin4 + rbin5 + rbin6 + rbin7,
  843. r_late = rbin8 + rbin9 + rbin10
  844. ) %>%
  845. select(-starts_with("lbin"), -starts_with("rbin"))
  846. }
  847. # Fit logistic regression model for each group
  848. fit_model <- function(data) {
  849. data <- group_bins(data)
  850. glm(choice ~ l_early + l_middle + l_late + r_early + r_middle + r_late,
  851. family = binomial,
  852. data = data)
  853. }
  854. models <- df %>%
  855. group_by(group_verbal) %>%
  856. nest() %>%
  857. mutate(model = map(data, fit_model))
  858. # Extract coefficients
  859. coef_data <- models %>%
  860. mutate(coef = map(model, tidy)) %>%
  861. select(group_verbal, coef) %>%
  862. unnest(coef) %>%
  863. filter(term != "(Intercept)") %>%
  864. mutate(
  865. side = ifelse(grepl("^l", term), "left", "right"),
  866. bin = gsub("^[lr]_", "", term),
  867. weight = estimate * ifelse(side == "left", -1, 1)
  868. ) %>%
  869. group_by(group_verbal, bin) %>%
  870. reframe(avg_weight = mean(weight),se_avg_weight = sqrt(sum(std.error^2)) / 2) # Using propagation of error method)
  871. coef_data$bin <- factor(coef_data$bin, levels = c("early", "middle", "late"),ordered = TRUE)
  872. coef_data$bin_num <- as.numeric(coef_data$bin)
  873. p2<-
  874. ggplot(coef_data, aes(x = bin_num, y = avg_weight, color = group_verbal),size=.5) +
  875. # geom_ribbon(aes(ymin=weight-std.error,ymax=weight+std.error,fill=group_verbal),alpha=.5) +
  876. geom_point(size=.5) +
  877. geom_errorbar(aes(ymin=avg_weight-se_avg_weight,ymax=avg_weight+se_avg_weight,color=group_verbal),size=.5,width=.2) +
  878. geom_smooth(span=1,se=F,size=.5)+
  879. ylim(0.25,1)+
  880. labs(x = "Flash count", y = TeX("$\\beta$"), color = "") +
  881. theme_bw()+
  882. scale_x_discrete(limits = c("early", "middle", "late")) +
  883. scale_color_manual(values = c("ASD" = "orange", "ASD (low verbal)" = "red", "TD" = "blue")) +
  884. scale_fill_manual(values = c("ASD" = "orange", "ASD (low verbal)" = "red", "TD" = "blue"))+
  885. theme(legend.position = "None",
  886. panel.grid.major = element_blank(),
  887. panel.grid.minor = element_blank())
  888. ggarrange(p1,p2, nrow=3,ncol = 3,common.legend = TRUE,legend = "top",align = "v",
  889. labels = c("B","C")) %>% ggexport(filename = "results/figures/geodems/flash_weights.pdf")
  890. # remove variables, functions and objects
  891. rm(df,coef_data,models,fit_model)
  892. rm(p1,p2)
  893. rm(group_bins)
  894. ########### Figure 3 (heatmaps) ############################
  895. df = alldata[game_version == "ft" & srs_normal == 1 & above_crit == 1 & group_verbal != "flag" & group_verbal != "flag_low_td"]
  896. df = df[,c("subID","binary_group","group_verbal","trial","score","dflash","tflash","choice","numflashright","numflashleft")]
  897. # Function to create heatmap for a single group
  898. create_heatmap <- function(data, group_name) {
  899. # Aggregate data
  900. agg_data <- data[, .(
  901. accuracy = mean(score),
  902. n = .N
  903. ), by = .(numflashright, numflashleft)]
  904. # Keep only upper triangular matrix (including diagonal)
  905. agg_data <- agg_data[numflashright > numflashleft]
  906. # agg_data <- agg_data[numflashright >0]
  907. # agg_data <- agg_data[numflashleft >0]
  908. # Create heatmap
  909. ggplot(agg_data, aes(x = numflashleft, y = numflashright, fill = accuracy)) +
  910. geom_tile() +
  911. scale_fill_gradient(low = "white", #mid = "white",
  912. high = "mediumblue",
  913. # midpoint = 0.75,
  914. limits = c(0.6, 1),
  915. name = "Accuracy") +
  916. labs(title = paste(group_name),
  917. x = "Smaller Number",
  918. y = "Larger Number") +
  919. theme_bw() +
  920. theme(plot.title = element_text(hjust = 0.5)) +
  921. # geom_text(aes(label = sprintf("%.2f\n(n=%d)", accuracy, n)),
  922. # color = "black", size = 3) +
  923. scale_x_continuous(breaks = seq(0, max(data$numflashleft)-1, by = 1)) +
  924. scale_y_continuous(breaks = seq(1, max(data$numflashright), by = 1)) +
  925. coord_fixed(ratio = 1) +
  926. theme(legend.position = c(.8, .3),
  927. panel.grid.major = element_blank(),
  928. panel.grid.minor = element_blank())
  929. # + # Make the plot square
  930. # geom_abline(intercept = 0, slope = 1, linetype = "dashed", color = "black") # Add diagonal line
  931. }
  932. # Create heatmaps for each group
  933. asd_heatmap <- create_heatmap(df[binary_group == "ASD"], "ASD")
  934. td_heatmap <- create_heatmap(df[binary_group == "TD"], "TD")
  935. # asd_low_heatmap <- create_heatmap(df[group_verbal == "ASD (low verbal)"], "ASD (low verbal)")
  936. # Display heatmaps
  937. print(asd_heatmap)
  938. print(td_heatmap)
  939. # print(asd_low_heatmap)
  940. agg_data_td <- df[binary_group=="TD", .(
  941. accuracy = mean(score),
  942. n = .N
  943. ), by = .(numflashright, numflashleft)]
  944. agg_data_asd <- df[binary_group=="ASD", .(
  945. accuracy = mean(score),
  946. n = .N
  947. ), by = .(numflashright, numflashleft)]
  948. agg_data <- merge(agg_data_td,agg_data_asd,by=c("numflashleft","numflashright"))
  949. agg_data$accuracy = agg_data$accuracy.x - agg_data$accuracy.y
  950. agg_data <- agg_data[numflashright > numflashleft]
  951. td_asd_heatmap <- ggplot(agg_data, aes(x = numflashleft, y = numflashright, fill = accuracy)) +
  952. geom_tile() +
  953. scale_fill_gradient2(low = "red", mid = "white",
  954. high = "mediumblue",
  955. # midpoint = 0.75,
  956. limits = c(-0.15, .15),
  957. name = "Accuracy") +
  958. labs(title = paste("TD - ASD"),
  959. x = "Smaller Number",
  960. y = "Larger Number") +
  961. theme_bw() +
  962. theme(plot.title = element_text(hjust = 0.5)) +
  963. # geom_text(aes(label = sprintf("%.2f\n(n=%d)", accuracy, n)),
  964. # color = "black", size = 3) +
  965. scale_x_continuous(breaks = seq(0, max(agg_data$numflashleft)-1, by = 1)) +
  966. scale_y_continuous(breaks = seq(1, max(agg_data$numflashright), by = 1)) +
  967. coord_fixed(ratio = 1) +
  968. theme(legend.position = c(.8, .3),
  969. panel.grid.major = element_blank(),
  970. panel.grid.minor = element_blank())
  971. # + # Make the plot square
  972. # geom_abline(i
  973. ggarrange(td_heatmap, asd_heatmap, nrow=2,ncol = 3,
  974. common.legend = TRUE,legend = "top",
  975. align = "h",
  976. labels = c("D","E","F")) %>% ggexport(filename = "results/figures/geodems/heatmaps.pdf")
  977. ggarrange(td_asd_heatmap, nrow=2,ncol = 3,
  978. common.legend = TRUE,legend = "right",
  979. align = "h",
  980. labels = c("F")) %>% ggexport(filename = "results/figures/geodems/diff_heatmaps.pdf")
  981. rm(df,asd_heatmap,td_heatmap,asd_low_heatmap,create_heatmap,agg_data,agg_data_asd,agg_data_asd_low,agg_data_td)
  982. ############### Figure 4 #################################
  983. library(optimx)
  984. library(boot)
  985. df = alldata[game_version == "ft" & srs_normal == 1 & above_crit == 1 & group_verbal != "flag" & group_verbal != "flag_low_td"]
  986. df = df[,c("subID","binary_group","group_verbal","trial","score","dflash","tflash","choice","numflashright","numflashleft","winstay","loseswitch")]
  987. df = df[trial %in% 11:200]
  988. df$r <- df$numflashright
  989. df$l <- df$numflashleft
  990. df <- as.data.frame(df[,c("subID","binary_group","choice","r","l","trial","winstay","loseswitch")])
  991. sdt_mdl <- function(params,data){
  992. p_choose_l = rep(0, length(data$choice))
  993. p_choose_r = rep(0, length(data$choice))
  994. if (min(params)<0){
  995. aic<-1000000
  996. }
  997. else{
  998. for (i in 1:length(data$r)){
  999. k0=params[1]
  1000. k1=params[2]
  1001. bs=params[3]
  1002. b1=params[4]
  1003. b2=params[5]
  1004. p=pnorm(0,mean=data$r[i]-data$l[i]+bs+b1*data$winstay[i]+b2*data$loseswitch[i],
  1005. sd=sqrt(k1*(data$r[i]^2+data$l[i]^2)+k0))
  1006. if (length(p)==0){
  1007. p=0
  1008. }
  1009. p_choose_l[i]=p
  1010. p_choose_r[i]<-1-p_choose_l[i]
  1011. }
  1012. lhd<-rep(NA,length(data$choice))
  1013. lhd[data$choice==1]<-p_choose_r[data$choice==1]
  1014. lhd[data$choice==0]<-p_choose_l[data$choice==0]
  1015. lhd[lhd==0]<-.0001
  1016. log.lhd<-log(lhd)
  1017. aic<- -2*sum(log.lhd)+2*length(params) # since there are 5 params
  1018. }
  1019. return(aic)
  1020. }
  1021. fit_model <- function(data) {
  1022. result <- optim(par = c(0.5, 0.5, 0, 0, 0),
  1023. fn = sdt_mdl,
  1024. data = data,
  1025. method = "SANN",
  1026. hessian = TRUE
  1027. )
  1028. }
  1029. # Fit model for ASD group
  1030. asd_data <- df[df$binary_group == "ASD", ]
  1031. asd_model <- fit_model(asd_data)
  1032. asd_se <- sqrt(diag(solve(asd_model$hessian)))
  1033. # Fit model for TD group
  1034. td_data <- df[df$binary_group == "TD", ]
  1035. td_model <- fit_model(td_data)
  1036. td_se <- sqrt(diag(solve(td_model$hessian)))
  1037. params_df <- data.frame(
  1038. group = rep(c("ASD", "TD"), each = 5),
  1039. parameter = rep(c("k0", "k1","bs", "b1", "b2"), 2),
  1040. value = c(asd_model$par, td_model$par),
  1041. se = c(asd_se, td_se)
  1042. )
  1043. params_df$parameter <- factor(params_df$parameter,
  1044. levels = c("k0","k1","bs", "b1", "b2"))
  1045. # p<-
  1046. # ggplot(params_df, aes(x = group, y = value)) +
  1047. # geom_point(aes(color = group), size = 2) +
  1048. # geom_errorbar(aes(ymin = value - se, ymax = value + se,color = group), width = 0.2) +
  1049. # # geom_bar(stat = "identity", position = "dodge") +
  1050. # facet_wrap(~ parameter, scales = "free_y", nrow = 3,ncol=3,shrink = T,dir="v") +
  1051. # theme_bw() +
  1052. # labs(title = "SDT Model Parameters by Group", x = "Group", y = "Value",color="") +
  1053. # theme(axis.text.x = element_text(angle = 45, hjust = 1))+
  1054. # theme(#legend.position = c(.8, .3),
  1055. # panel.grid.major = element_blank(),
  1056. # panel.grid.minor = element_blank())
  1057. # ggexport(p, filename = "results/figures/geodems/sdt_model_params.pdf")
  1058. # Generate predictions for ASD group
  1059. mdl_func <- function(params,data){
  1060. p_choose_l = rep(0, length(data$choice))
  1061. p_choose_r = rep(0, length(data$choice))
  1062. for (i in 1:length(data$r)){
  1063. k0=params[1]
  1064. k1=params[2]
  1065. bs=params[3]
  1066. b1=params[4]
  1067. b2=params[5]
  1068. p=pnorm(0,mean=data$r[i]-data$l[i]+bs+b1*data$winstay[i]+b2*data$loseswitch[i],
  1069. sd=sqrt(k1*(data$r[i]^2+data$l[i]^2)+k0))
  1070. if (length(p)==0){
  1071. p=0
  1072. }
  1073. p_choose_l[i]=p
  1074. p_choose_r[i]<-1-p_choose_l[i]
  1075. }
  1076. return(p_choose_r)
  1077. }
  1078. asd_pred <- mdl_func(asd_model$par, asd_data)
  1079. td_pred <- mdl_func(td_model$par, td_data)
  1080. fit_data <- rbind(
  1081. data.frame(subID = asd_data$subID,group = "ASD", observed = asd_data$choice, predicted = asd_pred, dflash=asd_data$r-asd_data$l),
  1082. data.frame(subID = td_data$subID,group = "TD", observed = td_data$choice, predicted = td_pred, dflash=td_data$r-td_data$l)
  1083. )
  1084. df_sub_avg <- fit_data %>%
  1085. group_by(subID,group, dflash) %>%
  1086. summarise(
  1087. avg_choice = mean(observed, na.rm = TRUE),
  1088. avg_pred = mean(predicted, na.rm=TRUE),
  1089. .groups = "drop")
  1090. df_group_avg <- df_sub_avg %>%
  1091. group_by(group, dflash) %>%
  1092. summarise(
  1093. mean_choice = mean(avg_choice, na.rm = TRUE),
  1094. mean_pred = mean(avg_pred, na.rm = TRUE),
  1095. choice_se = sd(avg_choice, na.rm = TRUE) / sqrt(n()),
  1096. pred_se = sd(avg_pred, na.rm = TRUE) / sqrt(n()),
  1097. .groups = "drop"
  1098. )
  1099. # p_fit <-
  1100. # ggplot(df_group_avg, aes(x=dflash,color = group)) +
  1101. # geom_point(aes(y=mean_choice),alpha = 0.5) +
  1102. # geom_errorbar(aes(y=mean_choice, ymin = mean_choice - choice_se, ymax = mean_choice + choice_se), width = 0.5) +
  1103. # geom_line(aes(y=mean_pred),alpha = 1) +
  1104. # # geom_smooth(method = "glm",method.args = list(family = "binomial") ,se = FALSE) +
  1105. # # facet_wrap(~group,nrow = 3,ncol=3) +
  1106. # theme_bw() +
  1107. # labs(title = "Model Fit vs Original Data", x = "Cue difference", y = "%Went right",color="")+
  1108. # theme(legend.position = "None",
  1109. # panel.grid.major = element_blank(),
  1110. # panel.grid.minor = element_blank())
  1111. # ggarrange(p_fit, nrow=3, ncol=3) %>%
  1112. # ggexport(p_fit, filename = "results/figures/geodems/sdt_model_fit.pdf")
  1113. # rm(p,p_fit)
  1114. setDT(df_group_avg)
  1115. p1 <-
  1116. ggplot(df_group_avg[group=="ASD"],aes(x=dflash))+
  1117. geom_point(aes(y=mean_choice),alpha = 0.5,color="red") +
  1118. geom_errorbar(aes(y=mean_choice, ymin = mean_choice - choice_se, ymax = mean_choice + choice_se), width = 1,color="red") +
  1119. geom_line(aes(y=mean_pred),alpha = 1,color="black") +
  1120. theme_bw() +
  1121. labs(title = "Model Fit: ASD", x = "Cue difference", y = "%Went right",color="")+
  1122. theme(legend.position = "None",
  1123. panel.grid.major = element_blank(),
  1124. panel.grid.minor = element_blank())
  1125. p2 <-
  1126. ggplot(df_group_avg[group=="TD"],aes(x=dflash))+
  1127. geom_point(aes(y=mean_choice),alpha = 0.5,color="blue") +
  1128. geom_errorbar(aes(y=mean_choice, ymin = mean_choice - choice_se, ymax = mean_choice + choice_se), width = 1,color="blue") +
  1129. geom_line(aes(y=mean_pred),alpha = 1,color="black") +
  1130. theme_bw() +
  1131. labs(title = "TD", x = "Cue difference", y = "%Went right",color="")+
  1132. theme(legend.position = "None",
  1133. panel.grid.major = element_blank(),
  1134. panel.grid.minor = element_blank())
  1135. ggarrange(p1,p2, nrow=3,ncol = 3,common.legend = TRUE,legend = "top",align = "h",
  1136. labels = c("E","F")) %>% ggexport(filename = "results/figures/geodems/sdt_model_fit.pdf")
  1137. rm(p1,p2)
  1138. setDT(params_df)
  1139. p1 <-
  1140. ggplot(params_df[parameter=="k0"], aes(x = group, y = value)) +
  1141. geom_point(aes(color = group), size = 2) +
  1142. geom_errorbar(aes(ymin = value - se, ymax = value + se,color = group), width = 0.2) +
  1143. # facet_wrap(~ parameter, scales = "free_y", nrow = 3,ncol=3,shrink = T,dir="v") +
  1144. theme_bw() +
  1145. ylim(0.5,2) +
  1146. labs(title = "", x="", y = TeX("$k_0$"),color="") +
  1147. theme(axis.text.x = element_text(angle = 45, hjust = 1))+
  1148. theme(#legend.position = c(.8, .3),
  1149. panel.grid.major = element_blank(),
  1150. panel.grid.minor = element_blank())
  1151. p2 <-
  1152. ggplot(params_df[parameter=="k1"], aes(x = group, y = value)) +
  1153. geom_point(aes(color = group), size = 2) +
  1154. geom_errorbar(aes(ymin = value - se, ymax = value + se,color = group), width = 0.2) +
  1155. theme_bw() +
  1156. ylim(0.05,0.15) +
  1157. labs(x = "", y = TeX("$k_1$"),color="") +
  1158. theme(axis.text.x = element_text(angle = 45, hjust = 1))+
  1159. theme(
  1160. panel.grid.major = element_blank(),
  1161. panel.grid.minor = element_blank())
  1162. p3 <-
  1163. ggplot(params_df[parameter=="bs"], aes(x = group, y = value)) +
  1164. geom_point(aes(color = group), size = 2) +
  1165. geom_errorbar(aes(ymin = value - se, ymax = value + se,color = group), width = 0.2) +
  1166. theme_bw() +
  1167. ylim(-0.1,0.1) +
  1168. labs(x = "", y = TeX("$bs$"),color="") +
  1169. theme(axis.text.x = element_text(angle = 45, hjust = 1))+
  1170. theme(
  1171. panel.grid.major = element_blank(),
  1172. panel.grid.minor = element_blank())
  1173. p4 <-
  1174. ggplot(params_df[parameter=="b1"], aes(x = group, y = value)) +
  1175. geom_point(aes(color = group), size = 2) +
  1176. geom_errorbar(aes(ymin = value - se, ymax = value + se,color = group), width = 0.2) +
  1177. theme_bw() +
  1178. ylim(0,0.3) +
  1179. labs(x = "", y = TeX("$b1$"),color="") +
  1180. theme(axis.text.x = element_text(angle = 45, hjust = 1))+
  1181. theme(
  1182. panel.grid.major = element_blank(),
  1183. panel.grid.minor = element_blank())
  1184. p5 <-
  1185. ggplot(params_df[parameter=="b2"], aes(x = group, y = value)) +
  1186. geom_point(aes(color = group), size = 2) +
  1187. geom_errorbar(aes(ymin = value - se, ymax = value + se,color = group), width = 0.2) +
  1188. theme_bw() +
  1189. ylim(0.3,0.7) +
  1190. labs(x = "", y = TeX("$b2$"),color="") +
  1191. theme(axis.text.x = element_text(angle = 45, hjust = 1))+
  1192. theme(
  1193. panel.grid.major = element_blank(),
  1194. panel.grid.minor = element_blank())
  1195. ggarrange(p1,p2,p3,p4,p5, nrow=3,ncol = 3,common.legend = TRUE,legend = "top",align = "v",
  1196. labels = c("G","H","I","J","K")) %>% ggexport(filename = "results/figures/geodems/sdt_model_params.pdf")
  1197. ################ Supplementary figure 4 #################################
  1198. sdt_mdl_linear <- function(params,data){
  1199. p_choose_l = rep(0, length(data$choice))
  1200. p_choose_r = rep(0, length(data$choice))
  1201. if (min(params)<0){
  1202. aic<-1000000
  1203. }
  1204. else{
  1205. for (i in 1:length(data$r)){
  1206. k0=params[1]
  1207. k1=params[2]
  1208. bs=params[3]
  1209. b1=params[4]
  1210. b2=params[5]
  1211. p=pnorm(0,mean=data$r[i]-data$l[i]+bs+b1*data$winstay[i]+b2*data$loseswitch[i],
  1212. sd=sqrt(k1*(data$r[i]+data$l[i])+k0))
  1213. if (length(p)==0){
  1214. p=0
  1215. }
  1216. p_choose_l[i]=p
  1217. p_choose_r[i]<-1-p_choose_l[i]
  1218. }
  1219. lhd<-rep(NA,length(data$choice))
  1220. lhd[data$choice==1]<-p_choose_r[data$choice==1]
  1221. lhd[data$choice==0]<-p_choose_l[data$choice==0]
  1222. lhd[lhd==0]<-.0001
  1223. log.lhd<-log(lhd)
  1224. aic<- -2*sum(log.lhd)+2*length(params) # since there are 5 params
  1225. }
  1226. return(aic)
  1227. }
  1228. fit_model_linear <- function(data) {
  1229. result <- optim(par = c(0.5, 0.5, 0, 0, 0),
  1230. fn = sdt_mdl_linear,
  1231. data = data,
  1232. method = "SANN",
  1233. hessian = TRUE
  1234. )
  1235. }
  1236. asd_model_linear <- fit_model_linear(asd_data)
  1237. asd_se_linear <- sqrt(diag(solve(asd_model_linear$hessian)))
  1238. td_model_linear <- fit_model_linear(td_data)
  1239. td_se_linear <- sqrt(diag(solve(td_model_linear$hessian)))
  1240. params_df_linear <- data.frame(
  1241. group = rep(c("ASD", "TD"), each = 5),
  1242. parameter = rep(c("k0", "k1","bs", "b1", "b2"), 2),
  1243. value = c(asd_model_linear$par, td_model_linear$par),
  1244. se = c(asd_se_linear, td_se_linear)
  1245. )
  1246. params_df_linear$parameter <- factor(params_df_linear$parameter,
  1247. levels = c("k0","k1","bs", "b1", "b2"))
  1248. mdl_func_linear <- function(params,data){
  1249. p_choose_l = rep(0, length(data$choice))
  1250. p_choose_r = rep(0, length(data$choice))
  1251. for (i in 1:length(data$r)){
  1252. k0=params[1]
  1253. k1=params[2]
  1254. bs=params[3]
  1255. b1=params[4]
  1256. b2=params[5]
  1257. p=pnorm(0,mean=data$r[i]-data$l[i]+bs+b1*data$winstay[i]+b2*data$loseswitch[i],
  1258. sd=sqrt(k1*(data$r[i]+data$l[i])+k0))
  1259. if (length(p)==0){
  1260. p=0
  1261. }
  1262. p_choose_l[i]=p
  1263. p_choose_r[i]<-1-p_choose_l[i]
  1264. }
  1265. return(p_choose_r)
  1266. }
  1267. asd_pred_linear <- mdl_func_linear(asd_model_linear$par, asd_data)
  1268. td_pred_linear <- mdl_func_linear(td_model_linear$par, td_data)
  1269. fit_data_linear <- rbind(
  1270. data.frame(subID = asd_data$subID,group = "ASD", observed = asd_data$choice, predicted = asd_pred_linear, dflash=asd_data$r-asd_data$l),
  1271. data.frame(subID = td_data$subID,group = "TD", observed = td_data$choice, predicted = td_pred_linear, dflash=td_data$r-td_data$l)
  1272. )
  1273. df_sub_avg_linear <- fit_data_linear %>%
  1274. group_by(subID,group, dflash) %>%
  1275. summarise(
  1276. avg_choice = mean(observed, na.rm = TRUE),
  1277. avg_pred = mean(predicted, na.rm=TRUE),
  1278. .groups = "drop")
  1279. df_group_avg_linear <- df_sub_avg_linear %>%
  1280. group_by(group, dflash) %>%
  1281. summarise(
  1282. mean_choice = mean(avg_choice, na.rm = TRUE),
  1283. mean_pred = mean(avg_pred, na.rm = TRUE),
  1284. choice_se = sd(avg_choice, na.rm = TRUE) / sqrt(n()),
  1285. pred_se = sd(avg_pred, na.rm = TRUE) / sqrt(n()),
  1286. .groups = "drop"
  1287. )
  1288. tmp <- data.frame(
  1289. group = rep(c("ASD", "TD"), each = 2),
  1290. mdl = rep(c("scalar noise","linear noise"), 2),
  1291. aic = c(asd_model$value,
  1292. asd_model_linear$value,
  1293. td_model$value,
  1294. td_model_linear$value)
  1295. )
  1296. p <-
  1297. ggplot(tmp, aes(x = group, y = log(aic),color = mdl)) +
  1298. geom_point(position = position_dodge(width = 0.5),size=.5) +
  1299. theme_bw() +
  1300. ylim(9,11)+
  1301. scale_color_manual(values = c("scalar noise" = "darkgreen", "linear noise" = "maroon")) +
  1302. geom_text(position = position_dodge(width = 1),aes(label = round(log(aic), 3)), vjust = -0.5,size=2) +
  1303. labs(title = "Model comparison: Scalar vs Linear Noise", x = "", y = "log(AIC)",color="") +
  1304. theme(axis.text.x = element_text(angle = 45, hjust = 1))+
  1305. theme(legend.position = "bottom",
  1306. # panel.grid.major = element_blank(),
  1307. # panel.grid.minor = element_blank()
  1308. )
  1309. ggarrange(p, nrow=3, ncol=3,align="h") %>% ggexport(filename = "results/figures/geodems/sdt_model_fit_linear.pdf")
  1310. setDT(params_df_linear)
  1311. p1 <-
  1312. ggplot(params_df_linear[parameter=="k0"], aes(x = group, y = value)) +
  1313. geom_point(aes(color = group), size = 2) +
  1314. geom_errorbar(aes(ymin = value - se, ymax = value + se,color = group), width = 0.2) +
  1315. # facet_wrap(~ parameter, scales = "free_y", nrow = 3,ncol=3,shrink = T,dir="v") +
  1316. theme_bw() +
  1317. ylim(0,1.5) +
  1318. labs(title = "", x="", y = TeX("$k_0$"),color="") +
  1319. theme(axis.text.x = element_text(angle = 45, hjust = 1))+
  1320. theme(#legend.position = c(.8, .3),
  1321. panel.grid.major = element_blank(),
  1322. panel.grid.minor = element_blank())
  1323. p2 <-
  1324. ggplot(params_df_linear[parameter=="k1"], aes(x = group, y = value)) +
  1325. geom_point(aes(color = group), size = 2) +
  1326. geom_errorbar(aes(ymin = value - se, ymax = value + se,color = group), width = 0.2) +
  1327. theme_bw() +
  1328. ylim(0.5,.9) +
  1329. labs(x = "", y = TeX("$k_1$"),color="") +
  1330. theme(axis.text.x = element_text(angle = 45, hjust = 1))+
  1331. theme(
  1332. panel.grid.major = element_blank(),
  1333. panel.grid.minor = element_blank())
  1334. p3 <-
  1335. ggplot(params_df_linear[parameter=="bs"], aes(x = group, y = value)) +
  1336. geom_point(aes(color = group), size = 2) +
  1337. geom_errorbar(aes(ymin = value - se, ymax = value + se,color = group), width = 0.2) +
  1338. theme_bw() +
  1339. ylim(-0.1,0.1) +
  1340. labs(x = "", y = TeX("$bs$"),color="") +
  1341. theme(axis.text.x = element_text(angle = 45, hjust = 1))+
  1342. theme(
  1343. panel.grid.major = element_blank(),
  1344. panel.grid.minor = element_blank())
  1345. p4 <-
  1346. ggplot(params_df_linear[parameter=="b1"], aes(x = group, y = value)) +
  1347. geom_point(aes(color = group), size = 2) +
  1348. geom_errorbar(aes(ymin = value - se, ymax = value + se,color = group), width = 0.2) +
  1349. theme_bw() +
  1350. ylim(0,0.3) +
  1351. labs(x = "", y = TeX("$b1$"),color="") +
  1352. theme(axis.text.x = element_text(angle = 45, hjust = 1))+
  1353. theme(
  1354. panel.grid.major = element_blank(),
  1355. panel.grid.minor = element_blank())
  1356. p5 <-
  1357. ggplot(params_df_linear[parameter=="b2"], aes(x = group, y = value)) +
  1358. geom_point(aes(color = group), size = 2) +
  1359. geom_errorbar(aes(ymin = value - se, ymax = value + se,color = group), width = 0.2) +
  1360. theme_bw() +
  1361. ylim(0.3,0.7) +
  1362. labs(x = "", y = TeX("$b2$"),color="") +
  1363. theme(axis.text.x = element_text(angle = 45, hjust = 1))+
  1364. theme(
  1365. panel.grid.major = element_blank(),
  1366. panel.grid.minor = element_blank())
  1367. ggarrange(p1,p2,p3,p4,p5, nrow=3,ncol = 3,common.legend = TRUE,legend = "top",align = "v") %>% ggexport(filename = "results/figures/geodems/sdt_model_params_linear.pdf")
  1368. # p <-
  1369. # ggplot(params_df_linear, aes(x = group, y = value)) +
  1370. # geom_point(aes(color = group), size = 2) +
  1371. # geom_errorbar(aes(ymin = value - se, ymax = value + se,color = group), width = 0.2) +
  1372. # # geom_bar(stat = "identity", position = "dodge") +
  1373. # facet_wrap(~ parameter, scales = "free_y", nrow = 3,ncol=3,shrink = T,dir="v") +
  1374. # theme_bw() +
  1375. # labs(title = "Linear Noise Model Parameters by Group", x = "Group", y = "Value",color="") +
  1376. # theme(axis.text.x = element_text(angle = 45, hjust = 1))+
  1377. # theme(#legend.position = c(.8, .3),
  1378. # panel.grid.major = element_blank(),
  1379. # panel.grid.minor = element_blank())
  1380. # ggexport(p, filename = "results/figures/geodems/sdt_model_params_linear.pdf")
  1381. # setDT(df_group_avg_linear)
  1382. # setDT(df_group_avg)
  1383. # df_group_avg$mdl <- "scalar noise"
  1384. # df_group_avg_linear$mdl <- "linear noise"
  1385. # tmp <- rbind(df_group_avg,df_group_avg_linear)
  1386. # setDT(tmp)
  1387. # # p_fit_linear <-
  1388. # ggplot(tmp[group=="ASD"],aes(x=dflash,color=mdl)) +
  1389. # geom_point(aes(y=mean_choice),color="red",alpha = 0.5,size=.5) +
  1390. # geom_errorbar(aes(y=mean_choice,ymin = mean_choice - choice_se, ymax = mean_choice + choice_se), color="red",width = 0.5) +
  1391. # geom_line(aes(y=mean_pred),alpha = 1) +
  1392. # theme_bw() +
  1393. # scale_color_manual(values = c("scalar noise" = "darkgreen", "linear noise" = "maroon")) +
  1394. # labs(title = "Scalar vs Linear Noise: ASD", x = "Cue difference", y = "%Went right",color="")+
  1395. # theme(#legend.position = "None",
  1396. # panel.grid.major = element_blank(),
  1397. # panel.grid.minor = element_blank())
  1398. # ggarrange(p, nrow=3, ncol=3,align="h") %>%
  1399. # ggexport(filename = "results/figures/geodems/sdt_model_fit_linear.pdf")
  1400. rm(p_fit_linear,p_fit,params_df_linear,asd_model_linear,asd_se_linear,td_model_linear,td_se_linear,fit_data_linear,df_sub_avg_linear,df_group_avg_linear)
  1401. rm(asd_data,td_data,asd_model,asd_se,td_model,td_se,fit_data,df_sub_avg,df_group_avg)
  1402. rm(tmp,p,fit_model_linear,sdt_mdl_linear,params_df_linear,mdl_func_linear,asd_pred_linear,td_pred_linear)
  1403. rm(fit_model,sdt_mdl,params_df,mdl_func,asd_pred,td_pred,fit_model_linear,sdt_mdl_linear,params_df_linear,mdl_func_linear,asd_pred_linear,td_pred_linear)
  1404. rm(df,df_avg)
  1405. ### Figure 5 ########
  1406. library(rstan)
  1407. load("results/files/params_expdata_svnarrow.RData")
  1408. asdhbmfit <- rstan::extract(fitasd)
  1409. tdhbmfit <- rstan::extract(fittd)
  1410. params_df <- data.frame(group=rep(c("ASD","TD"),each=length(asdhbmfit$mu0)),
  1411. mu0=c(asdhbmfit$mu0,tdhbmfit$mu0),
  1412. mu1=c(asdhbmfit$mu1,tdhbmfit$mu1),
  1413. sig0=c(asdhbmfit$sig0,tdhbmfit$sig0),
  1414. sig1=c(asdhbmfit$sig1,tdhbmfit$sig1),
  1415. k0=c(asdhbmfit$k0,tdhbmfit$k0))
  1416. p1<-ggplot(params_df, aes(color = group, x = k0,fill=group)) +
  1417. geom_histogram(aes(y=after_stat(density)), alpha=0.5,
  1418. position="identity")+
  1419. labs(fill="",color="")+
  1420. xlab(TeX("$k_{0}$"))+
  1421. geom_density(alpha=.2,show.legend = F)+
  1422. theme_bw()+
  1423. theme(legend.position =c(.8,.7),
  1424. panel.grid.major = element_blank(),
  1425. panel.grid.minor = element_blank())
  1426. p2<-ggplot(params_df, aes(color = group, x = mu0,fill=group)) +
  1427. geom_histogram(aes(y=after_stat(density)), alpha=0.5,
  1428. position="identity")+
  1429. labs(fill="",color="")+
  1430. xlab(TeX("$\\mu_{0}$"))+
  1431. geom_density(alpha=.2,show.legend = F)+
  1432. theme_bw()+
  1433. theme(legend.position =c(.8,.7),
  1434. panel.grid.major = element_blank(),
  1435. panel.grid.minor = element_blank())
  1436. p3<-ggplot(params_df, aes(color = group, x = sig0,fill=group)) +
  1437. geom_histogram(aes(y=after_stat(density)), alpha=0.5,
  1438. position="identity")+
  1439. labs(fill="",color="")+
  1440. xlab(TeX("$\\sigma_{0}$"))+
  1441. geom_density(alpha=.2,show.legend = F)+
  1442. theme_bw()+
  1443. theme(legend.position =c(.8,.7),
  1444. panel.grid.major = element_blank(),
  1445. panel.grid.minor = element_blank())
  1446. df = alldata[srs_normal == 1 & above_crit == 1 & game_version == "ft"]
  1447. df1 = data.frame("k1"=colMeans(asdhbmfit$k1),
  1448. "b1"=colMeans(asdhbmfit$b1),
  1449. "b2"=colMeans(asdhbmfit$b2),
  1450. "bs"=colMeans(asdhbmfit$bs),
  1451. "gr"="ASD")
  1452. df2 = data.frame("k1"=colMeans(tdhbmfit$k1),
  1453. "b1"=colMeans(tdhbmfit$b1),
  1454. "b2"=colMeans(tdhbmfit$b2),
  1455. "bs"=colMeans(tdhbmfit$bs),
  1456. "gr"="TD")
  1457. dt1 = df[binary_group=="ASD"]
  1458. dt2 = df[binary_group=="TD"]
  1459. dt1$subj_id<-as.numeric(factor(dt1$username,
  1460. levels=unique(dt1$username)))
  1461. dt2$subj_id<-as.numeric(factor(dt2$username,
  1462. levels=unique(dt2$username)))
  1463. df1$subj_id = as.numeric(c(1:length(unique(dt1$subj_id))))
  1464. df2$subj_id = as.numeric(c(1:length(unique(dt2$subj_id))))
  1465. dt1 = merge(dt1,df1,by="subj_id")
  1466. dt2 = merge(dt2,df2,by="subj_id")
  1467. dt = rbind(dt1,dt2)
  1468. dt = unique(dt[,.(username,k1,b1,b2,bs,binary_group)])
  1469. rm(dt1,dt2,df1,df2,df)
  1470. params_df2 <- dt
  1471. colnames(params_df2)[colnames(params_df2)=="binary_group"] <- "group"
  1472. # p4<-
  1473. ggplot(params_df2, aes(color = group, y = k1, x = group)) +
  1474. geom_boxplot(outlier.shape = NA,show.legend = F,alpha=1,width=.5)+
  1475. geom_jitter(width=.1,show.legend = FALSE,size=.5)+
  1476. labs(color="")+xlab("")+
  1477. ylab(TeX("$k_{1}$"))+
  1478. theme_bw()+
  1479. theme(legend.position ="None",
  1480. panel.grid.major = element_blank(),
  1481. panel.grid.minor = element_blank())
  1482. p5<-ggplot(params_df, aes(color = group, x = mu1,fill=group)) +
  1483. geom_histogram(aes(y=after_stat(density)), alpha=0.5,
  1484. position="identity")+
  1485. labs(fill="",color="")+
  1486. xlab(TeX("$\\mu_{1}$"))+
  1487. geom_density(alpha=.2,show.legend = F)+
  1488. theme_bw()+
  1489. theme(legend.position =c(.8,.7),
  1490. panel.grid.major = element_blank(),
  1491. panel.grid.minor = element_blank())
  1492. p6<-ggplot(params_df, aes(color = group, x = sig1,fill=group)) +
  1493. geom_histogram(aes(y=after_stat(density)), alpha=0.5,
  1494. position="identity")+
  1495. labs(fill="",color="")+
  1496. xlab(TeX("$\\sigma_{1}$"))+
  1497. geom_density(alpha=.2,show.legend = F)+
  1498. theme_bw()+
  1499. theme(legend.position =c(.8,.7),
  1500. panel.grid.major = element_blank(),
  1501. panel.grid.minor = element_blank())
  1502. p7<-
  1503. ggplot(params_df2, aes(color = group, y = abs(bs), x = group)) +
  1504. geom_boxplot(outlier.shape = NA,show.legend = F,alpha=1,width=.5)+
  1505. geom_jitter(width=.1,show.legend = FALSE,size=.5)+
  1506. labs(color="")+xlab("")+
  1507. ylab(TeX("Side-bias"))+
  1508. theme_bw()+
  1509. theme(legend.position ="None",
  1510. panel.grid.major = element_blank(),
  1511. panel.grid.minor = element_blank())
  1512. p8<-
  1513. ggplot(params_df2, aes(color = group, y = abs(b1), x = group)) +
  1514. geom_boxplot(outlier.shape = NA,show.legend = F,alpha=1,width=.5)+
  1515. geom_jitter(width=.1,show.legend = FALSE,size=.5)+
  1516. labs(color="")+xlab("")+
  1517. ylab(TeX("WS"))+
  1518. theme_bw()+
  1519. theme(legend.position ="None",
  1520. panel.grid.major = element_blank(),
  1521. panel.grid.minor = element_blank())
  1522. p9<-
  1523. ggplot(params_df2, aes(color = group, y = abs(b2), x = group)) +
  1524. geom_boxplot(outlier.shape = NA,show.legend = F,alpha=1,width=.5)+
  1525. geom_jitter(width=.1,show.legend = FALSE,size=.5)+
  1526. labs(color="")+xlab("")+
  1527. ylab(TeX("LS"))+
  1528. theme_bw()+
  1529. theme(legend.position ="None",
  1530. panel.grid.major = element_blank(),
  1531. panel.grid.minor = element_blank())
  1532. ggarrange(p1,p2,p3,p4,p5,p6,p7,p8,p9, nrow=3,ncol = 3,common.legend = TRUE,legend = "top",align = "v") %>%
  1533. ggexport(filename = "results/figures/geodems/hbm_params.pdf")
  1534. t.test(params_df2$k1~params_df2$group,var.equal = F)
  1535. t.test(abs(params_df2$bs)~params_df2$group)
  1536. t.test(abs(params_df2$b1)~params_df2$group)
  1537. t.test(abs(params_df2$b2)~params_df2$group)
  1538. summary2 = summary[game_version=="ft"]
  1539. # rename column binary_group to group
  1540. colnames(summary2)[colnames(summary2)=="binary_group"] <- "group"
  1541. summary2 = merge(summary2,params_df2,by=c("username","group"),all=TRUE)
  1542. colnames(summary2)
  1543. head(summary2)
  1544. df = alldata[above_crit == 1 & srs_normal == 1 & game_version == "ft"] %>%
  1545. group_by(username) %>% summarise(mean_rt = mean(rt,na.rm = TRUE),
  1546. mean_iti = mean(iti,na.rm = TRUE),
  1547. .groups = "drop")
  1548. summary2 = merge(summary2,df,by="username",all=TRUE)
  1549. rm(df,dt)
  1550. rm(p1,p2,p3,p4,p5,p6,p7,p8,p9,params_df,params_df2,asdhbmfit,tdhbmfit)
  1551. rm(asd_subid,td_subid,fitasd,fittd)
  1552. ###################### individual glms ###############################################
  1553. # Fit a GLM model to the choice data for each subject
  1554. df <- alldata[trial>10 & srs_normal == 1 & above_crit == 1 & game_version == "ft",c("username","subID","choice","dflash","winstay","loseswitch","binary_group","group_verbal")]
  1555. # Function to fit a GLM model to the choice data for a single subject
  1556. fit_subject_glm <- function(subject_data) {
  1557. glm(choice ~ dflash + winstay + loseswitch,
  1558. data = subject_data,
  1559. family = binomial(link = "logit"))
  1560. }
  1561. # fit the model for each subject
  1562. model_fits <- df %>%
  1563. group_by(subID) %>%
  1564. nest() %>%
  1565. mutate(model = map(data, fit_subject_glm))
  1566. # get the parameter estimates for each subject
  1567. param_estimates <- model_fits %>%
  1568. mutate(params = map(model, broom::tidy)) %>%
  1569. unnest(params) %>%
  1570. select(subID, term, estimate) %>%
  1571. pivot_wider(names_from = term, values_from = estimate)
  1572. # Merge with subject group information
  1573. param_estimates <- param_estimates %>%
  1574. left_join(df %>% distinct(subID, binary_group), by = "subID")
  1575. # # Plot the parameter estimates by group
  1576. # param_long <- param_estimates %>%
  1577. # pivot_longer(cols = c(dflash, winstay, loseswitch),
  1578. # names_to = "parameter",
  1579. # values_to = "estimate")
  1580. colnames(param_estimates) <- c("subID","Sidebias","Slope","WS","LS","group")
  1581. p2<-
  1582. ggplot(param_estimates,
  1583. aes(x = group, y = abs(Sidebias), color = group)) +
  1584. geom_boxplot(outlier.shape = NA,show.legend = F,alpha=1,width=.5)+
  1585. geom_jitter(width=.1,show.legend = FALSE,size=.5)+
  1586. labs(x = "", y = "Side-bias") +
  1587. theme_bw()+
  1588. theme(legend.position ="None",
  1589. panel.grid.major = element_blank(),
  1590. panel.grid.minor = element_blank())
  1591. p1<-
  1592. ggplot(param_estimates,
  1593. aes(x = group, y = Slope, color = group)) +
  1594. geom_boxplot(outlier.shape = NA,show.legend = F,alpha=1,width=.5)+
  1595. geom_jitter(width=.1,show.legend = FALSE,size=.5)+
  1596. labs(x = "", y = "Slope") +
  1597. theme_bw()+
  1598. theme(legend.position ="None",
  1599. panel.grid.major = element_blank(),
  1600. panel.grid.minor = element_blank())
  1601. p3<-
  1602. ggplot(param_estimates,
  1603. aes(x = group, y = abs(WS), color = group)) +
  1604. geom_boxplot(outlier.shape = NA,show.legend = F,alpha=1,width=.5)+
  1605. geom_jitter(width=.1,show.legend = FALSE,size=.5)+
  1606. labs(x = "", y = "WS") +
  1607. theme_bw()+
  1608. theme(legend.position ="None",
  1609. panel.grid.major = element_blank(),
  1610. panel.grid.minor = element_blank())
  1611. p4<-
  1612. ggplot(param_estimates,
  1613. aes(x = group, y = abs(LS), color = group)) +
  1614. geom_boxplot(outlier.shape = NA,show.legend = F,alpha=1,width=.5)+
  1615. geom_jitter(width=.1,show.legend = FALSE,size=.5)+
  1616. labs(x = "", y = "LS") +
  1617. theme_bw()+
  1618. theme(legend.position ="None",
  1619. panel.grid.major = element_blank(),
  1620. panel.grid.minor = element_blank())
  1621. summary2 = merge(summary2,param_estimates,by=c("subID","group"),all=TRUE)
  1622. colnames(summary2)
  1623. p5 <-
  1624. ggplot(summary2, aes(x = Slope, y = k1, color = group)) +
  1625. geom_point(size=.5) +
  1626. labs(x = "Psychometric slope", y = "Perceptual noise",color="") +
  1627. theme_bw() +
  1628. theme(legend.position ="None",
  1629. panel.grid.major = element_blank(),
  1630. panel.grid.minor = element_blank())
  1631. cor.test(summary2$Slope,summary2$k1)
  1632. cor.test(summary2$k1,abs(summary2$WS))
  1633. cor.test(summary2$k1,abs(summary2$LS))
  1634. cor.test(summary2$k1,abs(summary2$Sidebias))
  1635. ggplot(summary2, aes(x = k1, y = abs(LS), color = group)) +
  1636. geom_point(size=.5) +
  1637. labs(x = "Psychometric slope", y = "Perceptual noise",color="") +
  1638. theme_bw() +
  1639. theme(legend.position ="None",
  1640. panel.grid.major = element_blank(),
  1641. panel.grid.minor = element_blank())
  1642. # p6 <-
  1643. # ggplot(summary2, aes(x = abs(LS), y = abs(b2), color = group)) +
  1644. # geom_point() +
  1645. # # labs(title = "Sidebias vs. Performance", x = "Sidebias", y = "Performance",color="") +
  1646. # theme_bw() +
  1647. # theme(legend.position ="None",
  1648. # panel.grid.major = element_blank(),
  1649. # panel.grid.minor = element_blank())
  1650. ggarrange(p1,p2,p3,p4,p5, nrow=3,ncol = 3,common.legend = TRUE,legend = "top",align = "v",
  1651. labels = c("A","B","C","D","E")) %>%
  1652. ggexport(filename = "results/figures/geodems/indglm_params.pdf")
  1653. t.test(param_estimates$Slope~param_estimates$group)
  1654. t.test(abs(param_estimates$Sidebias)~param_estimates$group)
  1655. t.test(abs(param_estimates$WS)~param_estimates$group)
  1656. t.test(abs(param_estimates$LS)~param_estimates$group)
  1657. rm(p1,p2,p3,p4,p5,param_estimates,model_fits,fit_subject_glm,df)
  1658. # # Merge with the summary data
  1659. # summary <- summary %>%
  1660. # right_join(param_estimates, by = c("subID","binary_group"))
  1661. # # Plot the fitted choice probabilities for selected subjects
  1662. # plot_subject_fit <- function(subject_id, group_data, full_data) {
  1663. # subject_model <- group_data %>%
  1664. # filter(subject_id == !!subject_id) %>%
  1665. # pull(model) %>%
  1666. # .[[1]]
  1667. # subject_data <- full_data %>% filter(subject_id == !!subject_id)
  1668. # new_data <- expand.grid(
  1669. # dflash = seq(min(subject_data$dflash), max(subject_data$dflash), length.out = 100),
  1670. # winstay = mean(subject_data$winstay),
  1671. # loseswitch = mean(subject_data$loseswitch)
  1672. # )
  1673. # new_data$predicted_prob <- predict(subject_model, newdata = new_data, type = "response")
  1674. # ggplot(subject_data, aes(x = dflash, y = choice)) +
  1675. # geom_point(alpha = 0.5) +
  1676. # geom_line(data = new_data, aes(y = predicted_prob), color = "red") +
  1677. # labs(title = paste("Subject", subject_id, "(", group_data$group[1], ")"),
  1678. # x = "dflash", y = "Choice Probability") +
  1679. # theme_minimal()
  1680. # }
  1681. # # Select two subjects from each group
  1682. # asd_subjects <- param_estimates %>% filter(binary_group == "ASD") %>% slice_sample(n = 1) %>% pull(subID)
  1683. # td_subjects <- param_estimates %>% filter(binary_group == "TD") %>% slice_sample(n = 1) %>% pull(subID)
  1684. # # Plot for selected subjects
  1685. # selected_plots <- map(c(asd_subjects, td_subjects),
  1686. # ~plot_subject_fit(.x, model_fits, df))
  1687. # # Arrange plots in a grid
  1688. # grid.arrange(grobs = selected_plots, ncol = 2)
  1689. # rm(model_fits, param_estimates, param_long, selected_plots,plot_subject_fit,asd_subjects,td_subjects,fit_subject_glm)
  1690. #### correlations ####
  1691. library(Hmisc)
  1692. library(tidyverse)
  1693. library(gridExtra)
  1694. # Select the columns of interest
  1695. df1 <- summary2[#binary_group == "ASD"
  1696. ,c("perf","trainperf","Slope","k1","bs","b1","b2","mean_rt","mean_iti")]
  1697. df1$bs <- abs(df1$bs)
  1698. df1$b1 <- abs(df1$b1)
  1699. df1$b2 <- abs(df1$b2)
  1700. df2 <- summary2[,c("vinabcstd","srs2total","bistotal","aasp_low_reg_raw","aasp_sen_seek_raw","aasp_sen_ses_raw","aasp_sen_avoid_raw","group_verbal","group")]
  1701. df <- cbind(df1,df2)
  1702. df <- df[group_verbal != "flag" & group_verbal!="flag_low_td",]
  1703. # cor_test <- cor.test(df$vinabcstd, df$perf)
  1704. df$group_verbal <- factor(df$group_verbal, levels = c("TD", "ASD", "ASD (low verbal)"),ordered = T)
  1705. p1<-
  1706. ggplot(df,aes(x=vinabcstd,y=perf*100))+
  1707. geom_point(aes(color=group_verbal),size=.5)+
  1708. geom_smooth(method = "lm",se = T,color="black",size=.5)+
  1709. stat_cor(method = "pearson",label.y.npc = "bottom")+
  1710. labs(color="",x="VABS-3",y="%Correct")+
  1711. scale_color_manual(values = c("TD" = "blue", "ASD" = "orange", "ASD (low verbal)" = "red"))+
  1712. # annotate("text", x = 75, y = .5,
  1713. # label = paste("r =", round(cor_test$estimate, 2),
  1714. # ", p =", round(cor_test$p.value, 6)),
  1715. # hjust = 1, vjust = 1)+
  1716. theme_bw()+
  1717. theme(legend.position ="None",
  1718. panel.grid.major = element_blank(),
  1719. panel.grid.minor = element_blank())
  1720. p2 <-
  1721. ggplot(df,aes(x=vinabcstd,y=Slope))+
  1722. geom_point(aes(color=group_verbal),size=.5)+
  1723. geom_smooth(method = "lm",se = T,color="black",size=.5)+
  1724. stat_cor(method = "pearson",label.y.npc = "top")+
  1725. labs(color="",x="VABS-3",y="Psychometric Slope")+
  1726. scale_color_manual(values = c("TD" = "blue", "ASD" = "orange", "ASD (low verbal)" = "red"))+
  1727. theme_bw()+
  1728. theme(legend.position ="None",
  1729. panel.grid.major = element_blank(),
  1730. panel.grid.minor = element_blank())
  1731. p3 <-
  1732. ggplot(df,aes(x=vinabcstd,y=k1))+
  1733. geom_point(aes(color=group_verbal),size=.5)+
  1734. geom_smooth(method = "lm",se = T,color="black",size=.5)+
  1735. stat_cor(method = "pearson",label.y.npc = "top")+
  1736. labs(color="",x="VABS-3",y="Perceptual Noise")+
  1737. scale_color_manual(values = c("TD" = "blue", "ASD" = "orange", "ASD (low verbal)" = "red"))+
  1738. theme_bw()+
  1739. theme(legend.position ="None",
  1740. panel.grid.major = element_blank(),
  1741. panel.grid.minor = element_blank())
  1742. p4 <-
  1743. ggplot(df,aes(x=srs2total,y=perf*100))+
  1744. geom_point(aes(color=group_verbal),size=.5)+
  1745. geom_smooth(method = "lm",se = T,color="black",size=.5)+
  1746. stat_cor(method = "pearson",label.y.npc = "bottom")+
  1747. labs(color="",x="SRS-2",y="%Correct")+
  1748. scale_color_manual(values = c("TD" = "blue", "ASD" = "orange", "ASD (low verbal)" = "red"))+
  1749. theme_bw()+
  1750. theme(legend.position ="None",
  1751. panel.grid.major = element_blank(),
  1752. panel.grid.minor = element_blank())
  1753. p5 <-
  1754. ggplot(df,aes(x=srs2total,y=Slope))+
  1755. geom_point(aes(color=group_verbal),size=.5)+
  1756. geom_smooth(method = "lm",se = T,color="black",size=.5)+
  1757. stat_cor(method = "pearson",label.y.npc = "top")+
  1758. labs(color="",x="SRS-2",y="Psychometric Slope")+
  1759. scale_color_manual(values = c("TD" = "blue", "ASD" = "orange", "ASD (low verbal)" = "red"))+
  1760. theme_bw()+
  1761. theme(legend.position ="None",
  1762. panel.grid.major = element_blank(),
  1763. panel.grid.minor = element_blank())
  1764. p6 <-
  1765. ggplot(df,aes(x=srs2total,y=k1))+
  1766. geom_point(aes(color=group_verbal),size=.5)+
  1767. geom_smooth(method = "lm",se = T,color="black",size=.5)+
  1768. stat_cor(method = "pearson",label.y.npc = "top")+
  1769. labs(color="",x="SRS-2",y="Perceptual Noise")+
  1770. scale_color_manual(values = c("TD" = "blue", "ASD" = "orange", "ASD (low verbal)" = "red"))+
  1771. theme_bw()+
  1772. theme(legend.position ="None",
  1773. panel.grid.major = element_blank(),
  1774. panel.grid.minor = element_blank())
  1775. ggarrange(p1,p2,p3,p4,p5,p6, nrow=3,ncol = 3,common.legend = TRUE,legend = "top",align = "v",
  1776. labels = c("A","B","C","D","E","F")) %>%
  1777. ggexport(filename = "results/figures/geodems/correlations.pdf")
  1778. rm(df)
  1779. ###### cross correlation ######
  1780. library(corrplot)
  1781. library(psych)
  1782. library(coin)
  1783. library(reshape2)
  1784. df <- summary2
  1785. # df <- df[complete.cases(df),]
  1786. df1 <- df[#binary_group == "ASD"
  1787. ,c("perf","trainperf","Slope","k1","bs","b1","b2","mean_rt","mean_iti")]
  1788. df1$bs <- abs(df1$bs)
  1789. df1$b1 <- abs(df1$b1)
  1790. df1$b2 <- abs(df1$b2)
  1791. df2 <- df[,c("vinabcstd","srs2total","bistotal","aasp_low_reg_raw","aasp_sen_seek_raw","aasp_sen_ses_raw","aasp_sen_avoid_raw","group_verbal","group")]
  1792. df1 <- df1[,-c("trainperf","mean_rt")]
  1793. df2 <- df2[,-c("group_verbal","group")]
  1794. # df1 <- as.matrix(df1)
  1795. # df2 <- as.matrix(df2)
  1796. impute_mean <- function(df) {
  1797. for (col in names(df)) {
  1798. df[[col]][is.na(df[[col]])] <- median(df[[col]], na.rm = TRUE)
  1799. }
  1800. return(df)
  1801. }
  1802. # Impute missing values in both datasets
  1803. df1_imputed <- impute_mean(df1)
  1804. df2_imputed <- impute_mean(df2)
  1805. colnames(df1_imputed) <- c("Accuracy","Psychometric Slope","Perceptual Noise","Side-bias","WS","LS","ITI")
  1806. colnames(df2_imputed) <- c("VABS-3","SRS-2","BIS","AASP (Low Registration)","AASP (Sensation Seeking)","AASP (Sensory Sensitivity)","AASP (Sensation Avoiding)")
  1807. observed_corr <- corr.test(df1_imputed, df2_imputed, method = "spearman", adjust = "none")
  1808. set.seed(123) # For reproducibility
  1809. permutation_test <- function(A, B, n_perm = 1000) {
  1810. m <- ncol(A)
  1811. n <- ncol(B)
  1812. p_values <- matrix(NA, nrow = m, ncol = n)
  1813. for (i in 1:m) {
  1814. for (j in 1:n) {
  1815. observed_cor <- cor(A[[i]], B[[j]], method = "spearman", use = "complete.obs")
  1816. permuted_cors <- numeric(n_perm)
  1817. for (k in 1:n_perm) {
  1818. permuted_A <- A[sample(nrow(A)), ] # Permute entire rows of A
  1819. permuted_cors[k] <- cor(permuted_A[[i]], B[[j]], method = "spearman", use = "complete.obs")
  1820. }
  1821. p_values[i, j] <- mean(abs(permuted_cors) >= abs(observed_cor))
  1822. }
  1823. }
  1824. return(p_values)
  1825. }
  1826. # Run permutation test
  1827. p_values <- permutation_test(df1_imputed, df2_imputed)
  1828. p_adjusted <- p.adjust(p_values, method = "fdr")
  1829. sig_mask <- p_adjusted < 0.05
  1830. sig_corr <- observed_corr$r * sig_mask
  1831. melted_corr <- melt(sig_corr)
  1832. colnames(melted_corr) <- c("VarA", "VarB", "Correlation")
  1833. # p1 <-
  1834. ggplot(melted_corr[melted_corr$VarB %in% c("VABS-3","SRS-2","BIS"),], aes(x=VarA, VarB, fill = Correlation)) +
  1835. geom_tile(color="black",lwd = .2,
  1836. linetype = 2) +
  1837. scale_fill_gradient2(low = "yellow", high = "red", mid = "white",
  1838. midpoint = 0,
  1839. limit = c(-.35, .35), space = "Lab",
  1840. name="Spearman\nCorrelation") +
  1841. theme_bw() +
  1842. # coord_fixed()+
  1843. theme(axis.text.x = element_text(angle = 90, vjust = 1, hjust = 1)) +
  1844. coord_fixed() +
  1845. labs(x = "", y = "",
  1846. title = "Significant correlations between game measures and survey scores")
  1847. # p2 <-
  1848. ggplot(melted_corr[!(melted_corr$VarB %in% c("VABS-3","SRS-2","BIS")),], aes(x=VarA, VarB, fill = Correlation)) +
  1849. geom_tile(color="black",lwd = .2,
  1850. linetype = 2) +
  1851. scale_fill_gradient2(low = "yellow", high = "red", mid = "white",
  1852. midpoint = 0,
  1853. limit = c(-.35, .35), space = "Lab",
  1854. name="Spearman\nCorrelation") +
  1855. theme_bw() +
  1856. # coord_fixed()+
  1857. theme(axis.text.x = element_text(angle = 90, vjust = 1, hjust = 1)) +
  1858. coord_fixed() +
  1859. labs(x = "", y = "")
  1860. ggarrange(p1,p2, nrow=2,ncol = 2,common.legend = TRUE,legend = "top",align = "h",labels = c("G")) %>%
  1861. ggexport(filename = "results/figures/geodems/cross_correlations.pdf")
  1862. ############################################################################################################
  1863. # Load necessary libraries
  1864. # library(tidyverse)
  1865. # library(FactoMineR)
  1866. # library(factoextra)
  1867. # install.packages("leaps", type = "source")
  1868. # Assuming your data frame is called 'df'
  1869. # and has a column 'group' with values 'ASD' or 'TD'
  1870. data <- summary2[,c("group","srs2total","vinabcstd","bistotal","aasp_low_reg_raw","aasp_sen_seek_raw","aasp_sen_ses_raw","aasp_sen_avoid_raw",
  1871. "perf","trainperf","Slope","k1",
  1872. # "Sidebias","WS","LS",
  1873. "bs","b1","b2",
  1874. "mean_iti"
  1875. # ,"mean_rt"
  1876. )]
  1877. data <- data[complete.cases(data),]
  1878. # data$WS <- abs(data$WS)
  1879. # data$LS <- abs(data$LS)
  1880. # data$Sidebias <- abs(data$Sidebias)
  1881. data$bs <- abs(data$bs)
  1882. data$b1 <- abs(data$b1)
  1883. data$b2 <- abs(data$b2)
  1884. setDT(data)
  1885. colnames(data) <- c("group","SRS-2","VABS-3","BIS","AASP (Low Registration)","AASP (Sensation Seeking)","AASP (Sensory Sensitivity)","AASP (Sensation Avoiding)",
  1886. "Accuracy","Training Accuracy","Psychometric Slope","Perceptual Noise","Side-bias","WS","LS","Mean ITI"
  1887. # ,"Mean RT"
  1888. )
  1889. # Separate the group column and features
  1890. features <- data[,-c("group")]
  1891. scaled_features <- scale(features)
  1892. condition <- data$group
  1893. # 2. Perform PCA
  1894. pca_result <- prcomp(features, #center = TRUE,
  1895. scale. = TRUE)
  1896. # 3. Summary of PCA results
  1897. summary(pca_result)
  1898. print(pca_result$rotation)
  1899. # 4. Calculate variance explained by each component
  1900. var_explained <- pca_result$sdev^2 / sum(pca_result$sdev^2)
  1901. scree_data <- data.frame(PC = 1:length(var_explained),
  1902. VarExplained = var_explained)
  1903. p1<-
  1904. ggplot(scree_data, aes(x = PC*100, y = VarExplained)) +
  1905. geom_line() +
  1906. geom_point() +
  1907. theme_bw() +
  1908. labs(x = "Principal Component",
  1909. y = "%Variance Explained",
  1910. title = "Scree Plot") +
  1911. scale_x_continuous(breaks = 1:length(var_explained))
  1912. pc_scores <- pca_result$x
  1913. pca_data <- data.frame(pc_scores, condition = condition)
  1914. p2 <-
  1915. ggplot(pca_data, aes(x = PC1, y = PC2, color = condition)) +
  1916. geom_point(alpha = 0.7,size=.5) +
  1917. theme_bw() +
  1918. labs(x = "PC1", y = "PC2", color="",
  1919. title = "PC1 vs PC2 colored by group") +
  1920. scale_color_manual(values = c("ASD" = "red", "TD" = "blue")) +
  1921. # stat_ellipse(level = 0.95)+
  1922. theme(legend.position =c(.8,.8),
  1923. panel.grid.major = element_blank(),
  1924. panel.grid.minor = element_blank())
  1925. loadings <- data.frame(pca_result$rotation)
  1926. # loadings$feature <- rownames(loadings)
  1927. loadings_sorted <- loadings %>%
  1928. arrange(desc(abs(PC1)))
  1929. # loadings_sorted$PC1<- factor(loadings_sorted$PC1, levels = loadings_sorted$PC1)
  1930. loadings_sorted$feature <- rownames(loadings_sorted)
  1931. loadings_sorted$feature <- factor(loadings_sorted$feature,levels = loadings_sorted$feature)
  1932. # loadings_long <- tidyr::pivot_longer(loadings,
  1933. # cols = c(PC1, PC2),
  1934. # names_to = "PC",
  1935. # values_to = "loading")
  1936. loadings_long <- loadings_sorted %>%
  1937. pivot_longer(cols = c(PC1, PC2), names_to = "PC", values_to = "loading")
  1938. # ggplot(loadings_long, aes(x = feature, y = Loading, fill = PC)) +
  1939. # geom_col(position = "dodge") +
  1940. # coord_flip() +
  1941. # labs(x = "Variables", y = "Loadings") +
  1942. # theme_minimal()
  1943. p3 <-
  1944. ggplot(loadings_long, aes(x = feature, y = abs(loading), fill = PC)) +
  1945. geom_bar(stat = "identity", position = "dodge") +
  1946. theme_bw() +
  1947. labs(x = "Features", y = "Loading", fill="",
  1948. title = "Loadings for PC1 and PC2") +
  1949. theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  1950. scale_fill_manual(values = c("PC1" = "#F8766D", "PC2" = "#00BFC4"))+
  1951. theme(legend.position =c(.8,.8),
  1952. panel.grid.major = element_blank(),
  1953. panel.grid.minor = element_blank())
  1954. biplot_data <- data.frame(pca_result$x[,1:2])
  1955. biplot_data$condition <- condition
  1956. loadings_scaled <- pca_result$rotation[,1:2] *
  1957. (max(abs(pca_result$x[,1:2])) /
  1958. max(abs(pca_result$rotation[,1:2])))
  1959. p4 <-
  1960. ggplot(biplot_data, aes(x = PC1, y = PC2, color = condition)) +
  1961. geom_point(alpha = 0.7,size=.5) +
  1962. geom_segment(data = data.frame(loadings_scaled),
  1963. aes(x = 0, y = 0, xend = PC1, yend = PC2),
  1964. arrow = arrow(length = unit(0.1, "cm")),
  1965. linewidth = 0.5,
  1966. color = "black") +
  1967. xlim(-10,10) + ylim(-10,10) +
  1968. geom_text(data = data.frame(loadings_scaled),
  1969. aes(x = PC1, y = PC2, label = rownames(loadings_scaled)),
  1970. color = "black", vjust = 1, hjust = 1,size=2) +
  1971. theme_bw() +
  1972. labs(x = "PC1", y = "PC2", title = "Biplot",color="") +
  1973. scale_color_manual(values = c("ASD" = "red", "TD" = "blue")) +
  1974. theme(legend.position ="None",
  1975. panel.grid.major = element_blank(),
  1976. panel.grid.minor = element_blank())
  1977. t_test_pc1 <- t.test(abs(pc_scores[,1]) ~ condition)
  1978. print("T-test results for PC1:")
  1979. print(t_test_pc1)
  1980. t_test_pc2 <- t.test(abs(pc_scores[,2]) ~ condition)
  1981. print("T-test results for PC2:")
  1982. print(t_test_pc2)
  1983. ggarrange(p1,p2,p3, nrow=2,ncol = 2,common.legend = FALSE,legend = "top",
  1984. align = "v",labels = c("A","B","C")) %>%
  1985. ggexport(filename = "results/figures/geodems/pca.pdf")
  1986. rm(p1,p2,p3,p4,scree_data,pc_scores,pca_data,loadings,loadings_long,biplot_data,loadings_scaled)
  1987. rm(data,features,condition,pca_result,var_explained)
  1988. ##### classification analysis #####
  1989. library(caret)
  1990. library(pROC)
  1991. library(reshape2)
  1992. set.seed(123)
  1993. data <- summary2[,c("group","perf","trainperf","Slope","k1","mean_iti",
  1994. "bs","b1","b2",
  1995. # "WS","LS","Sidebias",
  1996. "aasp_low_reg_raw",
  1997. "aasp_sen_seek_raw","aasp_sen_ses_raw","aasp_sen_avoid_raw"
  1998. )]
  1999. data <- data[complete.cases(data),]
  2000. data$group <- as.factor(data$group)
  2001. data$bs <- abs(data$bs)
  2002. data$b1 <- abs(data$b1)
  2003. data$b2 <- abs(data$b2)
  2004. # data$WS <- abs(data$WS)
  2005. # data$LS <- abs(data$b1)
  2006. # data$b2 <- abs(data$b2)
  2007. data2 <- data[,c("group","perf","trainperf","Slope","k1","mean_iti",
  2008. # "WS","LS","Sidebias"
  2009. "bs","b1","b2"
  2010. )]
  2011. data3 <- data[,c("group",
  2012. "aasp_low_reg_raw",
  2013. "aasp_sen_seek_raw","aasp_sen_ses_raw","aasp_sen_avoid_raw"
  2014. )]
  2015. ctrl <- trainControl(method = "cv", number = 5, classProbs = TRUE, summaryFunction = twoClassSummary,savePredictions = "all")
  2016. mdl_name = "rf"
  2017. # "rf" - Random Forest
  2018. # "svmRadial" - Support Vector Machines with Radial Basis Function Kernel
  2019. # "glm" - Generalized Linear Model (including Logistic Regression)
  2020. # "knn" - k-Nearest Neighbors
  2021. # "nb" - Naive Bayes
  2022. # "nnet" - Neural Network
  2023. # "gbm" - Gradient Boosting Machine
  2024. # "rpart" - Decision Trees
  2025. # "lda" - Linear Discriminant Analysis
  2026. # "glmnet" - Regularized Generalized Linear Models
  2027. mdl <- train(group ~ .,
  2028. data = data,
  2029. method = mdl_name,
  2030. trControl = ctrl,
  2031. metric = "ROC")
  2032. mdl2 <- train(group ~ .,
  2033. data = data2,
  2034. method = mdl_name,
  2035. trControl = ctrl,
  2036. metric = "ROC")
  2037. mdl3 <- train(group ~ .,
  2038. data = data3,
  2039. method = mdl_name,
  2040. trControl = ctrl,
  2041. metric = "ROC")
  2042. roc_obj <- roc(mdl$pred$obs, mdl$pred$ASD)
  2043. roc_obj2 <- roc(mdl2$pred$obs, mdl2$pred$ASD)
  2044. roc_obj3 <- roc(mdl3$pred$obs, mdl3$pred$ASD)
  2045. roc_data1 = data.frame(
  2046. specificity = roc_obj$specificities,
  2047. sensitivity = roc_obj$sensitivities
  2048. )
  2049. # list(
  2050. roc_data2 = data.frame(
  2051. specificity = roc_obj2$specificities,
  2052. sensitivity = roc_obj2$sensitivities
  2053. )
  2054. # auc = auc(roc_obj2)
  2055. # )
  2056. # list(
  2057. roc_data3 = data.frame(
  2058. specificity = roc_obj3$specificities,
  2059. sensitivity = roc_obj3$sensitivities
  2060. )#,
  2061. # auc = auc(roc_obj3)
  2062. # )
  2063. # Combine ROC data
  2064. roc_data_all <- rbind(
  2065. cbind(roc_data1, set = "Game + AASP"),
  2066. cbind(roc_data2, set = "Game"),
  2067. cbind(roc_data3, set = "AASP")
  2068. )
  2069. aucs <- c(auc(roc_obj),auc(roc_obj2),auc(roc_obj3))
  2070. print(aucs)
  2071. # Plot ROC curves
  2072. p4<- ggplot(roc_data_all, aes(x = 1 - specificity, y = sensitivity, color = set)) +
  2073. geom_line(size = 1) +
  2074. geom_abline(intercept = 0, slope = 1, linetype = "dashed", color = "gray") +
  2075. coord_equal() +
  2076. scale_color_brewer(palette = "Set1") +
  2077. labs(x = "False Positive Rate", y = "True Positive Rate",
  2078. # title = "ROC Curves for Different Feature Sets (10-fold CV)",
  2079. color = "") +
  2080. theme_bw() +
  2081. # annotate("text", x = 0.75, y = 0.25,
  2082. # label = sprintf("AUC Set 1: %.3f\nAUC Set 2: %.3f\nAUC Set 3: %.3f",
  2083. # aucs[1], aucs[2], aucs[3]),
  2084. # hjust = 0, vjust = 0, size = 3) +
  2085. theme(legend.position ="bottom",
  2086. panel.grid.major = element_blank(),
  2087. panel.grid.minor = element_blank())
  2088. ggarrange(p4, nrow=2,ncol = 2,common.legend = FALSE,legend = "top",
  2089. align = "hv",labels = c("D")) %>%
  2090. ggexport(filename = "results/figures/geodems/classify.pdf")
  2091. ##### multiple sessions #####
  2092. df1 = summary[game_version=="ft" & above_crit ==1,c("username","subID","binary_group","perf")]
  2093. df2 = summary[game_version!="ft" & above_crit ==1,c("username","subID","binary_group","perf")]
  2094. df2[, subID := gsub("_2$|_3$", "", subID)]
  2095. common_ids = intersect(df1$subID, df2$subID)
  2096. df1 = df1[subID %in% common_ids]
  2097. df2 = df2[subID %in% common_ids]
  2098. df2 = df2[!duplicated(df2$subID),]
  2099. df2$perf2 <- df2$perf
  2100. subs <- union(df1$username,df2$username)
  2101. df1 <- df1[,-c("username")]
  2102. df2 <- df2[,-c("username","perf")]
  2103. df <- merge(df1,df2,by=c("subID","binary_group"),all=TRUE)
  2104. rm(df1,df2)
  2105. p1 <-
  2106. ggplot(df,aes(x=perf,y=perf2))+
  2107. geom_point(aes(color=binary_group),size=.5)+
  2108. geom_smooth(method = "lm",se = T,color="black",size=.5,linetype="dashed")+
  2109. stat_cor(method = "pearson",label.y.npc = "bottom")+
  2110. labs(title="Accuracy",color="",x="First attempt",y="Second attempt")+
  2111. scale_color_manual(values = c("ASD" = "red", "TD" = "blue"))+
  2112. theme_bw()+
  2113. theme(legend.position ="None",
  2114. panel.grid.major = element_blank(),
  2115. panel.grid.minor = element_blank())
  2116. df = alldata[username %in% subs] %>%
  2117. group_by(username,subID,binary_group,game_version) %>% summarise(mean_rt = mean(rt,na.rm = TRUE),
  2118. .groups = "drop")
  2119. df$mean_rt <- df$mean_rt/1000
  2120. df1 = df[df$game_version=="ft",]
  2121. df2 = df[df$game_version!="ft",]
  2122. setDT(df1)
  2123. setDT(df2)
  2124. df2[, subID := gsub("_2$|_3$", "", subID)]
  2125. # rename column mean_rt to mean_rt2
  2126. colnames(df2)[colnames(df2)=="mean_rt"] <- "mean_rt2"
  2127. df1 <- df1[,-c("username","game_version")]
  2128. df2 <- df2[,-c("username","game_version")]
  2129. # merge the two dataframes on subID and binary_group, remove username
  2130. df = merge(df1,df2,by=c("subID","binary_group"),all=TRUE)
  2131. rm(df1,df2)
  2132. p2 <-
  2133. ggplot(df,aes(x=mean_rt,y=mean_rt2))+
  2134. geom_point(aes(color=binary_group),size=.5)+
  2135. geom_smooth(method = "lm",se = T,color="black",size=.5,linetype="dashed")+
  2136. stat_cor(method = "pearson",label.y.npc = "bottom")+
  2137. labs(title="RT (s)",color="",x="First attempt",y="Second attempt")+
  2138. scale_color_manual(values = c("ASD" = "red", "TD" = "blue"))+
  2139. theme_bw()+
  2140. theme(legend.position ="None",
  2141. panel.grid.major = element_blank(),
  2142. panel.grid.minor = element_blank())
  2143. # make common x and y axis in the grid
  2144. rm(df)
  2145. ggarrange(p1,p2, nrow=3,ncol = 3,common.legend = TRUE,legend = "top",align = "v",labels = c("A")) %>%
  2146. ggexport(filename = "results/figures/geodems/multi_session.pdf")

data_analysis_final.R at commit 906c9be, under CC0-1.0 · at the source

Overview

  1. Psychological & Brain Sciences, Boston University, Boston, MA, USA
  2. Graduate Program for Neuroscience, Boston University, Boston, MA, USA
  3. Center for Systems Neuroscience, Boston University, Boston, MA, USA
Institutions: Boston University (United States)
Journal: Science advances, volume 12, issue 37, article eaec9291
Dates: received 8 October 2025; accepted 3 August 2026; published online 11 September 2026; in print September 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1126/sciadv.aec9291 · PMID 42726881 · PMCID PMC13564830 · OpenAlex W7212241919
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: human (organism), autism (population), cognitive (subfield)
Methods: Connectivity, Statistics, Smoothing, state filtering, decompositions, Machine learning, fMRI & imaging
MeSH: Autism Spectrum Disorder*, Autistic Disorder*, Conditioning, Operant*, Video Games*, Visual Perception*, Adolescent, Animals, Female, Humans, Learning, Male (* major topic)
Topic: Autism Spectrum Disorder Research (Cognitive Neuroscience, Neuroscience), according to OpenAlex
Funding: Simons Foundation (874568)
Citations: not cited yet (Europe PMC); 64 references in the paper

Abstract

Altered perception is a hallmark of autism spectrum disorder (ASD), yet its underlying mechanisms remain unclear, in part because individuals with greater impairments are often excluded from psychophysics research. Here, we introduce GEODE (gathering evidence to optimize decisions), an online video game designed to measure visual perception across the full autism spectrum. GEODE incorporates feedback-based shaping from animal training to teach task mechanics, customized graphics and storyline elements to sustain engagement, and touchscreen compatibility for broader accessibility, enabling participation from adolescents including those with profound autism. Across a large, heterogeneous cohort, autistic participants successfully played GEODE but showed slower learning and reduced accuracy compared with typically developing siblings. These deficits were best explained by increased noise in the integration of sensory evidence, which scaled nonlinearly with stimulus complexity and tracked standardized survey measures of adaptive functioning more closely than diagnostic category alone. Injecting equivalent noise into artificial neural networks produced agents that recapitulated ASD-like patterns of learning and decision-making. These findings implicate noisy evidence integration as a computational mechanism for altered visual perception across the autism spectrum and establish a scalable, accessible platform for probing altered perception and learning in neurodevelopmental and neuropsychiatric disorders.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repositories

Its files are read in the Code ↔ Paper reader above, with 9 matches between paragraphs and lines of code.

Zenodo 16998131

License: CC0-1.0
State: the link answers, verified on 26 September 2026
Evidence: files inventoried
Size: 1 file
Software Heritage: not checked
Found in: “Data, code, and materials availability:”
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: afex (1 file), broom (1 file), caret (1 file), data.table (1 file), emmeans (1 file), ggplot2 (1 file), ggpubr (1 file), glmnet (1 file), lme4 (1 file), pROC (1 file), psych (1 file), reshape2 (1 file), rstatix (1 file), Stan (1 file), tidyverse (1 file)
Availability: 1 check, the latest on 26 September 2026: the link answers (HTTP 200)
  • 26 September 2026: the link answers (HTTP 200)
3 files

Zenodo 20545339

License: CC-BY-4.0
State: the link answers, verified on 26 September 2026
Evidence: files inventoried
Size: 1 file
Software Heritage: not checked
Found in: “Data, code, and materials availability:”
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Availability: 1 check, the latest on 26 September 2026: the link answers (HTTP 200)
  • 26 September 2026: the link answers (HTTP 200)

ratacad/geode_data

License: CC0-1.0
State: the link answers, verified on 26 September 2026
Evidence: files inventoried
Commit: 906c9beeee3eab63dda2abd50f918cc940316ab7, 20 August 2026
Languages: R (9)
Size: 14 files, 9 scripts
Software Heritage: not archived
Found in: the Zenodo archive record
Holds: README, license file
Not found: CITATION.cff, environment file, tests, continuous integration, documentation
Tools: data.table (8 files), tidyverse (5 files), ggplot2 (3 files), ggpubr (3 files), lme4 (3 files), broom (2 files), nlme (2 files), Stan (2 files), afex (1 file), car (1 file), caret (1 file), easystats (1 file), emmeans (1 file), glmnet (1 file), patchwork (1 file), pROC (1 file), psych (1 file), reshape2 (1 file), rstatix (1 file)
Availability: 1 check, the latest on 26 September 2026: the link answers
  • 26 September 2026: the link answers
11 files

ratacad/geode

License: none: the authors keep all their rights
State: the link answers, verified on 26 September 2026
Evidence: files inventoried
Commit: 2a95c90a647fe2e4e920f688e5ed22c5f41541aa, 4 June 2026
Languages: JavaScript (366), TypeScript (8)
Size: 1,353 files, 374 scripts
Software Heritage: not archived
Found in: the Zenodo archive record
Holds: CITATION.cff
Not found: README, license file, environment file, tests, continuous integration, documentation
Availability: 1 check, the latest on 26 September 2026: the link answers
  • 26 September 2026: the link answers
374 files, not copied: shown from their source

OSCR keeps no copy of these files: this repository has no license that allows it. The reader above shows each one from its source, fetched by your browser at commit 2a95c90, when its fingerprint is the one OSCR verified. How this works.

The paper's code and data availability statement is in the Data section.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

  • 4 repositories of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 384 scripts, each with its path and the digest of its content;
  • 9 matches between paragraphs of the paper and lines of the code (method lexical-v1);
  • neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

No dataset and no data link were found in the paper.

Data, code, and materials availability

All data and code needed to evaluate and reproduce the results in the paper are present in the paper and/or the Supplementary Materials. All deidentified data and relevant analysis code can be found here: doi.org/10.5281/zenodo.16998131 (http://dx.doi.org/10.5281/zenodo.16998131). The source code for the game is available here: https://doi.org/10.5281/zenodo.20545339. This study did not generate new materials.

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 3, 28 September 2026

  • Funding: added Simons Foundation Autism Research Initiative: 874568

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 7 authors, 11 MeSH terms, 54 references.

Cite

This paper

Chakravarty, S., Do, Q., Li, Y., Torres-Lacarra, V., Tager-Flusberg, H., McGuire, J. T., & Scott, B. B. (2026). Computational mechanisms of perception in autism revealed using games inspired by rodent operant tasks. Science advances, 12(37), eaec9291. https://doi.org/10.1126/sciadv.aec9291

BibTeX

@article{chakravarty2026computational,
author = {Chakravarty, Sucheta and Do, Quan and Li, Yutong and Torres-Lacarra, Vanessa and Tager-Flusberg, Helen and McGuire, Joseph T and Scott, Benjamin B},
title = {{Computational mechanisms of perception in autism revealed using games inspired by rodent operant tasks}},
journal = {Science advances},
year = {2026},
month = sep,
volume = {12},
number = {37},
pages = {eaec9291},
publisher = {American Association for the Advancement of Science},
issn = {2375-2548},
doi = {10.1126/sciadv.aec9291},
url = {https://doi.org/10.1126/sciadv.aec9291},
pmid = {42726881},
pmcid = {PMC13564830}
}

RIS

TY - JOUR
AU - Chakravarty, Sucheta
AU - Do, Quan
AU - Li, Yutong
AU - Torres-Lacarra, Vanessa
AU - Tager-Flusberg, Helen
AU - McGuire, Joseph T
AU - Scott, Benjamin B
TI - Computational mechanisms of perception in autism revealed using games inspired by rodent operant tasks
T2 - Science advances
J2 - Sci Adv
PY - 2026
DA - 2026/09/11
VL - 12
IS - 37
SP - eaec9291
SN - 2375-2548
PB - American Association for the Advancement of Science
DO - 10.1126/sciadv.aec9291
UR - https://doi.org/10.1126/sciadv.aec9291
LA - en
ER -

CSL-JSON

{
"id": "10.1126/sciadv.aec9291",
"type": "article-journal",
"title": "Computational mechanisms of perception in autism revealed using games inspired by rodent operant tasks",
"container-title": "Science advances",
"author": [
{
"family": "Chakravarty",
"given": "Sucheta"
},
{
"family": "Do",
"given": "Quan"
},
{
"family": "Li",
"given": "Yutong"
},
{
"family": "Torres-Lacarra",
"given": "Vanessa"
},
{
"family": "Tager-Flusberg",
"given": "Helen"
},
{
"family": "McGuire",
"given": "Joseph T"
},
{
"family": "Scott",
"given": "Benjamin B"
}
],
"container-title-short": "Sci Adv",
"volume": "12",
"issue": "37",
"page": "eaec9291",
"DOI": "10.1126/sciadv.aec9291",
"PMID": "42726881",
"PMCID": "PMC13564830",
"ISSN": "2375-2548",
"publisher": "American Association for the Advancement of Science",
"URL": "https://doi.org/10.1126/sciadv.aec9291",
"language": "en",
"issued": {
"date-parts": [
[
2026,
9,
11
]
]
}
}

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1073/pnas.2606871123 [code]
Oxytocin modulates the neurocomputational mechanisms engaged in learning rank relationships in social networks.
Journal: Proceedings of the National Academy of Sciences of the United States of America
In common: nlme, psych, easystats, 10 other tools
[2] doi:10.1038/s42003-026-10957-8 [code]
Brain defence by the extracellular matrix protein Cochlin.
Journal: Communications biology
In common: Stan, glmnet, nlme, 9 other tools
[3] doi:10.1038/s41398-026-04010-9 [code]
Bullying victimization and brain development: a longitudinal structural magnetic resonance imaging study from adolescence to early adulthood.
Journal: Translational psychiatry
In common: nlme, psych, rstatix, 9 other tools
[4] doi:10.1093/bioinformatics/btag592 [code]
Network-based stratification of allele-specific expression reveals patient subgroups in Huntington's disease.
Journal: Bioinformatics (Oxford, England)
In common: pROC, rstatix, car, 9 other tools
[5] doi:10.1038/s41467-026-71415-x [code]
Regional BOLD variability reflects microstructural maturation and neuronal ensheathment in the preterm infant cortex.
Journal: Nature communications
In common: glmnet, nlme, psych, 8 other tools
[6] doi:10.1093/braincomms/fcag146 [code]
Convergent structural brain alterations in chronic pain: a multi-metric individual participant data meta-analysis.
Journal: Brain communications
In common: pROC, glmnet, caret, 8 other tools
[7] doi:10.1073/pnas.2603114123 [code]
The human hippocampus can pattern separate memories by meaning.
Journal: Proceedings of the National Academy of Sciences of the United States of America
In common: afex, psych, rstatix, 8 other tools, cognitive
[8] doi:10.1093/neuonc/noag128 [code]
Spatially-resolved single-cell imaging of melanoma brain metastases identifies localized immune patterns predictive of immune checkpoint blockade response.
Journal: Neuro-oncology
In common: pROC, glmnet, caret, 8 other tools
[9] doi:10.1192/bjp.2026.10664 [code]
Early effects of a novel 5-HT&lt;sub&gt;4&lt;/sub&gt;R agonist (PF-04995274) and the SSRI citalopram on emotional cognition in unmedicated depression: RESTAND study.
Journal: The British journal of psychiatry : the journal of mental science
In common: afex, rstatix, easystats, 8 other tools, cognitive
[10] doi:10.1016/j.celrep.2026.117505 [code]
Impaired spatial coding and neuronal hyperactivity in the medial entorhinal cortex of aged APP knock-in mice.
Journal: Cell reports
In common: Stan, nlme, easystats, 8 other tools

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.