Home > SurfStat > SurfStatP.m

SurfStatP

PURPOSE ^

Corrected P-values for vertices and clusters.

SYNOPSIS ^

function [ pval, peak, clus, clusid ] = SurfStatP( slm, mask, clusthresh );

DESCRIPTION ^

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.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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

Generated on Fri 26-Sep-2008 14:05:29 by m2html © 2003