Mastering Apply: From Matrices to Multidimensional Neuroimaging Data

Mastering Apply: From Matrices to Multidimensional Neuroimaging Data Dr. Mowinckel 2025-11-01 In my first apply post, I briefly showed apply() converting columns to character, but I glossed over what makes apply() special. Here’s the thing: apply() is actually quite different from sapply() and mapply(). While those work on lists and vectors, apply() is designed specifically for matrices and arrays. Pattern recognition Use apply() when you have: • Matrix or array data (not data.frame - use dplyr for those) • Need to collapse one or more dimensions by applying a function • Want to apply the same operation across rows, columns, or higher dimensions • Working with multidimensional scientific data (neuroimaging, climate models, genomics, etc.) The syntax is always: apply(X, MARGIN, FUN, ...) Where: - Xis your matrix or array - MARGIN specifies which dimension(s) to preserve (1 = rows, 2 = columns, 3+ = higher dimensions) - FUN is the function to apply - ... are additional arguments passed to FUN Think of MARGIN as telling R “keep these dimensions, collapse everything else.” That’s where the magic happens. I’m going to be honest - when I first encountered apply(), the MARGIN parameter confused the heck out of me. “Why do I need to specify margins? What does that even mean?” But once it clicked, I realized how powerful this function is for working with structured data. Understanding apply() properly will save you from writing many unnecessary loops when working with matrix-like data. 1 The MARGIN parameter explained The MARGIN parameter tells apply() which dimension to “collapse” by applying your function. I know that sounds abstract, so let me show you what I mean. It’s much easier to understand with examples than with explanations. Note: The MARGIN parameter can be either numeric (e.g., 1for rows, 2for columns) or a character vector specifying the names of the dimensions (if your array or matrix has named dimensions). For example, if your array has dimnames, you can use MARGIN = "Subject" instead of its numeric index. Let’s start with a simple matrix - just numbers arranged in rows and columns: # Create a 4x3 matrix (4 rows, 3 columns) test_matrix <- matrix(1:12,nrow = 4,ncol = 3) test_matrix [,1] [,2] [,3] [1,] 1 5 9 [2,] 2 6 10 [3,] 3 7 11 [4,] 4 8 12 We’ve got a matrix with 4 rows and 3 columns. Each row could represent a different observation, and each column could be a different variable. Now, let’s see what happens when we apply sum() across different margins: # MARGIN = 1: apply across rows (each row becomes one value) apply(test_matrix, MARGIN = 1,FUN = sum) [1] 15 18 21 24 Notice what happened? We got 4 values - one for each row. MARGIN = 1 means “for each row, apply the function and give me one result.” So R took each row, summed up all the values in that row, and gave us back a vector of 4 sums. Now let’s try the other direction: # MARGIN = 2: apply across columns (each column becomes one value) apply(test_matrix, MARGIN = 2,FUN = sum) 2 [1] 10 26 42 This time we got 3 values - one for each column! MARGIN = 2 means “for each column, apply the function and give me one result.” Here’s how I think about it: - MARGIN = 1: “Process each row” → Number of rows determines output length - MARGIN = 2: “Process each column” → Number of columns determines output length The margin you specify is the dimension that gets preserved in your output. Everything else gets collapsed by your function. Real-world example: Student grades Okay, let’s move beyond toy examples. Imagine Each row is a student, each column is an assignment. This is actually a really common data structure in many fields - observations in rows, variables in columns. # Rows = students, Columns = assignments # styler: off grades <- matrix( c( 85,92,78,88,# Student 1 90,85,95,87,# Student 2 78,83,80,85,# Student 3 95,88,92,90,# Student 4 82,79,84,81 # Student 5 ), nrow = 5, byrow = TRUE ) colnames(grades) <- c("Quiz1","Quiz2","Midterm","Final") rownames(grades) <- paste("Student",1:5) grades Quiz1 Quiz2 Midterm Final Student 1 85 92 78 88 Student 2 90 85 95 87 Student 3 78 83 80 85 Student 4 95 88 92 90 Student 5 82 79 84 81 3 Now we can easily answer questions like “What’s each student’s average?” and “Which assignment was hardest?” To get each student’s average, we want to work across columns for each row (sum up all their assignments and divide by the number of assignments): # Student averages (across columns for each row) student_averages <- apply(grades, MARGIN = 1,FUN = mean) student_averages Student 1 Student 2 Student 3 Student 4 Student 5 85.75 89.25 81.50 91.25 81.50 See how clean that is? No loops, no indexing nightmares. Just “hey R, for each student (row), calculate the mean.” If you want to use the names of the dimensions instead of the numbers, you will need to give the dimensions names first. We have so far given each column and each row a name, but we haven’t named the row and column themselves. This is done by setting the dimnames attribute of the matrix: dimnames(grades) [[1]] [1] "Student 1" "Student 2" "Student 3" "Student 4" "Student 5" [[2]] [1] "Quiz1" "Quiz2" "Midterm" "Final" dimnames(grades) <- list( Student = rownames(grades), Assignment = colnames(grades) ) dimnames(grades) $Student [1] "Student 1" "Student 2" "Student 3" "Student 4" "Student 5" $Assignment [1] "Quiz1" "Quiz2" "Midterm" "Final" 4 apply(grades, MARGIN = "Student",FUN = mean) Student 1 Student 2 Student 3 Student 4 Student 5 85.75 89.25 81.50 91.25 81.50 What about assignment difficulty? We can look at the average score for each assignment by working down the columns: # Assignment averages (across rows for each column) assignment_averages <- apply(grades, MARGIN = 2,FUN = mean) assignment_averages Quiz1 Quiz2 Midterm Final 86.0 85.4 85.8 86.2 Looks like Quiz 2 was the toughest (or maybe everyone was hungover that week?). The midterm was actually easier than the quizzes on average. This is so much cleaner than writing loops! Compare this to what you’d have to do with a for loop - you’d need to initialize empty vectors, loop through indices, store results… it’s just messy. Beyond mean: useful functions for apply The function you pass to apply() can be anything that makes sense for your data. Let’s explore some other common operations you might want to do with gradebook data. Want to know how consistent each student’s performance is? Standard deviation tells you that: # Standard deviation for each student (higher = more variable) apply(grades, 1, sd) Student 1 Student 2 Student 3 Student 4 Student 5 5.909033 4.349329 3.109126 2.986079 2.081666 Student 1 has the highest variability - maybe they’re really good at some things and struggle with others? Student 4 is super consistent. What about the range of scores for each assignment? 5 # Range (min to max) for each assignment apply(grades, 2, range) Assignment Quiz1 Quiz2 Midterm Final [1,] 78 79 78 81 [2,] 95 92 95 90 This gives us a matrix where the first row is minimums and second row is maximums for each assignment. Here’s a fun one - which assignment was each student’s worst? apply(grades, 1, which.min) Student 1 Student 2 Student 3 Student 4 Student 5 32122 These are column indices. We can make this more readable by looking up the actual assignment names: # Get the assignment names instead of indices colnames(grades)[apply(grades, 1, which.min)] [1] "Midterm" "Quiz2" "Quiz1" "Quiz2" "Quiz2" Now that’s informative! Three students struggled most with Quiz2. Custom functions with apply You’re not limited to built-in functions. You can pass your own functions too, which is where things get really powerful. Let’s say we want to calculate letter grades based on each student’s average. We need a function that takes a vector of scores, calculates the mean, and returns a letter: 6 get_letter_grade <- function(scores) { avg <- mean(scores) if (avg >= 90) { return("A") } if (avg >= 80) { return("B") } if (avg >= 70) { return("C") } if (avg >= 60) { return("D") } return("F") } # Letter grade for each student apply(grades, 1, get_letter_grade) Student 1 Student 2 Student 3 Student 4 Student 5 "B" "B" "B" "A" "B" One A and 4 B’s, impressive! Here’s something more interesting - let’s figure out if students are improving over time. We can calculate a simple correlation between time (assignment order) and scores: calculate_trend <- function(scores) { # Simple linear trend (positive = improving, negative = declining) x <- seq_along(scores) trend <- cor(x, scores) return(trend) } # Trend for each student (are they improving over time?) trends <- apply(grades, 1, calculate_trend) trends Student 1 Student 2 Student 3 Student 4 Student 5 -0.10923907 0.02968261 0.74740932 -0.47557147 0.12403473 7 Student 3 has the strongest positive trend (0.7474093) - they started rough but ended strong! Good for them. Student 4 has a negative trend (-0.4755715) - started great but faded. Maybe they need some extra help? The cool thing about custom functions is you can make them as complex as you need. Statistical tests, transformations, summarizations - whatever makes sense for your analysis. Working with 3D arrays Now we’re getting into territory where apply() really starts to shine. Let’s say instead of one class, you’re teaching multiple sections. Or maybe you’re running a longitudinal study with multiple measurement occasions. Suddenly you need a third dimension. Let me simulate something realistic: test scores for 5 students across 4 assignments in 3 different classes: # 3D array: 5 students × 4 assignments × 3 classes set.seed(42) scores_3d <- array( round(rnorm(60,mean = 85,sd = 10)), dim = c(5,4,3), dimnames = list( Student = paste("S",1:5,sep = ""), Assignment = c("Quiz1","Quiz2","Midterm","Final"), Class = c("Math","Science","English") ) ) # Look at the whole structure scores_3d , , Class = Math Assignment Student Quiz1 Quiz2 Midterm Final S1 99 84 98 91 S2 79 100 108 82 S3 89 84 71 58 S4 91 105 82 61 S5 89 84 84 98 , , Class = Science 8 Assignment Student Quiz1 Quiz2 Midterm Final S1 82 81 90 68 S2 67 82 92 77 S3 83 67 95 76 S4 97 90 79 61 S5 104 79 90 85 , , Class = English Assignment Student Quiz1 Quiz2 Midterm Final S1 87 89 88 88 S2 81 77 77 92 S3 93 99 101 86 S4 78 81 91 55 S5 71 92 86 88 This is now a 3-dimensional array. Think of it like a cube of numbers. We can slice it to look at just one class: # Just the Math class scores_3d[,, "Math"] Assignment Student Quiz1 Quiz2 Midterm Final S1 99 84 98 91 S2 79 100 108 82 S3 89 84 71 58 S4 91 105 82 61 S5 89 84 84 98 Now, with apply(), we can answer all sorts of questions by specifying which dimensions to collapse. Want each student’s overall average across all assignments and all classes? We keep dimension 1 (students) and collapse dimensions 2 (assignments) and 3 (classes): # Student averages across all assignments and classes # Keep dimension 1, collapse dimensions 2 and 3 apply(scores_3d, MARGIN = "Student",FUN = mean) 9 If you scan multiple subjects, you can imagine stacking these time series for each person, forming a “plane” of cubes—one row per subject, each row a sequence of 3D brain volumes over time. This is a helpful way to visualize multidimensional neuroimaging data, even if the actual storage in software may differ. 16 Figure 4: An illustration titled “5D MRI Data: Multi-Subject Time Series” showing a vertical stack of three rows. Each row represents a different Subject (the 5th dimension) and consists of a horizontal line of four translucent blue brain cubes (composed of voxels). Arrows connect the cubes within each row, representing the progression of Time (the 4th dimension). This arrangement visualizes 5 dimensions: X, Y, Z (spatial), Time, and Subject. So we have a cube of voxels changing over time, for multiple subjects. Each voxel is represented 17 as a cell in our multidimensional array. Analyzing this data without apply() would be a nightmare of nested loops and index juggling. In real life, we would do lots of stuff to normalize and smooth the data to reduce the immense noise inherent in this data. MRI analyses are complex and nuanced, but let me show you how apply() works with this kind of structure using actual data from a BIDS (Brain Imaging Data Standard) dataset I grabbed from OpenNeuro. The data were obtained by Mueckstein (2024), and are openly available for anyone to use (with credit, of course). This dataset comes from a dual-task training study where participants performed tasks requiring either compatible (visual+manual, auditory+vocal) or incompatible (visual+vocal, auditory+manual) modality pairings. It has pre and post training sessions, though we’ll only work with the pre-session data here. For our purposes, we’re using the anatomical scans and a simple auditory task to demonstrate how apply() works with real neuroimaging data. I’ve reduced the data to just 10 subjects to keep things manageable. Single subject Let’s start simple with loading in a single subject, and have a look at the MRI data. library(oro.nifti) # Path to our BIDS dataset bids_dir <- here::here("content/blog/2025/11-01_apply/data") Let’s load one subject’s functional run: # Load a functional run for subject 01 func_file <- file.path( bids_dir, "sub-101/ses-pre/func", "sub-101_ses-pre_task-DT_run-01_bold.nii.gz" ) fmri_data <- readNIfTI(func_file, reorient = FALSE) dim(fmri_data) [1] 64 64 37 133 This is a 4D array: 64 voxels in x, 64 in y, 37 in z, and 133 timepoints. That’s 64x64x37x133 = 298 individual numbers. You definitely don’t want to write nested loops for this! 18 One common analysis is calculating the mean activation for each voxel across time. Of course, real fMRI analysis is more complex than this, but this can illustrate how apply() works with 4D data. This tells you the average signal intensity at each location in the brain: # Calculate mean activation across time for each voxel # Keep x, y, z dimensions; collapse time (dimension 4) mean_activation <- apply( fmri_data, MARGIN = c(1,2,3), FUN = mean ) dim(mean_activation) [1] 64 64 37 # Visualize a slice of the mean activation image( mean_activation[,,20], # Slice at z=20 col = viridis::viridis(256), main = "Mean Activation - Subject 01 (Slice z=20)", xlab = "X-axis", ylab = "Y-axis" ) 19 Figure 5: Mean Activation across Voxels for Subject 01. Shows average signal intensity at each voxel location. Clear contour of a brain is visible. We went from 4D to 3D - we collapsed the time dimension and kept the spatial dimensions. Each voxel now has one number: its average signal across the whole scan. Similarly, we might want to know how variable the signal is at each location. High variability might indicate functional activity (or noise, depending on what you’re looking for): # Calculate temporal standard deviation for each voxel temporal_sd <- apply( fmri_data, MARGIN = c(1,2,3), FUN = sd ) 20 # Coefficient of variation (CV) - normalized measure of variability cv <- temporal_sd /mean_activation # Visualize a slice of the CV image( cv[,,20], # Slice at z=20 col = viridis::viridis(256), main = "Coefficient of Variation - Subject 01 (Slice z=20)", xlab = "X-axis", ylab = "Y-axis" ) 21 Figure 6: Coefficient of Variation across Voxels for Subject 01. Shows normalized variability of signal intensity at each voxel. Higher values indicate more variability relative to mean signal.Shows clear contour of a brain. The coefficient of variation is useful because it accounts for the baseline signal intensity. A voxel with high signal and high variability has a different meaning than a voxel with low signal and high variability. Sometimes you want to summarize across space instead of time. The “global signal” is the average intensity across all voxels at each timepoint: # Subject-level summary: mean signal across all voxels over time # Keep time dimension; collapse spatial dimensions (1, 2, 3) global_signal <- apply( fmri_data, 22 MARGIN = 4, FUN = mean ) # Plot the global signal library(ggplot2) data.frame( time = 1:length(global_signal), signal = global_signal )|> ggplot(aes(x = time, y = signal)) + geom_line(color = "steelblue",linewidth = 1.2)+ labs( title = "Global Signal - Subject 01", x = "Time (TRs)", y = "Mean Signal" )+ theme_minimal() 23 Figure 7: Global Signal over Time for Subject 01. Shows overall mean brain signal intensity at each timepoint. This gives you a sense of overall signal drift and scanner stability. You can see some slow drift over time, which is pretty typical in fMRI. Multiple subjects Now let’s scale up to multiple subjects. This is where neuroimaging analysis gets computationally intensive, but apply() keeps it manageable: fmri_files <- list.files( bids_dir, "ses-pre_task-DT_run-01_bold.nii.gz", 24 recursive = TRUE, full.names = TRUE ) # Extract subject IDs from filenames subjects <- gsub( ".*sub-(\\d+)_ses-.*", "\\1", basename(fmri_files) ) subjects [1] "101" "104" "105" "108" "109" "110" "111" "113" "114" "115" # Read in imaging files fmri_list <- lapply( fmri_files, readNIfTI, reorient = FALSE ) # Combine into 5D array: x, y, z, time, subjects fmri_array <- array( unlist(fmri_list), dim = c(dim(fmri_list[[1]]), length(fmri_list)) ) # Add dimension names dimnames(fmri_array) <- list( X = NULL, Y = NULL, Z = NULL, time = paste0("TR",seq_len(dim(fmri_array)[4])), Subject = paste0("sub-", subjects) ) dim(fmri_array) [1] 64 64 37 133 10 25