tmaxGRP() - 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 permutation test based on the
tmax statistic (e.g., Hemmelmann et al., 2004). 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 tmaxGND.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, crit_t]=tmaxGRP(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}
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 - 2D matrix of pairs of time values specifying the beginning
and end of one or more time windows in ms (e.g.,
[160 180; 350 550]). 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(s) 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}
time_block_dur- [integers] A number or numbers (in milliseconds)
specifying a duration of time blocks in which the time
windows specified by time_wind will be divided. For
example, if time_wind=[300 600] and time_block_dur=100,
the 300 to 600 ms time window will be sub-divided into
100 ms windows (i.e., time_wind will equal [300 396;
400 496; 500 596]). This is an easy way to break up
larger time windows of interest into smaller windows
for mean window analysis (if mean_wind option is set to
'yes'). If you specify multiple time windows with
time_wind, you can break them up using durations of
different lengths. For example, if time_wind=[150 250;
400 900] and time_block_dur=[25 100], the first time
window (150 to 250 ms) will be broken up into 25 ms
windows and the second window (400 to 900 ms) will be
broken up into 100 ms windows. {default: not used}
mean_wind - ['yes' or 'no'] If 'yes', the permutation test will be
performed on the mean amplitude within each 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 topographies of the mean
voltages/effects in each 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.
crit_t - The critical t-score(s) for the test. Any t-scores
that are more extreme than the critical t-score(s)
significantly deviate from 0.
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).
References:
Manly, B.F.J. (1997) Randomization, bootstrap, and Monte Carlo methods in
biology. 2nd ed. Chapmn and Hall, London.
Hemmelmann, et al. (2004) Multivariate tests for the evaluation of
high-dimensional EEG data. Journal of Neuroscience Methods.
Author:
David Groppe
Kutaslab, 3/20100001 % tmaxGRP() - 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 permutation test based on the 0004 % tmax statistic (e.g., Hemmelmann et al., 2004). This function 0005 % requires individual subject ERPs to be stored in a "GRP" 0006 % structure and outputs the test results in a number of graphical 0007 % and text formats. For analogous within-subject comparisons use 0008 % the function tmaxGND.m. 0009 % Note, when applied to a bin that is the mean ERP across two 0010 % groups (i.e., NOT a difference wave), a one-sample test is 0011 % executed. 0012 % 0013 % Usage: 0014 % >> [GRP, prm_pval, data_t, crit_t]=tmaxGRP(GRP_or_fname,bin,varargin) 0015 % 0016 % Required Inputs: 0017 % GRP_or_fname - A GRP structure variable or the filename of a GRP 0018 % structure that has been saved to disk. A GRP variable 0019 % is based on GND variables. To create a GRP variable from 0020 % GND variables use GNDs2GRP.m. See Mass Univariate ERP 0021 % Toolbox documentation for detailed information about the 0022 % format of a GRP variable. If you specifiy a filename be 0023 % sure to include the file's path, unless the file is in 0024 % the current working directory. 0025 % bin - [integer] The bin to contrast against the mean of the 0026 % null hypothesis. Use the function headinfo.m to see what 0027 % bins are stored in a GRP variable. Use the function 0028 % bin_opGRP.m to create a difference wave between two bins 0029 % whose significance you can test with this function. 0030 % 0031 % Optional Inputs: 0032 % tail - [1 | 0 | -1] An integer specifying the tail of the 0033 % hypothesis test. "1" means upper-tailed (i.e., alternative 0034 % hypothesis is that the ERP/difference wave is greater 0035 % than mu). "0" means two-tailed (i.e., alternative hypothesis 0036 % is that the ERP/difference wave is not equal to mu). 0037 % "-1" means lower-tailed (i.e., alternative hypothesis 0038 % is that the ERP/difference wave is less than mu). 0039 % {default: 0} 0040 % alpha - A number between 0 and 1 specifying the family-wise 0041 % alpha level of the test. {default: 0.05} 0042 % n_perm - The number of permutations to use in the test. As this 0043 % value increases, the test takes longer to compute and 0044 % the results become more reliable. Manly (1997) suggests 0045 % using at least 1000 permutations for an alpha level of 0046 % 0.05 and at least 5000 permutations for an alpha level 0047 % of 0.01. {default: 2500} 0048 % time_wind - 2D matrix of pairs of time values specifying the beginning 0049 % and end of one or more time windows in ms (e.g., 0050 % [160 180; 350 550]). Every single time point in 0051 % the time window will be individually tested (i.e., 0052 % maximal temporal resolution) if mean_wind option is 0053 % NOT used. Note, boundaries of time window(s) may not 0054 % exactly correspond to desired time window boundaries 0055 % because of temporal digitization (i.e., you only have 0056 % samples every so many ms). {default: 0 ms to the end of 0057 % the epoch} 0058 % time_block_dur- [integers] A number or numbers (in milliseconds) 0059 % specifying a duration of time blocks in which the time 0060 % windows specified by time_wind will be divided. For 0061 % example, if time_wind=[300 600] and time_block_dur=100, 0062 % the 300 to 600 ms time window will be sub-divided into 0063 % 100 ms windows (i.e., time_wind will equal [300 396; 0064 % 400 496; 500 596]). This is an easy way to break up 0065 % larger time windows of interest into smaller windows 0066 % for mean window analysis (if mean_wind option is set to 0067 % 'yes'). If you specify multiple time windows with 0068 % time_wind, you can break them up using durations of 0069 % different lengths. For example, if time_wind=[150 250; 0070 % 400 900] and time_block_dur=[25 100], the first time 0071 % window (150 to 250 ms) will be broken up into 25 ms 0072 % windows and the second window (400 to 900 ms) will be 0073 % broken up into 100 ms windows. {default: not used} 0074 % mean_wind - ['yes' or 'no'] If 'yes', the permutation test will be 0075 % performed on the mean amplitude within each time window 0076 % specified by time_wind. This sacrifices temporal 0077 % resolution to increase test power by reducing the number 0078 % of comparisons. If 'no', every single time point within 0079 % time_wind's time windows will be tested individually. 0080 % {default: 'no'} 0081 % null_mean - [number] The mean of the null hypothesis (in units of 0082 % microvolts). {default: 0} 0083 % exclude_chans - A cell array of channel labels to exclude from the 0084 % permutation test (e.g., {'A2','lle','rhe'}). This option 0085 % sacrifices spatial resolution to increase test power by 0086 % reducing the number of comparisons. Use headinfo.m to see 0087 % the channel labels stored in the GRP variable. You cannot 0088 % use both this option and 'include_chans' (below).{default: 0089 % not used, all channels included in test} 0090 % include_chans - A cell array of channel labels to use in the permutation 0091 % test (e.g., {'A2','lle','rhe'}). All other channels will 0092 % be ignored. This option sacrifices spatial resolution to 0093 % increase test power by reducing the number of comparisons. 0094 % Use headinfo.m to see the channel labels stored in the GRP 0095 % variable. You cannot use both this option and 0096 % 'exclude_chans' (above). {default: not used, all channels 0097 % included in test} 0098 % verblevel - An integer specifiying the amount of information you want 0099 % this function to provide about what it is doing during runtime. 0100 % Options are: 0101 % 0 - quiet, only show errors, warnings, and EEGLAB reports 0102 % 1 - stuff anyone should probably know 0103 % 2 - stuff you should know the first time you start working 0104 % with a data set {default value} 0105 % 3 - stuff that might help you debug (show all 0106 % reports) 0107 % plot_gui - ['yes' or 'no'] If 'yes', a GUI is created for 0108 % visualizing the results of the permutation test using the 0109 % function gui_erp.m. The GUI vizualizes the grand average 0110 % ERPs in each bin via various stats (uV, t-scores), shows 0111 % topographies at individual time points, and illustrates 0112 % which electrodes significantly differ from the null 0113 % hypothesis. This option does not work if mean_wind 0114 % option is set to 'yes.' This GUI can be reproduced using 0115 % the function gui_erp.m. {default: 'yes'} 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 topographies of the mean 0127 % voltages/effects in each 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: none} 0144 % save_GRP - ['yes' or 'no'] If 'yes', the GRP variable will be 0145 % saved to disk after the permutation test is completed 0146 % and added to it. User will first be prompted to verify 0147 % file name and path. {default: 'yes'} 0148 % reproduce_test- [integer] The number of the permutation test stored in 0149 % the GRP variable to reproduce. For example, if 0150 % 'reproduce_test' equals 2, the second t-test 0151 % stored in the GRP variable (i.e., GRP.t_tests(2)) will 0152 % be reproduced. Reproduction is accomplished by setting 0153 % the random number generator used in the permutation test 0154 % to the same initial state it was in when the permutation 0155 % test was first applied. 0156 % 0157 % Outputs: 0158 % GRP - GRP structure variable. This is the same as 0159 % the input GRP variable with one addition: the 0160 % field GRP.t_tests will contain the results of the 0161 % permutation test and the test parameters. 0162 % prm_pval - A two-dimensional matrix (channel x time) of the 0163 % p-values of each comparison. 0164 % data_t - A two-dimensional matrix (channel x time) of the 0165 % t-scores of each comparison. 0166 % crit_t - The critical t-score(s) for the test. Any t-scores 0167 % that are more extreme than the critical t-score(s) 0168 % significantly deviate from 0. 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 GRP variable, use the function "bin_opGRP". 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., 3 0186 % participants in each group). When you have such a small sample size, the 0187 % limited number of p-values may make the test overly 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 % References: 0192 % Manly, B.F.J. (1997) Randomization, bootstrap, and Monte Carlo methods in 0193 % biology. 2nd ed. Chapmn and Hall, London. 0194 % 0195 % Hemmelmann, et al. (2004) Multivariate tests for the evaluation of 0196 % high-dimensional EEG data. Journal of Neuroscience Methods. 0197 % 0198 % 0199 % Author: 0200 % David Groppe 0201 % Kutaslab, 3/2010 0202 0203 %%%%%%%%%%%%%%%% REVISION LOG %%%%%%%%%%%%%%%%% 0204 % 3/27/2010-'exclude_chans' field of GND.perm_tests removed (now there's 0205 % just a field for "include chans" 0206 % 0207 % 4/1/2010-Checks for compatibility of time points and channels between GND 0208 % and GRP variables. One-sample t-test performed on bins that are the average 0209 % of two groups. Added 'null_mean' option. 0210 % 0211 % 5/6/2010-Made compatible with FDR controls (e.g., GND.perm_tests turned 0212 % into GND.t_tests) 0213 % 0214 % 4/14/2011-Now reports min AND max significant p-values to command line 0215 % 0216 0217 function [GRP, prm_pval, data_t, crit_t]=tmaxGRP(GRP_or_fname,bin,varargin) 0218 0219 global VERBLEVEL; 0220 0221 p=inputParser; 0222 p.addRequired('GRP_or_fname',@(x) ischar(x) || isstruct(x)); 0223 p.addRequired('bin',@(x) isnumeric(x) && (length(x)==1) && (x>0)); 0224 p.addParamValue('tail',0,@(x) isnumeric(x) && (length(x)==1)); 0225 p.addParamValue('alpha',0.05,@(x) isnumeric(x) && (x>0) && (x<1)); 0226 p.addParamValue('time_wind',[],@(x) isnumeric(x) && (size(x,2)==2)); 0227 p.addParamValue('time_block_dur',[],@isnumeric); 0228 p.addParamValue('mean_wind','no',@(x) ischar(x) && (strcmpi(x,'yes') || strcmpi(x,'no'))); 0229 p.addParamValue('n_perm',2500,@(x) isnumeric(x) && (length(x)==1)); 0230 p.addParamValue('null_mean',0,@(x) isnumeric(x) && (length(x)==1)); 0231 p.addParamValue('verblevel',[],@(x) isnumeric(x) && (length(x)==1)); 0232 p.addParamValue('exclude_chans',[],@(x) ischar(x) || iscell(x)); 0233 p.addParamValue('include_chans',[],@(x) ischar(x) || iscell(x)); 0234 p.addParamValue('plot_gui','yes',@(x) ischar(x) && (strcmpi(x,'yes') || strcmpi(x,'no'))); 0235 p.addParamValue('plot_raster','yes',@(x) ischar(x) && (strcmpi(x,'yes') || strcmpi(x,'no'))); 0236 p.addParamValue('plot_mn_topo',[],@(x) ischar(x) && (strcmpi(x,'yes') || strcmpi(x,'no'))); 0237 p.addParamValue('output_file',[],@ischar); 0238 p.addParamValue('reproduce_test',[],@(x) isnumeric(x) && (length(x)==1)); 0239 p.addParamValue('save_GRP','yes',@(x) ischar(x) && (strcmpi(x,'yes') || strcmpi(x,'no'))); 0240 0241 p.parse(GRP_or_fname,bin,varargin{:}); 0242 0243 0244 if isempty(p.Results.verblevel), 0245 if isempty(VERBLEVEL), 0246 VERBLEVEL=2; 0247 end 0248 else 0249 VERBLEVEL=p.Results.verblevel; 0250 end 0251 0252 mean_wind=str2bool(p.Results.mean_wind); 0253 0254 %Load GRP struct 0255 if ischar(GRP_or_fname), 0256 fprintf('Loading GRP struct from file %s.\n',GRP_or_fname); 0257 load(GRP_or_fname,'-MAT'); 0258 else 0259 GRP=GRP_or_fname; 0260 clear GRP_or_fname; 0261 end 0262 [n_chan, n_pt, n_bin]=size(GRP.grands); 0263 n_group=length(GRP.GND_fnames); 0264 VerbReport(sprintf('Experiment: %s',GRP.exp_desc),2,VERBLEVEL); 0265 0266 if (bin>n_bin), 0267 error('There is no Bin %d in this GRP variable.',bin); 0268 end 0269 grpA=GRP.bin_info(bin).groupA; 0270 grpB=GRP.bin_info(bin).groupB; 0271 0272 0273 %% Figure out which channels to ignore if any 0274 %But first make sure exclude & include options were not both used. 0275 if ~isempty(p.Results.include_chans) && ~isempty(p.Results.exclude_chans) 0276 error('You cannot use BOTH ''include_chans'' and ''exclude_chans'' options.'); 0277 end 0278 if ischar(p.Results.exclude_chans), 0279 exclude_chans{1}=p.Results.exclude_chans; 0280 elseif isempty(p.Results.exclude_chans) 0281 exclude_chans=[]; 0282 else 0283 exclude_chans=p.Results.exclude_chans; 0284 end 0285 if ischar(p.Results.include_chans), 0286 include_chans{1}=p.Results.include_chans; 0287 elseif isempty(p.Results.include_chans) 0288 include_chans=[]; 0289 else 0290 include_chans=p.Results.include_chans; 0291 end 0292 0293 % exclude and include chan options 0294 if ~isempty(exclude_chans), 0295 ignore_chans=zeros(1,length(exclude_chans)); %preallocate mem 0296 ct=0; 0297 for x=1:length(exclude_chans), 0298 found=0; 0299 for c=1:n_chan, 0300 if strcmpi(exclude_chans{x},GRP.chanlocs(c).labels), 0301 found=1; 0302 ct=ct+1; 0303 ignore_chans(ct)=c; 0304 end 0305 end 0306 if ~found, 0307 watchit(sprintf('I attempted to exclude %s. However no such electrode was found in GRP variable.', ... 0308 exclude_chans{x})); 0309 end 0310 end 0311 ignore_chans=ignore_chans(1:ct); 0312 use_chans=setdiff(1:n_chan,ignore_chans); 0313 elseif ~isempty(include_chans), 0314 use_chans=zeros(1,length(include_chans)); %preallocate mem 0315 ct=0; 0316 for x=1:length(include_chans), 0317 found=0; 0318 for c=1:n_chan, 0319 if strcmpi(include_chans{x},GRP.chanlocs(c).labels), 0320 found=1; 0321 ct=ct+1; 0322 use_chans(ct)=c; 0323 end 0324 end 0325 if ~found, 0326 watchit(sprintf('I attempted to include %s. However no such electrode was found in GRP variable.', ... 0327 include_chans{x})); 0328 end 0329 end 0330 use_chans=use_chans(1:ct); 0331 else 0332 use_chans=1:n_chan; 0333 end 0334 0335 0336 %% Find time points 0337 if isempty(p.Results.time_wind), 0338 time_wind=[0 GRP.time_pts(end)]; %default time window 0339 else 0340 time_wind=p.Results.time_wind; 0341 end 0342 time_wind=sort(time_wind,2); %first make sure earlier of each pair of time points is first 0343 time_wind=sort(time_wind,1); %next sort time windows from earliest to latest onset 0344 n_wind=size(time_wind,1); 0345 if ~isempty(p.Results.time_block_dur), 0346 tbd=p.Results.time_block_dur; 0347 n_block=length(tbd); 0348 if (n_block>1) && (n_block~=n_wind), 0349 error('If you specify more than one time block duration you need to provide exactly one duration for every time window (in this case %d).',n_wind); 0350 elseif (n_block==1) && (n_wind>1), 0351 tbd=repmat(tbd,1,n_wind); 0352 end 0353 0354 t_winds=[]; 0355 one_step=1000/GRP.srate; 0356 for a=1:n_wind, 0357 wind_strt=time_wind(a,1); 0358 wind_stop=time_wind(a,2); 0359 new_winds=[wind_strt:tbd(a):wind_stop]; 0360 for b=1:(length(new_winds)-1), 0361 t_winds=[t_winds; new_winds(b) new_winds(b+1)-one_step]; 0362 end 0363 end 0364 time_wind=t_winds; 0365 clear t_winds; 0366 n_wind=size(time_wind,1); 0367 end 0368 0369 if mean_wind, 0370 use_tpts=cell(1,n_wind); 0371 else 0372 use_tpts=[]; 0373 end 0374 for a=1:n_wind, 0375 VerbReport(sprintf('Time Window #%d:',a),1,VERBLEVEL); 0376 VerbReport(sprintf('Attempting to use time boundaries of %d to %d ms for hypothesis test.',time_wind(a,1),time_wind(a,2)), ... 0377 1,VERBLEVEL); 0378 start_tpt=find_tpt(time_wind(a,1),GRP.time_pts); 0379 end_tpt=find_tpt(time_wind(a,2),GRP.time_pts); 0380 if mean_wind, 0381 use_tpts{a}=[start_tpt:end_tpt]; 0382 else 0383 use_tpts=[use_tpts [start_tpt:end_tpt]]; 0384 end 0385 %replace desired time points with closest matches 0386 time_wind(a,1)=GRP.time_pts(start_tpt); 0387 time_wind(a,2)=GRP.time_pts(end_tpt); 0388 VerbReport(sprintf('Exact window boundaries are %d to %d ms (that''s from time point %d to %d).', ... 0389 time_wind(a,1),time_wind(a,2),start_tpt,end_tpt),1,VERBLEVEL); 0390 end 0391 if ~mean_wind, 0392 use_tpts=unique(use_tpts); %sorts time points and gets rid of any redundant time points 0393 end 0394 0395 0396 %% Compile data from two groups of participants 0397 %pre-allocate memory 0398 grp_ct=0; 0399 for grp=[grpA grpB], 0400 grp_ct=grp_ct+1; 0401 if grp_ct==1, 0402 %Group A's bin 0403 ur_bin=GRP.bin_info(bin).source_binA; 0404 else 0405 %Group B's bin 0406 ur_bin=GRP.bin_info(bin).source_binB; 0407 end 0408 0409 %Load GND variable for this group 0410 load(GRP.GND_fnames{grp},'-MAT'); 0411 VerbReport(sprintf('Loading individual participant ERPs from Bin %d (%s) of GND variable from %s.', ... 0412 ur_bin,GND.bin_info(ur_bin).bindesc,GRP.GND_fnames{grp}),2,VERBLEVEL); 0413 0414 %Check to make sure time points are still compatible across GND and GRP 0415 %variables 0416 if ~isequal(GND.time_pts,GRP.time_pts) 0417 error('The time points in the GND variable from file %s are NOT the same as those in your GRP variable.', ... 0418 GRP.GND_fnames{grp}); 0419 end 0420 0421 %Derive channel indexes for this GND variable in case GND.chanlocs differs from 0422 %GRP.chanlocs (in order or in identity of channels) 0423 n_chanGND=length(GND.chanlocs); 0424 use_chansGND=[]; 0425 for a=use_chans, 0426 found=0; 0427 for b=1:n_chanGND, 0428 if strcmpi(GRP.chanlocs(a).labels,GND.chanlocs(b).labels), 0429 found=1; 0430 use_chansGND=[use_chansGND b]; 0431 break; 0432 end 0433 end 0434 if ~found, 0435 error('GND variable in file %s is missing channel %s.',GRP.GND_fnames{grp},GRP.chanlocs(a).labels); 0436 end 0437 end 0438 0439 %Use only subs with data in relevant bin(s) 0440 use_subs=find(GND.indiv_bin_ct(:,ur_bin)); 0441 n_sub=length(use_subs); 0442 0443 if mean_wind, 0444 %Take mean amplitude in time blocks and then test 0445 if grp_ct==1, 0446 %Group A's ERPs 0447 n_subA=n_sub; 0448 use_subsA=use_subs; 0449 erpsA=zeros(length(use_chans),n_wind,n_sub); 0450 for a=1:n_wind, 0451 for sub=1:n_sub, 0452 erpsA(:,a,sub)=mean(GND.indiv_erps(use_chansGND,use_tpts{a},ur_bin,use_subs(sub)),2); 0453 end 0454 end 0455 else 0456 %Group B's ERPs 0457 n_subB=n_sub; 0458 use_subsB=use_subs; 0459 erpsB=zeros(length(use_chans),n_wind,n_sub); 0460 for a=1:n_wind, 0461 for sub=1:n_sub, 0462 erpsB(:,a,sub)=mean(GND.indiv_erps(use_chansGND,use_tpts{a},ur_bin,use_subs(sub)),2); 0463 end 0464 end 0465 end 0466 else 0467 %Use every single time point in time window(s) 0468 if grp_ct==1, 0469 %Group A's ERPs 0470 n_subA=n_sub; 0471 use_subsA=use_subs; 0472 n_use_tpts=length(use_tpts); 0473 erpsA=zeros(length(use_chans),n_use_tpts,n_sub); 0474 for sub=1:n_sub, 0475 erpsA(:,:,sub)=GND.indiv_erps(use_chansGND,use_tpts,ur_bin,use_subs(sub)); 0476 end 0477 else 0478 %Group B's ERPs 0479 n_subB=n_sub; 0480 use_subsB=use_subs; 0481 n_use_tpts=length(use_tpts); 0482 erpsB=zeros(length(use_chans),n_use_tpts,n_sub); 0483 for sub=1:n_sub, 0484 erpsB(:,:,sub)=GND.indiv_erps(use_chansGND,use_tpts,ur_bin,use_subs(sub)); 0485 end 0486 end 0487 end 0488 end 0489 df=n_subA+n_subB-2; 0490 0491 %% Report tail of test & alpha levels 0492 VerbReport(sprintf('Testing null hypothesis that the grand average ERPs in GRP variable''s Bin %d (%s) have a mean of %f microvolts.',bin, ... 0493 GRP.bin_info(bin).bindesc,p.Results.null_mean),1,VERBLEVEL); 0494 if p.Results.tail==0 0495 VerbReport(sprintf('Alternative hypothesis is that the ERPs differ from %f (i.e., two-tailed test).',p.Results.null_mean), ... 0496 1,VERBLEVEL); 0497 elseif p.Results.tail<0, 0498 VerbReport(sprintf('Alternative hypothesis is that the ERPs are less than %f (i.e., lower-tailed test).',p.Results.null_mean), ... 0499 1,VERBLEVEL); 0500 else 0501 VerbReport(sprintf('Alternative hypothesis is that the ERPs are greater than %f (i.e., upper-tailed test).',p.Results.null_mean), ... 0502 1,VERBLEVEL); 0503 end 0504 0505 %% Optionally reset random number stream to reproduce a previous test 0506 if isempty(p.Results.reproduce_test), 0507 seed_state=[]; 0508 else 0509 if p.Results.reproduce_test>length(GRP.t_tests), 0510 error('Value of argument ''reproduce_test'' is too high. You only have %d permutation tests stored with this GRP variable.',length(GRP.t_tests)); 0511 else 0512 if isnan(GRP.t_tests(p.Results.reproduce_test).n_perm) 0513 error('t-test set %d is NOT a permutation test. You don''t need to seed the random number generator to reproduce it.', ... 0514 p.Results.reproduce_test); 0515 else 0516 seed_state=GRP.t_tests(p.Results.reproduce_test).seed_state; 0517 end 0518 end 0519 end 0520 0521 %Compute the permutation test 0522 if strcmpi(GRP.bin_info(bin).op,'(A+B)/n'), 0523 %one sample t-test 0524 VerbReport('Performing one sample/repeated measures t-tests.',1,VERBLEVEL); 0525 erpsAB=zeros(length(use_chans),size(erpsA,2),n_subA+n_subB); 0526 erpsAB(:,:,1:n_subA)=erpsA-p.Results.null_mean; 0527 erpsAB(:,:,(n_subA+1):(n_subA+n_subB))=erpsB-p.Results.null_mean; 0528 [prm_pval, data_t, crit_t, seed_state, est_alpha]=mxt_perm1(erpsAB,p.Results.n_perm,p.Results.alpha,p.Results.tail,VERBLEVEL,seed_state); 0529 else 0530 %independent samples t-test 0531 VerbReport('Performing independent samples t-tests.',1,VERBLEVEL); 0532 [prm_pval, data_t, crit_t, seed_state, est_alpha]=mxt_perm2(erpsA-p.Results.null_mean,erpsB,p.Results.n_perm,p.Results.alpha,p.Results.tail,VERBLEVEL,seed_state); 0533 end 0534 0535 %Command line summary of results 0536 VerbReport(['Critical t-score(s):' num2str(crit_t)],1,VERBLEVEL); 0537 if p.Results.tail 0538 %one-tailed test 0539 tw_alpha=1-cdf('t',max(abs(crit_t)),df); 0540 else 0541 %two-tailed test 0542 tw_alpha=2-cdf('t',max(crit_t),df)-cdf('t',-min(crit_t),df); 0543 end 0544 VerbReport(sprintf('That corresponds to a test-wise alpha level of %f.',tw_alpha),1,VERBLEVEL); 0545 VerbReport(sprintf('Bonferroni test-wise alpha would be %f.',p.Results.alpha/(size(prm_pval,1)* ... 0546 size(prm_pval,2))),1,VERBLEVEL); 0547 sig_tpts=find(sum(prm_pval<p.Results.alpha)); 0548 if isempty(sig_tpts), 0549 fprintf('ERPs are NOT significantly different from zero (alpha=%f) at any time point/window analyzed.\n', ... 0550 p.Results.alpha); 0551 fprintf('All p-values>=%f\n',min(min(prm_pval))); 0552 else 0553 fprintf('Significant differences from zero (in order of earliest to latest):\n'); 0554 max_sig_p=0; 0555 min_sig_p=2; 0556 for t=sig_tpts, 0557 if mean_wind 0558 %time windows instead of time points 0559 fprintf('%d to %d ms, electrode(s): ',GRP.time_pts(use_tpts{t}(1)), ... 0560 GRP.time_pts(use_tpts{t}(end))); 0561 else 0562 fprintf('%d ms, electrode(s): ',GRP.time_pts(use_tpts(t))); 0563 end 0564 sig_elec=find(prm_pval(:,t)<p.Results.alpha); 0565 ct=0; 0566 for c=sig_elec', 0567 ct=ct+1; 0568 if prm_pval(c,t)>max_sig_p, 0569 max_sig_p=prm_pval(c,t); 0570 end 0571 if prm_pval(c,t)<min_sig_p, 0572 min_sig_p=prm_pval(c,t); 0573 end 0574 0575 if ct==length(sig_elec), 0576 fprintf('%s.\n',GRP.chanlocs(use_chans(c)).labels); 0577 else 0578 fprintf('%s, ',GRP.chanlocs(use_chans(c)).labels); 0579 end 0580 end 0581 end 0582 fprintf('All significant corrected p-values are between %f and %f\n',max_sig_p,min_sig_p); 0583 end 0584 0585 0586 %Add permutation results to GRP struct 0587 n_t_tests=length(GRP.t_tests); 0588 neo_test=n_t_tests+1; 0589 GRP.t_tests(neo_test).bin=bin; 0590 GRP.t_tests(neo_test).time_wind=time_wind; 0591 GRP.t_tests(neo_test).used_tpt_ids=use_tpts; 0592 n_use_chans=length(use_chans); 0593 include_chans=cell(1,n_use_chans); 0594 for a=1:n_use_chans, 0595 include_chans{a}=GRP.chanlocs(use_chans(a)).labels; 0596 end 0597 GRP.t_tests(neo_test).include_chans=include_chans; 0598 GRP.t_tests(neo_test).used_chan_ids=use_chans; 0599 GRP.t_tests(neo_test).mult_comp_method='tmax perm test'; 0600 GRP.t_tests(neo_test).n_perm=p.Results.n_perm; 0601 GRP.t_tests(neo_test).desired_alphaORq=p.Results.alpha; 0602 GRP.t_tests(neo_test).estimated_alpha=est_alpha; 0603 GRP.t_tests(neo_test).null_mean=p.Results.null_mean; 0604 if mean_wind, 0605 GRP.t_tests(neo_test).data_t=data_t; 0606 GRP.t_tests(neo_test).mean_wind='yes'; 0607 else 0608 GRP.t_tests(neo_test).data_t='See GRP.grands_t'; 0609 GRP.t_tests(neo_test).mean_wind='no'; 0610 end 0611 GRP.t_tests(neo_test).crit_t=crit_t; 0612 GRP.t_tests(neo_test).df=df; 0613 GRP.t_tests(neo_test).adj_pval=prm_pval; 0614 GRP.t_tests(neo_test).fdr_rej=NaN; 0615 GRP.t_tests(neo_test).seed_state=seed_state; 0616 GRP.t_tests(neo_test).clust_info=NaN; 0617 GRP.t_tests(neo_test).chan_hood=NaN; 0618 0619 if strcmpi(p.Results.plot_raster,'yes'), 0620 sig_raster(GRP,neo_test,'verblevel',0,'use_color','rgb'); 0621 end 0622 0623 if mean_wind, 0624 if strcmpi(p.Results.plot_mn_topo,'yes') || isempty(p.Results.plot_mn_topo), 0625 sig_topo(GRP,neo_test,'units','t','verblevel',0); %t-score topographies 0626 sig_topo(GRP,neo_test,'units','uV','verblevel',0); %microvolt topographies 0627 end 0628 else 0629 %plot_pvals(GRP,p.Results.alpha,use_chans,use_tpts,prm_pval,bin); 0630 %%improve plot_pvals and make an option ?? 0631 if strcmpi(p.Results.plot_gui,'yes'), 0632 gui_erp('initialize','GNDorGRP',GRP,'t_test',neo_test,'stat','t', ... 0633 'verblevel',1); 0634 end 0635 end 0636 0637 if ~isempty(p.Results.output_file) 0638 [fid msg]=fopen(p.Results.output_file,'w'); 0639 if fid==-1, 0640 error('Cannot create file %s for writing. According to fopen.m: %s.', ... 0641 p.Results.file,msg); 0642 else 0643 %Write header column of times 0644 % Leave first column blank for channel labels 0645 if mean_wind, 0646 for t=1:n_wind 0647 fprintf(fid,' %d-%d',GRP.time_pts(use_tpts{t}(1)), ... 0648 GRP.time_pts(use_tpts{t}(end))); 0649 end 0650 0651 %write a couple spaces and then write header for t-scores 0652 fprintf(fid,' '); 0653 for t=1:n_wind 0654 fprintf(fid,' %d-%d',GRP.time_pts(use_tpts{t}(1)), ... 0655 GRP.time_pts(use_tpts{t}(end))); 0656 end 0657 else 0658 for t=use_tpts, 0659 fprintf(fid,' %d',GRP.time_pts(t)); 0660 end 0661 end 0662 fprintf(fid,' Milliseconds\n'); 0663 0664 % Write channel labels and p-values 0665 chan_ct=0; 0666 for c=use_chans, 0667 chan_ct=chan_ct+1; 0668 fprintf(fid,'%s',GRP.chanlocs(c).labels); 0669 for t=1:length(use_tpts), 0670 fprintf(fid,' %f',prm_pval(chan_ct,t)); 0671 end 0672 fprintf(fid,' p-value'); 0673 0674 if mean_wind, 0675 %write a couple spaces and then write t-scores if mean amp 0676 %in time windows used 0677 fprintf(fid,' '); 0678 for t=1:n_wind 0679 fprintf(fid,' %f',data_t(chan_ct,t)); 0680 end 0681 fprintf(fid,' t-score \n'); 0682 else 0683 fprintf(fid,'\n'); 0684 end 0685 end 0686 0687 % Write permutation test details 0688 fprintf(fid,'Experiment: %s\n',GRP.exp_desc); 0689 fprintf(fid,'Test_of_null_hypothesis_that_Bin_%d_equals: %f\n',bin,p.Results.null_mean); 0690 fprintf(fid,'#_of_time_windows: %d\n',n_wind); 0691 fprintf(fid,'#_of_permutations: %d\n',p.Results.n_perm); 0692 fprintf(fid,'Tail_of_test: '); 0693 if ~p.Results.tail, 0694 fprintf(fid,'Two_tailed\n'); 0695 fprintf(fid,'Critical_t_scores: %f %f\n',crit_t(1),crit_t(2)); 0696 elseif p.Results.tail>0 0697 fprintf(fid,'Upper_tailed\n'); 0698 fprintf(fid,'Critical_t_score: %f\n',crit_t(1)); 0699 else 0700 fprintf(fid,'Lower_tailed\n'); 0701 fprintf(fid,'Critical_t_score: %f\n',crit_t(1)); 0702 end 0703 fprintf(fid,'Degrees_of_freedom: %d\n',df); 0704 fprintf(fid,'Alpha_level: %f\n',p.Results.alpha); 0705 0706 for grp=[grpA grpB], 0707 % # of participants and filenames 0708 if grp==grpA, 0709 use_subs=use_subsA; 0710 fprintf(fid,'GND_fname_GroupA: %s\n',GRP.GND_fnames{grp}); 0711 fprintf(fid,'#_of_participants_GroupA: %d\n',n_subA); 0712 fprintf(fid,'Participant_names_GroupA: \n'); 0713 else 0714 use_subs=use_subsB; 0715 fprintf(fid,'GND_fname_GroupB: %s\n',GRP.GND_fnames{grp}); 0716 fprintf(fid,'#_of_participants_GroupB: %d\n',n_subB); 0717 fprintf(fid,'Participant_names_GroupB: \n'); 0718 end 0719 for s=1:length(use_subs), 0720 fprintf(fid,'%s\n',GRP.indiv_subnames{grp}{use_subs(s)}); 0721 end 0722 end 0723 end 0724 fclose(fid); 0725 end 0726 0727 if ~strcmpi(p.Results.save_GRP,'no'), 0728 GRP=save_matmk(GRP); 0729 end 0730 0731 % 0732 %% %%%%%%%%%%%%%%%%%%%%% function str2bool() %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0733 % 0734 function bool=str2bool(str) 0735 %function bool=str2bool(str) 0736 0737 if ischar(str), 0738 if strcmpi(str,'yes') || strcmpi(str,'y') 0739 bool=1; 0740 else 0741 bool=0; 0742 end 0743 else 0744 bool=str; 0745 end 0746