OSCR

Building bridges between brain and behavior: An open-source toolbox for joint modeling with fMRI.

Code ↔ Paper

15 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 15 matches
  1. [1] § Application › Convolving with the HRF ↔ R/MRI.R, lines 158–216 · score 0.79 · hrf_model, convolve_design_matrix, high_pass_model, spm, cosine, timeseries
  2. [2] § Application › Convolving with the HRF ↔ R/MRI.R, lines 158–216 · score 0.76 · convolved design matrix, HRF model, high pass filter, spm, cosine, mri
  3. [3] § Workflow › Design specification ↔ R/MRI.R, lines 380–464 · score 0.74 · cosine basis functions, low frequency drifts, high_pass_model, polynomial, filter
  4. [4] § Cognitive Modeling ↔ R/model_RDM.R, lines 274–358 · score 0.73 · response threshold, evidence accumulation, cognitive models, speeded, drift rate, respond
  5. [5] § Application › Model criticism ↔ R/plot_data.R, lines 1018–1058 · score 0.69 · Defective cumulative distributions, posterior predictive CDF, credible intervals, aggregating
  6. [6] § Workflow › Design specification ↔ R/MRI.R, lines 1131–1252 · score 0.66 · high_pass_model, high pass filtering, poly, design matrix, regressors
  7. [7] § Workflow › Design specification ↔ R/MRI.R, lines 218–269 · score 0.65 · double gamma, hrf_model, Glover, SPM, onset
  8. [8] § Cognitive Modeling ↔ R/model_LBA.R, lines 130–215 · score 0.64 · response threshold, evidence accumulation models, drift rate, respond, latent, stimulus
  9. [9] § Other Considerations › Hypothesis tests for correlations ↔ R/s3_funcs.R, lines 784–810 · score 0.64 · Savage Dickey ratio, Bayes factors, hypothesis, posterior, model
  10. [10] § Application › Model criticism ↔ R/plot_data.R, lines 216–256 · score 0.63 · plot_cdf, defective_factor, post_predict, fun
  11. [11] § Application › Convolving with the HRF ↔ R/MRI.R, lines 567–616 · score 0.63 · design_fmri, design_matrix, MRI_AR1, model
  12. [12] § Other Considerations › Hypothesis tests for correlations ↔ R/bridge_sampling.R, lines 187–276 · score 0.58 · Bridge sampling, Bayes factors, mathematical, selectively, posterior, model
  13. [13] § Application › Convolving with the HRF ↔ R/MRI.R, lines 775–860 · score 0.58 · plot_design_fmri, convolved design matrix, TRs, event
  14. [14] § Application › Prior specification ↔ R/MRI.R, lines 699–770 · score 0.55 · MRI_AR1, standard deviation, rho, transformed, sd, log
  15. [15] § Application › Model estimation ↔ R/design.R, lines 1419–1486 · score 0.54 · sampled_pars, par_names, joint design, grepl, Model

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 · 1,304 lines · 51 KB · GPL-3.0 · 8 matches

  1. apply_contrasts <- function(events, contrast = NULL, cell_coding = FALSE, remove_intercept = TRUE, levels = NULL) {
  2. factor_name <- events$factor[1]
  3. colnames(events)[colnames(events) == "event_type"] <- factor_name
  5. # If a contrast is provided, use it; otherwise let R default to its default contrasts.
  6. if(!is.null(contrast)){
  7. if(is.matrix(contrast)){
  8. if(!is.null(rownames(contrast))){
  9. events[[factor_name]] <- factor(events[[factor_name]], levels = rownames(contrast))
  10. } else {
  11. ## use levels provided
  12. events[[factor_name]] <- factor(events[[factor_name]], levels = levels)
  13. }
  14. stats::contrasts(events[[factor_name]], how.many = ncol(contrast)) <- contrast
  15. } else {
  16. events[[factor_name]] <- factor(events[[factor_name]], levels = levels)
  17. stats::contrasts(events[[factor_name]]) <- do.call(contrast, list(n = length(unique(events[[factor_name]]))))
  18. }
  19. } else {
  20. events[[factor_name]] <- factor(events[[factor_name]], levels = levels)
  21. # R's default contrasts will be used.
  22. }
  24. if(length(unique(events[[factor_name]])) == 1){
  25. design <- matrix(1, nrow = nrow(events))
  26. colnames(design) <- factor_name
  27. } else if(cell_coding){
  28. design <- model.matrix(as.formula(paste0("~ 0 + ", factor_name)), events)
  29. } else {
  30. design <- model.matrix(as.formula(paste0("~ ", factor_name)), events)
  31. colnames(design)[1] <- paste0(factor_name, "0")
  32. if(remove_intercept) design <- design[, -1, drop = FALSE]
  33. }
  34. events$factor <- NULL
  35. events_design <- cbind(design, events)
  37. long_events <- reshape(events_design,
  38. direction = "long",
  39. varying = colnames(design),
  40. v.names = "modulation",
  41. timevar = "regressor",
  42. times = colnames(design))
  43. rownames(long_events) <- NULL
  44. long_events <- long_events[long_events$modulation != 0, ]
  45. long_events <- long_events[, !(colnames(long_events) %in% c("id", factor_name))]
  46. long_events <- long_events[order(long_events$onset), ]
  47. return(long_events)
  48. }
  51. #' Reshape events data for fMRI analysis
  52. #'
  53. #' This function reshapes event data into a format suitable for fMRI analysis by
  54. #' converting specified event_types into separate event types with appropriate modulation values.
  55. #'
  56. #' @param events A data frame containing event information with required columns 'subjects', 'run', and 'onset'
  57. #' @param event_types A character vector of column names in the events data frame to be treated as event_types
  58. #' @param duration Either a single numeric value (applied to all event_types), a list with named elements
  59. #' corresponding to event_types, or a function that takes the events data frame and returns durations
  60. #' @param modulation Either a list with named elements corresponding to event_types, or a function that takes
  61. #' the events data frame and returns durations
  62. #' @return A data frame with columns 'subjects', 'onset', 'run', 'modulation', 'duration', and 'event_type'
  63. #' @export
  64. #' @examples
  65. #' # Create a simple events data frame
  66. #' events <- data.frame(
  67. #' subjects = rep(1, 10),
  68. #' run = rep(1, 10),
  69. #' onset = seq(0, 90, by = 10),
  70. #' condition = rep(c("A", "B"), 5),
  71. #' rt = runif(10, 0.5, 1.5),
  72. #' accuracy = sample(0:1, 10, replace = TRUE)
  73. #' )
  74. #'
  75. #' # Reshape with default duration
  76. #' reshaped1 <- reshape_events(events, event_types = c("condition", "accuracy"))
  77. #'
  78. #' # Reshape with custom duration for each event_type
  79. #' reshaped2 <- reshape_events(events,
  80. #' event_types = c("condition", "accuracy", "rt"),
  81. #' duration = list(condition = 0.5,
  82. #' accuracy = 0.2,
  83. #' rt = function(x) x$rt))
  84. reshape_events <- function(events, event_types, duration = 0.001, modulation = NULL){
  85. if(!(all(c("onset", "run", "subjects") %in% colnames(events)))){
  86. stop("Expected columns: subjects, duration, onset, run")
  87. }
  90. # First check if only 1 numeric entry is present
  91. if(length(duration) == 1 && !is.list(duration)){
  92. duration <- rep(duration, length(event_types))
  93. duration <- lapply(duration, function(x) return(x)) # and make it into a list
  94. } else if(is.list(duration) & any(names(duration) %in% event_types)){
  95. duration_tmp <- replicate(length(event_types), list(0.001)) # Fill in the default spike function
  96. for(i in 1:length(event_types)){
  97. if(any(names(duration) %in% event_types[i])){
  98. duration_tmp[[i]] <- duration[[event_types[i]]]
  99. }
  100. }
  101. duration <- duration_tmp
  102. } else if(length(duration) != length(event_types)){
  103. stop("Length of duration must be 1 or equal to the number of event_types")
  104. }
  105. if(!is.list(duration)){
  106. duration <- lapply(duration, function(x) return(x))
  107. }
  108. out <- list()
  109. for(i in 1:length(event_types)){
  110. fact <- event_types[i]
  111. tmp <- events[,c('subjects', 'run', 'onset', fact)]
  112. if(is.character(tmp[,fact]) || is.factor(tmp[,fact])){
  113. tmp[,fact] <- paste0(fact, "_", tmp[,fact])
  114. colnames(tmp)[4] <- "event_type"
  115. if(is.null(modulation[[fact]])){
  116. tmp$modulation <- 1
  117. } else{
  118. if(is.function(modulation[[fact]])){
  119. tmp$modulation <- modulation[[fact]](events)
  120. } else{
  121. tmp$modulation <- modulation[[fact]]
  122. }
  123. }
  124. } else{
  125. if(!is.null(modulation[[fact]])){
  126. if(is.function(modulation[[fact]])){
  127. tmp[,4] <- modulation[[fact]](events)
  128. } else{
  129. tmp[,4] <- modulation[[fact]]
  130. }
  131. }
  132. colnames(tmp)[4] <- "modulation"
  133. tmp$event_type <- fact
  134. }
  136. if(is.function(duration[[i]])){
  137. tmp$duration <- duration[[i]](events)
  138. } else{
  139. tmp$duration <- duration[[i]]
  140. }
  142. out[[fact]] <- tmp
  143. }
  144. out <- do.call(rbind, out)
  145. rownames(out) <- NULL
  146. out <- out[order(out$subjects, out$run, out$onset),]
  147. return(out)
  148. }
  150. map_MRI <- function(cur_draws){
  151. cur_draws[grepl("sd", rownames(cur_draws)),,] <- exp(cur_draws[grepl("sd", rownames(cur_draws)),,])
  152. cur_draws[grepl("rho", rownames(cur_draws)),,] <- pnorm(cur_draws[grepl("rho", rownames(cur_draws)),,])
  153. return(cur_draws)
  154. }
  158. #' Convolve Events with HRF to Construct Design Matrices
  159. #'
  160. #' This function convolves events with the HRF to construct design matrices for fMRI analysis.
  161. #'
  162. #' @param timeseries A data frame containing fMRI time series data with columns 'subjects', 'run', 'time', and at least one ROI column
  163. #' @param events A data frame containing event information with required columns `subjects`, `run`, `onset`, `duration`, `event_type`, and `modulation`
  164. #' @param factors A named list mapping factor names to event types
  165. #' @param contrasts A named list of contrast matrices for each factor
  166. #' @param covariates A character vector of event types to include as covariates
  167. #' @param add_constant A boolean specifying whether a 1 should be included to the design matrix post convolution
  168. #' @param hrf_model A character string specifying the HRF model to use ('glover', 'spm', 'glover + derivative', or 'spm + derivative')
  169. #' @param cell_coding A character vector of factor names to use cell coding for
  170. #' @param scale A boolean indicating whether to scale the design matrix.
  171. #' @param high_pass Logical indicating whether to apply high-pass filtering.
  172. #' Alternatively, specifying 'add' adds the regressors to the design matrix
  173. #' @param high_pass_model Character indicating which type of high-pass filtering to apply ('cosine', 'poly')
  174. #' @param cut_off A numeric value specifying the cutoff for the high-pass filter
  175. #'
  176. #' @return A list containing the design matrices
  177. #' @export
  178. #' @examples
  179. #' # Generate a simple example timeseries
  180. #' ts <- data.frame(
  181. #' subjects = rep(1, 100),
  182. #' run = rep(1, 100),
  183. #' time = seq(0, 99),
  184. #' ROI1 = rnorm(100)
  185. #' )
  186. #'
  187. #' # Generate example events
  188. #' events <- data.frame(
  189. #' subjects = rep(1, 4),
  190. #' run = rep(1, 4),
  191. #' onset = c(10, 30, 50, 70),
  192. #' duration = rep(0.5, 4),
  193. #' event_type = c("hard", "easy", "hard", "easy"),
  194. #' modulation = c(1, 1, 1, 1)
  195. #' )
  196. #'
  197. #' # Build design matrices
  198. #' design_matrices <- convolve_design_matrix(
  199. #' timeseries = ts,
  200. #' events = events,
  201. #' factors = list(difficulty = c("hard", "easy")),
  202. #' contrasts = list(difficulty = matrix(c(-1, 1)))
  203. #' )
  204. convolve_design_matrix <- function(timeseries, events, factors = NULL, contrasts = NULL,
  205. covariates = NULL, add_constant = TRUE,
  206. hrf_model = 'glover', cell_coding = NULL,
  207. scale = TRUE, high_pass = TRUE,
  208. high_pass_model = "cosine", cut_off = 1e-12) {
  210. if(!(all(c("onset", "run", "subjects", "modulation", "duration") %in% colnames(events)))){
  211. stop("Expected columns in events: subjects, duration, onset, run, modulation, duration")
  212. }
  214. if(!(all(c("run", "subjects", "time") %in% colnames(timeseries)))){
  215. stop("Expected columns in frame_times: run, subjects, time")
  216. }
  218. if(!setequal(unique(timeseries$subjects), unique(events$subjects))){
  219. stop("please make sure timeseries and events have the same subjects")
  220. }
  222. if(!is.null(cell_coding) && !cell_coding %in% names(factors)) stop("Cell coded factors must have same name as factors argument")
  223. # Define double-gamma hyperparameters
  224. if(grepl("glover", hrf_model)){
  225. undershoot <- 12 # When does the negative time point occur
  226. dispersion <- .9 # Width of positive gamma
  227. u_dispersion <- .9 # Width of negative gamma
  228. ratio <- .48 # Relative size of undershoot compared to overshoot
  229. } else{ # Different settings for SPM models
  230. undershoot <- 16
  231. dispersion <- 1
  232. u_dispersion <- 1
  233. ratio <- 1/6
  234. }
  236. if(!is.data.frame(events)) events <- do.call(rbind, events)
  237. subjects <- unique(timeseries$subjects)
  238. # Holders for filtered design matrix
  239. all_dms <- list()
  240. for(subject in subjects){
  241. ev_sub <- events[events$subjects == subject, ]
  242. ts_sub <- timeseries[timeseries$subjects == subject, ]
  243. runs <- unique(ts_sub$run)
  244. # Define subject-wise new design matrix and timeseries
  245. dms_sub <- vector("list", length = length(runs))
  246. for(run in runs){
  247. ev_run <- ev_sub[ev_sub$run == run, ]
  248. ts_run <- ts_sub[ts_sub$run == run, ]
  249. ev_tmp <- data.frame()
  251. for(fact in names(factors)){
  252. idx <- ev_run$event_type %in% factors[[fact]]
  253. if(!any(idx)) next
  254. ev_run$factor[idx] <- fact
  255. tmp <- ev_run[idx, ]
  256. new_tmp <- apply_contrasts(tmp, contrast = contrasts[[fact]],
  257. cell_coding = fact %in% cell_coding,
  258. levels = factors[[fact]])
  259. rownames(new_tmp) <- NULL
  260. new_tmp <- cbind(event_type = new_tmp$regressor, new_tmp)
  261. ev_tmp <- rbind(ev_tmp, new_tmp)
  262. }
  264. for(cov in covariates){
  265. idx <- ev_run$event_type %in% cov
  266. tmp <- ev_run[idx,c('event_type', 'subjects', 'run', 'onset', 'duration', 'modulation')]
  267. tmp$regressor <- cov
  268. ev_tmp <- rbind(ev_tmp, tmp)
  269. }
  271. ev_tmp <- ev_tmp[order(ev_tmp$onset), ]
  272. if((run == runs[1]) & (subject == subjects[1])){
  273. round_ev <- ev_tmp
  274. round_ev$duration <- round(scale(round_ev$duration))
  275. round_ev$modulation <- round(scale(round_ev$modulation))
  276. unq_idx <- !duplicated(round_ev[, !colnames(round_ev) %in% c("onset", "subjects")])
  277. print(ev_tmp[unq_idx,])
  278. }
  279. # Event_type was only included for printing
  280. ev_tmp <- ev_tmp[,colnames(ev_tmp) != "event_type"]
  281. dm <- construct_design_matrix(ts_run$time,
  282. events = ev_tmp,
  283. has_derivative = grepl("derivative", hrf_model),
  284. time_length = 32, # total hrf duration
  285. min_onset = -24, # Grid computation start time for oversampling
  286. oversampling = 50, # How many timepoints per tr
  287. onset = 0, # accounts for shifts in HRF
  288. delay = 6, # Time to peak of initial bump
  289. undershoot = undershoot,
  290. dispersion = dispersion,
  291. u_dispersion = u_dispersion,
  292. ratio = ratio,
  293. add_intercept = FALSE)
  294. if((run == runs[1]) & (subject == subjects[1]) & isTRUE(high_pass)){
  295. message("Filtering out high_pass noise, make sure you also use high_pass_filter(<timeseries>)")
  296. }
  297. if(!isFALSE(high_pass)){
  298. dm <- high_pass_filter(dm, high_pass_model, frame_times = ts_run$time, add=(high_pass == "add"))
  299. }
  300. if(add_constant) dm$constant <- 1
  301. dms_sub[[as.character(run)]] <- dm
  302. }
  303. dms_sub <- Filter(Negate(is.null), dms_sub)
  304. dm_cols <- unique(unlist(lapply(dms_sub, colnames), use.names = FALSE))
  305. dms_sub <- lapply(dms_sub, function(dm) {
  306. for(col in setdiff(dm_cols, colnames(dm))) dm[[col]] <- 0
  307. dm[, dm_cols, drop = FALSE]
  308. })
  309. dms_sub <- do.call(rbind, dms_sub)
  310. dms_sub[abs(dms_sub) < cut_off] <- 0
  311. rownames(dms_sub) <- NULL
  312. all_dms[[as.character(subject)]] <- dms_sub
  313. }
  314. dm_cols <- unique(unlist(lapply(all_dms, colnames), use.names = FALSE))
  315. all_dms <- lapply(all_dms, function(dm) {
  316. for(col in setdiff(dm_cols, colnames(dm))) dm[[col]] <- 0
  317. dm[, dm_cols, drop = FALSE]
  318. })
  319. if(scale){
  320. full_dm <- do.call(rbind, all_dms)
  321. maxs <- apply(full_dm, 2, max)
  322. all_dms <- lapply(all_dms, function(x){
  323. for(i in 1:ncol(x)){
  324. x[,i] <- x[,i]/maxs[i]
  325. }
  326. return(x)
  327. })
  328. }
  329. return(all_dms)
  330. }
  332. #' Split fMRI Timeseries Data by ROI Columns
  333. #'
  334. #' This function splits a timeseries data frame containing multiple ROI columns into a list
  335. #' of data frames, where each data frame contains the common columns (subjects, run, time)
  336. #' and one ROI column.
  337. #'
  338. #' @param timeseries A data frame containing fMRI timeseries data with required columns
  339. #' 'subjects', 'run', and 'time', plus one or more ROI columns.
  340. #' @param columns A character vector specifying which columns to split by. If NULL (default),
  341. #' all columns except 'subjects', 'run', and 'time' will be used.
  342. #'
  343. #' @return A named list of data frames, where each data frame contains the common columns
  344. #' (subjects, run, time) and one ROI column. The names of the list elements correspond
  345. #' to the ROI column names.
  346. #'
  347. #' @export
  348. #'
  349. #' @examples
  350. #' # Create a simple example timeseries with multiple ROIs
  351. #' set.seed(123)
  352. #' n_frames <- 100
  353. #'
  354. #' # Create a data frame with multiple ROIs
  355. #' timeseries <- data.frame(
  356. #' subjects = rep(1, n_frames),
  357. #' run = rep(1, n_frames),
  358. #' time = seq(0, n_frames-1),
  359. #' ROI1 = rnorm(n_frames),
  360. #' ROI2 = rnorm(n_frames),
  361. #' ROI3 = rnorm(n_frames)
  362. #' )
  363. #'
  364. #' # Split the timeseries by all ROI columns
  365. #' split_data <- split_timeseries(timeseries)
  366. split_timeseries <- function(timeseries, columns = NULL){
  367. if(!(all(c("run", "subjects", "time") %in% colnames(timeseries)))){
  368. stop("Expected columns in timeseries: run, subjects, time")
  369. }
  370. if(is.null(columns)) columns <- colnames(timeseries)[!colnames(timeseries) %in% c('run', 'subjects', 'time')]
  371. out <- list()
  372. for(col in columns){
  373. if(!col %in% colnames(timeseries)) stop("Please ensure selected columns are in timeseries")
  374. out[[col]] <- timeseries[,c('subjects', 'run', 'time', col)]
  375. }
  376. return(out)
  377. }
  380. #' Apply High-Pass Filtering to fMRI Data
  381. #'
  382. #' This function applies high-pass filtering to fMRI data to remove low-frequency noise
  383. #' and drift. It supports two filtering methods: cosine basis functions and polynomial
  384. #' regressors.
  385. #'
  386. #' @param X A data frame or matrix containing the data to be filtered. If it contains
  387. #' columns 'subjects' and 'run', the function will apply filtering separately for
  388. #' each subject-run combination.
  389. #' @param high_pass_model A character string specifying the high-pass filtering method.
  390. #' Options are 'cosine' (default) or 'poly' for polynomial regressors.
  391. #' @param frame_times A numeric vector of time points for each frame. If NULL, the
  392. #' function will attempt to extract this from a 'time' column in X.
  393. #' @param ... Additional arguments passed to the function.
  394. #'
  395. #' @return A data frame or matrix with the same structure as X, but with high-frequency
  396. #' components removed from the data columns.
  397. #'
  398. #' @export
  399. #'
  400. #' @examples
  401. #' # Create a simple example data frame with drift
  402. #' set.seed(123)
  403. #' n_frames <- 100
  404. #' time <- seq(0, 99)
  405. #'
  406. #' # Create a signal with low-frequency drift
  407. #' drift <- 0.1 * time
  408. #' signal <- sin(2 * pi * 0.1 * time) + drift
  409. #' noise <- rnorm(n_frames, 0, 0.5)
  410. #' data <- signal + noise
  411. #'
  412. #' # Create a data frame
  413. #' df <- data.frame(
  414. #' time = time,
  415. #' signal = data
  416. #' )
  417. #'
  418. #' # Apply high-pass filtering using cosine basis functions
  419. #' filtered_df <- high_pass_filter(df, high_pass_model = "cosine")
  420. high_pass_filter <- function(X, high_pass_model = 'cosine', frame_times = NULL, ...){
  421. if(is.null(frame_times)){
  422. if(!'time' %in% colnames(X)){
  423. stop("no column named 'time' for frame_times present, please separately provide")
  424. } else{
  425. frame_times <- X[,'time']
  426. }
  427. }
  428. if('subjects' %in% colnames(X) && is.null(list(...)$recursive)){
  429. out <- list()
  430. k <- 0
  431. message("Make sure you also high_pass_filter your events (set high_pass = TRUE in convolve_design_matrix)")
  432. for(sub in unique(X[,'subjects'])){
  433. tmp <- X[X[,'subjects'] == sub,]
  434. for(run in unique(tmp[,'run'])){
  435. k <- k + 1
  436. tmp_run <- tmp[tmp[,'run'] == run,]
  437. tmp_run <- high_pass_filter(tmp_run, recursive = TRUE)
  438. out[[k]] <- tmp_run
  439. }
  440. }
  441. return(do.call(rbind, out))
  442. }
  443. if(high_pass_model == "cosine"){
  444. nuisance <- cosine_drift(frame_times)
  445. } else if(high_pass_model == "poly"){
  446. nuisance <- poly_drift(frame_times)
  447. } else{
  448. stop("Only poly and cosine are supported as high_pass_model")
  449. }
  450. if(!is.null(list(...)$add)){
  451. if(list(...)$add){
  452. return(cbind(X, nuisance))
  453. }
  454. }
  456. gets_filter <- !colnames(X) %in% c('subjects', 'run', 'time')
  457. for(i in 1:ncol(X)){
  458. if(gets_filter[i]){
  459. fit <- lm(X[,i] ~ nuisance - 1)
  460. X[,i] <- residuals(fit)
  461. }
  462. }
  463. return(X)
  464. }
  466. poly_drift <- function(frame_times, order = 3) {
  467. # Ensure that 'order' is an integer
  468. order <- as.integer(order)
  469. n <- length(frame_times)
  470. # Compute maximum of frame_times (will be used to scale the time vector)
  471. tmax <- max(frame_times)
  473. # Create a matrix where column k corresponds to (frame_times/tmax)^k,
  474. # for k = 0, 1, ..., order; note that the 0th power yields a constant.
  475. pol <- sapply(0:order, function(k) (frame_times / tmax) ^ k)
  476. # 'pol' is now an n x (order+1) matrix
  478. # Orthogonalize the columns using QR decomposition.
  479. # The function qr.Q returns an orthonormal basis for the columns.
  480. pol_orth <- qr.Q(qr(pol))
  482. # Drop the constant
  483. result <- pol_orth[, -1, drop = FALSE]
  484. colnames(result) <- paste0("pol_", 1:ncol(result))
  485. return(result)
  486. }
  489. cosine_drift <- function(frame_times, high_pass = .01) {
  490. n_frames <- length(frame_times)
  491. n_times <- 0:(n_frames - 1)
  492. # Compute the time interval (dt) between frames.
  493. dt <- (frame_times[n_frames] - frame_times[1]) / (n_frames - 1)
  495. # Check if the product high_pass * dt is too high, issuing a warning if so.
  496. if (high_pass * dt >= 0.5) {
  497. warning(sprintf("High-pass filter will span all accessible frequencies and saturate the design matrix. You may want to reduce the high_pass value. The provided value is %.3f Hz", high_pass))
  498. }
  500. # Determine the number of cosine basis functions (excluding the constant).
  501. # This follows: order = min(n_frames - 1, floor(2 * n_frames * high_pass * dt))
  502. order <- min(n_frames - 1, floor(2 * n_frames * high_pass * dt))
  504. # Create a matrix to hold the cosine basis functions and a constant column.
  505. # The result will have (order + 1) columns.
  506. cosine_drift <- matrix(0, nrow = n_frames, ncol = order + 1)
  507. normalizer <- sqrt(2.0 / n_frames)
  509. # Fill the first 'order' columns with the cosine functions.
  510. # For each k = 1, 2, ..., order we compute:
  511. # normalizer * cos( (pi/n_frames) * (n_times + 0.5) * k )
  512. for (k in 1:order) {
  513. cosine_drift[, k] <- normalizer * cos((pi / n_frames) * (n_times + 0.5) * k)
  514. }
  515. cosine_drift <- cosine_drift[,-ncol(cosine_drift)]
  517. # Set the last column to a constant of 1.
  518. colnames(cosine_drift) <- paste0("drift_", 1:ncol(cosine_drift))
  520. return(cosine_drift)
  521. }
  524. #' Create fMRI Design for EMC2 Sampling
  525. #'
  526. #' This function takes the output from convolve_design_matrix and transforms it into a design
  527. #' suitable for sampling with EMC2. It properly configures parameter types, bounds, and transformations
  528. #' for the specified model.
  529. #'
  530. #' @param design_matrix A list of design matrices, the output from convolve_design_matrix
  531. #' @param model A function that returns a model specification, options are MRI or MRI_AR1
  532. #' @param ... Additional arguments passed to the model
  533. #'
  534. #' @return An object of class 'emc.design' suitable for EMC2 sampling
  535. #' @export
  536. #'
  537. #' @examples
  538. #' # Generate a simple example timeseries
  539. #' ts <- data.frame(
  540. #' subjects = rep(1, 100),
  541. #' run = rep(1, 100),
  542. #' time = cumsum(rep(1.38, 100)),
  543. #' ROI1 = rnorm(100)
  544. #' )
  545. #'
  546. #' # Generate example events
  547. #' events <- data.frame(
  548. #' subjects = rep(1, 4),
  549. #' run = rep(1, 4),
  550. #' onset = c(10, 30, 50, 70),
  551. #' duration = rep(0.5, 4),
  552. #' event_type = c("A", "B", "A", "B"),
  553. #' modulation = c(1, 1, 1, 1)
  554. #' )
  555. #'
  556. #' # Create convolved design matrix
  557. #' design_matrix <- convolve_design_matrix(
  558. #' timeseries = ts,
  559. #' events = events,
  560. #' factors = list(condition = c("A", "B")),
  561. #' hrf_model = "glover"
  562. #' )
  563. #'
  564. #' # Create fMRI design for EMC2
  565. #' fmri_design <- design_fmri(design_matrix, model = MRI_AR1)
  567. design_fmri <- function(design_matrix,
  568. model = MRI_AR1, ...) {
  569. dots <- list(...)
  570. betas <- colnames(design_matrix[[1]])
  571. subjects <- names(design_matrix)
  572. model_list <- model()
  573. # Fill in new p_types
  574. p_not_beta <- model_list$p_types[names(model_list$p_types) != "beta"]
  575. model_list$p_types <- c(setNames(rep(model_list$p_types["beta"], length(betas)), betas), p_not_beta)
  576. # Fill in new bound
  577. new_mm <- do.call(cbind, rep(list(model_list$bound$minmax[,"beta"]), length(betas)))
  578. colnames(new_mm) <- betas
  579. model_list$bound$minmax <- cbind(model_list$bound$minmax, new_mm)
  580. model_list$bound$minmax <- model_list$bound$minmax[,colnames(model_list$bound$minmax) != "beta"]
  582. # Fill in new transforms
  583. new_t <- setNames(rep(model_list$transform$func["beta"], length(betas)), betas)
  584. model_list$transform$func <- c(model_list$transform$func, new_t)
  585. model_list$transform$func <- model_list$transform$func[names(model_list$transform$func) != "beta"]
  586. # Make pre transforms
  587. par_names <- c(betas, names(p_not_beta))
  588. model_list$pre_transform$func <- setNames(rep("identity", length(par_names)), par_names)
  589. model <- function() {return(model_list)}
  590. n_pars <- length(par_names)
  592. # Fill up final results
  593. model_list$transform <- fill_transform(dots$transform,model)
  594. model_list$bound <- fill_bound(dots$bound,model)
  595. model_list$pre_transform <- fill_transform(dots$pre_transform, model = model, p_vector = model_list$p_types, is_pre = TRUE)
  596. model <- function(){return(model_list)}
  600. Flist <- vector("list", n_pars)
  601. for(i in 1:n_pars){
  602. Flist[[i]] <- as.formula(paste0(par_names[i], "~1"))
  603. }
  605. design <- list(Flist = Flist, model = model, Ffactors = list(subjects = subjects))
  606. attr(design, "design_matrix") <- lapply(design_matrix, FUN=function(x) {
  607. y <- x[,colnames(x) != 'subjects']
  608. DM_tmp <- data.matrix(y)
  609. rownames(DM_tmp) <- NULL
  610. return(DM_tmp)
  611. })
  612. par_names <- setNames(numeric(length(par_names)), par_names)
  613. attr(design, "p_vector") <- par_names
  614. class(design) <- 'emc.design'
  615. return(design)
  616. }
  619. #' GLM model for fMRI data
  620. #'
  621. #' Creates a model specification for fMRI data using a normal distribution.
  622. #' This model assumes that the observed BOLD signal follows a normal distribution
  623. #' with a mean determined by the design matrix and betas, and a standard deviation
  624. #' parameter for noise.
  625. #'
  626. #' @return A list containing model specification
  627. #'
  628. #' @details
  629. #' The model uses a normal distribution to model fMRI BOLD signals.
  630. #' Beta parameters represent the effect sizes for different conditions,
  631. #' and the sd parameter represents the standard deviation of the noise.
  632. #'
  633. #' The log-likelihood function centers the predicted values by subtracting
  634. #' the mean, which helps with model identifiability.
  635. #'
  636. #' @export
  637. #' @examples
  638. #' # Create a normal MRI model specification
  639. #' model_spec <- MRI()
  640. #'
  641. #' # Access model parameters
  642. #' model_spec$p_types
  643. MRI <- function(){
  644. return(
  645. list(
  646. type="MRI",
  647. c_name = "MRI",
  648. p_types=c("beta" = 0, "sd" = log(1)),
  649. transform=list(func=c(beta = "identity", sd = "exp")),
  650. bound=list(minmax=cbind(beta=c(-Inf,Inf),sd=c(0.001,Inf))),
  651. Ttransform = function(pars, dadm) return(pars),
  652. rfun=function(pars){
  653. # - Each row corresponds to an observation
  654. # - All columns except the last are betas (already multiplied by the design matrix)
  655. # - The last column is sigma (the noise standard deviation)
  656. # Extract sigma and betas
  657. sigma <- pars[, ncol(pars)]
  658. betas <- pars[, -ncol(pars)]
  659. # Compute the predicted mean for each observation as the sum of its betas
  660. y_hat <- rowSums(betas)
  661. # Generate simulated data: for each observation, add noise drawn from a normal distribution
  662. # with mean 0 and standard deviation sigma.
  663. y_sim <- y_hat + rnorm(n = length(y_hat), mean = 0, sd = sigma)
  664. return(y_sim)
  665. },
  666. log_likelihood=function(pars, dadm, model, min_ll=log(1e-10)){
  667. # Here pars already contains the mapped design contributions.
  668. y <- as.matrix(dadm[,!colnames(dadm) %in% c("subjects", 'run', 'time', "trials")])
  669. # grab the right parameters
  670. sigma <- pars[,ncol(pars)]
  671. betas <- pars[,-ncol(pars)]
  672. y_hat <- rowSums(betas)
  673. ll <- sum(pmax(dnorm(as.matrix(y), mean = y_hat, sd = sigma, log = T), min_ll))
  674. return(ll)
  675. }
  676. )
  677. )
  678. }
  680. #' Create an AR(1) GLM model for fMRI data
  681. #'
  682. #' This function creates a model specification for MRI data with an AR(1) error structure.
  683. #' The model includes beta parameters for the design matrix, a rho parameter for the
  684. #' autocorrelation, and a standard deviation parameter for the noise.
  685. #'
  686. #' The AR(1) model accounts for temporal autocorrelation in the data, where each timepoint
  687. #' is correlated with the previous timepoint according to the rho parameter.
  688. #'
  689. #' @return A list containing the model specifications
  690. #'
  691. #' @export
  692. #' @examples
  693. #' # Create an AR(1) GLM model for fMRI data
  694. #' model_spec <- MRI_AR1()
  695. #'
  696. #' # Access model parameters
  697. #' model_spec$p_types
  699. MRI_AR1 <- function(){
  700. return(
  701. list(
  702. type="MRI_AR1",
  703. c_name = "MRI_AR1",
  704. p_types=c("beta" = 0, "rho" = pnorm(0.001), "sd" = log(1)),
  705. transform=list(func=c(beta = "identity", rho = "pnorm", sd = "exp")),
  706. bound=list(minmax=cbind(beta=c(-Inf,Inf),sd=c(0.001,Inf), rho = c(0.0001, 1)),
  707. exception=c(rho=0)),
  708. Ttransform = function(pars, dadm) return(pars),
  709. rfun=function(pars){
  710. n <- nrow(pars)
  711. m <- ncol(pars)
  713. # - betas: columns 1 to (m-2)
  714. # - rho: column (m-1)
  715. # - sigma: column m (stationary standard deviation)
  716. betas <- pars[,1:(m-2), drop = FALSE]
  717. rho <- pars[,m-1]
  718. sigma <- pars[,m]
  720. # Compute the linear predictor (sum of beta contributions) and center it.
  721. y_hat <- rowSums(betas)
  722. # y_hat <- y_hat - mean(y_hat)
  724. # Allocate a vector for simulated data
  725. y_sim <- numeric(n)
  727. # Simulate the first observation
  728. y_sim[1] <- y_hat[1] + rnorm(1, mean = 0, sd = sigma[1])
  730. # Loop through time for the remaining observations
  731. for (t in 2:n) {
  732. # The conditional mean for observation t
  733. cond_mean <- y_hat[t] + rho[t] * (y_sim[t - 1] - y_hat[t - 1])
  734. # The conditional standard deviation for observation t
  735. cond_sd <- sigma[t] * sqrt(1 - rho[t]^2)
  736. # Simulate
  737. y_sim[t] <- cond_mean + rnorm(1, mean = 0, sd = cond_sd)
  738. }
  739. return(y_sim)
  740. },
  741. log_likelihood = function(pars, dadm, model, min_ll = log(1e-10)) {
  742. # Here pars already contains the mapped design contributions.
  744. # Extract observed data (as a vector)
  745. y <- as.vector(as.matrix(dadm[, !colnames(dadm) %in% c("subjects", "run", "time", "trials")]))
  746. n <- length(y)
  748. m <- ncol(pars)
  749. betas <- pars[, 1:(m - 2), drop = FALSE]
  750. rho <- pars[, m - 1]
  751. sigma <- pars[, m]
  753. y_hat <- rowSums(betas)
  754. # y_hat <- y_hat - mean(y_hat)
  756. # Log-likelihood for the first observation
  757. ll <- numeric(n)
  758. ll[1] <- dnorm(y[1], mean = y_hat[1], sd = sigma[1], log = TRUE)
  760. # For observations t = 2:n, compute conditional means
  761. cond_mean <- y_hat[-1] + rho[-1] * (y[-n] - y_hat[-n])
  762. cond_sd <- sigma[-1] * sqrt(1 - rho[-1]^2)
  763. ll[-1] <- dnorm(y[-1], mean = cond_mean, sd = cond_sd, log = TRUE)
  765. ll <- pmax(ll, min_ll)
  766. return(sum(ll))
  767. }
  768. )
  769. )
  770. }
  773. # Plotting functions ------------------------------------------------------
  775. #' Plot fMRI Design Matrix
  776. #'
  777. #' This function creates a visualization of an fMRI design matrix, showing the temporal
  778. #' evolution of regressors over time. It can handle various input formats and provides
  779. #' options to customize the visualization.
  780. #'
  781. #' @param design_matrix A design matrix for fMRI analysis. Can be a data frame, matrix,
  782. #' list of matrices, or an object of class 'emc.design'.
  783. #' @param TRs The number of time points (TRs) to plot. Default is 100.
  784. #' @param events A character vector specifying which regressors to plot. If NULL,
  785. #' all non-nuisance regressors will be plotted.
  786. #' @param remove_nuisance Logical indicating whether to remove nuisance regressors
  787. #' (drift terms, polynomial terms, derivatives) from the plot. Default is TRUE.
  788. #' @param subject The subject number to plot. Only applies for list of design matrices. Default is 1.
  789. #' @param legend_pos Position of the legend. Default is "bottomleft".
  790. #' @param ... Additional graphical parameters passed to matplot and legend.
  791. #'
  792. #' @return A plot showing the design matrix regressors over time.
  793. #'
  794. #' @export
  795. #'
  796. #' @examples
  797. #' # Example time series
  798. #' ts <- data.frame(
  799. #' subjects = rep(1, 100),
  800. #' run = rep(1, 100),
  801. #' time = seq(0, 99),
  802. #' ROI = rnorm(100)
  803. #' )
  804. #' # Create a simple events data frame
  805. #' events <- data.frame(
  806. #' subjects = rep(1, 10),
  807. #' run = rep(1, 10),
  808. #' onset = seq(0, 90, by = 10),
  809. #' condition = rep(c("A", "B"), 5),
  810. #' rt = runif(10, 0.5, 1.5),
  811. #' accuracy = sample(0:1, 10, replace = TRUE)
  812. #' )
  813. #' # Reshape with custom duration for each event_type
  814. #' reshaped <- reshape_events(events,
  815. #' event_types = c("condition", "accuracy", "rt"),
  816. #' duration = list(condition = 0.5,
  817. #' accuracy = 0.2,
  818. #' rt = function(x) x$rt))
  819. #' design_matrices <- convolve_design_matrix(
  820. #' timeseries = ts,
  821. #' events = reshaped,
  822. #' covariates = c('accuracy', 'rt'),
  823. #' factors = list(cond = c("condition_A", "condition_B")),
  824. #' contrasts = list(cond = matrix(c(-1, 1))))
  825. #'
  826. #' # Plot the design matrix
  827. #' plot_design_fmri(design_matrices)
  828. plot_design_fmri <- function(design_matrix, TRs = 100, events = NULL, remove_nuisance = TRUE, subject = 1,
  829. legend_pos = "bottomleft", ...){
  830. if(is(design_matrix,"emc.design")){
  831. design_matrix <- attr(design_matrix, "design_matrix")
  832. }
  833. if(!is.data.frame(design_matrix) && !is.matrix(design_matrix)){
  834. if(is.list(design_matrix)) design_matrix <- design_matrix[[subject]]
  835. }
  836. enames <- colnames(design_matrix)
  837. if(remove_nuisance & is.null(events)){
  838. is_nuisance <- grepl("drift", enames) | grepl("poly", enames) | grepl("derivative", enames) | apply(design_matrix, 2, sd) == 0
  839. design_matrix <- design_matrix[,!is_nuisance, drop = F]
  840. }
  841. enames <- colnames(design_matrix)
  842. if(is.null(events)){
  843. events <- enames
  844. } else{
  845. if(any(!events %in% enames)){
  846. stop("events not in colnames design matrix")
  847. }
  848. }
  849. distinct_colors <- c(
  850. "#E6194B", "#3CB44B", "#0082C8", "#F58231", "#911EB4",
  851. "#46F0F0", "#F032E6", "#D2F53C", "#FABEBE", "#008080",
  852. "#E6BEFF", "#AA6E28", "#FFFAC8", "#800000", "#FFD8B1"
  853. )
  854. dots <- add_defaults(list(...), col = distinct_colors, lwd = 2, lty = 1, main = NULL, xlab = "TRs", ylab = "Amplitude")
  856. TRs <- min(c(TRs, nrow(design_matrix)))
  857. design_matrix <- design_matrix[1:TRs, events, drop = F]
  858. do.call(matplot, c(list(design_matrix, type = "l"), fix_dots_plot(dots)))
  859. do.call(legend, c(list(legend_pos, legend = events, bty = "n"), fix_dots(dots, legend)))
  860. }
  864. # -------------------------------------------------------------------------
  865. # Filter events by event_type and bin/categorize the modulation
  866. # -------------------------------------------------------------------------
  867. prepare_event_groups <- function(events, event_type, n_bins = 4) {
  868. ev_sub <- events[events$event_type == event_type, ]
  869. if (nrow(ev_sub) == 0) {
  870. return(ev_sub) # empty
  871. }
  872. if (!"modulation" %in% names(ev_sub)) {
  873. stop("The 'events' data frame must have a 'modulation' column.")
  874. }
  876. # Round to reduce floating precision issues
  877. ev_sub$modulation <- round(ev_sub$modulation, 6)
  878. mod_vals <- unique(ev_sub$modulation)
  879. n_unique <- length(mod_vals)
  881. if (n_unique > 6) {
  882. # Calculate quartile breakpoints for the 'modulation' data
  883. quartile_breaks <- quantile(ev_sub$modulation, probs = seq(0, 1, length.out = n_bins + 1), na.rm = TRUE)
  885. # Bin the data into quartiles using these breakpoints
  886. ev_sub$mod_group <- cut(ev_sub$modulation, breaks = quartile_breaks, include.lowest = TRUE, labels = 1:n_bins)
  887. attr(ev_sub, "binned") <- TRUE
  889. } else {
  890. # treat as categorical
  891. ev_sub$mod_group <- factor(ev_sub$modulation, levels = sort(mod_vals))
  892. attr(ev_sub, "binned") <- FALSE
  893. }
  894. return(ev_sub)
  895. }
  897. # -----------------------------------------------------------------------------
  898. # Fast FIR utilities -----------------------------------------------------------
  899. # -----------------------------------------------------------------------------
  901. # 1) Build an FIR design matrix for one set of onsets --------------------------
  902. # frame_times : numeric vector of acquisition times (s) **monotone**
  903. # onsets : numeric vector of event onsets (s) (single modulation group)
  904. # pre, post : seconds before/after onset to model (positive numbers)
  905. #
  906. # Returns a dense matrix (length(frame_times) x K) where
  907. # K = floor((pre+post)/dt)+1 (dt = median diff(frame_times))
  908. # Column k corresponds to lag t = -pre + (k-1)*dt.
  909. # ---------------------------------------------------------------------------
  910. # Build an FIR design with an optional externally‑specified TR --------------
  911. # ---------------------------------------------------------------------------
  912. ###############################################################################
  913. ## Fast FIR utilities + plot_fmri (robust version, 3 May 2025) ##
  914. ###############################################################################
  916. # ─────────────────────────────────────────────────────────────────────────────
  918. # ─────────────────────────────────────────────────────────────────────────────
  919. # 1) FIR design – *global lag grid* ensured by `fixed_dt` argument
  920. # ─────────────────────────────────────────────────────────────────────────────
  921. build_fir_design <- function(frame_times, onsets, durations,
  922. pre, post, fixed_dt, weights = NULL) {
  924. if (is.null(weights)) weights <- rep(1, length(onsets))
  925. if (is.null(durations)) durations <- rep(0, length(onsets)) # impulse
  926. stopifnot(length(durations) == length(onsets))
  928. dt <- fixed_dt
  929. lags <- seq(-pre, post, by = dt)
  930. K <- length(lags)
  931. nTR <- length(frame_times)
  933. X <- matrix(0, nrow = nTR, ncol = K,
  934. dimnames = list(NULL, sprintf("lag_%0.3f", lags)))
  936. ## pre‑compute frame centres for fast look‑up
  937. fc <- frame_times
  939. for (j in seq_along(onsets)) {
  940. onset <- onsets[j]
  941. dur <- durations[j]
  942. w <- weights[j]
  944. ## frames whose centre is within the event window
  945. in_evt <- which(fc >= onset & fc < onset + dur)
  946. if (!length(in_evt)) { # shorter than first half‑TR → impulse
  947. in_evt <- which.min(abs(fc - onset))
  948. }
  950. ## update FIR cols for every lag
  951. for (k in seq_len(K)) {
  952. rows <- in_evt + round(lags[k] / dt)
  953. rows <- rows[rows >= 1 & rows <= nTR]
  954. if(length(rows))
  955. X[rows, k] <- X[rows, k] + w
  956. }
  957. }
  958. attr(X, "lag_seconds") <- lags
  959. X
  960. }
  963. # ─────────────────────────────────────────────────────────────────────────────
  964. # 2) AR(1) whitening
  965. # ─────────────────────────────────────────────────────────────────────────────
  966. estimate_rho <- function(y, clip = TRUE) {
  967. if (length(y) < 2L) return(0)
  969. y <- y - mean(y, na.rm = TRUE) # 1 de‑mean
  970. rho <- sum(y[-1] * y[-length(y)]) /
  971. (sum(y[-length(y)]^2) + 1e-12)
  973. if (clip) rho <- max(-0.99, min(0.99, rho)) # 2 safe bounds
  974. rho
  975. }
  977. whiten_series <- function(mat_or_vec, rho) {
  978. if (abs(rho) < 1e-6) return(mat_or_vec)
  980. if (is.vector(mat_or_vec))
  981. return(c(mat_or_vec[1],
  982. mat_or_vec[-1] - rho * mat_or_vec[-length(mat_or_vec)]))
  984. rbind(mat_or_vec[1, , drop = FALSE],
  985. mat_or_vec[-1, , drop = FALSE] -
  986. rho * mat_or_vec[-nrow(mat_or_vec), , drop = FALSE])
  987. }
  990. # ─────────────────────────────────────────────────────────────────────────────
  991. # 4) Quick FIR GLM fit (OLS + AR(1) whitening, baseline correction)
  992. # ─────────────────────────────────────────────────────────────────────────────
  993. fit_fir_glm <- function(y, X, lags_sec, target) {
  994. rho <- estimate_rho(y)
  995. yw <- whiten_series(y, rho)
  996. Xw <- whiten_series(X, rho)
  998. beta <- as.vector(solve(crossprod(Xw), crossprod(Xw, yw)))
  999. beta <- beta[target]
  1000. ## Baseline correction (pre‑stimulus lags)
  1001. idx_pre <- lags_sec < 0
  1002. if(any(idx_pre, na.rm = TRUE))
  1003. beta <- beta - mean(beta[idx_pre], na.rm = TRUE)
  1004. beta
  1005. }
  1007. # ─────────────────────────────────────────────────────────────────────────────
  1008. # 5) Main extractor – returns long data.frame for plotting
  1009. # ─────────────────────────────────────────────────────────────────────────────
  1010. # ─────────────────────────────────────────────────────────────────────────────
  1011. # 5) Main extractor – FIR for *all* events, returns long data‑frame
  1012. # ─────────────────────────────────────────────────────────────────────────────
  1013. get_fir_lines <- function(timeseries, events, event_type,
  1014. pre = 2, post = 18, n_bins = 4,
  1015. high_pass = TRUE, high_pass_model = "cosine") {
  1017. ## basic checks
  1018. stopifnot(all(c("run", "subjects", "time") %in% names(timeseries)))
  1019. ROI_col <- setdiff(names(timeseries),
  1020. c("run", "subjects", "time", "postn"))
  1021. if(length(ROI_col) != 1L) stop("Exactly one ROI column required")
  1023. ## global lag grid
  1024. global_dt <- median(diff(sort(unique(timeseries$time))))
  1025. global_lags <- seq(-pre, post, by = global_dt)
  1026. K <- length(global_lags)
  1028. ## TARGET events (binned)
  1029. ev_sub <- prepare_event_groups(events, event_type, n_bins)
  1030. if(nrow(ev_sub) == 0)
  1031. stop("No events found for event_type = ", event_type)
  1032. groups <- sort(unique(ev_sub$mod_group))
  1034. ## NUISANCE events (everything else)
  1035. ev_nuis <- events[events$event_type != event_type, ]
  1036. nuis_types <- sort(unique(ev_nuis$event_type)) # might be empty
  1038. betas <- setNames(vector("list", length(groups)), groups)
  1040. ts_split <- split(timeseries, list(timeseries$subjects,
  1041. timeseries$run), drop = TRUE)
  1043. for(chunk in ts_split) {
  1044. ft <- chunk$time
  1045. y <- chunk[[ROI_col]]
  1046. sid <- chunk$subjects[1]
  1047. rid <- chunk$run[1]
  1049. ## 1) build FIR blocks for *each nuisance type* in this run
  1050. X_nuis <- NULL
  1051. if(length(nuis_types)) {
  1052. for(et in nuis_types) {
  1053. ev_tmp <- ev_nuis[ev_nuis$event_type == et &
  1054. ev_nuis$subjects == sid &
  1055. ev_nuis$run == rid, ]
  1056. if(!nrow(ev_tmp)) next
  1058. if(sd(ev_tmp$modulation) == 0 & sd(ev_tmp$duration) == 0) next
  1061. ## centre the modulator so ‘main effect’ and modulation are orthogonal
  1062. if (sd(ev_tmp$modulation) > 0) {
  1063. ev_tmp$modulation <- ev_tmp$modulation - mean(ev_tmp$modulation)
  1064. } else{
  1065. ev_tmp$modulation <- ev_tmp$duration - mean(ev_tmp$duration)
  1066. }
  1067. Xet <- build_fir_design(ft,
  1068. ev_tmp$onset,
  1069. ev_tmp$duration,
  1070. pre, post,
  1071. fixed_dt = global_dt,
  1072. weights = ev_tmp$modulation)
  1073. colnames(Xet) <- paste0(et, "_", colnames(Xet))
  1074. X_nuis <- if(is.null(X_nuis)) Xet else cbind(X_nuis, Xet)
  1075. }
  1076. }
  1078. ## 2) loop over modulation groups of the TARGET event
  1079. for(g in groups) {
  1080. ev_g <- ev_sub[ev_sub$mod_group == g &
  1081. ev_sub$subjects == sid &
  1082. ev_sub$run == rid, ]
  1083. if(!nrow(ev_g)) next
  1085. ## NOTE: weights = NULL → un‑scaled FIR for the curve we plot
  1086. X_tar <- build_fir_design(ft,
  1087. ev_g$onset,
  1088. ev_g$duration, # durations for the target
  1089. pre, post,
  1090. fixed_dt = global_dt,
  1091. weights = NULL) # no amplitude scaling
  1093. X_full <- if(is.null(X_nuis)) X_tar else cbind(X_tar, X_nuis)
  1095. if(!isFALSE(high_pass)){
  1096. X_full <- high_pass_filter(X_full, high_pass_model, frame_times = ft, add=(high_pass == "add"))
  1097. }
  1098. b <- fit_fir_glm(y, X_full, global_lags, target = seq_len(K))
  1100. if(all(is.na(b))) next
  1101. betas[[g]] <- rbind(betas[[g]],
  1102. cbind(t(b), nrow(ev_g)))
  1103. }
  1104. }
  1106. ## 3) weight‑averaged β → long data.frame
  1107. out <- lapply(names(betas), function(g) {
  1108. mat <- betas[[g]]; if(is.null(mat)) return(NULL)
  1109. B <- mat[, 1:K, drop = FALSE]; w <- mat[, K+1L]
  1110. if(!is.matrix(B)) { B <- matrix(B, nrow = 1); w <- w[1] }
  1112. beta_avg <- colSums(B * w) / sum(w)
  1114. df <- data.frame(
  1115. mod_group = rep(g, K),
  1116. time = global_lags,
  1117. avg_signal = beta_avg
  1118. )
  1119. if(attr(ev_sub, "binned")) {
  1120. ev_all <- ev_sub[ev_sub$mod_group == g, ]
  1121. df$min_mod <- rep(min(ev_all$modulation), K)
  1122. df$max_mod <- rep(max(ev_all$modulation), K)
  1123. }
  1124. df
  1125. })
  1126. rownames_out <- NULL
  1127. do.call(rbind, out)
  1128. }
  1131. #' Plot fMRI peri-stimulus time courses
  1132. #'
  1133. #' This function plots average BOLD response around specified events for a single ROI
  1134. #' by using FIR based event estimation, all event_types in events are taken into account in the FIR.
  1135. #' Posterior predictives can be overlaid via the `post_predict` argument.
  1136. #'
  1137. #' @param timeseries A data frame with columns 'subjects', 'run', 'time', and one ROI measurement column.
  1138. #' @param events A data frame with columns 'subjects', 'run', 'onset', 'duration', 'event_type', and 'modulation'.
  1139. #' @param event_type Character string specifying which `event_type` in `events` to plot.
  1140. #' @param high_pass Logical indicating whether to apply high-pass filtering.
  1141. #' Alternatively, specifying 'add' adds the regressors to the design matrix in the FIR.
  1142. #' The choice here should be the same as the choice for `convolve_design_matrix`
  1143. #' @param high_pass_model Character indicating which type of high-pass filtering to apply ('cosine', 'poly')
  1144. #' @param post_predict Optional posterior predictive samples data frame (not shown in examples).
  1145. #' @param posterior_args Named list of graphical parameters for posterior predictive lines.
  1146. #' @param legend_pos Position of the legend. Default: "topleft".
  1147. #' @param layout Panel layout matrix for multiple modulation groups. NULL leaves current layout
  1148. #' @param n_cores Number of cores to calculate FIR across subjects with.
  1149. #' @param ... Additional graphical parameters passed to plotting functions (e.g., col, lwd, lty).
  1150. #'
  1151. #' @return NULL. Produces plots as a side-effect.
  1152. #' @export
  1153. #'
  1154. #' @examples
  1155. #' ts <- data.frame(
  1156. #' subjects = rep(1, 100),
  1157. #' run = rep(1, 100),
  1158. #' time = seq(0, 99),
  1159. #' ROI = rnorm(100)
  1160. #' )
  1161. #' events <- data.frame(
  1162. #' subjects = rep(1, 5),
  1163. #' run = rep(1, 5),
  1164. #' onset = c(10, 30, 50, 70, 90),
  1165. #' event_type = rep("A", 5),
  1166. #' modulation = rep(1, 5),
  1167. #' duration = rep(0.5, 5)
  1168. #' )
  1169. #' plot_fmri(ts, events = events, event_type = "A")
  1170. plot_fmri <- function(timeseries, post_predict = NULL, events, event_type,
  1171. high_pass = TRUE, high_pass_model = "cosine",
  1172. posterior_args = list(),
  1173. legend_pos = "topleft", layout = NA,
  1174. n_cores = 1, ...) {
  1176. posterior_args <- add_defaults(posterior_args, col = "darkgreen", lwd = 2)
  1177. plot_args <- add_defaults(list(...), col = "black", lwd = 2,
  1178. xlab = "time (s)", ylab = "BOLD response")
  1180. ROI_col <- setdiff(names(timeseries),
  1181. c("run", "subjects", "time", "postn"))
  1183. ts <- get_fir_lines(timeseries, events, event_type,
  1184. high_pass = high_pass,
  1185. high_pass_model = high_pass_model)
  1187. ## posterior predictive overlay
  1188. if(!is.null(post_predict)) {
  1189. pp_split <- split(post_predict, post_predict$postn)
  1190. pp_lines <- auto_mclapply(pp_split, get_fir_lines, events, event_type,
  1191. high_pass = high_pass,
  1192. high_pass_model = high_pass_model, mc.cores = n_cores)
  1193. pp_df <- do.call(rbind, pp_lines)
  1195. qfun <- function(x) quantile(x, c(.025, .5, .975))
  1196. agg <- aggregate(avg_signal ~ mod_group + time,
  1197. data = pp_df, FUN = qfun)
  1198. res_df <- data.frame(agg[, 1:2], agg$avg_signal)
  1199. names(res_df)[3:5] <- c("p025", "p50", "p975")
  1200. }
  1202. ## layout
  1203. n_plots <- length(unique(ts$mod_group))
  1204. if(n_plots == 1){
  1205. insert <- ""
  1206. } else{
  1207. insert <- "Q"
  1208. }
  1210. if(!is.null(layout)){
  1211. oldpar <- par(no.readonly = TRUE)
  1212. on.exit(par(oldpar))
  1213. }
  1214. if (any(is.na(layout))) {
  1215. par(mfrow = coda_setmfrow(Nchains = 1, Nparms = n_plots, nplots = 1))
  1216. } else {
  1217. par(mfrow = layout)
  1218. }
  1220. ## loop over modulation groups
  1221. for(mod in unique(ts$mod_group)){
  1222. if(!is.null(post_predict)){
  1223. ylim <- c(min(ts$avg_signal, res_df$p025), max(ts$avg_signal, res_df$p975))
  1224. } else{
  1225. ylim <- range(ts$avg_signal)
  1226. }
  1227. idx_ts <- ts$mod_group == mod
  1228. tmp_plot_args <- add_defaults(plot_args, ylim = ylim,
  1229. main = paste0(ROI_col, " - ", event_type, ifelse(n_plots == 1, "", paste0(": ", insert, mod))))
  1230. do.call(plot, c(list(ts$time[idx_ts], ts$avg_signal[idx_ts], type = "l"), fix_dots_plot(tmp_plot_args)))
  1231. abline(h = 0, lty = 2)
  1232. if(!is.null(post_predict)){
  1233. idx_pp <- res_df$mod_group == mod
  1234. do.call(lines, c(list(res_df$time[idx_pp], res_df$p50[idx_pp]), fix_dots_plot(posterior_args)))
  1235. polygon_args <- posterior_args
  1236. polygon_args$col <- adjustcolor(polygon_args$col, alpha.f = .2)
  1237. do.call(polygon, c(list(
  1238. x = c(res_df$time[idx_pp], rev(res_df$time[idx_pp])),
  1239. y = c(res_df$p025[idx_pp], rev(res_df$p975[idx_pp])), border = NA), fix_dots_plot(polygon_args)))
  1240. }
  1241. dots_legend <- plot_args
  1242. dots_legend$col <- c(plot_args$col, posterior_args$col)
  1243. dots_legend$lty <- c(plot_args$lty, posterior_args$lty)
  1244. dots_legend$lwd <- c(plot_args$lwd, posterior_args$lwd)
  1245. if(is.null(post_predict)){
  1246. legend_in <- "data"
  1247. } else{
  1248. legend_in <- c("data", "posterior")
  1249. }
  1250. do.call(legend, c(list(legend_pos, legend = legend_in, bty = "n"), fix_dots(dots_legend, legend)))
  1251. }
  1252. }
  1256. # Utility functions -------------------------------------------------------
  1260. make_mri_sampling_design <- function(design, sampled_p_names){
  1261. out <- list()
  1262. expand <- 1:nrow(design)
  1263. rownames(design) <- NULL
  1264. for(i in 1:ncol(design)){
  1265. out[[i]] <- design[,i, drop = F]
  1266. attr(out[[i]], "expand") <- expand
  1267. }
  1268. names(out) <- colnames(design)
  1269. no_design <- sampled_p_names[!sampled_p_names %in% colnames(design)]
  1270. for (nam in no_design){
  1271. mat <- matrix(1, 1, 1)
  1272. colnames(mat) <- nam
  1273. out[[nam]] <- mat
  1274. attr(out[[nam]], "expand") <- rep(1, nrow(design))
  1275. }
  1276. return(out)
  1277. }
  1280. make_data_wrapper_MRI <- function(parameters, data, design){
  1281. data_list <- list()
  1282. for(i in 1:nrow(parameters)){
  1283. data_tmp <- data[data$subjects == unique(data$subjects)[i],]
  1284. data_tmp$subjects <- factor(data_tmp$subjects)
  1285. attr(data_tmp, "designs") <- design$fMRI_design[[i]]
  1287. data_list[[i]] <- make_data_fMRI(parameters[i,,drop = F], design$model, data_tmp, design)
  1288. }
  1289. return(do.call(rbind, data_list))
  1290. }
  1292. make_data_fMRI <- function(parameters, model, data, design, ...){
  1293. # if(is.null(attr(design, "design_matrix"))){
  1294. # stop("for fMRI simulation the original design needs to be passed to the simulation function")
  1295. # }
  1296. pars <- get_pars_matrix_oo(parameters, data, model())
  1297. data[, !colnames(data) %in% c("subjects", "run", "time", "trials")] <- model()$rfun(pars)
  1298. return(data)
  1299. }
  1301. add_design_fMRI_predict <- function(design, emc){
  1302. design$fMRI_design <- lapply(emc[[1]]$data, function(x) return(attr(x,"designs")))
  1303. return(design)
  1304. }

MRI.R at commit beab948, under GPL-3.0 · at the source

Overview

Authors: Niek Stevenson1, Steven Miletić1,2, Birte U. Forstmann1
ORCID iDs: Niek Stevenson
  1. University of Amsterdam, Amsterdam, The Netherlands
  2. Leiden University, Leiden, The Netherlands
Institutions: University of Amsterdam (Netherlands); Leiden University (Netherlands)
Journal: Imaging neuroscience (Cambridge, Mass.), volume 4, article IMAG.a.1272
Dates: received 8 July 2025; accepted 8 May 2026; published online 15 June 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1162/imag.a.1272 · PMID 42312084 · PMCID PMC13271149 · OpenAlex W7162085102
Open access: diamond, a free copy (OpenAlex)
Status: code verified
Categories: fMRI (modality), methods / tools (subfield)
Methods: Connectivity, Preprocessing, fMRI & imaging, Statistics
Keywords: joint modeling, fMRI, hierarchical Bayes, evidence accumulation models, diffusion decision model
Journal subjects: Software Toolbox
Topic: Functional Brain Connectivity Studies (Cognitive Neuroscience, Neuroscience), according to OpenAlex
Funding: EC | Horizon Europe | Excellent Science | HORIZON EUROPE European Research Council (ERC) (864750); Nederlandse Organisatie voor Wetenschappelijk Onderzoek (016.vici.185.052)
Citations: not cited yet (Europe PMC); 86 references in the paper
Research resources: RRID:SCR_001362, RRID:SCR_001847, RRID:SCR_002438, which is based on Nipype 1.7.0 RRID:SCR_002502, RRID:SCR_002823, distributed with ANTs 2.3.3 RRID:SCR_004757, RRID:SCR_005927, RRID:SCR_008796, RRID:SCR_016216

Abstract

Understanding how neural activity relates to behavior remains a central challenge in cognitive neuroscience. Joint modeling offers a principled method by simultaneously fitting behavioral and fMRI data and estimating the relations between them, accounting for measurement error, inter-individual variability, and shared uncertainty. Here, we present a joint modeling toolbox built in the R package, EMC2, which streamlines the behavioral, neural, and joint model estimation. We apply the toolbox to a perceptual decision-making task performed in an fMRI scanner, demonstrating how to specify behavioral and neural design matrices, set priors, construct brain–behavioral links, and conduct group-level tests and individual-difference analyses. Through the hierarchical modeling framework, joint estimation stabilizes weakly identified parameters and mitigates attenuation bias in brain–behavior correlations. Our toolbox, accompanied by practical code examples, provides an easily accessible introduction for researchers to adopt joint models and to derive richer, more reliable insights into the cognitive and neural mechanisms underlying human behavior.

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 15 matches between paragraphs and lines of code.

ampl-psych/EMC2

License: GPL-3.0
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: beab948d283cfff25139de7f1a2dae11839cfddd, 14 July 2026
Languages: R (77), C++ (20), C/C++ (20), JavaScript (13)
Size: 429 files, 130 scripts
Software Heritage: not archived
Found in: “Data and Code Availability”
Holds: README, license file, environment (DESCRIPTION), tests, continuous integration, documentation, 5 notebooks
Not found: CITATION.cff
Tools: psych (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
132 files
At the source: github.com/ampl-psych/EMC2

OSF jxhcp

License: none: the authors keep all their rights
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Size: 1 file
Software Heritage: not checked
Found in: “Data and Code Availability”
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Availability: 1 check, the latest on 27 September 2026: the link answers (HTTP 200)
  • 27 September 2026: the link answers (HTTP 200)
At the source: osf.io/jxhcp/

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:

  • 2 repositories of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 130 scripts, each with its path and the digest of its content;
  • 15 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

Datasets cited

Data and Code Availability

The source code for the EMC2 package can be found at: https://github.com/ampl-psych/EMC2. The accompanying tutorial code can be found at: https://osf.io/jxhcp/. The archival data from Pereira et al. (2020) were retrieved from: https://openneuro.org/datasets/ds002158/versions/1.0.2.

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 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, pages, dates, 3 authors, 5 keywords, 2 funders, 86 references, 9 RRIDs.

Cite

This paper

Stevenson, N., Miletić, S., & Forstmann, B. U. (2026). Building bridges between brain and behavior: An open-source toolbox for joint modeling with fMRI. Imaging neuroscience (Cambridge, Mass.), 4, IMAG.a.1272. https://doi.org/10.1162/imag.a.1272

BibTeX

@article{stevenson2026building,
author = {Stevenson, Niek and Miletić, Steven and Forstmann, Birte U.},
title = {{Building bridges between brain and behavior: An open-source toolbox for joint modeling with fMRI}},
journal = {Imaging neuroscience (Cambridge, Mass.)},
year = {2026},
month = jun,
volume = {4},
pages = {IMAG.a.1272},
publisher = {MIT Press},
issn = {2837-6056},
doi = {10.1162/imag.a.1272},
url = {https://doi.org/10.1162/imag.a.1272},
pmid = {42312084},
pmcid = {PMC13271149}
}

RIS

TY - JOUR
AU - Stevenson, Niek
AU - Miletić, Steven
AU - Forstmann, Birte U.
TI - Building bridges between brain and behavior: An open-source toolbox for joint modeling with fMRI
T2 - Imaging neuroscience (Cambridge, Mass.)
J2 - Imaging Neurosci (Camb)
PY - 2026
DA - 2026/06/15
VL - 4
SP - IMAG.a.1272
SN - 2837-6056
PB - MIT Press
DO - 10.1162/imag.a.1272
UR - https://doi.org/10.1162/imag.a.1272
LA - en
ER -

CSL-JSON

{
"id": "10.1162/imag.a.1272",
"type": "article-journal",
"title": "Building bridges between brain and behavior: An open-source toolbox for joint modeling with fMRI",
"container-title": "Imaging neuroscience (Cambridge, Mass.)",
"author": [
{
"family": "Stevenson",
"given": "Niek"
},
{
"family": "Miletić",
"given": "Steven"
},
{
"family": "Forstmann",
"given": "Birte U."
}
],
"container-title-short": "Imaging Neurosci (Camb)",
"volume": "4",
"page": "IMAG.a.1272",
"DOI": "10.1162/imag.a.1272",
"PMID": "42312084",
"PMCID": "PMC13271149",
"ISSN": "2837-6056",
"publisher": "MIT Press",
"URL": "https://doi.org/10.1162/imag.a.1272",
"language": "en",
"issued": {
"date-parts": [
[
2026,
6,
15
]
]
}
}

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.1523/eneuro.0370-25.2026 [code]
Pretraining for Large-Scale Functional Connectome Fingerprinting Supports Generalization and Transfer Learning in Functional Neuroimaging.
Journal: eNeuro
In common: fMRI, methods / tools, 14 references
[2] doi:10.1371/journal.pbio.3003684 [code]
The retrieval of previously learned motor memories is facilitated by the reinstatement of default mode network manifold structures.
Journal: PLoS biology
In common: fMRI, 13 references
[3] doi:10.7554/elife.103956 [code]
Human brain-wide activation of sleep rhythms.
Journal: eLife
In common: 13 references
[4] doi:10.7554/elife.107955
Multiple event segmentation mechanisms in the human brain.
Journal: eLife
In common: 13 references
[5] doi:10.1016/j.celrep.2026.117454 [code]
Model-based and model-free valuation signals in the human brain vary markedly in relation to individual differences in behavioral control.
Journal: Cell reports
In common: 13 references
[6] doi:10.1038/s41597-026-07323-y [code]
A densely sampled fMRI dataset for investigating food valuation.
Journal: Scientific data
In common: fMRI, methods / tools, 12 references
[7] doi:10.1371/journal.pbio.3003755 [code]
Action information is integrated into entorhinal representations of conceptual space and is reflected in eye movements.
Journal: PLoS biology
In common: 12 references
[8] doi:10.1523/jneurosci.2178-25.2026 [code]
Increased Attentive Use Is Linked to More Idiosyncratic Functional Connections.
Journal: The Journal of neuroscience : the official journal of the Society for Neuroscience
In common: fMRI, 11 references
[9] doi:10.1162/imag.a.1209
Naturalistic movie viewing is an effective functional localizer of the fusiform face area in adolescents with and without autism.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: fMRI, 10 references
[10] doi:10.1038/s41597-026-07540-5
A multimodal epilepsy dataset of paired 3-Tesla and 7-Tesla MRI and intracranial EEG.
Journal: Scientific data
In common: methods / tools, 10 references

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.