Corrected P-values for vertices and clusters. Usage: [ pval, peak, clus, clusid ] = SurfStatP( slm [, mask, [, clusthresh] ] ); slm.t = l x v matrix of test statistics, v=#vertices; the first row slm.t(1,:) is the test statistic, and the other rows are used to calculate cluster resels if k>1. See SurfStatF for the precise definition of the extra rows. slm.df = degrees of freedom. slm.dfs = 1 x v vector of optional effective degrees of freedom. slm.k = k=#variates. slm.resl = e x k matrix of sum over observations of squares of differences of normalized residuals along each edge. slm.tri = 3 x t matrix of triangle indices, 1-based, t=#triangles. or slm.lat = nx x ny x nz matrix, 1=in, 0=out, [nx,ny,nz]=size(volume). mask = 1 x v logical vector, 1=inside, 0=outside, v=#vertices, = ones(1,v), i.e. the whole surface, by default. clusthresh = P-value threshold or statistic threshold for defining clusters, 0.001 by default. pval.P = 1 x v vector of corrected P-values for vertices. pval.C = 1 x v vector of corrected P-values for clusters. pval.mask = copy of mask. peak.t = np x 1 vector of peaks (local maxima). peak.vertid = np x 1 vector of vertex id's (1-based). peak.clusid = np x 1 vector of cluster id's that contain the peak. peak.P = np x 1 vector of corrected P-values for the peak. clus.clusid = nc x 1 vector of cluster id numbers clus.nverts = nc x 1 vector of number of vertices in the cluster. clus.resels = nc x 1 vector of resels in the cluster. clus.P = nc x 1 vector of corrected P-values for the cluster. clusid = 1 x v vector of cluster id's for each vertex. Reference: Worsley, K.J., Andermann, M., Koulis, T., MacDonald, D. & Evans, A.C. (1999). Detecting changes in nonisotropic images. Human Brain Mapping, 8:98-101.
0001 function [ pval, peak, clus, clusid ] = SurfStatP( slm, mask, clusthresh ); 0002 0003 %Corrected P-values for vertices and clusters. 0004 % 0005 % Usage: [ pval, peak, clus, clusid ] = 0006 % SurfStatP( slm [, mask, [, clusthresh] ] ); 0007 % 0008 % slm.t = l x v matrix of test statistics, v=#vertices; the first 0009 % row slm.t(1,:) is the test statistic, and the other 0010 % rows are used to calculate cluster resels if k>1. See 0011 % SurfStatF for the precise definition of the extra rows. 0012 % slm.df = degrees of freedom. 0013 % slm.dfs = 1 x v vector of optional effective degrees of freedom. 0014 % slm.k = k=#variates. 0015 % slm.resl = e x k matrix of sum over observations of squares of 0016 % differences of normalized residuals along each edge. 0017 % slm.tri = 3 x t matrix of triangle indices, 1-based, t=#triangles. 0018 % or 0019 % slm.lat = nx x ny x nz matrix, 1=in, 0=out, [nx,ny,nz]=size(volume). 0020 % mask = 1 x v logical vector, 1=inside, 0=outside, v=#vertices, 0021 % = ones(1,v), i.e. the whole surface, by default. 0022 % clusthresh = P-value threshold or statistic threshold for 0023 % defining clusters, 0.001 by default. 0024 % 0025 % pval.P = 1 x v vector of corrected P-values for vertices. 0026 % pval.C = 1 x v vector of corrected P-values for clusters. 0027 % pval.mask = copy of mask. 0028 % peak.t = np x 1 vector of peaks (local maxima). 0029 % peak.vertid = np x 1 vector of vertex id's (1-based). 0030 % peak.clusid = np x 1 vector of cluster id's that contain the peak. 0031 % peak.P = np x 1 vector of corrected P-values for the peak. 0032 % clus.clusid = nc x 1 vector of cluster id numbers 0033 % clus.nverts = nc x 1 vector of number of vertices in the cluster. 0034 % clus.resels = nc x 1 vector of resels in the cluster. 0035 % clus.P = nc x 1 vector of corrected P-values for the cluster. 0036 % clusid = 1 x v vector of cluster id's for each vertex. 0037 % 0038 % Reference: Worsley, K.J., Andermann, M., Koulis, T., MacDonald, D. 0039 % & Evans, A.C. (1999). Detecting changes in nonisotropic images. 0040 % Human Brain Mapping, 8:98-101. 0041 0042 [l,v]=size(slm.t); 0043 if nargin<2 | isempty(mask) 0044 mask=logical(ones(1,v)); 0045 end 0046 if nargin<3 0047 clusthresh=0.001; 0048 end 0049 0050 df=zeros(2); 0051 ndf=length(slm.df); 0052 df(1,1:ndf)=slm.df; 0053 df(2,1:2)=slm.df(ndf); 0054 if isfield(slm,'dfs') 0055 df(1,ndf)=mean(slm.dfs(mask>0)); 0056 end 0057 0058 if v==1 0059 pval.P=stat_threshold(0,1,0,df,[10 slm.t(1)],[],[],[],slm.k,[],[],0); 0060 pval.P=pval.P(2); 0061 peak=[]; 0062 clus=[]; 0063 clusid=[]; 0064 return 0065 end 0066 0067 if clusthresh<1 0068 thresh=stat_threshold(0,1,0,df,clusthresh,[],[],[],slm.k,[],[],0); 0069 else 0070 thresh=clusthresh; 0071 end 0072 0073 [resels,reselspvert,edg]=SurfStatResels(slm,mask); 0074 0075 if max(slm.t(1,mask))<thresh 0076 pval.P=stat_threshold(resels,v,1,df,[10 slm.t],[],[],[],slm.k,[],[],0); 0077 pval.P=pval.P((1:v)+1); 0078 peak=[]; 0079 clus=[]; 0080 clusid=[]; 0081 else 0082 [peak,clus,clusid]=SurfStatPeakClus(slm,mask,thresh,reselspvert,edg); 0083 [pp,clpval]=stat_threshold(resels,v,1,df,... 0084 [10 peak.t' slm.t(1,:)],thresh,[10; clus.resels],[],slm.k,[],[],0); 0085 peak.P=pp((1:length(peak.t))+1)'; 0086 pval.P=pp(length(peak.t)+(1:v)+1); 0087 if slm.k>1 0088 j=(slm.k-1):-2:0; 0089 sphere=zeros(1,slm.k); 0090 sphere(j+1)=exp((j+1)*log(2)+(j/2)*log(pi)+gammaln((slm.k+1)/2)- ... 0091 gammaln(j+1)-gammaln((slm.k+1-j)/2)); 0092 sphere=sphere.*(4*log(2)).^(-(0:(slm.k-1))/2)/ndf; 0093 [pp,clpval]=stat_threshold(conv(resels,sphere),Inf,1,df,... 0094 [],thresh,[10; clus.resels],[],[],[],[],0); 0095 end 0096 clus.P=clpval(2:length(clpval)); 0097 pval.C=interp1([0; clus.clusid],[1; clus.P],clusid); 0098 end 0099 tlim=stat_threshold(resels,v,1,df,[0.5 1],[],[],[],slm.k,[],[],0); 0100 tlim=tlim(2); 0101 pval.P=pval.P.*(slm.t(1,:)>tlim)+(slm.t(1,:)<=tlim); 0102 pval.mask=mask; 0103 0104 return 0105 end 0106