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