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