[pval, t_orig, clust_info, seed_state, est_alpha]=clust_perm2(dataA,dataB,chan_hood,n_perm,fwer,tail,thresh_p,verblevel,seed_state,freq_domain)
clust_perm2-Independent samples cluster-based permutation test using the "cluster
mass" statistic and a null hypothesis of a mean of zero. This
function can handle multiple electrodes and time points/frequencies.
This test was originally proposed for MRI data by Bullmore et al.
(1999) and for EEG/MEG analysis by Maris & Oostenveld (2007).
Required Inputs:
dataA - 3D matrix of ERP data (Channel x Time x Participant)
dataB - 3D matrix of ERP data (Channel x Time x Participant)
chan_hood - 2D symmetric binary matrix that indicates which channels are
considered neighbors of other channels. E.g., if
chan_hood(2,10)=1, then Channel 2 and Channel 10 are
nieghbors. You can produce a chan_hood matrix using the
function spatial_neighbors.m.
Optional Inputs:
n_perm - Number of permutations {default=2000}. Manly (1997) suggests
using at least 1000 permutations for an alpha level of 0.05 and
at least 5000 permutations for an alpha level of 0.01.
fwer - Desired family-wise error rate (i.e., alpha level) {default=.05}
tail - [1 | 0 | -1] If tail=1, the alternative hypothesis is that the
mean of the data is greater than 0 (upper tailed test). If tail=0,
the alternative hypothesis is that the mean of the data is different
than 0 (two tailed test). If tail=-1, the alternative hypothesis
is that the mean of the data is less than 0 (lower tailed test).
{default: 0}
thresh_p - The test-wise p-value threshold for cluster inclusion. If
a channel/time-point has a t-score that corresponds to an
uncorrected p-value greater than thresh_p, it is assigned
a p-value of 1 and not considered for clustering. Note
that thresh_p automatically takes into account the tail of
the test (e.g., you will get positive and negative t-score
thresholds for a two-tailed test).
verblevel - An integer specifiying the amount of information you want
this function to provide about what it is doing during runtime.
Options are:
0 - quiet, only show errors, warnings, and EEGLAB reports
1 - stuff anyone should probably know
2 - stuff you should know the first time you start working
with a data set {default value}
3 - stuff that might help you debug (show all
reports)
seed_state - The initial state of the random number generating stream
(see MATLAB documentation for "randstream"). If you pass
a value from a previous run of this function, it should
reproduce the exact same values.
freq_domain - If 0, command line report will be given in temporal units
(e.g. time points). Otherwise, the report will be given
in frequency domain units (e.g., frequencies).
{default: 0}
Outputs:
pval - p-value at each time point and electrode (corrected for
multiple comparisons via the permutation test)
t_orig - t-score at each time point and electrode
clust_info - A struct variable containing information about the
clusters found. Depending on the tail of the test it will
be composed of all or half of the following fields:
pos_clust_pval: p-values of the positive clusters
pos_clust_mass: t-score mass of the positive clusters
pos_clust_ids: channel x time point matrix that
indicated which positive cluster each channel/time point
pair belongs to. E.g., if pos_clust_ids(1,2)=4, then
the second time point of the first channel is a
member of the fourth cluster. 0 indicates that the
channel/time point is not a member of any positive
cluster.
neg_clust_pval: p-values of the negative clusters
neg_clust_mass: t-score mass of the negative clusters
neg_clust_ids: channel x time point matrix that
indicated which negative cluster each channel/time point
pair belongs to. E.g., if neg_clust_ids(1,2)=4, then
the second time point of the first channel is a
member of the fourth cluster. 0 indicates that the
channel/time point is not a member of any negative
cluster.
seed_state - The initial state of the random number generating stream
(see MATLAB documentation for "randstream") used to
generate the permutations. You can use this to reproduce
the output of this function.
est_alpha - The estimated family-wise alpha level of the test. With
permutation tests, a finite number of p-values are possible.
This function tries to use an alpha level that is as close
as possible to the desired alpha level. However, if the
sample size is small, a very limited number of p-values are
possible and the desired family-wise alpha level may be
impossible to acheive.
Notes:
-The null hypothesis of this test is that the data in the two groups are
"exchangeable" (i.e., the data from a participant in Group A was just as
likely to have come from Group B). If the distributions of the two
populations you're comparing differ in *any* way (i.e., not just in
central tendency), the null hypothesis will be violated. For example, if
the ERPs of the two groups differ in their variability, but not their
average value (e.g., one set of ERPs is from a group of children and the
other adults), you'll be more likely than alpha to reject the null
hypothesis.
-Unlike a parametric test (e.g., an ANOVA), a discrete set of p-values
are possible (at most the number of possible permutations). Since the
number of possible permutations grows rapdily with the number of
participants, this is only an issue for small sample sizes (e.g., 3
participants in each group). When you have such a small sample size, the
limited number of p-values may make the test overly conservative (e.g.,
you might be forced to use an alpha level of .0286 since it is the biggest
possible alpha level less than .05).
-When doing two-tailed tests it is possible to get p-values greater than 1.
These can be treated equivalent to a p-value of 1. The reason for
this is that the p-values for positive and negative clusters are computed
from two different distributions of permuted values. For a non-cluster
based permutation test, there is only one distribution of permuted
values.
Author:
David Groppe
May, 2011
Kutaslab, San Diego
References:
Bullmore, E. T., Suckling, J., Overmeyer, S., Rabe-Hesketh, S., Taylor,
E., & Brammer, M. J. (1999). Global, voxel, and cluster tests, by theory
and permutation, for a difference between two groups of structural MR
images of the brain. IEEE Transactions on Medical Imaging, 18(1), 32-42.
doi:10.1109/42.750253
Maris, E., & Oostenveld, R. (2007). Nonparametric statistical testing of
EEG- and MEG-data. Journal of Neuroscience Methods, 164(1), 177-190.
doi:10.1016/j.jneumeth.2007.03.024
Manly, B.F.J. (1997) Randomization, bootstrap, and Monte Carlo methods in
biology. 2nd ed. Chapmn and Hall, London.0001 %[pval, t_orig, clust_info, seed_state, est_alpha]=clust_perm2(dataA,dataB,chan_hood,n_perm,fwer,tail,thresh_p,verblevel,seed_state,freq_domain) 0002 % 0003 % clust_perm2-Independent samples cluster-based permutation test using the "cluster 0004 % mass" statistic and a null hypothesis of a mean of zero. This 0005 % function can handle multiple electrodes and time points/frequencies. 0006 % This test was originally proposed for MRI data by Bullmore et al. 0007 % (1999) and for EEG/MEG analysis by Maris & Oostenveld (2007). 0008 % 0009 % Required Inputs: 0010 % dataA - 3D matrix of ERP data (Channel x Time x Participant) 0011 % dataB - 3D matrix of ERP data (Channel x Time x Participant) 0012 % chan_hood - 2D symmetric binary matrix that indicates which channels are 0013 % considered neighbors of other channels. E.g., if 0014 % chan_hood(2,10)=1, then Channel 2 and Channel 10 are 0015 % nieghbors. You can produce a chan_hood matrix using the 0016 % function spatial_neighbors.m. 0017 % 0018 % Optional Inputs: 0019 % n_perm - Number of permutations {default=2000}. Manly (1997) suggests 0020 % using at least 1000 permutations for an alpha level of 0.05 and 0021 % at least 5000 permutations for an alpha level of 0.01. 0022 % fwer - Desired family-wise error rate (i.e., alpha level) {default=.05} 0023 % tail - [1 | 0 | -1] If tail=1, the alternative hypothesis is that the 0024 % mean of the data is greater than 0 (upper tailed test). If tail=0, 0025 % the alternative hypothesis is that the mean of the data is different 0026 % than 0 (two tailed test). If tail=-1, the alternative hypothesis 0027 % is that the mean of the data is less than 0 (lower tailed test). 0028 % {default: 0} 0029 % thresh_p - The test-wise p-value threshold for cluster inclusion. If 0030 % a channel/time-point has a t-score that corresponds to an 0031 % uncorrected p-value greater than thresh_p, it is assigned 0032 % a p-value of 1 and not considered for clustering. Note 0033 % that thresh_p automatically takes into account the tail of 0034 % the test (e.g., you will get positive and negative t-score 0035 % thresholds for a two-tailed test). 0036 % verblevel - An integer specifiying the amount of information you want 0037 % this function to provide about what it is doing during runtime. 0038 % Options are: 0039 % 0 - quiet, only show errors, warnings, and EEGLAB reports 0040 % 1 - stuff anyone should probably know 0041 % 2 - stuff you should know the first time you start working 0042 % with a data set {default value} 0043 % 3 - stuff that might help you debug (show all 0044 % reports) 0045 % seed_state - The initial state of the random number generating stream 0046 % (see MATLAB documentation for "randstream"). If you pass 0047 % a value from a previous run of this function, it should 0048 % reproduce the exact same values. 0049 % freq_domain - If 0, command line report will be given in temporal units 0050 % (e.g. time points). Otherwise, the report will be given 0051 % in frequency domain units (e.g., frequencies). 0052 % {default: 0} 0053 % 0054 % Outputs: 0055 % pval - p-value at each time point and electrode (corrected for 0056 % multiple comparisons via the permutation test) 0057 % t_orig - t-score at each time point and electrode 0058 % clust_info - A struct variable containing information about the 0059 % clusters found. Depending on the tail of the test it will 0060 % be composed of all or half of the following fields: 0061 % pos_clust_pval: p-values of the positive clusters 0062 % pos_clust_mass: t-score mass of the positive clusters 0063 % pos_clust_ids: channel x time point matrix that 0064 % indicated which positive cluster each channel/time point 0065 % pair belongs to. E.g., if pos_clust_ids(1,2)=4, then 0066 % the second time point of the first channel is a 0067 % member of the fourth cluster. 0 indicates that the 0068 % channel/time point is not a member of any positive 0069 % cluster. 0070 % neg_clust_pval: p-values of the negative clusters 0071 % neg_clust_mass: t-score mass of the negative clusters 0072 % neg_clust_ids: channel x time point matrix that 0073 % indicated which negative cluster each channel/time point 0074 % pair belongs to. E.g., if neg_clust_ids(1,2)=4, then 0075 % the second time point of the first channel is a 0076 % member of the fourth cluster. 0 indicates that the 0077 % channel/time point is not a member of any negative 0078 % cluster. 0079 % seed_state - The initial state of the random number generating stream 0080 % (see MATLAB documentation for "randstream") used to 0081 % generate the permutations. You can use this to reproduce 0082 % the output of this function. 0083 % est_alpha - The estimated family-wise alpha level of the test. With 0084 % permutation tests, a finite number of p-values are possible. 0085 % This function tries to use an alpha level that is as close 0086 % as possible to the desired alpha level. However, if the 0087 % sample size is small, a very limited number of p-values are 0088 % possible and the desired family-wise alpha level may be 0089 % impossible to acheive. 0090 % 0091 % Notes: 0092 % -The null hypothesis of this test is that the data in the two groups are 0093 % "exchangeable" (i.e., the data from a participant in Group A was just as 0094 % likely to have come from Group B). If the distributions of the two 0095 % populations you're comparing differ in *any* way (i.e., not just in 0096 % central tendency), the null hypothesis will be violated. For example, if 0097 % the ERPs of the two groups differ in their variability, but not their 0098 % average value (e.g., one set of ERPs is from a group of children and the 0099 % other adults), you'll be more likely than alpha to reject the null 0100 % hypothesis. 0101 % 0102 % -Unlike a parametric test (e.g., an ANOVA), a discrete set of p-values 0103 % are possible (at most the number of possible permutations). Since the 0104 % number of possible permutations grows rapdily with the number of 0105 % participants, this is only an issue for small sample sizes (e.g., 3 0106 % participants in each group). When you have such a small sample size, the 0107 % limited number of p-values may make the test overly conservative (e.g., 0108 % you might be forced to use an alpha level of .0286 since it is the biggest 0109 % possible alpha level less than .05). 0110 % 0111 % -When doing two-tailed tests it is possible to get p-values greater than 1. 0112 % These can be treated equivalent to a p-value of 1. The reason for 0113 % this is that the p-values for positive and negative clusters are computed 0114 % from two different distributions of permuted values. For a non-cluster 0115 % based permutation test, there is only one distribution of permuted 0116 % values. 0117 % 0118 % Author: 0119 % David Groppe 0120 % May, 2011 0121 % Kutaslab, San Diego 0122 % 0123 % References: 0124 % Bullmore, E. T., Suckling, J., Overmeyer, S., Rabe-Hesketh, S., Taylor, 0125 % E., & Brammer, M. J. (1999). Global, voxel, and cluster tests, by theory 0126 % and permutation, for a difference between two groups of structural MR 0127 % images of the brain. IEEE Transactions on Medical Imaging, 18(1), 32-42. 0128 % doi:10.1109/42.750253 0129 % 0130 % Maris, E., & Oostenveld, R. (2007). Nonparametric statistical testing of 0131 % EEG- and MEG-data. Journal of Neuroscience Methods, 164(1), 177-190. 0132 % doi:10.1016/j.jneumeth.2007.03.024 0133 % 0134 % Manly, B.F.J. (1997) Randomization, bootstrap, and Monte Carlo methods in 0135 % biology. 2nd ed. Chapmn and Hall, London. 0136 0137 %%%%%%%%%%%%%%%% REVISION LOG %%%%%%%%%%%%%%%%% 0138 % 0139 % 12/11/2011-Now uses Amy Guthormsen's recursiveless find_clusters.m. 0140 % 0141 % 3/17/2013-Randomization is now compatible with Matlab v13. Thanks to 0142 % Aaron Newman for the fix. 0143 % 0144 0145 % 0146 %This code has an appropriate false positive rate when run on Gaussian 0147 %noise (26 channels, 25 time points, 16 subs, two-tailed test). Done with 0148 %sim_test_clust2.m 0149 % 0150 0151 %%%%%%%%%%%%%%%% FUTURE IMPROVEMENTS %%%%%%%%%%%%%%%%% 0152 % 0153 % -Test for equality of trials or variance? 0154 % -Add alternate test statistics? 0155 0156 function [pval, t_orig, clust_info, seed_state, est_alpha]=clust_perm2(dataA,dataB,chan_hood,n_perm,fwer,tail,thresh_p,verblevel,seed_state,freq_domain) 0157 0158 if nargin<2, 0159 error('You need to provide data for two groups of participants.'); 0160 end 0161 0162 if nargin<3, 0163 error('You need to provide a chan_hood matrix.'); 0164 end 0165 0166 if nargin<4, 0167 n_perm=2000; 0168 end 0169 0170 if nargin<5, 0171 fwer=.05; 0172 elseif (fwer>=1) || (fwer<=0) 0173 error('Argument ''fwer'' needs to be between 0 and 1.'); 0174 end 0175 0176 if fwer<=.01 && n_perm<5000, 0177 watchit(sprintf('You are probably using too few permutations for an alpha level of %f. Type ">>help clust_perm2" for more info.',fwer)); 0178 elseif fwer<=.05 && n_perm<1000, 0179 watchit(sprintf('You are probably using too few permutations for an alpha level of %f. Type ">>help clust_perm2" for more info.',fwer)); 0180 end 0181 0182 if nargin<6, 0183 tail=0; 0184 elseif (tail~=0) && (tail~=1) && (tail~=-1), 0185 error('Argument ''tail'' needs to be 0,1, or -1.'); 0186 end 0187 0188 if nargin<7, 0189 thresh_p=.05; 0190 elseif thresh_p<=0 || thresh_p>1, 0191 error('Argument thresh_p needs to take a value between 0 and 1'); 0192 end 0193 0194 if nargin<8, 0195 verblevel=2; 0196 end 0197 0198 %get random # generator state 0199 if verLessThan('matlab','8.1') 0200 defaultStream=RandStream.getDefaultStream; 0201 else 0202 defaultStream=RandStream.getGlobalStream; 0203 end 0204 if (nargin<9) || isempty(seed_state), 0205 %Store state of random number generator 0206 seed_state=defaultStream.State; 0207 else 0208 defaultStream.State=seed_state; %reset random number generator to saved state 0209 end 0210 0211 if (nargin<10), 0212 freq_domain=0; 0213 end 0214 0215 if length(size(dataA))~=3 || length(size(dataB))~=3 0216 error('dataA and dataB need to be three dimensional (chan x time x participant)') 0217 end 0218 [n_chan, n_tpt, n_subsA]=size(dataA); 0219 [n_chan2, n_tpt2, n_subsB]=size(dataB); 0220 0221 if verblevel, 0222 warning('off','all'); %for large # of subjects, nchoosek warns that its result is approximate 0223 n_psbl_prms=nchoosek(n_subsA+n_subsB,n_subsA); 0224 if n_psbl_prms<100, 0225 watchit(sprintf(['Due to the very limited number of participants in each group,' ... 0226 ' the total number of possible permutations is small.\nThus only a limited number of p-values (at most %d) are possible and the test might be overly conservative.'], ... 0227 n_psbl_prms)); 0228 end 0229 warning('on','all'); 0230 end 0231 0232 if n_chan~=n_chan2, 0233 error('The number of channels in Group 1 (dataA) and Group 2 (dataB) need to be equal.'); 0234 elseif n_tpt~=n_tpt2, 0235 error('The number of time points in Group 1 (dataA) and Group 2 (dataB) need to be equal.'); 0236 end 0237 %combine data 0238 total_subs=n_subsA+n_subsB; 0239 data=dataA; 0240 data(:,:,(n_subsA+1):total_subs)=dataB; 0241 0242 if verblevel~=0, 0243 fprintf('clust_perm2: Number of channels: %d\n',n_chan); 0244 if freq_domain, 0245 fprintf('clust_perm2: Number of frequencies: %d\n',n_tpt); 0246 else 0247 fprintf('clust_perm2: Number of time points: %d\n',n_tpt); 0248 end 0249 fprintf('clust_perm2: Total # of comparisons: %d\n',n_chan*n_tpt); 0250 fprintf('clust_perm2: Number of participants in Group A: %d\n',n_subsA); 0251 fprintf('clust_perm2: Number of participants in Group B: %d\n',n_subsB); 0252 fprintf('t-score degrees of freedom: %d\n',total_subs-2); 0253 end 0254 VerbReport(sprintf('Executing permutation test with %d permutations...',n_perm),2,verblevel); 0255 if (verblevel>=2), 0256 fprintf('Permutations completed: '); 0257 end 0258 0259 % Factors that are used to compute t-scores. Saves time to compute them 0260 % now rather than to compute them anew for each permutation. 0261 df=n_subsA+n_subsB-2; 0262 mult_fact=(n_subsA+n_subsB)/(n_subsA*n_subsB); 0263 if tail 0264 %one tailed test 0265 thresh_t=tinv(thresh_p,df); %note unless thresh_p is greater than .5, thresh_t will be negative 0266 else 0267 %two tailed test 0268 thresh_t=tinv(thresh_p/2,df); 0269 end 0270 0271 if tail==0 0272 mx_clust_massNEG=zeros(1,n_perm); 0273 mx_clust_massPOS=zeros(1,n_perm); 0274 else 0275 mx_clust_mass=zeros(1,n_perm); 0276 end 0277 for perm=1:n_perm 0278 if ~rem(perm,100) 0279 if (verblevel>=2), 0280 if ~rem(perm-100,1000) 0281 fprintf('%d',perm); 0282 else 0283 fprintf(', %d',perm); 0284 end 0285 if ~rem(perm,1000) 0286 fprintf('\n'); 0287 end 0288 end 0289 end 0290 %randomly assign participants to conditions 0291 r=randperm(total_subs); 0292 grp1=r(1:n_subsA); 0293 grp2=r((n_subsA+1):total_subs); 0294 0295 %compute t-scores 0296 all_t=tmax2(data,grp1,grp2,n_subsA,n_subsB,df,mult_fact); 0297 0298 %form t-scores into clusters 0299 if tail==0, 0300 %two-tailed test 0301 0302 %positive clusters 0303 [clust_ids, n_clust]=find_clusters(all_t,-thresh_t,chan_hood,1); %note, thresh_t should be negative by default 0304 mx_clust_massPOS(perm)=find_mx_mass(clust_ids,all_t,n_clust,1); 0305 0306 %negative clusters 0307 [clust_ids, n_clust]=find_clusters(all_t,thresh_t,chan_hood,-1); 0308 mx_clust_massNEG(perm)=find_mx_mass(clust_ids,all_t,n_clust,-1); 0309 0310 elseif tail>0 0311 %upper tailed test 0312 [clust_ids, n_clust]=find_clusters(all_t,-thresh_t,chan_hood,1); %note, thresh_t should be negative by default 0313 mx_clust_mass(perm)=find_mx_mass(clust_ids,all_t,n_clust,1); 0314 else 0315 %lower tailed test 0316 [clust_ids, n_clust]=find_clusters(all_t,thresh_t,chan_hood,-1); %note, thresh_t should be negative by default 0317 mx_clust_mass(perm)=find_mx_mass(clust_ids,all_t,n_clust,-1); 0318 end 0319 end 0320 0321 %End of permutations, print carriage return if it hasn't already been done 0322 %(i.e., perm is NOT a multiple of 1000) 0323 if (verblevel>=2) && rem(perm,1000) 0324 fprintf('\n'); 0325 end 0326 0327 %Compute critical t's 0328 if tail==0, 0329 %two-tailed, test statistic is biggest absolute value of all t's 0330 mx_clust_mass=abs([mx_clust_massPOS mx_clust_massNEG]); 0331 tmx_ptile(2)=prctile(mx_clust_mass,100-100*fwer); 0332 tmx_ptile(1)=-tmx_ptile(2); 0333 est_alpha=mean(mx_clust_mass>=tmx_ptile(2)); 0334 elseif tail==1, 0335 %upper tailed 0336 tmx_ptile=prctile(mx_clust_mass,100-100*fwer); 0337 est_alpha=mean(mx_clust_mass>=tmx_ptile); 0338 else 0339 %tail=-1, lower tailed 0340 tmx_ptile=prctile(mx_clust_mass,fwer*100); 0341 est_alpha=mean(mx_clust_mass<=tmx_ptile); 0342 end 0343 if verblevel~=0, 0344 fprintf('Desired family-wise alpha level: %f\n',fwer); 0345 fprintf('Estimated actual family-wise alpha level: %f\n',est_alpha); 0346 end 0347 0348 %Compute t-scores of actual observations 0349 t_orig=tmax2(data,1:n_subsA,(n_subsA+1):total_subs,n_subsA,n_subsB,df,mult_fact); 0350 0351 %compute p-values 0352 pval=ones(n_chan,n_tpt); 0353 if tail==0, 0354 %two-tailed test 0355 pval=pval*2; %default p-value for channel/time point pairs is 2 0356 0357 %positive clusters 0358 [clust_ids, n_clust]=find_clusters(t_orig,-thresh_t,chan_hood,1); %note thresh_t is negative by default 0359 clust_info.pos_clust_pval=ones(1,n_clust); 0360 clust_info.pos_clust_mass=zeros(1,n_clust); 0361 clust_info.pos_clust_ids=clust_ids; 0362 for a=1:n_clust, 0363 use_ids=find(clust_ids==a); 0364 clust_mass=sum(t_orig(use_ids)); 0365 clust_p=mean(mx_clust_massPOS>=clust_mass)*2; %multiply by two to correct for doing an upper AND lower tailed test 0366 pval(use_ids)=clust_p; 0367 clust_info.pos_clust_pval(a)=clust_p; 0368 clust_info.pos_clust_mass(a)=clust_mass; 0369 end 0370 0371 %negative clusters 0372 [clust_ids, n_clust]=find_clusters(t_orig,thresh_t,chan_hood,-1); %note thresh_t is negative by default 0373 clust_info.neg_clust_pval=ones(1,n_clust); 0374 clust_info.neg_clust_mass=zeros(1,n_clust); 0375 clust_info.neg_clust_ids=clust_ids; 0376 for a=1:n_clust, 0377 use_ids=find(clust_ids==a); 0378 clust_mass=sum(t_orig(use_ids)); 0379 clust_p=mean(mx_clust_massNEG<=clust_mass)*2; %multiply by two to correct for doing an upper AND lower tailed test 0380 pval(use_ids)=clust_p; 0381 clust_info.neg_clust_pval(a)=clust_p; 0382 clust_info.neg_clust_mass(a)=clust_mass; 0383 end 0384 elseif tail>0 0385 %positive clusters 0386 [clust_ids, n_clust]=find_clusters(t_orig,-thresh_t,chan_hood,1); %note thresh_t is negative by default 0387 clust_info.pos_clust_pval=ones(1,n_clust); 0388 clust_info.pos_clust_mass=zeros(1,n_clust); 0389 clust_info.pos_clust_ids=clust_ids; 0390 for a=1:n_clust, 0391 use_ids=find(clust_ids==a); 0392 clust_mass=sum(t_orig(use_ids)); 0393 clust_p=mean(mx_clust_mass>=clust_mass); 0394 pval(use_ids)=clust_p; 0395 clust_info.pos_clust_pval(a)=clust_p; 0396 clust_info.pos_clust_mass(a)=clust_mass; 0397 end 0398 else 0399 %negative clusters 0400 [clust_ids, n_clust]=find_clusters(t_orig,thresh_t,chan_hood,-1); %note thresh_t is negative by default 0401 clust_info.neg_clust_pval=ones(1,n_clust); 0402 clust_info.neg_clust_mass=zeros(1,n_clust); 0403 clust_info.neg_clust_ids=clust_ids; 0404 for a=1:n_clust, 0405 use_ids=find(clust_ids==a); 0406 clust_mass=sum(t_orig(use_ids)); 0407 clust_p=mean(mx_clust_mass<=clust_mass); 0408 pval(use_ids)=clust_p; 0409 clust_info.neg_clust_pval(a)=clust_p; 0410 clust_info.neg_clust_mass(a)=clust_mass; 0411 end 0412 end 0413 0414 0415 %%% End of Main Function %%% 0416 0417 function all_t=tmax2(dat,grp1,grp2,n_subsA,n_subsB,df,mult_fact) 0418 % might make this faster by moving it into the code 0419 0420 x1=dat(:,:,grp1); 0421 x2=dat(:,:,grp2); 0422 0423 sm1=sum(x1,3); 0424 mn1=sm1/n_subsA; 0425 ss1=sum(x1.^2,3)-(sm1.^2)/n_subsA; 0426 0427 sm2=sum(x2,3); 0428 mn2=sm2/n_subsB; 0429 ss2=sum(x2.^2,3)-(sm2.^2)/n_subsB; 0430 0431 pooled_var=(ss1+ss2)/df; 0432 stder=sqrt(pooled_var*mult_fact); 0433 0434 all_t=(mn1-mn2)./stder; 0435 0436 0437 function mx_clust_mass=find_mx_mass(clust_ids,data_t,n_clust,tail) 0438 0439 mx_clust_mass=0; 0440 if tail<0 0441 %looking for most negative cluster mass 0442 for z=1:n_clust, 0443 use_ids=(clust_ids==z); 0444 use_mass=sum(data_t(use_ids)); 0445 if use_mass<mx_clust_mass, 0446 mx_clust_mass=use_mass; 0447 end 0448 end 0449 elseif tail>0, 0450 %looking for most positive cluster mass 0451 for z=1:n_clust, 0452 use_ids=(clust_ids==z); 0453 use_mass=sum(data_t(use_ids)); 0454 if use_mass>mx_clust_mass, 0455 mx_clust_mass=use_mass; 0456 end 0457 end 0458 end 0459 0460 0461