Home > chronux > spectral_analysis > statistical_tests > two_group_test_coherence.m

two_group_test_coherence

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:

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)
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:
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 %
0131 x=norminv([p/2 1-p/2],0,1);
0132 indx=find(dz>=x(1) & dz<=x(2));
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;