


function [dz,vdz,Adz]=two_group_test_spectrum(J1,J2,p)
Test the null hypothesis (H0) that data sets J1, J2 in
two conditions c1,c2 have equal population spectrum
Usage:
[dz,vdz,Adz]=two_sample_test_spectrum(J1,J2,p)
Inputs:
J1 tapered fourier transform in condition 1
J2 tapered fourier transform 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: J1: frequencies x number of samples in condition 1
J2: 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_spectrum(J1,J2,p,plt,f) 0002 % function [dz,vdz,Adz]=two_group_test_spectrum(J1,J2,p) 0003 % Test the null hypothesis (H0) that data sets J1, J2 in 0004 % two conditions c1,c2 have equal population spectrum 0005 % 0006 % Usage: 0007 % [dz,vdz,Adz]=two_sample_test_spectrum(J1,J2,p) 0008 % 0009 % Inputs: 0010 % J1 tapered fourier transform in condition 1 0011 % J2 tapered fourier transform in condition 2 0012 % p p value for test (default: 0.05) 0013 % plt 'y' for plot and 'n' for no plot 0014 % f frequencies (useful for plotting) 0015 % 0016 % 0017 % Dimensions: J1: frequencies x number of samples in condition 1 0018 % J2: 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 % 0027 % Note: all outputs are functions of frequency 0028 % 0029 % References: Arvesen, Jackkknifing U-statistics, Annals of Mathematical 0030 % Statisitics, vol 40, no. 6, pg 2076-2100 (1969) 0031 0032 if nargin < 2; error('Need four sets of Fourier transforms'); end; 0033 if nargin < 4 || isempty(plt); plt='n'; end; 0034 % 0035 % Test for matching dimensionalities 0036 % 0037 m1=size(J1,2); % number of samples, condition 1 0038 m2=size(J2,2); % number of samples, condition 2 0039 dof1=m1; % degrees of freedom, condition 1 0040 dof2=m2; % degrees of freedom, condition 2 0041 if nargin < 5 || isempty(f); f=size(J1c1,1); end; 0042 if nargin < 3; p=0.05; end; % set the default p value 0043 0044 % 0045 % Compute the individual condition spectra, coherences 0046 % 0047 S1=conj(J1).*J1; % spectrum, condition 1 0048 S2=conj(J2).*J2; % spectrum, condition 2 0049 0050 Sm1=squeeze(mean(S1,2)); % mean spectrum, condition 1 0051 Sm2=squeeze(mean(S2,2)); % mean spectrum, condition 2 0052 % 0053 % Compute the statistic dz, and the probability of observing the value dz 0054 % given an N(0,1) distribution i.e. under the null hypothesis 0055 % 0056 bias1=psi(dof1)-log(dof1); bias2=psi(dof2)-log(dof2); % bias from Thomson & Chave 0057 var1=psi(1,dof1); var2=psi(1,dof2); % variance from Thomson & Chave 0058 z1=log(Sm1)-bias1; % Bias-corrected Fisher z, condition 1 0059 z2=log(Sm2)-bias2; % Bias-corrected Fisher z, condition 2 0060 dz=(z1-z2)/sqrt(var1+var2); % z statistic 0061 pdz=normpdf(dz,0,1); % probability of observing value dz 0062 % 0063 % The remaining portion of the program computes Jackknife estimates of the mean (mdz) and variance (vdz) of dz 0064 % 0065 samples1=[1:m1]; 0066 samples2=[1:m2]; 0067 % 0068 % Leave one out of one sample 0069 % 0070 bias11=psi(dof1-1)-log(dof1-1); var11=psi(1,dof1-1); 0071 for i=1:m1; 0072 ikeep=setdiff(samples1,i); % all samples except i 0073 Sm1=squeeze(mean(S1(:,ikeep),2)); % 1 drop mean spectrum, data 1, condition 1 0074 z1i(:,i)=log(Sm1)-bias11; % 1 drop, bias-corrected Fisher z, condition 1 0075 dz1i(:,i)=(z1i(:,i)-z2)/sqrt(var11+var2); % 1 drop, z statistic, condition 1 0076 ps1(:,i)=m1*dz-(m1-1)*dz1i(:,i); 0077 end; 0078 ps1m=mean(ps1,2); 0079 bias21=psi(dof2-1)-log(dof2-1); var21=psi(1,dof2-1); 0080 for j=1:m2; 0081 jkeep=setdiff(samples2,j); % all samples except j 0082 Sm2=squeeze(mean(S2(:,jkeep),2)); % 1 drop mean spectrum, data 2, condition 2 0083 z2j(:,j)=log(Sm2)-bias21; % 1 drop, bias-corrected Fisher z, condition 2 0084 dz2j(:,j)=(z1-z2j(:,j))/sqrt(var1+var21); % 1 drop, z statistic, condition 2 0085 ps2(:,j)=m2*dz-(m2-1)*dz2j(:,j); 0086 end; 0087 % 0088 % Leave one out, both samples 0089 % and pseudo values 0090 % for i=1:m1; 0091 % for j=1:m2; 0092 % dzij(:,i,j)=(z1i(:,i)-z2j(:,j))/sqrt(var11+var21); 0093 % dzpseudoval(:,i,j)=m1*m2*dz-(m1-1)*m2*dz1i(:,i)-m1*(m2-1)*dz2j(:,j)+(m1-1)*(m2-1)*dzij(:,i,j); 0094 % end; 0095 % end; 0096 % 0097 % Jackknife mean and variance 0098 % 0099 % dzah=sum(sum(dzpseudoval,3),2)/(m1*m2); 0100 ps2m=mean(ps2,2); 0101 % dzar=(sum(ps1,2)+sum(ps2,2))/(m1+m2); 0102 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)); 0103 % vdzah=sum(sum((dzpseudoval-dzah(:,ones(1,m1),ones(1,m2))).*(dzpseudoval-dzah(:,ones(1,m1),ones(1,m2))),3),2)/(m1*m2); 0104 % 0105 % Test whether H0 is accepted at the specified p value 0106 % 0107 Adz=zeros(size(dz)); 0108 x=norminv([p/2 1-p/2],0,1); 0109 indx=find(dz>=x(1) & dz<=x(2)); 0110 Adz(indx)=1; 0111 0112 % Adzar=zeros(size(dzar)); 0113 % indx=find(dzar>=x(1) & dzar<=x(2)); 0114 % Adzar(indx)=1; 0115 % 0116 % Adzah=zeros(size(dzah)); 0117 % indx=find(dzah>=x(1) & dzah<=x(2)); 0118 % Adzah(indx)=1; 0119 if strcmp(plt,'y'); 0120 if isempty(f) || nargin < 6; 0121 f=linspace(0,1,length(dz)); 0122 end; 0123 % 0124 % Compute the coherences 0125 % 0126 S1=squeeze(mean(conj(J1).*J1,2)); 0127 S2=squeeze(mean(conj(J2).*J2,2)); 0128 % 0129 % Plot the coherence 0130 % 0131 subplot(311); 0132 plot(f,S1,f,S2); legend('Data 1','Data 2'); 0133 set(gca,'FontName','Times New Roman','Fontsize', 16); 0134 ylabel('Spectra'); 0135 title('Two group test for coherence'); 0136 subplot(312); 0137 plot(f,dz); 0138 set(gca,'FontName','Times New Roman','Fontsize', 16); 0139 ylabel('Test statistic'); 0140 conf=norminv(1-p/2,0,1); 0141 line(get(gca,'xlim'),[conf conf]); 0142 line(get(gca,'xlim'),[-conf -conf]); 0143 subplot(313); 0144 plot(f,vdz); 0145 set(gca,'FontName','Times New Roman','Fontsize', 16); 0146 xlabel('frequency'); ylabel('Jackknifed variance'); 0147 end;