private function to compute statistics for SPM2 design FORMAT [con stat Ps Pc] = pr_stat_compute(SPM,Ic) Input SPM - SPM design structure Ic - indices into contrast structure (xCon in SPM) Output con - contrast value (ess for F test) stat - statistic value Ps - uncorrected p value Pc - P value Bonferroni corrected for number of columns analyzed Based on: @(#)spm_contrasts.m 2.3 Andrew Holmes, Karl Friston & Jean-Baptiste Poline 02/12/30 $Id$
0001 function [con,stat,Ps,Pc] = pr_stat_compute(SPM,Ic) 0002 % private function to compute statistics for SPM2 design 0003 % FORMAT [con stat Ps Pc] = pr_stat_compute(SPM,Ic) 0004 % 0005 % Input 0006 % SPM - SPM design structure 0007 % Ic - indices into contrast structure (xCon in SPM) 0008 % 0009 % Output 0010 % con - contrast value (ess for F test) 0011 % stat - statistic value 0012 % Ps - uncorrected p value 0013 % Pc - P value Bonferroni corrected for number of columns analyzed 0014 % 0015 % Based on: 0016 % @(#)spm_contrasts.m 2.3 Andrew Holmes, Karl Friston & Jean-Baptiste Poline 02/12/30 0017 % 0018 % $Id$ 0019 0020 %-Get contrast definitions (if available) 0021 %----------------------------------------------------------------------- 0022 try 0023 xCon = SPM.xCon; 0024 catch 0025 xCon = []; 0026 end 0027 0028 %-set all contrasts by default 0029 %----------------------------------------------------------------------- 0030 if nargin < 2 0031 Ic = 1:length(xCon); 0032 end 0033 if any(Ic > length(xCon)) 0034 error('Indices too large for contrast structure'); 0035 end 0036 0037 % OLS estimators and error variance estimate 0038 %---------------------------------------------------------------- 0039 betas = SPM.betas; 0040 Hp = SPM.ResidualMS; 0041 0042 %-Compute contrast and statistic parameters 0043 %======================================================================= 0044 df = [NaN SPM.xX.erdf]; 0045 for i = 1:length(Ic) 0046 0047 %-Canonicalise contrast structure with required fields 0048 %------------------------------------------------------------------- 0049 ic = Ic(i); 0050 X1o = spm_FcUtil('X1o',xCon(ic),SPM.xX.xKXs); 0051 [trMV,trMVMV] = spm_SpUtil('trMV',X1o,SPM.xX.V); 0052 df(1) = trMV^2/trMVMV; % eidf 0053 0054 switch(xCon(ic).STAT) 0055 0056 case {'T'} %-Implement contrast as sum of betas 0057 0058 con(i,:) = xCon(ic).c'*betas; 0059 VcB = xCon(ic).c'*SPM.xX.Bcov*xCon(ic).c; 0060 stat(i,:) = con(i,:)./sqrt(Hp*VcB); 0061 Ps(i,:) = 1 - spm_Tcdf(stat(i,:),df(2)); 0062 0063 case 'F' %-Implement ESS 0064 0065 %-Residual (in parameter space) forming mtx 0066 %----------------------------------------------------------- 0067 h = spm_FcUtil('Hsqr',xCon(ic),SPM.xX.xKXs); 0068 con(i,:) = sum((h * betas).^2, 1); 0069 stat(i,:) = con(i,:) ./ Hp / trMV; 0070 Ps(i,:) = (1 - spm_Fcdf(stat(i,:),df)); 0071 0072 otherwise 0073 %--------------------------------------------------------------- 0074 error(['unknown STAT "',xCon(ic).STAT,'"']) 0075 0076 end % (switch(xCon...) 0077 end 0078 0079 % Compute corrected Bonferroni (corrected for number of regions) 0080 n = size(betas, 2); 0081 Pc = 1-(1-Ps).^n; 0082