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