% The first version is done by Wenlin % \$Id: spmd_stat.m,v 1.9 2008/02/20 12:19:02 nichols Exp \$ function [stat,p,coef] = spmd_stat(action,res,X,para,df) % calculate % 1. Dubin Watson statistic % 2. Cook Weisberg score statistic % 3. Shapiro-Wilk satatistic % 4. Proportion of outlier % 5. Cumulative periodogram statistic % ______________________________________________________________________ % FORMAT spmd_stat(action,res,x,para) % action - option % res - residual data % X - design matrix % Para - parameters % df - degree of freedom of error % ____________________ Function called _________________________________ % % spm_detrend % spm_Xcdf % spm_invNcdf % spm_Ncdf % spm_Icdf % ______________________________________________________________________ % @(#)spmd_stat.m 1.4 Wen-Lin Luo 04/06/16 %load SPM %df = SPM.xX.erdf; [nScan nVox] = size(res); if strmatch(action,'sw'), else if size(X,1) ~= nScan, error(['The length of the residuals and design matrix must be the' ... ' same!']); end end ResSS = sum(res.^2); switch lower(action) case 'dw' %---------------------------------------------------------------------------- % Durbin Watson statistics %---------------------------------------------------------------------------- if isempty(para) %-------------------------------------------------------------------------- % Calculate the approximate distribution parameters of D %-------------------------------------------------------------------------- nDWBeta = size(X,2); DWD1X = diff(X); DWD2X = diff(DWD1X); DWDP = (2*nScan-1)-trace(DWD1X'*DWD1X*(inv(X'*X))); DWDQ = 2*(3*nScan-4) - 2*trace(DWD2X'*DWD2X*(inv(X'*X))) + ... trace((DWD1X'*DWD1X*(inv(X'*X)))^2); DWED = DWDP/(nScan-(nDWBeta)-1); DWVarD= 2/((nScan-(nDWBeta)-1)*(nScan-(nDWBeta)+1))*(DWDQ-DWDP*DWED); DWp = (DWED*(4-DWED)/DWVarD-1)*DWED/4; DWq = (DWED*(4-DWED)/DWVarD)-1-DWp; coef = [DWp;DWq]; else coef = para; end %----------------------------------------------------------------- %-Durbin-Watson statistic %----------------------------------------------------------------- DW = sum((diff(res)).^2)./ResSS; %-Durbin-Watson statistic stat = coef(1)*(4-DW)./(coef(2)*DW); p = 1-betainc(1-DW/4,coef(2), coef(1)); %-DW probability p = p+(p==0)*eps; p = -log10(p); stat = {stat}; case 'score' %----------------------------------------------------------------------- % Score test / Cook Weisberg score statistic %------------------------------------------------------------------------- SSr = nScan*(res.^2)*(spdiags((ResSS.^(-1))',0,nVox,nVox)); meanSSr = mean(SSr); [nrowSSr ncolSSr] = size(meanSSr); if ncolSSr > 1 SYY = sum((SSr-repmat(meanSSr,nScan,1)).^2); %-scaled squared residual else SYY = sum((SSr-repmat(meanSSr(1,1),nScan,1)).^2); %-scaled squared residual end coef = para; %----------------------------------------------------------------- % regress the scaled squared residual on possible regressors %----------------------------------------------------------------- if isempty(coef), Rres = SSr-X*pinv(X)*SSr; %def = size(X,2)-1; def = max(1,size(X,2)-1); %add Tom's change elseif coef == 1, X = spm_detrend(X); SSr = spm_detrend(SSr); sumX2 = sum(X.^2)'; % define [xindex, yindex] = size(sumX2); for i=1:xindex for j=1:yindex if abs(sumX2(i,j)-0)<1e-7, sumX2(i,j)=0.000001; end end end %Rres = SSr - X*spdiags(sum(X.*SSr)'./sum(X.^2)', ... %0,nVox,nVox); Rres = SSr - X*spdiags(sum(X.*SSr)'./sumX2,0,nVox,nVox); def = 1; end stat = (SYY-sum(Rres.^2))./2; %-Score statistics p = 1-spm_Xcdf(stat,def); %-Score probability p = p+(p==0)*eps; p = -log10(p); stat = {stat}; case 'sw' %------------------------------------------------------------------------- %-For Shapiro-Wilks Normality Test %------------------------------------------------------------------------- if isempty(para) %----------------------------------------------------------------- % Calculate necessary coefficient for Shapiro-Wilks test %----------------------------------------------------------------- % Calculate expected normal order statistics %----------------------------------------------------------------- eprob = [0.5./nScan:1./nScan:(nScan-0.5)./nScan]; enorm = spm_invNcdf(eprob,0,1); % Calculate coefficient for W statistic %----------------------------------------------------------------- one = [0 ones(1,nScan-2) 0]; SWa = 2*enorm; SWa1s = (6.*nScan+7)./(6.*nScan+13).*... ((exp(1)./(nScan+2).*((nScan+1)./(nScan+2)).^(nScan-2)).^0.5); SWa(:,1) = -((SWa1s./(1-2.*SWa1s))*(one*(SWa.^2)')).^0.5; SWa(:,nScan) = ((SWa1s./(1-2.*SWa1s))*(one*(SWa.^2)')).^0.5; SWa = SWa./(SWa*SWa').^0.5; % % Calculate some constants for Approximation for Significance level % of W: lambda, mean and SD for the transformation of W % poly = log(nScan)-5; % for sample size greater than 20, d=5 %----------------------------------------------------------------- SWcoeff1 = [0.002989646 0.00879701 -0.0241665 0 0.318828 0.480385]; SWcoeff2 = [-0.01504614 -0.03513666 0.1066339 -0.04183209... -1.37888 -1.91487]; SWcoeff3 = [0.003852646 -0.03215018 -0.01638782 0.1773538... -0.331885 -1.015807 -3.73538]; SWlambda = polyval(SWcoeff1,log(nScan)-5); SWmu = exp(polyval(SWcoeff2,log(nScan)-5)); SWsd = exp(polyval(SWcoeff3,log(nScan)-5)); coef = {SWa SWlambda SWmu SWsd}; else coef = para; end %----------------------------------------------------------------- %-Shapiro-Wilks statistics %----------------------------------------------------------------- sRes = sort(res); WRes = (coef{1}*sRes).^2./sum((res-repmat(mean(res),nScan,1)).^2); TWRes = (1-WRes).^coef{2}; %-approximate normalizing % transformation of W stat = (TWRes-coef{3})./coef{4}; %-standard normal deviate for % transformation of W p = 1-spm_Ncdf(stat,0,1); %-Shapiro-Wilks probability p = p+(p==0)*eps; p = -log10(p); stat = {stat}; case 'outl' %----------------------------------------------------------------- %-outlier count %----------------------------------------------------------------- if isempty(para), coef = 3; else coef = para; end MSE = ResSS./df; %-Mean squared error Outh = X*pinv(X); %-leverage Outh = diag(Outh); % I_null= (1-Outh) coef); %-index of outliers %-spatial outlier count %--------------------------------------------------------------- p = 1-spm_Icdf(sum(Iout),nScan,2*(1-spm_Ncdf(coef))); p = p+(p==0)*eps; p = -log10(p); %-temporal outlier count stat = {sum(Iout) sum(Iout,2)}; case 'cp' %----------------------------------------------------------------- %-cumulative periodogram %----------------------------------------------------------------- nVar = rank(X); if isempty(para), %-------------------------------------------------------------- %-Calculate the vector for uncorrelating the residuals %-------------------------------------------------------------- if nVar == 0, coef = 0; else Xb0 = X(1:nVar,:); Xb1 = X(nVar+1:nScan,:); [V,D] = eig(Xb0*pinv(X'*X)*Xb0'); D = spdiags(D,0); V = V./(repmat(sqrt(sum(V.*V)),nVar,1)); H = sum(D==1); blus1 = zeros(nVar,nVar); for i = H+1:nVar blus1 = blus1 + D(i)/(1+D(i))*V(:,i)*V(:,i)'; end coef = Xb1*pinv(Xb0)*blus1; end else coef = para; end %-Cumulative periodogram %----------------------------------------------------------------- power = zeros(floor((nScan-nVar-1)/2),nVox); for i = 1:nVox %- Calculate BLUS residual vector %--------------------------------------------------------------- Rblus = res(nVar+1:nScan,i)-coef*res(1:nVar,i); %- Calculate periodogram %--------------------------------------------------------------- Spec = abs(fft(Rblus,nScan-nVar)).^2; % Select first half if rem(nScan-nVar,2), % nfft odd select = (2:(nScan-nVar+1)/2)'; %-don't include DC or Nyquist % components else select = (2:(nScan-nVar)/2)'; end freq = (select-1)*2/(nScan-nVar); % Calculate the single-sided spectrum which includes the full power power(:,i) = [2*Spec(select)]; end %-Calculate cumulative periodogram CPseries = cumsum(power); Flength = size(CPseries,1); CPS = CPseries./(repmat(sum(power),Flength,1)); %-setup the upper and lower statistic for the Cumulative periodogram MaxA1 = abs(CPS(1:Flength-1,:)-... (repmat([1:Flength-1]',1,nVox))/(Flength)); MaxB1 = abs(CPS(1:Flength-1,:)-... (repmat([1:Flength-1]',1,nVox)-1)/(Flength)); %-Calculate the Cumulative periodogram test statistic stat = max(max(MaxA1,MaxB1)); %-Calculate the coresponding p-value for the CPD CPZ = stat*(sqrt(Flength-1)+0.12+0.11/sqrt(Flength-1)); p = 1-sqrt(2*pi)*(sum(exp(-((2*repmat([1:30]',1,nVox)... -1)*pi).^2./(8*repmat(CPZ,size([1:30],2),1).^2))))./CPZ; p = p+(p==0)*eps; p = -log10(p); %-Calculate sum of periodogram stat = {stat sum(power,2) freq}; otherwise, warning('Unknown action string'); end;