# two_group_test_coherence

## PURPOSE

function [dz,vdz,Adz]=two_group_test_coherence(J1c1,J2c1,J1c2,J2c2,p,plt,f)

## SYNOPSIS

function [dz,vdz,Adz]=two_group_test_coherence(J1c1,J2c1,J1c2,J2c2,p,plt,f)

## DESCRIPTION

     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)

## CROSS-REFERENCE INFORMATION

This function calls:

This function is called by:

## SOURCE CODE

    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;