tfdr_specGND() - Tests the null hypothesis that the grand average power
of a bin is mu using a one sample t-test with FDR
correction for multiple comparisons. 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 t-test (a Gaussian population)
probably isn't valid for bins that are NOT the
difference between two other bins (unless your sample
size is quite large). For analogous between-subject
comparisons use the function tfdr_specGRP.m.
Usage:
>> [specGND, pval, data_t, crit_t]=tfdr_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}
q - A number between 0 and 1 specifying the family-wise
q level of the test. q is the upper bound on the
expected proportion of rejected null hypotheses that are
false rejections (i.e., the FDR). {default: 0.05}
method - ['bh', 'by', or 'bky'] The procedure used to control
the FDR. 'bh' is the classic Benjamini & Hochberg (1995)
procedure, which is guaranteed to control FDR when the
tests are independent or positively dependent (e.g.,
positively correlated Gaussians). 'by' is a much more
conservative version of 'bh' that always controls FDR
(regardless of the dependency structure of the tests--
Benjamini & Yekutieli, 2001). 'bky' is a "two-stage"
version of 'bh' that is more powerful than 'bh' when a
lot of the null hypotheses are false (Benjamini, Krieger, &
Yekutieli, 2006). 'bky' is guaranteed to control FDR when the
tests are independent and tends to be slightly less
powerful than 'bh' when few or no null hypothese are
false. {default: 'bh'}
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 temporal
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., subject names,
FDR method, 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'}
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.
pval - A two-dimensional matrix (channel x time) of the
p-values of each comparison (no correction for multiple
comparisons).
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".
References:
Benjamini, Y., & Hochberg, Y. (1995). Controlling the false discovery
rate: A practical and powerful approach to multiple testing. Journal of
the Royal Statistical Society. Series B (Methodological), 57(1), 289-300.
Benjamini, Y., Krieger, A. M., & Yekutieli, D. (2006). Adaptive linear
step-up procedures that control the false discovery rate. Biometrika,
93(3), 491-507.
Benjamini, Y., & Yekutieli, D. (2001). The control of the false
discovery rate in multiple testing under dependency. The Annals of
Statistics, 29(4), 1165-1188.
Author:
David Groppe
Kutaslab, 12/20100001 % tfdr_specGND() - Tests the null hypothesis that the grand average power 0002 % of a bin is mu using a one sample t-test with FDR 0003 % correction for multiple comparisons. Note, mu is 0004 % assumed to be 0 by default. This function requires 0005 % 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 t-test (a Gaussian population) 0011 % probably isn't valid for bins that are NOT the 0012 % difference between two other bins (unless your sample 0013 % size is quite large). For analogous between-subject 0014 % comparisons use the function tfdr_specGRP.m. 0015 % 0016 % 0017 % Usage: 0018 % >> [specGND, pval, data_t, crit_t]=tfdr_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 % q - A number between 0 and 1 specifying the family-wise 0045 % q level of the test. q is the upper bound on the 0046 % expected proportion of rejected null hypotheses that are 0047 % false rejections (i.e., the FDR). {default: 0.05} 0048 % method - ['bh', 'by', or 'bky'] The procedure used to control 0049 % the FDR. 'bh' is the classic Benjamini & Hochberg (1995) 0050 % procedure, which is guaranteed to control FDR when the 0051 % tests are independent or positively dependent (e.g., 0052 % positively correlated Gaussians). 'by' is a much more 0053 % conservative version of 'bh' that always controls FDR 0054 % (regardless of the dependency structure of the tests-- 0055 % Benjamini & Yekutieli, 2001). 'bky' is a "two-stage" 0056 % version of 'bh' that is more powerful than 'bh' when a 0057 % lot of the null hypotheses are false (Benjamini, Krieger, & 0058 % Yekutieli, 2006). 'bky' is guaranteed to control FDR when the 0059 % tests are independent and tends to be slightly less 0060 % powerful than 'bh' when few or no null hypothese are 0061 % false. {default: 'bh'} 0062 % freq_band - 2D matrix of pairs of frequency values specifying the beginning 0063 % and end of one or more frequency bands in Hz (e.g., 0064 % [3 8; 8 12]). Every single frequency in 0065 % the frequency band will be individually tested (i.e., 0066 % maximal temporal resolution) if mean_band option is 0067 % NOT used. Note, boundaries of frequency band(s) may not 0068 % exactly correspond to desired frequency band boundaries 0069 % because of temporal digitization (i.e., you only have 0070 % samples every so many Hz). {default: 0 Hz to the 0071 % Nyquist frequency} 0072 % freq_block_dur- [integers] A number or numbers (in Hz) specifying a 0073 % duration of frequency blocks in which the frequency 0074 % windows specified by freq_band will be divided. For 0075 % example, if freq_band=[0 40] and freq_block_dur=10, 0076 % the 0 to 40 Hz frequency band will be sub-divided into 0077 % 10 Hz bands (i.e., freq_band will equal something like 0078 % [0 9; 10 19; 20 29; 30 39]--depending on your frequency 0079 % resolution). This is an easy way to break up larger 0080 % frequency bands of interest into smaller bands for 0081 % mean band analysis (if the mean_band option is set to 0082 % 'yes'). If you specify multiple frequency bands with 0083 % freq_band, you can break them up using durations of 0084 % different lengths. For example, if freq_band=[0 10; 0085 % 40 90] and freq_block_dur=[5 10], the first frequency 0086 % band (0 to 10 Hz) will be broken up into 5 Hz 0087 % bands and the second band (40 to 90 Hz) will be 0088 % broken up into 10 Hz bands. {default: not used} 0089 % mean_band - ['yes' or 'no'] If 'yes', the permutation test will be 0090 % performed on the mean amplitude within each frequency band 0091 % specified by freq_band. This sacrifices temporal 0092 % resolution to increase test power by reducing the number 0093 % of comparisons. If 'no', every single frequency within 0094 % freq_band's frequency bands will be tested individually. 0095 % {default: 'no'} 0096 % null_mean - [number] The mean of the null hypothesis (i.e., mu) in 0097 % units of decibel power (i.e., 10*log10((uV^2)/Hz)). 0098 % {default: 0} 0099 % exclude_chans - A cell array of channel labels to exclude from the 0100 % permutation test (e.g., {'A2','lle','rhe'}). This option 0101 % sacrifices spatial resolution to increase test power by 0102 % reducing the number of comparisons. Use headinfo_spec.m 0103 % to see the channel labels stored in the specGND variable. 0104 % You cannot use both this option and 'include_chans' 0105 % (below). {default: not used, all channels included in test} 0106 % include_chans - A cell array of channel labels to use in the permutation 0107 % test (e.g., {'A2','lle','rhe'}). All other channels will 0108 % be ignored. This option sacrifices spatial resolution to 0109 % increase test power by reducing the number of comparisons. 0110 % Use headinfo_spec.m to see the channel labels stored in the 0111 % specGND variable. You cannot use both this option and 0112 % 'exclude_chans' (above). {default: not used, all channels 0113 % included in test} 0114 % verblevel - An integer specifiying the amount of information you want 0115 % this function to provide about what it is doing during runtime. 0116 % Options are: 0117 % 0 - quiet, only show errors, warnings, and EEGLAB reports 0118 % 1 - stuff anyone should probably know 0119 % 2 - stuff you should know the first time you start working 0120 % with a data set {default value} 0121 % 3 - stuff that might help you debug (show all 0122 % reports) 0123 % plot_gui - ['yes' or 'no'] If 'yes', a GUI is created for 0124 % visualizing the results of the permutation test using the 0125 % function gui_pow.m. The GUI vizualizes the grand average 0126 % power in each bin in various statistics (e.g., t-scores), 0127 % shows topographies at individual frequencies, and illustrates 0128 % which electrodes significantly differ from the null 0129 % hypothesis. This option does not work if the mean_band 0130 % option is set to 'yes.' This GUI can be reproduced using 0131 % the function gui_pow.m. {default: 'yes'} 0132 % plot_raster - ['yes' or 'no'] If 'yes', a two-dimensional (frequency x channel) 0133 % binary "raster" diagram is created to illustrate the 0134 % results of the permutation test. Significant negative and 0135 % positive deviations from the null hypothesis are shown 0136 % as black and white rectangles respectively. Non- 0137 % significant comparisons are shown as gray rectangles. 0138 % Clicking on the rectangles will show you the 0139 % corresponding frequency and channel label for that 0140 % rectangle. This figure can be reproduced with the 0141 % function sig_raster.m. {default: 'yes'} 0142 % plot_mn_topo - ['yes' or 'no'] If 'yes', the topographies of the mean 0143 % voltages/effects in each frequency band are produced. More 0144 % specifically, two figures are produced: one showing the 0145 % topographies as decibel power the other in t-scores. 0146 % Significant/nonsignificant comparisons are shown as 0147 % white/black electrodes. Clicking on electrodes will show the 0148 % electrode's name. This figure can be reproduced with 0149 % the function sig_topo.m. This option has NO effect if 0150 % mean_band option is set to 'no'. {default: 'yes'} 0151 % output_file - A string indicating the name of a space delimited text 0152 % file to produce containing the p-values of all comparisons 0153 % and the details of the test (e.g., subject names, 0154 % FDR method, etc...). If mean_band option is 0155 % set to 'yes,' t-scores of each comparison are also 0156 % included since you cannot derive them from the t-scores 0157 % at each frequency/electrode in a simple way. When 0158 % importing this file into a spreadsheet be sure NOT to count 0159 % consecutive spaces as multiple delimiters. {default: 0160 % none} 0161 % save_specGND - ['yes' or 'no'] If 'yes', the specGND variable will be 0162 % saved to disk after the permutation test is completed 0163 % and added to it. User will first be prompted to verify 0164 % file name and path. {default: 'yes'} 0165 % 0166 % Outputs: 0167 % specGND - specGND structure variable. This is the same as 0168 % the input specGND variable with one addition: the 0169 % field specGND.t_tests will contain the results of the 0170 % permutation test and the test parameters. 0171 % pval - A two-dimensional matrix (channel x time) of the 0172 % p-values of each comparison (no correction for multiple 0173 % comparisons). 0174 % data_t - A two-dimensional matrix (channel x frequency) of the 0175 % t-scores of each comparison. 0176 % crit_t - The critical t-score(s) for the test. Any t-scores 0177 % that are more extreme than the critical t-score(s) 0178 % significantly deviate from 0. 0179 % 0180 % Note also that a great deal of information about the test is displayed 0181 % in the MATLAB command window. You can easiy record of all this 0182 % information to a text file using the MATLAB command "diary." 0183 % 0184 % Global Variables: 0185 % VERBLEVEL - MATLABmk level of verbosity (i.e., tells 0186 % functions how much to report about what they're doing during 0187 % runtime) set by the optional function argument 'verblevel' 0188 % 0189 % Notes: 0190 % -To add a difference wave to a specGND variable, use the function "bin_dif_spec". 0191 % 0192 % 0193 % References: 0194 % Benjamini, Y., & Hochberg, Y. (1995). Controlling the false discovery 0195 % rate: A practical and powerful approach to multiple testing. Journal of 0196 % the Royal Statistical Society. Series B (Methodological), 57(1), 289-300. 0197 % 0198 % Benjamini, Y., Krieger, A. M., & Yekutieli, D. (2006). Adaptive linear 0199 % step-up procedures that control the false discovery rate. Biometrika, 0200 % 93(3), 491-507. 0201 % 0202 % Benjamini, Y., & Yekutieli, D. (2001). The control of the false 0203 % discovery rate in multiple testing under dependency. The Annals of 0204 % Statistics, 29(4), 1165-1188. 0205 % 0206 % 0207 % Author: 0208 % David Groppe 0209 % Kutaslab, 12/2010 0210 0211 %%%%%%%%%%%%%%%% REVISION LOG %%%%%%%%%%%%%%%%% 0212 % ?/?/2010- 0213 0214 0215 function [specGND, pval, data_t, crit_t]=tfdr_specGND(specGND_or_fname,bin,varargin) 0216 0217 global VERBLEVEL; 0218 0219 p=inputParser; 0220 p.addRequired('specGND_or_fname',@(x) ischar(x) || isstruct(x)); 0221 p.addRequired('bin',@(x) isnumeric(x) && (length(x)==1) && (x>0)); 0222 p.addParamValue('q',0.05,@(x) isnumeric(x) && (x>0) && (x<1)); 0223 p.addParamValue('method','bh',@(x) strcmpi(x,'bh') || strcmpi(x,'by') || strcmpi(x,'bky')); 0224 p.addParamValue('tail',0,@(x) isnumeric(x) && (length(x)==1)); 0225 p.addParamValue('freq_band',[],@(x) isnumeric(x) && (size(x,2)==2)); 0226 p.addParamValue('freq_block_dur',[],@isnumeric); 0227 p.addParamValue('mean_band','no',@(x) ischar(x) && (strcmpi(x,'yes') || strcmpi(x,'no'))); 0228 p.addParamValue('null_mean',0,@(x) isnumeric(x) && (length(x)==1)); 0229 p.addParamValue('verblevel',[],@(x) isnumeric(x) && (length(x)==1)); 0230 p.addParamValue('exclude_chans',[],@(x) ischar(x) || iscell(x)); 0231 p.addParamValue('include_chans',[],@(x) ischar(x) || iscell(x)); 0232 p.addParamValue('plot_gui','yes',@(x) ischar(x) && (strcmpi(x,'yes') || strcmpi(x,'no'))); 0233 p.addParamValue('plot_raster','yes',@(x) ischar(x) && (strcmpi(x,'yes') || strcmpi(x,'no'))); 0234 p.addParamValue('plot_mn_topo',[],@(x) ischar(x) && (strcmpi(x,'yes') || strcmpi(x,'no'))); 0235 p.addParamValue('output_file',[],@ischar); 0236 p.addParamValue('save_specGND','yes',@(x) ischar(x) && (strcmpi(x,'yes') || strcmpi(x,'no'))); 0237 0238 p.parse(specGND_or_fname,bin,varargin{:}); 0239 0240 0241 if isempty(p.Results.verblevel), 0242 if isempty(VERBLEVEL), 0243 VERBLEVEL=2; 0244 end 0245 else 0246 VERBLEVEL=p.Results.verblevel; 0247 end 0248 0249 mean_band=str2bool(p.Results.mean_band); 0250 0251 %Load specGND struct 0252 if ischar(specGND_or_fname), 0253 fprintf('Loading specGND struct from file %s.\n',specGND_or_fname); 0254 load(specGND_or_fname,'-MAT'); 0255 if ~exist('specGND') 0256 error('File %s does not contain a specGND variable.',specGND_or_fname); 0257 end 0258 else 0259 specGND=specGND_or_fname; 0260 clear specGND_or_fname; 0261 fldnames=fieldnames(specGND); 0262 if ismember('group_desc',fldnames), 0263 error('You passed a specGRP or GRP variable to this function instead of a specGND variable.'); 0264 end 0265 end 0266 [n_chan, n_freq, n_bin, total_subs]=size(specGND.indiv_pow_dB); 0267 VerbReport(sprintf('Experiment: %s',specGND.exp_desc),2,VERBLEVEL); 0268 0269 if (bin>n_bin), 0270 error('There is no Bin %d in this specGND variable.',bin); 0271 end 0272 0273 if ~specGND.dif(bin), 0274 error('You can only apply this function to bins that are the difference between two other bins.'); 0275 end 0276 0277 %Use only subs with data in relevant bin(s) 0278 use_subs=find(specGND.indiv_bin_ct(:,p.Results.bin)); 0279 n_sub=length(use_subs); 0280 VerbReport(sprintf('%d out of %d participants have data in relevant bin.',n_sub,total_subs), ... 0281 1,VERBLEVEL); 0282 0283 0284 %% Figure out which channels to ignore if any 0285 %But first make sure exclude & include options were not both used. 0286 if ~isempty(p.Results.include_chans) && ~isempty(p.Results.exclude_chans) 0287 error('You cannot use BOTH ''include_chans'' and ''exclude_chans'' options.'); 0288 end 0289 if ischar(p.Results.exclude_chans), 0290 exclude_chans{1}=p.Results.exclude_chans; 0291 elseif isempty(p.Results.exclude_chans) 0292 exclude_chans=[]; 0293 else 0294 exclude_chans=p.Results.exclude_chans; 0295 end 0296 if ischar(p.Results.include_chans), 0297 include_chans{1}=p.Results.include_chans; 0298 elseif isempty(p.Results.include_chans) 0299 include_chans=[]; 0300 else 0301 include_chans=p.Results.include_chans; 0302 end 0303 0304 % exclude and include chan options 0305 if ~isempty(exclude_chans), 0306 ignore_chans=zeros(1,length(exclude_chans)); %preallocate mem 0307 ct=0; 0308 for x=1:length(exclude_chans), 0309 found=0; 0310 for c=1:n_chan, 0311 if strcmpi(exclude_chans{x},specGND.chanlocs(c).labels), 0312 found=1; 0313 ct=ct+1; 0314 ignore_chans(ct)=c; 0315 end 0316 end 0317 if ~found, 0318 watchit(sprintf('I attempted to exclude %s. However no such electrode was found in specGND variable.', ... 0319 exclude_chans{x})); 0320 end 0321 end 0322 ignore_chans=ignore_chans(1:ct); 0323 use_chans=setdiff(1:n_chan,ignore_chans); 0324 elseif ~isempty(include_chans), 0325 use_chans=zeros(1,length(include_chans)); %preallocate mem 0326 ct=0; 0327 for x=1:length(include_chans), 0328 found=0; 0329 for c=1:n_chan, 0330 if strcmpi(include_chans{x},specGND.chanlocs(c).labels), 0331 found=1; 0332 ct=ct+1; 0333 use_chans(ct)=c; 0334 end 0335 end 0336 if ~found, 0337 watchit(sprintf('I attempted to include %s. However no such electrode was found in specGND variable.', ... 0338 include_chans{x})); 0339 end 0340 end 0341 use_chans=use_chans(1:ct); 0342 else 0343 use_chans=1:n_chan; 0344 end 0345 0346 0347 %% Find frequencies to analyze 0348 if isempty(p.Results.freq_band), 0349 freq_band=[0 specGND.freqs(end)]; %default frequency band 0350 else 0351 freq_band=p.Results.freq_band; 0352 end 0353 freq_band=sort(freq_band,2); %first make sure earlier of each pair of frequencies is first 0354 freq_band=sort(freq_band,1); %next sort frequency bands from earliest to latest lower frequency boundary 0355 n_band=size(freq_band,1); 0356 if ~isempty(p.Results.freq_block_dur), 0357 tbd=p.Results.freq_block_dur; 0358 n_block=length(tbd); 0359 if (n_block>1) && (n_block~=n_band), 0360 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); 0361 elseif (n_block==1) && (n_band>1), 0362 tbd=repmat(tbd,1,n_band); 0363 end 0364 0365 f_bands=[]; 0366 one_step=specGND.freqs(2)-specGND.freqs(1); %DFT frequency interval 0367 for a=1:n_band, 0368 band_strt=freq_band(a,1); 0369 band_stop=freq_band(a,2); 0370 new_bands=[band_strt:tbd(a):band_stop]; 0371 for b=1:(length(new_bands)-1), 0372 f_bands=[f_bands; new_bands(b) new_bands(b+1)-one_step]; 0373 end 0374 end 0375 freq_band=f_bands; 0376 clear f_bands; 0377 n_band=size(freq_band,1); 0378 end 0379 0380 if mean_band, 0381 use_freqs=cell(1,n_band); 0382 else 0383 use_freqs=[]; 0384 end 0385 for a=1:n_band, 0386 VerbReport(sprintf('Frequency Band #%d:',a),1,VERBLEVEL); 0387 VerbReport(sprintf('Attempting to use frequency boundaries of %f to %f Hz for hypothesis test.',freq_band(a,1),freq_band(a,2)), ... 0388 1,VERBLEVEL); 0389 start_freq=find_tpt(freq_band(a,1),specGND.freqs); 0390 end_freq=find_tpt(freq_band(a,2),specGND.freqs); 0391 if mean_band, 0392 use_freqs{a}=[start_freq:end_freq]; 0393 else 0394 use_freqs=[use_freqs [start_freq:end_freq]]; 0395 end 0396 %replace desired frequencies with closest matches 0397 freq_band(a,1)=specGND.freqs(start_freq); 0398 freq_band(a,2)=specGND.freqs(end_freq); 0399 VerbReport(sprintf('Exact frequency boundaries are %f to %f Hz (that''s from frequency# %d to %d).', ... 0400 freq_band(a,1),freq_band(a,2),start_freq,end_freq),1,VERBLEVEL); 0401 end 0402 if ~mean_band, 0403 use_freqs=unique(use_freqs); %sorts frequencies and gets rid of any redundant frequencies 0404 end 0405 0406 %% Compile data 0407 if mean_band, 0408 %Take mean amplitude in frequency blocks and then test 0409 mean_pwr=zeros(length(use_chans),n_band,n_sub); 0410 for a=1:n_band, 0411 for sub=1:n_sub, 0412 mean_pwr(:,a,sub)=mean(specGND.indiv_pow_dB(use_chans,use_freqs{a},bin,use_subs(sub)),2); 0413 end 0414 end 0415 else 0416 %Use every single frequency in frequency band(s) 0417 n_use_freqs=length(use_freqs); 0418 mean_pwr=zeros(length(use_chans),n_use_freqs,n_sub); 0419 for sub=1:n_sub, 0420 mean_pwr(:,:,sub)=specGND.indiv_pow_dB(use_chans,use_freqs,bin,use_subs(sub)); 0421 end 0422 end 0423 0424 0425 %% Report tail of test 0426 VerbReport(sprintf('Testing null hypothesis that the grand average power in Bin %d (%s) has a mean of %f (uV^2)/Hz.',bin, ... 0427 specGND.bindesc{bin},p.Results.null_mean),1,VERBLEVEL); 0428 if p.Results.tail==0 0429 VerbReport(sprintf('Alternative hypothesis is that the grand average power differs from %f (i.e., two-tailed test).',p.Results.null_mean), ... 0430 1,VERBLEVEL); 0431 elseif p.Results.tail<0, 0432 VerbReport(sprintf('Alternative hypothesis is that the grand average power is less than %f (i.e., lower-tailed test).',p.Results.null_mean), ... 0433 1,VERBLEVEL); 0434 else 0435 VerbReport(sprintf('Alternative hypothesis is that the grand average power is greater than %f (i.e., upper-tailed test).',p.Results.null_mean), ... 0436 1,VERBLEVEL); 0437 end 0438 0439 %Compute t-tests 0440 [p_values, data_t]=fast_t1(mean_pwr-p.Results.null_mean,p.Results.tail,VERBLEVEL); 0441 0442 %FDR correction for multiple comparisons 0443 switch lower(p.Results.method) 0444 case 'bh' 0445 VerbReport('FDR control procedure: Benjamini & Hochberg (independent or positive dependency)',1,VERBLEVEL); 0446 [h_rej, crit_p]=fdr_bh(p_values,p.Results.q,'pdep','no'); 0447 case 'by' 0448 VerbReport('FDR control procedure: Benjamini & Yekutieli (arbitrary dependency)',1,VERBLEVEL); 0449 [h_rej, crit_p]=fdr_bh(p_values,p.Results.q,'dep','no'); 0450 case 'bky' 0451 VerbReport('FDR control procedure: Benjamini, Krieger, & Yekutieli (two-stage)',1,VERBLEVEL); 0452 [h_rej, crit_p]=fdr_bky(p_values,p.Results.q,'no'); 0453 end 0454 0455 %Compute the permutation test ?? DELETE 0456 %freq_domain=1; 0457 %[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); 0458 0459 0460 sig_freqs=find(sum(h_rej)); 0461 if isempty(sig_freqs), 0462 if p.Results.tail==0, 0463 crit_t=[NaN NaN]; 0464 else 0465 crit_t=NaN; 0466 end 0467 fprintf('Power spectra are NOT significantly different from %f (q=%f) at any time point/window analyzed.\n', ... 0468 p.Results.null_mean,p.Results.q); 0469 %fprintf('All p-values>=%f\n',min(min(p_values))); ?? add FDR-corrected 0470 %p-values someday? 0471 else 0472 crit_t=min(abs(data_t(crit_p==p_values))); 0473 if p.Results.tail<0, 0474 crit_t=-crit_t; 0475 elseif p.Results.tail==0, 0476 crit_t=[-crit_t crit_t]; 0477 end 0478 0479 %Command line summary of results 0480 VerbReport(['Critical t-score(s):' num2str(crit_t)],1,VERBLEVEL); 0481 if p.Results.tail 0482 %one-tailed test 0483 tw_alpha=1-cdf('t',max(abs(crit_t)),n_sub-1); 0484 else 0485 %two-tailed test 0486 tw_alpha=(1-cdf('t',max(abs(crit_t)),n_sub-1))*2; 0487 end 0488 VerbReport(sprintf('That corresponds to a test-wise alpha level of %f.',tw_alpha),1,VERBLEVEL); 0489 VerbReport(sprintf('Bonferroni test-wise alpha would be %f.',p.Results.q/numel(p_values)) ... 0490 ,1,VERBLEVEL); 0491 0492 fprintf('Total number of significant differences: %d\n',sum(sum(h_rej))); 0493 fprintf('Estimated upper bound on expected number of false discoveries: %.1f\n',sum(sum(h_rej))*p.Results.q); 0494 fprintf('Significant differences from zero (in order of earliest to latest):\n'); 0495 %max_sig_p=0; ?? report FDR corrected p-values someday? 0496 %time windows instead of time points 0497 for t=sig_freqs, 0498 if mean_band 0499 fprintf('%f to %f Hz, electrode(s): ',specGND.freqs(use_freqs{t}(1)), ... 0500 specGND.freqs(use_freqs{t}(end))); 0501 else 0502 fprintf('%f Hz, electrode(s): ',specGND.freqs(use_freqs(t))); 0503 end 0504 sig_elec=find(h_rej(:,t)); 0505 ct=0; 0506 for c=sig_elec', 0507 ct=ct+1; 0508 % if p_values(c,t)>max_sig_p, ?? report FDR 0509 % corrected p-values? 0510 % max_sig_p=p_values(c,t); 0511 % end 0512 if ct==length(sig_elec), 0513 fprintf('%s.\n',specGND.chanlocs(use_chans(c)).labels); 0514 else 0515 fprintf('%s, ',specGND.chanlocs(use_chans(c)).labels); 0516 end 0517 end 0518 end 0519 %fprintf('All significant p-values<=%f\n',max_sig_p); %?? make FDR 0520 %corrected? 0521 end 0522 0523 0524 %Add t-test results to specGND struct 0525 n_t_tests=length(specGND.t_tests); 0526 neo_test=n_t_tests+1; 0527 specGND.t_tests(neo_test).bin=bin; 0528 specGND.t_tests(neo_test).freq_band=freq_band; 0529 specGND.t_tests(neo_test).used_freq_ids=use_freqs; 0530 n_use_chans=length(use_chans); 0531 include_chans=cell(1,n_use_chans); 0532 for a=1:n_use_chans, 0533 include_chans{a}=specGND.chanlocs(use_chans(a)).labels; 0534 end 0535 specGND.t_tests(neo_test).include_chans=include_chans; 0536 specGND.t_tests(neo_test).used_chan_ids=use_chans; 0537 specGND.t_tests(neo_test).mult_comp_method=p.Results.method; 0538 specGND.t_tests(neo_test).n_perm=NaN; 0539 specGND.t_tests(neo_test).desired_alphaORq=p.Results.q; 0540 specGND.t_tests(neo_test).estimated_alpha=NaN; 0541 specGND.t_tests(neo_test).null_mean=p.Results.null_mean; 0542 if mean_band, 0543 specGND.t_tests(neo_test).data_t=data_t; 0544 specGND.t_tests(neo_test).mean_band='yes'; 0545 else 0546 specGND.t_tests(neo_test).data_t='See specGND.grands_t'; 0547 specGND.t_tests(neo_test).mean_band='no'; 0548 end 0549 specGND.t_tests(neo_test).crit_t=crit_t; 0550 specGND.t_tests(neo_test).df=length(use_subs)-1; 0551 specGND.t_tests(neo_test).adj_pval=NaN; % ?? make FDR adjusted p-vals someday? 0552 specGND.t_tests(neo_test).fdr_rej=h_rej; 0553 specGND.t_tests(neo_test).seed_state=NaN; 0554 0555 if strcmpi(p.Results.plot_raster,'yes'), 0556 sig_raster(specGND,neo_test,'verblevel',0); 0557 end 0558 0559 if mean_band, 0560 if strcmpi(p.Results.plot_mn_topo,'yes') || isempty(p.Results.plot_mn_topo), 0561 sig_topo(specGND,neo_test,'units','t','verblevel',0); %t-score topographies 0562 sig_topo(specGND,neo_test,'units','dB','verblevel',0); %microvolt topographies 0563 end 0564 else 0565 if strcmpi(p.Results.plot_gui,'yes'), 0566 gui_pow('initialize','specGND_or_specGRP',specGND,'t_test',neo_test,'stat','t', ... 0567 'verblevel',1); 0568 end 0569 end 0570 0571 if ~isempty(p.Results.output_file) 0572 [fid msg]=fopen(p.Results.output_file,'w'); 0573 if fid==-1, 0574 error('Cannot create file %s for writing. According to fopen.m: %s.', ... 0575 p.Results.file,msg); 0576 else 0577 %Write header column of frequencies 0578 % Leave first column blank for channel labels 0579 if mean_band, 0580 for t=1:n_band 0581 fprintf(fid,' %f-%f',specGND.freqs(use_freqs{t}(1)), ... 0582 specGND.freqs(use_freqs{t}(end))); 0583 end 0584 0585 %write a couple spaces and then write header for t-scores 0586 fprintf(fid,' '); 0587 for t=1:n_band 0588 fprintf(fid,' %f-%f',specGND.freqs(use_freqs{t}(1)), ... 0589 specGND.freqs(use_freqs{t}(end))); 0590 end 0591 else 0592 for t=use_freqs, 0593 fprintf(fid,' %f',specGND.freqs(t)); 0594 end 0595 end 0596 fprintf(fid,' Hz\n'); 0597 0598 % Write channel labels and p-values 0599 chan_ct=0; 0600 for c=use_chans, 0601 chan_ct=chan_ct+1; 0602 fprintf(fid,'%s',specGND.chanlocs(c).labels); 0603 for t=1:length(use_freqs), 0604 fprintf(fid,' %f',pval(chan_ct,t)); 0605 end 0606 fprintf(fid,' p-value(uncorrected)'); %make corrected some day? 0607 0608 if mean_band, 0609 %write a couple spaces and then write t-scores if mean amp 0610 %in frequency bands used 0611 fprintf(fid,' '); 0612 for t=1:n_band 0613 fprintf(fid,' %f',data_t(chan_ct,t)); 0614 end 0615 fprintf(fid,' t-score \n'); 0616 else 0617 fprintf(fid,'\n'); 0618 end 0619 end 0620 0621 % Write permutation test details 0622 fprintf(fid,'Experiment: %s\n',specGND.exp_desc); 0623 fprintf(fid,'Test_of_null_hypothesis_that_Bin_%d_equals_%f\n',bin,p.Results.null_mean); 0624 fprintf(fid,'#_of_frequency_bands: %d\n',n_band); 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,'Critical_p-value: %f\n',crit_p); 0637 fprintf(fid,'Total_#_of_significant_differences: %d\n',sum(sum(h_rej))); 0638 fprintf(fid,'Estimated_upper_bound_on_expected_#_of_false_discoveries: %.1f\n',sum(sum(h_rej))*p.Results.q); 0639 fprintf(fid,'Degrees_of_freedom: %d\n',length(use_subs)-1); 0640 fprintf(fid,'q_level: %f\n',p.Results.q); 0641 0642 % # of participants and filenames 0643 fprintf(fid,'#_of_participants: %d\n',length(use_subs)); 0644 fprintf(fid,'Participant_names: \n'); 0645 for s=1:length(use_subs), 0646 fprintf(fid,'%s\n',specGND.indiv_subnames{use_subs(s)}); 0647 end 0648 end 0649 fclose(fid); 0650 end 0651 0652 if ~strcmpi(p.Results.save_specGND,'no'), 0653 specGND=save_matmk(specGND); 0654 end 0655 0656 % 0657 %% %%%%%%%%%%%%%%%%%%%%% function str2bool() %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0658 % 0659 function bool=str2bool(str) 0660 %function bool=str2bool(str) 0661 0662 if ischar(str), 0663 if strcmpi(str,'yes') || strcmpi(str,'y') 0664 bool=1; 0665 else 0666 bool=0; 0667 end 0668 else 0669 bool=str; 0670 end 0671