rm_regress() - Performs repeated measures multiple regression analysis
using EEG as a response variable via the regression method
described in Lorch & Meyers (1990). The method works by
performing a separate regression analysis for each
participant. Then the regression coefficients are
averaged across participants and their significance
determined by performing a t-test with a null hypothesis
of a coefficient of 0. The regression analysis can be
performed across multiple electrodes. When multiple
electrodes are used, the function corrects for multiple
comparisons by using the Bonferroni-Holms procedure. Each
electrode is considered a distinct test (i.e., the total
number of electrodes equals the number of tests in the
family of tests).
Usage:
>> regress_results=rm_regress(gui_infiles_or_tmplt,time_wind,varargin);
Required Inputs:
gui_infiles_or_tmplt - ['gui', a cell array of strings, or a string template]
If 'gui', a GUI is created that allows you to select
which set files to import (this is probably the
easiest way to import files). If a cell array of
of strings, each element of the cell array should
contain the filename of an EEGLAB set file (e.g.,
{'visodbl01.set','visodbl02.set'}. Otherwise, this
input should be a filename template (i.e., a string with
# where the subject number should be--'visodbl#.set').
If you provide a template, you must use the option
'sub_ids' (see below) to specify the ID numbers of
each subject. Include the files' path unless the files
are in the current working directory.
time_wind - [pair of integers] The time window (in milliseconds)
across which to average EEG amplitude for the response
variable (e.g., [300 500]).
Optional Inputs:
predictors - [cell array of strings] A cell array that specifies the
non-intercept predictors you wish to include in the
regression analysis (an intercept predictor is automatically
added). The predictors must be fields of the EEG.epoch
field of the EEG struct variable stored in each
set file. If not specified, the list of possible
predictors will be presented on the command line and the
user will be asked to select from that list. {default:
not specified}
sub_ids - [integer vector] A set of integers specifying the
subject ID numbers to include in the grand average.
Only necessary if a filename template is given as
the input to gui_infiles_or_tmplt. Leading 0's are not
necessary for ID numbers less than 10 (e.g., If 7 is a
sub_id and visodbl#.set the template, the function will
look for visodbl07.set AND visodbl7.set).
alpha - [value between 0 and 1] Family-wise alpha level (i.e.,
alpha level with adjustment for multiple comparisons)
{default: 0.05}
use_bins - [integer vector] The bin number or numbers to use in the
analysis (i.e., all trials that fall in those bins will
be included in the analysis). {default: all bins}
bsln_wind - [pair of integers] The time window (in units of
milliseconds) that will be used to baseline each trial
of EEG. The mean voltage in the time window will be
removed from the entire EEG epoch. If the data have
already been baselined and ICA artifact correction was NOT
applied, you won't need to re-baseline the data. {default:
no additional baselining will be done}
plot_topos - [0 or 1] If 1, topographies of mean regression
coefficients and their corresponding t-scores will be
produced (if more than one channel was used in the analysis).
{default: 1}
tail - [1 | 0 | -1] If tail=1, the alternative hypothesis is
that mean regression coeffcients are greater than 0 (upper
tailed test). If tail=0, the alternative hypothesis is
that the mean regression coefficients are different than
0 (two tailed test). If tail=-1, the alternative
hypothesis is that the mean regression coefficients
are less than 0 (lower tailed test). {default: 0}
exclude_chans - A cell array of channel labels to exclude from the
analysis (e.g., {'A2','lle','rhe'}). This option
sacrifices spatial resolution to increase test power by
reducing the number of comparisons. 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 analysis
(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.
You cannot use both this option and 'exclude_chans'
(above). {default: not used, all channels included in test}
correction - The method used to correct for multiple comparisons
(i.e., tests at multiple channels). There are two options:
'perm' -within subject permutation test {default}
'Bonf-Holm' -Bonferroni-Holm procedure (faster than but
probably not as powerful as permutation test
method when number of participants is 7 or more).
'FDR BH' -Benjamini & Hochberg (1995) false
discovery rate control procedure.
Guaranteed to control FDR for independent
or positive regression dependent data.
'FDR BY' -Benjamini & Yekutieli (2001) false
discovery rate control procedure.
Always guaranteed to control FDR.
If the regression analysis is being performed on a
single channel, no correction for multiple comparisons
is applied.
n_perm - Number of permutations to use for multiple comparison
correction. 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. Has no
effect if a permutation test is NOT used to correct
for multiple comparisons. {default: 5000}
output_file - A string indicating the name of a tab delimited text
file to produce containing the results of the analysis.
freq_filt - [low_boundary high_boundary] a two element vector indicating the frequency
cut-offs for a 3rd order Butterworth filter that will be applied to each
trial of data. If low_boundary=0, the filter is a low pass filter. If
high_boundary=(sampling rate)/2, then the filter is a high pass filter. If both
boundaries are between 0 and (sampling rate)/2, then the filter is a band-pass filter.
If both boundaries are between 0 and -(sampling rate)/2, then the filter is a band-
stop filter (with boundaries equal to the absolute values of the low_boundary and
high_boundary). Note, using this option requires the signal processing toolbox
function butter.m and filtfilt.m. If 'freq_filt' is not
specified, no fitering is done. {default: no filtering}
n_pad - [integer] number of time points to pad each epoch of
data before filtering. For example, if n_pad is 512 and
you have 256 time points of data per epoch, each epoch of
data will be grown to 1280 time points (i.e., 512*2+256) by
adding 512 padding points before and after the 256 time points
of data. The values of the padding points form a line from the
value of the first and last data point (i.e., if you wrapped
the padded waveform around a circle the pre-padding and
post-padding points would form a line). This method of
padding is the way the Kutaslab function filter.c works.
If n_pad is not specified, butterfilt.m will follow the Kutaslab
convention of n_pad=2*(the number of time points per
epoch). If you wish to avoid zero padding, set n_pad to
0. Note, 'n_pad' has no effect if the optional argument
'freq_filt' is not used.
seed_state - The seed for the MATLAB random number generator. This
will be used to seed the permutation test (if you've
requested to use a permutation test to correct for
multiple comparisons). This argument is useful for
reproducing previous test results. If you pass
regress_results.perm_test_info.seed_state from a previous
run of rm_regress and keep all the other input parameters
the same, you should exactly reproduce the results of the previous run.
Output:
Here's a brief explanation as to what each field of regress_results contains:
-sub_coefs: The regression coefficient at each channel for each
predictor and subject (the matrix is channel x predictor x
subject)
-mn_coefs: The mean regression coefficient across subjects for each
channel and predictor
-se_coefs: The standard error of the mean regression coefficient across
subjects for each channel and predictor
-hi_ci: Upper bound of 95% confidence interval for the estimate of the
mean regression coefficient across subjects for each channel
and predictor. Confidence intervals assume a t-distribution
and are NOT corrected for multiple comparisons.
-low_ci: Upper bound of 95% confidence interval for the estimate of the
mean regression coefficient across subjects for each channel
and predictor.
-t_coefs: t-scores of mean coefficients (i.e., the mean divided by the
standard error)
-t_df: degrees of freedom for the t-scores in t_coefs
-d_coefs: Cohen's d of the mean coefficients (i.e., the mean of the
coefficients divided by the standard deviation). Cohen's
d is a standard measure of effect size.
-tail: The tail of the hypothesis tests
-p_coefs_uncor: Uncorrected p-values for the mean coefficients
-crctn_method: The method used to correct for multiple comparisons
(if more than one channel analyzed)
-p_coefs_cor: p-values after correction for multiple comparisons
-h_coefs_cor: Binary matrix (channel x predictor) indicating whether
or not the null hypothesis of a mean coefficient of 0
could be rejected. 1=the null hypothesis is rejected.
0=you failed to reject the null hypothesis.
-predictors: The names of the predictor variables
-R_collinear: Estimate of the degree of collinearity for each predictor
(see next section)
-family_alpha: Specified family-wise alpha level
-chanlocs: Channel labels and location information for the EEG
channels included in the analysis
-bins: The bins included in the analysis
-time_wind: The lower and upper bound of time window used for analysis
-participant_filenames: Names of the EEGLAB set files used in the analysis.
-perm_test_info: Information about the permutation test used to correct for
multiple comparisons (if used at all). This
includes the number of permutations (n_perm), the
critical t-score for signficance (tmx_ptile), the
seed state of the random number generator (useful
for reproducing the analysis), and the estimated
family-wise alpha level of the test.
References:
Lorch, R.F.,Jr, & Myers, J.L. 1990, Regression analyses of repeated measures
data in cognitive research. Journal of Experimental Psychology: Learning,
Memory, and Cognition, 16(1), 149-157.
Manly, B.F.J. (1997) Randomization, bootstrap, and Monte Carlo methods in
biology. 2nd ed. Chapmn and Hall, London.
Notes:
-To plot topographies of coefficients from the output of this function at
a later time, use the function sig_topo_coefs.m.
-If ICA has been applied to the data and some ICs have been labeled as
artifacts, this function will remove them and re-baseline the data before
performing the regression analysis. The function should report on the
command line how many ICs have been removed.
-If you get the message "Warning: Matrix is singular to working
precision." something is wrong with your data or predictors. One or both
of them is exhibiting perfect collinearity and the regression
coefficients are undetermined.
-If you get a warning that some trials could not be used because one of
predictors have NaN values for that trial, that means that when you
loaded information about that predictor into the set file a predictor
value for that trial was not found. "NaN" simply stands for "not a
number" and indicates that that trial can not be used in the regression
analysis.
Author:
David Groppe
Kutaslab, 7/20100001 % rm_regress() - Performs repeated measures multiple regression analysis 0002 % using EEG as a response variable via the regression method 0003 % described in Lorch & Meyers (1990). The method works by 0004 % performing a separate regression analysis for each 0005 % participant. Then the regression coefficients are 0006 % averaged across participants and their significance 0007 % determined by performing a t-test with a null hypothesis 0008 % of a coefficient of 0. The regression analysis can be 0009 % performed across multiple electrodes. When multiple 0010 % electrodes are used, the function corrects for multiple 0011 % comparisons by using the Bonferroni-Holms procedure. Each 0012 % electrode is considered a distinct test (i.e., the total 0013 % number of electrodes equals the number of tests in the 0014 % family of tests). 0015 % 0016 % Usage: 0017 % >> regress_results=rm_regress(gui_infiles_or_tmplt,time_wind,varargin); 0018 % 0019 % Required Inputs: 0020 % gui_infiles_or_tmplt - ['gui', a cell array of strings, or a string template] 0021 % If 'gui', a GUI is created that allows you to select 0022 % which set files to import (this is probably the 0023 % easiest way to import files). If a cell array of 0024 % of strings, each element of the cell array should 0025 % contain the filename of an EEGLAB set file (e.g., 0026 % {'visodbl01.set','visodbl02.set'}. Otherwise, this 0027 % input should be a filename template (i.e., a string with 0028 % # where the subject number should be--'visodbl#.set'). 0029 % If you provide a template, you must use the option 0030 % 'sub_ids' (see below) to specify the ID numbers of 0031 % each subject. Include the files' path unless the files 0032 % are in the current working directory. 0033 % time_wind - [pair of integers] The time window (in milliseconds) 0034 % across which to average EEG amplitude for the response 0035 % variable (e.g., [300 500]). 0036 % 0037 % Optional Inputs: 0038 % predictors - [cell array of strings] A cell array that specifies the 0039 % non-intercept predictors you wish to include in the 0040 % regression analysis (an intercept predictor is automatically 0041 % added). The predictors must be fields of the EEG.epoch 0042 % field of the EEG struct variable stored in each 0043 % set file. If not specified, the list of possible 0044 % predictors will be presented on the command line and the 0045 % user will be asked to select from that list. {default: 0046 % not specified} 0047 % sub_ids - [integer vector] A set of integers specifying the 0048 % subject ID numbers to include in the grand average. 0049 % Only necessary if a filename template is given as 0050 % the input to gui_infiles_or_tmplt. Leading 0's are not 0051 % necessary for ID numbers less than 10 (e.g., If 7 is a 0052 % sub_id and visodbl#.set the template, the function will 0053 % look for visodbl07.set AND visodbl7.set). 0054 % alpha - [value between 0 and 1] Family-wise alpha level (i.e., 0055 % alpha level with adjustment for multiple comparisons) 0056 % {default: 0.05} 0057 % use_bins - [integer vector] The bin number or numbers to use in the 0058 % analysis (i.e., all trials that fall in those bins will 0059 % be included in the analysis). {default: all bins} 0060 % bsln_wind - [pair of integers] The time window (in units of 0061 % milliseconds) that will be used to baseline each trial 0062 % of EEG. The mean voltage in the time window will be 0063 % removed from the entire EEG epoch. If the data have 0064 % already been baselined and ICA artifact correction was NOT 0065 % applied, you won't need to re-baseline the data. {default: 0066 % no additional baselining will be done} 0067 % plot_topos - [0 or 1] If 1, topographies of mean regression 0068 % coefficients and their corresponding t-scores will be 0069 % produced (if more than one channel was used in the analysis). 0070 % {default: 1} 0071 % tail - [1 | 0 | -1] If tail=1, the alternative hypothesis is 0072 % that mean regression coeffcients are greater than 0 (upper 0073 % tailed test). If tail=0, the alternative hypothesis is 0074 % that the mean regression coefficients are different than 0075 % 0 (two tailed test). If tail=-1, the alternative 0076 % hypothesis is that the mean regression coefficients 0077 % are less than 0 (lower tailed test). {default: 0} 0078 % exclude_chans - A cell array of channel labels to exclude from the 0079 % analysis (e.g., {'A2','lle','rhe'}). This option 0080 % sacrifices spatial resolution to increase test power by 0081 % reducing the number of comparisons. You cannot 0082 % use both this option and 'include_chans' (below).{default: 0083 % not used, all channels included in test} 0084 % include_chans - A cell array of channel labels to use in the analysis 0085 % (e.g., {'A2','lle','rhe'}). All other channels will 0086 % be ignored. This option sacrifices spatial resolution to 0087 % increase test power by reducing the number of comparisons. 0088 % You cannot use both this option and 'exclude_chans' 0089 % (above). {default: not used, all channels included in test} 0090 % correction - The method used to correct for multiple comparisons 0091 % (i.e., tests at multiple channels). There are two options: 0092 % 'perm' -within subject permutation test {default} 0093 % 'Bonf-Holm' -Bonferroni-Holm procedure (faster than but 0094 % probably not as powerful as permutation test 0095 % method when number of participants is 7 or more). 0096 % 'FDR BH' -Benjamini & Hochberg (1995) false 0097 % discovery rate control procedure. 0098 % Guaranteed to control FDR for independent 0099 % or positive regression dependent data. 0100 % 'FDR BY' -Benjamini & Yekutieli (2001) false 0101 % discovery rate control procedure. 0102 % Always guaranteed to control FDR. 0103 % If the regression analysis is being performed on a 0104 % single channel, no correction for multiple comparisons 0105 % is applied. 0106 % n_perm - Number of permutations to use for multiple comparison 0107 % correction. Manly (1997) suggests using at least 1000 0108 % permutations for an alpha level of 0.05 and at least 0109 % 5000 permutations for an alpha level of 0.01. Has no 0110 % effect if a permutation test is NOT used to correct 0111 % for multiple comparisons. {default: 5000} 0112 % output_file - A string indicating the name of a tab delimited text 0113 % file to produce containing the results of the analysis. 0114 % freq_filt - [low_boundary high_boundary] a two element vector indicating the frequency 0115 % cut-offs for a 3rd order Butterworth filter that will be applied to each 0116 % trial of data. If low_boundary=0, the filter is a low pass filter. If 0117 % high_boundary=(sampling rate)/2, then the filter is a high pass filter. If both 0118 % boundaries are between 0 and (sampling rate)/2, then the filter is a band-pass filter. 0119 % If both boundaries are between 0 and -(sampling rate)/2, then the filter is a band- 0120 % stop filter (with boundaries equal to the absolute values of the low_boundary and 0121 % high_boundary). Note, using this option requires the signal processing toolbox 0122 % function butter.m and filtfilt.m. If 'freq_filt' is not 0123 % specified, no fitering is done. {default: no filtering} 0124 % n_pad - [integer] number of time points to pad each epoch of 0125 % data before filtering. For example, if n_pad is 512 and 0126 % you have 256 time points of data per epoch, each epoch of 0127 % data will be grown to 1280 time points (i.e., 512*2+256) by 0128 % adding 512 padding points before and after the 256 time points 0129 % of data. The values of the padding points form a line from the 0130 % value of the first and last data point (i.e., if you wrapped 0131 % the padded waveform around a circle the pre-padding and 0132 % post-padding points would form a line). This method of 0133 % padding is the way the Kutaslab function filter.c works. 0134 % If n_pad is not specified, butterfilt.m will follow the Kutaslab 0135 % convention of n_pad=2*(the number of time points per 0136 % epoch). If you wish to avoid zero padding, set n_pad to 0137 % 0. Note, 'n_pad' has no effect if the optional argument 0138 % 'freq_filt' is not used. 0139 % seed_state - The seed for the MATLAB random number generator. This 0140 % will be used to seed the permutation test (if you've 0141 % requested to use a permutation test to correct for 0142 % multiple comparisons). This argument is useful for 0143 % reproducing previous test results. If you pass 0144 % regress_results.perm_test_info.seed_state from a previous 0145 % run of rm_regress and keep all the other input parameters 0146 % the same, you should exactly reproduce the results of the previous run. 0147 % 0148 % Output: 0149 % Here's a brief explanation as to what each field of regress_results contains: 0150 % -sub_coefs: The regression coefficient at each channel for each 0151 % predictor and subject (the matrix is channel x predictor x 0152 % subject) 0153 % -mn_coefs: The mean regression coefficient across subjects for each 0154 % channel and predictor 0155 % -se_coefs: The standard error of the mean regression coefficient across 0156 % subjects for each channel and predictor 0157 % -hi_ci: Upper bound of 95% confidence interval for the estimate of the 0158 % mean regression coefficient across subjects for each channel 0159 % and predictor. Confidence intervals assume a t-distribution 0160 % and are NOT corrected for multiple comparisons. 0161 % -low_ci: Upper bound of 95% confidence interval for the estimate of the 0162 % mean regression coefficient across subjects for each channel 0163 % and predictor. 0164 % -t_coefs: t-scores of mean coefficients (i.e., the mean divided by the 0165 % standard error) 0166 % -t_df: degrees of freedom for the t-scores in t_coefs 0167 % -d_coefs: Cohen's d of the mean coefficients (i.e., the mean of the 0168 % coefficients divided by the standard deviation). Cohen's 0169 % d is a standard measure of effect size. 0170 % -tail: The tail of the hypothesis tests 0171 % -p_coefs_uncor: Uncorrected p-values for the mean coefficients 0172 % -crctn_method: The method used to correct for multiple comparisons 0173 % (if more than one channel analyzed) 0174 % -p_coefs_cor: p-values after correction for multiple comparisons 0175 % -h_coefs_cor: Binary matrix (channel x predictor) indicating whether 0176 % or not the null hypothesis of a mean coefficient of 0 0177 % could be rejected. 1=the null hypothesis is rejected. 0178 % 0=you failed to reject the null hypothesis. 0179 % -predictors: The names of the predictor variables 0180 % -R_collinear: Estimate of the degree of collinearity for each predictor 0181 % (see next section) 0182 % -family_alpha: Specified family-wise alpha level 0183 % -chanlocs: Channel labels and location information for the EEG 0184 % channels included in the analysis 0185 % -bins: The bins included in the analysis 0186 % -time_wind: The lower and upper bound of time window used for analysis 0187 % -participant_filenames: Names of the EEGLAB set files used in the analysis. 0188 % -perm_test_info: Information about the permutation test used to correct for 0189 % multiple comparisons (if used at all). This 0190 % includes the number of permutations (n_perm), the 0191 % critical t-score for signficance (tmx_ptile), the 0192 % seed state of the random number generator (useful 0193 % for reproducing the analysis), and the estimated 0194 % family-wise alpha level of the test. 0195 % 0196 % References: 0197 % Lorch, R.F.,Jr, & Myers, J.L. 1990, Regression analyses of repeated measures 0198 % data in cognitive research. Journal of Experimental Psychology: Learning, 0199 % Memory, and Cognition, 16(1), 149-157. 0200 % 0201 % Manly, B.F.J. (1997) Randomization, bootstrap, and Monte Carlo methods in 0202 % biology. 2nd ed. Chapmn and Hall, London. 0203 % 0204 % Notes: 0205 % -To plot topographies of coefficients from the output of this function at 0206 % a later time, use the function sig_topo_coefs.m. 0207 % -If ICA has been applied to the data and some ICs have been labeled as 0208 % artifacts, this function will remove them and re-baseline the data before 0209 % performing the regression analysis. The function should report on the 0210 % command line how many ICs have been removed. 0211 % -If you get the message "Warning: Matrix is singular to working 0212 % precision." something is wrong with your data or predictors. One or both 0213 % of them is exhibiting perfect collinearity and the regression 0214 % coefficients are undetermined. 0215 % -If you get a warning that some trials could not be used because one of 0216 % predictors have NaN values for that trial, that means that when you 0217 % loaded information about that predictor into the set file a predictor 0218 % value for that trial was not found. "NaN" simply stands for "not a 0219 % number" and indicates that that trial can not be used in the regression 0220 % analysis. 0221 % 0222 % Author: 0223 % David Groppe 0224 % Kutaslab, 7/2010 0225 % 0226 0227 %%%%%%%%%%%%%%%% Future Work %%%%%%%%%%%%%%%%% 0228 % -diagnostics plots for each sub? 0229 % -add verblevel as input and make more use of it? 0230 % -could speed up function by baselining just used trials 0231 % -option: exclude garv flagged trials? 0232 % -reaction time as predictor should be checked 0233 % -store additional permutation test info in regress_results? 0234 0235 %%%%%%%%%%%%%%%% REVISION LOG %%%%%%%%%%%%%%%%% 0236 % 7/19/2010 Added: 0237 % -check that there's more than one subject with a usable number of trials 0238 % -check that there aren't more predictors than trials 0239 % 0240 % 7/23/2010 0241 % -Added random generator seed state optional input argument 0242 % -Added additional information about test to regress_results and to command 0243 % line/output file output 0244 % -Moved report as to min/max predictor values so that it did not include 0245 % trials excluded from analysis 0246 % 0247 % 10/5/2010-freq_filt and n_pad options added 0248 % 0249 % 10/6/2010-'predictors' argument now works. p-value summary added to 0250 % command line output 0251 % 0252 % 12/9/2010-First trial in selected bins was being used even if predictors 0253 % had an NaN value for that trial. That's been fixed. Code now reports 0254 % how many trials could not be used because they have NaN values. Code 0255 % also checks for empty EEG.epoch(x).eventfield values, warns if they're 0256 % found, and replaces them with NaN 0257 % 0258 % 2/1/2010-remove_art_ics.m replaced by remove_artifact_ics.m. Bugs in help 0259 % documentation fixed. 0260 % 0261 % 6/28/2011-Added FDR control as an option for correcting for multiple 0262 % comparisons 0263 0264 function regress_results=rm_regress_mass(gui_infiles_or_tmplt,time_wind,varargin) 0265 0266 % One optional input arguments has no effect. 0267 % It just included so that items_regress.m and rm_regress.m accept the same optional input arguments 0268 % This allows minF_regress.m to call them with the same varagin. 0269 p=inputParser; 0270 p.addRequired('gui_infiles_or_tmplt',@(x) ischar(x) || iscell(x)); 0271 p.addRequired('time_wind',@(x) isnumeric(x) && (length(x)==2) && isvector(x)); 0272 p.addParamValue('predictors',[],@(x) ischar(x) || iscell(x)); 0273 p.addParamValue('alpha',.05,@(x) isnumeric(x) && (length(x)==1) && isvector(x) && (x>0) && (x<1)); 0274 p.addParamValue('sub_ids',[],@isnumeric); 0275 p.addParamValue('use_bins',[],@(x) isnumeric(x) && isvector(x)); 0276 p.addParamValue('bsln_wind',[],@(x) isnumeric(x) && length(x)==2); 0277 p.addParamValue('freq_filt',[],@(x) isnumeric(x) && length(x)==2); 0278 p.addParamValue('n_pad',[],@(x) isnumeric(x) && length(x)==1); 0279 p.addParamValue('plot_raster',1,@isnumeric); 0280 p.addParamValue('exclude_chans',[],@(x) ischar(x) || iscell(x)); 0281 p.addParamValue('include_chans',[],@(x) ischar(x) || iscell(x)); 0282 p.addParamValue('output_file',[],@ischar); 0283 p.addParamValue('n_perm',5000,@(x) isscalar(x) && (x>0)); 0284 p.addParamValue('tail',0,@(x) isscalar(x) && ((x==0) || (x==-1) || (x==1))); 0285 p.addParamValue('correction','perm',@(x) ischar(x) && (strcmpi(x,'Bonf-Holm') ... 0286 || strcmpi(x,'perm') || strcmpi(x,'FDR BH') || strcmpi(x,'FDR BY'))); 0287 p.addParamValue('seed_state',[],@isnumeric); 0288 p.addParamValue('use_bins_and_evcode','no',@(x) ischar(x) && (strcmpi(x,'yes')) || (strcmpi(x,'no'))); %Has no effect. 0289 0290 p.parse(gui_infiles_or_tmplt,time_wind,varargin{:}); 0291 0292 0293 %Make sure exclude & include options were not both used. 0294 if ~isempty(p.Results.include_chans) && ~isempty(p.Results.exclude_chans) 0295 error('You cannot use BOTH ''include_chans'' and ''exclude_chans'' options.'); 0296 end 0297 if ischar(p.Results.exclude_chans), 0298 exclude_chans{1}=p.Results.exclude_chans; 0299 elseif isempty(p.Results.exclude_chans) 0300 exclude_chans=[]; 0301 else 0302 exclude_chans=p.Results.exclude_chans; 0303 end 0304 if ischar(p.Results.include_chans), 0305 include_chans{1}=p.Results.include_chans; 0306 elseif isempty(p.Results.include_chans) 0307 include_chans=[]; 0308 else 0309 include_chans=p.Results.include_chans; 0310 end 0311 flt=p.Results.freq_filt; 0312 0313 VERBLEVEL=2; % ?? This is currently barely used. 0314 global EEG; %necessary for remove_artifact_ics.m 0315 0316 %% Get *.set fileames 0317 infiles=get_set_infiles(gui_infiles_or_tmplt,p.Results.sub_ids); 0318 0319 %% Import Data 0320 n_sub=length(infiles); 0321 if n_sub==1, 0322 error('You need more than one participant for repeated measures regression.'); 0323 end 0324 0325 n_regress_trial=zeros(1,n_sub); 0326 for sub=1:n_sub, 0327 fprintf('Loading data for Participant %s\n',infiles{sub}); 0328 EEG=pop_loadset(infiles{sub}); 0329 %load(infiles{sub},'-MAT'); 0330 0331 % Get time points 0332 if time_wind(1)<EEG.times(1), 0333 error('Analysis time window (%d to %d ms) begins before earliest time point in file %s (%d ms).', ... 0334 time_wind(1),time_wind(2),infiles{sub},EEG.times(1)); 0335 else 0336 start_tpt=find_tpt(time_wind(1),EEG.times); 0337 end 0338 if time_wind(2)>EEG.times(end), 0339 error('Analysis time window (%d to %d ms) extends beyond latest time point in file %s (%d ms).', ... 0340 time_wind(1),time_wind(2),infiles{sub},EEG.times(end)); 0341 else 0342 end_tpt=find_tpt(time_wind(2),EEG.times); 0343 end 0344 time_ids=start_tpt:end_tpt; 0345 n_tpts=length(time_ids); 0346 if n_tpts<=1, 0347 error('This function requires a time window of interest that spans more than one time point.'); 0348 end 0349 0350 %Check for artifact ICs and remove them if there 0351 %[EEG, ICs_removed]=remove_art_ics(EEG,p.Results.bsln_wind,VERBLEVEL);?? old code 0352 ICs_removed=remove_artifact_ics(p.Results.bsln_wind,VERBLEVEL); 0353 0354 %Filter the data if requested 0355 if ~isempty(flt) 0356 %first error check 0357 if ~isempty(flt), 0358 if max(flt)>(EEG.srate/2), 0359 error('''filt'' parameters need to be less than or equal to (sampling rate)/2 (i.e., %f).',EEG.srate/2); 0360 elseif (flt(2)==(EEG.srate/2)) && (flt(1)==0), 0361 error('If second element of ''filt'' parameter is EEG.srate/2, then the first element must be greater than 0.'); 0362 elseif abs(flt(2))<=abs(flt(1)), 0363 error('Second element of ''filt'' parameters must be greater than first in absolute value.'); 0364 elseif (flt(1)<0) || (flt(2)<0), 0365 if (flt(1)>=0) || (flt(2)>=0), 0366 error('BOTH parameters of ''filt'' need to be greater than or equal to zero OR need to be negative.'); 0367 end 0368 if min(flt)<=(-EEG.srate/2), 0369 error('''filt'' parameters need to be greater than sampling rate/2 (i.e., -%f) when creating a stop band.',EEG.srate/2); 0370 end 0371 end 0372 0373 end 0374 EEG=butterfilt(EEG,flt,[],p.Results.n_pad); 0375 end 0376 0377 %Baseline EEG (unless already baselined by remove_art_ics 0378 if ~isempty(p.Results.bsln_wind) && isempty(ICs_removed), 0379 EEG=baseline_EEG(EEG,p.Results.bsln_wind,VERBLEVEL); 0380 end 0381 0382 n_pos_chan=EEG.nbchan; 0383 chan_ids=zeros(1,n_pos_chan); 0384 0385 if sub==1, 0386 %% Identify predictors 0387 if isempty(p.Results.predictors), 0388 %Print list of possible predictors on command line 0389 epoch_flds=fieldnames(EEG.epoch); 0390 n_epoch_flds=length(epoch_flds); 0391 poss_pred=cell(1,n_epoch_flds); 0392 poss_pred_fullname=cell(1,n_epoch_flds); 0393 0394 fprintf('\n**** Choose Predictor(s) for Regression Analysis ****\n'); 0395 ct=0; 0396 for fld_loop=1:n_epoch_flds, 0397 if ~strcmpi('event',epoch_flds{fld_loop}) && ~strcmpi('eventtype',epoch_flds{fld_loop}), 0398 ct=ct+1; 0399 poss_pred_fullname{ct}=epoch_flds{fld_loop}; 0400 if strcmpi('event',epoch_flds{fld_loop}(1:5)), 0401 poss_pred{ct}=epoch_flds{fld_loop}(6:end); 0402 else 0403 poss_pred{ct}=epoch_flds{fld_loop}; 0404 end 0405 fprintf('%d) %s\n',ct,poss_pred{ct}); 0406 end 0407 end 0408 %Is reaction time a field? 0409 EEG_flds=fieldnames(EEG); 0410 if ismember('rtmsec',EEG_flds) && ismember('rtbins',EEG_flds), 0411 ct=ct+1; 0412 poss_pred{ct}='Reaction Time'; 0413 poss_pred_fullname{ct}='Reaction Time'; 0414 fprintf('%d) %s\n',ct,poss_pred{ct}); 0415 end 0416 0417 pred_ids=[]; 0418 while isempty(pred_ids), 0419 pred_ids=input('Enter the number(s) of desired predictor(s): ','s'); 0420 fprintf('\n'); 0421 pred_ids=str2num(pred_ids); 0422 if max(pred_ids)>ct, 0423 fprintf('ERROR: There are no predictor indices greater than %d.\n',ct); 0424 pred_ids=[]; 0425 elseif min(pred_ids)<1, 0426 fprintf('ERROR: All predictor indices must be greater than 0.\n'); 0427 pred_ids=[]; 0428 end 0429 end 0430 predictors=poss_pred_fullname(unique(pred_ids)); 0431 predictor_names_neat=poss_pred(unique(pred_ids)); 0432 if ismember('Reaction Time',predictors), 0433 rt_selected=1; 0434 else 0435 rt_selected=0; 0436 end 0437 else 0438 %get possible predictors 0439 epoch_flds=fieldnames(EEG.epoch); 0440 n_epoch_flds=length(epoch_flds); 0441 poss_pred=cell(1,n_epoch_flds); 0442 poss_pred_fullname=cell(1,n_epoch_flds); 0443 ct=0; 0444 for fld_loop=1:n_epoch_flds, 0445 if ~strcmpi('event',epoch_flds{fld_loop}) && ~strcmpi('eventtype',epoch_flds{fld_loop}), 0446 ct=ct+1; 0447 poss_pred_fullname{ct}=epoch_flds{fld_loop}; 0448 if strcmpi('event',epoch_flds{fld_loop}(1:5)), 0449 poss_pred{ct}=epoch_flds{fld_loop}(6:end); 0450 else 0451 poss_pred{ct}=epoch_flds{fld_loop}; 0452 end 0453 end 0454 end 0455 %Is reaction time a field? 0456 EEG_flds=fieldnames(EEG); 0457 if ismember('rtmsec',EEG_flds) && ismember('rtbins',EEG_flds), 0458 ct=ct+1; 0459 poss_pred{ct}='Reaction Time'; 0460 poss_pred_fullname{ct}='Reaction Time'; 0461 end 0462 0463 %figure out which of the possible predictors have been 0464 %requested 0465 if ischar(p.Results.predictors) 0466 requested_pred{1}=p.Results.predictors; 0467 else 0468 requested_pred=p.Results.predictors; 0469 end 0470 pred_ids=zeros(1,length(requested_pred)); 0471 ct=0; 0472 n_poss=length(poss_pred); 0473 for dloop=1:length(requested_pred), 0474 found=0; 0475 for eloop=1:n_poss, 0476 if strcmpi(requested_pred{dloop},poss_pred{eloop}) || ... 0477 strcmpi(requested_pred{dloop},poss_pred_fullname{eloop}) 0478 found=1; 0479 ct=ct+1; 0480 pred_ids(ct)=eloop; 0481 break; 0482 end 0483 end 0484 if ~found, 0485 fprintf('*** ERROR!! ***\n'); 0486 fprintf('Requested predictor %s was not found in this EEG variable.\n', ... 0487 requested_pred{dloop}); 0488 fprintf('Possible predictors are:\n'); 0489 for floop=1:n_poss, 0490 if ~isempty(poss_pred{floop}), 0491 fprintf('%d) %s\n',floop,poss_pred{floop}); 0492 end 0493 end 0494 error('Requested predictor not found.'); 0495 end 0496 end 0497 predictors=poss_pred_fullname(unique(pred_ids)); 0498 predictor_names_neat=poss_pred(unique(pred_ids)); 0499 if ismember('Reaction Time',predictors), 0500 rt_selected=1; 0501 else 0502 rt_selected=0; 0503 end 0504 end 0505 n_pred=length(predictors); 0506 0507 %% Figure out which channels to ignore if any 0508 % exclude and include chan options 0509 use_chan_labels=cell(1,n_pos_chan); %preallocate mem 0510 ct=0; 0511 if ~isempty(exclude_chans), 0512 for x=1:n_pos_chan, 0513 if ~ismember(EEG.chanlocs(x).labels,exclude_chans), 0514 ct=ct+1; 0515 chan_ids(ct)=x; 0516 use_chan_labels{ct}=EEG.chanlocs(x).labels; 0517 end 0518 end 0519 elseif ~isempty(include_chans), 0520 for x=1:n_pos_chan, 0521 if ismember(EEG.chanlocs(x).labels,include_chans), 0522 ct=ct+1; 0523 chan_ids(ct)=x; 0524 use_chan_labels{ct}=EEG.chanlocs(x).labels; 0525 end 0526 end 0527 else 0528 for x=1:n_pos_chan, 0529 use_chan_labels{x}=EEG.chanlocs(x).labels; 0530 end 0531 chan_ids=1:n_pos_chan; 0532 ct=n_pos_chan; 0533 end 0534 chan_ids=chan_ids(1:ct); 0535 use_chan_labels(ct+1:end)=[]; 0536 n_chan=length(chan_ids); 0537 chanlocs=EEG.chanlocs(chan_ids); 0538 0539 %Preallocate memory for regression coefficient estimates 0540 sub_coefs=zeros(n_chan,n_tpts,n_pred+1,n_sub); %add one for intercept 0541 R_collinear=zeros(n_pred+1,n_sub)*NaN; %add one for intercept 0542 else 0543 %Not the first subject, grab the same channels in the same order and throw error if missing used channel(s) 0544 ct=0; 0545 for a=1:n_chan, 0546 found=0; 0547 for b=1:n_pos_chan, 0548 %a=channel index for analysis 0549 %b=channel index for this particular subject 0550 if strcmpi(EEG.chanlocs(b).labels,use_chan_labels{a}), 0551 found=1; 0552 ct=ct+1; 0553 chan_ids(ct)=b; 0554 break; 0555 end 0556 end 0557 if ~found, 0558 error('Missing channel %s for participant %s.',use_chan_labels{a}, ... 0559 infiles{sub}); 0560 end 0561 end 0562 chan_ids=chan_ids(1:ct); 0563 end 0564 0565 poss_prdct=fieldnames(EEG.epoch); 0566 if rt_selected, 0567 EEG_flds=fieldnames(EEG); 0568 if ismember('rtmsec',EEG_flds) && ismember('rtbins',EEG_flds), 0569 poss_prdct{1+length(poss_prdct)}='Reaction Time'; 0570 end 0571 end 0572 missing_prdct=setdiff(predictors,poss_prdct); 0573 if ~isempty(missing_prdct), 0574 fprintf('ERROR: File %s is missing the following predictor(s): \n', ... 0575 infiles{sub}); 0576 for a=1:length(missing_prdct), 0577 fprintf('%s \n',missing_prdct{a}); 0578 end 0579 error('Regression aborted.'); 0580 end 0581 0582 n_trial=length(EEG.epoch); 0583 use_trial=zeros(1,n_trial); 0584 pred_val=zeros(n_pred,n_trial); 0585 EEG_val=zeros(n_chan,n_tpts,n_trial); 0586 for trial=1:n_trial, 0587 n_type=length(EEG.epoch(trial).eventtype); 0588 0589 if isempty(p.Results.use_bins), 0590 use_trial(trial)=1; 0591 else 0592 %Find out if current trial falls into desired bin 0593 for typ=1:n_type, 0594 for bin=p.Results.use_bins, 0595 if strcmpi(['bin' num2str(bin)],EEG.epoch(trial).eventtype{typ}), 0596 use_trial(trial)=1; 0597 break; %leave bin loop 0598 end 0599 end 0600 if use_trial(trial), 0601 break; %leave typ loop 0602 end 0603 end 0604 end 0605 0606 %Record predictor and response var values if trial is being used 0607 if use_trial(trial), 0608 for pred=1:n_pred, 0609 if strcmpi(predictors{pred},'Reaction Time'), 0610 n_rt_bin=length(EEG.rtbins{trial}); 0611 bin_rts=zeros(1,n_rt_bin); 0612 bin_ct=0; 0613 for b=1:n_rt_bin, 0614 if ismember(EEG.rtbins{trial}(b),p.Results.use_bins), 0615 bin_ct=bin_ct+1; 0616 bin_rts(bin_ct)=EEG.rtmsec{trial}(b); 0617 end 0618 end 0619 bin_rts=bin_rts(1:bin_ct); 0620 uni_rt=unique(bin_rts); 0621 if length(uni_rt)>1 0622 error('Trial %d has more than one reaction time for the bins you''ve selected.',trial); 0623 elseif isempty(uni_rt), 0624 watchit(sprintf('Epoch %d has an empty value for reaction time. Please use NaN for missing predictor values.\n', ... 0625 trial)); 0626 pred_val(pred,trial)=NaN; 0627 else 0628 pred_val(pred,trial)=uni_rt; 0629 end 0630 else 0631 cmnd=['emptyval=isempty(EEG.epoch(trial).' predictors{pred} ');']; 0632 eval(cmnd); %for some reason calling isempty.m wasn't working in the function so I used this workaround 0633 if emptyval, 0634 watchit(sprintf('Epoch %d has an empty value for predictor %s. Please use NaN for missing predictor values.\n', ... 0635 trial,predictor_names_neat{pred})); 0636 cmnd=['pred_val(pred,trial)=NaN;']; %?? Do RT too 0637 eval(cmnd); 0638 else 0639 cmnd=['temp_pred=EEG.epoch(trial).' predictors{pred} ';']; 0640 eval(cmnd); 0641 if iscell(temp_pred), 0642 if all_same_cell_elements(temp_pred), 0643 pred_val(pred,trial)=temp_pred{1}; 0644 else 0645 error('Epoch %d has multiple values for predictor %s.', ... 0646 trial,predictors{pred}); 0647 end 0648 else 0649 pred_val(pred,trial)=temp_pred; 0650 end 0651 end 0652 end 0653 end 0654 for c=1:n_chan, 0655 EEG_val(c,:,trial)=EEG.data(chan_ids(c),time_ids,trial); 0656 end 0657 end 0658 end 0659 0660 %get rid of unused trials 0661 trial_ids=find(use_trial==1); 0662 if n_chan>1, 0663 EEG_val=EEG_val(:,:,trial_ids); 0664 else 0665 EEG_val=EEG_val(:,trial_ids); 0666 end 0667 if n_pred>1, 0668 pred_val=pred_val(:,trial_ids); 0669 else 0670 pred_val=pred_val(trial_ids); 0671 end 0672 0673 %get rid of NaN trials (if any) 0674 if n_pred>1, 0675 has_NaN=sum(isnan(pred_val)); 0676 non_NaN_ids=find(has_NaN==0); 0677 NaN_ids=find(has_NaN); 0678 pred_val=pred_val(:,non_NaN_ids); 0679 else 0680 non_NaN_ids=find(~isnan(pred_val)); 0681 NaN_ids=find(isnan(pred_val)); 0682 pred_val=pred_val(non_NaN_ids); 0683 end 0684 if n_chan>1, 0685 EEG_val=EEG_val(:,:,non_NaN_ids); 0686 else 0687 EEG_val=EEG_val(:,non_NaN_ids); 0688 end 0689 0690 n_nonNaN_trials=length(non_NaN_ids); 0691 n_selected_trials=length(trial_ids); 0692 if n_selected_trials~=n_nonNaN_trials, 0693 watchit(sprintf('%d of the selected trials have NaN predictor values for at least one predictor variable.', ... 0694 n_selected_trials-n_nonNaN_trials)); 0695 NaN_trial_ids=trial_ids(NaN_ids); 0696 if VERBLEVEL>2, 0697 disp(['The EEG.epoch indices of those trials are:' num2str(NaN_trial_ids)]); 0698 end 0699 end 0700 n_regress_trial(sub)=n_nonNaN_trials; 0701 fprintf('Number of usable trials: %d\n',n_nonNaN_trials); 0702 %Check if sufficient # of usable trials 0703 if ~n_nonNaN_trials, 0704 if isempty(p.Results.use_bins), 0705 error('Participant %s has no usable trials. Regression aborted.', ... 0706 infiles{sub}); 0707 else 0708 error('Participant %s has no trials for Bin(s) %s. Regression aborted.', ... 0709 infiles{sub},num2str(p.Results.use_bins)); 0710 end 0711 elseif n_nonNaN_trials<=n_pred, 0712 error('Participant %s has insufficient usable trials (i.e., %d) for the number of predictors in this analysis (i.e., %d plus an intercept term).', ... 0713 infiles{sub},n_nonNaN_trials,n_pred); 0714 end 0715 0716 %Report min and max values for each predictor 0717 for pred_loop=1:n_pred, 0718 fprintf('%s (min/max): %f/%f\n',predictor_names_neat{pred_loop}, ... 0719 min(pred_val(pred_loop,:)),max(pred_val(pred_loop,:))); 0720 end 0721 0722 %perform multiple regression at each channel 0723 for c=1:n_chan, 0724 for tp=1:n_tpts, 0725 [coeff, S_err, XTXI, R_sq, R_sq_adj, F_val, Coef_stats, y_hat, residuals]= ... 0726 mregress(squeeze(EEG_val(c,tp,:)),pred_val',1); %final argument "1" means that an intercept term is included in the model 0727 %Note that the first element of coeff and associated variables 0728 %correpsponds to the intercept term 0729 sub_coefs(c,tp,1:n_pred+1,sub)=coeff; 0730 end 0731 end 0732 if sum(isnan(coeff)), 0733 error('The data for participant %s exhibits perfect collinearity and the regression coefficients cannot be determined.',infiles{sub}); 0734 end 0735 %Quantify degree of collinearity with other predictors 0736 if n_pred>1, 0737 for a=1:n_pred, 0738 [coeffCOL, S_errCOL, XTXIcol, R_sqCOL]=mregress(pred_val(a,:)', ... 0739 pred_val(setdiff(1:n_pred,a),:)',1); 0740 R_collinear(a+1,sub)=R_sqCOL; 0741 end 0742 end 0743 fprintf('\n'); 0744 end 0745 0746 regress_results.sub_coefs=sub_coefs; 0747 regress_results.mn_coefs=mean(sub_coefs,4); 0748 regress_results.se_coefs=stderr(sub_coefs,4); 0749 t_df=n_sub-1; 0750 crit_t=tinv(.975,t_df); %for 95% confidence intervals 0751 regress_results.hi_ci=regress_results.mn_coefs+crit_t*regress_results.se_coefs; 0752 regress_results.low_ci=regress_results.mn_coefs-crit_t*regress_results.se_coefs; 0753 regress_results.t_coefs=regress_results.mn_coefs./regress_results.se_coefs; 0754 regress_results.t_df=t_df; 0755 regress_results.d_coefs=regress_results.mn_coefs./std(sub_coefs,0,4); %0 indicates that n-1 SD formula will be used 0756 regress_results.F_coefs=zeros(n_chan,n_pred)*NaN; 0757 regress_results.tail=p.Results.tail; 0758 regress_results.p_coefs_uncor=zeros(n_chan,n_tpts,n_pred)*NaN; 0759 if n_chan==1, 0760 regress_results.crctn_method='Not applicable'; 0761 else 0762 regress_results.crctn_method=p.Results.correction; 0763 end 0764 regress_results.p_coefs_cor=regress_results.p_coefs_uncor; 0765 regress_results.h_coefs_cor=regress_results.p_coefs_uncor; 0766 regress_results.predictor_names=cell(1,n_pred+1); 0767 regress_results.predictor_names{1}='Intercept'; 0768 %regress_results.predictor_names(2:n_pred+1)=predictors; <=These are 0769 %predictor names with event in front of it 0770 regress_results.predictor_names(2:n_pred+1)=predictor_names_neat; 0771 regress_results.R_collinear=R_collinear; 0772 regress_results.family_alpha=p.Results.alpha; 0773 regress_results.chanlocs=chanlocs; 0774 if isempty(p.Results.use_bins), 0775 regress_results.bins='All'; 0776 else 0777 regress_results.bins=p.Results.use_bins; 0778 end 0779 regress_results.time_wind=time_wind; 0780 regress_results.participant_filenames=infiles; 0781 regress_results.perm_test_info.n_perm=NaN; 0782 regress_results.perm_test_info.tmx_ptile=NaN; 0783 regress_results.perm_test_info.seed_state=NaN; 0784 regress_results.perm_test_info.est_alpha=NaN; 0785 for pred=1:n_pred+1, 0786 fprintf('Computing significance of mean coefficient(s) for predictor %s.\n', ... 0787 regress_results.predictor_names{pred}); 0788 if ~isempty(p.Results.seed_state) && strcmpi(p.Results.correction,'perm'); 0789 fprintf('Seeding permutation test with specified seed_state.\n'); 0790 end 0791 for c=1:n_chan, 0792 for tp=1:n_tpts, 0793 if ~p.Results.tail 0794 %two-tailed test 0795 regress_results.p_coefs_uncor(c,tp,pred)=2*(1-tcdf(abs(regress_results.t_coefs(c,tp,pred)),n_sub-1)); 0796 elseif p.Results.tail==1, 0797 %upper-tailed test 0798 regress_results.p_coefs_uncor(c,tp,pred)=(1-tcdf(regress_results.t_coefs(c,tp,pred),n_sub-1)); 0799 else 0800 %lower-tailed test 0801 regress_results.p_coefs_uncor(c,tp,pred)=tcdf(regress_results.t_coefs(c,tp,pred),n_sub-1); 0802 end 0803 regress_results.F_coefs(c,tp,pred)=finv(1-regress_results.p_coefs_uncor(c,tp,pred),1,regress_results.t_df); 0804 end 0805 end 0806 if strcmpi(p.Results.correction,'Bonf-Holm'), 0807 [regress_results.p_coefs_cor(:,:,pred), regress_results.h_coefs_cor(:,:,pred)]= ... 0808 bonf_holm(squeeze(regress_results.p_coefs_uncor(:,:,pred)),p.Results.alpha); 0809 elseif strcmpi(p.Results.correction,'FDR BH'), 0810 [regress_results.h_coefs_cor(:,:,pred) , dummy, regress_results.p_coefs_cor(:,:,pred)]= ... 0811 fdr_bh(squeeze(regress_results.p_coefs_uncor(:,:,pred)),p.Results.alpha,'pdep','yes'); 0812 elseif strcmpi(p.Results.correction,'FDR BY'), 0813 [regress_results.h_coefs_cor(:,:,pred) dummy regress_results.p_coefs_cor(:,:,pred)]= ... 0814 fdr_bh(squeeze(regress_results.p_coefs_uncor(:,:,pred)),p.Results.alpha,'dep','yes'); 0815 else 0816 if pred>1 || isempty(p.Results.seed_state), 0817 [pval, t_orig, tmx_ptile, seed_state, est_alpha]=mxt_perm1(squeeze(sub_coefs(:,:,pred,:)), ... 0818 p.Results.n_perm,p.Results.alpha,p.Results.tail,VERBLEVEL); 0819 else 0820 %Just seed the first permutation test if seed_state 0821 %specified 0822 [pval, t_orig, tmx_ptile, seed_state, est_alpha]=mxt_perm1(squeeze(sub_coefs(:,:,pred,:)), ... 0823 p.Results.n_perm,p.Results.alpha,p.Results.tail,VERBLEVEL,p.Results.seed_state); 0824 end 0825 regress_results.p_coefs_cor(:,:,pred)=pval; 0826 regress_results.h_coefs_cor(:,:,pred)=pval<p.Results.alpha; 0827 regress_results.perm_test_info.n_perm=p.Results.n_perm; 0828 regress_results.perm_test_info.tmx_ptile=tmx_ptile; 0829 if pred==1, 0830 regress_results.perm_test_info.seed_state=seed_state; 0831 end 0832 regress_results.perm_test_info.est_alpha=est_alpha; 0833 end 0834 end 0835 0836 %add error check if only one participant submitted? ?? 0837 fprintf('\n\n'); 0838 0839 %% Report Results to Command Line 0840 fprintf('**ANALYSIS PARAMETERS**\n'); 0841 fprintf('Number of participants: %d\n',n_sub); 0842 fprintf('Regression analysis t-score degrees of freedom: %d\n',t_df); 0843 if p.Results.tail>0, 0844 fprintf('Tail of tests: Upper\n'); 0845 elseif p.Results.tail<0, 0846 fprintf('Tail of tests: Lower\n'); 0847 else 0848 fprintf('Tail of tests: Two-Tailed\n'); 0849 end 0850 if n_chan==1, 0851 fprintf('Multiple comparison correction method: Not Applicable\n'); 0852 else 0853 if strcmpi(p.Results.correction,'perm'), 0854 fprintf('Multiple comparison correction method: Within-subjects permutation test\n'); 0855 elseif strcmpi(p.Results.correction,'FDR BH'), 0856 fprintf('Benjamini & Hochberg FDR control procedure\n'); 0857 elseif strcmpi(p.Results.correction,'FDR BY'), 0858 fprintf('Benjamini & Yekutieli FDR control procedure\n'); 0859 else 0860 fprintf('Multiple comparison correction method: Bonferroni-Holm procedure\n'); 0861 end 0862 end 0863 fprintf('Specified Family-Wise Alpha Level: %.3f\n',p.Results.alpha); 0864 if strcmpi(p.Results.correction,'perm'), 0865 fprintf('Estimated Perm Test Family-Wise Alpha Level: %.3f\n', ... 0866 regress_results.perm_test_info.est_alpha); 0867 fprintf('# of permutations used for permutation test: %d\n', ... 0868 regress_results.perm_test_info.n_perm); 0869 end 0870 fprintf('Number of channels analyzed: %d\n',n_chan); 0871 fprintf('Mean EEG amplitude measured from %d to %d ms.\n',time_wind(1), ... 0872 time_wind(2)); 0873 fprintf('Mean (SD) # of trials per participant: %.1f (%.1f)\n',mean(n_regress_trial),std(n_regress_trial)); 0874 fprintf('Minimum and maximum # of trials per participant: %.1f to %.1f\n',min(n_regress_trial),max(n_regress_trial)); 0875 fprintf('**RESULTS (Summarized across all channels)**\n'); 0876 if n_pred>1, 0877 fprintf('Predictor\t#_of_Sig_Tests\tMedian(IQR)_Collinearity\tp-values\n'); 0878 for a=1:n_pred+1, 0879 mn_p=min(min(regress_results.p_coefs_cor(:,:,a))); 0880 if mn_p>=p.Results.alpha, 0881 %no significant effects for this predictor 0882 if a==1, 0883 %intercept term 0884 fprintf('%s\t\t%d\t\t\tNA\t\t\tp>=%f\n',regress_results.predictor_names{a}, ... 0885 sum(sum(regress_results.h_coefs_cor(:,:,a))),mn_p); %intercept term has no collinearity measure 0886 else 0887 fprintf('%s\t\t%d\t\t\t%.2f(%.2f)\t\tp>=%f\n',regress_results.predictor_names{a}, ... 0888 sum(sum(regress_results.h_coefs_cor(:,:,a))),median(R_collinear(a,:)), ... 0889 iqr(R_collinear(a,:)),mn_p); 0890 end 0891 else 0892 %at least one significant effects for this predictor 0893 p_coefs_cor=regress_results.p_coefs_cor(:,:,a); 0894 sig_p_id=find(p_coefs_cor<p.Results.alpha); 0895 mx_p=max(p_coefs_cor(sig_p_id)); 0896 if a==1, 0897 %intercept term 0898 fprintf('%s\t\t%d\t\t\tNA\t\t%f>=sig p>=%f\n',regress_results.predictor_names{a}, ... 0899 sum(sum(regress_results.h_coefs_cor(:,:,a))),mx_p,mn_p); %intercept term has no collinearity measure 0900 else 0901 fprintf('%s\t\t%d\t\t\t%.2f(%.2f)\t%f>=sig p>=%f\n',regress_results.predictor_names{a}, ... 0902 sum(sum(regress_results.h_coefs_cor(:,:,a))),median(R_collinear(a,:)), ... 0903 iqr(R_collinear(a,:)),mx_p,mn_p); 0904 end 0905 end 0906 end 0907 else 0908 fprintf('Predictor\t#_of_Sig_Channels\t\tp-values\n'); 0909 for a=1:2, 0910 mn_p=min(min(regress_results.p_coefs_cor(:,:,a))); 0911 if mn_p>=p.Results.alpha, 0912 %no significant effects for this predictor 0913 fprintf('%s\t\t%d\t\t\tp>=%f\n',regress_results.predictor_names{a}, ... 0914 sum(sum(regress_results.h_coefs_cor(:,:,a))),mn_p); 0915 else 0916 %at least one significant effect for this predictor 0917 p_coefs_cor=regress_results.p_coefs_cor(:,:,a); 0918 sig_p_id=find(p_coefs_cor<p.Results.alpha); 0919 mx_p=max(p_coefs_cor(sig_p_id)); 0920 % sig_p_id=find(regress_results.p_coefs_cor(:,a)<p.Results.alpha); 0921 % mx_p=max(regress_results.p_coefs_cor(sig_p_id,a)); 0922 fprintf('%s\t\t%d\t\t%f>=sig p>=%f\n',regress_results.predictor_names{a}, ... 0923 sum(sum(regress_results.h_coefs_cor(:,:,a))),mx_p,mn_p); 0924 end 0925 end 0926 end 0927 0928 %% Print Results to Output File 0929 if ~isempty(p.Results.output_file), 0930 watchit('This function cannot produce an output text file yet.'); 0931 if 0 0932 [fid msg]=fopen(p.Results.output_file,'w'); 0933 if fid==-1, 0934 error('Cannot create file %s for writing. According to fopen.m: %s.', ... 0935 p.Results.file,msg); 0936 else 0937 % Header Information 0938 fprintf(fid,'**ANALYSIS_PARAMETERS**\n'); 0939 fprintf(fid,'Number of participants: %d\n',n_sub); 0940 fprintf(fid,'Regression analysis t-score degrees of freedom: %d\n',n_sub-1); 0941 if p.Results.tail>0, 0942 fprintf(fid,'Tail of tests: Upper\n'); 0943 elseif p.Results.tail<0, 0944 fprintf(fid,'Tail of tests: Lower\n'); 0945 else 0946 fprintf(fid,'Tail of tests: Two-Tailed\n'); 0947 end 0948 if n_chan==1, 0949 fprintf(fid,'Multiple comparison correction method: Not Applicable\n'); 0950 else 0951 if strcmpi(p.Results.correction,'perm'), 0952 fprintf(fid,'Multiple comparison correction method: Within-subjects permutation test\n'); 0953 elseif strcmpi(p.Results.correction,'FDR BH'), 0954 fprintf(fid,'Benjamini & Hochberg FDR control procedure\n'); 0955 elseif strcmpi(p.Results.correction,'FDR BY'), 0956 fprintf(fid,'Benjamini & Yekutieli FDR control procedure\n'); 0957 else 0958 fprintf(fid,'Multiple comparison correction method: Bonferroni-Holm procedure\n'); 0959 end 0960 end 0961 fprintf(fid,'Family-Wise Alpha Level: %.3f\n',p.Results.alpha); 0962 if strcmpi(p.Results.correction,'perm'), 0963 fprintf(fid,'Estimated Perm Test Family-Wise Alpha Level: %.3f\n', ... 0964 regress_results.perm_test_info.est_alpha); 0965 fprintf(fid,'# of permutations used for permutation test: %d\n', ... 0966 regress_results.perm_test_info.n_perm); 0967 end 0968 fprintf(fid,'Number of channels analyzed: %d\n',n_chan); 0969 fprintf(fid,'Mean EEG amplitude measured from %d to %d ms.\n',time_wind(1), ... 0970 time_wind(2)); 0971 fprintf(fid,'Mean (SD) # of trials per participant: %.1f (%.1f)\n',mean(n_regress_trial),std(n_regress_trial)); 0972 fprintf(fid,'Minimum and maximum # of trials per participant: %.1f to %.1f\n',min(n_regress_trial),max(n_regress_trial)); 0973 fprintf(fid,'*.set_files_used_in_analysis #_of_trials_from_setfile_used\n'); 0974 for a=1:n_sub, 0975 fprintf(fid,'%s %d\n',infiles{a},n_regress_trial(a)); 0976 end 0977 fprintf(fid,'**REGRESSION_RESULTS_SUMMARY**\n'); 0978 if n_pred>1, 0979 fprintf(fid,'Predictor\t#_of_Sig_Channels\tMedian(IQR)_Collinearity\tp-values\n'); 0980 for a=1:n_pred+1, 0981 mn_p=min(regress_results.p_coefs_cor(:,a)); 0982 if mn_p>=p.Results.alpha, 0983 %no significant effects for this predictor 0984 if a==1, 0985 %intercept term 0986 fprintf(fid,'%s\t\t%d\t\t\tNA\t\t\tp>=%f\n',regress_results.predictor_names{a}, ... 0987 sum(regress_results.h_coefs_cor(:,a)),mn_p); %intercept term has no collinearity measure 0988 else 0989 fprintf(fid,'%s\t\t%d\t\t\t%.2f(%.2f)\t\tp>=%f\n',regress_results.predictor_names{a}, ... 0990 sum(regress_results.h_coefs_cor(:,a)),median(R_collinear(a,:)), ... 0991 iqr(R_collinear(a,:)),mn_p); 0992 end 0993 else 0994 %at least one significant effects for this predictor 0995 sig_p_id=find(regress_results.p_coefs_cor(:,a)<p.Results.alpha); 0996 mx_p=max(regress_results.p_coefs_cor(sig_p_id,a)); 0997 if a==1, 0998 %intercept term 0999 fprintf(fid,'%s\t\t%d\t\t\tNA\t\t%f>=sig p>=%f\n',regress_results.predictor_names{a}, ... 1000 sum(regress_results.h_coefs_cor(:,a)),mx_p,mn_p); %intercept term has no collinearity measure 1001 else 1002 fprintf(fid,'%s\t\t%d\t\t\t%.2f(%.2f)\t%f>=sig p>=%f\n',regress_results.predictor_names{a}, ... 1003 sum(regress_results.h_coefs_cor(:,a)),median(R_collinear(a,:)), ... 1004 iqr(R_collinear(a,:)),mx_p,mn_p); 1005 end 1006 end 1007 end 1008 else 1009 fprintf(fid,'Predictor\t#_of_Sig_Channels\t\tp-values\n'); 1010 for a=1:2, 1011 mn_p=min(regress_results.p_coefs_cor(:,a)); 1012 if mn_p>=p.Results.alpha, 1013 %no significant effects for this predictor 1014 fprintf(fid,'%s\t\t%d\t\t\tp>=%f\n',regress_results.predictor_names{a}, ... 1015 sum(regress_results.h_coefs_cor(:,a)),mn_p); 1016 else 1017 %at least one significant effect for this predictor 1018 sig_p_id=find(regress_results.p_coefs_cor(:,a)<p.Results.alpha); 1019 mx_p=max(regress_results.p_coefs_cor(sig_p_id,a)); 1020 fprintf(fid,'%s\t\t%d\t\t%f>=sig p>=%f\n',regress_results.predictor_names{a}, ... 1021 sum(regress_results.h_coefs_cor(:,a)),mx_p,mn_p); 1022 end 1023 end 1024 end 1025 1026 fprintf(fid,'**REGRESSION_RESULTS_PER_CHANNEL**\n'); 1027 for c=1:n_chan, 1028 fprintf(fid,'Results_for_Channel: %s \n',regress_results.chanlocs(c).labels); 1029 fprintf(fid,'Predictor\tMean_Coefficient\t95%c_CI_lower_bound\t95%c_CI_upper_bound\tt-score\tp-value(Uncorrected)\tSignificant(Corrected)\tCohen''s_d\n', ... 1030 37,37); 1031 for pred=1:n_pred+1, 1032 %mean coefficient 1033 %95% CI 1034 %t-score 1035 %p-value 1036 %Cohen's d 1037 if pred==1, 1038 pred_name='Intercept'; 1039 else 1040 pred_name=predictors{pred-1}; 1041 end 1042 fprintf(fid,'%s \t%f \t%f \t%f \t%.3f \t%.3f \t%d \t%.2f\n', ... 1043 pred_name,regress_results.mn_coefs(c,pred), ... 1044 regress_results.low_ci(c,pred),regress_results.hi_ci(c,pred), ... 1045 regress_results.t_coefs(c,pred),regress_results.p_coefs_uncor(c,pred), ... 1046 regress_results.h_coefs_cor(c,pred),regress_results.d_coefs(c,pred)); 1047 end 1048 end 1049 fclose(fid); 1050 end 1051 end 1052 end 1053 1054 %% Plot Results 1055 if p.Results.plot_raster, 1056 for a=1:length(regress_results.predictor_names) 1057 sig_raster_coefs(regress_results,'predictor',a); 1058 end 1059 end 1060 %%% End of Main Function %%% 1061 1062 1063 function all_same=all_same_cell_elements(cell_array) 1064 %function to ensure that all elements a 1D cell array are the same 1065 1066 n=length(cell_array); 1067 all_same=1; 1068 for a=2:n, 1069 if ~(isnan(cell_array{1}) && isnan(cell_array{a})) 1070 if ~isequal(cell_array{1},cell_array{a}), 1071 all_same=0; 1072 break; 1073 end 1074 end 1075 end 1076 1077