Calculating Intersubject Correlations
Intersubject correlation (ISC) is a measure used to assess the degree of synchrony among participants—that is, how similarly their responses change over time. ISC has been widely applied in fMRI and EEG research, particularly in naturalistic contexts, to quantify shared neural activity. More recently, ISC has been extended to webcam-based eye-tracking, allowing researchers to examine shared gaze patterns during dynamic stimuli.
This vignette demonstrates how to calculate ISC using two commonly used approaches:
- Pairwise ISC Approach
In this method, each participant’s data is correlated with every other participant’s data, generating a correlation matrix across all pairs.
Correlations are transformed using Fisher’s z to r transformation and averaged and then transformed back into r resulting in one ISC for each participant.
- Leave-One-Out ISC Approach - Here, each participant’s data is correlated with the mean time series of all other participants, excluding their own data. - This results in one correlation value for each participant.
To compute ISC, your data should be structured as a p Ă— t matrix, where:
- p represents the number of participants
- t represents time points
Before highlighting the function calculate_isc() I want to break down each step that goes into calculating the ISC.
Pairwise ISC
Below, we generate a 10 Ă— 1000 matrix
library(webgazeR)
set.seed(123)
data_matrix <- matrix(rnorm(10 * 1000), nrow = 1000, ncol = 10)
colnames(data_matrix) <- paste0("P", 1:10)
data_matrix |>
head() P1 P2 P3 P4 P5 P6
[1,] -0.56047565 -0.99579872 -0.5116037 -0.15030748 0.1965498 -0.4941739
[2,] -0.23017749 -1.03995504 0.2369379 -0.32775713 0.6501132 1.1275935
[3,] 1.55870831 -0.01798024 -0.5415892 -1.44816529 0.6710042 -1.1469495
[4,] 0.07050839 -0.13217513 1.2192276 -0.69728458 -1.2841578 1.4810186
[5,] 0.12928774 -2.54934277 0.1741359 2.59849023 -2.0261096 0.9161912
[6,] 1.71506499 1.04057346 -0.6152683 -0.03741501 2.2053261 0.3351310
P7 P8 P9 P10
[1,] -0.6992281 -1.6180367 0.5110004 1.9315759
[2,] 0.9964515 0.3791812 1.8079928 -0.6164747
[3,] -0.6927454 1.9022505 -1.7026150 -0.5625675
[4,] -0.1034830 0.6018743 0.2874488 -0.9899631
[5,] 0.6038661 1.7323497 -0.2691142 2.7312276
[6,] -0.6080450 -0.1468696 -0.3795247 -0.7216662We use the cor() function to compute pairwise correlations between each participant’s time-binned data and that of every other participant.
We specify use = “pairwise.complete.obs” to ensure that:
If any missing values (NA) exist for a participant, those time points are excluded from the correlation for all comparisons involving that participant.
This allows us to retain all participants in the dataset while ensuring that correlations are computed using only available data.
correlation_matrix <- cor(data_matrix, use = "pairwise.complete.obs")Since the diagonal of the correlation matrix represents self-correlations (i.e., each participant correlated with themselves), we replace those values with NA, as they are not meaningful for intersubject correlation (ISC):
diag(correlation_matrix) <- NAcorrelation_matrix P1 P2 P3 P4 P5
P1 NA 0.086479441 -0.019329537 -0.0029947104 0.0322973377
P2 0.086479441 NA 0.026503334 -0.0070290755 0.0924646624
P3 -0.019329537 0.026503334 NA 0.0505608499 -0.0150195347
P4 -0.002994710 -0.007029076 0.050560850 NA -0.0004005193
P5 0.032297338 0.092464662 -0.015019535 -0.0004005193 NA
P6 -0.018267991 -0.055292983 -0.033003346 0.0438457353 0.0138042312
P7 -0.034340849 0.024316458 0.005282974 -0.0200855749 -0.0090534146
P8 -0.030156049 0.031077180 -0.025946806 -0.0055278737 -0.0777495779
P9 -0.009370688 0.022155213 0.058417535 -0.0064371559 0.0470798891
P10 0.022344953 0.022481121 0.029148322 0.0459001272 -0.0047210318
P6 P7 P8 P9 P10
P1 -0.0182679912 -0.034340849 -0.030156049 -0.0093706879 0.022344953
P2 -0.0552929834 0.024316458 0.031077180 0.0221552128 0.022481121
P3 -0.0330033458 0.005282974 -0.025946806 0.0584175347 0.029148322
P4 0.0438457353 -0.020085575 -0.005527874 -0.0064371559 0.045900127
P5 0.0138042312 -0.009053415 -0.077749578 0.0470798891 -0.004721032
P6 NA -0.016497452 -0.016411441 -0.0006214018 -0.010177844
P7 -0.0164974516 NA 0.037173511 0.0791431861 0.002523005
P8 -0.0164114415 0.037173511 NA -0.0422433679 0.010503101
P9 -0.0006214018 0.079143186 -0.042243368 NA -0.006174781
P10 -0.0101778437 0.002523005 0.010503101 -0.0061747810 NANext, we apply Fisher’s Z-transformation using atanh() to convert the correlation coefficients into Z-scores before averaging.
Averaging correlation scores directly is bad practice because:
- Correlations are bounded between -1 and 1, which skews their distribution.
- Fisher’s Z-transformation normalizes correlation values, making them more suitable for averaging.
You can also choose to use the median instead (if this is the case you do not need to apply Fisher transformation)
z_values <- atanh(correlation_matrix) # Fisher's Z transformationOnce transformed, we compute the average ISC per participant and transform back to r.
isc_values <- tanh(rowMeans(z_values, na.rm = TRUE)) # Convert back to rThe calculate_isc() function with method = "pairwise" does all this for us:
isc_pairwise <- calculate_isc(data_matrix, method = "pairwise")
isc_pairwise P1 P2 P3 P4 P5 P6
0.002985107 0.027061018 0.008523871 0.010880947 0.008761466 -0.010296001
P7 P8 P9 P10
0.007625601 -0.013271741 0.015797990 0.012429909 You can also get the full correlation matrix back with return_matrix = TRUE:
isc_matrix <- calculate_isc(data_matrix, method = "pairwise", return_matrix = TRUE)
isc_matrix P1 P2 P3 P4 P5
P1 NA 0.086479441 -0.019329537 -0.0029947104 0.0322973377
P2 0.086479441 NA 0.026503334 -0.0070290755 0.0924646624
P3 -0.019329537 0.026503334 NA 0.0505608499 -0.0150195347
P4 -0.002994710 -0.007029076 0.050560850 NA -0.0004005193
P5 0.032297338 0.092464662 -0.015019535 -0.0004005193 NA
P6 -0.018267991 -0.055292983 -0.033003346 0.0438457353 0.0138042312
P7 -0.034340849 0.024316458 0.005282974 -0.0200855749 -0.0090534146
P8 -0.030156049 0.031077180 -0.025946806 -0.0055278737 -0.0777495779
P9 -0.009370688 0.022155213 0.058417535 -0.0064371559 0.0470798891
P10 0.022344953 0.022481121 0.029148322 0.0459001272 -0.0047210318
P6 P7 P8 P9 P10
P1 -0.0182679912 -0.034340849 -0.030156049 -0.0093706879 0.022344953
P2 -0.0552929834 0.024316458 0.031077180 0.0221552128 0.022481121
P3 -0.0330033458 0.005282974 -0.025946806 0.0584175347 0.029148322
P4 0.0438457353 -0.020085575 -0.005527874 -0.0064371559 0.045900127
P5 0.0138042312 -0.009053415 -0.077749578 0.0470798891 -0.004721032
P6 NA -0.016497452 -0.016411441 -0.0006214018 -0.010177844
P7 -0.0164974516 NA 0.037173511 0.0791431861 0.002523005
P8 -0.0164114415 0.037173511 NA -0.0422433679 0.010503101
P9 -0.0006214018 0.079143186 -0.042243368 NA -0.006174781
P10 -0.0101778437 0.002523005 0.010503101 -0.0061747810 NABy default, Pearson correlations are used. You can switch to Spearman with the cor_method argument:
isc_spearman <- calculate_isc(data_matrix, method = "pairwise", cor_method = "spearman")
isc_spearman P1 P2 P3 P4 P5 P6
-0.003115793 0.024277967 0.002729341 0.008779896 0.010114176 -0.011083150
P7 P8 P9 P10
0.007600397 -0.018210706 0.011954421 0.010942500 Leave-One-Out ISC
In the leave-one-out approach, we compute the correlation between each participant’s time series and the average time series of all other participants, excluding their own data.
First, we extract the time series for a single participant:
# Extract the full time series for the current subject
subject_time_series <- data_matrix[, 1]and calculate the mean of the other subjects data points, excluding the current subject
leave_one_out_mean <- rowMeans(data_matrix[, -1, drop = FALSE], na.rm = TRUE)We then compute the correlation between the participant’s time series and the leave-one-out mean:
leave_one_out_isc <- cor(subject_time_series,
leave_one_out_mean,
use = "pairwise.complete.obs")Since this needs to be done for every participant, we use calculate_isc() to automate the process:
isc_loo <- calculate_isc(data_matrix, method = "leave-one-out")
isc_loo P1 P2 P3 P4 P5 P6
0.008891827 0.080509746 0.025626567 0.031021518 0.026597776 -0.029482014
P7 P8 P9 P10
0.023592539 -0.038321505 0.046317774 0.036036027 Time-Windowed ISC
The time_window_isc() function computes ISC within sliding time windows rather than across the full time series. This is useful when you want to examine how synchrony changes over time—for example, whether participants’ gaze patterns converge during specific portions of a stimulus.
The function uses the leave-one-out approach within each window: for each participant, it correlates their time series against the mean of all other participants within each window. Fisher z-transformed correlations are averaged across windows to produce a single ISC value per participant.
isc_windowed <- time_window_isc(data_matrix, window_size = 50, step = 25)
isc_windowed P1 P2 P3 P4 P5 P6
0.011415919 0.076909638 0.022717568 0.028245792 0.036268444 -0.021398453
P7 P8 P9 P10
0.002853639 -0.039911715 0.048524607 0.036630041 Key parameters:
window_size: The number of time points in each window (default: 10).step: The step size between consecutive windows (default: 1). Use a larger step for non-overlapping or less overlapping windows.min_overlap: Minimum number of non-NA paired observations required to compute a correlation within a window (default: 3).
To inspect ISC at each window, set return_per_window = TRUE:
isc_detail <- time_window_isc(data_matrix, window_size = 50, step = 25,
return_per_window = TRUE)
# Per-subject ISC summary
isc_detail$isc P1 P2 P3 P4 P5 P6
0.011415919 0.076909638 0.022717568 0.028245792 0.036268444 -0.021398453
P7 P8 P9 P10
0.002853639 -0.039911715 0.048524607 0.036630041 # Per-window Fisher z values
head(isc_detail$per_window) window_start window_end subject z
1 1 50 P1 -0.13361077
2 26 75 P1 -0.07715429
3 51 100 P1 -0.03100864
4 76 125 P1 -0.04409713
5 101 150 P1 -0.14190318
6 126 175 P1 0.06666306