tmax_specGRP() - Tests the null hypothesis that the grand average voltage
of a between-subject difference in power spectra is mu
(usually 0) using a permutation test based on the
tmax statistic (e.g., Hemmelmann et al., 2004). This function
requires individual subject power spectra to be stored in a
"specGRP" structure and outputs the test results in a number
of graphical and text formats. For analogous within-subject
comparisons use the function tmax_specGND.m. Note, when applied
to a bin that is the mean power spectra across two groups
(i.e., NOT a difference groups), a one-sample test is executed.
Usage:
>> [specGRP, prm_pval, data_t, crit_t]=tmax_specGRP(specGRP_or_fname,bin,varargin)
Required Inputs:
specGRP_or_fname - A specGRP structure variable or the filename of a specGRP
structure that has been saved to disk. A specGRP variable
is based on specGND variables. To create a specGRP variable from
specGND variables use specGNDs2specGRP.m. See MATLABmk
documentation for detailed information about the
format of a specGRP 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_spec.m to see what
bins are stored in a specGRP variable. Use the function
bin_op_specGRP.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 power difference is greater
than mu). "0" means two-tailed (i.e., alternative hypothesis
is that the power difference 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}
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 power 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 (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_spec.m to see
the channel labels stored in the specGRP 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 specGRP
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
ERPs in each bin via various stats (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 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 is produced. More
specifically, two figures are produced: one showing the
topographies in 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_specGRP - ['yes' or 'no'] If 'yes', the specGRP 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 specGRP variable to reproduce. For example, if
'reproduce_test' equals 2, the second t-test
stored in the specGRP variable (i.e., specGRP.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:
specGRP - specGRP structure variable. This is the same as
the input specGRP variable with one addition: the
field specGRP.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.
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 bin that is the difference in spectral power between two groups
to a specGRP variable, use the function "bin_op_specGRP".
-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, 12/20100001 % tmax_specGRP() - Tests the null hypothesis that the grand average voltage 0002 % of a between-subject difference in power spectra is mu 0003 % (usually 0) using a permutation test based on the 0004 % tmax statistic (e.g., Hemmelmann et al., 2004). This function 0005 % requires individual subject power spectra to be stored in a 0006 % "specGRP" structure and outputs the test results in a number 0007 % of graphical and text formats. For analogous within-subject 0008 % comparisons use the function tmax_specGND.m. Note, when applied 0009 % to a bin that is the mean power spectra across two groups 0010 % (i.e., NOT a difference groups), a one-sample test is executed. 0011 % 0012 % 0013 % Usage: 0014 % >> [specGRP, prm_pval, data_t, crit_t]=tmax_specGRP(specGRP_or_fname,bin,varargin) 0015 % 0016 % Required Inputs: 0017 % specGRP_or_fname - A specGRP structure variable or the filename of a specGRP 0018 % structure that has been saved to disk. A specGRP variable 0019 % is based on specGND variables. To create a specGRP variable from 0020 % specGND variables use specGNDs2specGRP.m. See MATLABmk 0021 % documentation for detailed information about the 0022 % format of a specGRP 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_spec.m to see what 0027 % bins are stored in a specGRP variable. Use the function 0028 % bin_op_specGRP.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 power difference is greater 0035 % than mu). "0" means two-tailed (i.e., alternative hypothesis 0036 % is that the power difference 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 % freq_band - 2D matrix of pairs of frequency values specifying the beginning 0049 % and end of one or more frequency bands in Hz (e.g., 0050 % [3 8; 8 12]). Every single frequency in 0051 % the frequency band will be individually tested (i.e., 0052 % maximal temporal resolution) if mean_band option is 0053 % NOT used. Note, boundaries of frequency band(s) may not 0054 % exactly correspond to desired frequency band boundaries 0055 % because of temporal digitization (i.e., you only have 0056 % samples every so many Hz). {default: 0 Hz to the 0057 % Nyquist frequency} 0058 % freq_block_dur- [integers] A number or numbers (in Hz) specifying a 0059 % duration of frequency blocks in which the frequency 0060 % windows specified by freq_band will be divided. For 0061 % example, if freq_band=[0 40] and freq_block_dur=10, 0062 % the 0 to 40 Hz frequency band will be sub-divided into 0063 % 10 Hz bands (i.e., freq_band will equal something like 0064 % [0 9; 10 19; 20 29; 30 39]--depending on your frequency 0065 % resolution). This is an easy way to break up larger 0066 % frequency bands of interest into smaller bands for 0067 % mean band analysis (if the 'mean_band' option is set to 0068 % 'yes'). If you specify multiple frequency bands with 0069 % 'freq_band', you can break them up using durations of 0070 % different lengths. For example, if freq_band=[0 10; 0071 % 40 90] and freq_block_dur=[5 10], the first frequency 0072 % band (0 to 10 Hz) will be broken up into 5 Hz 0073 % bands and the second band (40 to 90 Hz) will be 0074 % broken up into 10 Hz bands. {default: not used} 0075 % mean_band - ['yes' or 'no'] If 'yes', the permutation test will be 0076 % performed on the mean power within each frequency band 0077 % specified by freq_band. This sacrifices frequency 0078 % resolution to increase test power by reducing the number 0079 % of comparisons. If 'no', every single frequency within 0080 % freq_band's frequency bands will be tested individually. 0081 % {default: 'no'} 0082 % null_mean - [number] The mean of the null hypothesis (in units of 0083 % 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_spec.m to see 0088 % the channel labels stored in the specGRP 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_spec.m to see the channel labels stored in the specGRP 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_pow.m. The GUI vizualizes the grand average 0111 % ERPs in each bin via various stats (e.g., t-scores), shows 0112 % topographies at individual frequencies, and illustrates 0113 % which electrodes significantly differ from the null 0114 % hypothesis. This option does not work if mean_band 0115 % option is set to 'yes.' This GUI can be reproduced using 0116 % the function gui_pow.m. {default: 'yes'} 0117 % plot_raster - ['yes' or 'no'] If 'yes', a two-dimensional (frequency 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 frequency 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 frequency band is produced. More 0129 % specifically, two figures are produced: one showing the 0130 % topographies in decibel power 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_band 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_band 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 frequency/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: none} 0145 % save_specGRP - ['yes' or 'no'] If 'yes', the specGRP variable will be 0146 % saved to disk after the permutation test is completed 0147 % and added to it. User will first be prompted to verify 0148 % file name and path. {default: 'yes'} 0149 % reproduce_test- [integer] The number of the permutation test stored in 0150 % the specGRP variable to reproduce. For example, if 0151 % 'reproduce_test' equals 2, the second t-test 0152 % stored in the specGRP variable (i.e., specGRP.t_tests(2)) will 0153 % be reproduced. Reproduction is accomplished by setting 0154 % the random number generator used in the permutation test 0155 % to the same initial state it was in when the permutation 0156 % test was first applied. 0157 % 0158 % Outputs: 0159 % specGRP - specGRP structure variable. This is the same as 0160 % the input specGRP variable with one addition: the 0161 % field specGRP.t_tests will contain the results of the 0162 % permutation test and the test parameters. 0163 % prm_pval - A two-dimensional matrix (channel x frequency) of the 0164 % p-values of each comparison. 0165 % data_t - A two-dimensional matrix (channel x frequency) of the 0166 % t-scores of each comparison. 0167 % crit_t - The critical t-score(s) for the test. Any t-scores 0168 % that are more extreme than the critical t-score(s) 0169 % significantly deviate from 0. 0170 % 0171 % Note also that a great deal of information about the test is displayed 0172 % in the MATLAB command window. You can easiy record of all this 0173 % information to a text file using the MATLAB command "diary." 0174 % 0175 % Global Variables: 0176 % VERBLEVEL = MATLABmk level of verbosity (i.e., tells functions how much 0177 % to report about what they're doing during runtime) set by 0178 % the optional function argument 'verblevel' 0179 % 0180 % Notes: 0181 % -To add a bin that is the difference in spectral power between two groups 0182 % to a specGRP variable, use the function "bin_op_specGRP". 0183 % 0184 % -Unlike a parametric test (e.g., an ANOVA), a discrete set of p-values 0185 % are possible (at most the number of possible permutations). Since the 0186 % number of possible permutations grows rapdily with the number of 0187 % participants, this is only issue for small sample sizes (e.g., 3 0188 % participants in each group). When you have such a small sample size, the 0189 % limited number of p-values may make the test overly conservative (e.g., 0190 % you might be forced to use an alpha level of .0286 since it is the biggest 0191 % possible alpha level less than .05). 0192 % 0193 % References: 0194 % Manly, B.F.J. (1997) Randomization, bootstrap, and Monte Carlo methods in 0195 % biology. 2nd ed. Chapmn and Hall, London. 0196 % 0197 % Hemmelmann, et al. (2004) Multivariate tests for the evaluation of 0198 % high-dimensional EEG data. Journal of Neuroscience Methods. 0199 % 0200 % 0201 % Author: 0202 % David Groppe 0203 % Kutaslab, 12/2010 0204 0205 %%%%%%%%%%%%%%%% REVISION LOG %%%%%%%%%%%%%%%%% 0206 % ?/?/2010-?? 0207 % 0208 0209 function [specGRP, prm_pval, data_t, crit_t]=tmax_specGRP(specGRP_or_fname,bin,varargin) 0210 0211 global VERBLEVEL; 0212 0213 p=inputParser; 0214 p.addRequired('specGRP_or_fname',@(x) ischar(x) || isstruct(x)); 0215 p.addRequired('bin',@(x) isnumeric(x) && (length(x)==1) && (x>0)); 0216 p.addParamValue('tail',0,@(x) isnumeric(x) && (length(x)==1)); 0217 p.addParamValue('alpha',0.05,@(x) isnumeric(x) && (x>0) && (x<1)); 0218 p.addParamValue('freq_band',[],@(x) isnumeric(x) && (size(x,2)==2)); 0219 p.addParamValue('freq_block_dur',[],@isnumeric); 0220 p.addParamValue('mean_band','no',@(x) ischar(x) && (strcmpi(x,'yes') || strcmpi(x,'no'))); 0221 p.addParamValue('n_perm',2500,@(x) isnumeric(x) && (length(x)==1)); 0222 p.addParamValue('null_mean',0,@(x) isnumeric(x) && (length(x)==1)); 0223 p.addParamValue('verblevel',[],@(x) isnumeric(x) && (length(x)==1)); 0224 p.addParamValue('exclude_chans',[],@(x) ischar(x) || iscell(x)); 0225 p.addParamValue('include_chans',[],@(x) ischar(x) || iscell(x)); 0226 p.addParamValue('plot_gui','yes',@(x) ischar(x) && (strcmpi(x,'yes') || strcmpi(x,'no'))); 0227 p.addParamValue('plot_raster','yes',@(x) ischar(x) && (strcmpi(x,'yes') || strcmpi(x,'no'))); 0228 p.addParamValue('plot_mn_topo',[],@(x) ischar(x) && (strcmpi(x,'yes') || strcmpi(x,'no'))); 0229 p.addParamValue('output_file',[],@ischar); 0230 p.addParamValue('reproduce_test',[],@(x) isnumeric(x) && (length(x)==1)); 0231 p.addParamValue('save_specGRP','yes',@(x) ischar(x) && (strcmpi(x,'yes') || strcmpi(x,'no'))); 0232 0233 p.parse(specGRP_or_fname,bin,varargin{:}); 0234 0235 0236 if isempty(p.Results.verblevel), 0237 if isempty(VERBLEVEL), 0238 VERBLEVEL=2; 0239 end 0240 else 0241 VERBLEVEL=p.Results.verblevel; 0242 end 0243 0244 mean_band=str2bool(p.Results.mean_band); 0245 0246 %Load specGRP struct 0247 if ischar(specGRP_or_fname), 0248 fprintf('Loading specGRP struct from file %s.\n',specGRP_or_fname); 0249 load(specGRP_or_fname,'-MAT'); 0250 else 0251 specGRP=specGRP_or_fname; 0252 clear specGRP_or_fname; 0253 end 0254 [n_chan, n_pt, n_bin]=size(specGRP.grands_pow_dB); 0255 n_group=length(specGRP.specGND_fnames); 0256 VerbReport(sprintf('Experiment: %s',specGRP.exp_desc),2,VERBLEVEL); 0257 0258 if (bin>n_bin), 0259 error('There is no Bin %d in this specGRP variable.',bin); 0260 end 0261 grpA=specGRP.bin_info(bin).groupA; 0262 grpB=specGRP.bin_info(bin).groupB; 0263 0264 0265 %% Figure out which channels to ignore if any 0266 %But first make sure exclude & include options were not both used. 0267 if ~isempty(p.Results.include_chans) && ~isempty(p.Results.exclude_chans) 0268 error('You cannot use BOTH ''include_chans'' and ''exclude_chans'' options.'); 0269 end 0270 if ischar(p.Results.exclude_chans), 0271 exclude_chans{1}=p.Results.exclude_chans; 0272 elseif isempty(p.Results.exclude_chans) 0273 exclude_chans=[]; 0274 else 0275 exclude_chans=p.Results.exclude_chans; 0276 end 0277 if ischar(p.Results.include_chans), 0278 include_chans{1}=p.Results.include_chans; 0279 elseif isempty(p.Results.include_chans) 0280 include_chans=[]; 0281 else 0282 include_chans=p.Results.include_chans; 0283 end 0284 0285 % exclude and include chan options 0286 if ~isempty(exclude_chans), 0287 ignore_chans=zeros(1,length(exclude_chans)); %preallocate mem 0288 ct=0; 0289 for x=1:length(exclude_chans), 0290 found=0; 0291 for c=1:n_chan, 0292 if strcmpi(exclude_chans{x},specGRP.chanlocs(c).labels), 0293 found=1; 0294 ct=ct+1; 0295 ignore_chans(ct)=c; 0296 end 0297 end 0298 if ~found, 0299 watchit(sprintf('I attempted to exclude %s. However no such electrode was found in specGRP variable.', ... 0300 exclude_chans{x})); 0301 end 0302 end 0303 ignore_chans=ignore_chans(1:ct); 0304 use_chans=setdiff(1:n_chan,ignore_chans); 0305 elseif ~isempty(include_chans), 0306 use_chans=zeros(1,length(include_chans)); %preallocate mem 0307 ct=0; 0308 for x=1:length(include_chans), 0309 found=0; 0310 for c=1:n_chan, 0311 if strcmpi(include_chans{x},specGRP.chanlocs(c).labels), 0312 found=1; 0313 ct=ct+1; 0314 use_chans(ct)=c; 0315 end 0316 end 0317 if ~found, 0318 watchit(sprintf('I attempted to include %s. However no such electrode was found in specGRP variable.', ... 0319 include_chans{x})); 0320 end 0321 end 0322 use_chans=use_chans(1:ct); 0323 else 0324 use_chans=1:n_chan; 0325 end 0326 0327 0328 %% Find frequency ids 0329 if isempty(p.Results.freq_band), 0330 freq_band=[0 specGRP.freqs(end)]; %default frequency band 0331 else 0332 freq_band=p.Results.freq_band; 0333 end 0334 freq_band=sort(freq_band,2); %first make sure earlier of each pair of frequencies is first 0335 freq_band=sort(freq_band,1); %next sort frequency bands from earliest to latest onset 0336 n_band=size(freq_band,1); 0337 if ~isempty(p.Results.freq_block_dur), 0338 fbd=p.Results.freq_block_dur; 0339 n_block=length(fbd); 0340 if (n_block>1) && (n_block~=n_band), 0341 error('If you specify more than one frequency block duration you need to provide exactly one duration for every frequency (in this case %d).',n_band); 0342 elseif (n_block==1) && (n_band>1), 0343 fbd=repmat(fbd,1,n_band); 0344 end 0345 0346 f_bands=[]; 0347 one_step=1000/specGRP.srate; 0348 for a=1:n_band, 0349 band_strt=freq_band(a,1); 0350 band_stop=freq_band(a,2); 0351 new_bands=[band_strt:fbd(a):band_stop]; 0352 for b=1:(length(new_bands)-1), 0353 f_bands=[f_bands; new_bands(b) new_bands(b+1)-one_step]; 0354 end 0355 end 0356 freq_band=f_bands; 0357 clear f_bands; 0358 n_band=size(freq_band,1); 0359 end 0360 0361 if mean_band, 0362 use_frq_ids=cell(1,n_band); 0363 else 0364 use_frq_ids=[]; 0365 end 0366 for a=1:n_band, 0367 VerbReport(sprintf('Frequency Band #%d:',a),1,VERBLEVEL); 0368 VerbReport(sprintf('Attempting to use frequency boundaries of %d to %d Hz for hypothesis test.',freq_band(a,1),freq_band(a,2)), ... 0369 1,VERBLEVEL); 0370 start_freq=find_tpt(freq_band(a,1),specGRP.freqs); 0371 end_freq=find_tpt(freq_band(a,2),specGRP.freqs); 0372 if mean_band, 0373 use_frq_ids{a}=[start_freq:end_freq]; 0374 else 0375 use_frq_ids=[use_frq_ids [start_freq:end_freq]]; 0376 end 0377 %replace desired frequencies with closest matches 0378 freq_band(a,1)=specGRP.freqs(start_freq); 0379 freq_band(a,2)=specGRP.freqs(end_freq); 0380 VerbReport(sprintf('Exact frequency boundaries are %d to %d Hz (that''s from frequency number %d to %d).', ... 0381 freq_band(a,1),freq_band(a,2),start_freq,end_freq),1,VERBLEVEL); 0382 end 0383 if ~mean_band, 0384 use_frq_ids=unique(use_frq_ids); %sorts frequency ids and gets rid of any redundant frequencies 0385 end 0386 0387 0388 %% Compile data from two groups of participants 0389 %pre-allocate memory 0390 grp_ct=0; 0391 for grp=[grpA grpB], 0392 grp_ct=grp_ct+1; 0393 if grp_ct==1, 0394 %Group A's bin 0395 ur_bin=specGRP.bin_info(bin).source_binA; 0396 else 0397 %Group B's bin 0398 ur_bin=specGRP.bin_info(bin).source_binB; 0399 end 0400 0401 %Load specGND variable for this group 0402 load(specGRP.specGND_fnames{grp},'-MAT'); 0403 VerbReport(sprintf('Loading individual participant power spectra from Bin %d (%s) of specGND variable from %s.', ... 0404 ur_bin,specGND.bindesc{ur_bin},specGRP.specGND_fnames{grp}),2,VERBLEVEL); 0405 0406 %Check to make sure frequencies are still compatible across specGND and specGRP 0407 %variables 0408 if ~isequal(specGND.freqs,specGRP.freqs) 0409 error('The frequencies in the specGND variable from file %s are NOT the same as those in your specGRP variable.', ... 0410 specGRP.specGND_fnames{grp}); 0411 end 0412 0413 %Derive channel indexes for this specGND variable in case specGND.chanlocs differs from 0414 %specGRP.chanlocs (in order or in identity of channels) 0415 n_chan_specGND=length(specGND.chanlocs); 0416 use_chans_specGND=[]; 0417 for a=use_chans, 0418 found=0; 0419 for b=1:n_chan_specGND, 0420 if strcmpi(specGRP.chanlocs(a).labels,specGND.chanlocs(b).labels), 0421 found=1; 0422 use_chans_specGND=[use_chans_specGND b]; 0423 break; 0424 end 0425 end 0426 if ~found, 0427 error('specGND variable in file %s is missing channel %s.',specGRP.specGND_fnames{grp},specGRP.chanlocs(a).labels); 0428 end 0429 end 0430 0431 %Use only subs with data in relevant bin(s) 0432 use_subs=find(specGND.indiv_bin_ct(:,ur_bin)); 0433 n_sub=length(use_subs); 0434 if mean_band, 0435 %Take mean amplitude in frequency band and then test 0436 if grp_ct==1, 0437 %Group A's Power Spectra 0438 n_subA=n_sub; 0439 use_subsA=use_subs; 0440 powA=zeros(length(use_chans),n_band,n_sub); 0441 for a=1:n_band, 0442 for sub=1:n_sub, 0443 powA(:,a,sub)=mean(specGND.indiv_pow_dB(use_chans_specGND,use_frq_ids{a},ur_bin,use_subs(sub)),2); 0444 end 0445 end 0446 else 0447 %Group B's Power Spectra 0448 n_subB=n_sub; 0449 use_subsB=use_subs; 0450 powB=zeros(length(use_chans),n_band,n_sub); 0451 for a=1:n_band, 0452 for sub=1:n_sub, 0453 powB(:,a,sub)=mean(specGND.indiv_pow_dB(use_chans_specGND,use_frq_ids{a},ur_bin,use_subs(sub)),2); 0454 end 0455 end 0456 end 0457 else 0458 %Use every single frequency in frequency band(s) 0459 if grp_ct==1, 0460 %Group A's Power Spectra 0461 n_subA=n_sub; 0462 use_subsA=use_subs; 0463 n_use_frq_ids=length(use_frq_ids); 0464 powA=zeros(length(use_chans),n_use_frq_ids,n_sub); 0465 for sub=1:n_sub, 0466 powA(:,:,sub)=specGND.indiv_pow_dB(use_chans_specGND,use_frq_ids,ur_bin,use_subs(sub)); 0467 end 0468 else 0469 %Group B's Power Spectra 0470 n_subB=n_sub; 0471 use_subsB=use_subs; 0472 n_use_frq_ids=length(use_frq_ids); 0473 powB=zeros(length(use_chans),n_use_frq_ids,n_sub); 0474 for sub=1:n_sub, 0475 powB(:,:,sub)=specGND.indiv_pow_dB(use_chans_specGND,use_frq_ids,ur_bin,use_subs(sub)); 0476 end 0477 end 0478 end 0479 end 0480 df=n_subA+n_subB-2; 0481 0482 %% Report tail of test & alpha levels 0483 VerbReport(sprintf('Testing null hypothesis that the grand average power spectra in specGRP variable''s Bin %d (%s) have a mean of %f 10*log10((uV^2)/Hz).',bin, ... 0484 specGRP.bin_info(bin).bindesc,p.Results.null_mean),1,VERBLEVEL); 0485 if p.Results.tail==0 0486 VerbReport(sprintf('Alternative hypothesis is that the power spectra differ from %f (i.e., two-tailed test).',p.Results.null_mean), ... 0487 1,VERBLEVEL); 0488 elseif p.Results.tail<0, 0489 VerbReport(sprintf('Alternative hypothesis is that the power spectra are less than %f (i.e., lower-tailed test).',p.Results.null_mean), ... 0490 1,VERBLEVEL); 0491 else 0492 VerbReport(sprintf('Alternative hypothesis is that the power spectra are greater than %f (i.e., upper-tailed test).',p.Results.null_mean), ... 0493 1,VERBLEVEL); 0494 end 0495 0496 %% Optionally reset random number stream to reproduce a previous test 0497 if isempty(p.Results.reproduce_test), 0498 seed_state=[]; 0499 else 0500 if p.Results.reproduce_test>length(specGRP.t_tests), 0501 error('Value of argument ''reproduce_test'' is too high. You only have %d permutation tests stored with this specGRP variable.',length(specGRP.t_tests)); 0502 else 0503 if isnan(specGRP.t_tests(p.Results.reproduce_test).n_perm) 0504 error('t-test set %d is NOT a permutation test. You don''t need to seed the random number generator to reproduce it.', ... 0505 p.Results.reproduce_test); 0506 else 0507 seed_state=specGRP.t_tests(p.Results.reproduce_test).seed_state; 0508 end 0509 end 0510 end 0511 0512 %Compute the permutation test 0513 freq_domain=1; 0514 if strcmpi(specGRP.bin_info(bin).op,'(A+B)/n'), 0515 %one sample t-test 0516 VerbReport('Performing one sample/repeated measures t-tests.',1,VERBLEVEL); 0517 powAB=zeros(length(use_chans),size(powA,2),n_subA+n_subB); 0518 powAB(:,:,1:n_subA)=powA-p.Results.null_mean; 0519 powAB(:,:,(n_subA+1):(n_subA+n_subB))=powB-p.Results.null_mean; 0520 [prm_pval, data_t, crit_t, seed_state, est_alpha]=mxt_perm1(powAB,p.Results.n_perm,p.Results.alpha,p.Results.tail,VERBLEVEL,seed_state,freq_domain); 0521 else 0522 %independent samples t-test 0523 VerbReport('Performing independent samples t-tests.',1,VERBLEVEL); 0524 [prm_pval, data_t, crit_t, seed_state, est_alpha]=mxt_perm2(powA-p.Results.null_mean,powB,p.Results.n_perm,p.Results.alpha,p.Results.tail,VERBLEVEL,seed_state,freq_domain); 0525 end 0526 0527 %Command line summary of results 0528 VerbReport(['Critical t-score(s):' num2str(crit_t)],1,VERBLEVEL); 0529 if p.Results.tail 0530 %one-tailed test 0531 tw_alpha=1-cdf('t',max(abs(crit_t)),df); 0532 else 0533 %two-tailed test 0534 tw_alpha=2-cdf('t',max(crit_t),df)-cdf('t',-min(crit_t),df); 0535 end 0536 VerbReport(sprintf('That corresponds to a test-wise alpha level of %f.',tw_alpha),1,VERBLEVEL); 0537 VerbReport(sprintf('Bonferroni test-wise alpha would be %f.',p.Results.alpha/(size(prm_pval,1)* ... 0538 size(prm_pval,2))),1,VERBLEVEL); 0539 sig_freq_ids=find(sum(prm_pval<p.Results.alpha)); 0540 if isempty(sig_freq_ids), 0541 fprintf('Spectral power is NOT significantly different from zero (alpha=%f) at any frequency/band analyzed.\n', ... 0542 p.Results.alpha); 0543 fprintf('All p-values>=%f\n',min(min(prm_pval))); 0544 else 0545 fprintf('Significant differences from zero (in order of lowest to highest frequency):\n'); 0546 max_sig_p=0; 0547 %frequency bands instead of individual frequencies 0548 for t=sig_freq_ids, 0549 if mean_band 0550 fprintf('%d to %d Hz, electrode(s): ',specGRP.freqs(use_frq_ids{t}(1)), ... 0551 specGRP.freqs(use_frq_ids{t}(end))); 0552 else 0553 fprintf('%d Hz, electrode(s): ',specGRP.freqs(use_frq_ids(t))); 0554 end 0555 sig_elec=find(prm_pval(:,t)<p.Results.alpha); 0556 ct=0; 0557 for c=sig_elec', 0558 ct=ct+1; 0559 if prm_pval(c,t)>max_sig_p, 0560 max_sig_p=prm_pval(c,t); 0561 end 0562 if ct==length(sig_elec), 0563 fprintf('%s.\n',specGRP.chanlocs(use_chans(c)).labels); 0564 else 0565 fprintf('%s, ',specGRP.chanlocs(use_chans(c)).labels); 0566 end 0567 end 0568 end 0569 fprintf('All significant p-values<=%f\n',max_sig_p); 0570 end 0571 0572 0573 %Add permutation results to specGRP struct 0574 n_t_tests=length(specGRP.t_tests); 0575 neo_test=n_t_tests+1; 0576 specGRP.t_tests(neo_test).bin=bin; 0577 specGRP.t_tests(neo_test).freq_band=freq_band; 0578 specGRP.t_tests(neo_test).used_freq_ids=use_frq_ids; 0579 n_use_chans=length(use_chans); 0580 include_chans=cell(1,n_use_chans); 0581 for a=1:n_use_chans, 0582 include_chans{a}=specGRP.chanlocs(use_chans(a)).labels; 0583 end 0584 specGRP.t_tests(neo_test).include_chans=include_chans; 0585 specGRP.t_tests(neo_test).used_chan_ids=use_chans; 0586 specGRP.t_tests(neo_test).mult_comp_method='tmax perm test'; 0587 specGRP.t_tests(neo_test).n_perm=p.Results.n_perm; 0588 specGRP.t_tests(neo_test).desired_alphaORq=p.Results.alpha; 0589 specGRP.t_tests(neo_test).estimated_alpha=est_alpha; 0590 specGRP.t_tests(neo_test).null_mean=p.Results.null_mean; 0591 if mean_band, 0592 specGRP.t_tests(neo_test).data_t=data_t; 0593 specGRP.t_tests(neo_test).mean_band='yes'; 0594 else 0595 specGRP.t_tests(neo_test).data_t='See specGRP.grands_pow_dB_t'; 0596 specGRP.t_tests(neo_test).mean_band='no'; 0597 end 0598 specGRP.t_tests(neo_test).crit_t=crit_t; 0599 specGRP.t_tests(neo_test).df=df; 0600 specGRP.t_tests(neo_test).adj_pval=prm_pval; 0601 specGRP.t_tests(neo_test).fdr_rej=NaN; 0602 specGRP.t_tests(neo_test).seed_state=seed_state; 0603 0604 0605 if strcmpi(p.Results.plot_raster,'yes'), 0606 sig_raster(specGRP,neo_test,'verblevel',0); 0607 end 0608 0609 if mean_band, 0610 if strcmpi(p.Results.plot_mn_topo,'yes') || isempty(p.Results.plot_mn_topo), 0611 sig_topo(specGRP,neo_test,'units','t','verblevel',0); %t-score topographies 0612 sig_topo(specGRP,neo_test,'units','dB','verblevel',0); %10*log10((uV^2)/Hz) topographies 0613 end 0614 else 0615 %plot_pvals(specGRP,p.Results.alpha,use_chans,use_frq_ids,prm_pval,bin); 0616 %%improve plot_pvals and make an option ?? 0617 if strcmpi(p.Results.plot_gui,'yes'), 0618 gui_pow('initialize','specGND_or_specGRP',specGRP,'t_test',neo_test,'stat','t', ... 0619 'verblevel',1); 0620 end 0621 end 0622 0623 if ~isempty(p.Results.output_file) 0624 [fid msg]=fopen(p.Results.output_file,'w'); 0625 if fid==-1, 0626 error('Cannot create file %s for writing. According to fopen.m: %s.', ... 0627 p.Results.file,msg); 0628 else 0629 %Write header column of frequencies 0630 % Leave first column blank for channel labels 0631 if mean_band, 0632 for t=1:n_band 0633 fprintf(fid,' %d-%d',specGRP.freqs(use_frq_ids{t}(1)), ... 0634 specGRP.freqs(use_frq_ids{t}(end))); 0635 end 0636 0637 %write a couple spaces and then write header for t-scores 0638 fprintf(fid,' '); 0639 for t=1:n_band 0640 fprintf(fid,' %d-%d',specGRP.freqs(use_frq_ids{t}(1)), ... 0641 specGRP.freqs(use_frq_ids{t}(end))); 0642 end 0643 else 0644 for t=use_frq_ids, 0645 fprintf(fid,' %d',specGRP.freqs(t)); 0646 end 0647 end 0648 fprintf(fid,' Hz\n'); 0649 0650 % Write channel labels and p-values 0651 chan_ct=0; 0652 for c=use_chans, 0653 chan_ct=chan_ct+1; 0654 fprintf(fid,'%s',specGRP.chanlocs(c).labels); 0655 for t=1:length(use_frq_ids), 0656 fprintf(fid,' %f',prm_pval(chan_ct,t)); 0657 end 0658 fprintf(fid,' p-value'); 0659 0660 if mean_band, 0661 %write a couple spaces and then write t-scores if mean amp 0662 %in frequency bands used 0663 fprintf(fid,' '); 0664 for t=1:n_band 0665 fprintf(fid,' %f',data_t(chan_ct,t)); 0666 end 0667 fprintf(fid,' t-score \n'); 0668 else 0669 fprintf(fid,'\n'); 0670 end 0671 end 0672 0673 % Write permutation test details 0674 fprintf(fid,'Experiment: %s\n',specGRP.exp_desc); 0675 fprintf(fid,'Test_of_null_hypothesis_that_Bin_%d_equals_%f\n',bin,p.Results.null_mean); 0676 fprintf(fid,'#_of_frequency_bands: %d\n',n_band); 0677 fprintf(fid,'#_of_permutations: %d\n',p.Results.n_perm); 0678 fprintf(fid,'Tail_of_test: '); 0679 if ~p.Results.tail, 0680 fprintf(fid,'Two_tailed\n'); 0681 fprintf(fid,'Critical_t_scores: %f %f\n',crit_t(1),crit_t(2)); 0682 elseif p.Results.tail>0 0683 fprintf(fid,'Upper_tailed\n'); 0684 fprintf(fid,'Critical_t_score: %f\n',crit_t(1)); 0685 else 0686 fprintf(fid,'Lower_tailed\n'); 0687 fprintf(fid,'Critical_t_score: %f\n',crit_t(1)); 0688 end 0689 fprintf(fid,'Degrees_of_freedom: %d\n',df); 0690 fprintf(fid,'Alpha_level: %f\n',p.Results.alpha); 0691 0692 for grp=[grpA grpB], 0693 % # of participants and filenames 0694 if grp==grpA, 0695 use_subs=use_subsA; 0696 fprintf(fid,'specGND_fname_GroupA: %s\n',specGRP.specGND_fnames{grp}); 0697 fprintf(fid,'#_of_participants_GroupA: %d\n',n_subA); 0698 fprintf(fid,'Participant_names_GroupA: \n'); 0699 else 0700 use_subs=use_subsB; 0701 fprintf(fid,'specGND_fname_GroupB: %s\n',specGRP.specGND_fnames{grp}); 0702 fprintf(fid,'#_of_participants_GroupB: %d\n',n_subB); 0703 fprintf(fid,'Participant_names_GroupB: \n'); 0704 end 0705 for s=1:length(use_subs), 0706 fprintf(fid,'%s\n',specGRP.indiv_subnames{grp}{use_subs(s)}); 0707 end 0708 end 0709 end 0710 fclose(fid); 0711 end 0712 0713 if ~strcmpi(p.Results.save_specGRP,'no'), 0714 specGRP=save_matmk(specGRP); 0715 end 0716 0717 % 0718 %% %%%%%%%%%%%%%%%%%%%%% function str2bool() %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0719 % 0720 function bool=str2bool(str) 0721 %function bool=str2bool(str) 0722 0723 if ischar(str), 0724 if strcmpi(str,'yes') || strcmpi(str,'y') 0725 bool=1; 0726 else 0727 bool=0; 0728 end 0729 else 0730 bool=str; 0731 end 0732