


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

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