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