korn_fd1-Korn et al. (2004) permutation based procedure for
controlling the "generalized family-wise error rate" (GFWER) of a
family of repeated-measure/one-sample t-tests. This
function implements the less computationally intensive/more
conservative GFWER procedure described by Korn and colleagues
(it is called Procedure A* in the appendix of their paper).
Usage:
>> [p_adj, p_raw, t_raw, seed_state]=korn_fd1(data,n_perm,n_fd,tail,verblevel,seed_state);
Required Input:
data - A 3D matrix of ERP or difference wave data (channel x time point x
participant)
Optional Inputs:
n_perm - Number of permutations to use to estimate the null
hypothesis distribution of all possible permutations.
{default: 2500}
n_fd - The lower bound on the number of false discoveries
whose probability you want to control. For example, if
n_fd is 1 and you reject all null hypotheses with
adjusted p-values less than .05, then the probably that the
number of false discoveries exceeds 1 is 5% or less.
{default: 1}
tail - [1 | 0 | -1] If tail=1, the alternative hypothesis is that the
mean of the data is greater than 0 (upper tailed test). If tail=0,
the alternative hypothesis is that the mean of the data is different
than 0 (two tailed test). If tail=-1, the alternative hypothesis
is that the mean of the data is less than 0 (lower tailed test).
{default: 0}
verblevel - An integer specifying the amount of information you want
this function to provide about what it is doing during runtime.
Options are:
0 - quiet, only show errors, warnings, and EEGLAB reports
1 - stuff anyone should probably know
2 - stuff you should know the first time you start working
with a data set {default value}
3 - stuff that might help you debug (show all
reports)
seed_state - The initial state of the random number generating stream
(see MATLAB documentation for "randstream"). If you pass
a value from a previous run of this function, it should
reproduce the exact same values.
Outputs:
p_adj - p-values corresponding to each dimension of the hypothesis
test after correcting for multiple comparisons
p_raw - p-values corresponding to each dimension of the hypothesis
test WITHOUT correction for multiple comparisons
t_raw - t-scores corresponding to each dimension of the hypothesis
test
seed_state - The initial state of the random number generating stream
(see MATLAB documentation for "randstream") used to
generate the permutations. You can use this to reproduce
the output of this function.
Author:
David Groppe
August, 2009
Kutaslab, San Diego
Reference:
Korn, E. L., Troendle, J. F., McShane, L. M., & Simon, R. (2004).
Controlling the number of false discoveries: Application to high-
dimensional genomic data. Journal of Statistical Planning and Inference,
124, 379-398.0001 % korn_fd1-Korn et al. (2004) permutation based procedure for 0002 % controlling the "generalized family-wise error rate" (GFWER) of a 0003 % family of repeated-measure/one-sample t-tests. This 0004 % function implements the less computationally intensive/more 0005 % conservative GFWER procedure described by Korn and colleagues 0006 % (it is called Procedure A* in the appendix of their paper). 0007 % 0008 % Usage: 0009 % >> [p_adj, p_raw, t_raw, seed_state]=korn_fd1(data,n_perm,n_fd,tail,verblevel,seed_state); 0010 % 0011 % 0012 % Required Input: 0013 % data - A 3D matrix of ERP or difference wave data (channel x time point x 0014 % participant) 0015 % 0016 % Optional Inputs: 0017 % n_perm - Number of permutations to use to estimate the null 0018 % hypothesis distribution of all possible permutations. 0019 % {default: 2500} 0020 % n_fd - The lower bound on the number of false discoveries 0021 % whose probability you want to control. For example, if 0022 % n_fd is 1 and you reject all null hypotheses with 0023 % adjusted p-values less than .05, then the probably that the 0024 % number of false discoveries exceeds 1 is 5% or less. 0025 % {default: 1} 0026 % tail - [1 | 0 | -1] If tail=1, the alternative hypothesis is that the 0027 % mean of the data is greater than 0 (upper tailed test). If tail=0, 0028 % the alternative hypothesis is that the mean of the data is different 0029 % than 0 (two tailed test). If tail=-1, the alternative hypothesis 0030 % is that the mean of the data is less than 0 (lower tailed test). 0031 % {default: 0} 0032 % verblevel - An integer specifying the amount of information you want 0033 % this function to provide about what it is doing during runtime. 0034 % Options are: 0035 % 0 - quiet, only show errors, warnings, and EEGLAB reports 0036 % 1 - stuff anyone should probably know 0037 % 2 - stuff you should know the first time you start working 0038 % with a data set {default value} 0039 % 3 - stuff that might help you debug (show all 0040 % reports) 0041 % seed_state - The initial state of the random number generating stream 0042 % (see MATLAB documentation for "randstream"). If you pass 0043 % a value from a previous run of this function, it should 0044 % reproduce the exact same values. 0045 % 0046 % Outputs: 0047 % p_adj - p-values corresponding to each dimension of the hypothesis 0048 % test after correcting for multiple comparisons 0049 % p_raw - p-values corresponding to each dimension of the hypothesis 0050 % test WITHOUT correction for multiple comparisons 0051 % t_raw - t-scores corresponding to each dimension of the hypothesis 0052 % test 0053 % seed_state - The initial state of the random number generating stream 0054 % (see MATLAB documentation for "randstream") used to 0055 % generate the permutations. You can use this to reproduce 0056 % the output of this function. 0057 % 0058 % Author: 0059 % David Groppe 0060 % August, 2009 0061 % Kutaslab, San Diego 0062 % 0063 % Reference: 0064 % Korn, E. L., Troendle, J. F., McShane, L. M., & Simon, R. (2004). 0065 % Controlling the number of false discoveries: Application to high- 0066 % dimensional genomic data. Journal of Statistical Planning and Inference, 0067 % 124, 379-398. 0068 % 0069 0070 0071 %%%%%%%%%%%%%%%% REVISION LOG %%%%%%%%%%%%%%%%% 0072 % 0073 % ?/?/2010- 0074 0075 function [p_adj, p_raw, t_raw, seed_state]=korn_fd1(data,n_perm,n_fd,tail,verblevel,seed_state) 0076 0077 if nargin<1, 0078 error('You need to provide data.'); 0079 end 0080 0081 if nargin<2. 0082 n_perm=2500; 0083 end 0084 0085 if nargin<3, 0086 n_fd=1; 0087 end 0088 0089 if nargin<4, 0090 tail=0; 0091 elseif (tail~=0) && (tail~=1) && (tail~=-1), 0092 error('Argument ''tail'' needs to be 0,1, or -1.'); 0093 end 0094 0095 if nargin<5, 0096 verblevel=2; 0097 end 0098 0099 if (nargin<6) || isempty(seed_state), 0100 seed_state=[]; 0101 end 0102 0103 [n_chan, n_tpt, n_sub]=size(data); 0104 0105 % if n_sub<7, 0106 % n_psbl_prms=2^n_sub; 0107 % watchit(sprintf(['Due to the very limited number of participants,' ... 0108 % ' the total number of possible permutations is small.\nThus only a limited number of p-values (at most %d) are possible and the test might be overly conservative.'], ... 0109 % n_psbl_prms)); 0110 % end 0111 0112 k=n_chan*n_tpt; %# of variables 0113 0114 %Step 0 <-Note, Step #'s follow those in Korn et al.'s appendix 0115 count=ones(1,k); 0116 0117 %Step 1 0118 p_tilda=zeros(1,k); 0119 0120 %Perform t-tests on n_perm permutations of the data (Step 2) 0121 [p_raw, t_raw, p_perms, seed_state]=p_perm1(data,n_perm,tail,verblevel,seed_state); 0122 0123 %Sort the p-values from permutations of the data and reshape them to be a 0124 %variable x permutation matrix 0125 p_perms=sort(reshape(p_perms,k,n_perm),1); 0126 0127 %Sort p-values of observations from smallest to largest 0128 [srtd_p, sort_ids]=sort(reshape(p_raw,1,n_chan*n_tpt)); 0129 [dummy, unsort_ids]=sort(sort_ids); %indices to return sorted_p to p_raw order 0130 0131 %Pre-allocate memory 0132 p_adj=-ones(1,k); 0133 p_adj_srtd=p_tilda; %should be all zeros 0134 0135 for r=n_fd+1:k, 0136 %Steps 3-5: Compute p_tilda 0137 perm_id=n_fd+1; 0138 count(r)=1+sum(p_perms(perm_id,:)<=srtd_p(r)); 0139 p_tilda(r)=count(r)/(n_perm+1); 0140 0141 %Note: p_tilda(r)=p_adj_srtd(r)=0 for all r<=n_fd 0142 0143 %Step 6: Compute adjusted p-value 0144 p_adj_srtd(r)=max(p_tilda(n_fd+1:r)); 0145 end 0146 p_adj=p_adj_srtd(unsort_ids); 0147 p_adj=reshape(p_adj,n_chan,n_tpt); 0148