pr_stat_compute

 calculates contrast value, stats and p values for contrasts
 FORMAT [Num, Stat, P, Pc] = pr_stat_compute(xCon, Xs, V, betas, ResidualMS);
 
 xCon      - contrast structure
 Xs        - design matrix
 V         - covariance matrix
 betas     - parameter estimates
 ResidualMS       - root mean square of residuals
 
 Output
 Num       - contrast value (ess for F test)
 Stat      - statistic value
 P         - uncorrected p value
 Pc        - P value corrected for number of columns analyzed
--------------------------------------------------------------------

 $Id$

0001 function [Num, Stat, P, Pc] = pr_stat_compute(xCon, Xs, V, betas, ResidualMS);
0002 % calculates contrast value, stats and p values for contrasts
0003 % FORMAT [Num, Stat, P, Pc] = pr_stat_compute(xCon, Xs, V, betas, ResidualMS);
0004 %
0005 % xCon      - contrast structure
0006 % Xs        - design matrix
0007 % V         - covariance matrix
0008 % betas     - parameter estimates
0009 % ResidualMS       - root mean square of residuals
0010 %
0011 % Output
0012 % Num       - contrast value (ess for F test)
0013 % Stat      - statistic value
0014 % P         - uncorrected p value
0015 % Pc        - P value corrected for number of columns analyzed
0016 %--------------------------------------------------------------------
0017 %
0018 % $Id$
0019 
0020 nROI = size(betas,2);
0021 nCon = length(xCon);
0022 
0023 xpxm = spm_sp('xpx-',Xs);
0024 xpVx = Xs.X'*V*Xs.X;
0025 Bcov = xpxm*xpVx*xpxm;
0026 
0027 [trRV trRVRV] = spm_SpUtil('trRV',Xs,V);
0028 erdf = trRV^2/trRVRV;
0029 RMS = sqrt(ResidualMS);
0030 
0031 T_indices = [];
0032 F_indices = [];
0033 dfnum     = [];
0034 
0035 Stat    = zeros(nCon, nROI);
0036 P    = zeros(nCon, nROI);
0037 check_Tvalue = zeros(nCon, nROI);
0038 
0039 for ii = 1:nCon
0040 
0041 %    [edf_tsp edf_Xsp] = spm_FcUtil('FconEdf', xCon(ii), Xs, V);
0042 
0043     switch(xCon(ii).STAT)
0044        case 'T'
0045         %----------- to check calculation with h -----------
0046         h       = spm_FcUtil('Hsqr', xCon(ii), Xs);
0047         [trMV trMVMV]= spm_SpUtil('trMV', ...
0048                 spm_FcUtil('X1o',xCon(ii),Xs), V);
0049 
0050         check_Tvalue(ii,:) = ((h*betas)/trMV)./RMS;
0051         check_dfnum(ii) = trMV;
0052 
0053         %----------- t value  -----------
0054 
0055         Num(ii,:) = xCon(ii).c'*betas;
0056             Stat(ii,:) = Num(ii,:) ./ ...
0057               (RMS .* sqrt((xCon(ii).c'*Bcov*xCon(ii).c)));
0058         P(ii,:) = 1 - spm_Tcdf(Stat(ii,:), erdf);
0059         T_indices = [T_indices ii];
0060 
0061 
0062            case 'F'  
0063         [trMV trMVMV]= spm_SpUtil('trMV', ...
0064                 spm_FcUtil('X1o',xCon(ii),Xs), V);
0065         dfnum   = [dfnum trMV^2/trMVMV];
0066         h       = spm_FcUtil('Hsqr', xCon(ii), Xs);
0067 
0068         Num(ii,:) = sum( (h*betas).^2, 1 );
0069         Stat(ii,:) = (Num(ii,:)/trMV) ./ (RMS.^2);
0070         check_Tvalue(ii,:) = Stat(ii,:) ;
0071         P(ii,:) = 1 - spm_Fcdf(Stat(ii,:), ...
0072                     [dfnum(end) erdf]);
0073         F_indices       = [F_indices ii];
0074         
0075 
0076        otherwise
0077             error(['unknown STAT "',xCon(ii).STAT,'"'])
0078     end
0079 end
0080 
0081 %- corrected P value for the number of ROI
0082 %--------------------------------------------------------------------
0083 Pc = ones(nCon, nROI) - (ones(nCon, nROI) - P).^nROI;