


function [dz,vdz,Adz]=two_group_test_coherence(J1c1,J2c1,J1c2,J2c2,p)
Test the null hypothesis (H0) that data sets J1c1,J2c1,J1c2,J2c2 in
two conditions c1,c2 have equal population coherence
Usage:
[dz,vdz,Adz]=two_sample_test_coherence(J1c1,J2c1,J1c2,J2c2,p)
Inputs:
J1c1 tapered fourier transform of dataset 1 in condition 1
J2c1 tapered fourier transform of dataset 1 in condition 1
J1c2 tapered fourier transform of dataset 1 in condition 2
J2c2 tapered fourier transform of dataset 1 in condition 2
p p value for test (default: 0.05)
Dimensions: J1c1,J2c2: frequencies x number of samples in condition 1
J1c2,J2c2: frequencies x number of samples in condition 2
number of samples = number of trials x number of tapers
Outputs:
dz test statistic (will be distributed as N(0,1) under H0
vdz Arvesen estimate of the variance of dz
Adz 1/0 for accept/reject null hypothesis of equal population
coherences based dz ~ N(0,1)
Note: all outputs are functions of frequency
References: Arvesen, Jackkknifing U-statistics, Annals of Mathematical
Statisitics, vol 40, no. 6, pg 2076-2100 (1969)

0001 function [dz,vdz,Adz]=two_group_test_coherence(J1c1,J2c1,J1c2,J2c2,p) 0002 % function [dz,vdz,Adz]=two_group_test_coherence(J1c1,J2c1,J1c2,J2c2,p) 0003 % Test the null hypothesis (H0) that data sets J1c1,J2c1,J1c2,J2c2 in 0004 % two conditions c1,c2 have equal population coherence 0005 % 0006 % Usage: 0007 % [dz,vdz,Adz]=two_sample_test_coherence(J1c1,J2c1,J1c2,J2c2,p) 0008 % 0009 % Inputs: 0010 % J1c1 tapered fourier transform of dataset 1 in condition 1 0011 % J2c1 tapered fourier transform of dataset 1 in condition 1 0012 % J1c2 tapered fourier transform of dataset 1 in condition 2 0013 % J2c2 tapered fourier transform of dataset 1 in condition 2 0014 % p p value for test (default: 0.05) 0015 % 0016 % 0017 % Dimensions: J1c1,J2c2: frequencies x number of samples in condition 1 0018 % J1c2,J2c2: frequencies x number of samples in condition 2 0019 % number of samples = number of trials x number of tapers 0020 % Outputs: 0021 % dz test statistic (will be distributed as N(0,1) under H0 0022 % vdz Arvesen estimate of the variance of dz 0023 % Adz 1/0 for accept/reject null hypothesis of equal population 0024 % coherences based dz ~ N(0,1) 0025 % 0026 % Note: all outputs are functions of frequency 0027 % 0028 % References: Arvesen, Jackkknifing U-statistics, Annals of Mathematical 0029 % Statisitics, vol 40, no. 6, pg 2076-2100 (1969) 0030 0031 0032 if nargin < 4; error('Need four sets of Fourier transforms'); end; 0033 % 0034 % Test for matching dimensionalities 0035 % 0036 if size(J1c1)~=size(J2c1) | size(J1c2)~=size(J2c2) | size(J1c1,1)~=size(J1c2,1); 0037 error('Need matching dimensionalities for the Fourier transforms: Check the help file for correct dimensionalities'); 0038 else; 0039 m1=size(J1c1,2); % number of samples, condition 1 0040 m2=size(J1c2,2); % number of samples, condition 2 0041 dof1=2*m1; % number of degrees of freedom in the first condition estimates 0042 dof2=2*m2; % number of degrees of freedom in the second condition estimates 0043 end; 0044 if nargin < 5; p=0.05; end; % set the default p value 0045 0046 % 0047 % Compute the individual condition spectra, coherences 0048 % 0049 S12c1=conj(J1c1).*J2c1; % individual sample cross-spectrum, condition 1 0050 S12c2=conj(J1c2).*J2c2; % individual sample cross-spectrum, condition 2 0051 S1c1=conj(J1c1).*J1c1; % individual sample spectrum, data 1, condition 1 0052 S2c1=conj(J2c1).*J2c1; % individual sample spectrum, data 2, condition 1 0053 S1c2=conj(J1c2).*J1c2; % individual sample spectrum, data 1, condition 2 0054 S2c2=conj(J2c2).*J2c2; % individual sample spectrum, data 2, condition 2 0055 0056 Sm12c1=squeeze(mean(S12c1,2)); % mean cross spectrum, condition 1 0057 Sm12c2=squeeze(mean(S12c2,2)); % mean cross spectrum, condition 2 0058 Sm1c1=squeeze(mean(S1c1,2)); % mean spectrum, data 1, condition 1 0059 Sm2c1=squeeze(mean(S2c1,2)); % mean spectrum, data 2, condition 1 0060 Sm1c2=squeeze(mean(S1c2,2)); % mean spectrum, data 1, condition 1 0061 Sm2c2=squeeze(mean(S2c2,2)); % mean spectrum, data 2, condition 1 0062 0063 Cm12c1=abs(Sm12c1./sqrt(Sm1c1.*Sm2c1)); % mean coherence, condition 1 0064 Cm12c2=abs(Sm12c2./sqrt(Sm1c2.*Sm2c2)); % mean coherence, condition 2 0065 0066 Ccm12c1=Cm12c1; % mean coherence saved for output 0067 Ccm12c2=Cm12c2; % mean coherence saved for output 0068 % 0069 % Compute the statistic dz, and the probability of observing the value dz 0070 % given an N(0,1) distribution i.e. under the null hypothesis 0071 % 0072 z1=atanh(Cm12c1)-1/(dof1-2); % Bias-corrected Fisher z, condition 1 0073 z2=atanh(Cm12c2)-1/(dof2-2); % Bias-corrected Fisher z, condition 2 0074 dz=(z1-z2)/sqrt(1/(dof1-2)+1/(dof2-2)); % z statistic 0075 % 0076 % The remaining portion of the program computes Jackknife estimates of the mean (mdz) and variance (vdz) of dz 0077 % 0078 samples1=[1:m1]; 0079 samples2=[1:m2]; 0080 % 0081 % Leave one out of one sample 0082 % 0083 for i=1:m1; 0084 ikeep=setdiff(samples1,i); % all samples except i 0085 Sm12c1=squeeze(mean(S12c1(:,ikeep),2)); % 1 drop mean cross-spectrum, condition 1 0086 Sm1c1=squeeze(mean(S1c1(:,ikeep),2)); % 1 drop mean spectrum, data 1, condition 1 0087 Sm2c1=squeeze(mean(S2c1(:,ikeep),2)); % 1 drop mean spectrum, data 2, condition 1 0088 Cm12c1(:,i)=abs(Sm12c1./sqrt(Sm1c1.*Sm2c1)); % 1 drop coherence, condition 1 0089 z1i(:,i)=atanh(Cm12c1(:,i))-1/(dof1-4); % 1 drop, bias-corrected Fisher z, condition 1 0090 dz1i(:,i)=(z1i(:,i)-z2)/sqrt(1/(dof1-4)+1/(dof2-2)); % 1 drop, z statistic, condition 1 0091 ps1(:,i)=m1*dz-(m1-1)*dz1i(:,i); 0092 % ps1(:,i)=dof1*dz-(dof1-2)*dz1i(:,i); 0093 end; 0094 ps1m=mean(ps1,2); 0095 for j=1:m2; 0096 jkeep=setdiff(samples2,j); % all samples except j 0097 Sm12c2=squeeze(mean(S12c2(:,jkeep),2)); % 1 drop mean cross-spectrum, condition 2 0098 Sm1c2=squeeze(mean(S1c2(:,jkeep),2)); % 1 drop mean spectrum, data 1, condition 2 0099 Sm2c2=squeeze(mean(S2c2(:,jkeep),2)); % 1 drop mean spectrum, data 2, condition 2 0100 Cm12c2(:,j)=abs(Sm12c2./sqrt(Sm1c2.*Sm2c2)); % 1 drop coherence, condition 2 0101 z2j(:,j)=atanh(Cm12c2(:,j))-1/(dof2-4); % 1 drop, bias-corrected Fisher z, condition 2 0102 dz2j(:,j)=(z1-z2j(:,j))/sqrt(1/(dof1-2)+1/(dof2-4)); % 1 drop, z statistic, condition 2 0103 ps2(:,j)=m2*dz-(m2-1)*dz2j(:,j); 0104 % ps2(:,j)=dof2*dz-(dof2-2)*dz2j(:,j); 0105 end; 0106 0107 % 0108 % Leave one out, both samples 0109 % and pseudo values 0110 % for i=1:m1; 0111 % for j=1:m2; 0112 % dzij(:,i,j)=(z1i(:,i)-z2j(:,j))/sqrt(1/(dof1-4)+1/(dof2-4)); 0113 % dzpseudoval(:,i,j)=m1*m2*dz-(m1-1)*m2*dz1i(:,i)-m1*(m2-1)*dz2j(:,j)+(m1-1)*(m2-1)*dzij(:,i,j); 0114 % % dzpseudoval(:,i,j)=dof1*dof2*dz-(dof1-2)*dof2*dz1i(:,i)-dof1*(dof2-2)*dz2j(:,j)+(dof1-2)*(dof2-2)*dzij(:,i,j); 0115 % end; 0116 % end; 0117 % dzah=sum(sum(dzpseudoval,3),2)/(m1*m2); 0118 ps2m=mean(ps2,2); 0119 % dzar=(sum(ps1,2)+sum(ps2,2))/(m1+m2); 0120 vdz=sum((ps1-ps1m(:,ones(1,m1))).*(ps1-ps1m(:,ones(1,m1))),2)/(m1*(m1-1))+sum((ps2-ps2m(:,ones(1,m2))).*(ps2-ps2m(:,ones(1,m2))),2)/(m2*(m2-1)); 0121 % vdzah=sum(sum((dzpseudoval-dzah(:,ones(1,m1),ones(1,m2))).*(dzpseudoval-dzah(:,ones(1,m1),ones(1,m2))),3),2)/(m1*m2); 0122 % 0123 % Test whether H0 is accepted at the specified p value 0124 % 0125 Adz=zeros(size(dz)); 0126 x=norminv([p/2 1-p/2],0,1); 0127 indx=find(dz>=x(1) & dz<=x(2)); 0128 Adz(indx)=1; 0129 0130 % Adzar=zeros(size(dzar)); 0131 % indx=find(dzar>=x(1) & dzar<=x(2)); 0132 % Adzar(indx)=1; 0133 % 0134 % Adzah=zeros(size(dzah)); 0135 % indx=find(dzah>=x(1) & dzah<=x(2)); 0136 % Adzah(indx)=1;