function FDRill(Timg,df,alpha,FDPplot) % Show two plots illustrating the FDR threshold % FDRill(Timg,df,alpha,FDPplot) % Timg - Filename of T or F statistic image % (If omitted, user is prompted) % df - Degrees of freedom; scalar for T, 2-vector for F statstic % If Timg isn't from SPM, df must be specified. % If df is a string, it can be 'P' or '1-P' to specify that % the image is P-values or one minus P-values, respectively. % alpha - Level at which to control FDR (0.05 default) % FDPplot- If 0, top plot is a histogram (default). % If 1, top plot is a cumulative plot, showing, for each % possible threshold, the observed and expected number of % detections and the estimated false discovery proportion (FDP). % % Alternate usage: % FDRill('FDPplot') or FDRill FDPplot % is equivalent to FDRill('',[],[],1) % % Top plot is a root-o-gram, a histogram where square-root counts (or % frequency) are plotted instead of counts. It has the effect of % magnifying the tails relative to the rest of the distribution. % % The bottom plot is the log-log pp-plot that defines the Benjamini & % Hochberg FDR rule. % % Works with either SPM99, SPM2, SPM5, or SPM8. % %________________________________________________________________________ % \$Id: FDRill.m,v 1.14 2008/05/27 11:58:46 nichols Exp \$ figure % 0 = Not Pvalues (t); 1 = Pvalues; -1 = 1-Pvalues isP=0; if nargin>1 & strcmp(Timg,'FDPplot') Timg = ''; FDPplot = 1; elseif nargin<4 FDPplot = 0; end % Get statistic image if nargin<1 | isempty(Timg) Timg = sf_spm_get(1,'spm*img','Select Statistic Image'); end V = spm_vol(Timg); % Get DF (hopefully from image comment field) if nargin<2 | isempty(df) df = sf_GetDF(V); if isempty(df) df = spm_input('Enter DF (df,[df1 df2], or 0 try SPM.mat)','+1','r'); if (df==0) warning('Assuming T statistic image; manually specify DF if F image'); SPM = sf_spm_get(1,'SPM.mat','Select SPM.mat'); SPM = load(SPM); if isfield(SPM,'SPM') % For SPM2 SPM = SPM.SPM; end df = SPM.xX.erdf; clear SPM end end else if isstr(df) if strcmp(df,'P') isP=1; elseif strcmp(df,'1-P') isP=-1; else error('Unknown DF code!') end df=100; end end if nargin<3 | isempty(alpha) alpha = 0.05; end if nargin<4 & isempty(FDPplot) FDPplot = 0; end % % Load and condition the T statistic image data % T = spm_read_vols(V); T=T(:); T(isnan(T))=[]; T(T==0)=[]; % % Get p-values, FDR threshold % if isP==0 if length(df)>1 Tp = 1-spm_Fcdf(T,df); STAT = 'F'; elseif finite(df) Tp = 1-spm_Tcdf(T,df); STAT = 'T'; elseif isinf(df) Tp = 1-spm_Ncdf(T); STAT = 'N'; end else if isP==1 Tp = T; elseif isP==-1 Tp = (1-T); end STAT = 'N'; T = spm_invNcdf(1-T); end Ts = sort(T); Tps = sort(Tp); iv = (1:length(T))/length(T); Tpt = myFDR(Tp,alpha); if isempty(Tpt) Tpt = NaN; Tt = NaN; else if length(df)>1 Tt = spm_invFcdf(1-Tpt,df); elseif finite(df) Tt = spm_invTcdf(1-Tpt,df); else Tt = spm_invNcdf(1-Tpt); end end V = length(Tps); % % Root-o-gram & PP plots % ax = []; h = []; ax = [subplot(2,1,1) subplot(2,1,2)]; axes(ax(1)); % ax = axes('position',[ 0.13 0.64 0.80 0.28]); [n,x]=hist(T,50);n=n/sum(n)/(x(2)-x(1)); if length(df)>1 NullPDF = spm_Fpdf(x,df); elseif finite(df) NullPDF = spm_Tpdf(x,df); else NullPDF = spm_Npdf(x); end h = [h bar(x,sqrt(n))]; hold on; h = [h plot(x,sqrt(NullPDF),'color','Green','LineWidth',2)]; hold off myabline('v',0) myabline('v',Tt,'LineStyle','-','Color',[0.3 0.3 0.3]); tmp = get(gca,'Ylim'); tmp = tmp(2)-diff(tmp)*0.15; if ~isnan(Tt) text(Tt+.25,tmp,sprintf('%c_{FDR}=%3.3g',STAT,Tt),'color',[0.3 0.3 0.3],'FontSize', 14) end ylabel('Sqrt-Frequency') if STAT=='T' str = sprintf('T_{%d}',df); elseif STAT=='F' str = sprintf('F_{%d,%d}',df); elseif STAT=='N' str = sprintf('N'); end title(['Null & Observed ' str ' dist^{n}s'],'FontSize',14) if FDPplot % % Plot of emperical FDR, the estimated FDP % % First, 'background' the histogram ylabel('') set(ax(1),'Ytick',[]) % Change histogram to light blue and gray % This will surely break in other old ML versions set(h(1),'FaceColor',[0.8 0.8 0.95],'EdgeColor',[0.8 0.8 0.8]) for i = 2:length(h) set(h(i),'Color',[0.8 0.8 0.8]+get(h(i),'color')*.2) end if length(df)>1 NullCDF = spm_Fcdf(Ts,df); elseif finite(df) NullCDF = spm_Tcdf(Ts,df); else NullCDF = spm_Ncdf(Ts); end nNulRej = V*flipud(Tps); nRej = flipud((1:V)'); FDPhat = min(1,nNulRej./nRej); ax2 = axes('position',get(ax(1),'position'),... 'color','none','Xtick',[]); h2 = semilogy(Ts,1-NullCDF,Ts,(V:-1:1)/V,Ts,FDPhat,'Linewidth',2); set(gca,'color','none',... 'Xtick',[],... 'Xlim',get(ax(1),'Xlim')) % set(gca,'Ylim',[0 1]) % Needed if not semilog set(h2(1),'color','Green'); set(h2(2),'color','Blue'); set(h2(3),'color',[1 .5 0]) ylabel('Percent') myabline('h',alpha,'color','red','linewidth',1,'linestyle','--') hl=legend('% Detected, Ho-expected',... '% Detected, actual',... 'False Discovery Proportion'); set(hl,'Color',[1 1 1]); end %ax = [ax axes('position',[ 0.13 0.11 0.80 0.30])]; subplot(2,1,2) h = [h loglog(iv,Tps,'-o')]; axis image set(get(gca,'children'),'LineWidth',2,'MarkerSize',4) myabline(0,1,'color','green','linestyle','-'); myabline(0,alpha,'LineStyle','--','color','red') if ~isnan(Tpt) myabline('h',Tpt,'LineStyle','-','Color',[0.3 0.3 0.3]); end title('Ordered p-value plots - loglog','FontSize',14) text(10^(-4.75),Tpt*4,sprintf('P_{FDR}=%3.3g',Tpt),'color',[0.3 0.3 0.3],'FontSize',14) xlabel('index/V'); ylabel('Ordered p-value') function [pID,pN] = myFDR(p,alpha) % % p - vector of p-values % alpha - False Discovery Rate level % % pID - p-value threshold based on independence or positive dependence % pN - "Nonparametric" (any covariance structure) p-value threshold %______________________________________________________________________________ % Based on FDR.m 1.4 Tom Nichols 02/07/02 p = sort(p(:)); V = length(p); I = (1:V)'; cVID = 1; cVN = sum(1./(1:V)); pID = p(max(find(p<=I/V*alpha/cVID))); pN = p(max(find(p<=I/V*alpha/cVN))); return function h=myabline(b,m,varargin) % FORMAT h = abline(b,m,...) % Plots y=a*x+b in dotted line % % ... Other graphics options % % Like Splus' abline % % To plot a horizontal line at y: abline('h',y) % To plot a vertical line at x: abline('v',x) % % @(#)abline.m 1.4 02/02/13 if (nargin==2) & isstr(b) b = lower(b); else if (nargin<1) b = 0; end if (nargin<2) m = 0; end end XX=get(gca,'Xlim'); YY=get(gca,'Ylim'); if isstr(b) & (b=='h') g=line(XX,[m m],'LineStyle',':',varargin{:}); elseif isstr(b) & (b=='v') g=line([m m],YY,'LineStyle',':',varargin{:}); else g=line(XX,m*XX+b,'LineStyle',':',varargin{:}); end if (nargout>0) h=g; end function df = sf_GetDF(V) % % Infer DF from comment field of statistic image % Fdf = sscanf(V.descrip,'SPM{F_[%f,%f]}'); Tdf = sscanf(V.descrip,'SPM{T_[%f]}'); if (length(Fdf) == 2) df = Fdf'; elseif (length(Tdf)==1) df = Tdf; else df = []; end return function P = sf_spm_get(n,wildcard,prompt) switch spm('ver') case 'SPM99' P = spm_get(n,wildcard,prompt); case 'SPM2' if max(findstr(wildcard,'*img')) == length(wildcard)-length('*img')+1 wildcard = [wildcard(1:end-4) '*IMAGE']; end P = spm_get(n,wildcard,prompt); case {'SPM5','SPM8'} wildcard = strrep(wildcard,'*','.*'); P = spm_select(n,'image',prompt,'',pwd,wildcard); end