Home > matlabmk > mxt_perm1old_matlab.m

mxt_perm1old_matlab

PURPOSE ^

function [pval, t_orig, tmx_ptile, seed_state, est_alpha]=mxt_perm1(data,n_perm,alph,tail,verblevel,seed_state,freq_domain)

SYNOPSIS ^

function [pval, t_orig, tmx_ptile, seed_state, est_alpha]=mxt_perm1(data,n_perm,alph,tail,verblevel,seed_state,freq_domain)

DESCRIPTION ^

function [pval, t_orig, tmx_ptile, seed_state, est_alpha]=mxt_perm1(data,n_perm,alph,tail,verblevel,seed_state,freq_domain)

 mxt_perm1-One sample permutation test based on a t-statistic and null
 hypothesis of a mean of zero.  Can handle multiple electrodes and time
 points/frequencies. According to Hemmelammnn et al. (2004), the t_max permutation
 test has relatively good power for multivariate data whose dimensions are
 highly correlated (like EEG/ERP data).  See Maris (2004) for an
 explanation of how permutation tests like this control for multiple
 comparisons.  Note, Maris uses Hotelling's multivariate T^2 statistic in
 his paper (instead of the t_max used here) and he advocated performing
 the tests over multiple time points (e.g., a 10 ms time window).
 However, our lab found that the T^2 statistic did not have much power when the
 sample size was near the number of electrodes and the results of the test
 could vary greatly depending on the size of the time window.

 Inputs:
  data   - 3D matrix of data (Channel x Time x Participant)

 Optional Inputs:
  n_perm      - number of permutations {default=2000}.  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.
  alph        - test alpha level {default=.05}
  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.
  freq_domain - If 0, command line report will be given in temporal units
                (e.g. time points).  Otherwise, the report will be given 
                in frequency domain units (e.g., frequencies). {default:
                0}
  

 Outputs:
  pval       - p-value at each time point and electrode (corrected for
               multiple comparisons via permutation test)
  t_orig     - t-score at each time point and electrode
  tmx_ptile  - critical t-scores for given alpha level
  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.
  est_alpha  - The estimated family-wise alpha level of the test.  With 
               permutation tests, a finite number of p-values are possible.
               This function tries to use an alpha level that is as close 
               as possible to the desired alpha level.  However, if the 
               sample size is small, a very limited number of p-values are 
               possible and the desired family-wise alpha level may be 
               impossible to approximately achieve.
  

 Notes:
 -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 rapidly with the number of
 participants, this is only an 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 overly 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).

 Author:
 David Groppe
 June, 2009
 Kutaslab, San Diego

 References:

 Hemmelmann, et al. (2004) Multivariate tests for the evaluation of
 high-dimensional EEG data. Journal of Neuroscience Methods.

 Maris, E. (2004) Randomization tests for ERP topographies and whole 
 spatiotemporal data matrices. Psychophysiology

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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 %function [pval, t_orig, tmx_ptile, seed_state, est_alpha]=mxt_perm1(data,n_perm,alph,tail,verblevel,seed_state,freq_domain)
0002 %
0003 % mxt_perm1-One sample permutation test based on a t-statistic and null
0004 % hypothesis of a mean of zero.  Can handle multiple electrodes and time
0005 % points/frequencies. According to Hemmelammnn et al. (2004), the t_max permutation
0006 % test has relatively good power for multivariate data whose dimensions are
0007 % highly correlated (like EEG/ERP data).  See Maris (2004) for an
0008 % explanation of how permutation tests like this control for multiple
0009 % comparisons.  Note, Maris uses Hotelling's multivariate T^2 statistic in
0010 % his paper (instead of the t_max used here) and he advocated performing
0011 % the tests over multiple time points (e.g., a 10 ms time window).
0012 % However, our lab found that the T^2 statistic did not have much power when the
0013 % sample size was near the number of electrodes and the results of the test
0014 % could vary greatly depending on the size of the time window.
0015 %
0016 % Inputs:
0017 %  data   - 3D matrix of data (Channel x Time x Participant)
0018 %
0019 % Optional Inputs:
0020 %  n_perm      - number of permutations {default=2000}.  Manly (1997) suggests
0021 %                using at least 1000 permutations for an alpha level of 0.05 and
0022 %                at least 5000 permutations for an alpha level of 0.01.
0023 %  alph        - test alpha level {default=.05}
0024 %  tail        - [1 | 0 | -1] If tail=1, the alternative hypothesis is that the
0025 %                mean of the data is greater than 0 (upper tailed test).  If tail=0,
0026 %                the alternative hypothesis is that the mean of the data is different
0027 %                than 0 (two tailed test).  If tail=-1, the alternative hypothesis
0028 %                is that the mean of the data is less than 0 (lower tailed test).
0029 %                {default: 0}
0030 %  verblevel   - An integer specifying the amount of information you want
0031 %                this function to provide about what it is doing during runtime.
0032 %                 Options are:
0033 %                    0 - quiet, only show errors, warnings, and EEGLAB reports
0034 %                    1 - stuff anyone should probably know
0035 %                    2 - stuff you should know the first time you start working
0036 %                        with a data set {default value}
0037 %                    3 - stuff that might help you debug (show all
0038 %                        reports)
0039 %  seed_state  - The initial state of the random number generating stream
0040 %                (see MATLAB documentation for "randstream"). If you pass
0041 %                a value from a previous run of this function, it should
0042 %                reproduce the exact same values.
0043 %  freq_domain - If 0, command line report will be given in temporal units
0044 %                (e.g. time points).  Otherwise, the report will be given
0045 %                in frequency domain units (e.g., frequencies). {default:
0046 %                0}
0047 %
0048 %
0049 % Outputs:
0050 %  pval       - p-value at each time point and electrode (corrected for
0051 %               multiple comparisons via permutation test)
0052 %  t_orig     - t-score at each time point and electrode
0053 %  tmx_ptile  - critical t-scores for given alpha level
0054 %  seed_state - The initial state of the random number generating stream
0055 %               (see MATLAB documentation for "randstream") used to
0056 %               generate the permutations. You can use this to reproduce
0057 %               the output of this function.
0058 %  est_alpha  - The estimated family-wise alpha level of the test.  With
0059 %               permutation tests, a finite number of p-values are possible.
0060 %               This function tries to use an alpha level that is as close
0061 %               as possible to the desired alpha level.  However, if the
0062 %               sample size is small, a very limited number of p-values are
0063 %               possible and the desired family-wise alpha level may be
0064 %               impossible to approximately achieve.
0065 %
0066 %
0067 % Notes:
0068 % -Unlike a parametric test (e.g., an ANOVA), a discrete set of p-values
0069 % are possible (at most the number of possible permutations).  Since the
0070 % number of possible permutations grows rapidly with the number of
0071 % participants, this is only an issue for small sample sizes (e.g., 6
0072 % participants).  When you have such a small sample size, the
0073 % limited number of p-values may make the test overly conservative (e.g.,
0074 % you might be forced to use an alpha level of .0286 since it is the biggest
0075 % possible alpha level less than .05).
0076 %
0077 % Author:
0078 % David Groppe
0079 % June, 2009
0080 % Kutaslab, San Diego
0081 %
0082 % References:
0083 %
0084 % Hemmelmann, et al. (2004) Multivariate tests for the evaluation of
0085 % high-dimensional EEG data. Journal of Neuroscience Methods.
0086 %
0087 % Maris, E. (2004) Randomization tests for ERP topographies and whole
0088 % spatiotemporal data matrices. Psychophysiology
0089 %
0090 % Manly, B.F.J. (1997) Randomization, bootstrap, and Monte Carlo methods in
0091 % biology. 2nd ed. Chapmn and Hall, London.
0092 
0093 %%%%%%%%%%%%%%%% REVISION LOG %%%%%%%%%%%%%%%%%
0094 %
0095 % 3/8/2010-rand_state optional input and output added
0096 %
0097 % 3/18/2010-modified to be much faster (time for-loop eliminated, some numbers needed to calculate
0098 % t are computed once for all t-tests; "computational formula" used instead of
0099 % conceptual formula).  New version gives exact same results as old version
0100 %
0101 % 4/1/2010-Estimated alpha level added
0102 %
0103 % 12/10/2010-Differences in the steps used for computing t-scores between
0104 % permutations and observed data could lead to rounding errors that would
0105 % produce inaccurate p-values for very small sample sizes (e.g., 3
0106 % participants).  That's been fixed.
0107 %
0108 % 12/14/2010-freq_domain optional input added
0109 
0110 %This code has an appropriate false positive rate when run on simulated
0111 %data at one time point (use sim_test1.m):
0112 %
0113 % n_sim=4000, n_perm=2000, data=2x1x16 Gaussian, alpha=.05
0114 % Two-tailed test: Empirical alpha rate: 0.047500
0115 % Upper tailed test: Empirical alpha rate: 0.048500
0116 % Lower tailed test: Empirical alpha rate: 0.045500
0117 
0118 function [pval, t_orig, tmx_ptile, seed_state, est_alpha]=mxt_perm1(data,n_perm,alph,tail,verblevel,seed_state,freq_domain)
0119 
0120 if nargin<1,
0121     error('You need to provide data.');
0122 end
0123 
0124 if nargin<2,
0125     n_perm=2000;
0126 end
0127 
0128 if nargin<3,
0129     alph=.05;
0130 elseif (alph>=1) || (alph<=0)
0131     error('Argument ''alph'' needs to be between 0 and 1.');
0132 end
0133 
0134 if alph<=.01 && n_perm<5000,
0135     watchit(sprintf('You are probably using too few permutations for an alpha level of %f. Type ">>help mxt_perm1" for more info.',alph));
0136 elseif alph<=.05 && n_perm<1000,
0137     watchit(sprintf('You are probably using too few permutations for an alpha level of %f. Type ">>help mxt_perm1" for more info.',alph));
0138 end
0139 
0140 if nargin<4,
0141     tail=0;
0142 elseif (tail~=0) && (tail~=1) && (tail~=-1),
0143     error('Argument ''tail'' needs to be 0,1, or -1.');
0144 end
0145 
0146 if nargin<5,
0147     verblevel=2;
0148 end
0149 
0150 if 0 % ?? for old versions of Matlab
0151     defaultStream=RandStream.getDefaultStream; %random # generator state
0152     if (nargin<6) || isempty(seed_state),
0153         %Store state of random number generator
0154         seed_state=defaultStream.State;
0155     else
0156         defaultStream.State=seed_state; %reset random number generator to saved state
0157     end
0158 else
0159    fprintf('Not seeding random number generator.\n'); 
0160    seed_state=NaN;
0161 end
0162 
0163 if (nargin<7),
0164    freq_domain=0; 
0165 end
0166 
0167 s=size(data);
0168 n_chan=s(1);
0169 n_pts=s(2);
0170 n_subs=s(3);
0171 if n_subs<2,
0172     error('You need data from at least two observations (e.g., participants) to perform a hypothesis test.')
0173 end
0174 
0175 if n_subs<7,
0176     n_psbl_prms=2^n_subs;
0177     watchit(sprintf(['Due to the very limited number of participants,' ...
0178         ' 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.'], ...
0179     n_psbl_prms));
0180 end
0181 
0182 if verblevel~=0,
0183     fprintf('mxt_perm1: Number of channels: %d\n',n_chan);
0184     if freq_domain,
0185         fprintf('mxt_perm1: Number of frequencies: %d\n',n_pts);
0186     else
0187         fprintf('mxt_perm1: Number of time points: %d\n',n_pts);
0188     end
0189     fprintf('mxt_perm1: Total # of comparisons: %d\n',n_pts*n_chan);
0190     fprintf('mxt_perm1: Number of participants: %d\n',n_subs);
0191     fprintf('t-score degrees of freedom: %d\n',n_subs-1);
0192 end
0193 
0194 
0195 VerbReport(sprintf('Executing permutation test with %d permutations...',n_perm),2,verblevel);
0196 if (verblevel>=2),
0197     fprintf('Permutations completed: ');
0198 end
0199 
0200 
0201 %Constant factor for computing t, speeds up computing t to precalculate
0202 %now
0203 sqrt_nXnM1=sqrt(n_subs*(n_subs-1));
0204 
0205 mxt=zeros(1,n_perm*2);
0206 sn=zeros(1,1,n_subs);
0207 for perm=1:n_perm
0208     if ~rem(perm,100)
0209         if (verblevel>=2),
0210             if ~rem(perm-100,1000)
0211                 fprintf('%d',perm);
0212             else
0213                 fprintf(', %d',perm);
0214             end
0215             if ~rem(perm,1000)
0216                 fprintf('\n');
0217             end
0218         end
0219     end
0220     %randomly set sign of each participant's data
0221     sn(1,1,1:n_subs)=(rand(1,n_subs)>.5)*2-1;
0222     sn_mtrx=repmat(sn,[n_chan n_pts 1]);
0223   
0224     d_perm=data.*sn_mtrx;
0225         
0226     %computes t-score of permuted data across all channels and time
0227     %points or frequencies
0228     sm=sum(d_perm,3);
0229     mn=sm/n_subs;
0230     sm_sqrs=sum(d_perm.^2,3)-(sm.^2)/n_subs;
0231     stder=sqrt(sm_sqrs)/sqrt_nXnM1;
0232     t=mn./stder;
0233     
0234     %get most extreme t-score (sign isn't immportant since we asumme
0235     %symmetric distribution of null hypothesis for one sample test)
0236     mxt(perm)=max(max(abs(t)));
0237 end
0238 mxt(n_perm+1:2*n_perm)=-mxt(1:n_perm); %add the negative of all values since we assumme
0239 %null hypothesis distribution is symmetric
0240 
0241 %End permutations completed line
0242 if (verblevel>=2) && rem(perm,1000)
0243     fprintf('\n');
0244 end
0245 
0246 if tail==0,
0247     %two-tailed
0248     tmx_ptile(1)=prctile(mxt,100*alph/2);
0249     tmx_ptile(2)=-tmx_ptile(1);
0250     est_alpha=mean(mxt>=tmx_ptile(2))*2;
0251 elseif tail==1,
0252     %upper tailed
0253     tmx_ptile=prctile(mxt,100-100*alph);
0254     est_alpha=mean(mxt>=tmx_ptile);
0255 else
0256     %tail=-1, lower tailed
0257     tmx_ptile=prctile(mxt,alph*100);
0258     est_alpha=mean(mxt<=tmx_ptile);
0259 end
0260 if verblevel~=0,
0261     fprintf('Desired family-wise alpha level: %f\n',alph);
0262     fprintf('Estimated actual family-wise alpha level: %f\n',est_alpha);
0263 end
0264 
0265 %computes t-scores of observations at all channels and time
0266 %points/frequencies
0267 sm=sum(data,3);
0268 mn=sm/n_subs;
0269 sm_sqrs=sum(data.^2,3)-(sm.^2)/n_subs;
0270 stder=sqrt(sm_sqrs)/sqrt_nXnM1;
0271 t_orig=mn./stder;
0272 
0273 %compute p-values
0274 pval=zeros(n_chan,n_pts);
0275 for t=1:n_pts,
0276     for c=1:n_chan,
0277         if tail==0,
0278             pval(c,t)=mean(mxt>=abs(t_orig(c,t)))*2;
0279         elseif tail==1,
0280             pval(c,t)=mean(mxt>=t_orig(c,t));
0281         elseif tail==-1,
0282             pval(c,t)=mean(mxt<=t_orig(c,t));
0283         end
0284     end
0285 end
0286

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