two_group_test_coherence

PURPOSE ^

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

SYNOPSIS ^

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

DESCRIPTION ^

 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)

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)
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;

Generated on Tue 28-Mar-2006 14:37:41 by m2html © 2003