two_group_test_coherence

 function [dz,vdz,Adz]=two_group_test_coherence(J1c1,J2c1,J1c2,J2c2,p,plt,f)
 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)
 plt    'y' for plot and 'n' for no plot
 f      frequencies (useful for plotting)


 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,plt,f)
0002 % function [dz,vdz,Adz]=two_group_test_coherence(J1c1,J2c1,J1c2,J2c2,p,plt,f)
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 % plt    'y' for plot and 'n' for no plot
0016 % f      frequencies (useful for plotting)
0017 %
0018 %
0019 % Dimensions: J1c1,J2c2: frequencies x number of samples in condition 1
0020 %              J1c2,J2c2: frequencies x number of samples in condition 2
0021 %              number of samples = number of trials x number of tapers
0022 % Outputs:
0023 % dz    test statistic (will be distributed as N(0,1) under H0
0024 % vdz   Arvesen estimate of the variance of dz
0025 % Adz   1/0 for accept/reject null hypothesis of equal population
0026 %       coherences based dz ~ N(0,1)
0027 %
0028 % Note: all outputs are functions of frequency
0029 %
0030 % References: Arvesen, Jackkknifing U-statistics, Annals of Mathematical
0031 % Statisitics, vol 40, no. 6, pg 2076-2100 (1969)
0032 
0033 
0034 if nargin < 4; error('Need four sets of Fourier transforms'); end;
0035 if nargin < 6 || isempty(plt); plt='n'; end;
0036 
0037 %
0038 % Test for matching dimensionalities
0039 %
0040 if size(J1c1)~=size(J2c1) | size(J1c2)~=size(J2c2) | size(J1c1,1)~=size(J1c2,1);
0041     error('Need matching dimensionalities for the Fourier transforms: Check the help file for correct dimensionalities');
0042 else;
0043     m1=size(J1c1,2); % number of samples, condition 1
0044     m2=size(J1c2,2); % number of samples, condition 2
0045     dof1=2*m1; % number of degrees of freedom in the first condition estimates
0046     dof2=2*m2; % number of degrees of freedom in the second condition estimates
0047 end;
0048 if nargin < 7 || isempty(f); f=size(J1c1,1); end;
0049 if nargin < 5 || isempty(p); p=0.05; end; % set the default p value
0050 
0051 %
0052 % Compute the individual condition spectra, coherences
0053 %
0054 S12c1=conj(J1c1).*J2c1; % individual sample cross-spectrum, condition 1
0055 S12c2=conj(J1c2).*J2c2; % individual sample cross-spectrum, condition 2
0056 S1c1=conj(J1c1).*J1c1; % individual sample spectrum, data 1, condition 1
0057 S2c1=conj(J2c1).*J2c1; % individual sample spectrum, data 2, condition 1
0058 S1c2=conj(J1c2).*J1c2; % individual sample spectrum, data 1, condition 2
0059 S2c2=conj(J2c2).*J2c2; % individual sample spectrum, data 2, condition 2
0060 
0061 Sm12c1=squeeze(mean(S12c1,2)); % mean cross spectrum, condition 1
0062 Sm12c2=squeeze(mean(S12c2,2)); % mean cross spectrum, condition 2
0063 Sm1c1=squeeze(mean(S1c1,2)); % mean spectrum, data 1, condition 1
0064 Sm2c1=squeeze(mean(S2c1,2)); % mean spectrum, data 2, condition 1
0065 Sm1c2=squeeze(mean(S1c2,2)); % mean spectrum, data 1, condition 1
0066 Sm2c2=squeeze(mean(S2c2,2)); % mean spectrum, data 2, condition 1
0067 
0068 Cm12c1=abs(Sm12c1./sqrt(Sm1c1.*Sm2c1)); % mean coherence, condition 1
0069 Cm12c2=abs(Sm12c2./sqrt(Sm1c2.*Sm2c2)); % mean coherence, condition 2
0070 
0071 Ccm12c1=Cm12c1; % mean coherence saved for output
0072 Ccm12c2=Cm12c2; % mean coherence saved for output
0073 %
0074 % Compute the statistic dz, and the probability of observing the value dz
0075 % given an N(0,1) distribution i.e. under the null hypothesis
0076 %
0077 z1=atanh(Cm12c1)-1/(dof1-2); % Bias-corrected Fisher z, condition 1
0078 z2=atanh(Cm12c2)-1/(dof2-2); % Bias-corrected Fisher z, condition 2
0079 dz=(z1-z2)/sqrt(1/(dof1-2)+1/(dof2-2)); % z statistic
0080 %
0081 % The remaining portion of the program computes Jackknife estimates of the mean (mdz) and variance (vdz) of dz
0082 %
0083 samples1=[1:m1];
0084 samples2=[1:m2];
0085 %
0086 % Leave one out of one sample
0087 %
0088 for i=1:m1;
0089     ikeep=setdiff(samples1,i); % all samples except i
0090     Sm12c1=squeeze(mean(S12c1(:,ikeep),2)); % 1 drop mean cross-spectrum, condition 1
0091     Sm1c1=squeeze(mean(S1c1(:,ikeep),2)); % 1 drop mean spectrum, data 1, condition 1
0092     Sm2c1=squeeze(mean(S2c1(:,ikeep),2)); % 1 drop mean spectrum, data 2, condition 1
0093     Cm12c1(:,i)=abs(Sm12c1./sqrt(Sm1c1.*Sm2c1)); % 1 drop coherence, condition 1
0094     z1i(:,i)=atanh(Cm12c1(:,i))-1/(dof1-4); % 1 drop, bias-corrected Fisher z, condition 1
0095     dz1i(:,i)=(z1i(:,i)-z2)/sqrt(1/(dof1-4)+1/(dof2-2)); % 1 drop, z statistic, condition 1
0096     ps1(:,i)=m1*dz-(m1-1)*dz1i(:,i);
0097 %      ps1(:,i)=dof1*dz-(dof1-2)*dz1i(:,i);
0098 end; 
0099 ps1m=mean(ps1,2);
0100 for j=1:m2;
0101     jkeep=setdiff(samples2,j); % all samples except j
0102     Sm12c2=squeeze(mean(S12c2(:,jkeep),2)); % 1 drop mean cross-spectrum, condition 2
0103     Sm1c2=squeeze(mean(S1c2(:,jkeep),2)); % 1 drop mean spectrum, data 1, condition 2
0104     Sm2c2=squeeze(mean(S2c2(:,jkeep),2)); % 1 drop mean spectrum, data 2, condition 2
0105     Cm12c2(:,j)=abs(Sm12c2./sqrt(Sm1c2.*Sm2c2)); % 1 drop coherence, condition 2
0106     z2j(:,j)=atanh(Cm12c2(:,j))-1/(dof2-4); % 1 drop, bias-corrected Fisher z, condition 2
0107     dz2j(:,j)=(z1-z2j(:,j))/sqrt(1/(dof1-2)+1/(dof2-4)); % 1 drop, z statistic, condition 2
0108     ps2(:,j)=m2*dz-(m2-1)*dz2j(:,j);
0109 %      ps2(:,j)=dof2*dz-(dof2-2)*dz2j(:,j);
0110 end;
0111 
0112 %
0113 % Leave one out, both samples
0114 % and pseudo values
0115 % for i=1:m1;
0116 %     for j=1:m2;
0117 %         dzij(:,i,j)=(z1i(:,i)-z2j(:,j))/sqrt(1/(dof1-4)+1/(dof2-4));
0118 %         dzpseudoval(:,i,j)=m1*m2*dz-(m1-1)*m2*dz1i(:,i)-m1*(m2-1)*dz2j(:,j)+(m1-1)*(m2-1)*dzij(:,i,j);
0119 % %         dzpseudoval(:,i,j)=dof1*dof2*dz-(dof1-2)*dof2*dz1i(:,i)-dof1*(dof2-2)*dz2j(:,j)+(dof1-2)*(dof2-2)*dzij(:,i,j);
0120 %     end;
0121 % end;
0122 % dzah=sum(sum(dzpseudoval,3),2)/(m1*m2);
0123 ps2m=mean(ps2,2);
0124 % dzar=(sum(ps1,2)+sum(ps2,2))/(m1+m2);
0125 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));
0126 % vdzah=sum(sum((dzpseudoval-dzah(:,ones(1,m1),ones(1,m2))).*(dzpseudoval-dzah(:,ones(1,m1),ones(1,m2))),3),2)/(m1*m2);
0127 %
0128 % Test whether H0 is accepted at the specified p value
0129 %
0130 Adz=zeros(size(dz));
0131 x=norminv([p/2 1-p/2],0,1);
0132 indx=find(dz>=x(1) & dz<=x(2)); 
0133 Adz(indx)=1;
0134 
0135 if strcmp(plt,'y');
0136     if isempty(f) || nargin < 6;
0137         f=linspace(0,1,length(dz));
0138     end;
0139     %
0140     % Compute the coherences
0141     %
0142     S121=mean(conj(J1c1).*J2c1,2);
0143     S122=mean(conj(J1c2).*J2c2,2);
0144     S111=mean(conj(J1c1).*J1c1,2);
0145     S221=mean(conj(J2c1).*J2c1,2);
0146     S112=mean(conj(J1c2).*J1c2,2);
0147     S222=mean(conj(J2c2).*J2c2,2);
0148     C121=abs(S121)./sqrt(S111.*S221);
0149     C122=abs(S122)./sqrt(S112.*S222);
0150     %
0151     % Plot the coherence
0152     %
0153     subplot(311); 
0154     plot(f,C121,f,C122); legend('Data 1','Data 2');
0155     set(gca,'FontName','Times New Roman','Fontsize', 16);
0156     ylabel('Coherence');
0157     title('Two group test for coherence');
0158     subplot(312);
0159     plot(f,dz);
0160     set(gca,'FontName','Times New Roman','Fontsize', 16);
0161     ylabel('Test statistic');
0162     conf=norminv(1-p/2,0,1);
0163     line(get(gca,'xlim'),[conf conf]);
0164     line(get(gca,'xlim'),[-conf -conf]);
0165     subplot(313);
0166     plot(f,vdz);
0167     set(gca,'FontName','Times New Roman','Fontsize', 16);
0168     xlabel('frequency'); ylabel('Jackknifed variance');
0169 end;
0170 % Adzar=zeros(size(dzar));
0171 % indx=find(dzar>=x(1) & dzar<=x(2));
0172 % Adzar(indx)=1;
0173 %
0174 % Adzah=zeros(size(dzah));
0175 % indx=find(dzah>=x(1) & dzah<=x(2));
0176 % Adzah(indx)=1;