Home > chronux_1_1 > statistical_tests > two_group_permutation_test.m

two_group_permutation_test

PURPOSE ^

Function to do a two group permutation test to determine the distribution

SYNOPSIS ^

function lstatperm=two_group_permutation_test(x,y1,y2,Nperm,family,nn,h)

DESCRIPTION ^

 Function to do a two group permutation test to determine the distribution
 of the statistic lstat based on locfit smoothing
 Usage: lstatperm=two_group_permutation_test(data1,data2,lstat,Nperm)
 Inputs:
 x: independent variable (at which the fits to y1 and y2 are computed
 using locfit) as a column vector
 y1, y2: data sets - BOTH struct arrays with length equal to number of trials
                     or BOTH matrices where data corresponding to each trial is
                     stored in a column - required
 lstat: the statistic to be computed - required
 Takes values 'd' for lstat=smooth(y1)-smooth(y2)
                 'ad' for lstat=|smooth(y1)-smooth(y2)|
                 'mad' for lstat=max(|smooth(y1)-smooth(y2)|)
                 'com' for lstat=differences in the centre of mass of y1
                       and y2
      smooth denotes a locfit operation - currently supports computing
      rates, densities, and doing a likelihood estimate
 Nperm: number of permutations - required
 family: for locfit - can take values that locfit.m allows e.g. rate, dens, 
                      poisson, binomial etc - required
         (currently density type estimates are only supported when data is
         a struct array and regression type estimates are supported when
         data is a matrix)
 nn : nearest neighbor bandwidth - required
 h : fixed bandwidth - required

 Outputs:
 lstatperm:cell array contains the statistics computed for each permutation
           (exchange of samples between the two groups) for the following 4 statistics 
                            lstat=smooth(y1)-smooth(y2) 
                            lstat=|smooth(y1)-smooth(y2)|
                            lstat=max(smooth(y1)-smooth(y2))
                            lstat=differences in the centre of mass of y1
                                  and y2

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function lstatperm=two_group_permutation_test(x,y1,y2,Nperm,family,nn,h)
0002 % Function to do a two group permutation test to determine the distribution
0003 % of the statistic lstat based on locfit smoothing
0004 % Usage: lstatperm=two_group_permutation_test(data1,data2,lstat,Nperm)
0005 % Inputs:
0006 % x: independent variable (at which the fits to y1 and y2 are computed
0007 % using locfit) as a column vector
0008 % y1, y2: data sets - BOTH struct arrays with length equal to number of trials
0009 %                     or BOTH matrices where data corresponding to each trial is
0010 %                     stored in a column - required
0011 % lstat: the statistic to be computed - required
0012 % Takes values 'd' for lstat=smooth(y1)-smooth(y2)
0013 %                 'ad' for lstat=|smooth(y1)-smooth(y2)|
0014 %                 'mad' for lstat=max(|smooth(y1)-smooth(y2)|)
0015 %                 'com' for lstat=differences in the centre of mass of y1
0016 %                       and y2
0017 %      smooth denotes a locfit operation - currently supports computing
0018 %      rates, densities, and doing a likelihood estimate
0019 % Nperm: number of permutations - required
0020 % family: for locfit - can take values that locfit.m allows e.g. rate, dens,
0021 %                      poisson, binomial etc - required
0022 %         (currently density type estimates are only supported when data is
0023 %         a struct array and regression type estimates are supported when
0024 %         data is a matrix)
0025 % nn : nearest neighbor bandwidth - required
0026 % h : fixed bandwidth - required
0027 %
0028 % Outputs:
0029 % lstatperm:cell array contains the statistics computed for each permutation
0030 %           (exchange of samples between the two groups) for the following 4 statistics
0031 %                            lstat=smooth(y1)-smooth(y2)
0032 %                            lstat=|smooth(y1)-smooth(y2)|
0033 %                            lstat=max(smooth(y1)-smooth(y2))
0034 %                            lstat=differences in the centre of mass of y1
0035 %                                  and y2
0036 x=x(:);
0037 if isstruct(y1) && isstruct(y2); 
0038    Ntr1=length(y1);
0039    Ntr2=length(y2);
0040    N=Ntr1+Ntr2;
0041    y=[y1 y2];
0042    lstatperm=cell(1,4);
0043    for n=1:Nperm;
0044        fprintf('%d\n',n);
0045        vv=randperm(N);
0046        v=vv(1:Ntr1);
0047        w=vv(Ntr1+1:end);
0048        dtmp1=y(v);
0049        dtmp2=y(w);
0050        %
0051        % Convert the struct arrays into 1d arrays
0052        %
0053        data1=[];
0054        for n=1:Ntr1;
0055            data1=[data1;dtmp1(n).times];
0056        end;
0057 
0058        data2=[];
0059        for n=1:Ntr2;
0060            data2=[data2;dtmp2(n).times];
0061        end;
0062        %
0063        % Compute the fits
0064        %
0065        fit1=locfit(data1,'nn',nn,'h',h,'family',family);
0066        fali1=fit1{4}{5};
0067        fit1=predict(fit1,x);
0068        fit1=invlink(fit1,fali1);
0069 
0070        fit2=locfit(data2,'nn',nn,'h',h,'family',family);
0071        fali2=fit2{4}{5};
0072        fit2=predict(fit2,x);
0073        fit2=invlink(fit2,fali2);
0074               
0075        if strcmp(family,'rate'); fit1=fit1/Ntr1; fit2=fit2/Ntr2; end;
0076        %
0077        % Compute the statistics
0078        %
0079        lstatperm{1}=[lstatperm{1} fit1-fit2];
0080        lstatperm{2}=[lstatperm{2} abs(fit1-fit2)];
0081        lstatperm{3}=[lstatperm{3} max(abs(fit1-fit2))];
0082 %        com1=trapz(x,x.*fit1)/trapz(x,fit1);
0083 %        com2=trapz(x,x.*fit2)/trapz(x,fit2);
0084        com1=trapz(x,x.*x.*fit1);
0085        com2=trapz(x,x.*x.*fit2);
0086        lstatperm{4}=[lstatperm{4} com1-com2];
0087    end;
0088 elseif ~isstruct(y1) && ~isstruct(y2);
0089     Ntr1=size(y1,2);
0090     Ntr2=size(y2,2);
0091     N=Ntr1+Ntr2;
0092     y=[y1 y2];
0093     x1=repmat(x,[1 Ntr1]);
0094     x2=repmat(x,[1 Ntr2]);
0095     x1=x1';
0096     x2=x2';
0097     x1=reshape(x1,[prod(size(x1)) 1]);
0098     x2=reshape(x2,[prod(size(x2)) 1]);
0099     lstatperm=cell(1,4);
0100     for n=1:Nperm;
0101         fprintf('%d\n',n);
0102         vv=randperm(N);
0103         v=vv(1:Ntr1);
0104         w=vv(Ntr1+1:end);
0105         data1=y(:,v)';
0106         data2=y(:,w)';
0107         data1=reshape(data1,[prod(size(data1)),1]);
0108         data2=reshape(data2,[prod(size(data2)),1]);
0109         %
0110         % Compute the fits
0111         %
0112         fit1=locfit(x1,data1,'nn',nn,'h',h,'family',family);
0113         fali1=fit1{4}{5};
0114         fit1=predict(fit1,x);
0115         fit1=invlink(fit1,fali1);
0116 
0117         fit2=locfit(x2,data2,'nn',nn,'h',h,'family',family);
0118         fali2=fit2{4}{5};
0119         fit2=predict(fit2,x);
0120         fit2=invlink(fit2,fali2);
0121         %
0122         % Compute the statistics
0123         %
0124        lstatperm{1}=[lstatperm{1} fit1-fit2];
0125        lstatperm{2}=[lstatperm{2} abs(fit1-fit2)];
0126        lstatperm{3}=[lstatperm{3} max(abs(fit1-fit2))];
0127 %        com1=trapz(x,x.*fit1)/trapz(x,fit1);
0128 %        com2=trapz(x,x.*fit2)/trapz(x,fit2);
0129        com1=trapz(x,x.*x.*fit1);
0130        com2=trapz(x,x.*x.*fit2);
0131        lstatperm{4}=[lstatperm{4} com1-com2];
0132     end;
0133 else;
0134     error('BOTH data sets should be struct arrays or matrices');
0135 end;
0136 
0137

Generated on Sun 13-Aug-2006 11:49:44 by m2html © 2003