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(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_topos',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 0346 %Check for artifact ICs and remove them if there 0347 %[EEG, ICs_removed]=remove_art_ics(EEG,p.Results.bsln_wind,VERBLEVEL);?? old code 0348 ICs_removed=remove_artifact_ics(p.Results.bsln_wind,VERBLEVEL); 0349 0350 %Filter the data if requested 0351 if ~isempty(flt) 0352 %first error check 0353 if ~isempty(flt), 0354 if max(flt)>(EEG.srate/2), 0355 error('''filt'' parameters need to be less than or equal to (sampling rate)/2 (i.e., %f).',EEG.srate/2); 0356 elseif (flt(2)==(EEG.srate/2)) && (flt(1)==0), 0357 error('If second element of ''filt'' parameter is EEG.srate/2, then the first element must be greater than 0.'); 0358 elseif abs(flt(2))<=abs(flt(1)), 0359 error('Second element of ''filt'' parameters must be greater than first in absolute value.'); 0360 elseif (flt(1)<0) || (flt(2)<0), 0361 if (flt(1)>=0) || (flt(2)>=0), 0362 error('BOTH parameters of ''filt'' need to be greater than or equal to zero OR need to be negative.'); 0363 end 0364 if min(flt)<=(-EEG.srate/2), 0365 error('''filt'' parameters need to be greater than sampling rate/2 (i.e., -%f) when creating a stop band.',EEG.srate/2); 0366 end 0367 end 0368 0369 end 0370 EEG=butterfilt(EEG,flt,[],p.Results.n_pad); 0371 end 0372 0373 %Baseline EEG (unless already baselined by remove_art_ics 0374 if ~isempty(p.Results.bsln_wind) && isempty(ICs_removed), 0375 EEG=baseline_EEG(EEG,p.Results.bsln_wind,VERBLEVEL); 0376 end 0377 0378 n_pos_chan=EEG.nbchan; 0379 chan_ids=zeros(1,n_pos_chan); 0380 0381 if sub==1, 0382 %% Identify predictors 0383 if isempty(p.Results.predictors), 0384 %Print list of possible predictors on command line 0385 epoch_flds=fieldnames(EEG.epoch); 0386 n_epoch_flds=length(epoch_flds); 0387 poss_pred=cell(1,n_epoch_flds); 0388 poss_pred_fullname=cell(1,n_epoch_flds); 0389 0390 fprintf('\n**** Choose Predictor(s) for Regression Analysis ****\n'); 0391 ct=0; 0392 for fld_loop=1:n_epoch_flds, 0393 if ~strcmpi('event',epoch_flds{fld_loop}) && ~strcmpi('eventtype',epoch_flds{fld_loop}), 0394 ct=ct+1; 0395 poss_pred_fullname{ct}=epoch_flds{fld_loop}; 0396 if strcmpi('event',epoch_flds{fld_loop}(1:5)), 0397 poss_pred{ct}=epoch_flds{fld_loop}(6:end); 0398 else 0399 poss_pred{ct}=epoch_flds{fld_loop}; 0400 end 0401 fprintf('%d) %s\n',ct,poss_pred{ct}); 0402 end 0403 end 0404 %Is reaction time a field? 0405 EEG_flds=fieldnames(EEG); 0406 if ismember('rtmsec',EEG_flds) && ismember('rtbins',EEG_flds), 0407 ct=ct+1; 0408 poss_pred{ct}='Reaction Time'; 0409 poss_pred_fullname{ct}='Reaction Time'; 0410 fprintf('%d) %s\n',ct,poss_pred{ct}); 0411 end 0412 0413 pred_ids=[]; 0414 while isempty(pred_ids), 0415 pred_ids=input('Enter the number(s) of desired predictor(s): ','s'); 0416 fprintf('\n'); 0417 pred_ids=str2num(pred_ids); 0418 if max(pred_ids)>ct, 0419 fprintf('ERROR: There are no predictor indices greater than %d.\n',ct); 0420 pred_ids=[]; 0421 elseif min(pred_ids)<1, 0422 fprintf('ERROR: All predictor indices must be greater than 0.\n'); 0423 pred_ids=[]; 0424 end 0425 end 0426 predictors=poss_pred_fullname(unique(pred_ids)); 0427 predictor_names_neat=poss_pred(unique(pred_ids)); 0428 if ismember('Reaction Time',predictors), 0429 rt_selected=1; 0430 else 0431 rt_selected=0; 0432 end 0433 else 0434 %get possible predictors 0435 epoch_flds=fieldnames(EEG.epoch); 0436 n_epoch_flds=length(epoch_flds); 0437 poss_pred=cell(1,n_epoch_flds); 0438 poss_pred_fullname=cell(1,n_epoch_flds); 0439 ct=0; 0440 for fld_loop=1:n_epoch_flds, 0441 if ~strcmpi('event',epoch_flds{fld_loop}) && ~strcmpi('eventtype',epoch_flds{fld_loop}), 0442 ct=ct+1; 0443 poss_pred_fullname{ct}=epoch_flds{fld_loop}; 0444 if strcmpi('event',epoch_flds{fld_loop}(1:5)), 0445 poss_pred{ct}=epoch_flds{fld_loop}(6:end); 0446 else 0447 poss_pred{ct}=epoch_flds{fld_loop}; 0448 end 0449 end 0450 end 0451 %Is reaction time a field? 0452 EEG_flds=fieldnames(EEG); 0453 if ismember('rtmsec',EEG_flds) && ismember('rtbins',EEG_flds), 0454 ct=ct+1; 0455 poss_pred{ct}='Reaction Time'; 0456 poss_pred_fullname{ct}='Reaction Time'; 0457 end 0458 0459 %figure out which of the possible predictors have been 0460 %requested 0461 if ischar(p.Results.predictors) 0462 requested_pred{1}=p.Results.predictors; 0463 else 0464 requested_pred=p.Results.predictors; 0465 end 0466 pred_ids=zeros(1,length(requested_pred)); 0467 ct=0; 0468 n_poss=length(poss_pred); 0469 for dloop=1:length(requested_pred), 0470 found=0; 0471 for eloop=1:n_poss, 0472 if strcmpi(requested_pred{dloop},poss_pred{eloop}) || ... 0473 strcmpi(requested_pred{dloop},poss_pred_fullname{eloop}) 0474 found=1; 0475 ct=ct+1; 0476 pred_ids(ct)=eloop; 0477 break; 0478 end 0479 end 0480 if ~found, 0481 fprintf('*** ERROR!! ***\n'); 0482 fprintf('Requested predictor %s was not found in this EEG variable.\n', ... 0483 requested_pred{dloop}); 0484 fprintf('Possible predictors are:\n'); 0485 for floop=1:n_poss, 0486 if ~isempty(poss_pred{floop}), 0487 fprintf('%d) %s\n',floop,poss_pred{floop}); 0488 end 0489 end 0490 error('Requested predictor not found.'); 0491 end 0492 end 0493 predictors=poss_pred_fullname(unique(pred_ids)); 0494 predictor_names_neat=poss_pred(unique(pred_ids)); 0495 if ismember('Reaction Time',predictors), 0496 rt_selected=1; 0497 else 0498 rt_selected=0; 0499 end 0500 end 0501 n_pred=length(predictors); 0502 0503 %% Figure out which channels to ignore if any 0504 % exclude and include chan options 0505 use_chan_labels=cell(1,n_pos_chan); %preallocate mem 0506 ct=0; 0507 if ~isempty(exclude_chans), 0508 for x=1:n_pos_chan, 0509 if ~ismember(EEG.chanlocs(x).labels,exclude_chans), 0510 ct=ct+1; 0511 chan_ids(ct)=x; 0512 use_chan_labels{ct}=EEG.chanlocs(x).labels; 0513 end 0514 end 0515 elseif ~isempty(include_chans), 0516 for x=1:n_pos_chan, 0517 if ismember(EEG.chanlocs(x).labels,include_chans), 0518 ct=ct+1; 0519 chan_ids(ct)=x; 0520 use_chan_labels{ct}=EEG.chanlocs(x).labels; 0521 end 0522 end 0523 else 0524 for x=1:n_pos_chan, 0525 use_chan_labels{x}=EEG.chanlocs(x).labels; 0526 end 0527 chan_ids=1:n_pos_chan; 0528 ct=n_pos_chan; 0529 end 0530 chan_ids=chan_ids(1:ct); 0531 use_chan_labels(ct+1:end)=[]; 0532 n_chan=length(chan_ids); 0533 chanlocs=EEG.chanlocs(chan_ids); 0534 0535 %Preallocate memory for regression coefficient estimates 0536 sub_coefs=zeros(n_chan,n_pred+1,n_sub); %add one for intercept 0537 R_collinear=zeros(n_pred+1,n_sub)*NaN; %add one for intercept 0538 else 0539 %Not the first subject, grab the same channels in the same order and throw error if missing used channel(s) 0540 ct=0; 0541 for a=1:n_chan, 0542 found=0; 0543 for b=1:n_pos_chan, 0544 %a=channel index for analysis 0545 %b=channel index for this particular subject 0546 if strcmpi(EEG.chanlocs(b).labels,use_chan_labels{a}), 0547 found=1; 0548 ct=ct+1; 0549 chan_ids(ct)=b; 0550 break; 0551 end 0552 end 0553 if ~found, 0554 error('Missing channel %s for participant %s.',use_chan_labels{a}, ... 0555 infiles{sub}); 0556 end 0557 end 0558 chan_ids=chan_ids(1:ct); 0559 end 0560 0561 poss_prdct=fieldnames(EEG.epoch); 0562 if rt_selected, 0563 EEG_flds=fieldnames(EEG); 0564 if ismember('rtmsec',EEG_flds) && ismember('rtbins',EEG_flds), 0565 poss_prdct{1+length(poss_prdct)}='Reaction Time'; 0566 end 0567 end 0568 missing_prdct=setdiff(predictors,poss_prdct); 0569 if ~isempty(missing_prdct), 0570 fprintf('ERROR: File %s is missing the following predictor(s): \n', ... 0571 infiles{sub}); 0572 for a=1:length(missing_prdct), 0573 fprintf('%s \n',missing_prdct{a}); 0574 end 0575 error('Regression aborted.'); 0576 end 0577 0578 n_trial=length(EEG.epoch); 0579 use_trial=zeros(1,n_trial); 0580 pred_val=zeros(n_pred,n_trial); 0581 EEG_val=zeros(n_chan,n_trial); 0582 for trial=1:n_trial, 0583 n_type=length(EEG.epoch(trial).eventtype); 0584 0585 if isempty(p.Results.use_bins), 0586 use_trial(trial)=1; 0587 else 0588 %Find out if current trial falls into desired bin 0589 for typ=1:n_type, 0590 for bin=p.Results.use_bins, 0591 if strcmpi(['bin' num2str(bin)],EEG.epoch(trial).eventtype{typ}), 0592 use_trial(trial)=1; 0593 break; %leave bin loop 0594 end 0595 end 0596 if use_trial(trial), 0597 break; %leave typ loop 0598 end 0599 end 0600 end 0601 0602 %Record predictor and response var values if trial is being used 0603 if use_trial(trial), 0604 for pred=1:n_pred, 0605 if strcmpi(predictors{pred},'Reaction Time'), 0606 n_rt_bin=length(EEG.rtbins{trial}); 0607 bin_rts=zeros(1,n_rt_bin); 0608 bin_ct=0; 0609 for b=1:n_rt_bin, 0610 if ismember(EEG.rtbins{trial}(b),p.Results.use_bins), 0611 bin_ct=bin_ct+1; 0612 bin_rts(bin_ct)=EEG.rtmsec{trial}(b); 0613 end 0614 end 0615 bin_rts=bin_rts(1:bin_ct); 0616 uni_rt=unique(bin_rts); 0617 if length(uni_rt)>1 0618 error('Trial %d has more than one reaction time for the bins you''ve selected.',trial); 0619 elseif isempty(uni_rt), 0620 watchit(sprintf('Epoch %d has an empty value for reaction time. Please use NaN for missing predictor values.\n', ... 0621 trial)); 0622 pred_val(pred,trial)=NaN; 0623 else 0624 pred_val(pred,trial)=uni_rt; 0625 end 0626 else 0627 cmnd=['emptyval=isempty(EEG.epoch(trial).' predictors{pred} ');']; 0628 eval(cmnd); %for some reason calling isempty.m wasn't working in the function so I used this workaround 0629 if emptyval, 0630 watchit(sprintf('Epoch %d has an empty value for predictor %s. Please use NaN for missing predictor values.\n', ... 0631 trial,predictor_names_neat{pred})); 0632 cmnd=['pred_val(pred,trial)=NaN;']; %?? Do RT too 0633 eval(cmnd); 0634 else 0635 cmnd=['temp_pred=EEG.epoch(trial).' predictors{pred} ';']; 0636 eval(cmnd); 0637 if iscell(temp_pred), 0638 if all_same_cell_elements(temp_pred), 0639 pred_val(pred,trial)=temp_pred{1}; 0640 else 0641 error('Epoch %d has multiple values for predictor %s.', ... 0642 trial,predictors{pred}); 0643 end 0644 else 0645 pred_val(pred,trial)=temp_pred; 0646 end 0647 end 0648 end 0649 end 0650 for c=1:n_chan, 0651 EEG_val(c,trial)=squeeze(mean(EEG.data(chan_ids(c),time_ids,trial),2)); 0652 end 0653 end 0654 end 0655 0656 %get rid of unused trials 0657 trial_ids=find(use_trial==1); 0658 if n_chan>1, 0659 EEG_val=EEG_val(:,trial_ids); 0660 else 0661 EEG_val=EEG_val(trial_ids); 0662 end 0663 if n_pred>1, 0664 pred_val=pred_val(:,trial_ids); 0665 else 0666 pred_val=pred_val(trial_ids); 0667 end 0668 0669 %get rid of NaN trials (if any) 0670 if n_pred>1, 0671 has_NaN=sum(isnan(pred_val)); 0672 non_NaN_ids=find(has_NaN==0); 0673 NaN_ids=find(has_NaN); 0674 pred_val=pred_val(:,non_NaN_ids); 0675 else 0676 non_NaN_ids=find(~isnan(pred_val)); 0677 NaN_ids=find(isnan(pred_val)); 0678 pred_val=pred_val(non_NaN_ids); 0679 end 0680 if n_chan>1, 0681 EEG_val=EEG_val(:,non_NaN_ids); 0682 else 0683 EEG_val=EEG_val(non_NaN_ids); 0684 end 0685 0686 n_nonNaN_trials=length(non_NaN_ids); 0687 n_selected_trials=length(trial_ids); 0688 if n_selected_trials~=n_nonNaN_trials, 0689 watchit(sprintf('%d of the selected trials have NaN predictor values for at least one predictor variable.', ... 0690 n_selected_trials-n_nonNaN_trials)); 0691 NaN_trial_ids=trial_ids(NaN_ids); 0692 if VERBLEVEL>2, 0693 disp(['The EEG.epoch indices of those trials are:' num2str(NaN_trial_ids)]); 0694 end 0695 end 0696 n_regress_trial(sub)=n_nonNaN_trials; 0697 fprintf('Number of usable trials: %d\n',n_nonNaN_trials); 0698 %Check if sufficient # of usable trials 0699 if ~n_nonNaN_trials, 0700 if isempty(p.Results.use_bins), 0701 error('Participant %s has no usable trials. Regression aborted.', ... 0702 infiles{sub}); 0703 else 0704 error('Participant %s has no trials for Bin(s) %s. Regression aborted.', ... 0705 infiles{sub},num2str(p.Results.use_bins)); 0706 end 0707 elseif n_nonNaN_trials<=n_pred, 0708 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).', ... 0709 infiles{sub},n_nonNaN_trials,n_pred); 0710 end 0711 0712 %Report min and max values for each predictor 0713 for pred_loop=1:n_pred, 0714 fprintf('%s (min/max): %f/%f\n',predictor_names_neat{pred_loop}, ... 0715 min(pred_val(pred_loop,:)),max(pred_val(pred_loop,:))); 0716 end 0717 0718 %perform multiple regression at each channel 0719 for c=1:n_chan, 0720 [coeff, S_err, XTXI, R_sq, R_sq_adj, F_val, Coef_stats, y_hat, residuals]= ... 0721 mregress(EEG_val(c,:)',pred_val',1); %final argument "1" means that an intercept term is included in the model 0722 %Note that the first element of coeff and associated variables 0723 %correpsponds to the intercept term 0724 sub_coefs(c,1:n_pred+1,sub)=coeff; 0725 end 0726 if sum(isnan(coeff)), 0727 error('The data for participant %s exhibits perfect collinearity and the regression coefficients cannot be determined.',infiles{sub}); 0728 end 0729 %Quantify degree of collinearity with other predictors 0730 if n_pred>1, 0731 for a=1:n_pred, 0732 [coeffCOL, S_errCOL, XTXIcol, R_sqCOL]=mregress(pred_val(a,:)', ... 0733 pred_val(setdiff(1:n_pred,a),:)',1); 0734 R_collinear(a+1,sub)=R_sqCOL; 0735 end 0736 end 0737 fprintf('\n'); 0738 end 0739 0740 regress_results.sub_coefs=sub_coefs; 0741 regress_results.mn_coefs=mean(sub_coefs,3); 0742 regress_results.se_coefs=stderr(sub_coefs,3); 0743 t_df=n_sub-1; 0744 crit_t=tinv(.975,t_df); %for 95% confidence intervals 0745 regress_results.hi_ci=regress_results.mn_coefs+crit_t*regress_results.se_coefs; 0746 regress_results.low_ci=regress_results.mn_coefs-crit_t*regress_results.se_coefs; 0747 regress_results.t_coefs=regress_results.mn_coefs./regress_results.se_coefs; 0748 regress_results.t_df=t_df; 0749 regress_results.d_coefs=regress_results.mn_coefs./std(sub_coefs,0,3); %0 indicates that n-1 SD formula will be used 0750 regress_results.F_coefs=zeros(n_chan,n_pred)*NaN; 0751 regress_results.tail=p.Results.tail; 0752 regress_results.p_coefs_uncor=zeros(n_chan,n_pred)*NaN; 0753 if n_chan==1, 0754 regress_results.crctn_method='Not applicable'; 0755 else 0756 regress_results.crctn_method=p.Results.correction; 0757 end 0758 regress_results.p_coefs_cor=regress_results.p_coefs_uncor; 0759 regress_results.h_coefs_cor=regress_results.p_coefs_uncor; 0760 regress_results.predictor_names=cell(1,n_pred+1); 0761 regress_results.predictor_names{1}='Intercept'; 0762 %regress_results.predictor_names(2:n_pred+1)=predictors; <=These are 0763 %predictor names with event in front of it 0764 regress_results.predictor_names(2:n_pred+1)=predictor_names_neat; 0765 regress_results.R_collinear=R_collinear; 0766 regress_results.family_alpha=p.Results.alpha; 0767 regress_results.chanlocs=chanlocs; 0768 if isempty(p.Results.use_bins), 0769 regress_results.bins='All'; 0770 else 0771 regress_results.bins=p.Results.use_bins; 0772 end 0773 regress_results.time_wind=time_wind; 0774 regress_results.participant_filenames=infiles; 0775 regress_results.perm_test_info.n_perm=NaN; 0776 regress_results.perm_test_info.tmx_ptile=NaN; 0777 regress_results.perm_test_info.seed_state=NaN; 0778 regress_results.perm_test_info.est_alpha=NaN; 0779 for pred=1:n_pred+1, 0780 fprintf('Computing significance of mean coefficient(s) for predictor %s.\n', ... 0781 regress_results.predictor_names{pred}); 0782 if ~isempty(p.Results.seed_state) && strcmpi(p.Results.correction,'perm'); 0783 fprintf('Seeding permutation test with specified seed_state.\n'); 0784 end 0785 for c=1:n_chan, 0786 if ~p.Results.tail 0787 %two-tailed test 0788 regress_results.p_coefs_uncor(c,pred)=2*(1-tcdf(abs(regress_results.t_coefs(c,pred)),n_sub-1)); 0789 elseif p.Results.tail==1, 0790 %upper-tailed test 0791 regress_results.p_coefs_uncor(c,pred)=(1-tcdf(regress_results.t_coefs(c,pred),n_sub-1)); 0792 else 0793 %lower-tailed test 0794 regress_results.p_coefs_uncor(c,pred)=tcdf(regress_results.t_coefs(c,pred),n_sub-1); 0795 end 0796 regress_results.F_coefs(c,pred)=finv(1-regress_results.p_coefs_uncor(c,pred),1,regress_results.t_df); 0797 end 0798 if n_chan==1, 0799 regress_results.p_coefs_cor(:,pred)=regress_results.p_coefs_uncor(:,pred); 0800 regress_results.h_coefs_cor(:,pred)=(regress_results.p_coefs_cor(:,pred)<p.Results.alpha); 0801 else 0802 if strcmpi(p.Results.correction,'Bonf-Holm'), 0803 [regress_results.p_coefs_cor(:,pred), regress_results.h_coefs_cor(:,pred)]= ... 0804 bonf_holm(regress_results.p_coefs_uncor(:,pred),p.Results.alpha); 0805 elseif strcmpi(p.Results.correction,'FDR BH'), 0806 [regress_results.h_coefs_cor(:,pred) dummy regress_results.p_coefs_cor(:,pred)]= ... 0807 fdr_bh(regress_results.p_coefs_uncor(:,pred),p.Results.alpha,'pdep','yes'); 0808 elseif strcmpi(p.Results.correction,'FDR BY'), 0809 [regress_results.h_coefs_cor(:,pred) dummy regress_results.p_coefs_cor(:,pred)]= ... 0810 fdr_bh(regress_results.p_coefs_uncor(:,pred),p.Results.alpha,'dep','yes'); 0811 else 0812 if pred>1 || isempty(p.Results.seed_state), 0813 [pval, t_orig, tmx_ptile, seed_state, est_alpha]=mxt_perm1(sub_coefs(:,pred,:), ... 0814 p.Results.n_perm,p.Results.alpha,p.Results.tail,VERBLEVEL); 0815 else 0816 %Just seed the first permutation test if seed_state 0817 %specified 0818 [pval, t_orig, tmx_ptile, seed_state, est_alpha]=mxt_perm1(sub_coefs(:,pred,:), ... 0819 p.Results.n_perm,p.Results.alpha,p.Results.tail,VERBLEVEL,p.Results.seed_state); 0820 end 0821 regress_results.p_coefs_cor(:,pred)=pval; 0822 regress_results.h_coefs_cor(:,pred)=pval<p.Results.alpha; 0823 regress_results.perm_test_info.n_perm=p.Results.n_perm; 0824 regress_results.perm_test_info.tmx_ptile=tmx_ptile; 0825 if pred==1, 0826 regress_results.perm_test_info.seed_state=seed_state; 0827 end 0828 regress_results.perm_test_info.est_alpha=est_alpha; 0829 end 0830 end 0831 end 0832 0833 %add error check if only one participant submitted? ?? 0834 fprintf('\n\n'); 0835 0836 %% Report Results to Command Line 0837 fprintf('**ANALYSIS PARAMETERS**\n'); 0838 fprintf('Number of participants: %d\n',n_sub); 0839 fprintf('Regression analysis t-score degrees of freedom: %d\n',t_df); 0840 if p.Results.tail>0, 0841 fprintf('Tail of tests: Upper\n'); 0842 elseif p.Results.tail<0, 0843 fprintf('Tail of tests: Lower\n'); 0844 else 0845 fprintf('Tail of tests: Two-Tailed\n'); 0846 end 0847 if n_chan==1, 0848 fprintf('Multiple comparison correction method: Not Applicable\n'); 0849 else 0850 if strcmpi(p.Results.correction,'perm'), 0851 fprintf('Multiple comparison correction method: Within-subjects permutation test\n'); 0852 elseif strcmpi(p.Results.correction,'FDR BH'), 0853 fprintf('Benjamini & Hochberg FDR control procedure\n'); 0854 elseif strcmpi(p.Results.correction,'FDR BY'), 0855 fprintf('Benjamini & Yekutieli FDR control procedure\n'); 0856 else 0857 fprintf('Multiple comparison correction method: Bonferroni-Holm procedure\n'); 0858 end 0859 end 0860 fprintf('Specified Family-Wise Alpha Level: %.3f\n',p.Results.alpha); 0861 if strcmpi(p.Results.correction,'perm'), 0862 fprintf('Estimated Perm Test Family-Wise Alpha Level: %.3f\n', ... 0863 regress_results.perm_test_info.est_alpha); 0864 fprintf('# of permutations used for permutation test: %d\n', ... 0865 regress_results.perm_test_info.n_perm); 0866 end 0867 fprintf('Number of channels analyzed: %d\n',n_chan); 0868 fprintf('Mean EEG amplitude measured from %d to %d ms.\n',time_wind(1), ... 0869 time_wind(2)); 0870 fprintf('Mean (SD) # of trials per participant: %.1f (%.1f)\n',mean(n_regress_trial),std(n_regress_trial)); 0871 fprintf('Minimum and maximum # of trials per participant: %.1f to %.1f\n',min(n_regress_trial),max(n_regress_trial)); 0872 fprintf('**RESULTS (Summarized across all channels)**\n'); 0873 if n_pred>1, 0874 fprintf('Predictor\t#_of_Sig_Channels\tMedian(IQR)_Collinearity\tp-values\n'); 0875 for a=1:n_pred+1, 0876 mn_p=min(regress_results.p_coefs_cor(:,a)); 0877 if mn_p>=p.Results.alpha, 0878 %no significant effects for this predictor 0879 if a==1, 0880 %intercept term 0881 fprintf('%s\t\t%d\t\t\tNA\t\t\tp>=%f\n',regress_results.predictor_names{a}, ... 0882 sum(regress_results.h_coefs_cor(:,a)),mn_p); %intercept term has no collinearity measure 0883 else 0884 fprintf('%s\t\t%d\t\t\t%.2f(%.2f)\t\tp>=%f\n',regress_results.predictor_names{a}, ... 0885 sum(regress_results.h_coefs_cor(:,a)),median(R_collinear(a,:)), ... 0886 iqr(R_collinear(a,:)),mn_p); 0887 end 0888 else 0889 %at least one significant effects for this predictor 0890 sig_p_id=find(regress_results.p_coefs_cor(:,a)<p.Results.alpha); 0891 mx_p=max(regress_results.p_coefs_cor(sig_p_id,a)); 0892 if a==1, 0893 %intercept term 0894 fprintf('%s\t\t%d\t\t\tNA\t\t%f>=sig p>=%f\n',regress_results.predictor_names{a}, ... 0895 sum(regress_results.h_coefs_cor(:,a)),mx_p,mn_p); %intercept term has no collinearity measure 0896 else 0897 fprintf('%s\t\t%d\t\t\t%.2f(%.2f)\t%f>=sig p>=%f\n',regress_results.predictor_names{a}, ... 0898 sum(regress_results.h_coefs_cor(:,a)),median(R_collinear(a,:)), ... 0899 iqr(R_collinear(a,:)),mx_p,mn_p); 0900 end 0901 end 0902 end 0903 else 0904 fprintf('Predictor\t#_of_Sig_Channels\t\tp-values\n'); 0905 for a=1:2, 0906 mn_p=min(regress_results.p_coefs_cor(:,a)); 0907 if mn_p>=p.Results.alpha, 0908 %no significant effects for this predictor 0909 fprintf('%s\t\t%d\t\t\tp>=%f\n',regress_results.predictor_names{a}, ... 0910 sum(regress_results.h_coefs_cor(:,a)),mn_p); 0911 else 0912 %at least one significant effect for this predictor 0913 sig_p_id=find(regress_results.p_coefs_cor(:,a)<p.Results.alpha); 0914 mx_p=max(regress_results.p_coefs_cor(sig_p_id,a)); 0915 fprintf('%s\t\t%d\t\t%f>=sig p>=%f\n',regress_results.predictor_names{a}, ... 0916 sum(regress_results.h_coefs_cor(:,a)),mx_p,mn_p); 0917 end 0918 end 0919 end 0920 0921 %% Print Results to Output File 0922 if ~isempty(p.Results.output_file), 0923 [fid msg]=fopen(p.Results.output_file,'w'); 0924 if fid==-1, 0925 error('Cannot create file %s for writing. According to fopen.m: %s.', ... 0926 p.Results.file,msg); 0927 else 0928 % Header Information 0929 fprintf(fid,'**ANALYSIS_PARAMETERS**\n'); 0930 fprintf(fid,'Number of participants: %d\n',n_sub); 0931 fprintf(fid,'Regression analysis t-score degrees of freedom: %d\n',n_sub-1); 0932 if p.Results.tail>0, 0933 fprintf(fid,'Tail of tests: Upper\n'); 0934 elseif p.Results.tail<0, 0935 fprintf(fid,'Tail of tests: Lower\n'); 0936 else 0937 fprintf(fid,'Tail of tests: Two-Tailed\n'); 0938 end 0939 if n_chan==1, 0940 fprintf(fid,'Multiple comparison correction method: Not Applicable\n'); 0941 else 0942 if strcmpi(p.Results.correction,'perm'), 0943 fprintf(fid,'Multiple comparison correction method: Within-subjects permutation test\n'); 0944 elseif strcmpi(p.Results.correction,'FDR BH'), 0945 fprintf(fid,'Benjamini & Hochberg FDR control procedure\n'); 0946 elseif strcmpi(p.Results.correction,'FDR BY'), 0947 fprintf(fid,'Benjamini & Yekutieli FDR control procedure\n'); 0948 else 0949 fprintf(fid,'Multiple comparison correction method: Bonferroni-Holm procedure\n'); 0950 end 0951 end 0952 fprintf(fid,'Family-Wise Alpha Level: %.3f\n',p.Results.alpha); 0953 if strcmpi(p.Results.correction,'perm'), 0954 fprintf(fid,'Estimated Perm Test Family-Wise Alpha Level: %.3f\n', ... 0955 regress_results.perm_test_info.est_alpha); 0956 fprintf(fid,'# of permutations used for permutation test: %d\n', ... 0957 regress_results.perm_test_info.n_perm); 0958 end 0959 fprintf(fid,'Number of channels analyzed: %d\n',n_chan); 0960 fprintf(fid,'Mean EEG amplitude measured from %d to %d ms.\n',time_wind(1), ... 0961 time_wind(2)); 0962 fprintf(fid,'Mean (SD) # of trials per participant: %.1f (%.1f)\n',mean(n_regress_trial),std(n_regress_trial)); 0963 fprintf(fid,'Minimum and maximum # of trials per participant: %.1f to %.1f\n',min(n_regress_trial),max(n_regress_trial)); 0964 fprintf(fid,'*.set_files_used_in_analysis #_of_trials_from_setfile_used\n'); 0965 for a=1:n_sub, 0966 fprintf(fid,'%s %d\n',infiles{a},n_regress_trial(a)); 0967 end 0968 fprintf(fid,'**REGRESSION_RESULTS_SUMMARY**\n'); 0969 if n_pred>1, 0970 fprintf(fid,'Predictor\t#_of_Sig_Channels\tMedian(IQR)_Collinearity\tp-values\n'); 0971 for a=1:n_pred+1, 0972 mn_p=min(regress_results.p_coefs_cor(:,a)); 0973 if mn_p>=p.Results.alpha, 0974 %no significant effects for this predictor 0975 if a==1, 0976 %intercept term 0977 fprintf(fid,'%s\t\t%d\t\t\tNA\t\t\tp>=%f\n',regress_results.predictor_names{a}, ... 0978 sum(regress_results.h_coefs_cor(:,a)),mn_p); %intercept term has no collinearity measure 0979 else 0980 fprintf(fid,'%s\t\t%d\t\t\t%.2f(%.2f)\t\tp>=%f\n',regress_results.predictor_names{a}, ... 0981 sum(regress_results.h_coefs_cor(:,a)),median(R_collinear(a,:)), ... 0982 iqr(R_collinear(a,:)),mn_p); 0983 end 0984 else 0985 %at least one significant effects for this predictor 0986 sig_p_id=find(regress_results.p_coefs_cor(:,a)<p.Results.alpha); 0987 mx_p=max(regress_results.p_coefs_cor(sig_p_id,a)); 0988 if a==1, 0989 %intercept term 0990 fprintf(fid,'%s\t\t%d\t\t\tNA\t\t%f>=sig p>=%f\n',regress_results.predictor_names{a}, ... 0991 sum(regress_results.h_coefs_cor(:,a)),mx_p,mn_p); %intercept term has no collinearity measure 0992 else 0993 fprintf(fid,'%s\t\t%d\t\t\t%.2f(%.2f)\t%f>=sig p>=%f\n',regress_results.predictor_names{a}, ... 0994 sum(regress_results.h_coefs_cor(:,a)),median(R_collinear(a,:)), ... 0995 iqr(R_collinear(a,:)),mx_p,mn_p); 0996 end 0997 end 0998 end 0999 else 1000 fprintf(fid,'Predictor\t#_of_Sig_Channels\t\tp-values\n'); 1001 for a=1:2, 1002 mn_p=min(regress_results.p_coefs_cor(:,a)); 1003 if mn_p>=p.Results.alpha, 1004 %no significant effects for this predictor 1005 fprintf(fid,'%s\t\t%d\t\t\tp>=%f\n',regress_results.predictor_names{a}, ... 1006 sum(regress_results.h_coefs_cor(:,a)),mn_p); 1007 else 1008 %at least one significant effect for this predictor 1009 sig_p_id=find(regress_results.p_coefs_cor(:,a)<p.Results.alpha); 1010 mx_p=max(regress_results.p_coefs_cor(sig_p_id,a)); 1011 fprintf(fid,'%s\t\t%d\t\t%f>=sig p>=%f\n',regress_results.predictor_names{a}, ... 1012 sum(regress_results.h_coefs_cor(:,a)),mx_p,mn_p); 1013 end 1014 end 1015 end 1016 1017 fprintf(fid,'**REGRESSION_RESULTS_PER_CHANNEL**\n'); 1018 for c=1:n_chan, 1019 fprintf(fid,'Results_for_Channel: %s \n',regress_results.chanlocs(c).labels); 1020 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', ... 1021 37,37); 1022 for pred=1:n_pred+1, 1023 %mean coefficient 1024 %95% CI 1025 %t-score 1026 %p-value 1027 %Cohen's d 1028 if pred==1, 1029 pred_name='Intercept'; 1030 else 1031 pred_name=predictors{pred-1}; 1032 end 1033 fprintf(fid,'%s \t%f \t%f \t%f \t%.3f \t%.3f \t%d \t%.2f\n', ... 1034 pred_name,regress_results.mn_coefs(c,pred), ... 1035 regress_results.low_ci(c,pred),regress_results.hi_ci(c,pred), ... 1036 regress_results.t_coefs(c,pred),regress_results.p_coefs_uncor(c,pred), ... 1037 regress_results.h_coefs_cor(c,pred),regress_results.d_coefs(c,pred)); 1038 end 1039 end 1040 fclose(fid); 1041 end 1042 end 1043 1044 %% Plot Results 1045 if p.Results.plot_topos && n_chan>1, 1046 sig_topo_coefs(regress_results); 1047 sig_topo_coefs(regress_results,'statistic','t'); 1048 end 1049 1050 %%% End of Main Function %%% 1051 1052 1053 function all_same=all_same_cell_elements(cell_array) 1054 %function to ensure that all elements a 1D cell array are the same 1055 1056 n=length(cell_array); 1057 all_same=1; 1058 for a=2:n, 1059 if ~(isnan(cell_array{1}) && isnan(cell_array{a})) 1060 if ~isequal(cell_array{1},cell_array{a}), 1061 all_same=0; 1062 break; 1063 end 1064 end 1065 end 1066 1067