Home > matlabmk > clustGND.m

clustGND

PURPOSE ^

clustGND() - Tests the null hypothesis that the grand average voltage

SYNOPSIS ^

function [GND, prm_pval, data_t]=clustGND(GND_or_fname,bin,varargin)

DESCRIPTION ^

 clustGND() - 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 a cluster-based 
             permutation test and the "cluster mass" statistic 
             (Bullmore et al., 1999; Maris & Oostenveld, 2007).  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 clustGRP.m.
             
             
 Usage:
  >> [GND, prm_pval, data_t]=clustGND(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
                  specify a filename be sure to include the file's path,
                  unless the file is in the current working directory.
   bin          - [integer] The bin to contrast against a voltage of 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}
   alpha         - A number between 0 and 1 specifying the family-wise
                   error rate (FWER) of the test. Note, cluster-based tests 
                   only weak control of FWER. {default: 0.05}
   thresh_p      - The test-wise p-value threshold for cluster inclusion. If
                   a channel/time-point has a t-score that corresponds to an
                   uncorrected p-value greater than thresh_p, it is assigned
                   a p-value of 1 and not considered for clustering. Note
                   that thresh_p automatically takes into account the tail of
                   the test (e.g., you will get positive and negative t-score
                   thresholds for a two-tailed test).
   chan_hood     - A scalar or a 2D symmetric binary matrix that indicates
                   which channels are considered neighbors of other 
                   channels. E.g., if chan_hood(2,10)=1, then Channel 2 
                   and Channel 10 are nieghbors. You can produce a 
                   chan_hood matrix using the function spatial_neighbors.m. 
                   If a scalar is provided, then all electrodes within that 
                   distance of a particular electrode are considered 
                   neighbors. Note, EEGLAB's electrode coordinates assume
                   the head has a radius of 1. See the help documentation 
                   of the function spatial_neighbors to see how you could
                   convert this distance threshold to centimeters. 
                   {default: 0.61}
   head_radius   - The radius of the head in whatever units the Cartesian
                   coordinates in GRP.chanlocs are in. This is used to
                   convert scalar values of chan_hood into centimeters.
                   {default: []}
   n_perm        - The number of permutations to use in the test.  As this
                   value increases, the test takes longer to compute and 
                   the results become more reliable.  Manly (1997) suggests 
                   using at least 1000 permutations for an alpha level of 
                   0.05 and at least 5000 permutations for an alpha level 
                   of 0.01. {default: 2500}
   time_wind     - Pair of time values specifying the beginning and
                   end of a time window in ms (e.g., [160 180]). 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
                   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}
   mean_wind     - ['yes' or 'no'] If 'yes', the permutation test will be
                   performed on the mean amplitude within the 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
                   permutation test (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 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.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_raster   - ['yes' or 'no'] If 'yes', a two-dimensional (time 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 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 topography of the mean
                   voltages/effects in the 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., number of permutations, 
                   family-wise alpha 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. {default:
                   none}
   save_GND      - ['yes' or 'no'] If 'yes', the GND variable will be
                   saved to disk after the permutation test is completed 
                   and added to it. User will first be prompted to verify 
                   file name and path. {default: 'yes'}
   reproduce_test- [integer] The number of the permutation test stored in
                   the GND variable to reproduce.  For example, if 
                   'reproduce_test' equals 2, the second t-test 
                   stored in the GND variable (i.e., GND.t_tests(2)) will 
                   be reproduced.  Reproduction is accomplished by setting
                   the random number generator used in the permutation test 
                   to the same initial state it was in when the permutation 
                   test was first applied.

 Outputs:
   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 
                   permutation test and the test parameters. 
   prm_pval      - A two-dimensional matrix (channel x time) of the
                   p-values of each comparison.  These p-values are
                   corrected for multiple comparisons by the permutation
                   test.
   data_t        - A two-dimensional matrix (channel x time) of the
                   t-scores of each comparison.

   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".

 -Unlike a parametric test (e.g., an ANOVA), a discrete set of p-values
 are possible (at most the number of possible permutations).  Since the
 number of possible permutations grows rapdily with the number of
 participants, this is only issue for small sample sizes (e.g., 6
 participants).  When you have such a small sample size, the
 limited number of p-values may make the test less conservative (e.g., 
 you might be forced to use an alpha level of .0286 since it is the biggest
 possible alpha level less than .05).

 -For two-tailed tests clustGND.m effectively performs an upper-tailed 
 test and a lower-tailed test and multiplies the resulting p-values by 2.
 In other words, it does two tests and then applies Bonferroni correction.
 This is faster than performing a "proper" two-tailed permutation test
 since you only have to form clusters for one polarity (e.g., if you form
 negative clusters for each permutation, the negative of that distribution
 will be valid for computing p-values for postive clusters) but can
 result in p-values greater than 1.

 Author:
 David Groppe
 May, 2011
 Kutaslab, San Diego

 References:
 Bullmore, E. T., Suckling, J., Overmeyer, S., Rabe-Hesketh, S., Taylor, 
 E., & Brammer, M. J. (1999). Global, voxel, and cluster tests, by theory 
 and permutation, for a difference between two groups of structural MR 
 images of the brain. IEEE Transactions on Medical Imaging, 18(1), 32-42. 
 doi:10.1109/42.750253

 Maris, E., & Oostenveld, R. (2007). Nonparametric statistical testing of 
 EEG- and MEG-data. Journal of Neuroscience Methods, 164(1), 177-190. 
 doi:10.1016/j.jneumeth.2007.03.024

 Manly, B.F.J. (1997) Randomization, bootstrap, and Monte Carlo methods in
 biology. 2nd ed. Chapmn and Hall, London.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 % clustGND() - 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 a cluster-based
0004 %             permutation test and the "cluster mass" statistic
0005 %             (Bullmore et al., 1999; Maris & Oostenveld, 2007).  Note, mu
0006 %             is assumed to be 0 by default. This function requires
0007 %             individual subject ERPs to be stored in a "GND" structure and
0008 %             outputs the test results in a number of graphical and text
0009 %             formats.  For analogous between-subject comparisons use the
0010 %             function clustGRP.m.
0011 %
0012 %
0013 % Usage:
0014 %  >> [GND, prm_pval, data_t]=clustGND(GND_or_fname,bin,varargin);
0015 %
0016 % Required Inputs:
0017 %   GND_or_fname - A GND structure variable or the filename of a
0018 %                  GND structure that has been saved to disk.  To
0019 %                  create a GND variable from Kutaslab ERP files (e.g.,
0020 %                  *.mas files) use avgs2GND.m.  To do the same from
0021 %                  EEGLAB *.set files use sets2GND.m.  See Mass
0022 %                  Univariate ERP Toolbox documentation for detailed
0023 %                  information about the format of a GND variable. If you
0024 %                  specify a filename be sure to include the file's path,
0025 %                  unless the file is in the current working directory.
0026 %   bin          - [integer] The bin to contrast against a voltage of mu
0027 %                  Use the function headinfo.m to see what bins are stored
0028 %                  in a GND variable.  Use the function bin_dif.m to create
0029 %                  a difference wave between two bins whose significance
0030 %                  you can test with this function.
0031 %
0032 % Optional Inputs:
0033 %   tail          - [1,0, or -1] An integer specifying the tail of the
0034 %                   hypothesis test.  "1" means upper-tailed (i.e., alternative
0035 %                   hypothesis is that the ERP/difference wave is greater
0036 %                   than the mean of the null hypothesis).  "0" means two-
0037 %                   tailed (i.e., alternative hypothesis is that the
0038 %                   ERP/difference wave is not equal to the mean of the null
0039 %                   hypothesis). "-1" means lower-tailed (i.e., alternative
0040 %                   hypothesis is that the ERP/difference wave is less than
0041 %                   the mean of the null hypothesis). {default: 0}
0042 %   alpha         - A number between 0 and 1 specifying the family-wise
0043 %                   error rate (FWER) of the test. Note, cluster-based tests
0044 %                   only weak control of FWER. {default: 0.05}
0045 %   thresh_p      - The test-wise p-value threshold for cluster inclusion. If
0046 %                   a channel/time-point has a t-score that corresponds to an
0047 %                   uncorrected p-value greater than thresh_p, it is assigned
0048 %                   a p-value of 1 and not considered for clustering. Note
0049 %                   that thresh_p automatically takes into account the tail of
0050 %                   the test (e.g., you will get positive and negative t-score
0051 %                   thresholds for a two-tailed test).
0052 %   chan_hood     - A scalar or a 2D symmetric binary matrix that indicates
0053 %                   which channels are considered neighbors of other
0054 %                   channels. E.g., if chan_hood(2,10)=1, then Channel 2
0055 %                   and Channel 10 are nieghbors. You can produce a
0056 %                   chan_hood matrix using the function spatial_neighbors.m.
0057 %                   If a scalar is provided, then all electrodes within that
0058 %                   distance of a particular electrode are considered
0059 %                   neighbors. Note, EEGLAB's electrode coordinates assume
0060 %                   the head has a radius of 1. See the help documentation
0061 %                   of the function spatial_neighbors to see how you could
0062 %                   convert this distance threshold to centimeters.
0063 %                   {default: 0.61}
0064 %   head_radius   - The radius of the head in whatever units the Cartesian
0065 %                   coordinates in GRP.chanlocs are in. This is used to
0066 %                   convert scalar values of chan_hood into centimeters.
0067 %                   {default: []}
0068 %   n_perm        - The number of permutations to use in the test.  As this
0069 %                   value increases, the test takes longer to compute and
0070 %                   the results become more reliable.  Manly (1997) suggests
0071 %                   using at least 1000 permutations for an alpha level of
0072 %                   0.05 and at least 5000 permutations for an alpha level
0073 %                   of 0.01. {default: 2500}
0074 %   time_wind     - Pair of time values specifying the beginning and
0075 %                   end of a time window in ms (e.g., [160 180]). Every
0076 %                   single time point in the time window will be individually
0077 %                   tested (i.e., maximal temporal resolution) if mean_wind
0078 %                   option is NOT used. Note, boundaries of time window
0079 %                   may not exactly correspond to desired time window
0080 %                   boundaries because of temporal digitization (i.e., you
0081 %                   only have samples every so many ms). {default: 0 ms to
0082 %                   the end of the epoch}
0083 %   mean_wind     - ['yes' or 'no'] If 'yes', the permutation test will be
0084 %                   performed on the mean amplitude within the time window
0085 %                   specified by time_wind.  This sacrifices temporal
0086 %                   resolution to increase test power by reducing the number
0087 %                   of comparisons.  If 'no', every single time point within
0088 %                   time_wind's time windows will be tested individually.
0089 %                   {default: 'no'}
0090 %   null_mean     - [number] The mean of the null hypothesis (i.e., mu) in
0091 %                   units of microvolts. {default: 0}
0092 %   exclude_chans - A cell array of channel labels to exclude from the
0093 %                   permutation test (e.g., {'A2','lle','rhe'}).  This option
0094 %                   sacrifices spatial resolution to increase test power by
0095 %                   reducing the number of comparisons. Use headinfo.m to see
0096 %                   the channel labels stored in the GND variable. You cannot
0097 %                   use both this option and 'include_chans' (below).{default:
0098 %                   not used, all channels included in test}
0099 %   include_chans - A cell array of channel labels to use in the permutation
0100 %                   test (e.g., {'A2','lle','rhe'}).  All other channels will
0101 %                   be ignored. This option sacrifices spatial resolution to
0102 %                   increase test power by reducing the number of comparisons.
0103 %                   Use headinfo.m to see the channel labels stored in the GND
0104 %                   variable. You cannot use both this option and
0105 %                   'exclude_chans' (above). {default: not used, all channels
0106 %                   included in test}
0107 %   verblevel     - An integer specifiying the amount of information you want
0108 %                   this function to provide about what it is doing during runtime.
0109 %                    Options are:
0110 %                      0 - quiet, only show errors, warnings, and EEGLAB reports
0111 %                      1 - stuff anyone should probably know
0112 %                      2 - stuff you should know the first time you start working
0113 %                          with a data set {default value}
0114 %                      3 - stuff that might help you debug (show all
0115 %                          reports)
0116 %   plot_raster   - ['yes' or 'no'] If 'yes', a two-dimensional (time x channel)
0117 %                   binary "raster" diagram is created to illustrate the
0118 %                   results of the permutation test.  Significant negative and
0119 %                   positive deviations from the null hypothesis are shown
0120 %                   as black and white rectangles respectively. Non-
0121 %                   significant comparisons are shown as gray rectangles.
0122 %                   Clicking on the rectangles will show you the
0123 %                   corresponding time and channel label for that
0124 %                   rectangle. This figure can be reproduced with the
0125 %                   function sig_raster.m. {default: 'yes'}
0126 %   plot_mn_topo  - ['yes' or 'no'] If 'yes', the topography of the mean
0127 %                   voltages/effects in the time window is produced.  More
0128 %                   specifically, two figures are produced: one showing the
0129 %                   topographies in uV the other in t-scores. Significant/
0130 %                   nonsignificant comparisons are shown as white/black
0131 %                   electrodes. Clicking on electrodes will show the
0132 %                   electrode's name.  This figure can be reproduced with
0133 %                   the function sig_topo.m.  This option has NO effect if
0134 %                   mean_wind option is set to 'no'. {default: 'yes'}
0135 %   output_file   - A string indicating the name of a space delimited text
0136 %                   file to produce containing the p-values of all comparisons
0137 %                   and the details of the test (e.g., number of permutations,
0138 %                   family-wise alpha level, etc...). If mean_wind option is
0139 %                   set to 'yes,' t-scores of each comparison are also
0140 %                   included since you cannot derive them from the t-scores
0141 %                   at each time point/electrode in a simple way. When
0142 %                   importing this file into a spreadsheet be sure NOT to count
0143 %                   consecutive spaces as multiple delimiters. {default:
0144 %                   none}
0145 %   save_GND      - ['yes' or 'no'] If 'yes', the GND variable will be
0146 %                   saved to disk after the permutation test is completed
0147 %                   and added to it. User will first be prompted to verify
0148 %                   file name and path. {default: 'yes'}
0149 %   reproduce_test- [integer] The number of the permutation test stored in
0150 %                   the GND variable to reproduce.  For example, if
0151 %                   'reproduce_test' equals 2, the second t-test
0152 %                   stored in the GND variable (i.e., GND.t_tests(2)) will
0153 %                   be reproduced.  Reproduction is accomplished by setting
0154 %                   the random number generator used in the permutation test
0155 %                   to the same initial state it was in when the permutation
0156 %                   test was first applied.
0157 %
0158 % Outputs:
0159 %   GND           - GND structure variable.  This is the same as
0160 %                   the input GND variable with one addition: the
0161 %                   field GND.t_tests will contain the results of the
0162 %                   permutation test and the test parameters.
0163 %   prm_pval      - A two-dimensional matrix (channel x time) of the
0164 %                   p-values of each comparison.  These p-values are
0165 %                   corrected for multiple comparisons by the permutation
0166 %                   test.
0167 %   data_t        - A two-dimensional matrix (channel x time) of the
0168 %                   t-scores of each comparison.
0169 %
0170 %   Note also that a great deal of information about the test is displayed
0171 %   in the MATLAB command window.  You can easiy record of all this
0172 %   information to a text file using the MATLAB command "diary."
0173 %
0174 % Global Variables:
0175 %   VERBLEVEL - Mass Univariate ERP Toolbox level of verbosity (i.e., tells
0176 %               functions how much to report about what they're doing during
0177 %               runtime) set by the optional function argument 'verblevel'
0178 %
0179 % Notes:
0180 % -To add a difference wave to a GND variable, use the function "bin_dif".
0181 %
0182 % -Unlike a parametric test (e.g., an ANOVA), a discrete set of p-values
0183 % are possible (at most the number of possible permutations).  Since the
0184 % number of possible permutations grows rapdily with the number of
0185 % participants, this is only issue for small sample sizes (e.g., 6
0186 % participants).  When you have such a small sample size, the
0187 % limited number of p-values may make the test less conservative (e.g.,
0188 % you might be forced to use an alpha level of .0286 since it is the biggest
0189 % possible alpha level less than .05).
0190 %
0191 % -For two-tailed tests clustGND.m effectively performs an upper-tailed
0192 % test and a lower-tailed test and multiplies the resulting p-values by 2.
0193 % In other words, it does two tests and then applies Bonferroni correction.
0194 % This is faster than performing a "proper" two-tailed permutation test
0195 % since you only have to form clusters for one polarity (e.g., if you form
0196 % negative clusters for each permutation, the negative of that distribution
0197 % will be valid for computing p-values for postive clusters) but can
0198 % result in p-values greater than 1.
0199 %
0200 % Author:
0201 % David Groppe
0202 % May, 2011
0203 % Kutaslab, San Diego
0204 %
0205 % References:
0206 % Bullmore, E. T., Suckling, J., Overmeyer, S., Rabe-Hesketh, S., Taylor,
0207 % E., & Brammer, M. J. (1999). Global, voxel, and cluster tests, by theory
0208 % and permutation, for a difference between two groups of structural MR
0209 % images of the brain. IEEE Transactions on Medical Imaging, 18(1), 32-42.
0210 % doi:10.1109/42.750253
0211 %
0212 % Maris, E., & Oostenveld, R. (2007). Nonparametric statistical testing of
0213 % EEG- and MEG-data. Journal of Neuroscience Methods, 164(1), 177-190.
0214 % doi:10.1016/j.jneumeth.2007.03.024
0215 %
0216 % Manly, B.F.J. (1997) Randomization, bootstrap, and Monte Carlo methods in
0217 % biology. 2nd ed. Chapmn and Hall, London.
0218 
0219 %%%%%%%%%%%%%%%% REVISION LOG %%%%%%%%%%%%%%%%%
0220 % 12/11/2011-Now uses Amy Guthormsen's recursiveless find_clusters.m.
0221 %
0222 
0223 function [GND, prm_pval, data_t]=clustGND(GND_or_fname,bin,varargin)
0224 
0225 global VERBLEVEL;
0226 
0227 p=inputParser;
0228 p.addRequired('GND_or_fname',@(x) ischar(x) || isstruct(x));
0229 p.addRequired('bin',@(x) isnumeric(x) && (length(x)==1) && (x>0));
0230 p.addParamValue('tail',0,@(x) isnumeric(x) && (length(x)==1));
0231 p.addParamValue('alpha',0.05,@(x) (length(x)==1) && isnumeric(x) && (x>0) && (x<1));
0232 p.addParamValue('thresh_p',0.05,@(x) (length(x)==1) && isnumeric(x) && (x>0) && (x<1));
0233 p.addParamValue('chan_hood',0.61,@isnumeric);
0234 p.addParamValue('head_radius',[],@isnumeric);
0235 p.addParamValue('time_wind',[],@(x) isnumeric(x) && (size(x,2)==2));
0236 p.addParamValue('mean_wind','no',@(x) ischar(x) && (strcmpi(x,'yes') || strcmpi(x,'no')));
0237 p.addParamValue('n_perm',2500,@(x) isnumeric(x) && (length(x)==1));
0238 p.addParamValue('null_mean',0,@(x) isnumeric(x) && (length(x)==1));
0239 p.addParamValue('verblevel',[],@(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('plot_raster','yes',@(x) ischar(x) && (strcmpi(x,'yes') || strcmpi(x,'no')));
0243 p.addParamValue('plot_gui','yes',@(x) ischar(x) && (strcmpi(x,'yes') || strcmpi(x,'no')));
0244 p.addParamValue('plot_mn_topo',[],@(x) ischar(x) && (strcmpi(x,'yes') || strcmpi(x,'no')));
0245 p.addParamValue('output_file',[],@ischar);
0246 p.addParamValue('reproduce_test',[],@(x) isnumeric(x) && (length(x)==1));
0247 p.addParamValue('save_GND','yes',@(x) ischar(x) && (strcmpi(x,'yes') || strcmpi(x,'no')));
0248 
0249 p.parse(GND_or_fname,bin,varargin{:});
0250 
0251 if isempty(p.Results.verblevel),
0252     if isempty(VERBLEVEL),
0253         VERBLEVEL=2;
0254     end
0255 else
0256    VERBLEVEL=p.Results.verblevel; 
0257 end
0258     
0259 mean_wind=str2bool(p.Results.mean_wind);
0260 
0261 %Load GND struct
0262 if ischar(GND_or_fname),
0263     fprintf('Loading GND struct from file %s.\n',GND_or_fname);
0264     load(GND_or_fname,'-MAT');
0265     if ~exist('GND')
0266         error('File %s does not contain a GND variable.',GND_or_fname);
0267     end
0268 else
0269     GND=GND_or_fname;
0270     clear GND_or_fname;
0271     fldnames=fieldnames(GND);
0272     if ismember('group_desc',fldnames),
0273        error('You passed a GRP variable to this function instead of a GND variable.'); 
0274     end
0275 end
0276 [n_chan, n_pt, n_bin, total_subs]=size(GND.indiv_erps);
0277 VerbReport(sprintf('Experiment: %s',GND.exp_desc),2,VERBLEVEL);
0278 
0279 if (bin>n_bin),
0280     error('There is no Bin %d in this GND variable.',bin); 
0281 end
0282 
0283 %Use only subs with data in relevant bin(s)
0284 use_subs=find(GND.indiv_bin_ct(:,p.Results.bin));
0285 n_sub=length(use_subs);
0286 VerbReport(sprintf('%d out of %d participants have data in relevant bin.',n_sub,total_subs), ...
0287     1,VERBLEVEL);
0288 
0289 
0290 %% Figure out which channels to ignore if any
0291 %But first make sure exclude & include options were not both used.
0292 if ~isempty(p.Results.include_chans) && ~isempty(p.Results.exclude_chans)
0293     error('You cannot use BOTH ''include_chans'' and ''exclude_chans'' options.');
0294 end
0295 if ischar(p.Results.exclude_chans),
0296     exclude_chans{1}=p.Results.exclude_chans;
0297 elseif isempty(p.Results.exclude_chans)
0298     exclude_chans=[];
0299 else
0300     exclude_chans=p.Results.exclude_chans;
0301 end
0302 if ischar(p.Results.include_chans),
0303     include_chans{1}=p.Results.include_chans;
0304 elseif isempty(p.Results.include_chans)
0305     include_chans=[];
0306 else
0307     include_chans=p.Results.include_chans;
0308 end
0309 
0310 % exclude and include chan options
0311 if ~isempty(exclude_chans),
0312     ignore_chans=zeros(1,length(exclude_chans)); %preallocate mem
0313     ct=0;
0314     for x=1:length(exclude_chans),
0315         found=0;
0316         for c=1:n_chan,
0317             if strcmpi(exclude_chans{x},GND.chanlocs(c).labels),
0318                 found=1;
0319                 ct=ct+1;
0320                 ignore_chans(ct)=c;
0321             end
0322         end
0323         if ~found,
0324             watchit(sprintf('I attempted to exclude %s.  However no such electrode was found in GND variable.', ...
0325                 exclude_chans{x}));
0326         end
0327     end
0328     ignore_chans=ignore_chans(1:ct);
0329     use_chans=setdiff(1:n_chan,ignore_chans);
0330 elseif ~isempty(include_chans),
0331     use_chans=zeros(1,length(include_chans)); %preallocate mem
0332     ct=0;
0333     for x=1:length(include_chans),
0334         found=0;
0335         for c=1:n_chan,
0336             if strcmpi(include_chans{x},GND.chanlocs(c).labels),
0337                 found=1;
0338                 ct=ct+1;
0339                 use_chans(ct)=c;
0340             end
0341         end
0342         if ~found,
0343             watchit(sprintf('I attempted to include %s.  However no such electrode was found in GND variable.', ...
0344                 include_chans{x}));
0345         end
0346     end
0347     use_chans=use_chans(1:ct);
0348 else
0349     use_chans=1:n_chan;
0350 end
0351 
0352 % Establish spatial neighborhood matrix for making clusters
0353 if isscalar(p.Results.chan_hood),
0354     chan_hood=spatial_neighbors(GND.chanlocs(use_chans),p.Results.chan_hood,p.Results.head_radius);
0355 else
0356     chan_hood=p.Results.chan_hood;
0357 end
0358 
0359 
0360 %% Find time points
0361 if isempty(p.Results.time_wind),
0362     time_wind=[0 GND.time_pts(end)]; %default time window
0363 else
0364     time_wind=p.Results.time_wind;
0365 end
0366 time_wind=sort(time_wind,2); %first make sure earlier of each pair of time points is first
0367 time_wind=sort(time_wind,1); %next sort time windows from earliest to latest onset
0368 n_wind=size(time_wind,1);
0369 if n_wind>1,
0370    error('clustGND.m can only handle a single mean time window at the moment.  Try tmaxGND.m or tfdrGND.m if you want to simultaneously test hypotheses in multiple mean time windows.'); 
0371 end
0372 
0373 if mean_wind,
0374     use_tpts=cell(1,n_wind);
0375 else
0376     use_tpts=[];
0377 end
0378 for a=1:n_wind,
0379     VerbReport(sprintf('Time Window #%d:',a),1,VERBLEVEL);
0380     VerbReport(sprintf('Attempting to use time boundaries of %d to %d ms for hypothesis test.',time_wind(a,1),time_wind(a,2)), ...
0381         1,VERBLEVEL);
0382     start_tpt=find_tpt(time_wind(a,1),GND.time_pts);
0383     end_tpt=find_tpt(time_wind(a,2),GND.time_pts);
0384     if mean_wind,
0385         use_tpts{a}=[start_tpt:end_tpt];
0386     else
0387         use_tpts=[use_tpts [start_tpt:end_tpt]];
0388     end
0389     %replace desired time points with closest matches
0390     time_wind(a,1)=GND.time_pts(start_tpt);
0391     time_wind(a,2)=GND.time_pts(end_tpt);
0392     VerbReport(sprintf('Exact window boundaries are %d to %d ms (that''s from time point %d to %d).', ...
0393         time_wind(a,1),time_wind(a,2),start_tpt,end_tpt),1,VERBLEVEL);
0394 end
0395 if ~mean_wind,
0396     use_tpts=unique(use_tpts); %sorts time points and gets rid of any redundant time points
0397 end
0398 
0399 
0400 %% Compile data
0401 if mean_wind,
0402     %Take mean amplitude in time blocks and then test
0403     erps=zeros(length(use_chans),n_wind,n_sub);
0404     for a=1:n_wind,
0405         for sub=1:n_sub,
0406             erps(:,a,sub)=mean(GND.indiv_erps(use_chans,use_tpts{a},bin,use_subs(sub)),2);
0407         end
0408     end
0409 else
0410     %Use every single time point in time window(s)
0411     n_use_tpts=length(use_tpts);
0412     erps=zeros(length(use_chans),n_use_tpts,n_sub);    
0413     for sub=1:n_sub,
0414         erps(:,:,sub)=GND.indiv_erps(use_chans,use_tpts,bin,use_subs(sub));
0415     end
0416 end
0417 
0418 
0419 %% Report tail of test & alpha levels
0420 VerbReport(sprintf('Testing null hypothesis that the grand average ERPs in Bin %d (%s) have a mean of %f microvolts.',bin, ...
0421     GND.bin_info(bin).bindesc,p.Results.null_mean),1,VERBLEVEL);
0422 if p.Results.tail==0
0423     VerbReport(sprintf('Alternative hypothesis is that the ERPs differ from %f (i.e., two-tailed test).',p.Results.null_mean), ...
0424         1,VERBLEVEL);
0425 elseif p.Results.tail<0,
0426     VerbReport(sprintf('Alternative hypothesis is that the ERPs are less than %f (i.e., lower-tailed test).',p.Results.null_mean), ...
0427         1,VERBLEVEL);
0428 else
0429     VerbReport(sprintf('Alternative hypothesis is that the ERPs are greater than %f (i.e., upper-tailed test).',p.Results.null_mean), ...
0430         1,VERBLEVEL);
0431 end
0432 
0433 
0434 %% Optionally reset random number stream to reproduce a previous test
0435 if isempty(p.Results.reproduce_test),
0436     seed_state=[];
0437 else
0438     if p.Results.reproduce_test>length(GND.t_tests),
0439         error('Value of argument ''reproduce_test'' is too high.  You only have %d t-tests stored with this GND variable.',length(GND.t_tests));
0440     else
0441         if isnan(GND.t_tests(p.Results.reproduce_test).n_perm)
0442             error('t-test set %d is NOT a permutation test. You don''t need to seed the random number generator to reproduce it.', ...
0443                 p.Results.reproduce_test);
0444         else
0445             seed_state=GND.t_tests(p.Results.reproduce_test).seed_state;
0446         end
0447     end
0448 end
0449 
0450 %% Compute the permutation test
0451 [prm_pval, data_t, clust_info, seed_state, est_alpha]=clust_perm1(erps-p.Results.null_mean, ...
0452     chan_hood,p.Results.n_perm,p.Results.alpha,p.Results.tail,p.Results.thresh_p, ...
0453     VERBLEVEL,seed_state,0);
0454 
0455 %% Command line summary of results
0456 if p.Results.tail>=0,
0457     %upper or two-tailed test
0458     n_pos=length(clust_info.pos_clust_pval);
0459     fprintf('# of positive clusters found: %d\n',n_pos);
0460     fprintf('# of significant positive clusters found: %d\n',sum(clust_info.pos_clust_pval<est_alpha));
0461     fprintf('Positive cluster p-values range from %g to %g.\n',min(clust_info.pos_clust_pval), ...
0462         max(clust_info.pos_clust_pval));
0463 end
0464 if p.Results.tail<=0,
0465     %lower or two-tailed test
0466     n_neg=length(clust_info.neg_clust_pval);
0467     fprintf('# of negative clusters found: %d\n',n_neg);
0468     fprintf('# of significant negative clusters found: %d\n',sum(clust_info.neg_clust_pval<est_alpha));
0469     fprintf('Negative cluster p-values range from %g to %g.\n',min(clust_info.neg_clust_pval), ...
0470         max(clust_info.neg_clust_pval));
0471 end
0472 
0473 sig_ids=find(prm_pval<p.Results.alpha);
0474 if isempty(sig_ids),
0475     fprintf('ERPs are NOT significantly different from zero (alpha=%f) at any time point/window analyzed.\n', ...
0476         p.Results.alpha);
0477     fprintf('All p-values>=%f\n',min(min(prm_pval)));
0478 else
0479     [dummy min_t_id]=min(abs(data_t(sig_ids)));
0480     min_t=data_t(sig_ids(min_t_id));
0481     VerbReport(['Smallest significant t-score(s):' num2str(min_t)],1,VERBLEVEL);
0482     if p.Results.tail
0483         %one-tailed test
0484         tw_alpha=1-cdf('t',max(abs(min_t)),n_sub-1);
0485     else
0486         %two-tailed test
0487         tw_alpha=(1-cdf('t',max(abs(min_t)),n_sub-1))*2;
0488     end
0489     
0490     VerbReport(sprintf('That corresponds to a test-wise alpha level of %f.',tw_alpha),1,VERBLEVEL);
0491     VerbReport(sprintf('Bonferroni test-wise alpha would be %f.',p.Results.alpha/(size(prm_pval,1)* ...
0492         size(prm_pval,2))),1,VERBLEVEL);
0493     sig_tpts=find(sum(prm_pval<p.Results.alpha));
0494     
0495     fprintf('Significant differences from zero (in order of earliest to latest):\n');
0496     max_sig_p=0;
0497     min_sig_p=2;
0498     for t=sig_tpts,
0499         if mean_wind
0500             %time windows instead of time points
0501             fprintf('%d to %d ms, electrode(s): ',GND.time_pts(use_tpts{t}(1)), ...
0502                 GND.time_pts(use_tpts{t}(end)));
0503         else
0504             fprintf('%d ms, electrode(s): ',GND.time_pts(use_tpts(t)));
0505         end
0506         sig_elec=find(prm_pval(:,t)<p.Results.alpha);
0507         ct=0;
0508         for c=sig_elec',
0509             ct=ct+1;
0510             if prm_pval(c,t)>max_sig_p,
0511                 max_sig_p=prm_pval(c,t);
0512             end
0513             if prm_pval(c,t)<min_sig_p,
0514                 min_sig_p=prm_pval(c,t);
0515             end
0516             
0517             if ct==length(sig_elec),
0518                 fprintf('%s.\n',GND.chanlocs(use_chans(c)).labels);
0519             else
0520                 fprintf('%s, ',GND.chanlocs(use_chans(c)).labels);
0521             end
0522         end
0523     end
0524     fprintf('All significant corrected p-values are between %f and %f\n',max_sig_p,min_sig_p);
0525 end
0526 
0527 
0528 %% Add permutation results to GND struct
0529 n_t_tests=length(GND.t_tests);
0530 neo_test=n_t_tests+1;
0531 GND.t_tests(neo_test).bin=bin;
0532 GND.t_tests(neo_test).time_wind=time_wind;
0533 GND.t_tests(neo_test).used_tpt_ids=use_tpts;
0534 n_use_chans=length(use_chans);
0535 include_chans=cell(1,n_use_chans);
0536 for a=1:n_use_chans,
0537    include_chans{a}=GND.chanlocs(use_chans(a)).labels; 
0538 end
0539 GND.t_tests(neo_test).include_chans=include_chans;
0540 GND.t_tests(neo_test).used_chan_ids=use_chans;
0541 GND.t_tests(neo_test).mult_comp_method='cluster mass perm test';
0542 GND.t_tests(neo_test).n_perm=p.Results.n_perm;
0543 GND.t_tests(neo_test).desired_alphaORq=p.Results.alpha;
0544 GND.t_tests(neo_test).estimated_alpha=est_alpha;
0545 GND.t_tests(neo_test).null_mean=p.Results.null_mean;
0546 if mean_wind,
0547     GND.t_tests(neo_test).data_t=data_t;
0548     GND.t_tests(neo_test).mean_wind='yes';
0549 else
0550     GND.t_tests(neo_test).data_t='See GND.grands_t';
0551     GND.t_tests(neo_test).mean_wind='no';
0552 end
0553 GND.t_tests(neo_test).crit_t=NaN;
0554 GND.t_tests(neo_test).df=length(use_subs)-1;
0555 GND.t_tests(neo_test).adj_pval=prm_pval;
0556 GND.t_tests(neo_test).fdr_rej=NaN;
0557 GND.t_tests(neo_test).seed_state=seed_state;
0558 GND.t_tests(neo_test).clust_info=clust_info;
0559 GND.t_tests(neo_test).chan_hood=chan_hood;
0560 
0561 if strcmpi(p.Results.plot_raster,'yes'),
0562     sig_raster(GND,neo_test,'verblevel',0,'use_color','rgb');
0563 end
0564 
0565 if mean_wind,
0566     if strcmpi(p.Results.plot_mn_topo,'yes') || isempty(p.Results.plot_mn_topo),
0567         sig_topo(GND,neo_test,'units','t','verblevel',0); %t-score topographies
0568         sig_topo(GND,neo_test,'units','uV','verblevel',0); %microvolt topographies
0569     end
0570 else
0571     if strcmpi(p.Results.plot_gui,'yes'),
0572         gui_erp('initialize','GNDorGRP',GND,'t_test',neo_test,'stat','t', ...
0573             'verblevel',1);
0574     end
0575 end
0576 
0577 %% Create text file output
0578 if ~isempty(p.Results.output_file)
0579     [fid msg]=fopen(p.Results.output_file,'w');
0580     if fid==-1,
0581         error('Cannot create file %s for writing.  According to fopen.m: %s.', ...
0582             p.Results.file,msg);
0583     else
0584         %Write header column of times
0585         % Leave first column blank for channel labels
0586         if mean_wind,
0587             for t=1:n_wind
0588                 fprintf(fid,' %d-%d',GND.time_pts(use_tpts{t}(1)), ...
0589                     GND.time_pts(use_tpts{t}(end)));
0590             end
0591             
0592             %write a couple spaces and then write header for t-scores
0593             fprintf(fid,'  ');
0594             for t=1:n_wind
0595                 fprintf(fid,' %d-%d',GND.time_pts(use_tpts{t}(1)), ...
0596                     GND.time_pts(use_tpts{t}(end)));
0597             end
0598         else
0599             for t=use_tpts,
0600                 fprintf(fid,' %d',GND.time_pts(t));
0601             end
0602         end
0603         fprintf(fid,' Milliseconds\n');
0604         
0605         % Write channel labels and p-values
0606         chan_ct=0;
0607         for c=use_chans,
0608             chan_ct=chan_ct+1;
0609             fprintf(fid,'%s',GND.chanlocs(c).labels);
0610             for t=1:length(use_tpts),
0611                 fprintf(fid,' %f',prm_pval(chan_ct,t));
0612             end
0613             fprintf(fid,' p-value');
0614             
0615             if mean_wind,
0616                 %write a couple spaces and then write t-scores if mean amp
0617                 %in time windows used
0618                 fprintf(fid,' ');
0619                 for t=1:n_wind
0620                     fprintf(fid,' %f',data_t(chan_ct,t));
0621                 end
0622                 fprintf(fid,' t-score \n');
0623             else
0624                 fprintf(fid,'\n');
0625             end
0626         end
0627         
0628         % Write permutation test details
0629         fprintf(fid,'Experiment: %s\n',GND.exp_desc);
0630         fprintf(fid,'Test_of_null_hypothesis_that_Bin_%d_equals: %f\n',bin,p.Results.null_mean);
0631         fprintf(fid,'#_of_time_windows: %d\n',n_wind);
0632         fprintf(fid,'#_of_permutations: %d\n',p.Results.n_perm);
0633         fprintf(fid,'Tail_of_test: ');
0634         if ~p.Results.tail,
0635             fprintf(fid,'Two_tailed\n');
0636         elseif p.Results.tail>0
0637             fprintf(fid,'Upper_tailed\n');
0638         else
0639             fprintf(fid,'Lower_tailed\n');
0640         end
0641         fprintf(fid,'Degrees_of_freedom: %d\n',length(use_subs)-1);
0642         fprintf(fid,'Family_wise_error_rate: %f\n',p.Results.alpha);
0643         fprintf(fid,'Cluster_inclusion_p_value_threshold: %f\n',p.Results.thresh_p);
0644         if isscalar(p.Results.chan_hood)
0645             fprintf(fid,'Max_spatial_neighbor_distance: %f\n',p.Results.chan_hood);
0646         end
0647         
0648         % # of participants and filenames
0649         fprintf(fid,'#_of_participants: %d\n',length(use_subs));
0650         fprintf(fid,'Participant_names: \n');
0651         for s=1:length(use_subs),
0652             fprintf(fid,'%s\n',GND.indiv_subnames{use_subs(s)});
0653         end
0654     end
0655     fclose(fid);
0656 end
0657 
0658 if ~strcmpi(p.Results.save_GND,'no'),
0659     GND=save_matmk(GND);
0660 end
0661 
0662 %
0663 %% %%%%%%%%%%%%%%%%%%%%% function str2bool() %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0664 %
0665 function bool=str2bool(str)
0666 %function bool=str2bool(str)
0667 
0668 if ischar(str),
0669     if strcmpi(str,'yes') || strcmpi(str,'y')
0670         bool=1;
0671     else
0672         bool=0;
0673     end
0674 else
0675    bool=str; 
0676 end
0677

Generated on Tue 10-May-2016 16:37:45 by m2html © 2005