tfdr_specGRP() - Tests the null hypothesis that the grand average voltage
of a between-subject difference in power spectra is mu
(usually 0) using t-tests based and one of various Benjamini
algorithms for control of the false discovery rate. 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 tfdr_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, p_values, data_t, crit_t]=tfdr_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}
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.
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,
false discovery rate, 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 t-tests are completed
and added to it. User will first be prompted to verify
file name and path. {default: 'yes'}
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
t-tests and the test parameters.
p_values - A two-dimensional matrix (channel x frequency) 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 bin that is the difference in spectral power between two groups
to a specGRP variable, use the function "bin_op_specGRP".
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_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 t-tests based and one of various Benjamini 0004 % algorithms for control of the false discovery rate. 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 tfdr_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, p_values, data_t, crit_t]=tfdr_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 % q - A number between 0 and 1 specifying the family-wise 0041 % q level of the test. q is the upper bound on the 0042 % expected proportion of rejected null hypotheses that are 0043 % false rejections (i.e., the FDR). {default: 0.05} 0044 % method - ['bh', 'by', or 'bky'] The procedure used to control 0045 % the FDR. 'bh' is the classic Benjamini & Hochberg (1995) 0046 % procedure, which is guaranteed to control FDR when the 0047 % tests are independent or positively dependent (e.g., 0048 % positively correlated Gaussians). 'by' is a much more 0049 % conservative version of 'bh' that always controls FDR 0050 % (regardless of the dependency structure of the tests-- 0051 % Benjamini & Yekutieli, 2001). 'bky' is a "two-stage" 0052 % version of 'bh' that is more powerful than 'bh' when a 0053 % lot of the null hypotheses are false (Benjamini, Krieger, & 0054 % Yekutieli, 2006). 'bky' is guaranteed to control FDR when the 0055 % tests are independent and tends to be slightly less 0056 % powerful than 'bh' when few or no null hypothese are 0057 % false. 0058 % freq_band - 2D matrix of pairs of frequency values specifying the 0059 % beginning and end of one or more frequency bands in Hz 0060 % (e.g., [3 8; 8 12]). Every single frequency in 0061 % the frequency band will be individually tested (i.e., 0062 % maximal temporal resolution) if mean_band option is 0063 % NOT used. Note, boundaries of frequency band(s) may not 0064 % exactly correspond to desired frequency band boundaries 0065 % because of temporal digitization (i.e., you only have 0066 % samples every so many Hz). {default: 0 Hz to the 0067 % Nyquist frequency} 0068 % freq_block_dur- [integers] A number or numbers (in Hz) specifying a 0069 % duration of frequency blocks in which the frequency 0070 % windows specified by freq_band will be divided. For 0071 % example, if freq_band=[0 40] and freq_block_dur=10, 0072 % the 0 to 40 Hz frequency band will be sub-divided into 0073 % 10 Hz bands (i.e., freq_band will equal something like 0074 % [0 9; 10 19; 20 29; 30 39]--depending on your frequency 0075 % resolution). This is an easy way to break up larger 0076 % frequency bands of interest into smaller bands for 0077 % mean band analysis (if the 'mean_band' option is set to 0078 % 'yes'). If you specify multiple frequency bands with 0079 % 'freq_band', you can break them up using durations of 0080 % different lengths. For example, if freq_band=[0 10; 0081 % 40 90] and freq_block_dur=[5 10], the first frequency 0082 % band (0 to 10 Hz) will be broken up into 5 Hz 0083 % bands and the second band (40 to 90 Hz) will be 0084 % broken up into 10 Hz bands. {default: not used} 0085 % mean_band - ['yes' or 'no'] If 'yes', the permutation test will be 0086 % performed on the mean power within each frequency band 0087 % specified by freq_band. This sacrifices frequency 0088 % resolution to increase test power by reducing the number 0089 % of comparisons. If 'no', every single frequency within 0090 % freq_band's frequency bands will be tested individually. 0091 % {default: 'no'} 0092 % null_mean - [number] The mean of the null hypothesis (in units of 0093 % microvolts). {default: 0} 0094 % exclude_chans - A cell array of channel labels to exclude from the 0095 % permutation test (e.g., {'A2','lle','rhe'}). This option 0096 % sacrifices spatial resolution to increase test power by 0097 % reducing the number of comparisons. Use headinfo_spec.m to see 0098 % the channel labels stored in the specGRP variable. You cannot 0099 % use both this option and 'include_chans' (below).{default: 0100 % not used, all channels included in test} 0101 % include_chans - A cell array of channel labels to use in the permutation 0102 % test (e.g., {'A2','lle','rhe'}). All other channels will 0103 % be ignored. This option sacrifices spatial resolution to 0104 % increase test power by reducing the number of comparisons. 0105 % Use headinfo_spec.m to see the channel labels stored in the specGRP 0106 % variable. You cannot use both this option and 0107 % 'exclude_chans' (above). {default: not used, all channels 0108 % included in test} 0109 % verblevel - An integer specifiying the amount of information you want 0110 % this function to provide about what it is doing during runtime. 0111 % Options are: 0112 % 0 - quiet, only show errors, warnings, and EEGLAB reports 0113 % 1 - stuff anyone should probably know 0114 % 2 - stuff you should know the first time you start working 0115 % with a data set {default value} 0116 % 3 - stuff that might help you debug (show all 0117 % reports) 0118 % plot_gui - ['yes' or 'no'] If 'yes', a GUI is created for 0119 % visualizing the results of the permutation test using the 0120 % function gui_pow.m. The GUI vizualizes the grand average 0121 % ERPs in each bin via various stats (e.g., t-scores), shows 0122 % topographies at individual frequencies, and illustrates 0123 % which electrodes significantly differ from the null 0124 % hypothesis. This option does not work if mean_band 0125 % option is set to 'yes.' This GUI can be reproduced using 0126 % the function gui_pow.m. {default: 'yes'} 0127 % plot_raster - ['yes' or 'no'] If 'yes', a two-dimensional (frequency x channel) 0128 % binary "raster" diagram is created to illustrate the 0129 % results of the permutation test. Significant negative and 0130 % positive deviations from the null hypothesis are shown 0131 % as black and white rectangles respectively. Non- 0132 % significant comparisons are shown as gray rectangles. 0133 % Clicking on the rectangles will show you the 0134 % corresponding frequency and channel label for that 0135 % rectangle. This figure can be reproduced with the 0136 % function sig_raster.m. {default: 'yes'} 0137 % plot_mn_topo - ['yes' or 'no'] If 'yes', the topographies of the mean 0138 % voltages/effects in each frequency band is produced. More 0139 % specifically, two figures are produced: one showing the 0140 % topographies in decibel power the other in t-scores. Significant/ 0141 % nonsignificant comparisons are shown as white/black 0142 % electrodes. Clicking on electrodes will show the 0143 % electrode's name. This figure can be reproduced with 0144 % the function sig_topo.m. This option has NO effect if 0145 % mean_band option is set to 'no'. {default: 'yes'} 0146 % output_file - A string indicating the name of a space delimited text 0147 % file to produce containing the p-values of all comparisons 0148 % and the details of the test (e.g., number of permutations, 0149 % false discovery rate, etc...). If mean_band option is 0150 % set to 'yes,' t-scores of each comparison are also 0151 % included since you cannot derive them from the t-scores 0152 % at each frequency/electrode in a simple way. When 0153 % importing this file into a spreadsheet be sure NOT to count 0154 % consecutive spaces as multiple delimiters. {default: none} 0155 % save_specGRP - ['yes' or 'no'] If 'yes', the specGRP variable will be 0156 % saved to disk after the t-tests are completed 0157 % and added to it. User will first be prompted to verify 0158 % file name and path. {default: 'yes'} 0159 % 0160 % Outputs: 0161 % specGRP - specGRP structure variable. This is the same as 0162 % the input specGRP variable with one addition: the 0163 % field specGRP.t_tests will contain the results of the 0164 % t-tests and the test parameters. 0165 % p_values - A two-dimensional matrix (channel x frequency) of the 0166 % p-values of each comparison (no correction for multiple 0167 % comparisons). 0168 % data_t - A two-dimensional matrix (channel x frequency) of the 0169 % t-scores of each comparison. 0170 % crit_t - The critical t-score(s) for the test. Any t-scores 0171 % that are more extreme than the critical t-score(s) 0172 % significantly deviate from 0. 0173 % 0174 % Note also that a great deal of information about the test is displayed 0175 % in the MATLAB command window. You can easiy record of all this 0176 % information to a text file using the MATLAB command "diary." 0177 % 0178 % Global Variables: 0179 % VERBLEVEL = MATLABmk level of verbosity (i.e., tells functions how much 0180 % to report about what they're doing during runtime) set by 0181 % the optional function argument 'verblevel' 0182 % 0183 % Notes: 0184 % -To add a bin that is the difference in spectral power between two groups 0185 % to a specGRP variable, use the function "bin_op_specGRP". 0186 % 0187 % 0188 % References: 0189 % Benjamini, Y., & Hochberg, Y. (1995). Controlling the false discovery 0190 % rate: A practical and powerful approach to multiple testing. Journal of 0191 % the Royal Statistical Society. Series B (Methodological), 57(1), 289-300. 0192 % 0193 % Benjamini, Y., Krieger, A. M., & Yekutieli, D. (2006). Adaptive linear 0194 % step-up procedures that control the false discovery rate. Biometrika, 0195 % 93(3), 491-507. 0196 % 0197 % Benjamini, Y., & Yekutieli, D. (2001). The control of the false 0198 % discovery rate in multiple testing under dependency. The Annals of 0199 % Statistics, 29(4), 1165-1188. 0200 % 0201 % 0202 % 0203 % Author: 0204 % David Groppe 0205 % Kutaslab, 12/2010 0206 0207 %%%%%%%%%%%%%%%% REVISION LOG %%%%%%%%%%%%%%%%% 0208 % ?/?/2010-?? 0209 % 0210 0211 function [specGRP, p_values, data_t, crit_t]=tfdr_specGRP(specGRP_or_fname,bin,varargin) 0212 0213 global VERBLEVEL; 0214 0215 p=inputParser; 0216 p.addRequired('specGRP_or_fname',@(x) ischar(x) || isstruct(x)); 0217 p.addRequired('bin',@(x) isnumeric(x) && (length(x)==1) && (x>0)); 0218 p.addParamValue('tail',0,@(x) isnumeric(x) && (length(x)==1)); 0219 p.addParamValue('q',0.05,@(x) isnumeric(x) && (x>0) && (x<1)); 0220 p.addParamValue('method','bh',@(x) strcmpi(x,'bh') || strcmpi(x,'by') || strcmpi(x,'bky')); 0221 p.addParamValue('freq_band',[],@(x) isnumeric(x) && (size(x,2)==2)); 0222 p.addParamValue('freq_block_dur',[],@isnumeric); 0223 p.addParamValue('mean_band','no',@(x) ischar(x) && (strcmpi(x,'yes') || strcmpi(x,'no'))); 0224 p.addParamValue('null_mean',0,@(x) isnumeric(x) && (length(x)==1)); 0225 p.addParamValue('verblevel',[],@(x) isnumeric(x) && (length(x)==1)); 0226 p.addParamValue('exclude_chans',[],@(x) ischar(x) || iscell(x)); 0227 p.addParamValue('include_chans',[],@(x) ischar(x) || iscell(x)); 0228 p.addParamValue('plot_gui','yes',@(x) ischar(x) && (strcmpi(x,'yes') || strcmpi(x,'no'))); 0229 p.addParamValue('plot_raster','yes',@(x) ischar(x) && (strcmpi(x,'yes') || strcmpi(x,'no'))); 0230 p.addParamValue('plot_mn_topo',[],@(x) ischar(x) && (strcmpi(x,'yes') || strcmpi(x,'no'))); 0231 p.addParamValue('output_file',[],@ischar); 0232 p.addParamValue('save_specGRP','yes',@(x) ischar(x) && (strcmpi(x,'yes') || strcmpi(x,'no'))); 0233 0234 p.parse(specGRP_or_fname,bin,varargin{:}); 0235 0236 0237 if isempty(p.Results.verblevel), 0238 if isempty(VERBLEVEL), 0239 VERBLEVEL=2; 0240 end 0241 else 0242 VERBLEVEL=p.Results.verblevel; 0243 end 0244 0245 mean_band=str2bool(p.Results.mean_band); 0246 0247 %Load specGRP struct 0248 if ischar(specGRP_or_fname), 0249 fprintf('Loading specGRP struct from file %s.\n',specGRP_or_fname); 0250 load(specGRP_or_fname,'-MAT'); 0251 else 0252 specGRP=specGRP_or_fname; 0253 clear specGRP_or_fname; 0254 end 0255 [n_chan, n_pt, n_bin]=size(specGRP.grands_pow_dB); 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 & q level 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 0497 %Compute the t-tests 0498 if strcmpi(specGRP.bin_info(bin).op,'(A+B)/n'), 0499 %one sample t-test 0500 VerbReport('Performing one sample/repeated measures t-tests.',1,VERBLEVEL); 0501 powAB=zeros(length(use_chans),size(powA,2),n_subA+n_subB); 0502 powAB(:,:,1:n_subA)=powA-p.Results.null_mean; 0503 powAB(:,:,(n_subA+1):(n_subA+n_subB))=powB-p.Results.null_mean; 0504 [p_values, data_t]=fast_t1(powAB,p.Results.tail,VERBLEVEL); 0505 else 0506 %independent samples t-test 0507 VerbReport('Performing independent samples t-tests.',1,VERBLEVEL); 0508 [p_values, data_t]=fast_t2(powA-p.Results.null_mean,powB,p.Results.tail,VERBLEVEL); 0509 end 0510 %Perform FDR control 0511 switch lower(p.Results.method) 0512 case 'bh' 0513 VerbReport('FDR control procedure: Benjamini & Hochberg (independent or positive dependency)',1,VERBLEVEL); 0514 [h_rej, crit_p]=fdr_bh(p_values,p.Results.q,'pdep','no'); 0515 case 'by' 0516 VerbReport('FDR control procedure: Benjamini & Yekutieli (arbitrary dependency)',1,VERBLEVEL); 0517 [h_rej, crit_p]=fdr_bh(p_values,p.Results.q,'dep','no'); 0518 case 'bky' 0519 VerbReport('FDR control procedure: Benjamini, Krieger, & Yekutieli (two-stage)',1,VERBLEVEL); 0520 [h_rej, crit_p]=fdr_bky(p_values,p.Results.q,'no'); 0521 end 0522 0523 sig_tpts=find(sum(h_rej)); 0524 if isempty(sig_tpts), 0525 if p.Results.tail==0, 0526 crit_t=[NaN NaN]; 0527 else 0528 crit_t=NaN; 0529 end 0530 fprintf('Power spectra are NOT significantly different from zero (q=%f) at any frequency/window analyzed.\n', ... 0531 p.Results.q); 0532 %fprintf('All p-values>=%f\n',min(min(p_values))); ?? add FDR-corrected 0533 %p-values someday? 0534 else 0535 crit_t=min(abs(data_t(crit_p==p_values))); 0536 if p.Results.tail<0, 0537 crit_t=-crit_t; 0538 elseif p.Results.tail==0, 0539 crit_t=[-crit_t crit_t]; 0540 end 0541 0542 %Command line summary of results 0543 VerbReport(['Critical t-score(s):' num2str(crit_t)],1,VERBLEVEL); 0544 if p.Results.tail 0545 %one-tailed test 0546 tw_alpha=1-cdf('t',max(abs(crit_t)),df); 0547 else 0548 %two-tailed test 0549 tw_alpha=2-cdf('t',max(crit_t),df)-cdf('t',-min(crit_t),df); 0550 end 0551 VerbReport(sprintf('That corresponds to a test-wise alpha level of %f.',tw_alpha),1,VERBLEVEL); 0552 VerbReport(sprintf('Bonferroni test-wise alpha would be %f.',p.Results.q/(size(p_values,1)* ... 0553 size(p_values,2))),1,VERBLEVEL); 0554 0555 fprintf('Significant differences from zero (in order of earliest to latest):\n'); 0556 %max_sig_p=0; ?? report FDR corrected p-values someday? 0557 %frequency bands instead of individual frequencies 0558 for t=sig_tpts, 0559 if mean_band 0560 fprintf('%d to %d Hz, electrode(s): ',specGRP.freqs(use_frq_ids{t}(1)), ... 0561 specGRP.freqs(use_frq_ids{t}(end))); 0562 else 0563 fprintf('%d Hz, electrode(s): ',specGRP.freqs(use_frq_ids(t))); 0564 end 0565 sig_elec=find(p_values(:,t)<p.Results.q); 0566 ct=0; 0567 for c=sig_elec', 0568 ct=ct+1; 0569 % if p_values(c,t)>max_sig_p, ?? report FDR corrected p-values someday? 0570 % max_sig_p=p_values(c,t); 0571 % end 0572 if ct==length(sig_elec), 0573 fprintf('%s.\n',specGRP.chanlocs(use_chans(c)).labels); 0574 else 0575 fprintf('%s, ',specGRP.chanlocs(use_chans(c)).labels); 0576 end 0577 end 0578 end 0579 %fprintf('All significant p-values<=%f\n',max_sig_p); ?? add FDR-corrected 0580 %p-values someday? 0581 end 0582 0583 %Add t-test results to specGRP struct 0584 n_t_tests=length(specGRP.t_tests); 0585 neo_test=n_t_tests+1; 0586 specGRP.t_tests(neo_test).bin=bin; 0587 specGRP.t_tests(neo_test).freq_band=freq_band; 0588 specGRP.t_tests(neo_test).used_freq_ids=use_frq_ids; 0589 n_use_chans=length(use_chans); 0590 include_chans=cell(1,n_use_chans); 0591 for a=1:n_use_chans, 0592 include_chans{a}=specGRP.chanlocs(use_chans(a)).labels; 0593 end 0594 specGRP.t_tests(neo_test).include_chans=include_chans; 0595 specGRP.t_tests(neo_test).used_chan_ids=use_chans; 0596 specGRP.t_tests(neo_test).mult_comp_method=p.Results.method; 0597 specGRP.t_tests(neo_test).n_perm=NaN; 0598 specGRP.t_tests(neo_test).desired_alphaORq=p.Results.q; 0599 specGRP.t_tests(neo_test).estimated_alpha=NaN; 0600 specGRP.t_tests(neo_test).null_mean=p.Results.null_mean; 0601 if mean_band, 0602 specGRP.t_tests(neo_test).data_t=data_t; 0603 specGRP.t_tests(neo_test).mean_band='yes'; 0604 else 0605 specGRP.t_tests(neo_test).data_t='See specGRP.grands_pow_dB_t'; 0606 specGRP.t_tests(neo_test).mean_band='no'; 0607 end 0608 specGRP.t_tests(neo_test).crit_t=crit_t; 0609 specGRP.t_tests(neo_test).df=df; 0610 specGRP.t_tests(neo_test).adj_pval=NaN; % ?? make FDR adjusted p-vals someday? 0611 specGRP.t_tests(neo_test).fdr_rej=h_rej; 0612 specGRP.t_tests(neo_test).seed_state=NaN; 0613 0614 if strcmpi(p.Results.plot_raster,'yes'), 0615 sig_raster(specGRP,neo_test,'verblevel',0); 0616 end 0617 0618 if mean_band, 0619 if strcmpi(p.Results.plot_mn_topo,'yes') || isempty(p.Results.plot_mn_topo), 0620 sig_topo(specGRP,neo_test,'units','t','verblevel',0); %t-score topographies 0621 sig_topo(specGRP,neo_test,'units','dB','verblevel',0); %10*log10((uV^2)/Hz) topographies 0622 end 0623 else 0624 %plot_pvals(specGRP,p.Results.alpha,use_chans,use_frq_ids,prm_pval,bin); 0625 %%improve plot_pvals and make an option ?? 0626 if strcmpi(p.Results.plot_gui,'yes'), 0627 gui_pow('initialize','specGND_or_specGRP',specGRP,'t_test',neo_test,'stat','t', ... 0628 'verblevel',1); 0629 end 0630 end 0631 0632 if ~isempty(p.Results.output_file) 0633 [fid msg]=fopen(p.Results.output_file,'w'); 0634 if fid==-1, 0635 error('Cannot create file %s for writing. According to fopen.m: %s.', ... 0636 p.Results.file,msg); 0637 else 0638 %Write header column of frequencies 0639 % Leave first column blank for channel labels 0640 if mean_band, 0641 for t=1:n_band 0642 fprintf(fid,' %d-%d',specGRP.freqs(use_frq_ids{t}(1)), ... 0643 specGRP.freqs(use_frq_ids{t}(end))); 0644 end 0645 0646 %write a couple spaces and then write header for t-scores 0647 fprintf(fid,' '); 0648 for t=1:n_band 0649 fprintf(fid,' %d-%d',specGRP.freqs(use_frq_ids{t}(1)), ... 0650 specGRP.freqs(use_frq_ids{t}(end))); 0651 end 0652 else 0653 for t=use_frq_ids, 0654 fprintf(fid,' %d',specGRP.freqs(t)); 0655 end 0656 end 0657 fprintf(fid,' Hz\n'); 0658 0659 % Write channel labels and p-values 0660 chan_ct=0; 0661 for c=use_chans, 0662 chan_ct=chan_ct+1; 0663 fprintf(fid,'%s',specGRP.chanlocs(c).labels); 0664 for t=1:length(use_tpts), 0665 fprintf(fid,' %f',p_values(chan_ct,t)); % ?? make FDR adjusted p-values someday? 0666 end 0667 fprintf(fid,' p-value'); 0668 0669 if mean_band, 0670 %write a couple spaces and then write t-scores if mean amp 0671 %in frequency bands used 0672 fprintf(fid,' '); 0673 for t=1:n_band 0674 fprintf(fid,' %f',data_t(chan_ct,t)); 0675 end 0676 fprintf(fid,' t-score \n'); 0677 else 0678 fprintf(fid,'\n'); 0679 end 0680 end 0681 0682 % Write t-test details 0683 fprintf(fid,'Experiment: %s\n',specGRP.exp_desc); 0684 fprintf(fid,'Test_of_null_hypothesis_that_Bin_%d_equals_%d\n',bin,p.Results.null_mean); 0685 fprintf(fid,'#_of_frequency_bands: %d\n',n_band); 0686 fprintf(fid,'FDR_method: %s\n',p.Results.method); 0687 fprintf(fid,'Tail_of_test: '); 0688 if ~p.Results.tail, 0689 fprintf(fid,'Two_tailed\n'); 0690 fprintf(fid,'Critical_t_scores: %f %f\n',crit_t(1),crit_t(2)); 0691 elseif p.Results.tail>0 0692 fprintf(fid,'Upper_tailed\n'); 0693 fprintf(fid,'Critical_t_score: %f\n',crit_t(1)); 0694 else 0695 fprintf(fid,'Lower_tailed\n'); 0696 fprintf(fid,'Critical_t_score: %f\n',crit_t(1)); 0697 end 0698 fprintf(fid,'Critical_p-value: %f\n',crit_p); 0699 fprintf(fid,'Total_#_of_significant_differences: %d\n',sum(sum(h_rej))); 0700 fprintf(fid,'Estimated_upper_bound_on_expected_#_of_false_discoveries: %.1f\n',sum(sum(h_rej))*p.Results.q); 0701 fprintf(fid,'Degrees_of_freedom: %d\n',df); 0702 fprintf(fid,'q_level: %f\n',p.Results.q); 0703 for grp=[grpA grpB], 0704 % # of participants and filenames 0705 if grp==grpA, 0706 use_subs=use_subsA; 0707 fprintf(fid,'specGRP_fname_GroupA: %s\n',specGRP.specGND_fnames{grp}); 0708 fprintf(fid,'#_of_participants_GroupA: %d\n',n_subA); 0709 fprintf(fid,'Participant_names_GroupA: \n'); 0710 else 0711 use_subs=use_subsB; 0712 fprintf(fid,'specGRP_fname_GroupB: %s\n',specGRP.specGND_fnames{grp}); 0713 fprintf(fid,'#_of_participants_GroupB: %d\n',n_subB); 0714 fprintf(fid,'Participant_names_GroupB: \n'); 0715 end 0716 for s=1:length(use_subs), 0717 fprintf(fid,'%s\n',specGRP.indiv_subnames{grp}{use_subs(s)}); 0718 end 0719 end 0720 end 0721 fclose(fid); 0722 end 0723 0724 if ~strcmpi(p.Results.save_specGRP,'no'), 0725 specGRP=save_erps(specGRP); 0726 end 0727 0728 0729 % 0730 %% %%%%%%%%%%%%%%%%%%%%% function str2bool() %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0731 % 0732 function bool=str2bool(str) 0733 %function bool=str2bool(str) 0734 0735 if ischar(str), 0736 if strcmpi(str,'yes') || strcmpi(str,'y') 0737 bool=1; 0738 else 0739 bool=0; 0740 end 0741 else 0742 bool=str; 0743 end 0744