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