Home > marsbar > @mardo_5 > private > pr_stat_compute_mv.m

pr_stat_compute_mv

PURPOSE ^

private function to compute mutlivariate statistics across ROIs

SYNOPSIS ^

function [MVres] = pr_stat_compute_mv(SPM,Ic)

DESCRIPTION ^

 private function to compute mutlivariate statistics across ROIs
 FORMAT [MVres] = pr_stat_compute_mv(SPM,Ic)
 
 Input
 SPM       - SPM design structure
 Ic        - indices into contrast structure (xCon in SPM)
  
 Output
 MVres     - mulitvariate result structure

 $Id: pr_stat_compute_mv.m 92 2004-01-07 08:37:16Z matthewbrett $

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function [MVres] = pr_stat_compute_mv(SPM,Ic)
0002 % private function to compute mutlivariate statistics across ROIs
0003 % FORMAT [MVres] = pr_stat_compute_mv(SPM,Ic)
0004 %
0005 % Input
0006 % SPM       - SPM design structure
0007 % Ic        - indices into contrast structure (xCon in SPM)
0008 %
0009 % Output
0010 % MVres     - mulitvariate result structure
0011 %
0012 % $Id: pr_stat_compute_mv.m 92 2004-01-07 08:37:16Z matthewbrett $
0013   
0014 %-Get contrast definitions (if available)
0015 %-----------------------------------------------------------------------
0016 try
0017   xCon  = SPM.xCon;
0018 catch
0019   xCon  = [];
0020 end
0021 
0022 %-set all contrasts by default
0023 %-----------------------------------------------------------------------
0024 if nargin < 2
0025   Ic    = 1:length(xCon);
0026 end
0027 if any(Ic > length(xCon))
0028   error('Indices too large for contrast structure');
0029 end
0030 
0031 
0032 % Get relevant fields from design
0033 xCon = xCon(Ic);
0034 Xs = SPM.xX.xKXs;
0035 V = SPM.xX.V;
0036 betas = SPM.betas;
0037 ResidualMS = SPM.ResidualMS;  
0038 Y = summary_data(SPM.marsY);
0039 
0040 % setup calculation
0041 [nBetas nROI]   = size(betas);
0042 nCon          = length(xCon);
0043 [trRV trRVRV] = spm_SpUtil('trRV',Xs,V);
0044 erdf = trRV^2/trRVRV;
0045 RMS = sqrt(ResidualMS);
0046 
0047 %--------------------------------------------------------------------
0048 %- Multivariate analysis
0049 %--------------------------------------------------------------------
0050 
0051 MVres = struct('y_pre',[], 'y_obs', [], 'Pf', [], 'u', [], 'ds', [] );
0052 
0053 if nCon == 1, return, end
0054 
0055 YpY     = Y'*Y;
0056 
0057 for ii = 1:nCon
0058 
0059      xC  = xCon(ii); 
0060 
0061      %--------------------------------------------------------------------
0062      [NF, nu, h, d, M12, XG, sXG] = sf_model_mlm(Xs, V, nROI, xC, erdf);
0063 
0064      %--------------------------------------------------------------------
0065      %- Compute svd
0066      %--------------------------------------------------------------------
0067      %- fprintf('%-40s\n','Computing Principal Components')
0068 
0069      Z    = ((NF*betas)./(ones(size(NF,1),1)*RMS));
0070      S    = Z*Z';
0071      S    = S/sum(nROI);
0072      [u s u] = svd(S,0);
0073      ds    = diag(s);
0074      clear  s;
0075 
0076 
0077      %--------------------------------------------------------------------
0078      %- STATISTICS if any ...
0079      %--------------------------------------------------------------------
0080      %- Fq : F values for the last q eigein values.
0081      %- P  : P values.for the last q eigein values.
0082 
0083      Fq = zeros(1,h);
0084      for q = 0:h-1;
0085               nu1    = d*(h-q);
0086              nu2    = d*nu - (d-1)*(4*(h-q)+2*nu)/(h-q+2);
0087              Fq(q+1)    = ((nu-2)/nu) * nu2/(nu2-2)*sum(ds(q+1:h))/(h-q);
0088      end
0089      Pf    = 1 - spm_Fcdf(Fq,nu1,nu2);
0090 
0091 
0092      %- fprintf('%-40s\n','Computing predicted and observed temporal reponse')
0093 %keyboard
0094 
0095      y_pre    = (pinv(XG)'* M12 * u)*diag(sqrt(ds)); % predicted temporal reponse
0096      
0097      gV = (diag(1./sqrt(ds))*Z)'*u;
0098      y_obs    = (Y./(ones(size(Y,1),1)*RMS)/nROI)*gV;
0099      
0100     %- save results for this constrast
0101     MVres(ii).y_pre  = y_pre;
0102     MVres(ii).y_obs  = y_obs;
0103     MVres(ii).Pf     = Pf;
0104     MVres(ii).u      = u;
0105     MVres(ii).ds     = ds;
0106     MVres(ii).df     = [nu1 nu2];    
0107     
0108 end
0109 
0110 
0111 
0112 
0113 
0114 
0115 %===================================================================
0116 function [NF,nu,h,d,M12,XG,sXG] = sf_model_mlm(Xs, V, nROI, xC, erdf);
0117 % Set sub-space of interest and the related matrix of normalisation.
0118 % FORMAT [NF,nu,h,d,M12,XG] = mm_model();
0119 %- nu, h, d : degrees of freedom
0120 %- NF : matrix of normalisation
0121 %===================================================================
0122 
0123 
0124 %--------------------------------------------------------------------
0125 %- SET, COMPUTE,NORMALIZE SPACES OF INTEREST
0126 %--------------------------------------------------------------------
0127 %- set X10 and XG
0128 %- XG= X -PG(X), PG projection operator on XG (cf. eq 1, 2)
0129 %--------------------------------------------------------------------
0130 sX1o    = spm_sp('set',spm_FcUtil('X1o',xC,Xs));
0131 sXG    = spm_sp('set',spm_FcUtil('X0',xC,Xs));
0132 X1o     = spm_sp('oP',sX1o,Xs.X);
0133 XG      = spm_sp('r',sXG,Xs.X);
0134 
0135 %- Compute Normalized effexts : M1/2=X'G*V*XG (cf eq 3)
0136 %--------------------------------------------------------------------
0137 % warning off;
0138 up    = spm_sp('ox',sX1o); ; %- PG=up*up'
0139 qi    = up'*Xs.X;
0140 sigma    = up'*V*up;
0141 M12    = (chol(sigma)*qi)';
0142 M_12    = pinv(M12);
0143 
0144 %- Compute NF : normalise factor (cf eq 4)
0145 %--------------------------------------------------------------------
0146 NF    = M_12*spm_sp('X',Xs)'*spm_sp('r',sXG,spm_sp('X',Xs));
0147 
0148 %- degrees of freedom
0149 %- nROI : number of ROI (corresponds to the number of Resels)
0150 %--------------------------------------------------------------------
0151 d    = nROI*(4*log(2)/pi)^(3/2);
0152 h    = sX1o.rk; %-rank of the sub-space of interest.
0153 nu    = erdf;

Generated on Thu 17-Mar-2016 01:14:46 by m2html © 2005