Repository navigation
Expand file tree
/
Copy pathisc.py
More file actions
1548 lines (1222 loc) · 62 KB
/
Copy pathisc.py
File metadata and controls
1548 lines (1222 loc) · 62 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
10BC8
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
# Copyright 2017 Intel Corporation
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
"""Intersubject correlation (ISC) analysis
Functions for computing intersubject correlation (ISC) and related
analyses (e.g., intersubject funtional correlations; ISFC), as well
as statistical tests designed specifically for ISC analyses.
The implementation is based on the work in [Hasson2004]_, [Kauppi2014]_,
[Simony2016]_, [Chen2016]_, and [Nastase2019]_.
.. [Chen2016] "Untangling the relatedness among correlations, part I:
nonparametric approaches to inter-subject correlation analysis at the
group level.", G. Chen, Y. W. Shin, P. A. Taylor, D. R. Glen, R. C.
Reynolds, R. B. Israel, R. W. Cox, 2016, NeuroImage, 142, 248-259.
https://doi.org/10.1016/j.neuroimage.2016.05.023
.. [Hasson2004] "Intersubject synchronization of cortical activity
during natural vision.", U. Hasson, Y. Nir, I. Levy, G. Fuhrmann,
R. Malach, 2004, Science, 303, 1634-1640.
https://doi.org/10.1126/science.1089506
.. [Kauppi2014] "A versatile software package for inter-subject
correlation based analyses of fMRI.", J. P. Kauppi, J. Pajula,
J. Tohka, 2014, Frontiers in Neuroinformatics, 8, 2.
https://doi.org/10.3389/fninf.2014.00002
.. [Simony2016] "Dynamic reconfiguration of the default mode network
during narrative comprehension.", E. Simony, C. J. Honey, J. Chen, O.
Lositsky, Y. Yeshurun, A. Wiesel, U. Hasson, 2016, Nature Communications,
7, 12141. https://doi.org/10.1038/ncomms12141
.. [Nastase2019] "Measuring shared responses across subjects using
intersubject correlation." S. A. Nastase, V. Gazzola, U. Hasson,
C. Keysers, 2019, Social Cognitive and Affective Neuroscience, 14,
667-685. https://doi.org/10.1093/scan/nsz037
"""
# Authors: Sam Nastase, Christopher Baldassano, Qihong Lu,
# Mai Nguyen, and Mor Regev
# Princeton University, 2018
import numpy as np
import logging
from scipy.spatial.distance import squareform
from itertools import combinations, permutations, product
from brainiak.fcma.util import compute_correlation
from brainiak.utils.utils import (array_correlation,
phase_randomize,
p_from_null,
_check_timeseries_input)
logger = logging.getLogger(__name__)
__all__ = [
"bootstrap_isc",
"compute_summary_statistic",
"isfc",
"isc",
"permutation_isc",
"phaseshift_isc",
"squareform_isfc",
"timeshift_isc",
]
MAX_RANDOM_SEED = 2**32 - 1
def isc(data, pairwise=False, summary_statistic=None, tolerate_nans=True):
"""Intersubject correlation
For each voxel or ROI, compute the Pearson correlation between each
subject's response time series and other subjects' response time series.
If pairwise is False (default), use the leave-one-out approach, where
correlation is computed between each subject and the average of the other
subjects. If pairwise is True, compute correlations between all pairs of
subjects. If summary_statistic is None, return N ISC values for N subjects
(leave-one-out) or N(N-1)/2 ISC values for each pair of N subjects,
corresponding to the upper triangle of the pairwise correlation matrix
(see scipy.spatial.distance.squareform). Alternatively, use either
'mean' or 'median' to compute summary statistic of ISCs (Fisher Z will
be applied if using mean). Input data should be a n_TRs by n_voxels by
n_subjects array (e.g., brainiak.image.MaskedMultiSubjectData) or a list
where each item is a n_TRs by n_voxels ndarray for a given subject.
Multiple input ndarrays must be the same shape. If a 2D array is supplied,
the last dimension is assumed to correspond to subjects. If only two
subjects are supplied, simply compute Pearson correlation (precludes
averaging in leave-one-out approach, and does not apply summary statistic).
When using leave-one-out approach, NaNs are ignored when computing mean
time series of N-1 subjects (default: tolerate_nans=True). Alternatively,
you may supply a float between 0 and 1 indicating a threshold proportion
of N subjects with non-NaN values required when computing the average time
series for a given voxel. For example, if tolerate_nans=.8, ISCs will be
computed for any voxel where >= 80% of subjects have non-NaN values,
while voxels with < 80% non-NaN values will be assigned NaNs. If set to
False, NaNs are not tolerated and voxels with one or more NaNs among the
N-1 subjects will be assigned NaN. Setting tolerate_nans to True or False
will not affect the pairwise approach; however, if a threshold float is
provided, voxels that do not reach this threshold will be excluded. Note
that accommodating NaNs may be notably slower than setting tolerate_nans to
False. Output is an ndarray where the first dimension is the number of
subjects or pairs and the second dimension is the number of voxels (or
ROIs). If only two subjects are supplied or a summary statistic is invoked,
the output is a ndarray n_voxels long.
The implementation is based on the work in [Hasson2004]_.
Parameters
----------
data : list or ndarray (n_TRs x n_voxels x n_subjects)
fMRI data for which to compute ISC
pairwise : bool, default: False
Whether to use pairwise (True) or leave-one-out (False) approach
summary_statistic : None or str, default: None
Return all ISCs or collapse using 'mean' or 'median'
tolerate_nans : bool or float, default: True
Accommodate NaNs (when averaging in leave-one-out approach)
Returns
-------
iscs : subjects or pairs by voxels ndarray
ISC for each subject or pair (or summary statistic) per voxel
"""
# Check response time series input format
data, n_TRs, n_voxels, n_subjects = _check_timeseries_input(data)
# No summary statistic if only two subjects
if n_subjects == 2:
logger.info("Only two subjects! Simply computing Pearson correlation.")
summary_statistic = None
# Check tolerate_nans input and use either mean/nanmean and exclude voxels
if tolerate_nans:
mean = np.nanmean
else:
mean = np.mean
data, mask = _threshold_nans(data, tolerate_nans)
# Compute correlation for only two participants
if n_subjects == 2:
# Compute correlation for each corresponding voxel
iscs_stack = array_correlation(data[..., 0],
data[..., 1])[np.newaxis, :]
# Compute pairwise ISCs using voxel loop and corrcoef for speed
elif pairwise:
# Swap axes for np.corrcoef
data = np.swapaxes(data, 2, 0)
# Loop through voxels
voxel_iscs = []
for v in np.arange(data.shape[1]):
voxel_data = data[:, v, :]
# Correlation matrix for all pairs of subjects (triangle)
iscs = squareform(np.corrcoef(voxel_data), checks=False)
voxel_iscs.append(iscs)
iscs_stack = np.column_stack(voxel_iscs)
# Compute leave-one-out ISCs
elif not pairwise:
# Loop through left-out subjects
iscs_stack = []
for s in np.arange(n_subjects):
# Correlation between left-out subject and mean of others
iscs_stack.append(array_correlation(
data[..., s],
mean(np.delete(data, s, axis=2), axis=2)))
iscs_stack = np.array(iscs_stack)
# Get ISCs back into correct shape after masking out NaNs
iscs = np.full((iscs_stack.shape[0], n_voxels), np.nan)
iscs[:, np.where(mask)[0]] = iscs_stack
# Summarize results (if requested)
if summary_statistic:
iscs = compute_summary_statistic(iscs,
summary_statistic=summary_statistic,
axis=0)[np.newaxis, :]
# Throw away first dimension if singleton
if iscs.shape[0] == 1:
iscs = iscs[0]
return iscs
def isfc(data, targets=None, pairwise=False, summary_statistic=None,
vectorize_isfcs=True, tolerate_nans=True):
"""Intersubject functional correlation (ISFC)
For each input voxel or ROI, compute the Pearson correlation between each
subject's response time series and all input voxels or ROIs in other
subjects. If a targets array is provided, instead compute ISFCs between
each input voxel time series and each voxel time series in targets across
subjects (resulting in asymmetric ISFC values). The targets array must have
the same number TRs and subjects as the input data. If pairwise is False
(default), use the leave-one-out approach, where correlation is computed
between each subject and the average of the other subjects. If pairwise is
True, compute correlations between all pairs of subjects. If a targets
array is provided, only the leave-one-out approach is supported. If
summary_statistic is None, return N ISFC values for N subjects (leave-one-
out) or N(N-1)/2 ISFC values for each pair of N subjects, corresponding to
the triangle of the correlation matrix (scipy.spatial.distance.squareform).
Alternatively, use either 'mean' or 'median' to compute summary statistic
of ISFCs (Fisher Z is applied if using mean). Input should be n_TRs by
n_voxels by n_subjects array (e.g., brainiak.image.MaskedMultiSubjectData)
or a list where each item is a n_TRs by n_voxels ndarray per subject.
Multiple input ndarrays must be the same shape. If a 2D array is supplied,
the last dimension is assumed to correspond to subjects. If only two
subjects are supplied, simply compute ISFC between these two subjects
(precludes averaging in leave-one-out approach, and does not apply summary
statistic). Returns vectorized upper triangle of ISFC matrices for each
subject or pair when vectorized_isfcs=True, or full (redundant) 2D ISFC
matrices when vectorized_isfcs=False. When using leave-one-out approach,
NaNs are ignored when computing mean time series of N-1 subjects (default:
tolerate_nans=True). Alternatively, you may supply a float between 0 and
1 indicating a threshold proportion of N subjects with non-NaN values
required when computing the average time series for a given voxel. For
example, if tolerate_nans=.8, ISCs will be computed for any voxel where
>= 80% of subjects have non-NaN values, while voxels with < 80% non-NaN
values will be assigned NaNs. If set to False, NaNs are not tolerated
and voxels with one or more NaNs among the N-1 subjects will be assigned
NaN. Setting tolerate_nans to True or False will not affect the pairwise
8
approa
10BC8
ch; however, if a threshold float is provided, voxels that do not
reach this threshold will be excluded. Note that accommodating NaNs may
be notably slower than setting tolerate_nans to False. Output is either
a tuple comprising condensed off-diagonal ISFC values and the diagonal
ISC values if vectorize_isfcs=True, or a single ndarray with shape
n_subjects (or n_pairs) by n_voxels by n_voxels 3D array if
vectorize_isfcs=False (see brainiak.isc.squareform_isfc). If targets array
is provided (yielding asymmetric ISFCs), output ISFCs are not vectorized,
resulting in an n_subjects by n_voxels by n_targets ISFC array. If
summary_statistic is supplied, output is collapsed along first dimension.
The implementation is based on the work in [Simony2016]_.
Parameters
----------
data : list or ndarray (n_TRs x n_voxels x n_subjects)
fMRI data for which to compute ISFC
targets : list or ndarray (n_TRs x n_voxels x n_subjects), optional
fMRI data to use as targets for ISFC
pairwise : bool, default: False
Whether to use pairwise (True) or leave-one-out (False) approach
summary_statistic : None or str, default: None
Return all ISFCs or collapse using 'mean' or 'median'
vectorize_isfcs : bool, default: True
Return tuple of condensed ISFCs and ISCs (True) or square (redundant)
ISFCs (False)
tolerate_nans : bool or float, default: True
Accommodate NaNs (when averaging in leave-one-out approach)
Returns
-------
isfcs : ndarray or tuple of ndarrays
ISFCs for each subject or pair (or summary statistic) per voxel pair
"""
# Check response time series input format
data, n_TRs, n_voxels, n_subjects = _check_timeseries_input(data)
# Check for optional targets input array
targets, t_n_TRs, t_n_voxels, t_n_subejcts, symmetric = (
_check_targets_input(targets, data))
if not symmetric:
pairwise = False
# Check tolerate_nans input and use either mean/nanmean and exclude voxels
if tolerate_nans:
mean = np.nanmean
else:
mean = np.mean
data, mask = _threshold_nans(data, tolerate_nans)
targets, targets_mask = _threshold_nans(targets, tolerate_nans)
# Handle just two subjects properly (for symmetric approach)
if symmetric and n_subjects == 2:
isfcs = compute_correlation(np.ascontiguousarray(data[..., 0].T),
np.ascontiguousarray(data[..., 1].T),
return_nans=True)
isfcs = (isfcs + isfcs.T) / 2
isfcs = isfcs[..., np.newaxis]
summary_statistic = None
logger.info("Only two subjects! Computing ISFC between them.")
# Compute all pairwise ISFCs (only for symmetric approach)
elif pairwise:
isfcs = []
for pair in combinations(np.arange(n_subjects), 2):
isfc_pair = compute_correlation(np.ascontiguousarray(
data[..., pair[0]].T),
np.ascontiguousarray(
targets[..., pair[1]].T),
return_nans=True)
if symmetric:
isfc_pair = (isfc_pair + isfc_pair.T) / 2
isfcs.append(isfc_pair)
isfcs = np.dstack(isfcs)
# Compute ISFCs using leave-one-out approach
else:
# Roll subject axis for loop
data = np.rollaxis(data, 2, 0)
targets = np.rollaxis(targets, 2, 0)
# Compute leave-one-out ISFCs
isfcs = [compute_correlation(np.ascontiguousarray(subject.T),
np.ascontiguousarray(mean(
np.delete(targets, s, axis=0),
axis=0).T),
return_nans=True)
for s, subject in enumerate(data)]
# Transpose and average ISFC matrices for both directions
isfcs = np.dstack([(isfc_matrix + isfc_matrix.T) / 2 if
symmetric else isfc_matrix for
isfc_matrix in isfcs])
# Get ISCs back into correct shape after masking out NaNs
isfcs_all = np.full((n_voxels, t_n_voxels, isfcs.shape[2]), np.nan)
isfcs_all[np.ix_(np.where(mask)[0], np.where(targets_mask)[0])] = isfcs
isfcs = np.moveaxis(isfcs_all, 2, 0)
# Summarize results (if requested)
if summary_statistic:
isfcs = compute_summary_statistic(isfcs,
summary_statistic=summary_statistic,
axis=0)
# Throw away first dimension if singleton
if isfcs.shape[0] == 1:
isfcs = isfcs[0]
# Optionally squareform to vectorize ISFC matrices (only if symmetric)
if vectorize_isfcs and symmetric:
isfcs, iscs = squareform_isfc(isfcs)
return isfcs, iscs
else:
return isfcs
def _check_isc_input(iscs, pairwise=False):
"""Checks ISC inputs for statistical tests
Input ISCs should be n_subjects (leave-one-out approach) or
n_pairs (pairwise approach) by n_voxels or n_ROIs array or a 1D
array (or list) of ISC values for a single voxel or ROI. This
function is only intended to be used internally by other
functions in this module (e.g., bootstrap_isc, permutation_isc).
Parameters
----------
iscs : ndarray or list
ISC values
Returns
-------
iscs : ndarray
Array of ISC values
n_subjects : int
Number of subjects
n_voxels : int
Number of voxels (or ROIs)
"""
# Standardize structure of input data
if type(iscs) == list:
iscs = np.array(iscs)[:, np.newaxis]
elif isinstance(iscs, np.ndarray):
if iscs.ndim == 1:
iscs = iscs[:, np.newaxis]
# Check if incoming pairwise matrix is vectorized triangle
if pairwise:
try:
test_square = squareform(iscs[:, 0], force='tomatrix')
n_subjects = test_square.shape[0]
except ValueError:
raise ValueError("For pairwise input, ISCs must be the "
"vectorized triangle of a square matrix.")
elif not pairwise:
n_subjects = iscs.shape[0]
# Infer subjects, voxels and print for user to check
n_voxels = iscs.shape[1]
logger.info("Assuming {0} subjects with and {1} "
"voxel(s) or ROI(s) in bootstrap ISC test.".format(n_subjects,
n_voxels))
return iscs, n_subjects, n_voxels
def _check_targets_input(targets, data):
"""Checks ISFC targets input array
For ISFC analysis, targets input array should either be a list
of n_TRs by n_targets arrays (where each array corresponds to
a subject), or an n_TRs by n_targets by n_subjects ndarray. This
function also checks the shape of the targets array against the
input data array.
Parameters
----------
data : list or ndarray (n_TRs x n_voxels x n_subjects)
fMRI data for which to compute ISFC
targets : list or ndarray (n_TRs x n_voxels x n_subjects)
fMRI data to use as targets for ISFC
Returns
-------
targets : ndarray (n_TRs x n_voxels x n_subjects)
ISFC targets with standadized structure
n_TRs : int
Number of time points (TRs) for targets array
n_voxels : int
Number of voxels (or ROIs) for targets array
n_subjects : int
Number of subjects for targets array
symmetric : bool
Indicator for symmetric vs. asymmetric
"""
if isinstance(targets, np.ndarray) or isinstance(targets, list):
targets, n_TRs, n_voxels, n_subjects = (
_check_timeseries_input(targets))
if data.shape[0] != n_TRs:
raise ValueError("Targets array must have same number of "
"TRs as input data")
if data.shape[2] != n_subjects:
raise ValueError("Targets array must have same number of "
"subjects as input data")
symmetric = False
else:
targets = data
n_TRs, n_voxels, n_subjects = data.shape
symmetric = True
return targets, n_TRs, n_voxels, n_subjects, symmetric
def compute_summary_statistic(iscs, summary_statistic='mean', axis=None):
"""Computes summary statistics for ISCs
Computes either the 'mean' or 'median' across a set of ISCs. In the
case of the mean, ISC values are first Fisher Z transformed (arctanh),
averaged, then inverse Fisher Z transformed (tanh).
The implementation is based on the work in [SilverDunlap1987]_.
.. [SilverDunlap1987] "Averaging correlation coefficients: should
Fisher's z transformation be used?", N. C. Silver, W. P. Dunlap, 1987,
Journal of Applied Psychology, 72, 146-148.
https://doi.org/10.1037/0021-9010.72.1.146
Parameters
----------
iscs : list or ndarray
ISC values
summary_statistic : str, default: 'mean'
Summary statistic, 'mean' or 'median'
axis : None or int or tuple of ints, optional
Axis or axes along which the means are computed. The default is to
compute the mean of the flattened array.
Returns
-------
statistic : float or ndarray
Summary statistic of ISC values
"""
if summary_statistic not in ('mean', 'median'):
raise ValueError("Summary statistic must be 'mean' or 'median'")
# Compute summary statistic
if summary_statistic == 'mean':
statistic = np.tanh(np.nanmean(np.arctanh(iscs), axis=axis))
elif summary_statistic == 'median':
statistic = np.nanmedian(iscs, axis=axis)
return statistic
def squareform_isfc(isfcs, iscs=None):
"""Converts square ISFCs to condensed ISFCs (and ISCs), and vice-versa
If input is a 2- or 3-dimensional array of square ISFC matrices, converts
this to the condensed off-diagonal ISFC values (i.e., the vectorized
triangle) and the diagonal ISC values. In this case, input must be a
single array of shape either n_voxels x n_voxels or n_subjects (or
n_pairs) x n_voxels x n_voxels. The condensed ISFC values are vectorized
according to scipy.spatial.distance.squareform, yielding n_voxels *
(n_voxels - 1) / 2 values comprising every voxel pair. Alternatively, if
input is an array of condensed off-diagonal ISFC values and an array of
diagonal ISC values, the square (redundant) ISFC values are returned.
This function mimics scipy.spatial.distance.squareform, but is intended
to retain the diagonal ISC values.
Parameters
----------
isfcs : ndarray
Either condensed or redundant ISFC values
iscs: ndarray, optional
Diagonal ISC values, required when input is condensed
Returns
-------
isfcs : ndarray or tuple of ndarrays
If condensed ISFCs are passed, a single redundant ISFC array is
returned; if redundant ISFCs are passed, both a condensed off-
diagonal ISFC array and the diagonal ISC values are returned
"""
# Check if incoming ISFCs are square (redundant)
if not type(iscs) == np.ndarray and isfcs.shape[-2] == isfcs.shape[-1]:
if isfcs.ndim == 2:
isfcs = isfcs[np.newaxis, ...]
if isfcs.ndim == 3:
iscs = np.diagonal(isfcs, axis1=1, axis2=2)
isfcs = np.vstack([squareform(isfc, checks=False)[np.newaxis, :]
for isfc in isfcs])
else:
raise ValueError("Square (redundant) ISFCs must be square "
"with multiple subjects or pairs of subjects "
"indexed by the first dimension")
if isfcs.shape[0] == iscs.shape[0] == 1:
isfcs, iscs = isfcs[0], iscs[0]
return isfcs, iscs
# Otherwise, convert from condensed to redundant
else:
if isfcs.ndim == iscs.ndim == 1:
isfcs, iscs = isfcs[np.newaxis, :], iscs[np.newaxis, :]
isfcs_stack = []
for isfc, isc in zip(isfcs, iscs):
isfc_sq = squareform(isfc, checks=False)
np.fill_diagonal(isfc_sq, isc)
isfcs_stack.append(isfc_sq[np.newaxis, ...])
isfcs = np.vstack(isfcs_stack)
if isfcs.shape[0] == 1:
isfcs = isfcs[0]
return isfcs
def _threshold_nans(data, tolerate_nans):
"""Thresholds data based on proportion of subjects with NaNs
Takes in data and a threshold value (float between 0.0 and 1.0) determining
the permissible proportion of subjects with non-NaN values. For example, if
threshold=.8, any voxel where >= 80% of subjects have non-NaN values will
be left unchanged, while any voxel with < 80% non-NaN values will be
assigned all NaN values and included in the nan_mask output. Note that the
output data has not been masked and will be same shape as the input data,
but may have a different number of NaNs based on the threshold.
Parameters
----------
data : ndarray (n_TRs x n_voxels x n_subjects)
fMRI time series data
tolerate_nans : bool or float (0.0 <= threshold <= 1.0)
Proportion of subjects with non-NaN values required to keep voxel
Returns
-------
data : ndarray (n_TRs x n_voxels x n_subjects)
fMRI time series data with adjusted NaNs
nan_mask : ndarray (n_voxels,)
Boolean mask array of voxels with too many NaNs based on threshold
"""
nans = np.all(np.any(np.isnan(data), axis=0), axis=1)
# Check tolerate_nans input and use either mean/nanmean and exclude voxels
if tolerate_nans is True:
logger.info("ISC computation will tolerate all NaNs when averaging")
elif type(tolerate_nans) is float:
if not 0.0 <= tolerate_nans <= 1.0:
raise ValueError("If threshold to tolerate NaNs is a float, "
"it must be between 0.0 and 1.0; got {0}".format(
tolerate_nans))
nans += ~(np.sum(~np.any(np.isnan(data), axis=0), axis=1) >=
data.shape[-1] * tolerate_nans)
logger.info("ISC computation will tolerate voxels with at least "
"{0} non-NaN values: {1} voxels do not meet "
"threshold".format(tolerate_nans,
np.sum(nans)))
else:
logger.info("ISC computation will not tolerate NaNs when averaging")
mask = ~nans
data = data[:, mask, :]
return data, mask
def bootstrap_isc(iscs, pairwise=False, summary_statistic='median',
n_bootstraps=1000, ci_percentile=95, side='right',
random_state=None):
"""One-sample group-level bootstrap hypothesis test for ISCs
For ISCs from one more voxels or ROIs, resample subjects with replacement
to construct a bootstrap distribution. Input is a list or ndarray of
ISCs for a single voxel/ROI, or an ISCs-by-voxels ndarray. ISC values
should be either N ISC values for N subjects in the leave-one-out appraoch
(pairwise=False), N(N-1)/2 ISC values for N subjects in the pairwise
approach (pairwise=True). In the pairwise approach, ISC values should
correspond to the vectorized upper triangle of a square correlation matrix
(see scipy.stats.distance.squareform). Shifts bootstrap distribution by
actual summary statistic (effectively to zero) for null hypothesis test
(Hall & Wilson, 1991). Uses subject-wise (not pair-wise) resampling in the
pairwise approach. Returns the observed ISC, the confidence interval, and
a p-value for the bootstrap hypothesis test, as well as the bootstrap
distribution of summary statistics. The p-value corresponds to either a
'two-sided', 'left'-, or 'right'-sided (default) test, as specified by
side. According to Chen et al., 2016, this is the preferred nonparametric
approach for controlling false positive rates (FPRs) for one-sample tests
in the pairwise approach. Note that the bootstrap hypothesis test may not
strictly control FPRs in the leave-one-out approach.
The implementation is based on the work in [Chen2016]_ and
[HallWilson1991]_.
.. [HallWilson1991] "Two guidelines for bootstrap hypothesis testing.",
P. Hall, S. R., Wilson, 1991, Biometrics, 757-762.
https://doi.org/10.2307/2532163
Parameters
----------
iscs : list or ndarray, ISCs by voxels array
ISC values for one or more voxels
pairwise : bool, default: False
Indicator of pairwise or leave-one-out, should match ISCs structure
summary_statistic : str, default: 'median'
Summary statistic, either 'median' (default) or 'mean'
n_bootstraps : int, default: 1000
Number of bootstrap samples (subject-level with replacement)
ci_percentile : int, default: 95
Percentile for computing confidence intervals
side : str
Perform one-sided ('left' or 'right') or 'two-sided' test
random_state = int or None, default: None
Initial random seed
Returns
-------
observed : float, median (or mean) ISC value
Summary statistic for actual ISCs
ci : tuple, bootstrap confidence intervals
Confidence intervals generated from bootstrap distribution
p : float, p-value
p-value based on bootstrap hypothesis test
distribution : ndarray, n_bootstraps by voxels
Bootstrap distribution
"""
# Standardize structure of input data
iscs, n_subjects, n_voxels = _check_isc_input(iscs, pairwise=pairwise)
# Check for valid summary statistic
if summary_statistic not in ('mean', 'median'):
raise ValueError("Summary statistic must be 'mean' or 'median'")
# Compute summary statistic for observed ISCs
observed = compute_summary_statistic(iscs,
summary_statistic=summary_statistic,
axis=0)
# Set up an empty list to build our bootstrap distribution
distribution = []
# Loop through n bootstrap iterations and populate distribution
for i in np.arange(n_bootstraps):
# Random seed to be deterministically re-randomized at each iteration
if isinstance(random_state, np.random.RandomState):
prng = random_state
else:
prng = np.random.RandomState(random_state)
# Randomly sample subject IDs with replacement
subject_sample = sorted(prng.choice(np.arange(n_subjects),
size=n_subjects))
# Squareform and shuffle rows/columns of pairwise ISC matrix to
# to retain correlation structure among ISCs, then get triangle
if pairwise:
# Loop through voxels
isc_sample = []
for voxel_iscs in iscs.T:
# Square the triangle and fill diagonal
voxel_iscs = squareform(voxel_iscs, force='tomatrix')
np.fill_diagonal(voxel_iscs, 1)
# Shuffle square correlation matrix and get triangle
voxel_sample = voxel_iscs[subject_sample, :][:, subject_sample]
voxel_sample = squareform(voxel_sample, checks=False)
# Censor off-diagonal 1s for same-subject pairs
voxel_sample[voxel_sample == 1.] = np.NaN
isc_sample.append(voxel_sample)
isc_sample = np.column_stack(isc_sample)
# Get simple bootstrap sample if not pairwise
elif not pairwise:
isc_sample = iscs[subject_sample, :]
# Compute summary statistic for bootstrap ISCs per voxel
# (alternatively could construct distribution for all voxels
# then compute statistics, but larger memory footprint)
distribution.append(compute_summary_statistic(
isc_sample,
summary_statistic=summary_statistic,
axis=0))
# Update random state for next iteration
random_state = np.random.RandomState(prng.randint(0, MAX_RANDOM_SEED))
# Convert distribution to numpy array
distribution = np.array(distribution)
# Compute CIs of median from bootstrap distribution (default: 95%)
ci = (np.percentile(distribution, (100 - ci_percentile)/2, axis=0),
np.percentile(distribution, ci_percentile + (100 - ci_percentile)/2,
axis=0))
# Shift bootstrap distribution to 0 for hypothesis test
shifted = distribution - observed
# Get p-value for actual median from shifted distribution
p = p_from_null(observed, shifted,
side=side, exact=False,
axis=0)
return observed, ci, p, distribution
def _check_group_assignment(group_assignment, n_subjects):
if type(group_assignment) == list:
pass
elif type(group_assignment) == np.ndarray:
group_assignment = group_assignment.tolist()
else:
logger.info("No group assignment provided, "
"performing one-sample test.")
if group_assignment and len(group_assignment) != n_subjects:
raise ValueError("Group assignments ({0}) "
"do not match number of subjects ({1})!".format(
len(group_assignment), n_subjects))
return group_assignment
def _get_group_parameters(group_assignment, n_subjects, pairwise=False):
# Set up dictionary to contain group info
group_parameters = {'group_assignment': group_assignment,
'n_subjects': n_subjects,
'group_labels': None, 'groups': None,
'sorter': None, 'unsorter': None,
'group_matrix': None, 'group_selector': None}
# Set up group selectors for two-group scenario
if group_assignment and len(np.unique(group_assignment)) == 2:
group_parameters['n_groups'] = 2
# Get group labels and counts
group_labels = np.unique(group_assignment)
groups = {group_labels[0]: group_assignment.count(group_labels[0]),
group_labels[1]: group_assignment.count(group_labels[1])}
# For two-sample pairwise approach set up selector from matrix
if pairwise:
# Sort the group_assignment variable if it came in shuffled
# so it's easier to build group assignment matrix
sorter = np.array(group_assignment).argsort()
unsorter = np.array(group_assignment).argsort().argsort()
# Populate a matrix with group assignments
upper_left = np.full((groups[group_labels[0]],
groups[group_labels[0]]),
group_labels[0])
upper_right = np.full((groups[group_labels[0]],
groups[group_labels[1]]),
np.nan)
lower_left = np.full((groups[group_labels[1]],
groups[group_labels[0]]),
np.nan)
lower_right = np.full((groups[group_labels[1]],
groups[group_labels[1]]),
group_labels[1])
group_matrix = np.vstack((np.hstack((upper_left, upper_right)),
np.hstack((lower_left, lower_right))))
np.fill_diagonal(group_matrix, np.nan)
group_parameters['group_matrix'] = group_matrix
# Unsort matrix and squareform to create selector
group_parameters['group_selector'] = squareform(
group_matrix[unsorter, :][:, unsorter],
checks=False)
group_parameters['sorter'] = sorter
group_parameters['unsorter'] = unsorter
# If leave-one-out approach, just user group assignment as selector
else:
group_parameters['group_selector'] = group_assignment
# Save these parameters for later
group_parameters['groups'] = groups
group_parameters['group_labels'] = group_labels
# Manage one-sample and incorrect group assignments
elif not group_assignment or len(np.unique(group_assignment)) == 1:
group_parameters['n_groups'] = 1
# If pairwise initialize matrix of ones for sign-flipping
if pairwise:
group_parameters['group_matrix'] = np.ones((
group_parameters['n_subjects'],
group_parameters['n_subjects']))
elif len(np.unique(group_assignment)) > 2:
raise ValueError("This test is not valid for more than "
"2 groups! (got {0})".format(
len(np.unique(group_assignment))))
else:
raise ValueError("Invalid group assignments!")
return group_parameters
def _permute_one_sample_iscs(iscs, group_parameters, i, pairwise=False,
summary_statistic='median', group_matrix=None,
exact_permutations=None, prng=None):
"""Applies one-sample permutations to ISC data
Input ISCs should be n_subjects (leave-one-out approach) or
n_pairs (pairwise approach) by n_voxels or n_ROIs array.
This function is only intended to be used internally by the
permutation_isc function in this module.
Parameters
----------
iscs : ndarray or list
ISC values
group_parameters : dict
Dictionary of group parameters
i : int
Permutation iteration
pairwise : bool, default: False
Indicator of pairwise or leave-one-out, should match ISCs variable
summary_statistic : str, default: 'median'
Summary statistic, either 'median' (default) or 'mean'
exact_permutations : list
List of permutations
prng = None or np.random.RandomState, default: None
Initial random seed
Returns
-------
isc_sample : ndarray
Array of permuted ISC values
"""
# Randomized sign-flips
if exact_permutations:
sign_flipper = np.array(exact_permutations[i])
else:
sign_flipper = prng.choice([-1, 1],
size=group_parameters['n_subjects'],
replace=True)
# If pairwise, apply sign-flips by rows and columns
if pairwise:
matrix_flipped = (group_parameters['group_matrix'] * sign_flipper
* sign_flipper[
:, np.newaxis])
sign_flipper = squareform(matrix_flipped, checks=False)
# Apply flips along ISC axis (same across voxels)
isc_flipped = iscs * sign_flipper[:, np.newaxis]
# Get summary statistics on sign-flipped ISCs
isc_sample = compute_summary_statistic(
isc_flipped,
summary_statistic=summary_statistic,
axis=0)
return isc_sample
def _permute_two_sample_iscs(iscs, group_parameters, i, pairwise=False,
summary_statistic='median',
exact_permutations=None, prng=None):
"""Applies two-sample permutations to ISC data
Input ISCs should be n_subjects (leave-one-out approach) or
n_pairs (pairwise approach) by n_voxels or n_ROIs array.
This function is only intended to be used internally by the
permutation_isc function in this module.
Parameters
----------
iscs : ndarray or list
ISC values
group_parameters : dict
Dictionary of group parameters
i : int
Permutation iteration
pairwise : bool, default: False
Indicator of pairwise or leave-one-out, should match ISCs variable
summary_statistic : str, default: 'median'
Summary statistic, either 'median' (default) or 'mean'
exact_permutations : list
List of permutations
prng = None or np.random.RandomState, default: None
Initial random seed