stats_tfrGND() - Performs a paired-samples (i.e., repeated measures) t-
test on the spectrogram in two bins of a MATLABmk tfrGND
variable with correction for multiple comparisons. The
analysis is done via the fieldtrip function
ft_freqstatistics.m.
Usage:
>> tfrGND=stats_tfrGND(tfrGND,binA,binB,varargin);
Required Input:
tfrGND - A MATLABmk tfrGND variable. These are produced by the
MATLABmk function tfrs2tfrGND.m. It contains spectrograms from
one or more participants, their grand averages, and
associated information.
binA - One of the bins in the tfrGND variable.
binB - One of the bins in the tfrGND variable. The analysis will be
done by subtracting binB from binA.
Optional Inputs:
'freq_limits' - [min max] Frequency limits (in Hz) of the frequency
band you wish to analyze. {default: lowest and highest
frequencies in the tfrGND variable}
'time_limits' - [min max] Time limits (in ms) of the window you wish
to analyze. {default: lowest and highest time
points in the tfrGND variable}
'tail' - [1,0, or -1] An integer specifying the tail of the
hypothesis test. "1" means upper-tailed (i.e., alternative
hypothesis is that the spectrogram difference is greater
than the mean of the null hypothesis). "0" means two-
tailed (i.e., alternative hypothesis is that the
spectrogram difference is not equal to the mean of the null
hypothesis). "-1" means lower-tailed (i.e., alternative
hypothesis is that the spectrogram difference is less
than the mean of the null hypothesis). {default: 0}
'alpha' - A number between 0 and 1 specifying the family-wise
alpha level of the test. {default: 0.05}
'correctm' - [string] The method used to correct for multiple
comparisons. Options are:
'max' = tmax permutation based method that
strongly controls the family-wise
error rate (FWER).
'cluster' = Cluster-based permutation test that
weakly controls FWER. See Maris &
Oostenveld (2007) for details.
{default}
'fdr_bh' = Benjamini & Hochberg (1995) method for
control of the false discovery rate.
'fdr_by' = Benjamini & Yekutieli (1999) method for
control of the false discovery rate.
'none' = No correction for multiple comparisons.
'n_perm' - The number of permutations to use in the test (if a
permutation test is used). As this value increases,
the test takes longer to compute and the results become
more reliable. Manly (1997) suggests using at least
1000 permutations for an alpha level of 0.05 and at
least 5000 permutations for an alpha level of 0.01.
{default: 1000}
'clusteralpha' - A number between 0 and 1 specifying the test-wise
alpha level threshold a variable needs to pass in order
to be included in a cluster. Test-wise alpha levels
are computed without correction for multiple comparisons.
The smaller the value of 'cluteralpha' the more focal the
clusters found by cluster-based test tend to be.
'clusteralpha' should probably not have a value greater
than 5% since it could lead to declaring differences
as significant that would not have been declared
significant if no correction for multiple comparisons
had been performed. This parameter only has an effect if
a cluster-based test is performed and does NOT affect
the family-wise error rate. {default: 0.05}
'minnbchan' - The minimum number of above threshold neighbors an
above threshold dependent variable has to have in order
to be included in a cluseter. Having a non-zero value
makes clusters more focal by preventing "bridges"
between clusters. A non-zero value will further bias the
test against finding narrowly distributed clusters.
This only affects cluster-based tests. {default: 2}
'neighbourdist' - All sensors within this distance are considered
neighbors for the purposes of the cluster-based test.
To get the approximate value of this parameter in cm,
mulitply it by 8.9. {default: 0.61 which correspons
approximately to 5.4 cm}
'clusterstatistic' - The statistic used for cluster based multiple
comparison correction. This is the statistic that
is extracted from each random permutation of the
data. Options are:
'maxsum' = The most extreme sum of the t-scores
across all clusters. {default}
'maxsize' = The size of the largest cluster.
'max' = The most extreme t-score across all
the clusters.
'units' - ['dB' or 'raw'] If 'dB' spectrogram color scale will
be in units of decibels 10*log10(uV^2/Hz). If 'raw'
spectrogram color scale will be in units of uV^2/Hz.
{default: whatever the units were when the tfrGND
variable was first created}
'bsln_type' - The type of baseline used to normalize each
participant's data. Default is to use whatever baseline
type was used when the tfrGND variable was first
created. Options are:
'relative' - Power at each time point is divided by
the mean power in the baseline window.
'absolute' - The mean power in the baseline window
is subtracted from the power at each
time point.
'bsln_wind' - [min max] Time window (in ms) used to baseline the
spectrograms {default: whatever baseline window was
originally used with the tfrGND variable}
'exclude_chans' - A cell array of channel labels (e.g., {'A2','lle'})
or a string channel label (e.g., 'A2') to exclude from the
analysis . This option sacrifices spatial resolution
to increase test power by reducing the number of
comparisons. Use headinfo_spec.m to see the channel
labels stored in the tfrGND variable. You cannot
use both this option and 'include_chans' (below).{default:
not used, all channels included in test}
'include_chans' - A cell array of channel labels (e.g., {'RLOc','LLOc'})
or a string channel label (e.g., 'RLOc') to use in
the analysis. All other channels will be ignored.
This option sacrifices spatial resolution to increase
test power by reducing the number of comparisons.
Use headinfo_spec.m to see the channel labels stored in the
tfrGND variable. You cannot use both this option and
'exclude_chans' (above). {default: not used, all channels
included in test}
'avgoverchan' - ['yes' or 'no'] Perform the analysis on the mean of
all included channels. Doing this sacrifices spatial
resolution to increase test power by reducing the number of
comparisons. {default: 'no'}
'avgovertime' - ['yes' or 'no'] Perform the analysis on the mean of
all included time windows. Doing this sacrifices temporal
resolution to increase test power by reducing the number of
comparisons. {default: 'no'}
'avgoverfreq' - ['yes' or 'no'] Perform the analysis on the mean of
all included frequencies. Doing this sacrifices temporal
resolution to increase test power by reducing the number of
comparisons. {default: 'no'}
Output:
tfrGND - MATLABmk tfrGND variable with fieldtrip formatted test
results stored in tfrGND.stats. Explanations of the various
fields of tfrGND.stats are below. The fields vary somewhat
depending on test parameters (e.g., two-tailed vs. one-tailed,
permutation test vs. FDR control):
prob - the p-value of each variable after being
adjusted for multiple comparisons
posclusters - the p-value and cluster statistic of each
positive cluster. Note, the p-value is
always the p-value for an UPPER TAILED test
regardless of the actual tail of
the test.
posclusterslabelmat - a matrix that specifies which positive cluster
each variable belongs to (or 0 if the variable
does not belong to a positive cluster).
posdistribution - the maximum cluster statistic for
each random permutation of the data
negclusters - the p-value and cluster statistic of each
negative cluster. Note, the p-value is
always the p-value for a LOWER TAILED test
regardless of the actual tail of
the test.
negclusterslabelmat - a matrix that specifies which negative cluster
each variable belongs to (or 0 if the variable
does not belong to a negative cluster).
negdistribution: - the minimum cluster statistic for
each random permutation of the data
mask - a matrix with 1 for each variable with a
significant mult comp corrected p-value,
0 otherwise.
stat - the t-score of each variable
df - the degrees of freedom for t-scores
critval - the critical t-score(s) for significance
ref - the mean t-scores from all the
random permutations
dimord - What each dimension of matrices like 'prob,'
'mask,' and 'stat' correspondes to (e.g.,
which dimension corresponds to time, frequency,
and channel)
label - Electrode names
freq - Specific frequencies analyzed
time - Center of moving time windows (in
SECONDS)
cfg - The parameters used to call the
fieldtrip function ft_freqstatistics.
correctm - Which method was used for multiple
comparison correction.
bins - The bins that were analyzed. A minus
sign will indicate which bin was
subtracted from which bin.
Notes:
-If you run out of memory trying to use this function, you can reduce
memory load by:
1) Using fewer permutations (if you're performing a permutation test)
2) Analyzing a smaller frequency band
3) Analyzing a shorter time window
4) Computing spectrograms with fewer frequencies and using fewer time
windows
-Unlike a parametric test (e.g., an ANOVA), a discrete set of p-values
are possible (at most the number of possible permutations). Since the
number of possible permutations grows rapdily with the number of
participants, this is only issue for small sample sizes (e.g., 6
participants). When you have such a small sample size, the
limited number of p-values may make the test less conservative (e.g.,
you might be forced to use an alpha level of .0286 since it is the biggest
possible alpha level less than .05).
-When doing two-tailed test with cluster-based permutation test it is
possible to get multiple comparison corrected p-values as big as 2. When
using FDR control, multiple comparisons corrected p-values (often called
q-values for FDR control) can also be bigger than 1.
References:
Benjamini, Y. & Hochberg, Y. (1995) Controlling the false discovery
rate: A practical and powerful approach to multiple testing. Journal
of the Royal Statistical Society, Series B (Methodological). 57(1),
289-300.
Benjamini, Y. & Yekutieli, D. (2001) The control of the false discovery
rate in multiple testing under dependency. The Annals of Statistics.
29(4), 1165-1188.
Manly, B.F.J. (1997) Randomization, Bootstrap, and Monte Carlo Methods in
Biology. 2nd ed. Chapman and Hall, London.
Maris, E., & Oostenveld, R. (2007). Nonparametric statistical testing of
EEG- and MEG-data. Journal of Neuroscience Methods, 164(1), 177-190.
Author:
David Groppe
Kutaslab, 1/20110001 % stats_tfrGND() - Performs a paired-samples (i.e., repeated measures) t- 0002 % test on the spectrogram in two bins of a MATLABmk tfrGND 0003 % variable with correction for multiple comparisons. The 0004 % analysis is done via the fieldtrip function 0005 % ft_freqstatistics.m. 0006 % 0007 % Usage: 0008 % >> tfrGND=stats_tfrGND(tfrGND,binA,binB,varargin); 0009 % 0010 % Required Input: 0011 % tfrGND - A MATLABmk tfrGND variable. These are produced by the 0012 % MATLABmk function tfrs2tfrGND.m. It contains spectrograms from 0013 % one or more participants, their grand averages, and 0014 % associated information. 0015 % binA - One of the bins in the tfrGND variable. 0016 % binB - One of the bins in the tfrGND variable. The analysis will be 0017 % done by subtracting binB from binA. 0018 % 0019 % Optional Inputs: 0020 % 'freq_limits' - [min max] Frequency limits (in Hz) of the frequency 0021 % band you wish to analyze. {default: lowest and highest 0022 % frequencies in the tfrGND variable} 0023 % 'time_limits' - [min max] Time limits (in ms) of the window you wish 0024 % to analyze. {default: lowest and highest time 0025 % points in the tfrGND variable} 0026 % 'tail' - [1,0, or -1] An integer specifying the tail of the 0027 % hypothesis test. "1" means upper-tailed (i.e., alternative 0028 % hypothesis is that the spectrogram difference is greater 0029 % than the mean of the null hypothesis). "0" means two- 0030 % tailed (i.e., alternative hypothesis is that the 0031 % spectrogram difference is not equal to the mean of the null 0032 % hypothesis). "-1" means lower-tailed (i.e., alternative 0033 % hypothesis is that the spectrogram difference is less 0034 % than the mean of the null hypothesis). {default: 0} 0035 % 'alpha' - A number between 0 and 1 specifying the family-wise 0036 % alpha level of the test. {default: 0.05} 0037 % 'correctm' - [string] The method used to correct for multiple 0038 % comparisons. Options are: 0039 % 'max' = tmax permutation based method that 0040 % strongly controls the family-wise 0041 % error rate (FWER). 0042 % 'cluster' = Cluster-based permutation test that 0043 % weakly controls FWER. See Maris & 0044 % Oostenveld (2007) for details. 0045 % {default} 0046 % 'fdr_bh' = Benjamini & Hochberg (1995) method for 0047 % control of the false discovery rate. 0048 % 'fdr_by' = Benjamini & Yekutieli (1999) method for 0049 % control of the false discovery rate. 0050 % 'none' = No correction for multiple comparisons. 0051 % 'n_perm' - The number of permutations to use in the test (if a 0052 % permutation test is used). As this value increases, 0053 % the test takes longer to compute and the results become 0054 % more reliable. Manly (1997) suggests using at least 0055 % 1000 permutations for an alpha level of 0.05 and at 0056 % least 5000 permutations for an alpha level of 0.01. 0057 % {default: 1000} 0058 % 'clusteralpha' - A number between 0 and 1 specifying the test-wise 0059 % alpha level threshold a variable needs to pass in order 0060 % to be included in a cluster. Test-wise alpha levels 0061 % are computed without correction for multiple comparisons. 0062 % The smaller the value of 'cluteralpha' the more focal the 0063 % clusters found by cluster-based test tend to be. 0064 % 'clusteralpha' should probably not have a value greater 0065 % than 5% since it could lead to declaring differences 0066 % as significant that would not have been declared 0067 % significant if no correction for multiple comparisons 0068 % had been performed. This parameter only has an effect if 0069 % a cluster-based test is performed and does NOT affect 0070 % the family-wise error rate. {default: 0.05} 0071 % 'minnbchan' - The minimum number of above threshold neighbors an 0072 % above threshold dependent variable has to have in order 0073 % to be included in a cluseter. Having a non-zero value 0074 % makes clusters more focal by preventing "bridges" 0075 % between clusters. A non-zero value will further bias the 0076 % test against finding narrowly distributed clusters. 0077 % This only affects cluster-based tests. {default: 2} 0078 % 'neighbourdist' - All sensors within this distance are considered 0079 % neighbors for the purposes of the cluster-based test. 0080 % To get the approximate value of this parameter in cm, 0081 % mulitply it by 8.9. {default: 0.61 which correspons 0082 % approximately to 5.4 cm} 0083 % 'clusterstatistic' - The statistic used for cluster based multiple 0084 % comparison correction. This is the statistic that 0085 % is extracted from each random permutation of the 0086 % data. Options are: 0087 % 'maxsum' = The most extreme sum of the t-scores 0088 % across all clusters. {default} 0089 % 'maxsize' = The size of the largest cluster. 0090 % 'max' = The most extreme t-score across all 0091 % the clusters. 0092 % 'units' - ['dB' or 'raw'] If 'dB' spectrogram color scale will 0093 % be in units of decibels 10*log10(uV^2/Hz). If 'raw' 0094 % spectrogram color scale will be in units of uV^2/Hz. 0095 % {default: whatever the units were when the tfrGND 0096 % variable was first created} 0097 % 'bsln_type' - The type of baseline used to normalize each 0098 % participant's data. Default is to use whatever baseline 0099 % type was used when the tfrGND variable was first 0100 % created. Options are: 0101 % 'relative' - Power at each time point is divided by 0102 % the mean power in the baseline window. 0103 % 'absolute' - The mean power in the baseline window 0104 % is subtracted from the power at each 0105 % time point. 0106 % 'bsln_wind' - [min max] Time window (in ms) used to baseline the 0107 % spectrograms {default: whatever baseline window was 0108 % originally used with the tfrGND variable} 0109 % 'exclude_chans' - A cell array of channel labels (e.g., {'A2','lle'}) 0110 % or a string channel label (e.g., 'A2') to exclude from the 0111 % analysis . This option sacrifices spatial resolution 0112 % to increase test power by reducing the number of 0113 % comparisons. Use headinfo_spec.m to see the channel 0114 % labels stored in the tfrGND variable. You cannot 0115 % use both this option and 'include_chans' (below).{default: 0116 % not used, all channels included in test} 0117 % 'include_chans' - A cell array of channel labels (e.g., {'RLOc','LLOc'}) 0118 % or a string channel label (e.g., 'RLOc') to use in 0119 % the analysis. All other channels will be ignored. 0120 % This option sacrifices spatial resolution to increase 0121 % test power by reducing the number of comparisons. 0122 % Use headinfo_spec.m to see the channel labels stored in the 0123 % tfrGND variable. You cannot use both this option and 0124 % 'exclude_chans' (above). {default: not used, all channels 0125 % included in test} 0126 % 'avgoverchan' - ['yes' or 'no'] Perform the analysis on the mean of 0127 % all included channels. Doing this sacrifices spatial 0128 % resolution to increase test power by reducing the number of 0129 % comparisons. {default: 'no'} 0130 % 'avgovertime' - ['yes' or 'no'] Perform the analysis on the mean of 0131 % all included time windows. Doing this sacrifices temporal 0132 % resolution to increase test power by reducing the number of 0133 % comparisons. {default: 'no'} 0134 % 'avgoverfreq' - ['yes' or 'no'] Perform the analysis on the mean of 0135 % all included frequencies. Doing this sacrifices temporal 0136 % resolution to increase test power by reducing the number of 0137 % comparisons. {default: 'no'} 0138 % 0139 % Output: 0140 % tfrGND - MATLABmk tfrGND variable with fieldtrip formatted test 0141 % results stored in tfrGND.stats. Explanations of the various 0142 % fields of tfrGND.stats are below. The fields vary somewhat 0143 % depending on test parameters (e.g., two-tailed vs. one-tailed, 0144 % permutation test vs. FDR control): 0145 % prob - the p-value of each variable after being 0146 % adjusted for multiple comparisons 0147 % posclusters - the p-value and cluster statistic of each 0148 % positive cluster. Note, the p-value is 0149 % always the p-value for an UPPER TAILED test 0150 % regardless of the actual tail of 0151 % the test. 0152 % posclusterslabelmat - a matrix that specifies which positive cluster 0153 % each variable belongs to (or 0 if the variable 0154 % does not belong to a positive cluster). 0155 % posdistribution - the maximum cluster statistic for 0156 % each random permutation of the data 0157 % negclusters - the p-value and cluster statistic of each 0158 % negative cluster. Note, the p-value is 0159 % always the p-value for a LOWER TAILED test 0160 % regardless of the actual tail of 0161 % the test. 0162 % negclusterslabelmat - a matrix that specifies which negative cluster 0163 % each variable belongs to (or 0 if the variable 0164 % does not belong to a negative cluster). 0165 % negdistribution: - the minimum cluster statistic for 0166 % each random permutation of the data 0167 % mask - a matrix with 1 for each variable with a 0168 % significant mult comp corrected p-value, 0169 % 0 otherwise. 0170 % stat - the t-score of each variable 0171 % df - the degrees of freedom for t-scores 0172 % critval - the critical t-score(s) for significance 0173 % ref - the mean t-scores from all the 0174 % random permutations 0175 % dimord - What each dimension of matrices like 'prob,' 0176 % 'mask,' and 'stat' correspondes to (e.g., 0177 % which dimension corresponds to time, frequency, 0178 % and channel) 0179 % label - Electrode names 0180 % freq - Specific frequencies analyzed 0181 % time - Center of moving time windows (in 0182 % SECONDS) 0183 % cfg - The parameters used to call the 0184 % fieldtrip function ft_freqstatistics. 0185 % correctm - Which method was used for multiple 0186 % comparison correction. 0187 % bins - The bins that were analyzed. A minus 0188 % sign will indicate which bin was 0189 % subtracted from which bin. 0190 % 0191 % Notes: 0192 % -If you run out of memory trying to use this function, you can reduce 0193 % memory load by: 0194 % 1) Using fewer permutations (if you're performing a permutation test) 0195 % 2) Analyzing a smaller frequency band 0196 % 3) Analyzing a shorter time window 0197 % 4) Computing spectrograms with fewer frequencies and using fewer time 0198 % windows 0199 % 0200 % -Unlike a parametric test (e.g., an ANOVA), a discrete set of p-values 0201 % are possible (at most the number of possible permutations). Since the 0202 % number of possible permutations grows rapdily with the number of 0203 % participants, this is only issue for small sample sizes (e.g., 6 0204 % participants). When you have such a small sample size, the 0205 % limited number of p-values may make the test less conservative (e.g., 0206 % you might be forced to use an alpha level of .0286 since it is the biggest 0207 % possible alpha level less than .05). 0208 % 0209 % -When doing two-tailed test with cluster-based permutation test it is 0210 % possible to get multiple comparison corrected p-values as big as 2. When 0211 % using FDR control, multiple comparisons corrected p-values (often called 0212 % q-values for FDR control) can also be bigger than 1. 0213 % 0214 % 0215 % 0216 % References: 0217 % Benjamini, Y. & Hochberg, Y. (1995) Controlling the false discovery 0218 % rate: A practical and powerful approach to multiple testing. Journal 0219 % of the Royal Statistical Society, Series B (Methodological). 57(1), 0220 % 289-300. 0221 % 0222 % Benjamini, Y. & Yekutieli, D. (2001) The control of the false discovery 0223 % rate in multiple testing under dependency. The Annals of Statistics. 0224 % 29(4), 1165-1188. 0225 % 0226 % Manly, B.F.J. (1997) Randomization, Bootstrap, and Monte Carlo Methods in 0227 % Biology. 2nd ed. Chapman and Hall, London. 0228 % 0229 % Maris, E., & Oostenveld, R. (2007). Nonparametric statistical testing of 0230 % EEG- and MEG-data. Journal of Neuroscience Methods, 164(1), 177-190. 0231 % 0232 % 0233 % Author: 0234 % David Groppe 0235 % Kutaslab, 1/2011 0236 0237 %%%%%%%%%%%%%%%% Revision History %%%%%%%%%%%%%%%%% 0238 % 0239 0240 %% Other Notes 0241 %-cfg.clusteralpha is the threshold for being included in a cluster. 0242 %Fieldtrip's code automatically figures out if a 1 or 2 sided test 0243 %should be computed when deriving the test-wise p-value for each test 0244 % 0245 %-Another possible value for the optional argument 'clusterstatistic' is 0246 %'wcm'. I have no idea what it means. It's not commented on the fieldtrip 0247 %website nor is it explained in the code. The function in which it is 0248 %actually used is clusterstat.m. 0249 % 0250 %-If you set cfg.method='montecarlo' you get test-wise p-values via a 0251 %permutation test and then (optionally) correct for multiple comparisons 0252 %with any of the following methods: 'no', 'max', cluster', 'bonferoni', 0253 %'holms', 'fdr' with the cfg parameter cfg.correctm. Right now, stats_tfrGND.m 0254 %only supports 'cluster' and 'max'. For doing FDR correction I think it's 0255 %better to use an analytic test because the repeated-measures t-test is 0256 %very robust to violations of normality and will perform better when sample 0257 %sizes (or the number of permutations) are small. 0258 % 0259 %-The tmax procedure is implemented using Fieldtrip's built in function, 0260 %which I don't think is as accurate as our MATLABmk code. It should only 0261 %be a problem for small sample sizes or numbers of permutations. 0262 % 0263 %-Fieldtrip stats flow of control is: 0264 % 1) Call ft_freqstatistics 0265 % 2) Go to statistics_wrapper.m 0266 % 3) Go to statistics_montecarlo (for permutation based methods) 0267 % 4) Go to clusterstat.m for cluster-based tests 0268 % 0269 %-Note, at some point between January and March of 2011, the p-values in 0270 % tfrGND.stats{neo_test}.posclusters(a).prob were changed. Originally, the 0271 % p-values were for one-tailed tests (regardless of the actual tail of the 0272 % test). Now, they accurately reflect the tail of the test. The change 0273 % was made by the FieldTrip people. 0274 % 0275 0276 0277 %%% Future Work: 0278 % 0279 0280 function tfrGND=stats_tfrGND(tfrGND,binA,binB,varargin) 0281 0282 p=inputParser; 0283 p.addRequired('tfrGND',@isstruct); 0284 p.addRequired('binA',@(x) isnumeric(x) && length(x)==1); 0285 p.addRequired('binB',@(x) isnumeric(x) && length(x)==1); 0286 p.addParamValue('freq_limits',[],@(x) isnumeric(x) && length(x)==2); 0287 p.addParamValue('time_limits',[],@(x) isnumeric(x) && length(x)==2); 0288 p.addParamValue('correctm','cluster',@ischar); 0289 p.addParamValue('clusterstatistic','maxsum',@(x) (strcmpi(x,'maxsum') || strcmpi(x,'maxsize') || strcmpi(x,'max'))); 0290 p.addParamValue('alpha',.05,@(x) isnumeric(x) && (length(x)==1) && (x<=1) && (x>0)); 0291 p.addParamValue('clusteralpha',.05,@(x) isnumeric(x) && (length(x)==1) && (x<=1) && (x>0)); 0292 p.addParamValue('tail',0,@(x) isnumeric(x) && length(x)==1); 0293 p.addParamValue('neighbourdist',.61,@(x) isnumeric(x) && length(x)==1); 0294 p.addParamValue('minnbchan',2,@(x) isnumeric(x) && length(x)==1 && (x>=0)); 0295 p.addParamValue('n_perm',1000,@(x) isnumeric(x) && length(x)==1 && (x>0)); 0296 p.addParamValue('units',[],@(x) strcmpi(x,'dB') || strcmpi(x,'raw')); 0297 p.addParamValue('bsln_wind',[],@(x) isnumeric(x) && length(x)==2); 0298 p.addParamValue('bsln_type',[],@(x) strcmpi(x,'relative') || strcmpi(x,'absolute')); 0299 p.addParamValue('exclude_chans',[],@(x) ischar(x) || iscell(x)); 0300 p.addParamValue('include_chans',[],@(x) ischar(x) || iscell(x)); 0301 p.addParamValue('avgoverchan','no',@(x) strcmpi(x,'yes') || strcmpi(x,'no')); 0302 p.addParamValue('avgovertime','no',@(x) strcmpi(x,'yes') || strcmpi(x,'no')); 0303 p.addParamValue('avgoverfreq','no',@(x) strcmpi(x,'yes') || strcmpi(x,'no')); 0304 0305 p.parse(tfrGND,binA,binB,varargin{:}); 0306 0307 %%%% I. Transform spectrograms to dB or re-baseline (if requested) 0308 if isempty(p.Results.bsln_wind), 0309 bsln_wind=tfrGND.bsln_wind; %note p.Results.bsln_wind should be in units of ms 0310 else 0311 bsln_wind=p.Results.bsln_wind; 0312 end 0313 if isempty(p.Results.bsln_type), 0314 bsln_type=tfrGND.bsln_type; 0315 else 0316 bsln_type=p.Results.bsln_type; 0317 end 0318 if isempty(p.Results.units), 0319 units=tfrGND.units; 0320 else 0321 units=p.Results.units; 0322 end 0323 if ~strcmpi(units,tfrGND.units) || ... 0324 ~isequal(bsln_wind,tfrGND.bsln_wind) || ... 0325 ~strcmpi(bsln_type,tfrGND.bsln_type) 0326 tfrGND=baseline_tfrGND(tfrGND,units,bsln_wind,bsln_type); 0327 end 0328 0329 0330 %%%% II. Generic testing paramters (common to all multiple comparison correction methods) 0331 cfg = []; 0332 0333 %% Channel selection 0334 if ~isempty(p.Results.include_chans) 0335 if ~isempty(p.Results.exclude_chans), 0336 error('You cannot use BOTH include_chans and exclude_chans options.'); 0337 end 0338 if ischar(p.Results.include_chans) 0339 cfg.channel{1}=p.Results.include_chans; %cfg.channel needs to be a cell array I believe 0340 else 0341 cfg.channel=p.Results.include_chans; 0342 end 0343 elseif ~isempty(p.Results.exclude_chans), 0344 if ischar(p.Results.exclude_chans) 0345 exclude{1}=p.Results.exclude_chans; %convert to cell array 0346 else 0347 exclude=p.Results.exclude_chans; 0348 end 0349 n_exclude=length(exclude); 0350 cfg.channel=cell(1,n_exclude+1); 0351 cfg.channel{1}='all'; 0352 for a=1:n_exclude, 0353 cfg.channel{a+1}=['-' exclude{a}]; 0354 end 0355 else 0356 %default 0357 cfg.channel = {'all'}; 0358 end 0359 0360 %% Other paramaters 0361 cfg.tail = p.Results.tail; 0362 if cfg.tail, 0363 cfg.alpha = p.Results.alpha; 0364 else 0365 cfg.alpha = p.Results.alpha/2; %for a two-tailed test FWER=cfg.alpha*2; presumably this just effects the mask 0366 %this appears to still be accurate despite the changes Fieldtrip made 0367 %to stats{#}.negclusters and stats{#}.posclusters 0368 end 0369 cfg.statistic = 'depsamplesT'; 0370 cfg.correctm = p.Results.correctm; 0371 cfg.correcttail = 'prob'; %this makes the p-values of ft_freqstatistics reflect the tail of the analysis 0372 % (i.e., the p-values in tfrGND.stats{neo_test}.prob). Note that the 0373 % p-values in tfrGND.stats{neo_test}.negclusters and 0374 % tfrGND.stats{neo_test}.posclusters are one-tailed p-values 0375 0376 cfg.avgoverchan=p.Results.avgoverchan; 0377 cfg.avgovertime=p.Results.avgovertime; 0378 cfg.avgoverfreq=p.Results.avgoverfreq; 0379 0380 if isempty(p.Results.time_limits), 0381 cfg.latency=[tfrGND.ftrip.time(1) tfrGND.ftrip.time(end)]; 0382 else 0383 cfg.latency=p.Results.time_limits/1000; %convert from ms to sec 0384 %error check 0385 if max(cfg.latency)>max(tfrGND.ftrip.time) 0386 error('Specified maximum time limit %g ms exceeds data (i.e., tfrGND) max time limit %g ms.', ... 0387 max(cfg.latency)*1000,max(tfrGND.ftrip.time)*1000); 0388 elseif min(cfg.latency)<min(tfrGND.ftrip.time) 0389 error('Specified minimum time limit %g ms exceeds data (i.e., tfrGND) min time limit %g ms.', ... 0390 min(cfg.latency)*1000,min(tfrGND.ftrip.time)*1000); 0391 end 0392 end 0393 if isempty(p.Results.freq_limits), 0394 cfg.frequency=[tfrGND.ftrip.freq(1) tfrGND.ftrip.freq(end)]; 0395 %error check 0396 if max(cfg.frequency)>max(tfrGND.ftrip.freq) 0397 error('Specified maximum time limit %g Hz exceeds data (i.e., tfrGND) max time limit %g Hz.', ... 0398 max(cfg.frequency),max(tfrGND.ftrip.freq)); 0399 elseif min(cfg.frequency)<min(tfrGND.ftrip.freq) 0400 error('Specified minimum time limit %g Hz exceeds data (i.e., tfrGND) min time limit %g Hz.', ... 0401 min(cfg.frequency),min(tfrGND.ftrip.freq)); 0402 end 0403 else 0404 cfg.frequency=p.Results.freq_limits; 0405 end 0406 0407 %% Parameters that vary depending on correction for multiple comparisons 0408 if strcmpi(cfg.correctm,'cluster'), 0409 cfg.method = 'montecarlo'; 0410 cfg.clusterstatistic = p.Results.clusterstatistic; 0411 cfg.clusteralpha = p.Results.clusteralpha; 0412 cfg.minnbchan = p.Results.minnbchan; 0413 cfg.numrandomization = p.Results.n_perm; %this massively eats up memory (I can get 500 to work though) 0414 cfg.neighbourdist=p.Results.neighbourdist; % This is approximately 5.4 cm (assuming a spherical head with 56 cm circumference) 0415 % It seems to work well for 26 channel cap 0416 fprintf('Minimal number of above threshold neighbors required for cluster membership (i.e., minnbchan parameter of ft_freqstatistics): %d\n', ... 0417 cfg.minnbchan); 0418 elseif strcmpi(cfg.correctm,'fdr_bh') || strcmpi(cfg.correctm,'fdr_by') || strcmpi(cfg.correctm,'none'), 0419 cfg.method = 'analytic'; %this won't correct for multiple comparisons, but we can do that later 0420 0421 %Alternatively, you could use the code below that uses a permutation 0422 %test to get test-wise p-values 0423 %cfg.method = 'montecarlo'; 0424 %cfg.correctm = 'fdr'; 0425 %cfg.numrandomization = p.Results.n_perm; 0426 elseif strcmpi(cfg.correctm,'max'), 0427 cfg.method = 'montecarlo'; 0428 cfg.numrandomization = p.Results.n_perm; 0429 else 0430 error('%s is not currently a supported method for multiple comparison correction.', ... 0431 p.Results.correctm); 0432 end 0433 0434 %% Create design matrix 0435 n_subs=size(tfrGND.indiv{1},1); 0436 design = zeros(2,2*n_subs); 0437 for i = 1:n_subs 0438 design(1,i) = i; 0439 end 0440 for i = 1:n_subs 0441 design(1,n_subs+i) = i; 0442 end 0443 design(2,1:n_subs) = 1; 0444 design(2,n_subs+1:2*n_subs) = 2; 0445 0446 cfg.design = design; 0447 cfg.uvar = 1; 0448 cfg.ivar = 2; 0449 0450 0451 %% Extract spectrograms from desired bins into dummy fieldtrip variables & 0452 % perform t-tests 0453 cond1=tfrGND.ftrip; 0454 cond1.dimord='subj_chan_freq_time'; 0455 cond1.powspctrm=tfrGND.indiv{binA}; 0456 0457 cond2=tfrGND.ftrip; 0458 cond2.dimord='subj_chan_freq_time'; 0459 cond2.powspctrm=tfrGND.indiv{binB}; 0460 n_tests=length(tfrGND.stats); 0461 neo_test=n_tests+1; 0462 tfrGND.stats{neo_test}=ft_freqstatistics(cfg,cond1,cond2); %perform the test 0463 tfrGND.stats{neo_test}.correctm=p.Results.correctm; 0464 tfrGND.stats{neo_test}.family_alpha=p.Results.alpha; 0465 tfrGND.stats{neo_test}.bins=[num2str(binA) '-' num2str(binB)]; 0466 0467 0468 %% Command line summary of results 0469 if strcmpi(cfg.correctm,'cluster'), 0470 if isfield(tfrGND.stats{neo_test},'posclusters') 0471 npc=length(tfrGND.stats{neo_test}.posclusters); 0472 sig_p=0; 0473 if npc, 0474 pc_pvals=ones(1,npc); 0475 end 0476 for a=1:npc, 0477 pc_pvals(a)=tfrGND.stats{neo_test}.posclusters(a).prob; 0478 sig_p=sig_p+(tfrGND.stats{neo_test}.posclusters(a).prob<=p.Results.alpha); 0479 end 0480 fprintf('Number of positive clusters found: %d\n',npc); 0481 fprintf('Number of *significant* positive clusters: %d\n',sig_p); 0482 if npc 0483 fprintf('Min/max positive cluster p-value: %g/%g\n',min(pc_pvals),max(pc_pvals)); 0484 end 0485 end 0486 0487 if isfield(tfrGND.stats{neo_test},'negclusters') 0488 nnc=length(tfrGND.stats{neo_test}.negclusters); 0489 sig_n=0; 0490 if nnc, 0491 nc_pvals=ones(1,nnc); 0492 end 0493 for a=1:nnc, 0494 sig_n=sig_n+(tfrGND.stats{neo_test}.negclusters(a).prob<=p.Results.alpha); 0495 nc_pvals(a)=tfrGND.stats{neo_test}.negclusters(a).prob; 0496 end 0497 fprintf('Number of negative clusters found: %d\n',nnc); 0498 fprintf('Number of *significant* negative clusters: %d\n',sig_n); 0499 if nnc 0500 fprintf('Min/max negative cluster p-value: %g/%g\n',min(nc_pvals),max(nc_pvals)); 0501 end 0502 end 0503 elseif strcmpi(cfg.correctm,'fdr_bh'), 0504 %correct p-values using Benjamini-Hochberg method. Note adjusted p-values may be 0505 %greater than 1. 0506 fprintf('NOW correcting for multiple comparisons with Benjamini & Hochberg FDR control algorithm.\n'); 0507 [h, crit_p, tfrGND.stats{neo_test}.prob]=fdr_bh(tfrGND.stats{neo_test}.prob); 0508 tfrGND.stats{neo_test}.mask=tfrGND.stats{neo_test}.prob<=p.Results.alpha; 0509 clear h crit_p; %not needed 0510 elseif strcmpi(cfg.correctm,'fdr_by'), 0511 %correct p-values using Benjamini-Yekutieli method. Note adjusted p-values may be 0512 %greater than 1. 0513 fprintf('NOW correcting for multiple comparisons with Benjamini & Yekutieli FDR control algorithm.\n'); 0514 [h, crit_p, tfrGND.stats{neo_test}.prob]=fdr_bh(tfrGND.stats{neo_test}.prob,.05,'dep'); 0515 tfrGND.stats{neo_test}.mask=tfrGND.stats{neo_test}.prob<=p.Results.alpha; 0516 clear h crit_p; %not needed 0517 end 0518 0519 0520 if ~strcmpi(cfg.correctm,'cluster'), 0521 if strcmpi(cfg.correctm(1:3),'fdr'), 0522 pORq='q'; 0523 else 0524 pORq='p'; 0525 end 0526 0527 %summarize results 0528 mn_p=min(min(min(tfrGND.stats{neo_test}.prob))); 0529 mx_p=max(max(max(tfrGND.stats{neo_test}.prob))); 0530 fprintf('Min/Max %c-values: %g/%g\n',pORq,mn_p,mx_p); 0531 sig_ids=find(tfrGND.stats{neo_test}.prob<=p.Results.alpha); 0532 if isempty(sig_ids), 0533 fprintf('No comparisons are significant (%c<=%g)\n',pORq,p.Results.alpha); 0534 else 0535 mn_p=min(min(min(tfrGND.stats{neo_test}.prob(sig_ids)))); 0536 mx_p=max(max(max(tfrGND.stats{neo_test}.prob(sig_ids)))); 0537 fprintf('Min/max significant %c-values: %g/%g\n',pORq,mn_p,mx_p); 0538 end 0539 end 0540 0541 %tag on baseline and unit information 0542 tfrGND.stats{neo_test}.units=units; 0543 tfrGND.stats{neo_test}.bsln_wind=bsln_wind; 0544 tfrGND.stats{neo_test}.bsln_type=bsln_type; 0545 0546 %Indicate tfrGND variable hasn't been saved since last change 0547 tfrGND.saved='no'; 0548