clustGND() - Tests the null hypothesis that the grand average voltage
of a bin is mu or that the grand average within-subject
difference between two bins is mu using a cluster-based
permutation test and the "cluster mass" statistic
(Bullmore et al., 1999; Maris & Oostenveld, 2007). Note, mu
is assumed to be 0 by default. This function requires
individual subject ERPs to be stored in a "GND" structure and
outputs the test results in a number of graphical and text
formats. For analogous between-subject comparisons use the
function clustGRP.m.
Usage:
>> [GND, prm_pval, data_t]=clustGND(GND_or_fname,bin,varargin);
Required Inputs:
GND_or_fname - A GND structure variable or the filename of a
GND structure that has been saved to disk. To
create a GND variable from Kutaslab ERP files (e.g.,
*.mas files) use avgs2GND.m. To do the same from
EEGLAB *.set files use sets2GND.m. See Mass
Univariate ERP Toolbox documentation for detailed
information about the format of a GND variable. If you
specify a filename be sure to include the file's path,
unless the file is in the current working directory.
bin - [integer] The bin to contrast against a voltage of mu
Use the function headinfo.m to see what bins are stored
in a GND variable. Use the function bin_dif.m to create
a difference wave between two bins whose significance
you can test with this function.
Optional Inputs:
tail - [1,0, or -1] An integer specifying the tail of the
hypothesis test. "1" means upper-tailed (i.e., alternative
hypothesis is that the ERP/difference wave is greater
than the mean of the null hypothesis). "0" means two-
tailed (i.e., alternative hypothesis is that the
ERP/difference wave is not equal to the mean of the null
hypothesis). "-1" means lower-tailed (i.e., alternative
hypothesis is that the ERP/difference wave is less than
the mean of the null hypothesis). {default: 0}
alpha - A number between 0 and 1 specifying the family-wise
error rate (FWER) of the test. Note, cluster-based tests
only weak control of FWER. {default: 0.05}
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).
chan_hood - A scalar or a 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.
If a scalar is provided, then all electrodes within that
distance of a particular electrode are considered
neighbors. Note, EEGLAB's electrode coordinates assume
the head has a radius of 1. See the help documentation
of the function spatial_neighbors to see how you could
convert this distance threshold to centimeters.
{default: 0.61}
head_radius - The radius of the head in whatever units the Cartesian
coordinates in GRP.chanlocs are in. This is used to
convert scalar values of chan_hood into centimeters.
{default: []}
n_perm - The number of permutations to use in the test. As this
value increases, the test takes longer to compute and
the results become more reliable. 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. {default: 2500}
time_wind - Pair of time values specifying the beginning and
end of a time window in ms (e.g., [160 180]). Every
single time point in the time window will be individually
tested (i.e., maximal temporal resolution) if mean_wind
option is NOT used. Note, boundaries of time window
may not exactly correspond to desired time window
boundaries because of temporal digitization (i.e., you
only have samples every so many ms). {default: 0 ms to
the end of the epoch}
mean_wind - ['yes' or 'no'] If 'yes', the permutation test will be
performed on the mean amplitude within the time window
specified by time_wind. This sacrifices temporal
resolution to increase test power by reducing the number
of comparisons. If 'no', every single time point within
time_wind's time windows will be tested individually.
{default: 'no'}
null_mean - [number] The mean of the null hypothesis (i.e., mu) in
units of microvolts. {default: 0}
exclude_chans - A cell array of channel labels to exclude from the
permutation test (e.g., {'A2','lle','rhe'}). This option
sacrifices spatial resolution to increase test power by
reducing the number of comparisons. Use headinfo.m to see
the channel labels stored in the GND variable. You cannot
use both this option and 'include_chans' (below).{default:
not used, all channels included in test}
include_chans - A cell array of channel labels to use in the permutation
test (e.g., {'A2','lle','rhe'}). All other channels will
be ignored. This option sacrifices spatial resolution to
increase test power by reducing the number of comparisons.
Use headinfo.m to see the channel labels stored in the GND
variable. You cannot use both this option and
'exclude_chans' (above). {default: not used, all channels
included in 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)
plot_raster - ['yes' or 'no'] If 'yes', a two-dimensional (time x channel)
binary "raster" diagram is created to illustrate the
results of the permutation test. Significant negative and
positive deviations from the null hypothesis are shown
as black and white rectangles respectively. Non-
significant comparisons are shown as gray rectangles.
Clicking on the rectangles will show you the
corresponding time and channel label for that
rectangle. This figure can be reproduced with the
function sig_raster.m. {default: 'yes'}
plot_mn_topo - ['yes' or 'no'] If 'yes', the topography of the mean
voltages/effects in the time window is produced. More
specifically, two figures are produced: one showing the
topographies in uV the other in t-scores. Significant/
nonsignificant comparisons are shown as white/black
electrodes. Clicking on electrodes will show the
electrode's name. This figure can be reproduced with
the function sig_topo.m. This option has NO effect if
mean_wind option is set to 'no'. {default: 'yes'}
output_file - A string indicating the name of a space delimited text
file to produce containing the p-values of all comparisons
and the details of the test (e.g., number of permutations,
family-wise alpha level, etc...). If mean_wind option is
set to 'yes,' t-scores of each comparison are also
included since you cannot derive them from the t-scores
at each time point/electrode in a simple way. When
importing this file into a spreadsheet be sure NOT to count
consecutive spaces as multiple delimiters. {default:
none}
save_GND - ['yes' or 'no'] If 'yes', the GND variable will be
saved to disk after the permutation test is completed
and added to it. User will first be prompted to verify
file name and path. {default: 'yes'}
reproduce_test- [integer] The number of the permutation test stored in
the GND variable to reproduce. For example, if
'reproduce_test' equals 2, the second t-test
stored in the GND variable (i.e., GND.t_tests(2)) will
be reproduced. Reproduction is accomplished by setting
the random number generator used in the permutation test
to the same initial state it was in when the permutation
test was first applied.
Outputs:
GND - GND structure variable. This is the same as
the input GND variable with one addition: the
field GND.t_tests will contain the results of the
permutation test and the test parameters.
prm_pval - A two-dimensional matrix (channel x time) of the
p-values of each comparison. These p-values are
corrected for multiple comparisons by the permutation
test.
data_t - A two-dimensional matrix (channel x time) of the
t-scores of each comparison.
Note also that a great deal of information about the test is displayed
in the MATLAB command window. You can easiy record of all this
information to a text file using the MATLAB command "diary."
Global Variables:
VERBLEVEL - Mass Univariate ERP Toolbox level of verbosity (i.e., tells
functions how much to report about what they're doing during
runtime) set by the optional function argument 'verblevel'
Notes:
-To add a difference wave to a GND variable, use the function "bin_dif".
-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 issue for small sample sizes (e.g., 6
participants). When you have such a small sample size, the
limited number of p-values may make the test less 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).
-For two-tailed tests clustGND.m effectively performs an upper-tailed
test and a lower-tailed test and multiplies the resulting p-values by 2.
In other words, it does two tests and then applies Bonferroni correction.
This is faster than performing a "proper" two-tailed permutation test
since you only have to form clusters for one polarity (e.g., if you form
negative clusters for each permutation, the negative of that distribution
will be valid for computing p-values for postive clusters) but can
result in p-values greater than 1.
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 % clustGND() - Tests the null hypothesis that the grand average voltage 0002 % of a bin is mu or that the grand average within-subject 0003 % difference between two bins is mu using a cluster-based 0004 % permutation test and the "cluster mass" statistic 0005 % (Bullmore et al., 1999; Maris & Oostenveld, 2007). Note, mu 0006 % is assumed to be 0 by default. This function requires 0007 % individual subject ERPs to be stored in a "GND" structure and 0008 % outputs the test results in a number of graphical and text 0009 % formats. For analogous between-subject comparisons use the 0010 % function clustGRP.m. 0011 % 0012 % 0013 % Usage: 0014 % >> [GND, prm_pval, data_t]=clustGND(GND_or_fname,bin,varargin); 0015 % 0016 % Required Inputs: 0017 % GND_or_fname - A GND structure variable or the filename of a 0018 % GND structure that has been saved to disk. To 0019 % create a GND variable from Kutaslab ERP files (e.g., 0020 % *.mas files) use avgs2GND.m. To do the same from 0021 % EEGLAB *.set files use sets2GND.m. See Mass 0022 % Univariate ERP Toolbox documentation for detailed 0023 % information about the format of a GND variable. If you 0024 % specify a filename be sure to include the file's path, 0025 % unless the file is in the current working directory. 0026 % bin - [integer] The bin to contrast against a voltage of mu 0027 % Use the function headinfo.m to see what bins are stored 0028 % in a GND variable. Use the function bin_dif.m to create 0029 % a difference wave between two bins whose significance 0030 % you can test with this function. 0031 % 0032 % Optional Inputs: 0033 % tail - [1,0, or -1] An integer specifying the tail of the 0034 % hypothesis test. "1" means upper-tailed (i.e., alternative 0035 % hypothesis is that the ERP/difference wave is greater 0036 % than the mean of the null hypothesis). "0" means two- 0037 % tailed (i.e., alternative hypothesis is that the 0038 % ERP/difference wave is not equal to the mean of the null 0039 % hypothesis). "-1" means lower-tailed (i.e., alternative 0040 % hypothesis is that the ERP/difference wave is less than 0041 % the mean of the null hypothesis). {default: 0} 0042 % alpha - A number between 0 and 1 specifying the family-wise 0043 % error rate (FWER) of the test. Note, cluster-based tests 0044 % only weak control of FWER. {default: 0.05} 0045 % thresh_p - The test-wise p-value threshold for cluster inclusion. If 0046 % a channel/time-point has a t-score that corresponds to an 0047 % uncorrected p-value greater than thresh_p, it is assigned 0048 % a p-value of 1 and not considered for clustering. Note 0049 % that thresh_p automatically takes into account the tail of 0050 % the test (e.g., you will get positive and negative t-score 0051 % thresholds for a two-tailed test). 0052 % chan_hood - A scalar or a 2D symmetric binary matrix that indicates 0053 % which channels are considered neighbors of other 0054 % channels. E.g., if chan_hood(2,10)=1, then Channel 2 0055 % and Channel 10 are nieghbors. You can produce a 0056 % chan_hood matrix using the function spatial_neighbors.m. 0057 % If a scalar is provided, then all electrodes within that 0058 % distance of a particular electrode are considered 0059 % neighbors. Note, EEGLAB's electrode coordinates assume 0060 % the head has a radius of 1. See the help documentation 0061 % of the function spatial_neighbors to see how you could 0062 % convert this distance threshold to centimeters. 0063 % {default: 0.61} 0064 % head_radius - The radius of the head in whatever units the Cartesian 0065 % coordinates in GRP.chanlocs are in. This is used to 0066 % convert scalar values of chan_hood into centimeters. 0067 % {default: []} 0068 % n_perm - The number of permutations to use in the test. As this 0069 % value increases, the test takes longer to compute and 0070 % the results become more reliable. Manly (1997) suggests 0071 % using at least 1000 permutations for an alpha level of 0072 % 0.05 and at least 5000 permutations for an alpha level 0073 % of 0.01. {default: 2500} 0074 % time_wind - Pair of time values specifying the beginning and 0075 % end of a time window in ms (e.g., [160 180]). Every 0076 % single time point in the time window will be individually 0077 % tested (i.e., maximal temporal resolution) if mean_wind 0078 % option is NOT used. Note, boundaries of time window 0079 % may not exactly correspond to desired time window 0080 % boundaries because of temporal digitization (i.e., you 0081 % only have samples every so many ms). {default: 0 ms to 0082 % the end of the epoch} 0083 % mean_wind - ['yes' or 'no'] If 'yes', the permutation test will be 0084 % performed on the mean amplitude within the time window 0085 % specified by time_wind. This sacrifices temporal 0086 % resolution to increase test power by reducing the number 0087 % of comparisons. If 'no', every single time point within 0088 % time_wind's time windows will be tested individually. 0089 % {default: 'no'} 0090 % null_mean - [number] The mean of the null hypothesis (i.e., mu) in 0091 % units of microvolts. {default: 0} 0092 % exclude_chans - A cell array of channel labels to exclude from the 0093 % permutation test (e.g., {'A2','lle','rhe'}). This option 0094 % sacrifices spatial resolution to increase test power by 0095 % reducing the number of comparisons. Use headinfo.m to see 0096 % the channel labels stored in the GND variable. You cannot 0097 % use both this option and 'include_chans' (below).{default: 0098 % not used, all channels included in test} 0099 % include_chans - A cell array of channel labels to use in the permutation 0100 % test (e.g., {'A2','lle','rhe'}). All other channels will 0101 % be ignored. This option sacrifices spatial resolution to 0102 % increase test power by reducing the number of comparisons. 0103 % Use headinfo.m to see the channel labels stored in the GND 0104 % variable. You cannot use both this option and 0105 % 'exclude_chans' (above). {default: not used, all channels 0106 % included in test} 0107 % verblevel - An integer specifiying the amount of information you want 0108 % this function to provide about what it is doing during runtime. 0109 % Options are: 0110 % 0 - quiet, only show errors, warnings, and EEGLAB reports 0111 % 1 - stuff anyone should probably know 0112 % 2 - stuff you should know the first time you start working 0113 % with a data set {default value} 0114 % 3 - stuff that might help you debug (show all 0115 % reports) 0116 % plot_raster - ['yes' or 'no'] If 'yes', a two-dimensional (time x channel) 0117 % binary "raster" diagram is created to illustrate the 0118 % results of the permutation test. Significant negative and 0119 % positive deviations from the null hypothesis are shown 0120 % as black and white rectangles respectively. Non- 0121 % significant comparisons are shown as gray rectangles. 0122 % Clicking on the rectangles will show you the 0123 % corresponding time and channel label for that 0124 % rectangle. This figure can be reproduced with the 0125 % function sig_raster.m. {default: 'yes'} 0126 % plot_mn_topo - ['yes' or 'no'] If 'yes', the topography of the mean 0127 % voltages/effects in the time window is produced. More 0128 % specifically, two figures are produced: one showing the 0129 % topographies in uV the other in t-scores. Significant/ 0130 % nonsignificant comparisons are shown as white/black 0131 % electrodes. Clicking on electrodes will show the 0132 % electrode's name. This figure can be reproduced with 0133 % the function sig_topo.m. This option has NO effect if 0134 % mean_wind option is set to 'no'. {default: 'yes'} 0135 % output_file - A string indicating the name of a space delimited text 0136 % file to produce containing the p-values of all comparisons 0137 % and the details of the test (e.g., number of permutations, 0138 % family-wise alpha level, etc...). If mean_wind option is 0139 % set to 'yes,' t-scores of each comparison are also 0140 % included since you cannot derive them from the t-scores 0141 % at each time point/electrode in a simple way. When 0142 % importing this file into a spreadsheet be sure NOT to count 0143 % consecutive spaces as multiple delimiters. {default: 0144 % none} 0145 % save_GND - ['yes' or 'no'] If 'yes', the GND variable will be 0146 % saved to disk after the permutation test is completed 0147 % and added to it. User will first be prompted to verify 0148 % file name and path. {default: 'yes'} 0149 % reproduce_test- [integer] The number of the permutation test stored in 0150 % the GND variable to reproduce. For example, if 0151 % 'reproduce_test' equals 2, the second t-test 0152 % stored in the GND variable (i.e., GND.t_tests(2)) will 0153 % be reproduced. Reproduction is accomplished by setting 0154 % the random number generator used in the permutation test 0155 % to the same initial state it was in when the permutation 0156 % test was first applied. 0157 % 0158 % Outputs: 0159 % GND - GND structure variable. This is the same as 0160 % the input GND variable with one addition: the 0161 % field GND.t_tests will contain the results of the 0162 % permutation test and the test parameters. 0163 % prm_pval - A two-dimensional matrix (channel x time) of the 0164 % p-values of each comparison. These p-values are 0165 % corrected for multiple comparisons by the permutation 0166 % test. 0167 % data_t - A two-dimensional matrix (channel x time) of the 0168 % t-scores of each comparison. 0169 % 0170 % Note also that a great deal of information about the test is displayed 0171 % in the MATLAB command window. You can easiy record of all this 0172 % information to a text file using the MATLAB command "diary." 0173 % 0174 % Global Variables: 0175 % VERBLEVEL - Mass Univariate ERP Toolbox level of verbosity (i.e., tells 0176 % functions how much to report about what they're doing during 0177 % runtime) set by the optional function argument 'verblevel' 0178 % 0179 % Notes: 0180 % -To add a difference wave to a GND variable, use the function "bin_dif". 0181 % 0182 % -Unlike a parametric test (e.g., an ANOVA), a discrete set of p-values 0183 % are possible (at most the number of possible permutations). Since the 0184 % number of possible permutations grows rapdily with the number of 0185 % participants, this is only issue for small sample sizes (e.g., 6 0186 % participants). When you have such a small sample size, the 0187 % limited number of p-values may make the test less conservative (e.g., 0188 % you might be forced to use an alpha level of .0286 since it is the biggest 0189 % possible alpha level less than .05). 0190 % 0191 % -For two-tailed tests clustGND.m effectively performs an upper-tailed 0192 % test and a lower-tailed test and multiplies the resulting p-values by 2. 0193 % In other words, it does two tests and then applies Bonferroni correction. 0194 % This is faster than performing a "proper" two-tailed permutation test 0195 % since you only have to form clusters for one polarity (e.g., if you form 0196 % negative clusters for each permutation, the negative of that distribution 0197 % will be valid for computing p-values for postive clusters) but can 0198 % result in p-values greater than 1. 0199 % 0200 % Author: 0201 % David Groppe 0202 % May, 2011 0203 % Kutaslab, San Diego 0204 % 0205 % References: 0206 % Bullmore, E. T., Suckling, J., Overmeyer, S., Rabe-Hesketh, S., Taylor, 0207 % E., & Brammer, M. J. (1999). Global, voxel, and cluster tests, by theory 0208 % and permutation, for a difference between two groups of structural MR 0209 % images of the brain. IEEE Transactions on Medical Imaging, 18(1), 32-42. 0210 % doi:10.1109/42.750253 0211 % 0212 % Maris, E., & Oostenveld, R. (2007). Nonparametric statistical testing of 0213 % EEG- and MEG-data. Journal of Neuroscience Methods, 164(1), 177-190. 0214 % doi:10.1016/j.jneumeth.2007.03.024 0215 % 0216 % Manly, B.F.J. (1997) Randomization, bootstrap, and Monte Carlo methods in 0217 % biology. 2nd ed. Chapmn and Hall, London. 0218 0219 %%%%%%%%%%%%%%%% REVISION LOG %%%%%%%%%%%%%%%%% 0220 % 12/11/2011-Now uses Amy Guthormsen's recursiveless find_clusters.m. 0221 % 0222 0223 function [GND, prm_pval, data_t]=clustGND(GND_or_fname,bin,varargin) 0224 0225 global VERBLEVEL; 0226 0227 p=inputParser; 0228 p.addRequired('GND_or_fname',@(x) ischar(x) || isstruct(x)); 0229 p.addRequired('bin',@(x) isnumeric(x) && (length(x)==1) && (x>0)); 0230 p.addParamValue('tail',0,@(x) isnumeric(x) && (length(x)==1)); 0231 p.addParamValue('alpha',0.05,@(x) (length(x)==1) && isnumeric(x) && (x>0) && (x<1)); 0232 p.addParamValue('thresh_p',0.05,@(x) (length(x)==1) && isnumeric(x) && (x>0) && (x<1)); 0233 p.addParamValue('chan_hood',0.61,@isnumeric); 0234 p.addParamValue('head_radius',[],@isnumeric); 0235 p.addParamValue('time_wind',[],@(x) isnumeric(x) && (size(x,2)==2)); 0236 p.addParamValue('mean_wind','no',@(x) ischar(x) && (strcmpi(x,'yes') || strcmpi(x,'no'))); 0237 p.addParamValue('n_perm',2500,@(x) isnumeric(x) && (length(x)==1)); 0238 p.addParamValue('null_mean',0,@(x) isnumeric(x) && (length(x)==1)); 0239 p.addParamValue('verblevel',[],@(x) isnumeric(x) && (length(x)==1)); 0240 p.addParamValue('exclude_chans',[],@(x) ischar(x) || iscell(x)); 0241 p.addParamValue('include_chans',[],@(x) ischar(x) || iscell(x)); 0242 p.addParamValue('plot_raster','yes',@(x) ischar(x) && (strcmpi(x,'yes') || strcmpi(x,'no'))); 0243 p.addParamValue('plot_gui','yes',@(x) ischar(x) && (strcmpi(x,'yes') || strcmpi(x,'no'))); 0244 p.addParamValue('plot_mn_topo',[],@(x) ischar(x) && (strcmpi(x,'yes') || strcmpi(x,'no'))); 0245 p.addParamValue('output_file',[],@ischar); 0246 p.addParamValue('reproduce_test',[],@(x) isnumeric(x) && (length(x)==1)); 0247 p.addParamValue('save_GND','yes',@(x) ischar(x) && (strcmpi(x,'yes') || strcmpi(x,'no'))); 0248 0249 p.parse(GND_or_fname,bin,varargin{:}); 0250 0251 if isempty(p.Results.verblevel), 0252 if isempty(VERBLEVEL), 0253 VERBLEVEL=2; 0254 end 0255 else 0256 VERBLEVEL=p.Results.verblevel; 0257 end 0258 0259 mean_wind=str2bool(p.Results.mean_wind); 0260 0261 %Load GND struct 0262 if ischar(GND_or_fname), 0263 fprintf('Loading GND struct from file %s.\n',GND_or_fname); 0264 load(GND_or_fname,'-MAT'); 0265 if ~exist('GND') 0266 error('File %s does not contain a GND variable.',GND_or_fname); 0267 end 0268 else 0269 GND=GND_or_fname; 0270 clear GND_or_fname; 0271 fldnames=fieldnames(GND); 0272 if ismember('group_desc',fldnames), 0273 error('You passed a GRP variable to this function instead of a GND variable.'); 0274 end 0275 end 0276 [n_chan, n_pt, n_bin, total_subs]=size(GND.indiv_erps); 0277 VerbReport(sprintf('Experiment: %s',GND.exp_desc),2,VERBLEVEL); 0278 0279 if (bin>n_bin), 0280 error('There is no Bin %d in this GND variable.',bin); 0281 end 0282 0283 %Use only subs with data in relevant bin(s) 0284 use_subs=find(GND.indiv_bin_ct(:,p.Results.bin)); 0285 n_sub=length(use_subs); 0286 VerbReport(sprintf('%d out of %d participants have data in relevant bin.',n_sub,total_subs), ... 0287 1,VERBLEVEL); 0288 0289 0290 %% Figure out which channels to ignore if any 0291 %But first make sure exclude & include options were not both used. 0292 if ~isempty(p.Results.include_chans) && ~isempty(p.Results.exclude_chans) 0293 error('You cannot use BOTH ''include_chans'' and ''exclude_chans'' options.'); 0294 end 0295 if ischar(p.Results.exclude_chans), 0296 exclude_chans{1}=p.Results.exclude_chans; 0297 elseif isempty(p.Results.exclude_chans) 0298 exclude_chans=[]; 0299 else 0300 exclude_chans=p.Results.exclude_chans; 0301 end 0302 if ischar(p.Results.include_chans), 0303 include_chans{1}=p.Results.include_chans; 0304 elseif isempty(p.Results.include_chans) 0305 include_chans=[]; 0306 else 0307 include_chans=p.Results.include_chans; 0308 end 0309 0310 % exclude and include chan options 0311 if ~isempty(exclude_chans), 0312 ignore_chans=zeros(1,length(exclude_chans)); %preallocate mem 0313 ct=0; 0314 for x=1:length(exclude_chans), 0315 found=0; 0316 for c=1:n_chan, 0317 if strcmpi(exclude_chans{x},GND.chanlocs(c).labels), 0318 found=1; 0319 ct=ct+1; 0320 ignore_chans(ct)=c; 0321 end 0322 end 0323 if ~found, 0324 watchit(sprintf('I attempted to exclude %s. However no such electrode was found in GND variable.', ... 0325 exclude_chans{x})); 0326 end 0327 end 0328 ignore_chans=ignore_chans(1:ct); 0329 use_chans=setdiff(1:n_chan,ignore_chans); 0330 elseif ~isempty(include_chans), 0331 use_chans=zeros(1,length(include_chans)); %preallocate mem 0332 ct=0; 0333 for x=1:length(include_chans), 0334 found=0; 0335 for c=1:n_chan, 0336 if strcmpi(include_chans{x},GND.chanlocs(c).labels), 0337 found=1; 0338 ct=ct+1; 0339 use_chans(ct)=c; 0340 end 0341 end 0342 if ~found, 0343 watchit(sprintf('I attempted to include %s. However no such electrode was found in GND variable.', ... 0344 include_chans{x})); 0345 end 0346 end 0347 use_chans=use_chans(1:ct); 0348 else 0349 use_chans=1:n_chan; 0350 end 0351 0352 % Establish spatial neighborhood matrix for making clusters 0353 if isscalar(p.Results.chan_hood), 0354 chan_hood=spatial_neighbors(GND.chanlocs(use_chans),p.Results.chan_hood,p.Results.head_radius); 0355 else 0356 chan_hood=p.Results.chan_hood; 0357 end 0358 0359 0360 %% Find time points 0361 if isempty(p.Results.time_wind), 0362 time_wind=[0 GND.time_pts(end)]; %default time window 0363 else 0364 time_wind=p.Results.time_wind; 0365 end 0366 time_wind=sort(time_wind,2); %first make sure earlier of each pair of time points is first 0367 time_wind=sort(time_wind,1); %next sort time windows from earliest to latest onset 0368 n_wind=size(time_wind,1); 0369 if n_wind>1, 0370 error('clustGND.m can only handle a single mean time window at the moment. Try tmaxGND.m or tfdrGND.m if you want to simultaneously test hypotheses in multiple mean time windows.'); 0371 end 0372 0373 if mean_wind, 0374 use_tpts=cell(1,n_wind); 0375 else 0376 use_tpts=[]; 0377 end 0378 for a=1:n_wind, 0379 VerbReport(sprintf('Time Window #%d:',a),1,VERBLEVEL); 0380 VerbReport(sprintf('Attempting to use time boundaries of %d to %d ms for hypothesis test.',time_wind(a,1),time_wind(a,2)), ... 0381 1,VERBLEVEL); 0382 start_tpt=find_tpt(time_wind(a,1),GND.time_pts); 0383 end_tpt=find_tpt(time_wind(a,2),GND.time_pts); 0384 if mean_wind, 0385 use_tpts{a}=[start_tpt:end_tpt]; 0386 else 0387 use_tpts=[use_tpts [start_tpt:end_tpt]]; 0388 end 0389 %replace desired time points with closest matches 0390 time_wind(a,1)=GND.time_pts(start_tpt); 0391 time_wind(a,2)=GND.time_pts(end_tpt); 0392 VerbReport(sprintf('Exact window boundaries are %d to %d ms (that''s from time point %d to %d).', ... 0393 time_wind(a,1),time_wind(a,2),start_tpt,end_tpt),1,VERBLEVEL); 0394 end 0395 if ~mean_wind, 0396 use_tpts=unique(use_tpts); %sorts time points and gets rid of any redundant time points 0397 end 0398 0399 0400 %% Compile data 0401 if mean_wind, 0402 %Take mean amplitude in time blocks and then test 0403 erps=zeros(length(use_chans),n_wind,n_sub); 0404 for a=1:n_wind, 0405 for sub=1:n_sub, 0406 erps(:,a,sub)=mean(GND.indiv_erps(use_chans,use_tpts{a},bin,use_subs(sub)),2); 0407 end 0408 end 0409 else 0410 %Use every single time point in time window(s) 0411 n_use_tpts=length(use_tpts); 0412 erps=zeros(length(use_chans),n_use_tpts,n_sub); 0413 for sub=1:n_sub, 0414 erps(:,:,sub)=GND.indiv_erps(use_chans,use_tpts,bin,use_subs(sub)); 0415 end 0416 end 0417 0418 0419 %% Report tail of test & alpha levels 0420 VerbReport(sprintf('Testing null hypothesis that the grand average ERPs in Bin %d (%s) have a mean of %f microvolts.',bin, ... 0421 GND.bin_info(bin).bindesc,p.Results.null_mean),1,VERBLEVEL); 0422 if p.Results.tail==0 0423 VerbReport(sprintf('Alternative hypothesis is that the ERPs differ from %f (i.e., two-tailed test).',p.Results.null_mean), ... 0424 1,VERBLEVEL); 0425 elseif p.Results.tail<0, 0426 VerbReport(sprintf('Alternative hypothesis is that the ERPs are less than %f (i.e., lower-tailed test).',p.Results.null_mean), ... 0427 1,VERBLEVEL); 0428 else 0429 VerbReport(sprintf('Alternative hypothesis is that the ERPs are greater than %f (i.e., upper-tailed test).',p.Results.null_mean), ... 0430 1,VERBLEVEL); 0431 end 0432 0433 0434 %% Optionally reset random number stream to reproduce a previous test 0435 if isempty(p.Results.reproduce_test), 0436 seed_state=[]; 0437 else 0438 if p.Results.reproduce_test>length(GND.t_tests), 0439 error('Value of argument ''reproduce_test'' is too high. You only have %d t-tests stored with this GND variable.',length(GND.t_tests)); 0440 else 0441 if isnan(GND.t_tests(p.Results.reproduce_test).n_perm) 0442 error('t-test set %d is NOT a permutation test. You don''t need to seed the random number generator to reproduce it.', ... 0443 p.Results.reproduce_test); 0444 else 0445 seed_state=GND.t_tests(p.Results.reproduce_test).seed_state; 0446 end 0447 end 0448 end 0449 0450 %% Compute the permutation test 0451 [prm_pval, data_t, clust_info, seed_state, est_alpha]=clust_perm1(erps-p.Results.null_mean, ... 0452 chan_hood,p.Results.n_perm,p.Results.alpha,p.Results.tail,p.Results.thresh_p, ... 0453 VERBLEVEL,seed_state,0); 0454 0455 %% Command line summary of results 0456 if p.Results.tail>=0, 0457 %upper or two-tailed test 0458 n_pos=length(clust_info.pos_clust_pval); 0459 fprintf('# of positive clusters found: %d\n',n_pos); 0460 fprintf('# of significant positive clusters found: %d\n',sum(clust_info.pos_clust_pval<est_alpha)); 0461 fprintf('Positive cluster p-values range from %g to %g.\n',min(clust_info.pos_clust_pval), ... 0462 max(clust_info.pos_clust_pval)); 0463 end 0464 if p.Results.tail<=0, 0465 %lower or two-tailed test 0466 n_neg=length(clust_info.neg_clust_pval); 0467 fprintf('# of negative clusters found: %d\n',n_neg); 0468 fprintf('# of significant negative clusters found: %d\n',sum(clust_info.neg_clust_pval<est_alpha)); 0469 fprintf('Negative cluster p-values range from %g to %g.\n',min(clust_info.neg_clust_pval), ... 0470 max(clust_info.neg_clust_pval)); 0471 end 0472 0473 sig_ids=find(prm_pval<p.Results.alpha); 0474 if isempty(sig_ids), 0475 fprintf('ERPs are NOT significantly different from zero (alpha=%f) at any time point/window analyzed.\n', ... 0476 p.Results.alpha); 0477 fprintf('All p-values>=%f\n',min(min(prm_pval))); 0478 else 0479 [dummy min_t_id]=min(abs(data_t(sig_ids))); 0480 min_t=data_t(sig_ids(min_t_id)); 0481 VerbReport(['Smallest significant t-score(s):' num2str(min_t)],1,VERBLEVEL); 0482 if p.Results.tail 0483 %one-tailed test 0484 tw_alpha=1-cdf('t',max(abs(min_t)),n_sub-1); 0485 else 0486 %two-tailed test 0487 tw_alpha=(1-cdf('t',max(abs(min_t)),n_sub-1))*2; 0488 end 0489 0490 VerbReport(sprintf('That corresponds to a test-wise alpha level of %f.',tw_alpha),1,VERBLEVEL); 0491 VerbReport(sprintf('Bonferroni test-wise alpha would be %f.',p.Results.alpha/(size(prm_pval,1)* ... 0492 size(prm_pval,2))),1,VERBLEVEL); 0493 sig_tpts=find(sum(prm_pval<p.Results.alpha)); 0494 0495 fprintf('Significant differences from zero (in order of earliest to latest):\n'); 0496 max_sig_p=0; 0497 min_sig_p=2; 0498 for t=sig_tpts, 0499 if mean_wind 0500 %time windows instead of time points 0501 fprintf('%d to %d ms, electrode(s): ',GND.time_pts(use_tpts{t}(1)), ... 0502 GND.time_pts(use_tpts{t}(end))); 0503 else 0504 fprintf('%d ms, electrode(s): ',GND.time_pts(use_tpts(t))); 0505 end 0506 sig_elec=find(prm_pval(:,t)<p.Results.alpha); 0507 ct=0; 0508 for c=sig_elec', 0509 ct=ct+1; 0510 if prm_pval(c,t)>max_sig_p, 0511 max_sig_p=prm_pval(c,t); 0512 end 0513 if prm_pval(c,t)<min_sig_p, 0514 min_sig_p=prm_pval(c,t); 0515 end 0516 0517 if ct==length(sig_elec), 0518 fprintf('%s.\n',GND.chanlocs(use_chans(c)).labels); 0519 else 0520 fprintf('%s, ',GND.chanlocs(use_chans(c)).labels); 0521 end 0522 end 0523 end 0524 fprintf('All significant corrected p-values are between %f and %f\n',max_sig_p,min_sig_p); 0525 end 0526 0527 0528 %% Add permutation results to GND struct 0529 n_t_tests=length(GND.t_tests); 0530 neo_test=n_t_tests+1; 0531 GND.t_tests(neo_test).bin=bin; 0532 GND.t_tests(neo_test).time_wind=time_wind; 0533 GND.t_tests(neo_test).used_tpt_ids=use_tpts; 0534 n_use_chans=length(use_chans); 0535 include_chans=cell(1,n_use_chans); 0536 for a=1:n_use_chans, 0537 include_chans{a}=GND.chanlocs(use_chans(a)).labels; 0538 end 0539 GND.t_tests(neo_test).include_chans=include_chans; 0540 GND.t_tests(neo_test).used_chan_ids=use_chans; 0541 GND.t_tests(neo_test).mult_comp_method='cluster mass perm test'; 0542 GND.t_tests(neo_test).n_perm=p.Results.n_perm; 0543 GND.t_tests(neo_test).desired_alphaORq=p.Results.alpha; 0544 GND.t_tests(neo_test).estimated_alpha=est_alpha; 0545 GND.t_tests(neo_test).null_mean=p.Results.null_mean; 0546 if mean_wind, 0547 GND.t_tests(neo_test).data_t=data_t; 0548 GND.t_tests(neo_test).mean_wind='yes'; 0549 else 0550 GND.t_tests(neo_test).data_t='See GND.grands_t'; 0551 GND.t_tests(neo_test).mean_wind='no'; 0552 end 0553 GND.t_tests(neo_test).crit_t=NaN; 0554 GND.t_tests(neo_test).df=length(use_subs)-1; 0555 GND.t_tests(neo_test).adj_pval=prm_pval; 0556 GND.t_tests(neo_test).fdr_rej=NaN; 0557 GND.t_tests(neo_test).seed_state=seed_state; 0558 GND.t_tests(neo_test).clust_info=clust_info; 0559 GND.t_tests(neo_test).chan_hood=chan_hood; 0560 0561 if strcmpi(p.Results.plot_raster,'yes'), 0562 sig_raster(GND,neo_test,'verblevel',0,'use_color','rgb'); 0563 end 0564 0565 if mean_wind, 0566 if strcmpi(p.Results.plot_mn_topo,'yes') || isempty(p.Results.plot_mn_topo), 0567 sig_topo(GND,neo_test,'units','t','verblevel',0); %t-score topographies 0568 sig_topo(GND,neo_test,'units','uV','verblevel',0); %microvolt topographies 0569 end 0570 else 0571 if strcmpi(p.Results.plot_gui,'yes'), 0572 gui_erp('initialize','GNDorGRP',GND,'t_test',neo_test,'stat','t', ... 0573 'verblevel',1); 0574 end 0575 end 0576 0577 %% Create text file output 0578 if ~isempty(p.Results.output_file) 0579 [fid msg]=fopen(p.Results.output_file,'w'); 0580 if fid==-1, 0581 error('Cannot create file %s for writing. According to fopen.m: %s.', ... 0582 p.Results.file,msg); 0583 else 0584 %Write header column of times 0585 % Leave first column blank for channel labels 0586 if mean_wind, 0587 for t=1:n_wind 0588 fprintf(fid,' %d-%d',GND.time_pts(use_tpts{t}(1)), ... 0589 GND.time_pts(use_tpts{t}(end))); 0590 end 0591 0592 %write a couple spaces and then write header for t-scores 0593 fprintf(fid,' '); 0594 for t=1:n_wind 0595 fprintf(fid,' %d-%d',GND.time_pts(use_tpts{t}(1)), ... 0596 GND.time_pts(use_tpts{t}(end))); 0597 end 0598 else 0599 for t=use_tpts, 0600 fprintf(fid,' %d',GND.time_pts(t)); 0601 end 0602 end 0603 fprintf(fid,' Milliseconds\n'); 0604 0605 % Write channel labels and p-values 0606 chan_ct=0; 0607 for c=use_chans, 0608 chan_ct=chan_ct+1; 0609 fprintf(fid,'%s',GND.chanlocs(c).labels); 0610 for t=1:length(use_tpts), 0611 fprintf(fid,' %f',prm_pval(chan_ct,t)); 0612 end 0613 fprintf(fid,' p-value'); 0614 0615 if mean_wind, 0616 %write a couple spaces and then write t-scores if mean amp 0617 %in time windows used 0618 fprintf(fid,' '); 0619 for t=1:n_wind 0620 fprintf(fid,' %f',data_t(chan_ct,t)); 0621 end 0622 fprintf(fid,' t-score \n'); 0623 else 0624 fprintf(fid,'\n'); 0625 end 0626 end 0627 0628 % Write permutation test details 0629 fprintf(fid,'Experiment: %s\n',GND.exp_desc); 0630 fprintf(fid,'Test_of_null_hypothesis_that_Bin_%d_equals: %f\n',bin,p.Results.null_mean); 0631 fprintf(fid,'#_of_time_windows: %d\n',n_wind); 0632 fprintf(fid,'#_of_permutations: %d\n',p.Results.n_perm); 0633 fprintf(fid,'Tail_of_test: '); 0634 if ~p.Results.tail, 0635 fprintf(fid,'Two_tailed\n'); 0636 elseif p.Results.tail>0 0637 fprintf(fid,'Upper_tailed\n'); 0638 else 0639 fprintf(fid,'Lower_tailed\n'); 0640 end 0641 fprintf(fid,'Degrees_of_freedom: %d\n',length(use_subs)-1); 0642 fprintf(fid,'Family_wise_error_rate: %f\n',p.Results.alpha); 0643 fprintf(fid,'Cluster_inclusion_p_value_threshold: %f\n',p.Results.thresh_p); 0644 if isscalar(p.Results.chan_hood) 0645 fprintf(fid,'Max_spatial_neighbor_distance: %f\n',p.Results.chan_hood); 0646 end 0647 0648 % # of participants and filenames 0649 fprintf(fid,'#_of_participants: %d\n',length(use_subs)); 0650 fprintf(fid,'Participant_names: \n'); 0651 for s=1:length(use_subs), 0652 fprintf(fid,'%s\n',GND.indiv_subnames{use_subs(s)}); 0653 end 0654 end 0655 fclose(fid); 0656 end 0657 0658 if ~strcmpi(p.Results.save_GND,'no'), 0659 GND=save_matmk(GND); 0660 end 0661 0662 % 0663 %% %%%%%%%%%%%%%%%%%%%%% function str2bool() %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0664 % 0665 function bool=str2bool(str) 0666 %function bool=str2bool(str) 0667 0668 if ischar(str), 0669 if strcmpi(str,'yes') || strcmpi(str,'y') 0670 bool=1; 0671 else 0672 bool=0; 0673 end 0674 else 0675 bool=str; 0676 end 0677